BioDiVinE: A Framework for Parallel Analysis of Biological Models

In this paper a novel tool BioDiVinEfor parallel analysis of biological models is presented. The tool allows analysis of biological models specified in terms of a set of chemical reactions. Chemical reactions are transformed into a system of multi-affine differential equations. BioDiVinE employs techniques for finite discrete abstraction of the continuous state space. At that level, parallel analysis algorithms based on model checking are provided. In the paper, the key tool features are described and their application is demonstrated by means of a case study

Robustness Analysis of Stochastic Biochemical Systems

We propose a new framework for rigorous robustness analysis of stochastic biochemical systems that is based on probabilistic model checking techniques. We adapt the general definition of robustness introduced by Kitano to the class of stochastic systems modelled as continuous time Markov Chains in order to extensively analyse and compare robustness of biological models with uncertain parameters. The framework utilises novel computational methods that enable to effectively evaluate the robustness of models with respect to quantitative temporal properties and parameters such as reaction rate constants and initial conditions. We have applied the framework to gene regulation as an example of a central biological mechanism where intrinsic and extrinsic stochasticity plays crucial role due to low numbers of DNA and RNA molecules. Using our methods we have obtained a comprehensive and precise analysis of stochastic dynamics under parameter uncertainty. Furthermore, we apply our framework to compare several variants of two-component signalling networks from the perspective of robustness with respect to intrinsic noise caused by low populations of signalling components. We have successfully extended previous studies performed on deterministic models (ODE) and showed that stochasticity may significantly affect obtained predictions. Our case studies demonstrate that the framework can provide deeper insight into the role of key parameters in maintaining the system functionality and thus it significantly contributes to formal methods in computational systems biology

Robustness Analysis of Stochastic Biochemical Systems

<div><p>We propose a new framework for rigorous robustness analysis of stochastic biochemical systems that is based on probabilistic model checking techniques. We adapt the general definition of robustness introduced by Kitano to the class of stochastic systems modelled as continuous time Markov Chains in order to extensively analyse and compare robustness of biological models with uncertain parameters. The framework utilises novel computational methods that enable to effectively evaluate the robustness of models with respect to quantitative temporal properties and parameters such as reaction rate constants and initial conditions. We have applied the framework to gene regulation as an example of a central biological mechanism where intrinsic and extrinsic stochasticity plays crucial role due to low numbers of DNA and RNA molecules. Using our methods we have obtained a comprehensive and precise analysis of stochastic dynamics under parameter uncertainty. Furthermore, we apply our framework to compare several variants of two-component signalling networks from the perspective of robustness with respect to intrinsic noise caused by low populations of signalling components. We have successfully extended previous studies performed on deterministic models (ODE) and showed that stochasticity may significantly affect obtained predictions. Our case studies demonstrate that the framework can provide deeper insight into the role of key parameters in maintaining the system functionality and thus it significantly contributes to formal methods in computational systems biology.</p></div

Model of a two-component signalling pathway.

<p>(A) Basic topology of the two-component signalling pathway. (B) Modified topology of the two-component signalling pathway, additionally, histidine kinase <i>H</i> catalyses dephosporylation of the response regulator <i>R</i>. (C) Reactions specifying the biochemical model of the two considered topologies of the two-component signalling pathway. Phosphorylation of the first component <i>H</i> catalysed by the input signal <i>S</i> and phosphorylation of the second component <i>R</i> are shared by both topologies, the only difference is in the second component’s dephophorylation. Additionally, we consider unregulated proteosynthesis/degradation reactions for both topology variants. Reaction topology in (A) and (B) was created using CellDesigner <a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0094553#pone.0094553-Funahashi1" target="_blank">[54]</a>.</p

Robustness analysis workflow.

<p>The robustness analysis framework considers several objects on the input side. In particular, stochastic kinetic model is supplied with the quantitative hypothesis and the perturbation space of interest. The robustness analysis procedure systematically traverses the perturbation space and explores the system’s functionality determined by the quantitative hypothesis. The output side of the framework provides the evaluation function describing the system’s functionality with respect to the perturbation space. A single value characterising the system robustness is computed by integrating the evaluation function over the perturbation space.</p

High signal region in Model 2.

<p>A magnification of the high signal region in Model 2, where increasing levels of regulation by the sigmoid coefficient <i>n</i> leads to a paradoxical increase of output response noise instead of a decrease. Even though the inaccuracy is large we consider the trend to be strong and thus real.</p

Noise in populations or <i>H</i> and <i>R</i> in both models.

<p>Noise in <i>H</i> (A) and <i>R</i> (B) in both models has been computed with respect to perturbations of signal <i>S</i> over five selected intervals and for three distinct levels of inherent production noise represented by sigmoid coefficient . We can see that in all cases, with increasing regulation by <i>n</i>, the intrinsic noise in the dynamics of each of the signalling components decreases.</p

Results of robustness analysis for hypothesis (1) using the until operator.

<p>Hypothesis (1) requires stabilisation of in the low concentration mode (). A CSL formula with the until operator is used in this case. Each of the curves represents the evaluation function over degradation obtained for a particular setting of . More precisely, the horizontal axis shows the perturbation of <i>pRB</i> degradation rate and the vertical axis shows the probability of the hypothesis to be satisfied. In the upper left corner, robustness values are shown for each of the curves. The values are displayed with the absolute error quantifying the precision of the approximate method. For comparison, the values are computed also on piece-wise affine approximations of the evaluation function. It can be seen that the robustness values are small which is due to the fact that fluctuations of molecular numbers cause frequent crossing of the required bound in the considered time horizon.</p

Comparison of models by <i>Rp</i> noise robustness.

<p>Robustness <i>Rp</i> noise in both models has been computed with respect to perturbations of signal <i>S</i> over five selected intervals of the input signal and for three distinct levels of the intrinsic noise in signalling component dynamics represented by sigmoid coefficient . Perturbations were not computed over the whole interval due to high computational demands. From the computed values of individual refined subspaces as well as the aggregated robustness values for each input signal interval we can see that for lower values of signal <i>S</i> (up-to 10), Model 2 embodies lower output response noise than Model 1 (spontaneous dephosphorylation). While the output response noise in Model 1 tends to converge to values between 8 and 10, Model 2 exhibits a permanent (almost linear) increase in the output response noise over most of the studied portion of the perturbation space. A super-linear increase of the noise is observed for strong input signals. Another interesting aspect is that, with increasing levels of gene regulation given by sigmoid coefficient <i>n</i>, the overall noise in <i>Rp</i> decreases over the whole interval of signal values for Model 1 and most of the interval for Model 2. However, there is an anomaly in Model 2 in the high signal region [19.0, 20.0], where with decreasing noise in <i>R</i> and <i>H</i> (see <a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0094553#pone-0094553-g015" target="_blank">Figure 15</a>) the noise in <i>Rp</i> increases.</p

Results of robustness analysis for hypothesis (1) using the reward operator.

<p>Hypothesis (1) requires stabilisation of in the low concentration mode (). A CSL formula with cumulative reward operator is used in this case. Each of the curves represents the evaluation function over degradation obtained for a particular setting of . More precisely, the horizontal axis shows the perturbation of <i>pRB</i> degradation rate and the vertical axis shows the probability of the hypothesis to be satisfied. In the upper left corner, robustness values are shown for each of the curves. The values are displayed with the absolute error quantifying the precision of the approximate method. For comparison, the values are computed also on piece-wise affine approximations of the evaluation function. It can be seen that the robustness values change rapidly with different settings of . This observation goes with the fact that with faster degradation of there is a higher probability that the positively self-regulated protein is locked in the stable mode of no production. The decrease of the value with increasing is due to the weakening effect of inhibition by <i>pRB</i>.</p