Benchmarks for identification of ordinary differential equations from time series data

SBbadger: biochemical reaction networks with definable degree distributions

MOTIVATION

An essential step in developing computational tools for the inference, optimization and simulation of biochemical reaction networks is gauging tool performance against earlier efforts using an appropriate set of benchmarks. General strategies for the assembly of benchmark models include collection from the literature, creation via subnetwork extraction and de novo generation. However, with respect to biochemical reaction networks, these approaches and their associated tools are either poorly suited to generate models that reflect the wide range of properties found in natural biochemical networks or to do so in numbers that enable rigorous statistical analysis.

RESULTS

In this work, we present SBbadger, a python-based software tool for the generation of synthetic biochemical reaction or metabolic networks with user-defined degree distributions, multiple available kinetic formalisms and a host of other definable properties. SBbadger thus enables the creation of benchmark model sets that reflect properties of biological systems and generate the kinetics and model structures typically targeted by computational analysis and inference software. Here, we detail the computational and algorithmic workflow of SBbadger, demonstrate its performance under various settings, provide sample outputs and compare it to currently available biochemical reaction network generation software.

AVAILABILITY AND IMPLEMENTATION

SBbadger is implemented in Python and is freely available at https://github.com/sys-bio/SBbadger and via PyPI at https://pypi.org/project/SBbadger/. Documentation can be found at https://SBbadger.readthedocs.io.

SUPPLEMENTARY INFORMATION

Supplementary data are available at Bioinformatics online.

Discrete Biochemical Systems Theory

Almost every biomedical systems analysis requires early decisions regarding the choice of the most suitable representations to be used. De facto the most prevalent choice is a system of ordinary differential equations (ODEs). This framework is very popular because it is flexible and fairly easy to use. It is also supported by an enormous array of stand-alone programs for analysis, including many distinct numerical solvers that are implemented in the main programming languages. Having selected ODEs, the modeler must then choose a mathematical format for the equations. This selection is not trivial as nearly unlimited options exist and there is seldom objective guidance. The typical choices include ad hoc representations, default models like mass-action or Lotka-Volterra equations, and generic approximations. Within the realm of approximations, linear models are typically successful for analyses of engineered systems, but they are not as appropriate for biomedical phenomena, which often display nonlinear features such as saturation, threshold effects or limit cycle oscillations, and possibly even chaos. Power-law approximations are simple but overcome these limitations. They are the key ingredient of Biochemical Systems Theory (BST), which uses ODEs exclusively containing power-law representations for all processes within a model. BST models cover a vast repertoire of nonlinear responses and, at the same time, have structural properties that are advantageous for a wide range of analyses. Nonetheless, as all ODE models, the BST approach has limitations. In particular, it is not always straightforward to account for genuine discreteness, time delays, and stochastic processes. As a new option, we therefore propose here an alternative to BST in the form of discrete Biochemical Systems Theory (dBST). dBST models have the same generality and practicality as their BST-ODE counterparts, but they are readily implemented even in situations where ODEs struggle. As a case study, we illustrate dBST applied to the dynamics of the aryl hydrocarbon receptor (AhR), a signal transduction system that simultaneously involves time delays and stochasticity.

The best models of metabolism

Biochemical systems are among of the oldest application areas of mathematical modeling. Spanning a time period of over one hundred years, the repertoire of options for structuring a model and for formulating reactions has been constantly growing, and yet, it is still unclear whether or to what degree some models are better than others and how the modeler is to choose among them. In fact, the variety of options has become overwhelming and difficult to maneuver for novices and experts alike. This review outlines the metabolic model design process and discusses the numerous choices for modeling frameworks and mathematical representations. It tries to be inclusive, even though it cannot be complete, and introduces the various modeling options in a manner that is as unbiased as that is feasible. However, the review does end with personal recommendations for the choices of default models. WIREs Syst Biol Med 2017, 9:e1391. doi: 10.1002/wsbm.1391 For further resources related to this article, please visit the WIREs website.

Nonparametric dynamic modeling

Challenging as it typically is, the estimation of parameter values seems to be an unavoidable step in the design and implementation of any dynamic model. Here, we demonstrate that it is possible to set up, diagnose, and simulate dynamic models without the need to estimate parameter values, if the situation is favorable. Specifically, it is possible to establish nonparametric models for nonlinear compartment models, including metabolic pathway models, if sufficiently many high-quality time series data are available that describe the biological phenomenon under investigation in an appropriate and representative manner. The proposed nonparametric strategy is a variant of the method of Dynamic Flux Estimation (DFE), which permits the estimation of numerical flux profiles from metabolic time series data. However, instead of attempting to formulate these numerical profiles as explicit functions and to optimize their parameter values, as it is done in DFE, the metabolite and flux profiles are used here directly as a scaffold for a library from which values are interpolated and retrieved for the simulation of the differential equations describing the model. Beyond simulations, the proposed methods render it possible to determine steady states from non-steady state data, perform sensitivity analyses, and estimate the Jacobian of the system at a steady state.

JuPOETs: a constrained multiobjective optimization approach to estimate biochemical model ensembles in the Julia programming language

Background

Ensemble modeling is a promising approach for obtaining robust predictions and coarse grained population behavior in deterministic mathematical models. Ensemble approaches address model uncertainty by using parameter or model families instead of single best-fit parameters or fixed model structures. Parameter ensembles can be selected based upon simulation error, along with other criteria such as diversity or steady-state performance. Simulations using parameter ensembles can estimate confidence intervals on model variables, and robustly constrain model predictions, despite having many poorly constrained parameters.

Results

In this software note, we present a multiobjective based technique to estimate parameter or models ensembles, the Pareto Optimal Ensemble Technique in the Julia programming language (JuPOETs). JuPOETs integrates simulated annealing with Pareto optimality to estimate ensembles on or near the optimal tradeoff surface between competing training objectives. We demonstrate JuPOETs on a suite of multiobjective problems, including test functions with parameter bounds and system constraints as well as for the identification of a proof-of-concept biochemical model with four conflicting training objectives. JuPOETs identified optimal or near optimal solutions approximately six-fold faster than a corresponding implementation in Octave for the suite of test functions. For the proof-of-concept biochemical model, JuPOETs produced an ensemble of parameters that gave both the mean of the training data for conflicting data sets, while simultaneously estimating parameter sets that performed well on each of the individual objective functions.

Conclusions

JuPOETs is a promising approach for the estimation of parameter and model ensembles using multiobjective optimization. JuPOETs can be adapted to solve many problem types, including mixed binary and continuous variable types, bilevel optimization problems and constrained problems without altering the base algorithm. JuPOETs is open source, available under an MIT license, and can be installed using the Julia package manager from the JuPOETs GitHub repository

Modelling and analysis of an ensemble of eukaryotic translation initiation models

Programmed protein synthesis plays an important role in the cell cycle. Deregulated translation has been observed in several cancers. In this study, the authors constructed an ensemble of mathematical models describing the integration of growth factor signals with translation initiation. Using these models, the authors estimated critical structural features of the translation architecture. Sensitivity and robustness analysis with and without growth factors suggested that a balance between competing regulatory programmes governed translation initiation. Proteins such as Akt and mTor promoted initiation by integrating growth factor signals with the assembly of the 80S initiation complex. However, negative regulators such as PTEN and 4EBP1 restrained initiation in the absence of stimulation. Other proteins such as eIF4E were also found to be structurally critical as deletion of amplification of these components resulted in a network incapable of nominal operation. These findings could help understand the molecular basis of translation deregulation observed in cancer and perhaps lead to new anti-cancer therapeutic strategies. [Includes supplementary material].

BioPreDyn-bench: a suite of benchmark problems for dynamic modelling in systems biology

Background

Dynamic modelling is one of the cornerstones of systems biology. Many research efforts are currently being invested in the development and exploitation of large-scale kinetic models. The associated problems of parameter estimation (model calibration) and optimal experimental design are particularly challenging. The community has already developed many methods and software packages which aim to facilitate these tasks. However, there is a lack of suitable benchmark problems which allow a fair and systematic evaluation and comparison of these contributions.

Results

Here we present BioPreDyn-bench, a set of challenging parameter estimation problems which aspire to serve as reference test cases in this area. This set comprises six problems including medium and large-scale kinetic models of the bacterium E. coli, baker’s yeast S. cerevisiae, the vinegar fly D. melanogaster, Chinese Hamster Ovary cells, and a generic signal transduction network. The level of description includes metabolism, transcription, signal transduction, and development. For each problem we provide (i) a basic description and formulation, (ii) implementations ready-to-run in several formats, (iii) computational results obtained with specific solvers, (iv) a basic analysis and interpretation.

Conclusions

This suite of benchmark problems can be readily used to evaluate and compare parameter estimation methods. Further, it can also be used to build test problems for sensitivity and identifiability analysis, model reduction and optimal experimental design methods. The suite, including codes and documentation, can be freely downloaded from the BioPreDyn-bench website, https://sites.google.com/site/biopredynbenchmarks/.

Electronic supplementary material

The online version of this article (doi:10.1186/s12918-015-0144-4) contains supplementary material, which is available to authorized users.

Automated analysis of routinely generated preclinical pharmacokinetic and pharmacodynamic data

Model-based analysis of routinely generated pharmacokinetic and pharmacodynamic (PK-PD) data is a key component of preclinical drug discovery. The work process of such analyses can be automated by properly designed computer programs that reduce the number of manual steps, resulting in time saving and significantly fewer errors. Critical decisions can still be made by modelers. Using concrete animal data examples this paper illustrates when, and demonstrates how, automated PK-PD approaches can be used and what benefits they offer to the modeling and simulation community. Specifically, we describe two compound optimization case studies from drug discovery projects, and also demonstrate how a subsequent optimization step to predict the human dose can be coupled to an automated approach.

ODEion--a software module for structural identification of ordinary differential equations

In the systems biology field, algorithms for structural identification of ordinary differential equations (ODEs) have mainly focused on fixed model spaces like S-systems and/or on methods that require sufficiently good data so that derivatives can be accurately estimated. There is therefore a lack of methods and software that can handle more general models and realistic data. We present ODEion, a software module for structural identification of ODEs. Main characteristic features of the software are: • The model space is defined by arbitrary user-defined functions that can be nonlinear in both variables and parameters, such as for example chemical rate reactions. • ODEion implements computationally efficient algorithms that have been shown to efficiently handle sparse and noisy data. It can run a range of realistic problems that previously required a supercomputer. • ODEion is easy to use and provides SBML output. We describe the mathematical problem, the ODEion system itself, and provide several examples of how the system can be used. Available at: http://www.odeidentification.org.

Biochemical Systems Theory: A Review

Biochemical systems theory (BST) is the foundation for a set of analytical andmodeling tools that facilitate the analysis of dynamic biological systems. This paper depicts major developments in BST up to the current state of the art in 2012. It discusses its rationale, describes the typical strategies and methods of designing, diagnosing, analyzing, and utilizing BST models, and reviews areas of application. The paper is intended as a guide for investigators entering the fascinating field of biological systems analysis and as a resource for practitioners and experts.

Multiscale models of breast cancer progression

Breast cancer initiation, invasion and metastasis span multiple length and time scales. Molecular events at short length scales lead to an initial tumorigenic population, which left unchecked by immune action, acts at increasingly longer length scales until eventually the cancer cells escape from the primary tumor site. This series of events is highly complex, involving multiple cell types interacting with (and shaping) the microenvironment. Multiscale mathematical models have emerged as a powerful tool to quantitatively integrate the convective-diffusion-reaction processes occurring on the systemic scale, with the molecular signaling processes occurring on the cellular and subcellular scales. In this study, we reviewed the current state of the art in cancer modeling across multiple length scales, with an emphasis on the integration of intracellular signal transduction models with pro-tumorigenic chemical and mechanical microenvironmental cues. First, we reviewed the underlying biomolecular origin of breast cancer, with a special emphasis on angiogenesis. Then, we summarized the development of tissue engineering platforms which could provide high-fidelity ex vivo experimental models to identify and validate multiscale simulations. Lastly, we reviewed top-down and bottom-up multiscale strategies that integrate subcellular networks with the microenvironment. We present models of a variety of cancers, in addition to breast cancer specific models. Taken together, we expect as the sophistication of the simulations increase, that multiscale modeling and bottom-up agent-based models in particular will become an increasingly important platform technology for basic scientific discovery, as well as the identification and validation of potentially novel therapeutic targets.

Estimation of dynamic flux profiles from metabolic time series data

Background

Advances in modern high-throughput techniques of molecular biology have enabled top-down approaches for the estimation of parameter values in metabolic systems, based on time series data. Special among them is the recent method of dynamic flux estimation (DFE), which uses such data not only for parameter estimation but also for the identification of functional forms of the processes governing a metabolic system. DFE furthermore provides diagnostic tools for the evaluation of model validity and of the quality of a model fit beyond residual errors. Unfortunately, DFE works only when the data are more or less complete and the system contains as many independent fluxes as metabolites. These drawbacks may be ameliorated with other types of estimation and information. However, such supplementations incur their own limitations. In particular, assumptions must be made regarding the functional forms of some processes and detailed kinetic information must be available, in addition to the time series data.

Results

The authors propose here a systematic approach that supplements DFE and overcomes some of its shortcomings. Like DFE, the approach is model-free and requires only minimal assumptions. If sufficient time series data are available, the approach allows the determination of a subset of fluxes that enables the subsequent applicability of DFE to the rest of the flux system. The authors demonstrate the procedure with three artificial pathway systems exhibiting distinct characteristics and with actual data of the trehalose pathway in Saccharomyces cerevisiae.

Conclusions

The results demonstrate that the proposed method successfully complements DFE under various situations and without a priori assumptions regarding the model representation. The proposed method also permits an examination of whether at all, to what degree, or within what range the available time series data can be validly represented in a particular functional format of a flux within a pathway system. Based on these results, further experiments may be designed to generate data points that genuinely add new information to the structure identification and parameter estimation tasks at hand.

Modeling gene regulation networks using ordinary differential equations

Gene regulation networks are composed of transcription factors, their interactions, and targets. It is of great interest to reconstruct and study these regulatory networks from genomics data. Ordinary differential equations (ODEs) are popular tools to model the dynamic system of gene regulation networks. Although the form of ODEs is often provided based on expert knowledge, the values for ODE parameters are seldom known. It is a challenging problem to infer ODE parameters from gene expression data, because the ODEs do not have analytic solutions and the time-course gene expression data are usually sparse and associated with large noise. In this chapter, we review how the generalized profiling method can be applied to obtain estimates for ODE parameters from the time-course gene expression data. We also summarize the consistency and asymptotic normality results for the generalized profiling estimates.

SimBioNeT: a simulator of biological network topology

Studying biological networks at topological level is a major issue in computational biology studies and simulation is often used in this context, either to assess reverse engineering algorithms or to investigate how topological properties depend on network parameters. In both contexts, it is desirable for a topology simulator to reproduce the current knowledge on biological networks, to be able to generate a number of networks with the same properties and to be flexible with respect to the possibility to mimic networks of different organisms. We propose a biological network topology simulator, SimBioNeT, in which module structures of different type and size are replicated at different level of network organization and interconnected, so to obtain the desired degree distribution, e.g., scale free, and a clustering coefficient constant with the number of nodes in the network, a typical characteristic of biological networks. Empirical assessment of the ability of the simulator to reproduce characteristic properties of biological network and comparison with E. coli and S. cerevisiae transcriptional networks demonstrates the effectiveness of our proposal.

Modeling and analysis of retinoic acid induced differentiation of uncommitted precursor cells

Manipulation of differentiation programs has therapeutic potential in a spectrum of human cancers and neurodegenerative disorders. In this study, we integrated computational and experimental methods to unravel the response of a lineage uncommitted precursor cell-line, HL-60, to Retinoic Acid (RA). HL-60 is a human myeloblastic leukemia cell-line used extensively to study human differentiation programs. Initially, we focused on the role of the BLR1 receptor in RA-induced differentiation and G1/0-arrest in HL-60. BLR1, a putative G protein-coupled receptor expressed following RA exposure, is required for RA-induced cell-cycle arrest and differentiation and causes persistent MAPK signaling. A mathematical model of RA-induced cell-cycle arrest and differentiation was formulated and tested against BLR1 wild-type (wt) knock-out and knock-in HL-60 cell-lines with and without RA. The current model described the dynamics of 729 proteins and protein complexes interconnected by 1356 interactions. An ensemble strategy was used to compensate for uncertain model parameters. The ensemble of HL-60 models recapitulated the positive feedback between BLR1 and MAPK signaling. The ensemble of models also correctly predicted Rb and p47phox regulation and the correlation between p21-CDK4-cyclin D formation and G1/0-arrest following exposure to RA. Finally, we investigated the robustness of the HL-60 network architecture to structural perturbations and generated experimentally testable hypotheses for future study. Taken together, the model presented here was a first step toward a systematic framework for analysis of programmed differentiation. These studies also demonstrated that mechanistic network modeling can help prioritize experimental directions by generating falsifiable hypotheses despite uncertainty.

Automated analysis of biological oscillator models using mode decomposition

MOTIVATION

Oscillating signals produced by biological systems have shapes, described by their Fourier spectra, that can potentially reveal the mechanisms that generate them. Extracting this information from measured signals is interesting for the validation of theoretical models, discovery and classification of interaction types, and for optimal experiment design.

RESULTS

An automated workflow is described for the analysis of oscillating signals. A software package is developed to match signal shapes to hundreds of a priori viable model structures defined by a class of first-order differential equations. The package computes parameter values for each model by exploiting the mode decomposition of oscillating signals and formulating the matching problem in terms of systems of simultaneous polynomial equations. On the basis of the computed parameter values, the software returns a list of models consistent with the data. In validation tests with synthetic datasets, it not only shortlists those model structures used to generate the data but also shows that excellent fits can sometimes be achieved with alternative equations. The listing of all consistent equations is indicative of how further invalidation might be achieved with additional information. When applied to data from a microarray experiment on mice, the procedure finds several candidate model structures to describe interactions related to the circadian rhythm. This shows that experimental data on oscillators is indeed rich in information about gene regulation mechanisms.

AVAILABILITY

The software package is available at http://babylone.ulb.ac.be/autoosc/.

Continuous-time modeling of cell fate determination in Arabidopsis flowers

Background

The genetic control of floral organ specification is currently being investigated by various approaches, both experimentally and through modeling. Models and simulations have mostly involved boolean or related methods, and so far a quantitative, continuous-time approach has not been explored.

Results

We propose an ordinary differential equation (ODE) model that describes the gene expression dynamics of a gene regulatory network that controls floral organ formation in the model plant Arabidopsis thaliana. In this model, the dimerization of MADS-box transcription factors is incorporated explicitly. The unknown parameters are estimated from (known) experimental expression data. The model is validated by simulation studies of known mutant plants.

Conclusions

The proposed model gives realistic predictions with respect to independent mutation data. A simulation study is carried out to predict the effects of a new type of mutation that has so far not been made in Arabidopsis, but that could be used as a severe test of the validity of the model. According to our predictions, the role of dimers is surprisingly important. Moreover, the functional loss of any dimer leads to one or more phenotypic alterations.

Ensembles of signal transduction models using Pareto Optimal Ensemble Techniques (POETs)

Mathematical modeling of complex gene expression programs is an emerging tool for understanding disease mechanisms. However, identification of large models sometimes requires training using qualitative, conflicting or even contradictory data sets. One strategy to address this challenge is to estimate experimentally constrained model ensembles using multiobjective optimization. In this study, we used Pareto Optimal Ensemble Techniques (POETs) to identify a family of proof-of-concept signal transduction models. POETs integrate Simulated Annealing (SA) with Pareto optimality to identify models near the optimal tradeoff surface between competing training objectives. We modeled a prototypical-signaling network using mass-action kinetics within an ordinary differential equation (ODE) framework (64 ODEs in total). The true model was used to generate synthetic immunoblots from which the POET algorithm identified the 117 unknown model parameters. POET generated an ensemble of signaling models, which collectively exhibited population-like behavior. For example, scaled gene expression levels were approximately normally distributed over the ensemble following the addition of extracellular ligand. Also, the ensemble recovered robust and fragile features of the true model, despite significant parameter uncertainty. Taken together, these results suggest that experimentally constrained model ensembles could capture qualitatively important network features without exact parameter information.

Analysis of the molecular networks in androgen dependent and independent prostate cancer revealed fragile and robust subsystems

Androgen ablation therapy is currently the primary treatment for metastatic prostate cancer. Unfortunately, in nearly all cases, androgen ablation fails to permanently arrest cancer progression. As androgens like testosterone are withdrawn, prostate cancer cells lose their androgen sensitivity and begin to proliferate without hormone growth factors. In this study, we constructed and analyzed a mathematical model of the integration between hormone growth factor signaling, androgen receptor activation, and the expression of cyclin D and Prostate-Specific Antigen in human LNCaP prostate adenocarcinoma cells. The objective of the study was to investigate which signaling systems were important in the loss of androgen dependence. The model was formulated as a set of ordinary differential equations which described 212 species and 384 interactions, including both the mRNA and protein levels for key species. An ensemble approach was chosen to constrain model parameters and to estimate the impact of parametric uncertainty on model predictions. Model parameters were identified using 14 steady-state and dynamic LNCaP data sets taken from literature sources. Alterations in the rate of Prostatic Acid Phosphatase expression was sufficient to capture varying levels of androgen dependence. Analysis of the model provided insight into the importance of network components as a function of androgen dependence. The importance of androgen receptor availability and the MAPK/Akt signaling axes was independent of androgen status. Interestingly, androgen receptor availability was important even in androgen-independent LNCaP cells. Translation became progressively more important in androgen-independent LNCaP cells. Further analysis suggested a positive synergy between the MAPK and Akt signaling axes and the translation of key proliferative markers like cyclin D in androgen-independent cells. Taken together, the results support the targeting of both the Akt and MAPK pathways. Moreover, the analysis suggested that direct targeting of the translational machinery, specifically eIF4E, could be efficacious in androgen-independent prostate cancers.