Reconstructing nonlinear dynamic models of gene regulation using stochastic sampling

Abstract

Background: 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.

Results: 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 #3 data from the DREAM 2 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.

Conclusions: 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.

Figure 1

Regulation functions. Hill functions f_ij for different Hill coefficients m = 1, 3, 5. The left plot shows an activation, the right plot an inhibitory effect. The threshold θ_j was chosen equal to 3 for both plots, at this concentration of the regulating gene j, half the maximum effect on gene i is achieved.

Figure 2

L_q Prior. Plot of the two-dimensional L_q prior p(β₁, β₂):= L_q(β₁; q, s)·L_q(β₂; q, s) for q = 0.5 and s = 1. It can clearly be seen how this prior favors points (β₁, β₂) where one of the two components is approximately zero over points at the same distance from the origin with both β₁, β₂ ≠ 0.

Figure 3

Gold Standard Topology and Simulated Data. (a) True network of the DREAM 2 challenge #3 five gene time series data, showing the bio-engineered interactions between the five genes artificially inserted into yeast. (b) Time course of simulation with model in arbitrary time and concentration units, for the simulated five gene model. Different numbers of equidistant time points from this data were used for network reconstruction in the simulation study. The time courses of gene 2 and gene 3 are almost the same.

Figure 4

AUC Values for Simulated Data. AUC values for different noise levels and different numbers of time points used for network reconstruction. The standard deviation of the noise was varied from σ = 0 to σ = 0.3, the number of time points from T = 10 to T = 200. The plots show (a) AUC values under the ROC curve and (b) AUC values for PR curves for varying T and σ. The blue surface indicates the AUC_ROC and AUC_PR values that would follow for random guessing.

Figure 5

Learned Dynamics of DREAM 2 Challenge Data. Plot of the experimental data from the DREAM 2 challenge, in comparison to time courses simulated with reconstructed model parameters. Shown in black is the smoothed data.