Models Of Biochemical Reaction Networks Research Articles

A query-response causal analysis of reaction events in biochemical reaction networks

The stochastic kinetics of biochemical reaction networks is described by a chemical master equation (CME) and the underlying laws of mass action. Assuming network-free simulations of the rule-based models of biochemical reaction networks (BRNs), this paper departs from the usual analysis of network dynamics as the time-dependent distributions of chemical species counts, and instead considers statistically evaluating the sequences of reaction events generated from the stochastic simulations. The reaction event-time series can be used for reaction clustering, identifying rare events, and recognizing the periods of increased or steady-state activity. However, the main aim of this paper is to device an effective method for identifying causally and anti-causally related sub-sequences of reaction events using their empirical probabilities. This allows discovering some of the causal dynamics of BRNs as well as uncovering their short-term deterministic behaviors. In particular, it is proposed that the reaction sub-sequences that are conditionally nearly certain or nearly uncertain can be considered as being causally related. Moreover, since the time-ordering of reaction events is locally irrelevant, the reaction sub-sequences can be transformed into the reaction sets or multi-sets. The distance metrics can be then used to define the equivalences among the reaction events. The proposed method for identifying the causally related reaction sub-sequences has been implemented as a computationally efficient query-response mechanism. The method was evaluated for five models of genetic networks in seven defined numerical experiments. The models were simulated in BioNetGen using the open-source network-free simulator NFsim. This simulator had to be modified first to allow recording the traces of reaction events, and it is available in the Github repository, ploskot/nfsim_1.20. The generated event time-series were analyzed with Python and Matlab scripts. The whole process of data generation, analysis and visualization has been nearly fully automated using shell scripts. This demonstrates the opportunities for substantially increasing the research productivity by creating automated data generation and processing pipelines.

Efficient computation of adjoint sensitivities at steady-state in ODE models of biochemical reaction networks.

Dynamical models in the form of systems of ordinary differential equations have become a standard tool in systems biology. Many parameters of such models are usually unknown and have to be inferred from experimental data. Gradient-based optimization has proven to be effective for parameter estimation. However, computing gradients becomes increasingly costly for larger models, which are required for capturing the complex interactions of multiple biochemical pathways. Adjoint sensitivity analysis has been pivotal for working with such large models, but methods tailored for steady-state data are currently not available. We propose a new adjoint method for computing gradients, which is applicable if the experimental data include steady-state measurements. The method is based on a reformulation of the backward integration problem to a system of linear algebraic equations. The evaluation of the proposed method using real-world problems shows a speedup of total simulation time by a factor of up to 4.4. Our results demonstrate that the proposed approach can achieve a substantial improvement in computation time, in particular for large-scale models, where computational efficiency is critical.

The chemical Langevin equation for biochemical systems in dynamic environments.

Modeling and simulation of complex biochemical reaction networks form cornerstones of modern biophysics. Many of the approaches developed so far capture temporal fluctuations due to the inherent stochasticity of the biophysical processes, referred to as intrinsic noise. Stochastic fluctuations, however, predominantly stem from the interplay of the network with many other-and mostly unknown-fluctuating processes, as well as with various random signals arising from the extracellular world; these sources contribute extrinsic noise. Here, we provide a computational simulation method to probe the stochastic dynamics of biochemical systems subject to both intrinsic and extrinsic noise. We develop an extrinsic chemical Langevin equation (CLE)-a physically motivated extension of the CLE-to model intrinsically noisy reaction networks embedded in a stochastically fluctuating environment. The extrinsic CLE is a continuous approximation to the chemical master equation (CME) with time-varying propensities. In our approach, noise is incorporated at the level of the CME, and it can account for the full dynamics of the exogenous noise process, irrespective of timescales and their mismatches. We show that our method accurately captures the first two moments of the stationary probability density when compared with exact stochastic simulation methods while reducing the computational runtime by several orders of magnitude. Our approach provides a method that is practical, computationally efficient, and physically accurate to study systems that are simultaneously subject to a variety of noise sources.

SMGen: A Generator of Synthetic Models of Biochemical Reaction Networks

Several software tools for the simulation and analysis of biochemical reaction networks have been developed in the last decades; however, assessing and comparing their computational performance in executing the typical tasks of computational systems biology can be limited by the lack of a standardized benchmarking approach. To overcome these limitations, we propose here a novel tool, named SMGen, designed to automatically generate synthetic models of reaction networks that, by construction, are characterized by relevant features (e.g., system connectivity and reaction discreteness) and non-trivial emergent dynamics of real biochemical networks. The generation of synthetic models in SMGen is based on the definition of an undirected graph consisting of a single connected component that, generally, results in a computationally demanding task; to speed up the overall process, SMGen exploits a main–worker paradigm. SMGen is also provided with a user-friendly graphical user interface, which allows the user to easily set up all the parameters required to generate a set of synthetic models with any number of reactions and species. We analysed the computational performance of SMGen by generating batches of symmetric and asymmetric reaction-based models (RBMs) of increasing size, showing how a different number of reactions and/or species affects the generation time. Our results show that when the number of reactions is higher than the number of species, SMGen has to identify and correct a large number of errors during the creation process of the RBMs, a circumstance that increases the running time. Still, SMGen can generate synthetic models with hundreds of species and reactions in less than 7 s.

An automatic adaptive method to combine summary statistics in approximate Bayesian computation.

To infer the parameters of mechanistic models with intractable likelihoods, techniques such as approximate Bayesian computation (ABC) are increasingly being adopted. One of the main disadvantages of ABC in practical situations, however, is that parameter inference must generally rely on summary statistics of the data. This is particularly the case for problems involving high-dimensional data, such as biological imaging experiments. However, some summary statistics contain more information about parameters of interest than others, and it is not always clear how to weight their contributions within the ABC framework. We address this problem by developing an automatic, adaptive algorithm that chooses weights for each summary statistic. Our algorithm aims to maximize the distance between the prior and the approximate posterior by automatically adapting the weights within the ABC distance function. Computationally, we use a nearest neighbour estimator of the distance between distributions. We justify the algorithm theoretically based on properties of the nearest neighbour distance estimator. To demonstrate the effectiveness of our algorithm, we apply it to a variety of test problems, including several stochastic models of biochemical reaction networks, and a spatial model of diffusion, and compare our results with existing algorithms.

A Probabilistic Approach to Explore Signal Execution Mechanisms With Limited Experimental Data.

Mathematical models of biochemical reaction networks are central to the study of dynamic cellular processes and hypothesis generation that informs experimentation and validation. Unfortunately, model parameters are often not available and sparse experimental data leads to challenges in model calibration and parameter estimation. This can in turn lead to unreliable mechanistic interpretations of experimental data and the generation of poorly conceived hypotheses for experimental validation. To address this challenge, we evaluate whether a Bayesian-inspired probability-based approach, that relies on expected values for quantities of interest calculated from available information regarding the reaction network topology and parameters can be used to qualitatively explore hypothetical biochemical network execution mechanisms in the context of limited available data. We test our approach on a model of extrinsic apoptosis execution to identify preferred signal execution modes across varying conditions. Apoptosis signal processing can take place either through a mitochondria independent (Type I) mode or a mitochondria dependent (Type II) mode. We first show that in silico knockouts, represented by model subnetworks, successfully identify the most likely execution mode for specific concentrations of key molecular regulators. We then show that changes in molecular regulator concentrations alter the overall reaction flux through the network by shifting the primary route of signal flow between the direct caspase and mitochondrial pathways. Our work thus demonstrates that probabilistic approaches can be used to explore the qualitative dynamic behavior of model biochemical systems even with missing or sparse data.

Parameter estimation for biochemical reaction networks using Wasserstein distances

We present a method for estimating parameters in stochastic models of biochemical reaction networks by fitting steady-state distributions using Wasserstein distances. We simulate a reaction network at different parameter settings and train a Gaussian process to learn the Wasserstein distance between observations and the simulator output for all parameters. We then use Bayesian optimization to find parameters minimizing this distance based on the trained Gaussian process. The effectiveness of our method is demonstrated on the three-stage model of gene expression and a genetic feedback loop for which moment-based methods are known to perform poorly. Our method is applicable to any simulator model of stochastic reaction networks, including Brownian dynamics.

Efficient parameterization of large-scale dynamic models based on relative measurements.

Mechanistic models of biochemical reaction networks facilitate the quantitative understanding of biological processes and the integration of heterogeneous datasets. However, some biological processes require the consideration of comprehensive reaction networks and therefore large-scale models. Parameter estimation for such models poses great challenges, in particular when the data are on a relative scale. Here, we propose a novel hierarchical approach combining (i) the efficient analytic evaluation of optimal scaling, offset and error model parameters with (ii) the scalable evaluation of objective function gradients using adjoint sensitivity analysis. We evaluate the properties of the methods by parameterizing a pan-cancer ordinary differential equation model (>1000 state variables, >4000 parameters) using relative protein, phosphoprotein and viability measurements. The hierarchical formulation improves optimizer performance considerably. Furthermore, we show that this approach allows estimating error model parameters with negligible computational overhead when no experimental estimates are available, providing an unbiased way to weight heterogeneous data. Overall, our hierarchical formulation is applicable to a wide range of models, and allows for the efficient parameterization of large-scale models based on heterogeneous relative measurements. Supplementary code and data are available online at http://doi.org/10.5281/zenodo.3254429 and http://doi.org/10.5281/zenodo.3254441. Supplementary data are available at Bioinformatics online.

Simulation and inference algorithms for stochastic biochemical reaction networks: from basic concepts to state-of-the-art.

Stochasticity is a key characteristic of intracellular processes such as gene regulation and chemical signalling. Therefore, characterizing stochastic effects in biochemical systems is essential to understand the complex dynamics of living things. Mathematical idealizations of biochemically reacting systems must be able to capture stochastic phenomena. While robust theory exists to describe such stochastic models, the computational challenges in exploring these models can be a significant burden in practice since realistic models are analytically intractable. Determining the expected behaviour and variability of a stochastic biochemical reaction network requires many probabilistic simulations of its evolution. Using a biochemical reaction network model to assist in the interpretation of time-course data from a biological experiment is an even greater challenge due to the intractability of the likelihood function for determining observation probabilities. These computational challenges have been subjects of active research for over four decades. In this review, we present an accessible discussion of the major historical developments and state-of-the-art computational techniques relevant to simulation and inference problems for stochastic biochemical reaction network models. Detailed algorithms for particularly important methods are described and complemented with Matlab® implementations. As a result, this review provides a practical and accessible introduction to computational methods for stochastic models within the life sciences community.

Efficient computation of steady states in large-scale ODE models of biochemical reaction networks

In systems and computational biology, ordinary differential equations are used for the mechanistic modelling of biochemical networks. These models can easily have hundreds of states and parameters. Typically most parameters are unknown and estimated by fitting model output to observation. During parameter estimation the model needs to be solved repeatedly, sometimes millions of times. This can then be a computational bottleneck, and limits the employment of such models.In many situations the experimental data provides information about the steady state of the biochemical reaction network. In such cases one only needs to obtain the equilibrium state for a given set of model parameters. In this paper we exploit this fact and solve the steady state problem directly rather than integrating the ODE forward in time until steady state is reached. We use Newton’s method - like some previous studies - and develop several improvements to achieve robust convergence. To address the reliance of Newtons method on good initial guesses, we propose a continuation method. We show that the method works robustly in this setting and achieves a speed up of up to 100 compared to using ODE solves.

Optimal experimental design for machine learning using the Fisher information matrix

Optimal experimental design (OED) refers to a class of methods for selecting new data collection conditions that minimize the statistical uncertainty in the inferred parameter values of a model. The Fisher information matrix (FIM) gives an estimate of the relative uncertainty in and correlation among the model parameters based on the local curvature of the cost function. FIM-based approaches to OED allow for rapid assessment of many different experimental conditions (e.g., input data type, parameterizations, etc.). In machine learning models, accurate parameter estimates are often not a priority (nor even desirable) as they have no direct physical meaning. Instead, one would like to minimize the uncertainty in the model predictions for several quantities of interest. FIM approaches to OED can be generalized to minimize statistical variance, not in parameters, but in predictions of the quantities of interest. This approach has been applied, for example, to systems biology models of biochemical reaction networks [Transtrum and Qiu, BMC Bioinformatics 13(1), 181 (2012)]. Preliminary application of the FIM to optimize experimental design for source localization in an uncertain ocean environment is a first step towards an efficient machine learning algorithm that produces results with the least uncertainty in the quantities of interest.

Bayesian parameter estimation for biochemical reaction networks using region-based adaptive parallel tempering.

MotivationMathematical models have become standard tools for the investigation of cellular processes and the unraveling of signal processing mechanisms. The parameters of these models are usually derived from the available data using optimization and sampling methods. However, the efficiency of these methods is limited by the properties of the mathematical model, e.g. non-identifiabilities, and the resulting posterior distribution. In particular, multi-modal distributions with long valleys or pronounced tails are difficult to optimize and sample. Thus, the developement or improvement of optimization and sampling methods is subject to ongoing research.ResultsWe suggest a region-based adaptive parallel tempering algorithm which adapts to the problem-specific posterior distributions, i.e. modes and valleys. The algorithm combines several established algorithms to overcome their individual shortcomings and to improve sampling efficiency. We assessed its properties for established benchmark problems and two ordinary differential equation models of biochemical reaction networks. The proposed algorithm outperformed state-of-the-art methods in terms of calculation efficiency and mixing. Since the algorithm does not rely on a specific problem structure, but adapts to the posterior distribution, it is suitable for a variety of model classes.Availability and implementationThe code is available both as Supplementary Material and in a Git repository written in MATLAB.Supplementary information Supplementary data are available at Bioinformatics online.

Inferring a nonlinear biochemical network model from a heterogeneous single-cell time course data

Mathematical modeling and analysis of biochemical reaction networks are key routines in computational systems biology and biophysics; however, it remains difficult to choose the most valid model. Here, we propose a computational framework for data-driven and systematic inference of a nonlinear biochemical network model. The framework is based on the expectation-maximization algorithm combined with particle smoother and sparse regularization techniques. In this method, a “redundant” model consisting of an excessive number of nodes and regulatory paths is iteratively updated by eliminating unnecessary paths, resulting in an inference of the most likely model. Using artificial single-cell time-course data showing heterogeneous oscillatory behaviors, we demonstrated that this algorithm successfully inferred the true network without any prior knowledge of network topology or parameter values. Furthermore, we showed that both the regulatory paths among nodes and the optimal number of nodes in the network could be systematically determined. The method presented in this study provides a general framework for inferring a nonlinear biochemical network model from heterogeneous single-cell time-course data.

Parameter Estimation by Picard-Iteration for Biochemical Networks with Noisy Data

Mathematical models of biochemical reaction networks are widespread used to shade light and understand biological process. The validation and parametert estimation of such models is an important task to ensure the quality of the model. Many parameter estimation strategies for biochemical reaction networks exist by now. This work focuses on the use of an Picard iteration based parameter estimation approach for ODE-based dynamical models, which inherently exploits the model structure. In the frame of this work the influence of measurement errors and noise on the resulting parameter estimation is examined.

Methods of Model Reduction for Large-Scale Biological Systems: A Survey of Current Methods and Trends

Complex models of biochemical reaction systems have become increasingly common in the systems biology literature. The complexity of such models can present a number of obstacles for their practical use, often making problems difficult to intuit or computationally intractable. Methods of model reduction can be employed to alleviate the issue of complexity by seeking to eliminate those portions of a reaction network that have little or no effect upon the outcomes of interest, hence yielding simplified systems that retain an accurate predictive capacity. This review paper seeks to provide a brief overview of a range of such methods and their application in the context of biochemical reaction network models. To achieve this, we provide a brief mathematical account of the main methods including timescale exploitation approaches, reduction via sensitivity analysis, optimisation methods, lumping, and singular value decomposition-based approaches. Methods are reviewed in the context of large-scale systems biology type models, and future areas of research are briefly discussed.

Sensitivity estimation for stochastic models of biochemical reaction networks in the presence of extrinsic variability.

Determining the sensitivity of certain system states or outputs to variations in parameters facilitates our understanding of the inner working of that system and is an essential design tool for the de novo construction of robust systems. In cell biology, the output of interest is often the response of a certain reaction network to some input (e.g., stressors or nutrients) and one aims to quantify the sensitivity of this response in the presence of parameter heterogeneity. We argue that for such applications, parametric sensitivities in their standard form do not paint a complete picture of a system's robustness since one assumes that all cells in the population have the same parameters and are perturbed in the same way. Here, we consider stochasticreaction networks in which the parameters are randomly distributed over the population and propose a new sensitivity index that captures the robustness of system outputs upon changes in the characteristics of the parameter distribution, rather than the parameters themselves. Subsequently, we make use of Girsanov's likelihood ratio method to construct a Monte Carlo estimator of this sensitivity index. However, it turns out that this estimator has an exceedingly large variance. To overcome this problem, we propose a novel estimation algorithm that makes use of a marginalization of the path distribution of stochasticreaction networks and leads to Rao-Blackwellized estimators with reduced variance.

A combined model reduction algorithm for controlled biochemical systems

BackgroundSystems Biology continues to produce increasingly large models of complex biochemical reaction networks. In applications requiring, for example, parameter estimation, the use of agent-based modelling approaches, or real-time simulation, this growing model complexity can present a significant hurdle. Often, however, not all portions of a model are of equal interest in a given setting. In such situations methods of model reduction offer one possible approach for addressing the issue of complexity by seeking to eliminate those portions of a pathway that can be shown to have the least effect upon the properties of interest.MethodsIn this paper a model reduction algorithm bringing together the complementary aspects of proper lumping and empirical balanced truncation is presented. Additional contributions include the development of a criterion for the selection of state-variable elimination via conservation analysis and use of an ‘averaged’ lumping inverse. This combined algorithm is highly automatable and of particular applicability in the context of ‘controlled’ biochemical networks.ResultsThe algorithm is demonstrated here via application to two examples; an 11 dimensional model of bacterial chemotaxis in Escherichia coli and a 99 dimensional model of extracellular regulatory kinase activation (ERK) mediated via the epidermal growth factor (EGF) and nerve growth factor (NGF) receptor pathways. In the case of the chemotaxis model the algorithm was able to reduce the model to 2 state-variables producing a maximal relative error between the dynamics of the original and reduced models of only 2.8% whilst yielding a 26 fold speed up in simulation time. For the ERK activation model the algorithm was able to reduce the system to 7 state-variables, incurring a maximal relative error of 4.8%, and producing an approximately 10 fold speed up in the rate of simulation. Indices of controllability and observability are additionally developed and demonstrated throughout the paper. These provide insight into the relative importance of individual reactants in mediating a biochemical system’s input-output response even for highly complex networks.ConclusionsThrough application, this paper demonstrates that combined model reduction methods can produce a significant simplification of complex Systems Biology models whilst retaining a high degree of predictive accuracy. In particular, it is shown that by combining the methods of proper lumping and empirical balanced truncation it is often possible to produce more accurate reductions than can be obtained by the use of either method in isolation.

Scalable Parameter Estimation for Genome-Scale Biochemical Reaction Networks.

Mechanistic mathematical modeling of biochemical reaction networks using ordinary differential equation (ODE) models has improved our understanding of small- and medium-scale biological processes. While the same should in principle hold for large- and genome-scale processes, the computational methods for the analysis of ODE models which describe hundreds or thousands of biochemical species and reactions are missing so far. While individual simulations are feasible, the inference of the model parameters from experimental data is computationally too intensive. In this manuscript, we evaluate adjoint sensitivity analysis for parameter estimation in large scale biochemical reaction networks. We present the approach for time-discrete measurement and compare it to state-of-the-art methods used in systems and computational biology. Our comparison reveals a significantly improved computational efficiency and a superior scalability of adjoint sensitivity analysis. The computational complexity is effectively independent of the number of parameters, enabling the analysis of large- and genome-scale models. Our study of a comprehensive kinetic model of ErbB signaling shows that parameter estimation using adjoint sensitivity analysis requires a fraction of the computation time of established methods. The proposed method will facilitate mechanistic modeling of genome-scale cellular processes, as required in the age of omics.

An efficient finite-difference strategy for sensitivity analysis of stochastic models of biochemical systems

Sensitivity analysis characterizes the dependence of a model's behaviour on system parameters. It is a critical tool in the formulation, characterization, and verification of models of biochemical reaction networks, for which confident estimates of parameter values are often lacking. In this paper, we propose a novel method for sensitivity analysis of discrete stochastic models of biochemical reaction systems whose dynamics occur over a range of timescales. This method combines finite-difference approximations and adaptive tau-leaping strategies to efficiently estimate parametric sensitivities for stiff stochastic biochemical kinetics models, with negligible loss in accuracy compared with previously published approaches. We analyze several models of interest to illustrate the advantages of our method.

Regularization Techniques to Overcome Overparameterization of Complex Biochemical Reaction Networks.

Models of biochemical reaction networks commonly contain a large number of parameters while at the same time there is only a limited amount of (noisy) data available for their estimation. As such, the values of many parameters are not well known as nominal parameter values have to be determined from the open scientific literature and a significant number of the values may have been derived in different cell types or organisms than that which is modeled. There clearly is a need to estimate at least some of the parameter values from experimental data, however, the small amount of available data and the large number of parameters commonly found in these types of models, require the use of regularization techniques to avoid over fitting. A tutorial of regularization techniques, including parameter set selection, precedes a case study of estimating parameters in a signal transduction network. Cross validation rather than fitting results are presented to further emphasize the need for models that generalize well to new data instead of simply fitting the current data.