Identification of dynamic mass-action biochemical reaction networks using sparse Bayesian methods

Identifying the reactions that govern a dynamical biological system is a crucial but challenging task in systems biology. In this work, we present a data-driven method to infer the underlying biochemical reaction system governing a set of observed species concentrations over time. We formulate the problem as a regression over a large, but limited, mass-action constrained reaction space and utilize sparse Bayesian inference via the regularized horseshoe prior to produce robust, interpretable biochemical reaction networks, along with uncertainty estimates of parameters. The resulting systems of chemical reactions and posteriors inform the biologist of potentially several reaction systems that can be further investigated. We demonstrate the method on two examples of recovering the dynamics of an unknown reaction system, to illustrate the benefits of improved accuracy and information obtained.

Direct statistical inference for finite Markov jump processes via the matrix exponential

Given noisy, partial observations of a time-homogeneous, finite-statespace Markov chain, conceptually simple, direct statistical inference is available, in theory, via its rate matrix, or infinitesimal generator, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathsf {Q}}$$\end{document}Q, since \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\exp ({\mathsf {Q}}t)$$\end{document}exp(Qt) is the transition matrix over time t. However, perhaps because of inadequate tools for matrix exponentiation in programming languages commonly used amongst statisticians or a belief that the necessary calculations are prohibitively expensive, statistical inference for continuous-time Markov chains with a large but finite state space is typically conducted via particle MCMC or other relatively complex inference schemes. When, as in many applications \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathsf {Q}}$$\end{document}Q arises from a reaction network, it is usually sparse. We describe variations on known algorithms which allow fast, robust and accurate evaluation of the product of a non-negative vector with the exponential of a large, sparse rate matrix. Our implementation uses relatively recently developed, efficient, linear algebra tools that take advantage of such sparsity. We demonstrate the straightforward statistical application of the key algorithm on a model for the mixing of two alleles in a population and on the Susceptible-Infectious-Removed epidemic model.

Likelihood-free nested sampling for parameter inference of biochemical reaction networks

The development of mechanistic models of biological systems is a central part of Systems Biology. One major challenge in developing these models is the accurate inference of model parameters. In recent years, nested sampling methods have gained increased attention in the Systems Biology community due to the fact that they are parallelizable and provide error estimates with no additional computations. One drawback that severely limits the usability of these methods, however, is that they require the likelihood function to be available, and thus cannot be applied to systems with intractable likelihoods, such as stochastic models. Here we present a likelihood-free nested sampling method for parameter inference which overcomes these drawbacks. This method gives an unbiased estimator of the Bayesian evidence as well as samples from the posterior. We derive a lower bound on the estimators variance which we use to formulate a novel termination criterion for nested sampling. The presented method enables not only the reliable inference of the posterior of parameters for stochastic systems of a size and complexity that is challenging for traditional methods, but it also provides an estimate of the obtained variance. We illustrate our approach by applying it to several realistically sized models with simulated data as well as recently published biological data. We also compare our developed method with the two most popular other likeliood-free approaches: pMCMC and ABC-SMC. The C++ code of the proposed methods, together with test data, is available at the github web page https://github.com/Mijan/LFNS_paper.

The behaviour of mathematical models of biochemical reactions is governed by model parameters encoding for various reaction rates, molecule concentrations and other biochemical quantities. As the general purpose of these models is to reproduce and predict the true biological response to different stimuli, the inference of these parameters, given experimental observations, is a crucial part of Systems Biology. While plenty of methods have been published for the inference of model parameters, most of them require the availability of the likelihood function and thus cannot be applied to models that do not allow for the computation of the likelihood. Further, most established methods do not provide an estimate of the variance of the obtained estimator. In this paper, we present a novel inference method that accurately approximates the posterior distribution of parameters and does not require the evaluation of the likelihood function. Our method is based on the nested sampling algorithm and approximates the likelihood with a particle filter. We show that the resulting posterior estimates are unbiased and provide a way to estimate not just the posterior distribution, but also an error estimate of the final estimator. We illustrate our method on several stochastic models with simulated data as well as one model of transcription with real biological data.

Inferring Mechanistic Parameters from Amyloid Formation Kinetics by Approximate Bayesian Computation

Amyloid formation is implicated in a number of human diseases, and is thought to proceed via a nucleation-dependent polymerization mechanism. Experimenters often wish to relate changes in amyloid formation kinetics, for example, in response to small molecules to specific mechanistic steps along this pathway. However, fitting kinetic fibril formation data to a complex model including explicit rate constants results in an ill-posed problem with a vast number of potential solutions. The levels of uncertainty remaining in parameters calculated from these models, arising both from experimental noise and high levels of degeneracy or codependency in parameters, is often unclear. Here, we demonstrate that a combination of explicit mathematical models with an approximate Bayesian computation approach can be used to assign the mechanistic effects of modulators on amyloid fibril formation. We show that even when exact rate constants cannot be extracted, parameters derived from these rate constants can be recovered and used to assign mechanistic effects and their relative magnitudes with a great deal of confidence. Furthermore, approximate Bayesian computation provides a robust method for visualizing uncertainty remaining in the model parameters, regardless of its origin. We apply these methods to the problem of heparin-mediated tau polymerization, which displays complex kinetic behavior not amenable to analysis by more traditional methods. Our analysis indicates that the role of heparin cannot be explained by enhancement of nucleation alone, as has been previously proposed. The methods described here are applicable to a wide range of systems, as models can be easily adapted to account for new reactions and reversibility.

Bayesian inference of reaction kinetics from single-cell recordings across a heterogeneous cell population

Single-cell experimental techniques provide informative data to help uncover dynamical processes inside a cell. Making full use of such data requires dedicated computational methods to estimate biophysical process parameters and states in a model-based manner. In particular, the treatment of heterogeneity or cell-to-cell variability deserves special attention. The present article provides an introduction to one particular class of algorithms which employ marginalization in order to take heterogeneity into account. An overview of alternative approaches is provided for comparison. We treat two frequently encountered scenarios in single-cell experiments, namely, single-cell trajectory data and single-cell distribution data.