A Provably Convergent Control Closure Scheme for the Method of Moments of the Chemical Master Equation

In this article, we introduce a novel moment closure scheme based on concepts from model predictive control (MPC) to accurately describe the time evolution of the statistical moments of the solution of the chemical master equation (CME). The method of moments, a set of ordinary differential equations frequently used to calculate the first nm moments, is generally not closed since lower-order moments depend on higher-order moments. To overcome this limitation, we interpret the moment equations as a nonlinear dynamical system, where the first nm moments serve as states, and the closing moments serve as the control input. We demonstrate the efficacy of our approach using three example systems and show that it outperforms existing closure schemes. For polynomial systems, which encompass all mass-action systems, we provide probability bounds for the error between true and estimated moment trajectories. We achieve this by combining the convergence properties of a priori moment estimates from stochastic simulations with guarantees for nonlinear reference tracking MPC. Our proposed method offers an effective solution to accurately predict the time evolution of moments of the CME, which has wide-ranging implications for many fields, including biology, chemistry, and engineering.

The impossible challenge of estimating non-existent moments of the Chemical Master Equation

MOTIVATION

The Chemical Master Equation (CME) is a set of linear differential equations that describes the evolution of the probability distribution on all possible configurations of a (bio-)chemical reaction system. Since the number of configurations and therefore the dimension of the CME rapidly increases with the number of molecules, its applicability is restricted to small systems. A widely applied remedy for this challenge is moment-based approaches which consider the evolution of the first few moments of the distribution as summary statistics for the complete distribution. Here, we investigate the performance of two moment-estimation methods for reaction systems whose equilibrium distributions encounter fat-tailedness and do not possess statistical moments.

RESULTS

We show that estimation via stochastic simulation algorithm (SSA) trajectories lose consistency over time and estimated moment values span a wide range of values even for large sample sizes. In comparison, the method of moments returns smooth moment estimates but is not able to indicate the non-existence of the allegedly predicted moments. We furthermore analyze the negative effect of a CME solution's fat-tailedness on SSA run times and explain inherent difficulties. While moment-estimation techniques are a commonly applied tool in the simulation of (bio-)chemical reaction networks, we conclude that they should be used with care, as neither the system definition nor the moment-estimation techniques themselves reliably indicate the potential fat-tailedness of the CME's solution.

Periportal steatosis in mice affects distinct parameters of pericentral drug metabolism

Little is known about the impact of morphological disorders in distinct zones on metabolic zonation. It was described recently that periportal fibrosis did affect the expression of CYP proteins, a set of pericentrally located drug-metabolizing enzymes. Here, we investigated whether periportal steatosis might have a similar effect. Periportal steatosis was induced in C57BL6/J mice by feeding a high-fat diet with low methionine/choline content for either two or four weeks. Steatosis severity was quantified using image analysis. Triglycerides and CYP activity were quantified in photometric or fluorometric assay. The distribution of CYP3A4, CYP1A2, CYP2D6, and CYP2E1 was visualized by immunohistochemistry. Pharmacokinetic parameters of test drugs were determined after injecting a drug cocktail (caffeine, codeine, and midazolam). The dietary model resulted in moderate to severe mixed steatosis confined to periportal and midzonal areas. Periportal steatosis did not affect the zonal distribution of CYP expression but the activity of selected CYPs was associated with steatosis severity. Caffeine elimination was accelerated by microvesicular steatosis, whereas midazolam elimination was delayed in macrovesicular steatosis. In summary, periportal steatosis affected parameters of pericentrally located drug metabolism. This observation calls for further investigations of the highly complex interrelationship between steatosis and drug metabolism and underlying signaling mechanisms.

Quasi-Entropy Closure: a fast and reliable approach to close the moment equations of the Chemical Master Equation

MOTIVATION

The Chemical Master Equation is a stochastic approach to describe the evolution of a (bio)chemical reaction system. Its solution is a time-dependent probability distribution on all possible configurations of the system. As this number is typically large, the Master Equation is often practically unsolvable. The Method of Moments reduces the system to the evolution of a few moments, which are described by ordinary differential equations. Those equations are not closed, since lower order moments generally depend on higher order moments. Various closure schemes have been suggested to solve this problem. Two major problems with these approaches are first that they are open loop systems, which can diverge from the true solution, and second, some of them are computationally expensive.

RESULTS

Here we introduce Quasi-Entropy Closure, a moment-closure scheme for the Method of Moments. It estimates higher order moments by reconstructing the distribution that minimizes the distance to a uniform distribution subject to lower order moment constraints. Quasi-Entropy Closure can be regarded as an advancement of Zero-Information Closure, which similarly maximizes the information entropy. Results show that both approaches outperform truncation schemes. Quasi-Entropy Closure is computationally much faster than Zero-Information Closure, although both methods consider solutions on the space of configurations and hence do not completely overcome the curse of dimensionality. In addition, our scheme includes a plausibility check for the existence of a distribution satisfying a given set of moments on the feasible set of configurations. All results are evaluated on different benchmark problems.

SUPPLEMENTARY INFORMATION

Supplementary data are available at Bioinformatics online.

ROSIE: RObust Sparse ensemble for outlIEr detection and gene selection in cancer omics data

The extraction of novel information from omics data is a challenging task, in particular, since the number of features (e.g. genes) often far exceeds the number of samples. In such a setting, conventional parameter estimation leads to ill-posed optimization problems, and regularization may be required. In addition, outliers can largely impact classification accuracy.

Here we introduce ROSIE, an ensemble classification approach, which combines three sparse and robust classification methods for outlier detection and feature selection and further performs a bootstrap-based validity check. Outliers of ROSIE are determined by the rank product test using outlier rankings of all three methods, and important features are selected as features commonly selected by all methods.

We apply ROSIE to RNA-Seq data from The Cancer Genome Atlas (TCGA) to classify observations into Triple-Negative Breast Cancer (TNBC) and non-TNBC tissue samples. The pre-processed dataset consists of 16,600 genes and more than 1,000 samples. We demonstrate that ROSIE selects important features and outliers in a robust way. Identified outliers are concordant with the distribution of the commonly selected genes by the three methods, and results are in line with other independent studies. Furthermore, we discuss the association of some of the selected genes with the TNBC subtype in other investigations. In summary, ROSIE constitutes a robust and sparse procedure to identify outliers and important genes through binary classification. Our approach is ad hoc applicable to other datasets, fulfilling the overall goal of simultaneously identifying outliers and candidate disease biomarkers to the targeted in therapy research and personalized medicine frameworks.

Hepatectomy-Induced Alterations in Hepatic Perfusion and Function - Toward Multi-Scale Computational Modeling for a Better Prediction of Post-hepatectomy Liver Function

Liver resection causes marked perfusion alterations in the liver remnant both on the organ scale (vascular anatomy) and on the microscale (sinusoidal blood flow on tissue level). These changes in perfusion affect hepatic functions via direct alterations in blood supply and drainage, followed by indirect changes of biomechanical tissue properties and cellular function. Changes in blood flow impose compression, tension and shear forces on the liver tissue. These forces are perceived by mechanosensors on parenchymal and non-parenchymal cells of the liver and regulate cell-cell and cell-matrix interactions as well as cellular signaling and metabolism. These interactions are key players in tissue growth and remodeling, a prerequisite to restore tissue function after PHx. Their dysregulation is associated with metabolic impairment of the liver eventually leading to liver failure, a serious post-hepatectomy complication with high morbidity and mortality. Though certain links are known, the overall functional change after liver surgery is not understood due to complex feedback loops, non-linearities, spatial heterogeneities and different time-scales of events. Computational modeling is a unique approach to gain a better understanding of complex biomedical systems. This approach allows (i) integration of heterogeneous data and knowledge on multiple scales into a consistent view of how perfusion is related to hepatic function; (ii) testing and generating hypotheses based on predictive models, which must be validated experimentally and clinically. In the long term, computational modeling will (iii) support surgical planning by predicting surgery-induced perfusion perturbations and their functional (metabolic) consequences; and thereby (iv) allow minimizing surgical risks for the individual patient. Here, we review the alterations of hepatic perfusion, biomechanical properties and function associated with hepatectomy. Specifically, we provide an overview over the clinical problem, preoperative diagnostics, functional imaging approaches, experimental approaches in animal models, mechanoperception in the liver and impact on cellular metabolism, omics approaches with a focus on transcriptomics, data integration and uncertainty analysis, and computational modeling on multiple scales. Finally, we provide a perspective on how multi-scale computational models, which couple perfusion changes to hepatic function, could become part of clinical workflows to predict and optimize patient outcome after complex liver surgery.

The role of stochastic sequestration dynamics for intrinsic noise filtering in signaling network motifs

The relation between design principles of signaling network motifs and their robustness against intrinsic noise still remains illusive. In this work we investigate the role of cascading for coping with intrinsic noise due to stochasticity in molecular reactions. We use stochastic approaches to quantify fluctuations in the terminal kinase of phosphorylation-dephosphorylation cascade motifs and demonstrate that cascading highly affects these fluctuations. We show that this purely stochastic effect can be explained by time-varying sequestration of upstream kinase molecules. In particular, we discuss conditions on time scales and parameter regimes which lead to a reduction of output fluctuations. Our results are put into biological context by adapting rate parameters of our modeling approach to biologically feasible ranges for general binding-unbinding and phosphorylation-dephosphorylation mechanisms. Overall, this study reveals a novel role of stochastic sequestration for dynamic noise filtering in signaling cascade motifs.

The circuit-breaking algorithm for monotone systems

In earlier work, we have introduced the circuit-breaking algorithm (CBA) for the analysis of intracellular regulation networks. This algorithm uses the network topology to construct a one-dimensional circuit-characteristic whose zeros correspond to the fixed points of the system. In this study, we apply the CBA to monotone systems whose flow preserves a partial order with respect to some orthant cone. We consider relations between stability of fixed points and the derivative of the corresponding zeros of the circuit-characteristic. In particular, we derive sufficient conditions for instability in case of global asymptotic stability of the open-loop system. Furthermore, we fully characterize stability of the fixed points if in addition the system is monotone. Combined with the theory of monotone systems, our results are used to characterize the long-term behavior of two models for different intracellular regulation processes.

Hamiltonian Monte Carlo methods for efficient parameter estimation in steady state dynamical systems

Background

Parameter estimation for differential equation models of intracellular processes is a highly relevant bu challenging task. The available experimental data do not usually contain enough information to identify all parameters uniquely, resulting in ill-posed estimation problems with often highly correlated parameters. Sampling-based Bayesian statistical approaches are appropriate for tackling this problem. The samples are typically generated via Markov chain Monte Carlo, however such methods are computationally expensive and their convergence may be slow, especially if there are strong correlations between parameters. Monte Carlo methods based on Euclidean or Riemannian Hamiltonian dynamics have been shown to outperform other samplers by making proposal moves that take the local sensitivities of the system’s states into account and accepting these moves with high probability. However, the high computational cost involved with calculating the Hamiltonian trajectories prevents their widespread use for all but the smallest differential equation models. The further development of efficient sampling algorithms is therefore an important step towards improving the statistical analysis of predictive models of intracellular processes.

Results

We show how state of the art Hamiltonian Monte Carlo methods may be significantly improved for steady state dynamical models. We present a novel approach for efficiently calculating the required geometric quantities by tracking steady states across the Hamiltonian trajectories using a Newton-Raphson method and employing local sensitivity information. Using our approach, we compare both Euclidean and Riemannian versions of Hamiltonian Monte Carlo on three models for intracellular processes with real data and demonstrate at least an order of magnitude improvement in the effective sampling speed. We further demonstrate the wider applicability of our approach to other gradient based MCMC methods, such as those based on Langevin diffusions.

Conclusion

Our approach is strictly benefitial in all test cases. The Matlab sources implementing our MCMC methodology is available from https://github.com/a-kramer/ode_rmhmc.

Electronic supplementary material

The online version of this article (doi:10.1186/1471-2105-15-253) contains supplementary material, which is available to authorized users.

MCMC_CLIB-an advanced MCMC sampling package for ODE models

SUMMARY

We present a new C implementation of an advanced Markov chain Monte Carlo (MCMC) method for the sampling of ordinary differential equation (ode) model parameters. The software mcmc_clib uses the simplified manifold Metropolis-adjusted Langevin algorithm (SMMALA), which is locally adaptive; it uses the parameter manifold's geometry (the Fisher information) to make efficient moves. This adaptation does not diminish with MC length, which is highly advantageous compared with adaptive Metropolis techniques when the parameters have large correlations and/or posteriors substantially differ from multivariate Gaussians. The software is standalone (not a toolbox), though dependencies include the GNU scientific library and sundials libraries for ode integration and sensitivity analysis.

AVAILABILITY AND IMPLEMENTATION

The source code and binary files are freely available for download at http://a-kramer.github.io/mcmc_clib/. This also includes example files and data. A detailed documentation, an example model and user manual are provided with the software.

CONTACT

andrei.kramer@ist.uni-stuttgart.de.

iVUN: interactive Visualization of Uncertain biochemical reaction Networks

BACKGROUND

Mathematical models are nowadays widely used to describe biochemical reaction networks. One of the main reasons for this is that models facilitate the integration of a multitude of different data and data types using parameter estimation. Thereby, models allow for a holistic understanding of biological processes. However, due to measurement noise and the limited amount of data, uncertainties in the model parameters should be considered when conclusions are drawn from estimated model attributes, such as reaction fluxes or transient dynamics of biological species.

METHODS AND RESULTS

We developed the visual analytics system iVUN that supports uncertainty-aware analysis of static and dynamic attributes of biochemical reaction networks modeled by ordinary differential equations. The multivariate graph of the network is visualized as a node-link diagram, and statistics of the attributes are mapped to the color of nodes and links of the graph. In addition, the graph view is linked with several views, such as line plots, scatter plots, and correlation matrices, to support locating uncertainties and the analysis of their time dependencies. As demonstration, we use iVUN to quantitatively analyze the dynamics of a model for Epo-induced JAK2/STAT5 signaling.

CONCLUSION

Our case study showed that iVUN can be used to perform an in-depth study of biochemical reaction networks, including attribute uncertainties, correlations between these attributes and their uncertainties as well as the attribute dynamics. In particular, the linking of different visualization options turned out to be highly beneficial for the complex analysis tasks that come with the biological systems as presented here.

Modeling sphingomyelin synthase 1 driven reaction at the Golgi apparatus can explain data by inclusion of a positive feedback mechanism

Here we present a minimal mathematical model for the sphingomyelin synthase 1 (SMS1) driven conversion of ceramide to sphingomyelin based on chemical reaction kinetics. We demonstrate via mathematical analysis that this model is not able to qualitatively reproduce experimental measurements on lipid compositions after altering SMS1 activity. We prove that a positive feedback mechanism from the products to the reactants of the reaction is one possible model extension to explain these specific experimental data. The proposed mechanism in fact exists in vivo via protein kinase D and the ceramide transfer protein CERT. The model is further evaluated by additional observations from the literature.

Trajectory-oriented Bayesian experiment design versus Fisher A-optimal design: an in depth comparison study

Motivation: Experiment design strategies for biomedical models with the purpose of parameter estimation or model discrimination are in the focus of intense research. Experimental limitations such as sparse and noisy data result in unidentifiable parameters and render-related design tasks challenging problems. Often, the temporal resolution of data is a limiting factor and the amount of possible experimental interventions is finite. To address this issue, we propose a Bayesian experiment design algorithm to minimize the prediction uncertainty for a given set of experiments and compare it to traditional A-optimal design.

Results: In an in depth numerical study involving an ordinary differential equation model of the trans-Golgi network with 12 partly non-identifiable parameters, we minimized the prediction uncertainty efficiently for predefined scenarios. The introduced method results in twice the prediction precision as the same amount of A-optimal designed experiments while introducing a useful stopping criterion. The simulation intensity of the algorithm's major design step is thereby reasonably affordable. Besides smaller variances in the predicted trajectories compared with Fisher design, we could also achieve smaller parameter posterior distribution entropies, rendering this method superior to A-optimal Fisher design also in the parameter space.

Availability: Necessary software/toolbox information are available in the supplementary material. The project script including example data can be downloaded from http://www.ist.uni-stuttgart.de/%7eweber/BayesFisher2012.

Contact:patrick.weber@ist.uni-stuttgart.de

Supplementary Information:Supplementary data are available at Bioinformatics online.

Analyzing fixed points of intracellular regulation networks with interrelated feedback topology

Background

Modeling the dynamics of intracellular regulation networks by systems of ordinary differential equations has become a standard method in systems biology, and it has been shown that the behavior of these networks is often tightly connected to the network topology. We have recently introduced the circuit-breaking algorithm, a method that uses the network topology to construct a one-dimensional circuit-characteristic of the system. It was shown that this characteristic can be used for an efficient calculation of the system’s fixed points.

Results

Here we extend previous work and show several connections between the circuit-characteristic and the stability of fixed points. In particular, we derive a sufficient condition on the characteristic for a fixed point to be unstable for certain graph structures and demonstrate that the characteristic does not contain the information to decide whether a fixed point is asymptotically stable. All statements are illustrated on biological network models.

Conclusions

Single feedback circuits and their role for complex dynamic behavior of biological networks have extensively been investigated, but a transfer of most of these concepts to more complex topologies is difficult. In this context, our algorithm is a powerful new approach for the analysis of regulation networks that goes beyond single isolated feedback circuits.

Fixed point characterization of biological networks with complex graph topology

MOTIVATION

Feedback circuits are important motifs in biological networks and part of virtually all regulation processes that are needed for a reliable functioning of the cell. Mathematically, feedback is connected to complex behavior of the systems, which is often related to bifurcations of fixed points. Therefore, several approaches for the investigation of fixed points in biological networks have been developed in recent years. Many of them assume the fixed point coordinates to be known, and an efficient way to calculate the entire set of fixed points for interrelated feedback structures is highly desirable.

RESULTS

In this article, we consider regulatory network models, which are differential equations with an underlying directed graph that illustrates independencies among variables. We introduce the circuit-breaking algorithm (CBA), a method that constructs one-dimensional characteristics for these network models, which inherit important information about the system. In particular, fixed points are related to the zeros of these characteristics. The CBA operates on the graph topology, and results from graph theory are used in order to make calculations efficient. Our framework provides a general scheme for analyzing network models in terms of interrelated feedback circuits. The efficiency of the approach is demonstrated on a model for calcium oscillations based on experiments in hepatocytes, which consists of several interrelated feedback circuits.

Modeling feedback loops in the H-NS-mediated regulation of the Escherichia coli bgl operon

The histone-like nucleoid-associated protein H-NS is a global transcriptional repressor that controls approximately 5% of all genes in Escherichia coli and other enterobacteria. H-NS binds to DNA with low specificity. Nonetheless, repression of some loci is exceptionally specific. Experimental data for the E. coli bgl operon suggest that highly specific repression is caused by regulatory feedback loops. To analyze whether such feedback loops can account for the observed specificity of repression, here a model was built based on expression data. The model includes several regulatory interactions, which are synergy of repression by binding of H-NS to two regulatory elements, an inverse correlation of the rate of repression by H-NS and transcription, and a threshold for positive regulation by anti-terminator BglG, which is encoded within the operon. The latter two regulatory interactions represent feedback loops in the model. The resulting system of equations was solved for the expression level of the operon and analyzed with respect to different promoter activities. This analysis demonstrates that a small (3-fold) increase of the bgl promoter activity results in a strong (80-fold) enhancement of bgl operon expression. Thus, the parameters included into the model are sufficient to simulate specific repression by H-NS.

Systematic component selection for gene-network refinement

MOTIVATION

A quantitative description of interactions between cell components is a major challenge in Computational Biology. As a method of choice, differential equations are used for this purpose, because they provide a detailed insight into the dynamic behavior of the system. In most cases, the number of time points of experimental time series is usually too small to estimate the parameters of a model of a whole gene regulatory network based on differential equations, such that one needs to focus on subnetworks consisting of only a few components. For most approaches, the set of components of the subsystem is given in advance and only the structure has to be estimated. However, the set of components that influence the system significantly are not always known in advance, making a method desirable that determines both, the components that are included into the model and the parameters.

RESULTS

We have developed a method that uses gene expression data as well as interaction data between cell components to define a set of genes that we use for our modeling. In a subsequent step, we estimate the parameters of our model of piecewise linear differential equations and evaluate the results simulating the behavior of the system with our model. We have applied our method to the DNA repair system of Mycobacterium tuberculosis. Our analysis predicts that the gene Rv2719c plays an important role in this system.