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

Monte Carlo samplers for efficient network inference

Accessing information on an underlying network driving a biological process often involves interrupting the process and collecting snapshot data. When snapshot data are stochastic, the data's structure necessitates a probabilistic description to infer underlying reaction networks. As an example, we may imagine wanting to learn gene state networks from the type of data collected in single molecule RNA fluorescence in situ hybridization (RNA-FISH). In the networks we consider, nodes represent network states, and edges represent biochemical reaction rates linking states. Simultaneously estimating the number of nodes and constituent parameters from snapshot data remains a challenging task in part on account of data uncertainty and timescale separations between kinetic parameters mediating the network. While parametric Bayesian methods learn parameters given a network structure (with known node numbers) with rigorously propagated measurement uncertainty, learning the number of nodes and parameters with potentially large timescale separations remain open questions. Here, we propose a Bayesian nonparametric framework and describe a hybrid Bayesian Markov Chain Monte Carlo (MCMC) sampler directly addressing these challenges. In particular, in our hybrid method, Hamiltonian Monte Carlo (HMC) leverages local posterior geometries in inference to explore the parameter space; Adaptive Metropolis Hastings (AMH) learns correlations between plausible parameter sets to efficiently propose probable models; and Parallel Tempering takes into account multiple models simultaneously with tempered information content to augment sampling efficiency. We apply our method to synthetic data mimicking single molecule RNA-FISH, a popular snapshot method in probing transcriptional networks to illustrate the identified challenges inherent to learning dynamical models from these snapshots and how our method addresses them.

A quantitative model for human neurovascular coupling with translated mechanisms from animals

Neurons regulate the activity of blood vessels through the neurovascular coupling (NVC). A detailed understanding of the NVC is critical for understanding data from functional imaging techniques of the brain. Many aspects of the NVC have been studied both experimentally and using mathematical models; various combinations of blood volume and flow, local field potential (LFP), hemoglobin level, blood oxygenation level-dependent response (BOLD), and optogenetics have been measured and modeled in rodents, primates, or humans. However, these data have not been brought together into a unified quantitative model. We now present a mathematical model that describes all such data types and that preserves mechanistic behaviors between experiments. For instance, from modeling of optogenetics and microscopy data in mice, we learn cell-specific contributions; the first rapid dilation in the vascular response is caused by NO-interneurons, the main part of the dilation during longer stimuli is caused by pyramidal neurons, and the post-peak undershoot is caused by NPY-interneurons. These insights are translated and preserved in all subsequent analyses, together with other insights regarding hemoglobin dynamics and the LFP/BOLD-interplay, obtained from other experiments on rodents and primates. The model can predict independent validation-data not used for training. By bringing together data with complementary information from different species, we both understand each dataset better, and have a basis for a new type of integrative analysis of human data.

Fides: Reliable trust-region optimization for parameter estimation of ordinary differential equation models

Ordinary differential equation (ODE) models are widely used to study biochemical reactions in cellular networks since they effectively describe the temporal evolution of these networks using mass action kinetics. The parameters of these models are rarely known a priori and must instead be estimated by calibration using experimental data. Optimization-based calibration of ODE models on is often challenging, even for low-dimensional problems. Multiple hypotheses have been advanced to explain why biochemical model calibration is challenging, including non-identifiability of model parameters, but there are few comprehensive studies that test these hypotheses, likely because tools for performing such studies are also lacking. Nonetheless, reliable model calibration is essential for uncertainty analysis, model comparison, and biological interpretation.

We implemented an established trust-region method as a modular Python framework (fides) to enable systematic comparison of different approaches to ODE model calibration involving a variety of Hessian approximation schemes. We evaluated fides on a recently developed corpus of biologically realistic benchmark problems for which real experimental data are available. Unexpectedly, we observed high variability in optimizer performance among different implementations of the same mathematical instructions (algorithms). Analysis of possible sources of poor optimizer performance identified limitations in the widely used Gauss-Newton, BFGS and SR1 Hessian approximation schemes. We addressed these drawbacks with a novel hybrid Hessian approximation scheme that enhances optimizer performance and outperforms existing hybrid approaches. When applied to the corpus of test models, we found that fides was on average more reliable and efficient than existing methods using a variety of criteria. We expect fides to be broadly useful for ODE constrained optimization problems in biochemical models and to be a foundation for future methods development.

Adaptive Rejection of a Sinusoidal Disturbance with Unknown Frequency in a Flexible Rotor with Lubricated Journal Bearings

Frequency-dependent adaptive noise cancellation-based tracking controller (ANC-TC) is a known technique for the stabilization of several nonlinear dynamical systems. In recent years, this control strategy has been introduced and applied for the stabilization of a flexible rotor supported on full-lubricated journal bearings. This paper proposes a theoretical investigation, based on robust immersion and invariance (I&I) theory, of a novel ANC-frequency estimation (FE) technique designed to stabilize a flexible rotor shaft affected by self-generated sinusoidal disturbances, generalized to the case of unknown frequency. A structural proof, under assumptions on closed-loop output signals, shows that the sinusoidal disturbance rejection is exponential. Numerical simulations are presented to validate the mathematical results in silico. The iterative Inexact Newton method is applied to the disturbance frequency and phase estimation error point series. The data fitting confirms that the phase estimation succession has an exponential convergence behavior and that the asymptotical frequency estimation is a warm-up phase of the overall close-loop disturbance estimation process. In two different operating conditions, the orders of convergence obtained by phase and frequency estimate timeseries are pφ=1, pω,unc=0.9983 and pω,cav=1.005. Rejection of the rotor dynamic disturbance occurs approximately 76% before in the cavitated than in the uncavitated condition, 2 (s) and 8.5 (s), respectively.

A protocol for dynamic model calibration

Ordinary differential equation models are nowadays widely used for the mechanistic description of biological processes and their temporal evolution. These models typically have many unknown and nonmeasurable parameters, which have to be determined by fitting the model to experimental data. In order to perform this task, known as parameter estimation or model calibration, the modeller faces challenges such as poor parameter identifiability, lack of sufficiently informative experimental data and the existence of local minima in the objective function landscape. These issues tend to worsen with larger model sizes, increasing the computational complexity and the number of unknown parameters. An incorrectly calibrated model is problematic because it may result in inaccurate predictions and misleading conclusions. For nonexpert users, there are a large number of potential pitfalls. Here, we provide a protocol that guides the user through all the steps involved in the calibration of dynamic models. We illustrate the methodology with two models and provide all the code required to reproduce the results and perform the same analysis on new models. Our protocol provides practitioners and researchers in biological modelling with a one-stop guide that is at the same time compact and sufficiently comprehensive to cover all aspects of the problem.

Accelerated regression-based summary statistics for discrete stochastic systems via approximate simulators

BACKGROUND

Approximate Bayesian Computation (ABC) has become a key tool for calibrating the parameters of discrete stochastic biochemical models. For higher dimensional models and data, its performance is strongly dependent on having a representative set of summary statistics. While regression-based methods have been demonstrated to allow for the automatic construction of effective summary statistics, their reliance on first simulating a large training set creates a significant overhead when applying these methods to discrete stochastic models for which simulation is relatively expensive. In this τ work, we present a method to reduce this computational burden by leveraging approximate simulators of these systems, such as ordinary differential equations and τ-Leaping approximations.

RESULTS

We have developed an algorithm to accelerate the construction of regression-based summary statistics for Approximate Bayesian Computation by selectively using the faster approximate algorithms for simulations. By posing the problem as one of ratio estimation, we use state-of-the-art methods in machine learning to show that, in many cases, our algorithm can significantly reduce the number of simulations from the full resolution model at a minimal cost to accuracy and little additional tuning from the user. We demonstrate the usefulness and robustness of our method with four different experiments.

CONCLUSIONS

We provide a novel algorithm for accelerating the construction of summary statistics for stochastic biochemical systems. Compared to the standard practice of exclusively training from exact simulator samples, our method is able to dramatically reduce the number of required calls to the stochastic simulator at a minimal loss in accuracy. This can immediately be implemented to increase the overall speed of the ABC workflow for estimating parameters in complex systems.

High-speed parameter search of dynamic biological pathways from time-course transcriptomic profiles using high-level Petri net

Dynamic simulation promises a deeper understanding of complex molecular mechanisms of biological pathways. How to determine the reaction kinetic parameters which govern the simulation results is still an open question in the field of systems biology. (1) Background: To execute simulation experiments, it is an essential first step to search effective values of model parameters. The complexity of biological systems and the experimental measurement technology severely limit the acquirement of accurate kinetic parameters. Previously proposed genomic data assimilation (GDA) approach enables users to handle parameter estimation using time-course information. However, it highly depends on successive time points and costs massive computational resource; (2) Methods: To address this problem, we present a new high-speed parameter search method for estimating the kinetic parameters of quantitative biological pathways using time-course transcriptomic profiles. The key idea of our method is to interactively prune the search space by introducing Probabilistic Linear-time Temporal Logic (PLTL) based model checking into GDA. (3) Results and conclusion: We demonstrated the effectiveness of our method by comparing with GDA on Mus musculus transcription circuits modelled by hybrid functional Petri net with extension. As a result, our method works faster and more accurate than GDA for both time-course datasets with dense and sparse observed values.

Reoccurring neural stem cell divisions in the adult zebrafish telencephalon are sufficient for the emergence of aggregated spatiotemporal patterns

Regulation of quiescence and cell cycle entry is pivotal for the maintenance of stem cell populations. Regulatory mechanisms, however, are poorly understood. In particular, it is unclear how the activity of single stem cells is coordinated within the population or if cells divide in a purely random fashion. We addressed this issue by analyzing division events in an adult neural stem cell (NSC) population of the zebrafish telencephalon. Spatial statistics and mathematical modeling of over 80,000 NSCs in 36 brain hemispheres revealed weakly aggregated, nonrandom division patterns in space and time. Analyzing divisions at 2 time points allowed us to infer cell cycle and S-phase lengths computationally. Interestingly, we observed rapid cell cycle reentries in roughly 15% of newly born NSCs. In agent-based simulations of NSC populations, this redividing activity sufficed to induce aggregated spatiotemporal division patterns that matched the ones observed experimentally. In contrast, omitting redivisions leads to a random spatiotemporal distribution of dividing cells. Spatiotemporal aggregation of dividing stem cells can thus emerge solely from the cells’ history.

An interdisciplinary study of the rules governing cell divisions in a population of neural stem cells in the zebrafish brain reveals the existence of aggregated spatio-temporal division patterns of rapid cell cycles in stem cells, and shows that these patterns can be explained by a simple agent-based model relying solely on the cells‘ division history.