Inference for discretely observed stochastic kinetic networks with applications to epidemic modeling

Statistical Inference of Rate Constants in Chemical and Biochemical Reaction Networks Using an "Inverse" Event-Driven Kinetic Monte Carlo Method

The use of rate models for networks of stochastic reactions is frequently used to comprehend the macroscopically observed dynamic properties of finite size reactive systems as well as their relationship to the underlying molecular events. Τhis particular approach usually stumbles on parameter derivation associated with stochastic kinetics, a quite demanding procedure. The present study incorporates a novel algorithm, which infers kinetic parameters from the system's time evolution, manifested as changes in molecular species populations. The proposed methodology reconstructs distributions required to infer kinetic parameters of a stochastic process pertaining to either a simulation or experimental data. The suggested approach accurately replicates rate constants of the stochastic reaction networks, which have evolved over time by event-driven Monte Carlo (MC) simulations using the Gillespie algorithm. Furthermore, our approach has been successfully used to estimate rate constants of association and dissociation events between molecular species developing during molecular dynamics (MD) simulations. We certainly believe that our method will be remarkably helpful for considering the macroscopic characteristic molecular roots related to stochastic physical and biological processes.

The multilevel hierarchical data EM-algorithm. Applications to discrete-time Markov chain epidemic models

The theory of multilevel hierarchical data Expectation Maximization (EM)-algorithm is introduced via discrete time Markov chain (DTMC) epidemic models. A general model for a multilevel hierarchical discrete data is derived. The observed sample Y in the system is a stochastic incomplete data, and the missing data Z exhibits a multilevel hierarchical data structure. The EM-algorithm to find ML-estimates for parameters in the stochastic system is derived. Applications of the EM-algorithm are exhibited in the two DTMC models, to find ML-estimates of the system parameters. Numerical results are given for influenza epidemics in the state of Georgia (GA), USA.

A surrogate Bayesian framework for a SARS-CoV-2 data driven stochastic model

Dynamic compartmentalized data (DCD) and compartmentalized differential equations (CDEs) are key instruments for modeling transmission of pathogens such as the SARS-CoV-2 virus. We describe an effi-cient nowcasting algorithm for modeling transmission of SARS-CoV-2 with uncertainty quantification for the COVID-19 impact. A key concern for transmission of SARS-CoV-2 is under-reporting of cases, and this is addressed in our data-driven model by providing an estimate for the detection rate. Our novel top-down model is based on CDEs with stochastic constitutive parameters obtained from the DCD using Bayesian inference. We demonstrate the robustness of our algorithm for simulation studies using synthetic DCD, and nowcasting COVID-19 using real DCD from several regions across five continents.

COVID-19: Estimation of the transmission dynamics in Spain using a stochastic simulator and black-box optimization techniques

BACKGROUND AND OBJECTIVES

Epidemiological models of epidemic spread are an essential tool for optimizing decision-making. The current literature is very extensive and covers a wide variety of deterministic and stochastic models. However, with the increase in computing resources, new, more general, and flexible procedures based on simulation models can assess the effectiveness of measures and quantify the current state of the epidemic. This paper illustrates the potential of this approach to build a new dynamic probabilistic model to estimate the prevalence of SARS-CoV-2 infections in different compartments.

METHODS

We propose a new probabilistic model in which, for the first time in the epidemic literature, parameter learning is carried out using gradient-free stochastic black-box optimization techniques simulating multiple trajectories of the infection dynamics in a general way, solving an inverse problem that is defined employing the daily information from mortality records.

RESULTS

After the application of the new proposal in Spain in the first and successive waves, the result of the model confirms the accuracy to estimate the seroprevalence and allows us to know the real dynamics of the pandemic a posteriori to assess the impact of epidemiological measures by the Spanish government and to plan more efficiently the subsequent decisions with the prior knowledge obtained.

CONCLUSIONS

The model results allow us to estimate the daily patterns of COVID-19 infections in Spain retrospectively and examine the population's exposure to the virus dynamically in contrast to seroprevalence surveys. Furthermore, given the flexibility of our simulation framework, we can model situations -even using non-parametric distributions between the different compartments in the model- that other models in the existing literature cannot. Our general optimization strategy remains valid in these cases, and we can easily create other non-standard simulation epidemic models that incorporate more complex and dynamic structures.

Network reconstruction from infection cascades

Accessing the network through which a propagation dynamics diffuses is essential for understanding and controlling it. In a few cases, such information is available through direct experiments or thanks to the very nature of propagation data. In a majority of cases however, available information about the network is indirect and comes from partial observations of the dynamics, rendering the network reconstruction a fundamental inverse problem. Here we show that it is possible to reconstruct the whole structure of an interaction network and to simultaneously infer the complete time course of activation spreading, relying just on single epoch (i.e. snapshot) or time-scattered observations of a small number of activity cascades. The method that we present is built on a belief propagation approximation, that has shown impressive accuracy in a wide variety of relevant cases, and is able to infer interactions in the presence of incomplete time-series data by providing a detailed modelling of the posterior distribution of trajectories conditioned to the observations. Furthermore, we show by experiments that the information content of full cascades is relatively smaller than that of sparse observations or single snapshots.

Reverse engineering gene networks using global-local shrinkage rules

Inferring gene regulatory networks from high-throughput 'omics' data has proven to be a computationally demanding task of critical importance. Frequently, the classical methods break down owing to the curse of dimensionality, and popular strategies to overcome this are typically based on regularized versions of the classical methods. However, these approaches rely on loss functions that may not be robust and usually do not allow for the incorporation of prior information in a straightforward way. Fully Bayesian methods are equipped to handle both of these shortcomings quite naturally, and they offer the potential for improvements in network structure learning. We propose a Bayesian hierarchical model to reconstruct gene regulatory networks from time-series gene expression data, such as those common in perturbation experiments of biological systems. The proposed methodology uses global-local shrinkage priors for posterior selection of regulatory edges and relaxes the common normal likelihood assumption in order to allow for heavy-tailed data, which were shown in several of the cited references to severely impact network inference. We provide a sufficient condition for posterior propriety and derive an efficient Markov chain Monte Carlo via Gibbs sampling in the electronic supplementary material. We describe a novel way to detect multiple scales based on the corresponding posterior quantities. Finally, we demonstrate the performance of our approach in a simulation study and compare it with existing methods on real data from a T-cell activation study.

Survival dynamical systems: individual-level survival analysis from population-level epidemic models

In this paper, we show that solutions to ordinary differential equations describing the large-population limits of Markovian stochastic epidemic models can be interpreted as survival or cumulative hazard functions when analysing data on individuals sampled from the population. We refer to the individual-level survival and hazard functions derived from population-level equations as a survival dynamical system (SDS). To illustrate how population-level dynamics imply probability laws for individual-level infection and recovery times that can be used for statistical inference, we show numerical examples based on synthetic data. In these examples, we show that an SDS analysis compares favourably with a complete-data maximum-likelihood analysis. Finally, we use the SDS approach to analyse data from a 2009 influenza A(H1N1) outbreak at Washington State University.

Bayesian inference of distributed time delay in transcriptional and translational regulation

MOTIVATION

Advances in experimental and imaging techniques have allowed for unprecedented insights into the dynamical processes within individual cells. However, many facets of intracellular dynamics remain hidden, or can be measured only indirectly. This makes it challenging to reconstruct the regulatory networks that govern the biochemical processes underlying various cell functions. Current estimation techniques for inferring reaction rates frequently rely on marginalization over unobserved processes and states. Even in simple systems this approach can be computationally challenging, and can lead to large uncertainties and lack of robustness in parameter estimates. Therefore we will require alternative approaches to efficiently uncover the interactions in complex biochemical networks.

RESULTS

We propose a Bayesian inference framework based on replacing uninteresting or unobserved reactions with time delays. Although the resulting models are non-Markovian, recent results on stochastic systems with random delays allow us to rigorously obtain expressions for the likelihoods of model parameters. In turn, this allows us to extend MCMC methods to efficiently estimate reaction rates, and delay distribution parameters, from single-cell assays. We illustrate the advantages, and potential pitfalls, of the approach using a birth-death model with both synthetic and experimental data, and show that we can robustly infer model parameters using a relatively small number of measurements. We demonstrate how to do so even when only the relative molecule count within the cell is measured, as in the case of fluorescence microscopy.

AVAILABILITY AND IMPLEMENTATION

Accompanying code in R is available at https://github.com/cbskust/DDE_BD.

SUPPLEMENTARY INFORMATION

Supplementary data are available at Bioinformatics online.

Modeling epidemics: A primer and Numerus Model Builder implementation

Epidemiological models are dominated by compartmental models, of which SIR formulations are the most commonly used. These formulations can be continuous or discrete (in either the state-variable values or time), deterministic or stochastic, or spatially homogeneous or heterogeneous, the latter often embracing a network formulation. Here we review the continuous and discrete deterministic and discrete stochastic formulations of the SIR dynamical systems models, and we outline how they can be easily and rapidly constructed using Numerus Model Builder, a graphically-driven coding platform. We also demonstrate how to extend these models to a metapopulation setting using NMB network and mapping tools.

Beyond the Michaelis-Menten equation: Accurate and efficient estimation of enzyme kinetic parameters

Examining enzyme kinetics is critical for understanding cellular systems and for using enzymes in industry. The Michaelis-Menten equation has been widely used for over a century to estimate the enzyme kinetic parameters from reaction progress curves of substrates, which is known as the progress curve assay. However, this canonical approach works in limited conditions, such as when there is a large excess of substrate over enzyme. Even when this condition is satisfied, the identifiability of parameters is not always guaranteed, and often not verifiable in practice. To overcome such limitations of the canonical approach for the progress curve assay, here we propose a Bayesian approach based on an equation derived with the total quasi-steady-state approximation. In contrast to the canonical approach, estimates obtained with this proposed approach exhibit little bias for any combination of enzyme and substrate concentrations. Importantly, unlike the canonical approach, an optimal experiment to identify parameters with certainty can be easily designed without any prior information. Indeed, with this proposed design, the kinetic parameters of diverse enzymes with disparate catalytic efficiencies, such as chymotrypsin, fumarase, and urease, can be accurately and precisely estimated from a minimal amount of timecourse data. A publicly accessible computational package performing such accurate and efficient Bayesian inference for enzyme kinetics is provided.

Network-based analysis of a small Ebola outbreak

We present a method for estimating epidemic parameters in network-based stochastic epidemic models when the total number of infections is assumed to be small. We illustrate the method by reanalyzing the data from the 2014 Democratic Republic of the Congo (DRC) Ebola outbreak described in Maganga et al. (2014).

A Likelihood Approach for Real-Time Calibration of Stochastic Compartmental Epidemic Models

Stochastic transmission dynamic models are especially useful for studying the early emergence of novel pathogens given the importance of chance events when the number of infectious individuals is small. However, methods for parameter estimation and prediction for these types of stochastic models remain limited. In this manuscript, we describe a calibration and prediction framework for stochastic compartmental transmission models of epidemics. The proposed method, Multiple Shooting for Stochastic systems (MSS), applies a linear noise approximation to describe the size of the fluctuations, and uses each new surveillance observation to update the belief about the true epidemic state. Using simulated outbreaks of a novel viral pathogen, we evaluate the accuracy of MSS for real-time parameter estimation and prediction during epidemics. We assume that weekly counts for the number of new diagnosed cases are available and serve as an imperfect proxy of incidence. We show that MSS produces accurate estimates of key epidemic parameters (i.e. mean duration of infectiousness, R0, and Reff) and can provide an accurate estimate of the unobserved number of infectious individuals during the course of an epidemic. MSS also allows for accurate prediction of the number and timing of future hospitalizations and the overall attack rate. We compare the performance of MSS to three state-of-the-art benchmark methods: 1) a likelihood approximation with an assumption of independent Poisson observations; 2) a particle filtering method; and 3) an ensemble Kalman filter method. We find that MSS significantly outperforms each of these three benchmark methods in the majority of epidemic scenarios tested. In summary, MSS is a promising method that may improve on current approaches for calibration and prediction using stochastic models of epidemics.

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.

A stochastic transcriptional switch model for single cell imaging data

Gene expression is made up of inherently stochastic processes within single cells and can be modeled through stochastic reaction networks (SRNs). In particular, SRNs capture the features of intrinsic variability arising from intracellular biochemical processes. We extend current models for gene expression to allow the transcriptional process within an SRN to follow a random step or switch function which may be estimated using reversible jump Markov chain Monte Carlo (MCMC). This stochastic switch model provides a generic framework to capture many different dynamic features observed in single cell gene expression. Inference for such SRNs is challenging due to the intractability of the transition densities. We derive a model-specific birth–death approximation and study its use for inference in comparison with the linear noise approximation where both approximations are considered within the unifying framework of state-space models. The methodology is applied to synthetic as well as experimental single cell imaging data measuring expression of the human prolactin gene in pituitary cells.

Comparison of Parameter Estimation Methods in Stochastic Chemical Kinetic Models: Examples in Systems Biology

Stochastic chemical kinetics has become a staple for mechanistically modeling various phenomena in systems biology. These models, even more so than their deterministic counterparts, pose a challenging problem in the estimation of kinetic parameters from experimental data. As a result of the inherent randomness involved in stochastic chemical kinetic models, the estimation methods tend to be statistical in nature. Three classes of estimation methods are implemented and compared in this paper. The first is the exact method, which uses the continuous-time Markov chain representation of stochastic chemical kinetics and is tractable only for a very restricted class of problems. The next class of methods is based on Markov chain Monte Carlo (MCMC) techniques. The third method, termed conditional density importance sampling (CDIS), is a new method introduced in this paper. The use of these methods is demonstrated on two examples taken from systems biology, one of which is a new model of single-cell viral infection. The applicability, strengths and weaknesses of the three classes of estimation methods are discussed. Using simulated data for the two examples, some guidelines are provided on experimental design to obtain more information from a limited number of measurements.

Least squares estimation in stochastic biochemical networks

The paper presents results on the asymptotic properties of the least-squares estimates (LSEs) of the reaction constants in mass-action, stochastic, biochemical network models. LSEs are assumed to be based on the longitudinal data from partially observed trajectories of a stochastic dynamical system, modeled as a continuous-time, pure jump Markov process. Under certain regularity conditions on such a process, it is shown that the vector of LSEs is jointly consistent and asymptotically normal, with the asymptotic covariance structure given in terms of a system of ordinary differential equations (ODE). The derived asymptotic properties hold true as the biochemical network size (the total species number) increases, in which case the stochastic dynamical system converges to the deterministic mass-action ODE. An example is provided, based on synthetic as well as RT-PCR data from the retro-transcription network of the LINE1 gene.