Moment-based inference predicts bimodality in transient gene expression

Advanced methods for gene network identification and noise decomposition from single-cell data

Central to analyzing noisy gene expression systems is solving the Chemical Master Equation (CME), which characterizes the probability evolution of the reacting species' copy numbers. Solving CMEs for high-dimensional systems suffers from the curse of dimensionality. Here, we propose a computational method for improved scalability through a divide-and-conquer strategy that optimally decomposes the whole system into a leader system and several conditionally independent follower subsystems. The CME is solved by combining Monte Carlo estimation for the leader system with stochastic filtering procedures for the follower subsystems. We demonstrate this method with high-dimensional numerical examples and apply it to identify a yeast transcription system at the single-cell resolution, leveraging mRNA time-course experimental data. The identification results enable an accurate examination of the heterogeneity in rate parameters among isogenic cells. To validate this result, we develop a noise decomposition technique exploiting time-course data but requiring no supplementary components, e.g., dual-reporters.

Quantifying and correcting bias in transcriptional parameter inference from single-cell data

The snapshot distribution of mRNA counts per cell can be measured using single-molecule fluorescence in situ hybridization or single-cell RNA sequencing. These distributions are often fit to the steady-state distribution of the two-state telegraph model to estimate the three transcriptional parameters for a gene of interest: mRNA synthesis rate, the switching on rate (the on state being the active transcriptional state), and the switching off rate. This model assumes no extrinsic noise, i.e., parameters do not vary between cells, and thus estimated parameters are to be understood as approximating the average values in a population. The accuracy of this approximation is currently unclear. Here, we develop a theory that explains the size and sign of estimation bias when inferring parameters from single-cell data using the standard telegraph model. We find specific bias signatures depending on the source of extrinsic noise (which parameter is most variable across cells) and the mode of transcriptional activity. If gene expression is not bursty then the population averages of all three parameters are overestimated if extrinsic noise is in the synthesis rate; underestimation occurs if extrinsic noise is in the switching on rate; both underestimation and overestimation can occur if extrinsic noise is in the switching off rate. We find that some estimated parameters tend to infinity as the size of extrinsic noise approaches a critical threshold. In contrast when gene expression is bursty, we find that in all cases the mean burst size (ratio of the synthesis rate to the switching off rate) is overestimated while the mean burst frequency (the switching on rate) is underestimated. We estimate the size of extrinsic noise from the covariance matrix of sequencing data and use this together with our theory to correct published estimates of transcriptional parameters for mammalian genes.

Stochastic chemical kinetics of cell fate decision systems: From single cells to populations and back

Stochastic chemical kinetics is a widely used formalism for studying stochasticity of chemical reactions inside single cells. Experimental studies of reaction networks are generally performed with cells that are part of a growing population, yet the population context is rarely taken into account when models are developed. Models that neglect the population context lose their validity whenever the studied system influences traits of cells that can be selected in the population, a property that naturally arises in the complex interplay between single-cell and population dynamics of cell fate decision systems. Here, we represent such systems as absorbing continuous-time Markov chains. We show that conditioning on non-absorption allows one to derive a modified master equation that tracks the time evolution of the expected population composition within a growing population. This allows us to derive consistent population dynamics models from a specification of the single-cell process. We use this approach to classify cell fate decision systems into two types that lead to different characteristic phases in emerging population dynamics. Subsequently, we deploy the gained insights to experimentally study a recurrent problem in biology: how to link plasmid copy number fluctuations and plasmid loss events inside single cells to growth of cell populations in dynamically changing environments.

Exploring HERV-K (HML-2) Influence in Cancer and Prospects for Therapeutic Interventions

This review investigates the intricate role of human endogenous retroviruses (HERVs) in cancer development and progression, explicitly focusing on HERV-K (HML-2). This paper sheds light on the latest research advancements and potential treatment strategies by examining the historical context of HERVs and their involvement in critical biological processes such as embryonic development, immune response, and disease progression. This review covers computational modeling for drug-target binding assessment, systems biology modeling for simulating HERV-K viral cargo dynamics, and using antiviral drugs to combat HERV-induced diseases. The findings presented in this review contribute to our understanding of HERV-mediated disease mechanisms and provide insights into future therapeutic approaches. They emphasize why HERV-K holds significant promise as a biomarker and a target.

MaxCal can infer models from coupled stochastic trajectories of gene expression and cell division

Gene expression is inherently noisy due to small numbers of proteins and nucleic acids inside a cell. Likewise, cell division is stochastic, particularly when tracking at the level of a single cell. The two can be coupled when gene expression affects the rate of cell division. Single-cell time-lapse experiments can measure both fluctuations by simultaneously recording protein levels inside a cell and its stochastic division. These information-rich noisy trajectory data sets can be harnessed to learn about the underlying molecular and cellular details that are often not known a priori. A critical question is: How can we infer a model given data where fluctuations at two levels-gene expression and cell division-are intricately convoluted? We show the principle of maximum caliber (MaxCal)-integrated within a Bayesian framework-can be used to infer several cellular and molecular details (division rates, protein production, and degradation rates) from these coupled stochastic trajectories (CSTs). We demonstrate this proof of concept using synthetic data generated from a known model. An additional challenge in data analysis is that trajectories are often not in protein numbers, but in noisy fluorescence that depends on protein number in a probabilistic manner. We again show that MaxCal can infer important molecular and cellular rates even when data are in fluorescence, another example of CST with three confounding factors-gene expression noise, cell division noise, and fluorescence distortion-all coupled. Our approach will provide guidance to build models in synthetic biology experiments as well as general biological systems where examples of CSTs are abundant.

Inferring transcriptional bursting kinetics from single-cell snapshot data using a generalized telegraph model

Gene expression has inherent stochasticity resulting from transcription's burst manners. Single-cell snapshot data can be exploited to rigorously infer transcriptional burst kinetics, using mathematical models as blueprints. The classical telegraph model (CTM) has been widely used to explain transcriptional bursting with Markovian assumptions. However, growing evidence suggests that the gene-state dwell times are generally non-exponential, as gene-state switching is a multi-step process in organisms. Therefore, interpretable non-Markovian mathematical models and efficient statistical inference methods are urgently required in investigating transcriptional burst kinetics. We develop an interpretable and tractable model, the generalized telegraph model (GTM), to characterize transcriptional bursting that allows arbitrary dwell-time distributions, rather than exponential distributions, to be incorporated into the ON and OFF switching process. Based on the GTM, we propose an inference method for transcriptional bursting kinetics using an approximate Bayesian computation framework. This method demonstrates an efficient and scalable estimation of burst frequency and burst size on synthetic data. Further, the application of inference to genome-wide data from mouse embryonic fibroblasts reveals that GTM would estimate lower burst frequency and higher burst size than those estimated by CTM. In conclusion, the GTM and the corresponding inference method are effective tools to infer dynamic transcriptional bursting from static single-cell snapshot data.

Modeling Cellular Signaling Variability Based on Single-Cell Data: The TGFβ-SMAD Signaling Pathway

Nongenetic heterogeneity is key to cellular decisions, as even genetically identical cells respond in very different ways to the same external stimulus, e.g., during cell differentiation or therapeutic treatment of disease. Strong heterogeneity is typically already observed at the level of signaling pathways that are the first sensors of external inputs and transmit information to the nucleus where decisions are made. Since heterogeneity arises from random fluctuations of cellular components, mathematical models are required to fully describe the phenomenon and to understand the dynamics of heterogeneous cell populations. Here, we review the experimental and theoretical literature on cellular signaling heterogeneity, with special focus on the TGFβ/SMAD signaling pathway.

Approximating solutions of the Chemical Master equation using neural networks

The Chemical Master Equation (CME) provides an accurate description of stochastic biochemical reaction networks in well-mixed conditions, but it cannot be solved analytically for most systems of practical interest. Although Monte Carlo methods provide a principled means to probe system dynamics, the large number of simulations typically required can render the estimation of molecule number distributions and other quantities infeasible. In this article, we aim to leverage the representational power of neural networks to approximate the solutions of the CME and propose a framework for the Neural Estimation of Stochastic Simulations for Inference and Exploration (Nessie). Our approach is based on training neural networks to learn the distributions predicted by the CME from relatively few stochastic simulations. We show on biologically relevant examples that simple neural networks with one hidden layer can capture highly complex distributions across parameter space, thereby accelerating computationally intensive tasks such as parameter exploration and inference.

Inference and uncertainty quantification of stochastic gene expression via synthetic models

Estimating uncertainty in model predictions is a central task in quantitative biology. Biological models at the single-cell level are intrinsically stochastic and nonlinear, creating formidable challenges for their statistical estimation which inevitably has to rely on approximations that trade accuracy for tractability. Despite intensive interest, a sweet spot in this trade-off has not been found yet. We propose a flexible procedure for uncertainty quantification in a wide class of reaction networks describing stochastic gene expression including those with feedback. The method is based on creating a tractable coarse-graining of the model that is learned from simulations, a synthetic model, to approximate the likelihood function. We demonstrate that synthetic models can substantially outperform state-of-the-art approaches on a number of non-trivial systems and datasets, yielding an accurate and computationally viable solution to uncertainty quantification in stochastic models of gene expression.

Bayesian Inference of Stochastic Dynamic Models Using Early-Rejection Methods Based on Sequential Stochastic Simulations

Stochastic modelling is an important method to investigate the functions of noise in a wide range of biological systems. However, the parameter inference for stochastic models is still a challenging problem partially due to the large computing time required for stochastic simulations. To address this issue, we propose a novel early-rejection method by using sequential stochastic simulations. We first show that a large number of stochastic simulations are required to obtain reliable inference results. Instead of generating a large number of simulations for each parameter sample, we propose to generate these simulations in a number of stages. The simulation process will go to the next stage only if the accuracy of simulations at the current stage satisfies a given error criterion. We propose a formula to determine the error criterion and use a stochastic differential equation model to examine the effects of different criteria. Three biochemical network models are used to evaluate the efficiency and accuracy of the proposed method. Numerical results suggest the proposed early-rejection method achieves substantial improvement in the efficiency for the inference of stochastic models.

Systematic inference identifies a major source of heterogeneity in cell signaling dynamics: The rate-limiting step number

Identifying the sources of cell-to-cell variability in signaling dynamics is essential to understand drug response variability and develop effective therapeutics. However, it is challenging because not all signaling intermediate reactions can be experimentally measured simultaneously. This can be overcome by replacing them with a single random time delay, but the resulting process is non-Markovian, making it difficult to infer cell-to-cell heterogeneity in reaction rates and time delays. To address this, we developed an efficient and scalable moment-based Bayesian inference method (MBI) with a user-friendly computational package that infers cell-to-cell heterogeneity in the non-Markovian signaling process. We applied MBI to single-cell expression profiles from promoters responding to antibiotics and discovered a major source of cell-to-cell variability in antibiotic stress response: the number of rate-limiting steps in signaling cascades. This knowledge can help identify effective therapies that destroy all pathogenic or cancer cells, and the approach can be applied to precision medicine.

Using single-cell models to predict the functionality of synthetic circuits at the population scale

SignificanceAt the single-cell level, biochemical processes are inherently stochastic. For many natural systems, the resulting cell-to-cell variability is exploited by microbial populations. In synthetic biology, however, the interplay of cell-to-cell variability and population processes such as selection or growth often leads to circuits not functioning as predicted by simple models. Here we show how multiscale stochastic kinetic models that simultaneously track single-cell and population processes can be obtained based on an augmentation of the chemical master equation. These models enable us to quantitatively predict complex population dynamics of a yeast optogenetic differentiation system from a specification of the circuit's components and to demonstrate how cell-to-cell variability can be exploited to purposefully create unintuitive circuit functionality.

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 calibration of a stochastic, multiscale agent-based model for predicting in vitro tumor growth

Hybrid multiscale agent-based models (ABMs) are unique in their ability to simulate individual cell interactions and microenvironmental dynamics. Unfortunately, the high computational cost of modeling individual cells, the inherent stochasticity of cell dynamics, and numerous model parameters are fundamental limitations of applying such models to predict tumor dynamics. To overcome these challenges, we have developed a coarse-grained two-scale ABM (cgABM) with a reduced parameter space that allows for an accurate and efficient calibration using a set of time-resolved microscopy measurements of cancer cells grown with different initial conditions. The multiscale model consists of a reaction-diffusion type model capturing the spatio-temporal evolution of glucose and growth factors in the tumor microenvironment (at tissue scale), coupled with a lattice-free ABM to simulate individual cell dynamics (at cellular scale). The experimental data consists of BT474 human breast carcinoma cells initialized with different glucose concentrations and tumor cell confluences. The confluence of live and dead cells was measured every three hours over four days. Given this model, we perform a time-dependent global sensitivity analysis to identify the relative importance of the model parameters. The subsequent cgABM is calibrated within a Bayesian framework to the experimental data to estimate model parameters, which are then used to predict the temporal evolution of the living and dead cell populations. To this end, a moment-based Bayesian inference is proposed to account for the stochasticity of the cgABM while quantifying uncertainties due to limited temporal observational data. The cgABM reduces the computational time of ABM simulations by 93% to 97% while staying within a 3% difference in prediction compared to ABM. Additionally, the cgABM can reliably predict the temporal evolution of breast cancer cells observed by the microscopy data with an average error and standard deviation for live and dead cells being 7.61±2.01 and 5.78±1.13, respectively.

The calibration of agent-based models of tumor cell growth to experimental data remains a challenge in computational oncology. Besides the computational cost of modeling thousands of agents, the model’s intrinsic stochasticity demands numerous realizations of the simulations to accurately represent the statistical features of the model predictions. We developed a hybrid, multiscale, coarse-grain, agent-based model that captures the growth and decline of human breast carcinoma cells under different initial conditions. We determined the effects of coarse-graining the ABM on the multiscale model output and the number of repetitions necessary to capture the stochastic transitions present in the model. We identified the most influential parameters on the model prediction through a sensitivity analysis and selected which parameters can be fixed and which ones should be calibrated. Using Bayesian calibration, we show that the model can accurately represent the experimental data. The validation step indicates that our model can reliably predict the in vitro temporal data, depending on the choice of the training (calibration data) sets.

Cell-to-cell variability in JAK2/STAT5 pathway components and cytoplasmic volumes defines survival threshold in erythroid progenitor cells

Survival or apoptosis is a binary decision in individual cells. However, at the cell-population level, a graded increase in survival of colony-forming unit-erythroid (CFU-E) cells is observed upon stimulation with erythropoietin (Epo). To identify components of Janus kinase 2/signal transducer and activator of transcription 5 (JAK2/STAT5) signal transduction that contribute to the graded population response, we extended a cell-population-level model calibrated with experimental data to study the behavior in single cells. The single-cell model shows that the high cell-to-cell variability in nuclear phosphorylated STAT5 is caused by variability in the amount of Epo receptor (EpoR):JAK2 complexes and of SHP1, as well as the extent of nuclear import because of the large variance in the cytoplasmic volume of CFU-E cells. 24-118 pSTAT5 molecules in the nucleus for 120 min are sufficient to ensure cell survival. Thus, variability in membrane-associated processes is sufficient to convert a switch-like behavior at the single-cell level to a graded population-level response.

Tensor-train approximation of the chemical master equation and its application for parameter inference

In this work, we perform Bayesian inference tasks for the chemical master equation in the tensor-train format. The tensor-train approximation has been proven to be very efficient in representing high-dimensional data arising from the explicit representation of the chemical master equation solution. An additional advantage of representing the probability mass function in the tensor-train format is that parametric dependency can be easily incorporated by introducing a tensor product basis expansion in the parameter space. Time is treated as an additional dimension of the tensor and a linear system is derived to solve the chemical master equation in time. We exemplify the tensor-train method by performing inference tasks such as smoothing and parameter inference using the tensor-train framework. A very high compression ratio is observed for storing the probability mass function of the solution. Since all linear algebra operations are performed in the tensor-train format, a significant reduction in the computational time is observed as well.

Neural network aided approximation and parameter inference of non-Markovian models of gene expression

Non-Markovian models of stochastic biochemical kinetics often incorporate explicit time delays to effectively model large numbers of intermediate biochemical processes. Analysis and simulation of these models, as well as the inference of their parameters from data, are fraught with difficulties because the dynamics depends on the system's history. Here we use an artificial neural network to approximate the time-dependent distributions of non-Markovian models by the solutions of much simpler time-inhomogeneous Markovian models; the approximation does not increase the dimensionality of the model and simultaneously leads to inference of the kinetic parameters. The training of the neural network uses a relatively small set of noisy measurements generated by experimental data or stochastic simulations of the non-Markovian model. We show using a variety of models, where the delays stem from transcriptional processes and feedback control, that the Markovian models learnt by the neural network accurately reflect the stochastic dynamics across parameter space.

Critical Comparison of MaxCal and Other Stochastic Modeling Approaches in Analysis of Gene Networks

Learning the underlying details of a gene network with feedback is critical in designing new synthetic circuits. Yet, quantitative characterization of these circuits remains limited. This is due to the fact that experiments can only measure partial information from which the details of the circuit must be inferred. One potentially useful avenue is to harness hidden information from single-cell stochastic gene expression time trajectories measured for long periods of time—recorded at frequent intervals—over multiple cells. This raises the feasibility vs. accuracy dilemma while deciding between different models of mining these stochastic trajectories. We demonstrate that inference based on the Maximum Caliber (MaxCal) principle is the method of choice by critically evaluating its computational efficiency and accuracy against two other typical modeling approaches: (i) a detailed model (DM) with explicit consideration of multiple molecules including protein-promoter interaction, and (ii) a coarse-grain model (CGM) using Hill type functions to model feedback. MaxCal provides a reasonably accurate model while being significantly more computationally efficient than DM and CGM. Furthermore, MaxCal requires minimal assumptions since it is a top-down approach and allows systematic model improvement by including constraints of higher order, in contrast to traditional bottom-up approaches that require more parameters or ad hoc assumptions. Thus, based on efficiency, accuracy, and ability to build minimal models, we propose MaxCal as a superior alternative to traditional approaches (DM, CGM) when inferring underlying details of gene circuits with feedback from limited data.

Stochastic models for single-cell data: Current challenges and the way forward

Although the quantity and quality of single-cell data have progressed rapidly, making quantitative predictions with single-cell stochastic models remains challenging. The stochastic nature of cellular processes leads to at least three challenges in building models with single-cell data: (a) because variability in single-cell data can be attributed to multiple different sources, it is difficult to rule out conflicting mechanistic models that explain the same data equally well; (b) the distinction between interesting biological variability and experimental variability is sometimes ambiguous; (c) the nonstandard distributions of single-cell data can lead to violations of the assumption of symmetric errors in least-squares fitting. In this review, we first discuss recent studies that overcome some of the challenges or set up a promising direction and then introduce some powerful statistical approaches utilized in these studies. We conclude that applying and developing statistical approaches could lead to further progress in building stochastic models for single-cell data.

Promoters adopt distinct dynamic manifestations depending on transcription factor context

Cells respond to external signals and stresses by activating transcription factors (TF), which induce gene expression changes. Prior work suggests that signal-specific gene expression changes are partly achieved because different gene promoters exhibit distinct induction dynamics in response to the same TF input signal. Here, using high-throughput quantitative single-cell measurements and a novel statistical method, we systematically analyzed transcriptional responses to a large number of dynamic TF inputs. In particular, we quantified the scaling behavior among different transcriptional features extracted from the measured trajectories such as the gene activation delay or duration of promoter activity. Surprisingly, we found that even the same gene promoter can exhibit qualitatively distinct induction and scaling behaviors when exposed to different dynamic TF contexts. While it was previously known that promoters fall into distinct classes, here we show that the same promoter can switch between different classes depending on context. Thus, promoters can adopt context-dependent "manifestations". Our analysis suggests that the full complexity of signal processing by genetic circuits may be significantly underestimated when studied in only specific contexts.

Statistics of Nascent and Mature RNA Fluctuations in a Stochastic Model of Transcriptional Initiation, Elongation, Pausing, and Termination

Recent advances in fluorescence microscopy have made it possible to measure the fluctuations of nascent (actively transcribed) RNA. These closely reflect transcription kinetics, as opposed to conventional measurements of mature (cellular) RNA, whose kinetics is affected by additional processes downstream of transcription. Here, we formulate a stochastic model which describes promoter switching, initiation, elongation, premature detachment, pausing, and termination while being analytically tractable. We derive exact closed-form expressions for the mean and variance of nascent RNA fluctuations on gene segments, as well as of total nascent RNA on a gene. We also obtain exact expressions for the first two moments of mature RNA fluctuations and approximate distributions for total numbers of nascent and mature RNA. Our results, which are verified by stochastic simulation, uncover the explicit dependence of the statistics of both types of RNA on transcriptional parameters and potentially provide a means to estimate parameter values from experimental data.

Identifiability analysis for stochastic differential equation models in systems biology

Mathematical models are routinely calibrated to experimental data, with goals ranging from building predictive models to quantifying parameters that cannot be measured. Whether or not reliable parameter estimates are obtainable from the available data can easily be overlooked. Such issues of parameter identifiability have important ramifications for both the predictive power of a model, and the mechanistic insight that can be obtained. Identifiability analysis is well-established for deterministic, ordinary differential equation (ODE) models, but there are no commonly adopted methods for analysing identifiability in stochastic models. We provide an accessible introduction to identifiability analysis and demonstrate how existing ideas for analysis of ODE models can be applied to stochastic differential equation (SDE) models through four practical case studies. To assess structural identifiability, we study ODEs that describe the statistical moments of the stochastic process using open-source software tools. Using practically motivated synthetic data and Markov chain Monte Carlo methods, we assess parameter identifiability in the context of available data. Our analysis shows that SDE models can often extract more information about parameters than deterministic descriptions. All code used to perform the analysis is available on Github.

Optimal Design of Single-Cell Experiments within Temporally Fluctuating Environments

Modern biological experiments are becoming increasingly complex, and designing these experiments to yield the greatest possible quantitative insight is an open challenge. Increasingly, computational models of complex stochastic biological systems are being used to understand and predict biological behaviors or to infer biological parameters. Such quantitative analyses can also help to improve experiment designs for particular goals, such as to learn more about specific model mechanisms or to reduce prediction errors in certain situations. A classic approach to experiment design is to use the Fisher information matrix (FIM), which quantifies the expected information a particular experiment will reveal about model parameters. The Finite State Projection based FIM (FSP-FIM) was recently developed to compute the FIM for discrete stochastic gene regulatory systems, whose complex response distributions do not satisfy standard assumptions of Gaussian variations. In this work, we develop the FSP-FIM analysis for a stochastic model of stress response genes in S. cerevisae under time-varying MAPK induction. We verify this FSP-FIM analysis and use it to optimize the number of cells that should be quantified at particular times to learn as much as possible about the model parameters. We then extend the FSP-FIM approach to explore how different measurement times or genetic modifications help to minimize uncertainty in the sensing of extracellular environments, and we experimentally validate the FSP-FIM to rank single-cell experiments for their abilities to minimize estimation uncertainty of NaCl concentrations during yeast osmotic shock. This work demonstrates the potential of quantitative models to not only make sense of modern biological data sets, but to close the loop between quantitative modeling and experimental data collection.

Inheritance and variability of kinetic gene expression parameters in microbial cells: modeling and inference from lineage tree data

Motivation

Modern experimental technologies enable monitoring of gene expression dynamics in individual cells and quantification of its variability in isogenic microbial populations. Among the sources of this variability is the randomness that affects inheritance of gene expression factors at cell division. Known parental relationships among individually observed cells provide invaluable information for the characterization of this extrinsic source of gene expression noise. Despite this fact, most existing methods to infer stochastic gene expression models from single-cell data dedicate little attention to the reconstruction of mother–daughter inheritance dynamics.

Results

Starting from a transcription and translation model of gene expression, we propose a stochastic model for the evolution of gene expression dynamics in a population of dividing cells. Based on this model, we develop a method for the direct quantification of inheritance and variability of kinetic gene expression parameters from single-cell gene expression and lineage data. We demonstrate that our approach provides unbiased estimates of mother–daughter inheritance parameters, whereas indirect approaches using lineage information only in the post-processing of individual-cell parameters underestimate inheritance. Finally, we show on yeast osmotic shock response data that daughter cell parameters are largely determined by the mother, thus confirming the relevance of our method for the correct assessment of the onset of gene expression variability and the study of the transmission of regulatory factors.

Availability and implementation

Software code is available at https://github.com/almarguet/IdentificationWithARME. Lineage tree data is available upon request.

Supplementary information

Supplementary material is available at Bioinformatics online.

Time-dependent product-form Poisson distributions for reaction networks with higher order complexes

It is well known that stochastically modeled reaction networks that are complex balanced admit a stationary distribution that is a product of Poisson distributions. In this paper, we consider the following related question: supposing that the initial distribution of a stochastically modeled reaction network is a product of Poissons, under what conditions will the distribution remain a product of Poissons for all time? By drawing inspiration from Crispin Gardiner's "Poisson representation" for the solution to the chemical master equation, we provide a necessary and sufficient condition for such a product-form distribution to hold for all time. Interestingly, the condition is a dynamical "complex-balancing" for only those complexes that have multiplicity greater than or equal to two (i.e. the higher order complexes that yield non-linear terms to the dynamics). We term this new condition the "dynamical and restricted complex balance" condition (DR for short).

Maximum Entropy Framework for Predictive Inference of Cell Population Heterogeneity and Responses in Signaling Networks

Predictive models of signaling networks are essential for understanding cell population heterogeneity and designing rational interventions in disease. However, using computational models to predict heterogeneity of signaling dynamics is often challenging because of the extensive variability of biochemical parameters across cell populations. Here, we describe a maximum entropy-based framework for inference of heterogeneity in dynamics of signaling networks (MERIDIAN). MERIDIAN estimates the joint probability distribution over signaling network parameters that is consistent with experimentally measured cell-to-cell variability of biochemical species. We apply the developed approach to investigate the response heterogeneity in the EGFR/Akt signaling network. Our analysis demonstrates that a significant fraction of cells exhibits high phosphorylated Akt (pAkt) levels hours after EGF stimulation. Our findings also suggest that cells with high EGFR levels predominantly contribute to the subpopulation of cells with high pAkt activity. We also discuss how MERIDIAN can be extended to accommodate various experimental measurements.

Filtering and inference for stochastic oscillators with distributed delays

Motivation

The time evolution of molecular species involved in biochemical reaction networks often arises from complex stochastic processes involving many species and reaction events. Inference for such systems is profoundly challenged by the relative sparseness of experimental data, as measurements are often limited to a small subset of the participating species measured at discrete time points. The need for model reduction can be realistically achieved for oscillatory dynamics resulting from negative translational and transcriptional feedback loops by the introduction of probabilistic time-delays. Although this approach yields a simplified model, inference is challenging and subject to ongoing research. The linear noise approximation (LNA) has recently been proposed to address such systems in stochastic form and will be exploited here.

Results

We develop a novel filtering approach for the LNA in stochastic systems with distributed delays, which allows the parameter values and unobserved states of a stochastic negative feedback model to be inferred from univariate time-series data. The performance of the methods is tested for simulated data. Results are obtained for real data when the model is fitted to imaging data on Cry1, a key gene involved in the mammalian central circadian clock, observed via a luciferase reporter construct in a mouse suprachiasmatic nucleus.

Availability and implementation

Programmes are written in MATLAB and Statistics Toolbox Release 2016 b, The MathWorks, Inc., Natick, Massachusetts, USA. Sample code and Cry1 data are available on GitHub https://github.com/scalderazzo/FLNADD.

Supplementary information

Supplementary data are available at Bioinformatics online.

Phase separation provides a mechanism to reduce noise in cells

Expression of proteins inside cells is noisy, causing variability in protein concentration among identical cells. A central problem in cellular control is how cells cope with this inherent noise. Compartmentalization of proteins through phase separation has been suggested as a potential mechanism to reduce noise, but systematic studies to support this idea have been missing. In this study, we used a physical model that links noise in protein concentration to theory of phase separation to show that liquid droplets can effectively reduce noise. We provide experimental support for noise reduction by phase separation using engineered proteins that form liquid-like compartments in mammalian cells. Thus, phase separation can play an important role in biological signal processing and control.

BAYESIAN INFERENCE OF STOCHASTIC REACTION NETWORKS USING MULTIFIDELITY SEQUENTIAL TEMPERED MARKOV CHAIN MONTE CARLO

Stochastic reaction network models are often used to explain and predict the dynamics of gene regulation in single cells. These models usually involve several parameters, such as the kinetic rates of chemical reactions, that are not directly measurable and must be inferred from experimental data. Bayesian inference provides a rigorous probabilistic framework for identifying these parameters by finding a posterior parameter distribution that captures their uncertainty. Traditional computational methods for solving inference problems such as Markov Chain Monte Carlo methods based on classical Metropolis-Hastings algorithm involve numerous serial evaluations of the likelihood function, which in turn requires expensive forward solutions of the chemical master equation (CME). We propose an alternate approach based on a multifidelity extension of the Sequential Tempered Markov Chain Monte Carlo (ST-MCMC) sampler. This algorithm is built upon Sequential Monte Carlo and solves the Bayesian inference problem by decomposing it into a sequence of efficiently solved subproblems that gradually increase both model fidelity and the influence of the observed data. We reformulate the finite state projection (FSP) algorithm, a well-known method for solving the CME, to produce a hierarchy of surrogate master equations to be used in this multifidelity scheme. To determine the appropriate fidelity, we introduce a novel information-theoretic criteria that seeks to extract the most information about the ultimate Bayesian posterior from each model in the hierarchy without inducing significant bias. This novel sampling scheme is tested with high performance computing resources using biologically relevant problems.

Dynamic distribution decomposition for single-cell snapshot time series identifies subpopulations and trajectories during iPSC reprogramming

Recent high-dimensional single-cell technologies such as mass cytometry are enabling time series experiments to monitor the temporal evolution of cell state distributions and to identify dynamically important cell states, such as fate decision states in differentiation. However, these technologies are destructive, and require analysis approaches that temporally map between cell state distributions across time points. Current approaches to approximate the single-cell time series as a dynamical system suffer from too restrictive assumptions about the type of kinetics, or link together pairs of sequential measurements in a discontinuous fashion. We propose Dynamic Distribution Decomposition (DDD), an operator approximation approach to infer a continuous distribution map between time points. On the basis of single-cell snapshot time series data, DDD approximates the continuous time Perron-Frobenius operator by means of a finite set of basis functions. This procedure can be interpreted as a continuous time Markov chain over a continuum of states. By only assuming a memoryless Markov (autonomous) process, the types of dynamics represented are more general than those represented by other common models, e.g., chemical reaction networks, stochastic differential equations. Furthermore, we can a posteriori check whether the autonomy assumptions are valid by calculation of prediction error-which we show gives a measure of autonomy within the studied system. The continuity and autonomy assumptions ensure that the same dynamical system maps between all time points, not arbitrarily changing at each time point. We demonstrate the ability of DDD to reconstruct dynamically important cell states and their transitions both on synthetic data, as well as on mass cytometry time series of iPSC reprogramming of a fibroblast system. We use DDD to find previously identified subpopulations of cells and to visualise differentiation trajectories. Dynamic Distribution Decomposition allows interpretation of high-dimensional snapshot time series data as a low-dimensional Markov process, thereby enabling an interpretable dynamics analysis for a variety of biological processes by means of identifying their dynamically important cell states.

Accuracy of parameter estimation for auto-regulatory transcriptional feedback loops from noisy data

Bayesian and non-Bayesian moment-based inference methods are commonly used to estimate the parameters defining stochastic models of gene regulatory networks from noisy single cell or population snapshot data. However, a systematic investigation of the accuracy of the predictions of these methods remains missing. Here, we present the results of such a study using synthetic noisy data of a negative auto-regulatory transcriptional feedback loop, one of the most common building blocks of complex gene regulatory networks. We study the error in parameter estimation as a function of (i) number of cells in each sample; (ii) the number of time points; (iii) the highest-order moment of protein fluctuations used for inference; (iv) the moment-closure method used for likelihood approximation. We find that for sample sizes typical of flow cytometry experiments, parameter estimation by maximizing the likelihood is as accurate as using Bayesian methods but with a much reduced computational time. We also show that the choice of moment-closure method is the crucial factor determining the maximum achievable accuracy of moment-based inference methods. Common likelihood approximation methods based on the linear noise approximation or the zero cumulants closure perform poorly for feedback loops with large protein-DNA binding rates or large protein bursts; this is exacerbated for highly heterogeneous cell populations. By contrast, approximating the likelihood using the linear-mapping approximation or conditional derivative matching leads to highly accurate parameter estimates for a wide range of conditions.

Quantifying heterogeneity of stochastic gene expression

The heterogeneity of stochastic gene expression, which refers to the temporal fluctuation in a gene product and its cell-to-cell variation, has attracted considerable interest from biologists, physicists, and mathematicians. The dynamics of protein production and degradation have been modeled as random processes with transition probabilities. However, there is a gap between theory and phenomena, particularly in terms of analytical formulation and parameter estimation. In this study, we propose a theoretical framework in which we present a basic model of a gene regulatory system, derive a steady-state solution, and provide a Bayesian approach for estimating the model parameters from single-cell experimental data. The proposed framework is demonstrated to be applicable for various scales of single-cell experiments at both the mRNA and protein levels and is useful for comparing kinetic parameters across species, genomes, and cell strains.

Bayesian Estimation for Stochastic Gene Expression Using Multifidelity Models

The finite state projection (FSP) approach to solving the chemical master equation has enabled successful inference of discrete stochastic models to predict single-cell gene regulation dynamics. Unfortunately, the FSP approach is highly computationally intensive for all but the simplest models, an issue that is highly problematic when parameter inference and uncertainty quantification takes enormous numbers of parameter evaluations. To address this issue, we propose two new computational methods for the Bayesian inference of stochastic gene expression parameters given single-cell experiments. We formulate and verify an adaptive delayed acceptance Metropolis-Hastings (ADAMH) algorithm to utilize with reduced Krylov-basis projections of the FSP. We then introduce an extension of the ADAMH into a hybrid scheme that consists of an initial phase to construct a reduced model and a faster second phase to sample from the approximate posterior distribution determined by the constructed model. We test and compare both algorithms to an adaptive Metropolis algorithm with full FSP-based likelihood evaluations on three example models and simulated data to show that the new ADAMH variants achieve substantial speedup in comparison to the full FSP approach. By reducing the computational costs of parameter estimation, we expect the ADAMH approach to enable efficient data-driven estimation for more complex gene regulation models.

Control of Hox transcription factor concentration and cell-to-cell variability by an auto-regulatory switch

The variability in transcription factor concentration among cells is an important developmental determinant, yet how variability is controlled remains poorly understood. Studies of variability have focused predominantly on monitoring mRNA production noise. Little information exists about transcription factor protein variability, as this requires the use of quantitative methods with single-molecule sensitivity. Using Fluorescence Correlation Spectroscopy (FCS), we have characterized the concentration and variability of 14 endogenously tagged TFs in live Drosophila imaginal discs. For the Hox TF Antennapedia, we investigated whether protein variability results from random stochastic events or is developmentally regulated. We found that Antennapedia transitioned from low concentration/high variability early, to high concentration/low variability later, in development. FCS and temporally resolved genetic studies uncovered that Antennapedia itself is necessary and sufficient to drive a developmental regulatory switch from auto-activation to auto-repression, thereby reducing variability. This switch is controlled by progressive changes in relative concentrations of preferentially activating and repressing Antennapedia isoforms, which bind chromatin with different affinities. Mathematical modeling demonstrated that the experimentally supported auto-regulatory circuit can explain the increase of Antennapedia concentration and suppression of variability over time.

Summary: Preferentially repressing and activating isoforms of the Hox transcription factor Antennapedia elicit a developmental regulatory switch from auto-activation to auto-repression that increases concentration and suppresses cell-to-cell variability over time.

Intrinsic and extrinsic noise of gene expression in lineage trees

Cell-to-cell heterogeneity is driven by stochasticity in intracellular reactions and the population dynamics. While these sources are usually studied separately, we develop an agent-based framework that accounts for both factors while tracking every single cell of a growing population. Apart from the common intrinsic variability, the framework also predicts extrinsic noise without the need to introduce fluctuating rate constants. Instead, extrinsic fluctuations are explained by cell cycle fluctuations and differences in cell age. We provide explicit formulas to quantify mean molecule numbers, intrinsic and extrinsic noise statistics in two-colour experiments. We find that these statistics differ significantly depending on the experimental setup used to observe the cells. We illustrate this fact using (i) averages over an isolated cell lineage tracked over many generations as observed in the mother machine, (ii) population snapshots with known cell ages as recorded in time-lapse microscopy, and (iii) snapshots with unknown cell ages as measured from static images or flow cytometry. Applying the method to models of stochastic gene expression and feedback regulation elucidates that isolated lineages, as compared to snapshot data, can significantly overestimate the mean number of molecules, overestimate extrinsic noise but underestimate intrinsic noise and have qualitatively different sensitivities to cell cycle fluctuations.

Extreme Osmotolerance and Halotolerance in Food-Relevant Yeasts and the Role of Glycerol-Dependent Cell Individuality

Osmotolerance or halotolerance are used to describe resistance to sugars and salt, or only salt, respectively. Here, a comprehensive screen of more than 600 different yeast isolates revealed that osmosensitive species were equally affected by NaCl and glucose. However, the relative toxicity of salt became increasingly prominent in more osmoresistant species. We confirmed that growth inhibition by glucose in a laboratory strain of Saccharomyces cerevisiae occurred at a lower water activity (A_w) than by salt (NaCl), and pre-growth in high levels of glucose or salt gave enhanced cross-resistance to either. Salt toxicity was largely due to osmotic stress but with an additive enhancement due to effects of the relevant cation. Almost all of the yeast isolates from the screen were also noted to exhibit hetero-resistance to both salt and sugar, whereby high concentrations restricted growth to a small minority of cells within the clonal populations. Rare resistant colonies required growth for up to 28 days to become visible. This cell individuality was more marked with salt than sugar, a possible further reflection of the ion toxicity effect. In both cases, heteroresistance in S. cerevisiae was strikingly dependent on the GPD1 gene product, important for glycerol synthesis. In contrast, a tps1Δ deletant impaired for trehalose showed altered MIC but no change in heteroresistance. Effects on heteroresistance were evident in chronic (but not acute) salt or glucose stress, particularly relevant to growth on low A_w foods. The study reports diverse osmotolerance and halotolerance phenotypes and heteroresistance across an extensive panel of yeast isolates, and indicates that Gpd1-dependent glycerol synthesis is a key determinant enabling growth of rare yeast subpopulations at low A_w, brought about by glucose and in particular salt.

The rise of the distributions: why non-normality is important for understanding the transcriptome and beyond

The application of statistics has been instrumental in clarifying our understanding of the genome. While insights have been derived for almost all levels of genome function, most importantly, statistics has had the greatest impact on improving our knowledge of transcriptional regulation. But the drive to extract the most meaningful inferences from big data can often force us to overlook the fundamental role that statistics plays, and specifically, the basic assumptions that we make about big data. Normality is a statistical property that is often swept up into an assumption that we may or may not be consciously aware of making. This review highlights the inherent value of non-normal distributions to big data analysis by discussing use cases of non-normality that focus on gene expression data. Collectively, these examples help to motivate the premise of why at this stage, now more than ever, non-normality is important for learning about gene regulation, transcriptomics, and more.

The finite state projection based Fisher information matrix approach to estimate information and optimize single-cell experiments

Modern optical imaging experiments not only measure single-cell and single-molecule dynamics with high precision, but they can also perturb the cellular environment in myriad controlled and novel settings. Techniques, such as single-molecule fluorescence in-situ hybridization, microfluidics, and optogenetics, have opened the door to a large number of potential experiments, which begs the question of how to choose the best possible experiment. The Fisher information matrix (FIM) estimates how well potential experiments will constrain model parameters and can be used to design optimal experiments. Here, we introduce the finite state projection (FSP) based FIM, which uses the formalism of the chemical master equation to derive and compute the FIM. The FSP-FIM makes no assumptions about the distribution shapes of single-cell data, and it does not require precise measurements of higher order moments of such distributions. We validate the FSP-FIM against well-known Fisher information results for the simple case of constitutive gene expression. We then use numerical simulations to demonstrate the use of the FSP-FIM to optimize the timing of single-cell experiments with more complex, non-Gaussian fluctuations. We validate optimal simulated experiments determined using the FSP-FIM with Monte-Carlo approaches and contrast these to experiment designs chosen by traditional analyses that assume Gaussian fluctuations or use the central limit theorem. By systematically designing experiments to use all of the measurable fluctuations, our method enables a key step to improve co-design of experiments and quantitative models.

Estimation methods for heterogeneous cell population models in systems biology

Heterogeneity among individual cells is a characteristic and relevant feature of living systems. A range of experimental techniques to investigate this heterogeneity is available, and multiple modelling frameworks have been developed to describe and simulate the dynamics of heterogeneous populations. Measurement data are used to adjust computational models, which results in parameter and state estimation problems. Methods to solve these estimation problems need to take the specific properties of data and models into account. The aim of this review is to give an overview on the state of the art in estimation methods for heterogeneous cell population data and models. The focus is on models based on the population balance equation, but stochastic and individual-based models are also discussed. It starts with a brief discussion of common experimental approaches and types of measurement data that can be obtained in this context. The second part describes computational modelling frameworks for heterogeneous populations and the types of estimation problems occurring for these models. The third part starts with a discussion of observability and identifiability properties, after which the computational methods to solve the various estimation problems are described.

Solving Stochastic Reaction Networks with Maximum Entropy Lagrange Multipliers

The time evolution of stochastic reaction networks can be modeled with the chemical master equation of the probability distribution. Alternatively, the numerical problem can be reformulated in terms of probability moment equations. Herein we present a new alternative method for numerically solving the time evolution of stochastic reaction networks. Based on the assumption that the entropy of the reaction network is maximum, Lagrange multipliers are introduced. The proposed method derives equations that model the time derivatives of these Lagrange multipliers. We present detailed steps to transform moment equations to Lagrange multiplier equations. In order to demonstrate the method, we present examples of non-linear stochastic reaction networks of varying degrees of complexity, including multistable and oscillatory systems. We find that the new approach is as accurate and significantly more efficient than Gillespie’s original exact algorithm for systems with small number of interacting species. This work is a step towards solving stochastic reaction networks accurately and efficiently.

Identifiability and Reconstruction of Biochemical Reaction Networks from Population Snapshot Data

Inference of biochemical network models from experimental data is a crucial problem in systems and synthetic biology that includes parameter calibration but also identification of unknown interactions. Stochastic modelling from single-cell data is known to improve identifiability of reaction network parameters for specific systems. However, general results are lacking, and the advantage over deterministic, population-average approaches has not been explored for network reconstruction. In this work, we study identifiability and propose new reconstruction methods for biochemical interaction networks. Focusing on population-snapshot data and networks with reaction rates affine in the state, for parameter estimation, we derive general methods to test structural identifiability and demonstrate them in connection with practical identifiability for a reporter gene in silico case study. In the same framework, we next develop a two-step approach to the reconstruction of unknown networks of interactions. We apply it to compare the achievable network reconstruction performance in a deterministic and a stochastic setting, showing the advantage of the latter, and demonstrate it on population-snapshot data from a simulated example.

Moment-Based Parameter Estimation for Stochastic Reaction Networks in Equilibrium

Calibrating parameters is a crucial problem within quantitative modeling approaches to reaction networks. Existing methods for stochastic models rely either on statistical sampling or can only be applied to small systems. Here, we present an inference procedure for stochastic models in equilibrium that is based on a moment matching scheme with optimal weighting and that can be used with high-throughput data like the one collected by flow cytometry. Our method does not require an approximation of the underlying equilibrium probability distribution and, if reaction rate constants have to be learned, the optimal values can be computed by solving a linear system of equations. We discuss important practical issues such as the selection of the moments and evaluate the effectiveness of the proposed approach on three case studies.

Selected-node stochastic simulation algorithm

Stochastic simulations of biochemical networks are of vital importance for understanding complex dynamics in cells and tissues. However, existing methods to perform such simulations are associated with computational difficulties and addressing those remains a daunting challenge to the present. Here we introduce the selected-node stochastic simulation algorithm (snSSA), which allows us to exclusively simulate an arbitrary, selected subset of molecular species of a possibly large and complex reaction network. The algorithm is based on an analytical elimination of chemical species, thereby avoiding explicit simulation of the associated chemical events. These species are instead described continuously in terms of statistical moments derived from a stochastic filtering equation, resulting in a substantial speedup when compared to Gillespie's stochastic simulation algorithm (SSA). Moreover, we show that statistics obtained via snSSA profit from a variance reduction, which can significantly lower the number of Monte Carlo samples needed to achieve a certain performance. We demonstrate the algorithm using several biological case studies for which the simulation time could be reduced by orders of magnitude.

Linear approximations of global behaviors in nonlinear systems with moderate or strong noise

While many physical or chemical systems can be modeled by nonlinear Langevin equations (LEs), dynamical analysis of these systems is challenging in the cases of moderate and strong noise. Here we develop a linear approximation scheme, which can transform an often intractable LE into a linear set of binomial moment equations (BMEs). This scheme provides a feasible way to capture nonlinear behaviors in the sense of probability distribution and is effective even when the noise is moderate or big. Based on BMEs, we further develop a noise reduction technique, which can effectively handle tough cases where traditional small-noise theories are inapplicable. The overall method not only provides an approximation-based paradigm to analysis of the local and global behaviors of nonlinear noisy systems but also has a wide range of applications.