In this paper we develop a new method for numerically approximating sensitivities in parameter-dependent ordinary differential equations (ODEs). Our approach, intended for situations where the standard forward and adjoint sensitivity analyses become too computationally costly for practical purposes, is based on the Peano-Baker series from control theory. Using this series, we construct a representation of the sensitivity matrix $$\textbf{S}$$ S and, from this representation, a numerical method for approximating $$\textbf{S}$$ S . We prove that, under standard regularity assumptions, the error of our method scales as $$\mathcal {O}(\Delta t^2_{\textrm{max}})$$ O ( Δ t max 2 ) , where $$\Delta t_{\textrm{max}}$$ Δ t max is the largest time step used when numerically solving the ODE. We illustrate the performance of the method in several numerical experiments, taken from both the systems biology setting and more classical dynamical systems. The experiments show the sought-after improvement in running time of our method compared to the forward sensitivity approach. In experiments involving a random linear system, the forward approach requires roughly $$\sqrt{n}$$ n longer computational time, where n is the dimension of the parameter space, than our proposed method.

We present an R-package developed for modeling of biochemical reaction networks, uncertainty quantification (UQ) and sensitivity analysis (SA). Estimating parameters and quantifying their uncertainty (and resulting prediction uncertainty), is required for data-driven systems biology modeling. Sampling methods need to be efficient when confronted with high-dimensional, correlated parameter distributions. We have developed the UQSA package to be fast for this problem class and work well with other tools for modelling. We aim for simplicity, and part of that is our use of the SBtab format for the unified storage of model and data. Our tool-set is modular enough, that parts can be replaced. We use intermediate formats that are not hidden from the user to make this feasible. UQ is performed through Markov chain Monte Carlo (MCMC) sampling in an Approximate Bayesian Computation (ABC) setting. This can be followed by a variance-decomposition based global sensitivity analysis. If needed, complex parameter distributions can be described, evaluated, and sampled from, with the help of Vine-copulas that are available in R. This approach is especially useful when new experimental data become available, and a previously calibrated model needs to be updated. Implementation: R is a high-level language and allows the use of sophisticated statistical methods. The ode solver we used is written in C (gsl_odeiv2, interface to R is ours). We use the SBtab tabular format for the model description, as well as the data and an event system to be able to model inputs frequently encountered in systems biology and neuroscience. The code has been tested on one node with 256 cores of a computing cluster, but smaller examples are included in the repository that can be run on a laptop. Source code: https://github.com/icpm-kth/uqsa

Modeling in neuroscience occurs at the intersection of different points of view and approaches. Typically, hypothesis-driven modeling brings a question into focus so that a model is constructed to investigate a specific hypothesis about how the system works or why certain phenomena are observed. Data-driven modeling, on the other hand, follows a more unbiased approach, with model construction informed by the computationally intensive use of data. At the same time, researchers employ models at different biological scales and at different levels of abstraction. Combining these models while validating them against experimental data increases understanding of the multiscale brain. However, a lack of interoperability, transparency, and reusability of both models and the workflows used to construct them creates barriers for the integration of models representing different biological scales and built using different modeling philosophies. We argue that the same imperatives that drive resources and policy for data - such as the FAIR (Findable, Accessible, Interoperable, Reusable) principles - also support the integration of different modeling approaches. The FAIR principles require that data be shared in formats that are Findable, Accessible, Interoperable, and Reusable. Applying these principles to models and modeling workflows, as well as the data used to constrain and validate them, would allow researchers to find, reuse, question, validate, and extend published models, regardless of whether they are implemented phenomenologically or mechanistically, as a few equations or as a multiscale, hierarchical system. To illustrate these ideas, we use a classical synaptic plasticity model, the Bienenstock-Cooper-Munro rule, as an example due to its long history, different levels of abstraction, and implementation at many scales.

Neuroscience incorporates knowledge from a range of scales, from single molecules to brain wide neural networks. Modeling is a valuable tool in understanding processes at a single scale or the interactions between two adjacent scales and researchers use a variety of different software tools in the model building and analysis process. Here we focus on the scale of biochemical pathways, which is one of the main objects of study in systems biology. While systems biology is among the more standardized fields, conversion between different model formats and interoperability between various tools is still somewhat problematic. To offer our take on tackling these shortcomings and by keeping in mind the FAIR (findability, accessibility, interoperability, reusability) data principles, we have developed a workflow for building and analyzing biochemical pathway models, using pre-existing tools that could be utilized for the storage and refinement of models in all phases of development. We have chosen the SBtab format which allows the storage of biochemical models and associated data in a single file and provides a human readable set of syntax rules. Next, we implemented custom-made MATLAB® scripts to perform parameter estimation and global sensitivity analysis used in model refinement. Additionally, we have developed a web-based application for biochemical models that allows simulations with either a network free solver or stochastic solvers and incorporating geometry. Finally, we illustrate convertibility and use of a biochemical model in a biophysically detailed single neuron model by running multiscale simulations in NEURON. Using this workflow, we can simulate the same model in three different simulators, with a smooth conversion between the different model formats, enhancing the characterization of different aspects of the model.

We present a workflow for kinetic modeling of biocatalytic reactions which combines methods from Bayesian learning and uncertainty quantification for model calibration, model selection, evaluation, and model reduction in a consistent statistical framework. Our workflow is particularly tailored to sparse data settings in which a considerable variability of the parameters remains after the models have been adapted to available data, a ubiquitous problem in many real-world applications. Our workflow is exemplified on an enzyme-catalyzed two-substrate reaction mechanism describing the symmetric carboligation of 3,5-dimethoxy-benzaldehyde to (R)-3,3 ',5,5 '-tetramethoxybenzoin catalyzed by benzaldehyde lyase from Pseudomonas fluorescens. Results indicate a substrate-dependent inactivation of enzyme, which is in accordance with other recent studies.

Motivation Dynamical models describing intracellular phenomena are increasing in size and complexity as more information is obtained from experiments. These models are often over-parameterized with respect to the quantitative data used for parameter estimation, resulting in uncertainty in the individual parameter estimates as well as in the predictions made from the model. Here we combine Bayesian analysis with global sensitivity analysis in order to give better informed predictions; to point out weaker parts of the model that are important targets for further experiments, as well as give guidance on parameters that are essential in distinguishing different qualitative output behaviours. Results We used approximate Bayesian computation (ABC) to estimate the model parameters from experimental data, as well as to quantify the uncertainty in this estimation (inverse uncertainty quantification), resulting in a posterior distribution for the parameters. This parameter uncertainty was next propagated to a corresponding uncertainty in the predictions (forward uncertainty propagation), and a global sensitivity analysis was performed on the prediction using the posterior distribution as the possible values for the parameters. This methodology was applied on a relatively large and complex model relevant for synaptic plasticity, using experimental data from several sources. We could hereby point out those parameters that by themselves have the largest contribution to the uncertainty of the prediction as well as identify parameters important to separate between qualitatively different predictions.This approach is useful both for experimental design as well as model building.

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.

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.

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.

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.

We present a visual analytics system that supports an uncertainty-aware analysis of static and dynamic attributes of biochemical reaction networks (BRNs). These are often described by mathematical models, such as ordinary differential equations (ODEs), which enable the integration of a multitude of different data and data types using parameter estimation. Due to the limited amount of data, parameter estimation does not necessarily yield a single point in parameter space and many attributes of the model remain uncertain. Our system visualizes the model as a graph, where the statistics of the attributes are mapped to the color of edges and vertices. The graph view is combined with several linked views such as lineplots, scatterplots, and correlation matrices, to support the identification of uncertainties and the analysis of their mutual dependencies as well as their time dependencies. To assess the utility of the individual visualization approaches and multiple linked views, a qualitative user study with domain experts was performed. We found that all users were able to process analysis tasks using our system.

Biological measurements of intracellular regulation processes are typically noisy, and time resolution is low. In practice often only steady state measurements of perturbation experiments are available. Since data acquisition is expensive, a framework for experimental design that allows the inclusion of prior knowledge and takes uncertainty into account is highly desirable. We introduce a framework for the experimental design problem to infer parameters from steady state observations of intracellular networks. Our network model consists of (nonlinear) ordinary differential equations based on chemical reaction kinetics. We consider sets of structural perturbation experiments, that is, steady state measurements of the system subject to gene knockout or mutations. The model is stochastically embedded by introducing Gaussian measurement errors. This allows the application of a statistical Bayesian framework and usage of information-theoretic measures for experimental design. We propose to choose the optimal experiments with respect to identifiability of model parameters by maximizing the information content of the expected outcome, measured as the entropy of the posterior distributions. In this setting the posterior has no closed form and an analysis requires efficient sampling methods. We introduce a simulation-based experimental design framework for the identification of network parameters with an efficient entropy estimation approach. First results are shown on a network model for secretory pathway control. Secretion of proteins from cells involves the budding of vesicles at the Golgi. For this process PKD activity is central.

In this paper we consider the problem of parameter estimation for intracellular network models with statistical Bayesian approaches. We use systems of nonlinear differential equations in order to describe the dynamics of those networks. In this setting, the posterior distribution has to be investigated via Markov chain Monte Carlo sampling. An estimation of summary statistics of the posterior from these samples requires appropriate density estimation methods. We focus in this study particularly on the influence of kernel density estimators on the expected information content of the posterior. A new method for the calculation of this information content is introduced that uses directly the unnormalized posterior values at the sample points. We exemplarily show its superiority to kernel estimators on a model of secretory pathway control at the trans-Golgi network in mammalian cells.

We use a discrete‐time Markov process modeling approach to describe the stochastic dynamics of intracellular networks and focus on the inverse problem of network inference using experimental data. In addition to the intrinsic noise of the system, data are corrupted by measurement errors. Thus, the model describes two superimposed stochastic processes with different variances. We introduce a statistical framework for the separation of these two processes via analyzing the marginal likelihood function and apply this concept to small network examples.