BACKGROUND
since in  <dig> the human genome project released the complete human genome sequence, there is great interest in the complex interplay between different genes and proteins.
instead of focusing on individual cellular components, interest has shifted to the interplay between these components, introducing the view that a biological system is more than the sum of its parts.
one of the most difficult problems in systems biology is the reconstruction of gene regulatory networks from experimental data.
this difficulty arises from numerous sources, among them the combinatorial explosion of possible network topologies for a given number of genes, limited information content and high levels of noise in experimental data, limited amounts of data, and the complexity of regulatory processes in cells during transcription, translation and post-translation.
many approaches have been proposed to infer networks from data, good reviews are, for example,  <cit> .
a common method to represent dynamics in biochemical systems are differential equations  <cit> .
rich mathematical theory has been established for their solution and analysis, and can be exploited  <cit> .
linear differential equation models have been proposed to infer gene regulatory networks  <cit> .
these are attractive models due to the low number of parameters and their analytical tractability.
however, since biological networks are typically highly nonlinear, linear differential equations are usually not adequate to accurately capture a regulatory network's dynamic behavior  <cit> .
some authors argue that if a system is linearized around a specific point of interest, e.g., a steady state, one may describe the local behavior using linear models  <cit> .
to describe more complex dynamic behavior, nonlinear models are needed.
such models can describe nonlinear behavior such as oscillations, multi-stationarity and biochemical switches.
furthermore, by using differential equations which are based on chemical reaction kinetics, model parameters directly correspond to reaction rates, thus models and model parameters can be immediately interpreted biochemically  <cit> .
on the other hand, due to the high-dimensional search space, inference of nonlinear models from data is much more complex than linear system identification, and serious problems with over-fitting and non-identifiability arise.
nevertheless, nonlinear models are increasingly being used, and are very likely to play an important role in our ability to understand progressively more complex systems in the future.
bongard and lipson recently published a method that can be used to symbolically infer nonlinear systems without prior specification of a model class, which they applied to simulated data of a three-component model of the lac operon  <cit> .
while such a model-free approach is very interesting, it remains to be seen whether the methodology can be extended to larger networks.
making assumptions about the underlying model class, spieth et al.
used s-systems  <cit> , generalized linear models  <cit>  and so-called h-systems and inferred models with up to  <dig> genes from data, using different search strategies, including evolutionary algorithms  <cit> .
a cooperative, coevolutionary algorithm was used by kimura et al.
for the inference of s-system models of genetic networks  <cit> .
perkins et al.
used partial differential equations to reverse engineer the gap gene network in drosophila melanogaster  <cit> .
busch et al.
recently used an approach related to delay differential equations to infer the regulatory network underlying keratinocyte migration  <cit> .
while models based on differential equations provide a quantitative dynamic description of a system under consideration, they completely disregard the stochastic nature of biological data.
linear stochastic differential equations have been proposed for this reason  <cit> , but they still require strong assumptions, and it is unclear if larger, nonlinear stochastic differential equation models of genetic regulatory networks can successfully be inferred from experimental data.
a further difficulty with differential equation models is, that it is not straightforward to compute probability distributions over alternative models or model parameters.
this would be most useful in particular if several alternative models fit the data well, and could be used to design additional experiments.
furthermore, such information would make it possible to consider alternative scenarios also in simulation-based perturbation studies, e.g., when interest is on the effect of potential drug candidates.
the problems of over-fitting and non-identifiability of models typically encountered with nonlinear differential equation models can be addressed by regularization  <cit> , or by including additional biological knowledge in the inference process  <cit> .
however, the former requires setting a regularization parameter, which is often a nontrivial problem, whereas the latter approach requires a systematic way to include such information in parameter estimation.
both issues can nicely be addressed in a bayesian framework.
we therefore embed a nonlinear ordinary differential equation  model, based on chemical reaction kinetics, into a bayesian framework.
network inference then amounts to evaluating the posterior distribution over models and model parameters, given the experimental data.
a related idea has recently been pursued by steinke et al.
for linear models  <cit> .
in their paper, the authors use expectation propagation to evaluate the posterior distribution.
we combine a nonlinear differential equation model based on chemical reaction kinetics with a bayesian framework, and use a markov chain monte carlo  approach to sample from the posterior distribution.
a difficulty with this approach comes from the fact, that it requires solving the differential equations at every step in the markov chain.
to avoid this problem, we adapt the parameter estimation method proposed by ramsay et al.
<cit>  for markov chains.
the authors iteratively fit smoothing splines to the data, and then learn the parameters of the differential equations using a least squares procedure on the slope estimates of the splines and the differential equations.
the idea to carry out the optimization on slopes instead of concentrations was first suggested by varah in  <dig>  <cit> .
it was then improved by poyton et al., who proposed to iterate between spline interpolation and parameter estimation  <cit> .
ramsay et al.
further improved this approach by proposing profiled estimation  <cit> .
we adapt their objective function, and sample both model parameters and smoothing factor using a markov chain.
the combination of differential equations with a bayesian framework proposed in this work allows it to adequately describe the nonlinear dynamic behavior of gene regulatory networks, and to incorporate prior knowledge into the network inference process at the same time.
in contrast to simple optimization of the posterior distribution as we pursued in previous work  <cit> , the mcmc approach used here provides confidences on learned parameters and computes probability distributions over alternative network topologies.
this can be used to consider alternative future scenarios in simulation, and permits the design of additional, most informative experiments to improve the inference procedure.
methods
differential equations model
we represent gene regulatory networks as directed graphs g = , with vertices v = {v <dig> ..., vn} corresponding to genes and directed edges e corresponding to regulatory interactions.
an edge from gene i to gene j indicates that the product of gene i, xi, influences the expression of gene j either by activating or by inhibiting it.
we assume that different regulators act independently, such that the total effect on the expression of gene i can be written as the sum of the individual effects.
this clearly is a simplification, and can be generalized by considering products of effects from different genes.
our ode model is written as
where xi is the concentration of gene i at time t. furthermore, si and γi are basal synthesis and degradation rates for each gene i, which in the absence of regulations from other genes determine the dynamic behavior of component i. variable βij denotes the regulation strength of component xj on xi, and fij is the corresponding regulation function.
βij >  <dig> corresponds to an activation, βij <  <dig> to an inhibition, and βij =  <dig> means that there is no regulation from gene j to gene i. as regulation functions we use hill-type functions
where m denotes the hill coefficient and θj is related to the reaction rate by describing the concentration of xj needed to significantly activate or inhibit expression of xi.
to keep the number of parameters small, we use one joint hill coefficient m for all regulations, and use the same threshold parameter θj for all regulations out of the same gene xj.
the regulation functions  are obtained from chemical reaction kinetics by considering the binding process of a transcription factor to a promoter a reversible chemical reaction.
for a more detailed derivation, see, for example,  <cit> .
parameter estimation of ode systems
the estimation of model parameters from experimental time series data for differential equation models is typically carried out iteratively in two steps:  numerically solve the differential equations for the time interval of interest, and  compute an error between experimental data points and model prediction.
initial values and model parameters are then modified to minimize this error.
the disadvantage of this procedure is, that the differential equations have to be solved numerically in every iteration of the optimization, which is very time consuming.
as an alternative, varah proposed a two stage method  <cit> .
in step one, interpolating cubic splines z are fitted to the data d. thereafter, the squared difference between the differential equations and the slope estimates from the interpolation is minimized according to
here,  is the slope estimate from the cubic splines z, ω are the differential equation model parameters, t is the number of different time points tτ and n is the number of time series to be fitted, for example, the number of different genes in the network.
an obvious drawback of varah's approach is, that the quality of the parameter estimates can only be as good as the spline fit, which is particularly difficult in case of noisy data.
to address this problem, poyton et al.
developed a recursive method, called iterative principal differential analysis, where the two steps of varah's method are iterated, and the model predictions are fed back into the spline estimation  <cit> .
ramsay et al.
improved this method further using a generalization of profiled estimation to learn the parameters of interest  <cit> .
we adapt this iterative method by simultaneously estimating model parameters ω :=  and smoothing factor λ of the smoothing splines.
this could be done by minimization of the function
where diτ denotes the measured data, t <dig> is the number of time points in the experimental data and t <dig> denotes the number of points to be used in the squared error parameter fitting on the slopes.
to adequately describe the dynamics of a system using derivatives, a large number of slope estimates  is required, we will therefore usually have t <dig> >> t <dig>  we note that equation  requires concentrations to compute the derivatives, these are taken from the spline interpolation.
bayesian learning framework
we now address two further problems in parameter estimation, regarding the entire topology of the regulatory network, and variability in experimental measurements.
the network topology  can either be determined in a separate step prior to parameter estimation, or can be solved implicitly by assuming a fully connected network, and pruning edges with very small regulation strengths afterwards.
determination of the network topology in a separate step has the disadvantage, that edges not included in this prior step cannot be re-introduced in parameter estimation.
we therefore use the latter approach, with appropriate regularization to prune many edges during the inference process.
to account for noise in the experimental data, we embed the differential equations into a bayesian framework.
for this purpose, we assume that the measured data diτ is corrupted by independent mean zero gaussian noise with variance .
the assumption of normally distributed noise is clearly a simplification, which is made here to keep the model simple.
other noise models could be used.
we furthermore assume the differences between slope estimates and differential equation predictions to follow a normal distribution with mean zero and variance .
we note that the ratio between  and  corresponds to a parameter that weighs the two terms in  relative to one another.
the assumption of gaussian noise leads to the likelihood
which is equivalent to equation  up to log-transformation and scaling.
since we are interested in the probability distribution over the model parameters ω, the smoothing factor and the variances  and , given the experimental data d, we use bayes' theorem to obtain
where  is given by the likelihood ,  is a prior distribution on the model parameters, and p is a normalizing factor which is independent of ω, λ  and .
for simplicity, we treat the variance parameters  and  as user parameters, which are set in advance and not sampled.
inclusion of prior knowledge
the prior distribution p on the model parameters can be used to include available biological knowledge on the system under consideration into the network inference process, as demonstrated, for example, by wehrli and husmeier  <cit> .
it is a huge advantage of the bayesian framework that it allows the easy and systematic integration of such expert knowledge.
if no such detailed knowledge is available, one can resort to very general assumptions, such as sparsity of the interaction network  <cit>  or rough estimates of reasonable ranges for parameters.
we assume independent prior distributions for the different model parameters, and suggest to use gamma priors
for the synthesis and degradation rates si and γi, for the hill coefficient m and the threshold parameters θj.
the parameters a and r are scale and rate parameters of the gamma distribution, respectively, and Γ denotes the gamma function
this choice of prior ensures that the parameters are positive, and will not become too large.
since the smoothing factor λ ranges between zero and one, we use a beta prior
on λ. the parameters a' and b' are positive shape parameters of the beta distribution, and b is the beta function
to reflect the assumption of sparsity of gene regulatory networks, we use a prior based on the lq norm  <cit>  for the interaction strengths βij, i, j =  <dig> ..., n:
for β ∈ ℝ and q, s >  <dig>  where  is the normalizing factor
for q =  <dig>  equation  is a normal distribution, for q =  <dig> it corresponds to the laplace distribution.
values of q < <dig> enforce stronger sparsity constraints, as can be seen in figure  <dig> for the two-dimensional case with q =  <dig>  and s =  <dig>  in comparison with the prior proposed in  <cit>  and used in network inference in  <cit> , we avoid the numerical integration of the prior required in these publications, and obtain similar sparseness constraints.
mcmc sampling from the posterior
the posterior distribution  can now be maximized using, for example, gradient based methods, simulated annealing or genetic algorithms.
however, since multiple parameter combinations, corresponding to alternative network topologies, may explain the data equally well, we sample from p using markov chain monte carlo.
this way, full distributions over each parameter are available, and can be used, for example, to consider different likely topologies, and to design experiments that will resolve ambiguities.
this would not be possible with simple maximization approaches.
to sample from p we use an iterative approach.
first, we sample the model parameters ω using the hybrid monte carlo algorithm , with fixed smoothing factor λ. hmc has originally been proposed by duane et al.
for problems arising in quantum chromodynamics  <cit> , and has been introduced to the general bayesian statistical community by neal  <cit> .
the method samples points from a given n-dimensional distribution p by introducing momentum variables ρ =  with associated energy k, and iterative sampling for the momentum variables from k and following the hamiltonian dynamics of the system h:= -ln p + k in time.
doing so, hmc generates a sequence of points distributed according to p, and can avoid the random walk behavior typically observed with the metropolis hastings algorithm  <cit> .
as the second step, we sample the smoothing factor λ using metropolis hastings  <cit> , with ω fixed to the values sampled in the previous step.
pseudocode for our iterative sampling procedure is given in table  <dig>
evaluation of reconstructed networks
to evaluate reconstructed networks, we summarize the markov chains by using the mean value of the chain for each model parameter, after excluding points from the burn-in phase of the chain.
this is of course a very crude simplification, which we take to allow for an automated, quantitative evaluation of reconstructed networks.
obviously, in case of, for example, bimodal distributions, the mean will be located somewhere between the two modes, possibly in a region with very low probability mass.
we therefore emphasize here that the full set of points sampled can and should be analyzed in more detail.
to quantitatively evaluate inferred networks, we use receiver operator characteristic  and precision to recall  analysis, and summarize these using the area under the curve .
for two-class classification problems , roc curves plot sensitivity against 1-specificity for varying thresholds on the predictor , whereas pr curves similarly plot precision against recall .
the auc is then simply the area under the roc or pr curve, and on a scale from  <dig> to  <dig> provides a single number to measure the quality of a predictor.
we note that changing the threshold in roc and pr curves corresponds to different thresholds for edge pruning in reconstructed networks.
in our case, we want to distinguish between three classes - no edge, activation, or inhibition.
therefore, we map the three-class problem onto a two-class problem as indicated in table  <dig>  with this approach we calculate sensitivity, specificity and precision, to calculate the area under the roc and pr curves .
we point out that for our three-class problem, for a random network, the average expected auc value will not be  <dig>  as in the two-class case, but will vary between zero and  <dig>  for aucroc and between zero and  <dig>  for aucpr, depending on the number of non-existing edges in the reference network.
for a mathematical proof we refer to additional file  <dig>
mapping of three-class classification problem  onto two-class roc/pr evaluation.
tp: true positive, tn: true negative, fp: false positive, fn: false negative.
implementation
we implemented our algorithm in matlab, release 2008b , using the spline and statistics toolboxes.
computations were carried out on a linux cluster with dual-processor  <dig>  ghz xeon quadcore machines with  <dig> gb ram, running each markov chain in a single thread .
the matlab code is available on request from the authors.
RESULTS
simulated five gene network
we first evaluated our approach on a simulated five-gene network.
this allows it to systematically study the performance of network inference under varying levels of noise and dataset sizes, while the true network topology is known.
we simulated data using the network topology shown in figure 3a.
since this is the topology also underlying the experimental data used later, this allows a direct comparison of simulation results with inference results on real experimental data.
data was simulated using the differential equation model , with synthesis and degradation rates s =  and γ =  for the five genes.
the threshold parameter of the hill functions was set to θ = , with hill coefficient m =  <dig>  the parameters for the regulation strength were set to βij =  <dig> for activations, βij = - <dig> for inhibitions, and zero if there was no interaction between to genes, compare figure 3a.
data was simulated by numerical integration of the differential equations in matlab using the function ode <dig>  simulated data shows oscillations for all genes, see figure 3b.
to simulate the typical setting in network inference, where only a limited number of noisy measurements are available, we evaluated our network reconstruction approach using different numbers of time points subsampled equidistantly from the simulated data, and added mean-zero normally distributed noise with different standard deviations to the concentration values.
we then used our method to reconstruct the original network from this data.
for this purpose, we sampled  <dig>   <dig> points for  using the algorithm described in table  <dig>  with a burn-in of  <dig>   <dig> points.
parameters of the prior distributions were set to a =  <dig>  r =  <dig> for the gamma prior on synthesis and degradation rates, a =  <dig> , r =  <dig>  for the gamma prior on the hill coefficient m, a' =  <dig>  b' =  <dig> for the beta prior on λ, and q =  <dig> , s =  <dig> for the lq prior on the regulation strengths βij.
shape and scale parameters for the gamma priors on the θj for each gene j where chosen such that mean and variance of the priors correspond to mean and variance of the training data.
the number of slope estimates t <dig> is set to  <dig> and the corresponding variance  is set to  <dig>  furthermore, the variance  is set to t1/t <dig>  where t <dig> denotes the number of time points.
results on  <dig> time points
we first describe results for an ideal case with  <dig> time points and no noise.
in that case, mean values for inferred synthesis and degradation rates were s =  and γ = .
mean value for the hill coefficient m was  <dig> , means for the thresholds θj ranged from  <dig>  to  <dig>  and the mean smoothing factor λ was  <dig> .
inferred regulation strengths  are given in table  <dig>  the large standard deviations for some regulation strengths, e.g., the self-regulation on gene  <dig> or the regulation from gene  <dig> to gene  <dig>  indicate that there either are different network topologies which describe the data well, or that the dynamics of the system is insensitive to changes of this parameter.
learned regulation strength parameters β for the simulated dataset with  <dig> time points.
given are mean and standard deviation of the sampled interaction parameters.
true edges which are present in the reference network are indicated in bold.
precision to recall and receiver operator characteristic analysis of results yield auc values of aucpr =  <dig>   for precision to recall curves and aucroc =  <dig>   for sensitivity vs. 1-specificity curves.
to close the circle from the original concentration data over the reconstructed model back to dynamic simulation, we used the mean inferred model parameters to simulate gene concentrations.
this simulation shows an accurate match of simulated and experimental data .
it is well known that fitting models to oscillating data, and even more so, reconstructing networks from such data, are extremely hard problems, since models tend to learn a steady state  <cit> .
in spite of this, oscillations were captured with high precision by our approach .
effect of noise and dataset size
we next studied the effect of different dataset sizes  and different levels of noise in the data on the quality of network reconstruction.
for this purpose, we added mean zero gaussian noise with standard deviations  <dig> ,  <dig> ,  <dig> ,  <dig>  and  <dig>  to the simulated concentration data, and furthermore subsampled equidistantly from the time series to generate data sets with t <dig> =  <dig>   <dig>   <dig>   <dig>   <dig>   <dig>   <dig>   <dig>   <dig>   <dig> and  <dig> different time points for each of the five noise levels.
then the network reconstruction was performed as described in the methods section.
evaluation on experimental data: the dream  <dig>  challenge  <dig> dataset
to assess the performance of different reverse engineering approaches, stolovitzky et al.
fathered the dream  initiative  <cit> .
for this purpose, cantone et al.
provided in-vivo data on a small, bio-engineered five-gene network, which was posted as challenge # <dig> within dream  <dig>  <cit> .
this data was generated by inserting new promoter/gene combinations directly into the chromosomal dna of budding yeast.
two time series of gene expression of the five inserted genes after stimulation were measured using quantitative pcr, measuring  <dig> time points in  <dig> minute intervals in time series  <dig>  and  <dig> time points in  <dig> minute intervals in time series  <dig>
measurements in both time series consist of negative  log-ratios of the genes of interest to housekeeping genes, we therefore transformed the measured data to recover the original ratios.
we used the first time series  for network inference, i.e., t <dig> =  <dig>
the hyperparameters for the lq prior on the regulation strengths were set to q =  <dig> and s =  <dig>  parameters for the gamma priors on the threshold parameters were manually set to concentrate probability mass near the mean concentration value for each gene individually.
an evaluation of results using the mean of the markov chain for each parameter, as done for the simulated data, results in auc values that are equivalent to guessing .
this might indicate that either the level of noise present in the experimental data is too high, or that the posterior distribution has multiple modes, with the mean being an inappropriate summary statistic.
we therefore searched the points sampled for the maximum a-posteriori mode, and evaluate this mode further in the following.
clearly, data from the additional modes are available and can be studied similarly.
the regulation strength parameters for the maximum a-posteriori mode are shown in table  <dig>  the dynamic behavior and fit of the model prediction to the experimental data is depicted in figure  <dig>  obviously, the general dynamics of the original data is well represented, with a moderate amount of smoothing.
auc values of the reconstructed network are  <dig>   for sensitivity vs. 1-specificity, and  <dig>   for precision to recall.
reconstructed maximum-a-posteriori regulation strength parameters β for the dream  <dig> challenge # <dig> data.
true edges present in the reference topology are marked in bold.
to assess the quality of this result, we next compared performance of our approach to the performance of other approaches submitted to the dream  <dig> challenge.
for this purpose, we computed performance measures as were used in the dream  <dig> challenge for our inferred network, and show results in table  <dig>  auc values for this comparison were calculated as described in  <cit> .
since our method gives us a topology with positive and negative regulation strengths, we have transformed our results to be suitable for the evaluation method used in dream 2:
• by skipping the sign and dividing by the largest learned regulation strength for the directed-unsigned challenge.
• for the two directed-signed challenges we only took the regulation strengths with the appropriate sign and divide them by the highest absolute regulation strength.
results of the dream  <dig> challenge # <dig> data of our method  compared to submitted best results from  <cit>  .
the third column gives the auc values for our method when self-regulations are omitted.
our method clearly outperforms all submitted methods in the directed-signed-inhibitory challenge; furthermore, when neglecting self-regulations, we also beat the best submitted method in the directed-unsigned challenge.
our method outperforms all submitted results in the challenge directed-signed-inhibitory.
one difficulty we observed was, that our approach learned many strong self-regulations of genes, possibly because of an improper balancing of the priors on synthesis/degradation rates and regulation strengths.
since there are no self-regulations in the dream  <dig> challenge # <dig> data, we provide an additional evaluation when disregarding self-regulations, results are shown in the third column of table  <dig>  in this case, we not only outperform all submitted approaches in the directed-signed-inhibitory challenge, but also beat the best models in the directed-unsigned challenge.
discussion and
CONCLUSIONS
we have developed a novel methodological approach to reverse engineer gene regulatory networks from gene expression time series data, and evaluated this approach on both simulated and real gene expression data.
the combination of ordinary differential equations and the bayes' regularized inference technique can be used for the quantitative analysis of complex cellular processes.
non-linear differential equations are able to describe complex dynamic behavior, and a rich mathematical theory for analyzing them is well established.
our method combines these advantages of differential equations with the advantages of a bayesian framework, which is able to capture noise in data, makes it possible to include biological knowledge into the learning process, and allows the computation of probability distributions over model topologies and model parameters.
the latter is one of the main advantages of the mcmc approach.
the information about distributions can be used to make predictions of future states of the network together with confidence intervals on the predictions made.
this may allow it to take alternative future scenarios into account, and could be used to design additional most informative experiments that will help to distinguish between corresponding topologies or parameter sets.
we therefore think that our approach will be highly useful to elucidate regulatory networks in an iterative procedure with several rounds of experiment, network inference, and experiment design.
in contrast to the usual approach of minimizing an error between experimentally measured concentrations and model predictions, our likelihood function uses the difference between model slopes and experimental slopes.
we furthermore integrate a smoothing spline approximation into the likelihood, automatically performing an optimized tradeoff between an accurate representation of the experimental data, and smoothing out noise.
fitting of model parameters on slopes has the advantage that no numerical integration of the model is required in each step of the optimization or sampling.
instead, we must estimate smoothing splines and slopes.
this can be carried out much faster than numerical integration, enabling us to use a markov chain sampler on the posterior distribution instead of plain maximization.
we have evaluated our approach on simulated and on real experimental data from a synthetic gene regulatory network.
on the simulated data, we have shown that our approach can reconstruct the underlying topology with high accuracy.
as expected, performance deteriorates with increasing levels of noise and with decreasing number of different time points available.
we emphasize that the simulated example chosen is a difficult learning task due to the oscillations in the data.
it is obvious that sufficient data points are required to sample the full dynamics of the oscillating network, and that oscillations quickly break down in the presence of noise.
on the dream  <dig> data, our method yields superior results when compared to other approaches that were submitted to dream  <dig> in the directed-signed-inhibitory and directed-unsigned categories.
importantly, our analysis shows that there are multiple posterior modes that describe the data well, which may explain a surprising result of the original dream  <dig> challenge: as reported by stolovitzky et al., none of the submitted models were able to accurately reconstruct the original network topology from the synthetic data  <cit> .
our results indicate that this might be due to a dense population of local optima, in which network reconstruction approaches looking for a single optimal topology might get trapped and return suboptimal solutions.
an obvious conclusion is that further experiments are required to resolve ambiguities in network reconstruction.
this emphasizes the need for robust and efficient methods for optimum experimental design.
our sampling approach may be a good starting point for such experiment design, since it analyzes the full distribution over model parameters, and therefore yields information on alternative network topologies and confidence intervals on parameters, which are instrumental to design experiments that elucidate the network topology further.
authors' contributions
lk designed the research and methodology.
jm and dr wrote the software, analyzed the data, and drafted the paper.
gr and lk advised on algorithm design and data analysis, and wrote the final version of the paper.
all authors read and approved the final version of the paper.
supplementary material
additional file 1
supplementary material.
mathematical proof.
click here for file
acknowledgements
we would like to thank a. szabowski, h. busch, n. radde and r. eils for helpful comments and discussions.
furthermore, we thank m. p. cosma and d. di bernardo for the provision of the challenge # <dig> data from dream  <dig>  we are grateful for funding to the german ministry of education and research , grant number  <dig> .

 since in the human genome project released the complete human genome sequence, there is great interest in the complex interplay between different genes and proteins. instead of focusing on individual cellular components, interest has shifted to the interplay between these components, introducing the view that a biological system is more than the sum of its parts. one of the most difficult problems in systems biology is the reconstruction of gene regulatory networks from experimental data. this difficulty arises from numerous sources, among them the combinatorial explosion of possible network topologies for a given number of genes, limited information content and high levels of noise in experimental data, limited amounts of data, and the complexity of regulatory processes in cells during transcription, translation and post-translation. many approaches have been proposed to infer networks from data, good reviews are, for example,. a common method to represent dynamics in biochemical systems are differential equations. rich mathematical theory has been established for their solution and analysis, and can be exploited. linear differential equation models have been proposed to infer gene regulatory networks. these are attractive models due to the low number of parameters and their analytical tractability. however, since biological networks are typically highly nonlinear, linear differential equations are usually not adequate to accurately capture a regulatory network's dynamic behavior. some authors argue that if a system is linearized around a specific point of interest, e.g., a steady state, one may describe the local behavior using linear models. to describe more complex dynamic behavior, nonlinear models are needed. such models can describe nonlinear behavior such as oscillations, multi-stationarity and biochemical switches. furthermore, by using differential equations which are based on chemical reaction kinetics, model parameters directly correspond to reaction rates, thus models and model parameters can be immediately interpreted biochemically. on the other hand, due to the high-dimensional search space, inference of nonlinear models from data is much more complex than linear system identification, and serious problems with over-fitting and non-identifiability arise. nevertheless, nonlinear models are increasingly being used, and are very likely to play an important role in our ability to understand progressively more complex systems in the future. bongard and lipson recently published a method that can be used to symbolically infer nonlinear systems without prior specification of a model class, which they applied to simulated data of a three-component model of the lac operon. while such a model-free approach is very interesting, it remains to be seen whether the methodology can be extended to larger networks. making assumptions about the underlying model class, spieth et al. used s-systems, generalized linear models and so-called h-systems and inferred models with up to genes from data, using different search strategies, including evolutionary algorithms. a cooperative, coevolutionary algorithm was used by kimura et al. for the inference of s-system models of genetic networks. perkins et al. used partial differential equations to reverse engineer the gap gene network in drosophila melanogaster. busch et al. recently used an approach related to delay differential equations to infer the regulatory network underlying keratinocyte migration. while models based on differential equations provide a quantitative dynamic description of a system under consideration, they completely disregard the stochastic nature of biological data. linear stochastic differential equations have been proposed for this reason, but they still require strong assumptions, and it is unclear if larger, nonlinear stochastic differential equation models of genetic regulatory networks can successfully be inferred from experimental data. a further difficulty with differential equation models is, that it is not straightforward to compute probability distributions over alternative models or model parameters. this would be most useful in particular if several alternative models fit the data well, and could be used to design additional experiments. furthermore, such information would make it possible to consider alternative scenarios also in simulation-based perturbation studies, e.g., when interest is on the effect of potential drug candidates. the problems of over-fitting and non-identifiability of models typically encountered with nonlinear differential equation models can be addressed by regularization, or by including additional biological knowledge in the inference process. however, the former requires setting a regularization parameter, which is often a nontrivial problem, whereas the latter approach requires a systematic way to include such information in parameter estimation. both issues can nicely be addressed in a bayesian framework. we therefore embed a nonlinear ordinary differential equation model, based on chemical reaction kinetics, into a bayesian framework. network inference then amounts to evaluating the posterior distribution over models and model parameters, given the experimental data. a related idea has recently been pursued by steinke et al. for linear models. in their paper, the authors use expectation propagation to evaluate the posterior distribution. we combine a nonlinear differential equation model based on chemical reaction kinetics with a bayesian framework, and use a markov chain monte carlo approach to sample from the posterior distribution. a difficulty with this approach comes from the fact, that it requires solving the differential equations at every step in the markov chain. to avoid this problem, we adapt the parameter estimation method proposed by ramsay et al. for markov chains. the authors iteratively fit smoothing splines to the data, and then learn the parameters of the differential equations using a least squares procedure on the slope estimates of the splines and the differential equations. the idea to carry out the optimization on slopes instead of concentrations was first suggested by varah in. it was then improved by poyton et al., who proposed to iterate between spline interpolation and parameter estimation. ramsay et al. further improved this approach by proposing profiled estimation. we adapt their objective function, and sample both model parameters and smoothing factor using a markov chain. the combination of differential equations with a bayesian framework proposed in this work allows it to adequately describe the nonlinear dynamic behavior of gene regulatory networks, and to incorporate prior knowledge into the network inference process at the same time. in contrast to simple optimization of the posterior distribution as we pursued in previous work, the mcmc approach used here provides confidences on learned parameters and computes probability distributions over alternative network topologies. this can be used to consider alternative future scenarios in simulation, and permits the design of additional, most informative experiments to improve the inference procedure. methods differential equations model we represent gene regulatory networks as directed graphs g =, with vertices v = {v..., vn} corresponding to genes and directed edges e corresponding to regulatory interactions. an edge from gene i to gene j indicates that the product of gene i, xi, influences the expression of gene j either by activating or by inhibiting it. we assume that different regulators act independently, such that the total effect on the expression of gene i can be written as the sum of the individual effects. this clearly is a simplification, and can be generalized by considering products of effects from different genes. our ode model is written as where xi is the concentration of gene i at time t. furthermore, si and γi are basal synthesis and degradation rates for each gene i, which in the absence of regulations from other genes determine the dynamic behavior of component i. variable βij denotes the regulation strength of component xj on xi, and fij is the corresponding regulation function. βij > corresponds to an activation, βij < to an inhibition, and βij = means that there is no regulation from gene j to gene i. as regulation functions we use hill-type functions where m denotes the hill coefficient and θj is related to the reaction rate by describing the concentration of xj needed to significantly activate or inhibit expression of xi. to keep the number of parameters small, we use one joint hill coefficient m for all regulations, and use the same threshold parameter θj for all regulations out of the same gene xj. the regulation functions are obtained from chemical reaction kinetics by considering the binding process of a transcription factor to a promoter a reversible chemical reaction. for a more detailed derivation, see, for example,. parameter estimation of ode systems the estimation of model parameters from experimental time series data for differential equation models is typically carried out iteratively in two steps: numerically solve the differential equations for the time interval of interest, and compute an error between experimental data points and model prediction. initial values and model parameters are then modified to minimize this error. the disadvantage of this procedure is, that the differential equations have to be solved numerically in every iteration of the optimization, which is very time consuming. as an alternative, varah proposed a two stage method. in step one, interpolating cubic splines z are fitted to the data d. thereafter, the squared difference between the differential equations and the slope estimates from the interpolation is minimized according to here, is the slope estimate from the cubic splines z, ω are the differential equation model parameters, t is the number of different time points tτ and n is the number of time series to be fitted, for example, the number of different genes in the network. an obvious drawback of varah's approach is, that the quality of the parameter estimates can only be as good as the spline fit, which is particularly difficult in case of noisy data. to address this problem, poyton et al. developed a recursive method, called iterative principal differential analysis, where the two steps of varah's method are iterated, and the model predictions are fed back into the spline estimation. ramsay et al. improved this method further using a generalization of profiled estimation to learn the parameters of interest. we adapt this iterative method by simultaneously estimating model parameters ω:= and smoothing factor λ of the smoothing splines. this could be done by minimization of the function where diτ denotes the measured data, t is the number of time points in the experimental data and t denotes the number of points to be used in the squared error parameter fitting on the slopes. to adequately describe the dynamics of a system using derivatives, a large number of slope estimates is required, we will therefore usually have t >> t we note that equation requires concentrations to compute the derivatives, these are taken from the spline interpolation. bayesian learning framework we now address two further problems in parameter estimation, regarding the entire topology of the regulatory network, and variability in experimental measurements. the network topology can either be determined in a separate step prior to parameter estimation, or can be solved implicitly by assuming a fully connected network, and pruning edges with very small regulation strengths afterwards. determination of the network topology in a separate step has the disadvantage, that edges not included in this prior step cannot be re-introduced in parameter estimation. we therefore use the latter approach, with appropriate regularization to prune many edges during the inference process. to account for noise in the experimental data, we embed the differential equations into a bayesian framework. for this purpose, we assume that the measured data diτ is corrupted by independent mean zero gaussian noise with variance. the assumption of normally distributed noise is clearly a simplification, which is made here to keep the model simple. other noise models could be used. we furthermore assume the differences between slope estimates and differential equation predictions to follow a normal distribution with mean zero and variance. we note that the ratio between and corresponds to a parameter that weighs the two terms in relative to one another. the assumption of gaussian noise leads to the likelihood which is equivalent to equation up to log-transformation and scaling. since we are interested in the probability distribution over the model parameters ω, the smoothing factor and the variances and, given the experimental data d, we use bayes' theorem to obtain where is given by the likelihood, is a prior distribution on the model parameters, and p is a normalizing factor which is independent of ω, λ and. for simplicity, we treat the variance parameters and as user parameters, which are set in advance and not sampled. inclusion of prior knowledge the prior distribution p on the model parameters can be used to include available biological knowledge on the system under consideration into the network inference process, as demonstrated, for example, by wehrli and husmeier. it is a huge advantage of the bayesian framework that it allows the easy and systematic integration of such expert knowledge. if no such detailed knowledge is available, one can resort to very general assumptions, such as sparsity of the interaction network or rough estimates of reasonable ranges for parameters. we assume independent prior distributions for the different model parameters, and suggest to use gamma priors for the synthesis and degradation rates si and γi, for the hill coefficient m and the threshold parameters θj. the parameters a and r are scale and rate parameters of the gamma distribution, respectively, and Γ denotes the gamma function this choice of prior ensures that the parameters are positive, and will not become too large. since the smoothing factor λ ranges between zero and one, we use a beta prior on λ. the parameters a' and b' are positive shape parameters of the beta distribution, and b is the beta function to reflect the assumption of sparsity of gene regulatory networks, we use a prior based on the lq norm for the interaction strengths βij, i, j =..., n: for β ∈ ℝ and q, s > where is the normalizing factor for q = equation is a normal distribution, for q = it corresponds to the laplace distribution. values of q < enforce stronger sparsity constraints, as can be seen in figure for the two-dimensional case with q = and s = in comparison with the prior proposed in and used in network inference in, we avoid the numerical integration of the prior required in these publications, and obtain similar sparseness constraints. mcmc sampling from the posterior the posterior distribution can now be maximized using, for example, gradient based methods, simulated annealing or genetic algorithms. however, since multiple parameter combinations, corresponding to alternative network topologies, may explain the data equally well, we sample from p using markov chain monte carlo. this way, full distributions over each parameter are available, and can be used, for example, to consider different likely topologies, and to design experiments that will resolve ambiguities. this would not be possible with simple maximization approaches. to sample from p we use an iterative approach. first, we sample the model parameters ω using the hybrid monte carlo algorithm, with fixed smoothing factor λ. hmc has originally been proposed by duane et al. for problems arising in quantum chromodynamics, and has been introduced to the general bayesian statistical community by neal. the method samples points from a given n-dimensional distribution p by introducing momentum variables ρ = with associated energy k, and iterative sampling for the momentum variables from k and following the hamiltonian dynamics of the system h:= -ln p + k in time. doing so, hmc generates a sequence of points distributed according to p, and can avoid the random walk behavior typically observed with the metropolis hastings algorithm. as the second step, we sample the smoothing factor λ using metropolis hastings, with ω fixed to the values sampled in the previous step. pseudocode for our iterative sampling procedure is given in table evaluation of reconstructed networks to evaluate reconstructed networks, we summarize the markov chains by using the mean value of the chain for each model parameter, after excluding points from the burn-in phase of the chain. this is of course a very crude simplification, which we take to allow for an automated, quantitative evaluation of reconstructed networks. obviously, in case of, for example, bimodal distributions, the mean will be located somewhere between the two modes, possibly in a region with very low probability mass. we therefore emphasize here that the full set of points sampled can and should be analyzed in more detail. to quantitatively evaluate inferred networks, we use receiver operator characteristic and precision to recall analysis, and summarize these using the area under the curve. for two-class classification problems, roc curves plot sensitivity against 1-specificity for varying thresholds on the predictor, whereas pr curves similarly plot precision against recall. the auc is then simply the area under the roc or pr curve, and on a scale from to provides a single number to measure the quality of a predictor. we note that changing the threshold in roc and pr curves corresponds to different thresholds for edge pruning in reconstructed networks. in our case, we want to distinguish between three classes - no edge, activation, or inhibition. therefore, we map the three-class problem onto a two-class problem as indicated in table with this approach we calculate sensitivity, specificity and precision, to calculate the area under the roc and pr curves. we point out that for our three-class problem, for a random network, the average expected auc value will not be as in the two-class case, but will vary between zero and for aucroc and between zero and for aucpr, depending on the number of non-existing edges in the reference network. for a mathematical proof we refer to additional file mapping of three-class classification problem onto two-class roc/pr evaluation. tp: true positive, tn: true negative, fp: false positive, fn: false negative. implementation we implemented our algorithm in matlab, release 2008b, using the spline and statistics toolboxes. computations were carried out on a linux cluster with dual-processor ghz xeon quadcore machines with gb ram, running each markov chain in a single thread. the matlab code is available on request from the authors.  simulated five gene network we first evaluated our approach on a simulated five-gene network. this allows it to systematically study the performance of network inference under varying levels of noise and dataset sizes, while the true network topology is known. we simulated data using the network topology shown in figure 3a. since this is the topology also underlying the experimental data used later, this allows a direct comparison of simulation results with inference results on real experimental data. data was simulated using the differential equation model, with synthesis and degradation rates s = and γ = for the five genes. the threshold parameter of the hill functions was set to θ =, with hill coefficient m = the parameters for the regulation strength were set to βij = for activations, βij = - for inhibitions, and zero if there was no interaction between to genes, compare figure 3a. data was simulated by numerical integration of the differential equations in matlab using the function ode simulated data shows oscillations for all genes, see figure 3b. to simulate the typical setting in network inference, where only a limited number of noisy measurements are available, we evaluated our network reconstruction approach using different numbers of time points subsampled equidistantly from the simulated data, and added mean-zero normally distributed noise with different standard deviations to the concentration values. we then used our method to reconstruct the original network from this data. for this purpose, we sampled points for using the algorithm described in table with a burn-in of points. parameters of the prior distributions were set to a = r = for the gamma prior on synthesis and degradation rates, a =, r = for the gamma prior on the hill coefficient m, a' = b' = for the beta prior on λ, and q =, s = for the lq prior on the regulation strengths βij. shape and scale parameters for the gamma priors on the θj for each gene j where chosen such that mean and variance of the priors correspond to mean and variance of the training data. the number of slope estimates t is set to and the corresponding variance is set to furthermore, the variance is set to t1/t where t denotes the number of time points. results on time points we first describe results for an ideal case with time points and no noise. in that case, mean values for inferred synthesis and degradation rates were s = and γ =. mean value for the hill coefficient m was, means for the thresholds θj ranged from to and the mean smoothing factor λ was. inferred regulation strengths are given in table the large standard deviations for some regulation strengths, e.g., the self-regulation on gene or the regulation from gene to gene indicate that there either are different network topologies which describe the data well, or that the dynamics of the system is insensitive to changes of this parameter. learned regulation strength parameters β for the simulated dataset with time points. given are mean and standard deviation of the sampled interaction parameters. true edges which are present in the reference network are indicated in bold. precision to recall and receiver operator characteristic analysis of results yield auc values of aucpr = for precision to recall curves and aucroc = for sensitivity vs. 1-specificity curves. to close the circle from the original concentration data over the reconstructed model back to dynamic simulation, we used the mean inferred model parameters to simulate gene concentrations. this simulation shows an accurate match of simulated and experimental data. it is well known that fitting models to oscillating data, and even more so, reconstructing networks from such data, are extremely hard problems, since models tend to learn a steady state. in spite of this, oscillations were captured with high precision by our approach. effect of noise and dataset size we next studied the effect of different dataset sizes and different levels of noise in the data on the quality of network reconstruction. for this purpose, we added mean zero gaussian noise with standard deviations,,, and to the simulated concentration data, and furthermore subsampled equidistantly from the time series to generate data sets with t = and different time points for each of the five noise levels. then the network reconstruction was performed as described in the methods section. evaluation on experimental data: the dream challenge dataset to assess the performance of different reverse engineering approaches, stolovitzky et al. fathered the dream initiative. for this purpose, cantone et al. provided in-vivo data on a small, bio-engineered five-gene network, which was posted as challenge # within dream. this data was generated by inserting new promoter/gene combinations directly into the chromosomal dna of budding yeast. two time series of gene expression of the five inserted genes after stimulation were measured using quantitative pcr, measuring time points in minute intervals in time series and time points in minute intervals in time series measurements in both time series consist of negative log-ratios of the genes of interest to housekeeping genes, we therefore transformed the measured data to recover the original ratios. we used the first time series for network inference, i.e., t = the hyperparameters for the lq prior on the regulation strengths were set to q = and s = parameters for the gamma priors on the threshold parameters were manually set to concentrate probability mass near the mean concentration value for each gene individually. an evaluation of results using the mean of the markov chain for each parameter, as done for the simulated data, results in auc values that are equivalent to guessing. this might indicate that either the level of noise present in the experimental data is too high, or that the posterior distribution has multiple modes, with the mean being an inappropriate summary statistic. we therefore searched the points sampled for the maximum a-posteriori mode, and evaluate this mode further in the following. clearly, data from the additional modes are available and can be studied similarly. the regulation strength parameters for the maximum a-posteriori mode are shown in table the dynamic behavior and fit of the model prediction to the experimental data is depicted in figure obviously, the general dynamics of the original data is well represented, with a moderate amount of smoothing. auc values of the reconstructed network are for sensitivity vs. 1-specificity, and for precision to recall. reconstructed maximum-a-posteriori regulation strength parameters β for the dream challenge # data. true edges present in the reference topology are marked in bold. to assess the quality of this result, we next compared performance of our approach to the performance of other approaches submitted to the dream challenge. for this purpose, we computed performance measures as were used in the dream challenge for our inferred network, and show results in table auc values for this comparison were calculated as described in. since our method gives us a topology with positive and negative regulation strengths, we have transformed our results to be suitable for the evaluation method used in dream 2: • by skipping the sign and dividing by the largest learned regulation strength for the directed-unsigned challenge. • for the two directed-signed challenges we only took the regulation strengths with the appropriate sign and divide them by the highest absolute regulation strength. results of the dream challenge # data of our method compared to submitted best results from. the third column gives the auc values for our method when self-regulations are omitted. our method clearly outperforms all submitted methods in the directed-signed-inhibitory challenge; furthermore, when neglecting self-regulations, we also beat the best submitted method in the directed-unsigned challenge. our method outperforms all submitted results in the challenge directed-signed-inhibitory. one difficulty we observed was, that our approach learned many strong self-regulations of genes, possibly because of an improper balancing of the priors on synthesis/degradation rates and regulation strengths. since there are no self-regulations in the dream challenge # data, we provide an additional evaluation when disregarding self-regulations, results are shown in the third column of table in this case, we not only outperform all submitted approaches in the directed-signed-inhibitory challenge, but also beat the best models in the directed-unsigned challenge. discussion and  we have developed a novel methodological approach to reverse engineer gene regulatory networks from gene expression time series data, and evaluated this approach on both simulated and real gene expression data. the combination of ordinary differential equations and the bayes' regularized inference technique can be used for the quantitative analysis of complex cellular processes. non-linear differential equations are able to describe complex dynamic behavior, and a rich mathematical theory for analyzing them is well established. our method combines these advantages of differential equations with the advantages of a bayesian framework, which is able to capture noise in data, makes it possible to include biological knowledge into the learning process, and allows the computation of probability distributions over model topologies and model parameters. the latter is one of the main advantages of the mcmc approach. the information about distributions can be used to make predictions of future states of the network together with confidence intervals on the predictions made. this may allow it to take alternative future scenarios into account, and could be used to design additional most informative experiments that will help to distinguish between corresponding topologies or parameter sets. we therefore think that our approach will be highly useful to elucidate regulatory networks in an iterative procedure with several rounds of experiment, network inference, and experiment design. in contrast to the usual approach of minimizing an error between experimentally measured concentrations and model predictions, our likelihood function uses the difference between model slopes and experimental slopes. we furthermore integrate a smoothing spline approximation into the likelihood, automatically performing an optimized tradeoff between an accurate representation of the experimental data, and smoothing out noise. fitting of model parameters on slopes has the advantage that no numerical integration of the model is required in each step of the optimization or sampling. instead, we must estimate smoothing splines and slopes. this can be carried out much faster than numerical integration, enabling us to use a markov chain sampler on the posterior distribution instead of plain maximization. we have evaluated our approach on simulated and on real experimental data from a synthetic gene regulatory network. on the simulated data, we have shown that our approach can reconstruct the underlying topology with high accuracy. as expected, performance deteriorates with increasing levels of noise and with decreasing number of different time points available. we emphasize that the simulated example chosen is a difficult learning task due to the oscillations in the data. it is obvious that sufficient data points are required to sample the full dynamics of the oscillating network, and that oscillations quickly break down in the presence of noise. on the dream data, our method yields superior results when compared to other approaches that were submitted to dream in the directed-signed-inhibitory and directed-unsigned categories. importantly, our analysis shows that there are multiple posterior modes that describe the data well, which may explain a surprising result of the original dream challenge: as reported by stolovitzky et al., none of the submitted models were able to accurately reconstruct the original network topology from the synthetic data. our results indicate that this might be due to a dense population of local optima, in which network reconstruction approaches looking for a single optimal topology might get trapped and return suboptimal solutions. an obvious conclusion is that further experiments are required to resolve ambiguities in network reconstruction. this emphasizes the need for robust and efficient methods for optimum experimental design. our sampling approach may be a good starting point for such experiment design, since it analyzes the full distribution over model parameters, and therefore yields information on alternative network topologies and confidence intervals on parameters, which are instrumental to design experiments that elucidate the network topology further. authors' contributions lk designed the research and methodology. jm and dr wrote the software, analyzed the data, and drafted the paper. gr and lk advised on algorithm design and data analysis, and wrote the final version of the paper. all authors read and approved the final version of the paper. supplementary material additional file 1 supplementary material. mathematical proof. click here for file acknowledgements we would like to thank a. szabowski, h. busch, n. radde and r. eils for helpful comments and discussions. furthermore, we thank m. p. cosma and d. di bernardo for the provision of the challenge # data from dream we are grateful for funding to the german ministry of education and research, grant number.
we present a novel approach to infer such models from data, based on nonlinear differential equations, which we embed into a stochastic bayesian framework. our approach combines dynamic modeling using differential equations with a stochastic learning framework, thus bridging the gap between biophysical modeling and stochastic inference approaches. on the simulated data, we study effects of different levels of noise and dataset sizes. at the same time, quantitative, dynamic models, ideally with probability distributions over model topologies and parameters, are highly desirable. results on real data show that the dynamics and main regulatory interactions are correctly reconstructed. results show that the method can reap the advantages of both worlds, and allows the reconstruction of biophysically accurate dynamic models from noisy data. the reconstruction of gene regulatory networks from time series gene expression data is one of the most difficult problems in systems biology. this is due to several reasons, among them the combinatorial explosion of possible network topologies, limited information content of the experimental data with high levels of noise, and the complexity of gene regulation at the transcriptional, translational and post-translational levels. in addition, the stochastic learning framework used permits the computation of probability distributions over models and model parameters, which holds interesting prospects for experimental design purposes. using stochastic sampling from the bayes' posterior distribution, our approach can infer different likely network topologies and model parameters along with their respective probabilities from given data. furthermore, the bayesian framework allows it to easily integrate prior knowledge into the inference process. we thus address both the stochasticity of experimental data and the need for quantitative dynamic models.
 the reconstruction of gene regulatory networks from time series gene expression data is one of the most difficult problems in systems biology. this is due to several reasons, among them the combinatorial explosion of possible network topologies, limited information content of the experimental data with high levels of noise, and the complexity of gene regulation at the transcriptional, translational and post-translational levels. at the same time, quantitative, dynamic models, ideally with probability distributions over model topologies and parameters, are highly desirable.  we present a novel approach to infer such models from data, based on nonlinear differential equations, which we embed into a stochastic bayesian framework. we thus address both the stochasticity of experimental data and the need for quantitative dynamic models. furthermore, the bayesian framework allows it to easily integrate prior knowledge into the inference process. using stochastic sampling from the bayes' posterior distribution, our approach can infer different likely network topologies and model parameters along with their respective probabilities from given data. we evaluate our approach on simulated data and the challenge # data from the dream initiative. on the simulated data, we study effects of different levels of noise and dataset sizes. results on real data show that the dynamics and main regulatory interactions are correctly reconstructed.  our approach combines dynamic modeling using differential equations with a stochastic learning framework, thus bridging the gap between biophysical modeling and stochastic inference approaches. results show that the method can reap the advantages of both worlds, and allows the reconstruction of biophysically accurate dynamic models from noisy data. in addition, the stochastic learning framework used permits the computation of probability distributions over models and model parameters, which holds interesting prospects for experimental design purposes. 