Global parameter estimation methods for stochastic biochemical systems

Background  

The importance of stochasticity in cellular processes having low number of molecules has resulted in the development of stochastic models such as chemical master equation. As in other modelling frameworks, the accompanying rate constants are important for the end-applications like analyzing system properties (e.g. robustness) or predicting the effects of genetic perturbations. Prior knowledge of kinetic constants is usually limited and the model identification routine typically includes parameter estimation from experimental data. Although the subject of parameter estimation is well-established for deterministic models, it is not yet routine for the chemical master equation. In addition, recent advances in measurement technology have made the quantification of genetic substrates possible to single molecular levels. Thus, the purpose of this work is to develop practical and effective methods for estimating kinetic model parameters in the chemical master equation and other stochastic models from single cell and cell population experimental data.     

Accelerated maximum likelihood parameter estimation for stochastic biochemical systems

ABSTRACT: BACKGROUND: A prerequisite for the mechanistic simulation of a biochemical system is detailed knowledge of its kinetic parameters. Despite recent experimental advances, the estimation of unknown parameter values from observed data is still a bottleneck for obtaining accurate simulation results. Many methods exist for parameter estimation in deterministic biochemical systems; methods for discrete stochastic systems are less well developed. Given the probabilistic nature of stochastic biochemical models, a natural approach is to choose parameter values that maximize the probability of the observed data with respect to the unknown parameters, a.k.a. the maximum likelihood parameter estimates (MLEs). MLE computation for all but the simplest models requires the simulation of many system trajectories that are consistent with experimental data. For models with unknown parameters, this presents a computational challenge, as the generation of consistent trajectories can be an extremely rare occurrence. RESULTS: We have developed Monte Carlo Expectation-Maximization with Modified Cross-Entropy Method (MCEM2): an accelerated method for calculating MLEs that combines advances in rare event simulation with a computationally efficient version of the Monte Carlo expectation-maximization (MCEM) algorithm. Our method requires no prior knowledge regarding parameter values, and it automatically provides a multivariate parameter uncertainty estimate. We applied the method to five stochastic systems of increasing complexity, progressing from an analytically tractable pure-birth model to a computationally demanding model of yeast-polarization. Our results demonstrate that MCEM2 substantially accelerates MLE computation on all tested models when compared to a stand-alone version of MCEM. Additionally, we show how our method identifies parameter values for certain classes of models more accurately than two recently proposed computationally efficient methods. CONCLUSIONS: This work provides a novel, accelerated version of a likelihood-based parameter estimation method that can be readily applied to stochastic biochemical systems. In addition, our results suggest opportunities for added efficiency improvements that will further enhance our ability to mechanistically simulate biological processes.     

Robust parameter estimation techniques for stochastic within-host macroparasite models

We present a stochastic model of the within-host population dynamics of lymphatic filariasis, and use a simulated goodness-of-fit (GOF) method to estimate immunological parameters and their confidence intervals from experimental data. A variety of deterministic moment closure approximations to the stochastic system are explored and compared with simulation results. For the maximum GOF parameter estimates, none of the methods of closure accurately reproduce the behaviour of the stochastic model. However, direct analysis of the stochastic model demonstrates that the high levels of variation observed in the data can be reproduced without requiring parameters to vary between hosts. This indicates that the observed aggregation of parasite load may be dynamically generated by random variation in the development of an effective immune response against parasite larvae.     

A robust fed-batch feeding strategy for optimal parameter estimation for baker's yeast production

Parameter identification of structured models is often a problem in biotechnology, because the poor data situation and the number of unknown parameters only allow for inaccurate estimates. But often only a subset of all kinetic parameters of the model are of interest for production purposes, e.g. for fed-batch cultivation. These parameters should be estimated with a given accuracy. In addition, the experiments for information acquisition with respect to these parameters should be as simple as possible and should consider some practical restrictions. In this contribution a fed-batch feeding strategy is proposed to allow for an accurate estimation of yield and of critical growth rate of baker's yeast. The feeding also allows for economic and stereotyped use of staff and equipment and is therefore suitable for routine use in screening of strains and media. The overall pattern is similar to that one, usually used in production scale to minimize errors by limited model validity. After an initial phase for achieving a reproducible state three different growth rates are adjusted to cover the range of possible critical growth rates. From biomass and ethanol measurements yield and critical growth rate can be estimated with an accuracy of about 2.1%. The fermentation pattern ends up with a constant feeding rate to simulate a limited oxygen transfer rate and to allow for an uptake of residual sugar and ethanol before a dough test can be carried out. Beside experimental results simulations and sensitivity analyses are shown.List of Symbols P ethanol concentration - S substrate concentration - S _f substrate concentration in feed - T fermentation time - V fermenter volume - X biomass concentration - C measurement error covariance matrix - F Fisher information matrix - X state variables - Y output variables - X _p state sensitivity functions with respect to parameters - Y _p output sensitivity functions - e eigenvectors - k vector of limitation and inhibition parameters - n number of observations - q _in feeding stream - q _b stream for samples and ammonia feed - r vector of specific turnover rates - y vector of yields - specific weight - eigenvalues - specific growth rate - _set exponent in exponential feeding - standard deviation Dedicated to the 65th birthday of Professor Fritz Wagner.A. O. Ejiofor and B. O. Solomon are grateful to the Alexander von Humboldt Stiftung for granting them fellowships and to GBF for providing all the materials necessary for their successful research stay in Germany.     

Tools for simulating evolution of aligned genomic regions with integrated parameter estimation

Controlled simulations of genome evolution are useful for benchmarking tools. However, many simulators lack extensibility and cannot measure parameters directly from data. These issues are addressed by three new open-source programs: GSIMULATOR (for neutrally evolving DNA), SIMGRAM (for generic structured features) and SIMGENOME (for syntenic genome blocks). Each offers algorithms for parameter measurement and reconstruction of ancestral sequence. All three tools out-perform the leading neutral DNA simulator (DAWG) in benchmarks. The programs are available at .     

A parameter estimation technique for stochastic self-assembly systems and its application to human papillomavirus self-assembly

Virus capsid assembly has been a key model system for studies of complex self-assembly but it does pose some significant challenges for modeling studies. One important limitation is the difficulty of determining accurate rate parameters. The large size and rapid assembly of typical viruses make it infeasible to directly measure coat protein binding rates or deduce them from the relatively indirect experimental measures available. In this work, we develop a computational strategy to deduce coat-coat binding rate parameters for viral capsid assembly systems by fitting stochastic simulation trajectories to experimental measures of assembly progress. Our method combines quadratic response surface and quasi-gradient descent approximations to deal with the high computational cost of simulations, stochastic noise in simulation trajectories and limitations of the available experimental data. The approach is demonstrated on a light scattering trajectory for a human papillomavirus (HPV) in vitro assembly system, showing that the method can provide rate parameters that produce accurate curve fits and are in good concordance with prior analysis of the data. These fits provide an insight into potential assembly mechanisms of the in vitro system and give a basis for exploring how these mechanisms might vary between in vitro and in vivo assembly conditions.     

A cloud-based enhanced differential evolution algorithm for parameter estimation problems in computational systems biology

Metaheuristics are gaining increasing recognition in many research areas, computational systems biology among them. Recent advances in metaheuristics can be helpful in locating the vicinity of the global solution in reasonable computation times, with Differential Evolution (DE) being one of the most popular methods. However, for most realistic applications, DE still requires excessive computation times. With the advent of Cloud Computing effortless access to large number of distributed resources has become more feasible, and new distributed frameworks, like Spark, have been developed to deal with large scale computations on commodity clusters and cloud resources. In this paper we propose a parallel implementation of an enhanced DE using Spark. The proposal drastically reduces the execution time, by means of including a selected local search and exploiting the available distributed resources. The performance of the proposal has been thoroughly assessed using challenging parameter estimation problems from the domain of computational systems biology. Two different platforms have been used for the evaluation, a local cluster and the Microsoft Azure public cloud. Additionally, it has been also compared with other parallel approaches, another cloud-based solution (a MapReduce implementation) and a traditional HPC solution (a MPI implementation)     

Cytochrome b562 variants: a library for examining redox potential evolution

A general understanding of how cytochromes evolve within a fixed structure to optimize redox potential for specific bioenergetic processes does not exist. Toward this end, a library approach is used to investigate the range and distribution of redox potential which occurs when all sequence space available through mutation at two positions is examined within a fixed structural motif. Random mutation of Phe61 and Phe65 of cytochrome b562 (E. coli), and subsequent examination of a statistically significant sampling of this library, demonstrates that the redox potential can vary over 100 mV (>25% of the known accessible potential in native proteins with axial His-Met ligation) through mutation at these two positions. The redox potential of the wild-type protein occurs at an extremum of the distribution observed, indicating that Phe61 and Phe65 were most likely naturally selected to differentially stabilize the reduced state of the protein. At the other extremum, a compositionally conservative set of mutations (F61I, F65Y) leads to a 100 mV shift in the redox equilibrium toward the oxidized state. NMR analyses indicate that a charge-dipole interaction which results from mutation of phenylalanine to tyrosine at position 65 may be responsible.     

Analysis of the respiro-fermentative growth of Saccharomyces cerevisiae on glucose in a fed-batch fermentation strategy for accurate parameter estimation

Analysis of the respiro-fermentive growth of a strain of Saccharomyces cerevisiae, DSM 2155 on glucose, in a simulated 5-phase feeding strategy of fedbatch cultures executed on the Universal BIoprocess CONtrol (UBICON) system, was carried out. There was a good agreement between the estimated and the simulated values of specific growth rates. In this study, which was designed to span 0.20–0.23 h^–1 growth rates before returning to lower growth rates, the critical dilution rate at which the switch between purely oxidative and respiro-fermentative growth takes was not observed. The biomass yield, specific substrate uptake and O₂ consumption rates as well as the consistency of the data using both carbon and available electron balances were examined. A high average value of true biomass energetic yield, _max = 0.707, and a low value of maintenance coefficient, m_e = 0.0114 h^–1, were obtained indicating that the organism was in no danger from the ethanol produced as a high-density fermentation with a yeast concentration above 54 g 1^–1 was possible within a period of 24 h. The yeast produced also had good dough-leavening characteristics. Thus it is possible to operate a yeast plant without resorting to using respiratory quotient, which may be problematic, as the controlling parameter.     

CADLIVE optimizer: web-based parameter estimation for dynamic models

ABSTRACT: Computer simulation has been an important technique to capture the dynamics of biochemical networks. In most networks, however, few kinetic parameters have been measured in vivo because of experimental complexity. We develop a kinetic parameter estimation system, named the CADLIVE Optimizer, which comprises genetic algorithms-based solvers with a graphical user interface. This optimizer is integrated into the CADLIVE Dynamic Simulator to attain efficient simulation for dynamic models.     

Mutational order: a major stochastic process in evolution

Computer simulations in which selection acts on a quantitative character show that the randomness of mutations can contribute significantly to evolutionary divergence between populations. In different populations, different advantageous mutations occur, and are selected to fixation, so that the populations diverge even when they are initially identical, and are subject to identical selection. This stochastic process is distinct from random genetic drift. In some circumstances (large populations or strong selection, or both) mutational order can be greatly more important than random drift in bringing about divergence. It can generate a 'disconnection' between evolution at the phenotypic and genotypic levels, and can give rise to a rough 'molecular clock', albeit episodic, that is driven by selection. In the absence of selection, mutational order has little or no effect.     

SBML-PET: a Systems Biology Markup Language-based parameter estimation tool

The estimation of model parameters from experimental data remains a bottleneck for a major breakthrough in systems biology. We present a Systems Biology Markup Language (SBML) based Parameter Estimation Tool (SBML-PET). The tool is designed to enable parameter estimation for biological models including signaling pathways, gene regulation networks and metabolic pathways. SBML-PET supports import and export of the models in the SBML format. It can estimate the parameters by fitting a variety of experimental data from different experimental conditions. SBML-PET has a unique feature of supporting event definition in the SMBL model. SBML models can also be simulated in SBML-PET. Stochastic Ranking Evolution Strategy (SRES) is incorporated in SBML-PET for parameter estimation jobs. A classic ODE Solver called ODEPACK is used to solve the Ordinary Differential Equation (ODE) system. AVAILABILITY: http://sysbio.molgen.mpg.de/SBML-PET/. The website also contains detailed documentation for SBML-PET.     

Individual based modeling and parameter estimation for a Lotka-Volterra system

Stochastic component, inevitable in biological systems, makes problematic the estimation of the model parameters from a single sequence of measurements, despite the complete knowledge of the system. We studied the problem of parameter estimation using individual-based computer simulations of a 'Lotka-Volterra world'. Two kinds (species) of particles--X (preys) and Y (predators)--moved on a sphere according to deterministic rules and at the collision (interaction) of X and Y the particle X was changed to a new particle Y. Birth of preys and death of predators were simulated by addition of X and removal of Y, respectively, according to exponential probability distributions. With this arrangement of the system, the numbers of particles of each kind might be described by the Lotka-Volterra equations. The simulations of the system with low (200-400 particles on average) number of individuals showed unstable oscillations of the population size. In some simulation runs one of the species became extinct. Nevertheless, the oscillations had some generic properties (e.g. mean, in one simulation run, oscillation period, mean ratio of the amplitudes of the consecutive maxima of X and Y numbers, etc.) characteristic for the solutions of the Lotka-Volterra equations. This observation made it possible to estimate the four parameters of the Lotka-Volterra model with high accuracy and good precision. The estimation was performed using the integral form of the Lotka-Volterra equations and two parameter linear regression for each oscillation cycle separately. We conclude that in spite of the irregular time course of the number of individuals in each population due to stochastic intraspecies component, the generic features of the simulated system evolution can provide enough information for quantitative estimation of the system parameters.     

SParSE++: improved event-based stochastic parameter search

Background

Despite the increasing availability of high performance computing capabilities, analysis and characterization of stochastic biochemical systems remain a computational challenge. To address this challenge, the Stochastic Parameter Search for Events (SParSE) was developed to automatically identify reaction rates that yield a probabilistic user-specified event. SParSE consists of three main components: the multi-level cross-entropy method, which identifies biasing parameters to push the system toward the event of interest, the related inverse biasing method, and an optional interpolation of identified parameters. While effective for many examples, SParSE depends on the existence of a sufficient amount of intrinsic stochasticity in the system of interest. In the absence of this stochasticity, SParSE can either converge slowly or not at all.

Results

We have developed SParSE++, a substantially improved algorithm for characterizing target events in terms of system parameters. SParSE++ makes use of a series of novel parameter leaping methods that accelerate the convergence rate to the target event, particularly in low stochasticity cases. In addition, the interpolation stage is modified to compute multiple interpolants and to choose the optimal one in a statistically rigorous manner. We demonstrate the performance of SParSE++ on four example systems: a birth-death process, a reversible isomerization model, SIRS disease dynamics, and a yeast polarization model. In all four cases, SParSE++ shows significantly improved computational efficiency over SParSE, with the largest improvements resulting from analyses with the strictest error tolerances.

Conclusions

As researchers continue to model realistic biochemical systems, the need for efficient methods to characterize target events will grow. The algorithmic advancements provided by SParSE++ fulfill this need, enabling characterization of computationally intensive biochemical events that are currently resistant to analysis.

     

Superior feature-set ranking for small samples using bolstered error estimation

MOTIVATION: Ranking feature sets is a key issue for classification, for instance, phenotype classification based on gene expression. Since ranking is often based on error estimation, and error estimators suffer to differing degrees of imprecision in small-sample settings, it is important to choose a computationally feasible error estimator that yields good feature-set ranking. RESULTS: This paper examines the feature-ranking performance of several kinds of error estimators: resubstitution, cross-validation, bootstrap and bolstered error estimation. It does so for three classification rules: linear discriminant analysis, three-nearest-neighbor classification and classification trees. Two measures of performance are considered. One counts the number of the truly best feature sets appearing among the best feature sets discovered by the error estimator and the other computes the mean absolute error between the top ranks of the truly best feature sets and their ranks as given by the error estimator. Our results indicate that bolstering is superior to bootstrap, and bootstrap is better than cross-validation, for discovering top-performing feature sets for classification when using small samples. A key issue is that bolstered error estimation is tens of times faster than bootstrap, and faster than cross-validation, and is therefore feasible for feature-set ranking when the number of feature sets is extremely large.     

Libsequence: a C++ class library for evolutionary genetic analysis

A C++ class library is available to facilitate the implementation of software for genomics and sequence polymorphism analysis. The library implements methods for data manipulation and the calculation of several statistics commonly used to analyze SNP data. The object-oriented design of the library is intended to be extensible, allowing users to design custom classes for their own needs. In addition, routines are provided to process samples generated by a widely used coalescent simulation. AVAILABILITY: The source code (in C++) is available from http://www.molpopgen.org     

Computational approaches for RNA energy parameter estimation

Methods for efficient and accurate prediction of RNA structure are increasingly valuable, given the current rapid advances in understanding the diverse functions of RNA molecules in the cell. To enhance the accuracy of secondary structure predictions, we developed and refined optimization techniques for the estimation of energy parameters. We build on two previous approaches to RNA free-energy parameter estimation: (1) the Constraint Generation (CG) method, which iteratively generates constraints that enforce known structures to have energies lower than other structures for the same molecule; and (2) the Boltzmann Likelihood (BL) method, which infers a set of RNA free-energy parameters that maximize the conditional likelihood of a set of reference RNA structures. Here, we extend these approaches in two main ways: We propose (1) a max-margin extension of CG, and (2) a novel linear Gaussian Bayesian network that models feature relationships, which effectively makes use of sparse data by sharing statistical strength between parameters. We obtain significant improvements in the accuracy of RNA minimum free-energy pseudoknot-free secondary structure prediction when measured on a comprehensive set of 2518 RNA molecules with reference structures. Our parameters can be used in conjunction with software that predicts RNA secondary structures, RNA hybridization, or ensembles of structures. Our data, software, results, and parameter sets in various formats are freely available at http://www.cs.ubc.ca/labs/beta/Projects/RNA-Params.     

A general method for biological parameter estimation

Summary A general observer-based estimator method is developed and applied for process modelling and monitoring. This parameter estimation technique was successfully applied to a L-lysine fermentation process. It was a useful tool to detect the effect of major culture conditions on cell growth and product synthesis. It can also be used for the development of adaptive optimal control schemes.