Bayesian inference for dynamic transcriptional regulation; the Hes1 system as a case study

MOTIVATION: In this study, we address the problem of estimating the parameters of regulatory networks and provide the first application of Markov chain Monte Carlo (MCMC) methods to experimental data. As a case study, we consider a stochastic model of the Hes1 system expressed in terms of stochastic differential equations (SDEs) to which rigorous likelihood methods of inference can be applied. When fitting continuous-time stochastic models to discretely observed time series the lengths of the sampling intervals are important, and much of our study addresses the problem when the data are sparse. RESULTS: We estimate the parameters of an autoregulatory network providing results both for simulated and real experimental data from the Hes1 system. We develop an estimation algorithm using MCMC techniques which are flexible enough to allow for the imputation of latent data on a finer time scale and the presence of prior information about parameters which may be informed from other experiments as well as additional measurement error.     

A diffusion process to model generalized von Bertalanffy growth patterns: Fitting to real data

The von Bertalanffy growth curve has been commonly used for modeling animal growth (particularly fish). Both deterministic and stochastic models exist in association with this curve, the latter allowing for the inclusion of fluctuations or disturbances that might exist in the system under consideration which are not always quantifiable or may even be unknown. This curve is mainly used for modeling the length variable whereas a generalized version, including a new parameter b≥1, allows for modeling both length and weight for some animal species in both isometric (b=3) and allometric (b≠3) situations.In this paper a stochastic model related to the generalized von Bertalanffy growth curve is proposed. This model allows to investigate the time evolution of growth variables associated both with individual behaviors and mean population behavior. Also, with the purpose of fitting the above-mentioned model to real data and so be able to forecast and analyze particular characteristics, we study the maximum likelihood estimation of the parameters of the model. In addition, and regarding the numerical problems posed by solving the likelihood equations, a strategy is developed for obtaining initial solutions for the usual numerical procedures. Such strategy is validated by means of simulated examples. Finally, an application to real data of mean weight of swordfish is presented.     

A finite-state continuous-time approach for inferring regional migration and mortality rates from archival tagging and conventional tag-recovery experiments

Spatially structured population dynamics models are important management tools for harvested, highly mobile species and although conventional tag recovery experiments remain useful for estimation of various demographic parameters of these models, archival tagging experiments are becoming an important data source for analyzing migratory behavior of mobile marine species. We provide a likelihood-based approach for estimating the regional migration and mortality rate parameters intrinsic to these models that may use information obtained from conventional tag recovery and archival tagging experiments. Specifically, we assume that the regional location and survival of animals through time is a finite-state continuous-time stochastic process. The stochastic process is the basis of probability models for observations provided by the different types of tags. Results from application to simulated tagging experiments for western Atlantic bluefin tuna show that maximum likelihood estimators based on archival tagging observations and corresponding confidence intervals perform similar to conventional tagging observations for a given number of tag releases and releasing tags in each region can improve the behavior of maximum likelihood estimators regardless of tag type. We provide an example application with Atlantic bluefin tuna released with conventional tags in 1990-1992.     

Longitudinal Data Analysis with Event Time as a Covariate

We consider the estimation of a nonparametric smooth function of some event time in a semiparametric mixed effects model from repeatedly measured data when the event time is subject to right censoring. The within-subject correlation is captured by both cross-sectional and time-dependent random effects, where the latter is modeled by a nonhomogeneous Ornstein–Uhlenbeck stochastic process. When the censoring probability depends on other variables in the model, which often happens in practice, the event time data are not missing completely at random. Hence, the complete case analysis by eliminating all the censored observations may yield biased estimates of the regression parameters including the smooth function of the event time, and is less efficient. To remedy, we derive the likelihood function for the observed data by modeling the event time distribution given other covariates. We propose a two-stage pseudo-likelihood approach for the estimation of model parameters by first plugging an estimator of the conditional event time distribution into the likelihood and then maximizing the resulting pseudo-likelihood function. Empirical evaluation shows that the proposed method yields negligible biases while significantly reduces the estimation variability. This research is motivated by the project of hormone profile estimation around age at the final menstrual period for the cohort of women in the Michigan Bone Health and Metabolism Study.     

Parameter estimation in a Gompertzian stochastic model for tumor growth

The problem of estimating parameters in the drift coefficient when a diffusion process is observed continuously requires some specific assumptions. In this paper, we consider a stochastic version of the Gompertzian model that describes in vivo tumor growth and its sensitivity to treatment with antiangiogenic drugs. An explicit likelihood function is obtained, and we discuss some properties of the maximum likelihood estimator for the intrinsic growth rate of the stochastic Gompertzian model. Furthermore, we show some simulation results on the behavior of the corresponding discrete estimator. Finally, an application is given to illustrate the estimate of the model parameters using real data.     

Stochastic analogues of deterministic single-species population models

Although single-species deterministic difference equations have long been used in modeling the dynamics of animal populations, little attention has been paid to how stochasticity should be incorporated into these models. By deriving stochastic analogues to difference equations from first principles, we show that the form of these models depends on whether noise in the population process is demographic or environmental. When noise is demographic, we argue that variance around the expectation is proportional to the expectation. When noise is environmental the variance depends in a non-trivial way on how variation enters into model parameters, but we argue that if the environment affects the population multiplicatively then variance is proportional to the square of the expectation. We compare various stochastic analogues of the Ricker map model by fitting them, using maximum likelihood estimation, to data generated from an individual-based model and the weevil data of Utida. Our demographic models are significantly better than our environmental models at fitting noise generated by population processes where noise is mainly demographic. However, the traditionally chosen stochastic analogues to deterministic models--additive normally distributed noise and multiplicative lognormally distributed noise--generally fit all data sets well. Thus, the form of the variance does play a role in the fitting of models to ecological time series, but may not be important in practice as first supposed.     

A Stochastic Compartmental Model with Long Lasting Infusion

Most of the compartmental models in current use to model pharmacokinetic systems are deterministic. Stochastic formulations of pharmacokinetic compartmental models introduce stochasticity through either a probabilistic transfer mechanism or the randomization of the transfer rate constants. In this paper we consider a linear stochastic differential equation (LSDE) which represents a stochastic version of a one‐compartment linear model when input function undergoes random fluctuations. The solution of the LSDE, its mean value and covariance structure are derived. An explicit likelihood function is obtained either when the process is observed continuously over a period of time or when sampled data are available, as it is generally feasible. We discuss some asymptotic properties of the maximum likelihood estimators for the model parameters. Furthermore we develop expressions for two random variables of interest in pharmacokinetics: the area under the time‐concentration curve, M₀(T), and the plateau concentration, x_ss. Finally the estimation procedure is illustrated by an application to real data.     

Toward a stochastic formulation of microbial growth in relation to bioreactor performances: case study of an E. coli fed-batch process

A stochastic microbial growth model has been elaborated in the case of the culture of E. coli in fed-batch and scale-down reactors. This model is based on the stochastic determination of the generation time of the microbial cells. The determination of generation time is determined by choosing the appropriate value on a log-normal distribution. The appropriateness of such distribution is discussed and growth curves are obtained that show good agreement compared with the experimental results. The mean and the standard deviation of the log-normal distribution can be considered to be constant during the batch phase of the culture, but they vary when the fed-batch mode is started. It has been shown that the parameters related to the log-normal distribution are submitted to an exponential evolution. The aim of this study is to explore the bioreactor hydrodynamic effect on microbial growth. Thus, in a second time, the stochastic growth model has been reinforced by data coming from a previous stochastic bioreactor mixing model (1). The connection of these hydrodynamic data with the actual stochastic growth model has allowed us to explain the scale-down effect associated with the glucose concentration fluctuations. It is important to point out that the scale-down effect is induced differently according to the feeding strategy involved in the fed-batch experiments.     

Demography of Verreaux’s sifaka in a stochastic rainfall environment

In this study, we use deterministic and stochastic models to analyze the demography of Verreaux’s sifaka (Propithecus verreauxi verreauxi) in a fluctuating rainfall environment. The model is based on 16 years of data from Beza Mahafaly Special Reserve, southwest Madagascar. The parameters in the stage-classified life cycle were estimated using mark-recapture methods. Statistical models were evaluated using information-theoretic techniques and multi-model inference. The highest ranking model is time-invariant, but the averaged model includes rainfall-dependence of survival and breeding. We used a time-series model of rainfall to construct a stochastic demographic model. The time-invariant model and the stochastic model give a population growth rate of about 0.98. Bootstrap confidence intervals on the growth rates, both deterministic and stochastic, include 1. Growth rates are most elastic to changes in adult survival. Many demographic statistics show a nonlinear response to annual rainfall but are depressed when annual rainfall is low, or the variance in annual rainfall is high. Perturbation analyses from both the time-invariant and stochastic models indicate that recruitment and survival of older females are key determinants of population growth rate.     

Inference of Epidemiological Dynamics Based on Simulated Phylogenies Using Birth-Death and Coalescent Models

Quantifying epidemiological dynamics is crucial for understanding and forecasting the spread of an epidemic. The coalescent and the birth-death model are used interchangeably to infer epidemiological parameters from the genealogical relationships of the pathogen population under study, which in turn are inferred from the pathogen genetic sequencing data. To compare the performance of these widely applied models, we performed a simulation study. We simulated phylogenetic trees under the constant rate birth-death model and the coalescent model with a deterministic exponentially growing infected population. For each tree, we re-estimated the epidemiological parameters using both a birth-death and a coalescent based method, implemented as an MCMC procedure in BEAST v2.0. In our analyses that estimate the growth rate of an epidemic based on simulated birth-death trees, the point estimates such as the maximum a posteriori/maximum likelihood estimates are not very different. However, the estimates of uncertainty are very different. The birth-death model had a higher coverage than the coalescent model, i.e. contained the true value in the highest posterior density (HPD) interval more often (2–13% vs. 31–75% error). The coverage of the coalescent decreases with decreasing basic reproductive ratio and increasing sampling probability of infecteds. We hypothesize that the biases in the coalescent are due to the assumption of deterministic rather than stochastic population size changes. Both methods performed reasonably well when analyzing trees simulated under the coalescent. The methods can also identify other key epidemiological parameters as long as one of the parameters is fixed to its true value. In summary, when using genetic data to estimate epidemic dynamics, our results suggest that the birth-death method will be less sensitive to population fluctuations of early outbreaks than the coalescent method that assumes a deterministic exponentially growing infected population.     

Identifiability of parameters in MCMC Bayesian inference of phylogeny

Methods for Bayesian inference of phylogeny using DNA sequences based on Markov chain Monte Carlo (MCMC) techniques allow the incorporation of arbitrarily complex models of the DNA substitution process, and other aspects of evolution. This has increased the realism of models, potentially improving the accuracy of the methods, and is largely responsible for their recent popularity. Another consequence of the increased complexity of models in Bayesian phylogenetics is that these models have, in several cases, become overparameterized. In such cases, some parameters of the model are not identifiable; different combinations of nonidentifiable parameters lead to the same likelihood, making it impossible to decide among the potential parameter values based on the data. Overparameterized models can also slow the rate of convergence of MCMC algorithms due to large negative correlations among parameters in the posterior probability distribution. Functions of parameters can sometimes be found, in overparameterized models, that are identifiable, and inferences based on these functions are legitimate. Examples are presented of overparameterized models that have been proposed in the context of several Bayesian methods for inferring the relative ages of nodes in a phylogeny when the substitution rate evolves over time.     

A Bayesian perspective on a non-parsimonious parsimony model

Several stochastic models of character change, when implemented in a maximum likelihood framework, are known to give a correspondence between the maximum parsimony method and the method of maximum likelihood. One such model has an independently estimated branch-length parameter for each site and each branch of the phylogenetic tree. This model--the no-common-mechanism model--has many parameters, and, in fact, the number of parameters increases as fast as the alignment is extended. We take a Bayesian approach to the no-common-mechanism model and place independent gamma prior probability distributions on the branch-length parameters. We are able to analytically integrate over the branch lengths, and this allowed us to implement an efficient Markov chain Monte Carlo method for exploring the space of phylogenetic trees. We were able to reliably estimate the posterior probabilities of clades for phylogenetic trees of up to 500 sequences. However, the Bayesian approach to the problem, at least as implemented here with an independent prior on the length of each branch, does not tame the behavior of the branch-length parameters. The integrated likelihood appears to be a simple rescaling of the parsimony score for a tree, and the marginal posterior probability distribution of the length of a branch is dependent upon how the maximum parsimony method reconstructs the characters at the interior nodes of the tree. The method we describe, however, is of potential importance in the analysis of morphological character data and also for improving the behavior of Markov chain Monte Carlo methods implemented for models in which sites share a common branch-length parameter.     

Predicting variability in biological control of a plant-pathogen system using stochastic models.

A stochastic model for the dynamics of a plant-pathogen interaction is developed and fitted to observations of the fungal pathogen Rhizoctonia solani (Kühn) in radish (Raphanus sativus L.), in both the presence and absence of the antagonistic fungus Trichoderma viride (Pers ex Gray). The model incorporates parameters for primary and secondary infection mechanisms and for characterizing the time-varying susceptibility of the host population. A parameter likelihood is developed and used to fit the model to data from microcosm experiments. It is shown that the stochastic model accounts well for observed variability both within and between treatments. Moreover, it enables us to describe the time evolution of the probability distribution for the variability among replicate epidemics in terms of the underlying epidemiological parameters for primary and secondary infection and decay in susceptibility. Consideration of profile likelihoods for each parameter provides strong evidence that T. viride mainly affects primary infection. By using the stochastic model to study the dependence of the probability distribution of disease levels on the primary infection rate we are therefore able to predict the effectiveness of a widely used biological control agent.     

Integral equation methods for computing likelihoods and their derivatives in the stochastic integrate-and-fire model

We recently introduced likelihood-based methods for fitting stochastic integrate-and-fire models to spike train data. The key component of this method involves the likelihood that the model will emit a spike at a given time t. Computing this likelihood is equivalent to computing a Markov first passage time density (the probability that the model voltage crosses threshold for the first time at time t). Here we detail an improved method for computing this likelihood, based on solving a certain integral equation. This integral equation method has several advantages over the techniques discussed in our previous work: in particular, the new method has fewer free parameters and is easily differentiable (for gradient computations). The new method is also easily adaptable for the case in which the model conductance, not just the input current, is time-varying. Finally, we describe how to incorporate large deviations approximations to very small likelihoods. Action Editor: Barry J. Richmond     

A unifying framework for chaos and stochastic stability in discrete population models

 In this paper we propose a general framework for discrete time one-dimensional Markov population models which is based on two fundamental premises in population dynamics. We show that this framework incorporates both earlier population models, like the Ricker and Hassell models, and experimental observations concerning the structure of density dependence. The two fundamental premises of population dynamics are sufficient to guarantee that the model will exhibit chaotic behaviour for high values of the natural growth and the density-dependent feedback, and this observation is independent of the particular structure of the model. We also study these models when the environment of the population varies stochastically and address the question under what conditions we can find an invariant probability distribution for the population under consideration. The sufficient conditions for this stochastic stability that we derive are of some interest, since studying certain statistical characteristics of these stochastic population processes may only be possible if the process converges to such an invariant distribution. Received 15 May 1995; received in revised form 17 April 1996     

Predicted steady-state cell size distributions for various growth models

The question of how an individual bacterial cell grows during its life cycle remains controversial. In 1962 Collins and Richmond derived a very general expression relating the size distributions of newborn, dividing and extant cells in steady-state growth and their growth rate; it represents the most powerful framework currently available for the analysis of bacterial growth kinetics. The Collins-Richmond equation is in effect a statement of the conservation of cell numbers for populations in steady-state exponential growth. It has usually been used to calculate the growth rate from a measured cell size distribution under various assumptions regarding the dividing and newborn cell distributions, but can also be applied in reverse--to compute the theoretical cell size distribution from a specified growth law. This has the advantage that it is not limited to models in which growth rate is a deterministic function of cell size, such as in simple exponential or linear growth, but permits evaluation of far more sophisticated hypotheses. Here we employed this reverse approach to obtain theoretical cell size distributions for two exponential and six linear growth models. The former differ as to whether there exists in each cell a minimal size that does not contribute to growth, the latter as to when the presumptive doubling of the growth rate takes place: in the linear age models, it is taken to occur at a particular cell age, at a fixed time prior to division, or at division itself; in the linear size models, the growth rate is considered to double with a constant probability from cell birth, with a constant probability but only after the cell has reached a minimal size, or after the minimal size has been attained but with a probability that increases linearly with cell size. Each model contains a small number of adjustable parameters but no assumptions other than that all cells obey the same growth law. In the present article, the various growth laws are described and rigorous mathematical expressions developed to predict the size distribution of extant cells in steady-state exponential growth; in the following paper, these predictions are tested against high-quality experimental data.     

Parameter inference for stochastic biochemical models from perturbation experiments parallelised at the single cell level

Understanding and characterising biochemical processes inside single cells requires experimental platforms that allow one to perturb and observe the dynamics of such processes as well as computational methods to build and parameterise models from the collected data. Recent progress with experimental platforms and optogenetics has made it possible to expose each cell in an experiment to an individualised input and automatically record cellular responses over days with fine time resolution. However, methods to infer parameters of stochastic kinetic models from single-cell longitudinal data have generally been developed under the assumption that experimental data is sparse and that responses of cells to at most a few different input perturbations can be observed. Here, we investigate and compare different approaches for calculating parameter likelihoods of single-cell longitudinal data based on approximations of the chemical master equation (CME) with a particular focus on coupling the linear noise approximation (LNA) or moment closure methods to a Kalman filter. We show that, as long as cells are measured sufficiently frequently, coupling the LNA to a Kalman filter allows one to accurately approximate likelihoods and to infer model parameters from data even in cases where the LNA provides poor approximations of the CME. Furthermore, the computational cost of filtering-based iterative likelihood evaluation scales advantageously in the number of measurement times and different input perturbations and is thus ideally suited for data obtained from modern experimental platforms. To demonstrate the practical usefulness of these results, we perform an experiment in which single cells, equipped with an optogenetic gene expression system, are exposed to various different light-input sequences and measured at several hundred time points and use parameter inference based on iterative likelihood evaluation to parameterise a stochastic model of the system.     

Bayesian Analysis of Growth Curves Using Mixed Models Defined by Stochastic Differential Equations

Summary Growth curve data consist of repeated measurements of a continuous growth process over time in a population of individuals. These data are classically analyzed by nonlinear mixed models. However, the standard growth functions used in this context prescribe monotone increasing growth and can fail to model unexpected changes in growth rates. We propose to model these variations using stochastic differential equations (SDEs) that are deduced from the standard deterministic growth function by adding random variations to the growth dynamics. A Bayesian inference of the parameters of these SDE mixed models is developed. In the case when the SDE has an explicit solution, we describe an easily implemented Gibbs algorithm. When the conditional distribution of the diffusion process has no explicit form, we propose to approximate it using the Euler–Maruyama scheme. Finally, we suggest validating the SDE approach via criteria based on the predictive posterior distribution. We illustrate the efficiency of our method using the Gompertz function to model data on chicken growth, the modeling being improved by the SDE approach.