Genome-wide single-molecule analysis of long-read DNA methylation reveals heterogeneous patterns at heterochromatin that reflect nucleosome organisation.

High-throughput sequencing technology is central to our current understanding of the human methylome. The vast majority of studies use chemical conversion to analyse bulk-level patterns of DNA methylation across the genome from a population of cells. While this technology has been used to probe single-molecule methylation patterns, such analyses are limited to short reads of a few hundred basepairs. DNA methylation can also be directly detected using Nanopore sequencing which can generate reads measuring megabases in length. However, thus far these analyses have largely focused on bulk-level assessment of DNA methylation. Here, we analyse DNA methylation in single Nanopore reads from human lymphoblastoid cells, to show that bulk-level metrics underestimate large-scale heterogeneity in the methylome. We use the correlation in methylation state between neighbouring sites to quantify single-molecule heterogeneity and find that heterogeneity varies significantly across the human genome, with some regions having heterogeneous methylation patterns at the single-molecule level and others possessing more homogeneous methylation patterns. By comparing the genomic distribution of the correlation to epigenomic annotations, we find that the greatest heterogeneity in single-molecule patterns is observed within heterochromatic partially methylated domains (PMDs). In contrast, reads originating from euchromatic regions and gene bodies have more ordered DNA methylation patterns. By analysing the patterns of single molecules in more detail, we show the existence of a nucleosome-scale periodicity in DNA methylation that accounts for some of the heterogeneity we uncover in long single-molecule DNA methylation patterns. We find that this periodic structure is partially masked in bulk data and correlates with DNA accessibility as measured by nanoNOMe-seq, suggesting that it could be generated by nucleosomes. Our findings demonstrate the power of single-molecule analysis of long-read data to understand the structure of the human methylome.

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.

Solving stochastic gene-expression models using queueing theory: A tutorial review.

Stochastic models of gene expression are typically formulated using the chemical master equation, which can be solved exactly or approximately using a repertoire of analytical methods. Here, we provide a tutorial review of an alternative approach based on queueing theory that has rarely been used in the literature of gene expression. We discuss the interpretation of six types of infinite-server queues from the angle of stochastic single-cell biology and provide analytical expressions for the stationary and nonstationary distributions and/or moments of mRNA/protein numbers and bounds on the Fano factor. This approach may enable the solution of complex models that have hitherto evaded analytical solution.

Solving the time-dependent protein distributions for autoregulated bursty gene expression using spectral decomposition.

In this study, we obtain an exact time-dependent solution of the chemical master equation (CME) of an extension of the two-state telegraph model describing bursty or non-bursty protein expression in the presence of positive or negative autoregulation. Using the method of spectral decomposition, we show that the eigenfunctions of the generating function solution of the CME are Heun functions, while the eigenvalues can be determined by solving a continued fraction equation. Our solution generalizes and corrects a previous time-dependent solution for the CME of a gene circuit describing non-bursty protein expression in the presence of negative autoregulation [Ramos et al., Phys. Rev. E 83, 062902 (2011)]. In particular, we clarify that the eigenvalues are generally not real as previously claimed. We also investigate the relationship between different types of dynamic behavior and the type of feedback, the protein burst size, and the gene switching rate.

DelaySSAToolkit.jl: stochastic simulation of reaction systems with time delays in Julia.

SUMMARY: DelaySSAToolkit.jl is a Julia package for modelling reaction systems with non-Markovian dynamics, specifically those with time delays. These delays implicitly capture multiple intermediate reaction steps and hence serve as an effective model reduction technique for complex systems in biology, chemistry, ecology and genetics. The package implements a variety of exact formulations of the delay stochastic simulation algorithm. AVAILABILITY AND IMPLEMENTATION: The source code and documentation of DelaySSAToolkit.jl are available at https://github.com/palmtree2013/DelaySSAToolkit.jl.

Characterizing non-exponential growth and bimodal cell size distributions in fission yeast: An analytical approach.

Unlike many single-celled organisms, the growth of fission yeast cells within a cell cycle is not exponential. It is rather characterized by three distinct phases (elongation, septation, and reshaping), each with a different growth rate. Experiments also showed that the distribution of cell size in a lineage can be bimodal, unlike the unimodal distributions measured for the bacterium Escherichia coli. Here we construct a detailed stochastic model of cell size dynamics in fission yeast. The theory leads to analytic expressions for the cell size and the birth size distributions, and explains the origin of bimodality seen in experiments. In particular, our theory shows that the left peak in the bimodal distribution is associated with cells in the elongation phase, while the right peak is due to cells in the septation and reshaping phases. We show that the size control strategy, the variability in the added size during a cell cycle, and the fraction of time spent in each of the three cell growth phases have a strong bearing on the shape of the cell size distribution. Furthermore, we infer all the parameters of our model by matching the theoretical cell size and birth size distributions to those from experimental single-cell time-course data for seven different growth conditions. Our method provides a much more accurate means of determining the size control strategy (timer, adder or sizer) than the standard method based on the slope of the best linear fit between the birth and division sizes. We also show that the variability in added size and the strength of size control in fission yeast depend weakly on the temperature but strongly on the culture medium. More importantly, we find that stronger size homeostasis and larger added size variability are required for fission yeast to adapt to unfavorable environmental conditions.

Concentration fluctuations in growing and dividing cells: Insights into the emergence of concentration homeostasis.

Intracellular reaction rates depend on concentrations and hence their levels are often regulated. However classical models of stochastic gene expression lack a cell size description and cannot be used to predict noise in concentrations. Here, we construct a model of gene product dynamics that includes a description of cell growth, cell division, size-dependent gene expression, gene dosage compensation, and size control mechanisms that can vary with the cell cycle phase. We obtain expressions for the approximate distributions and power spectra of concentration fluctuations which lead to insight into the emergence of concentration homeostasis. We find that (i) the conditions necessary to suppress cell division-induced concentration oscillations are difficult to achieve; (ii) mRNA concentration and number distributions can have different number of modes; (iii) two-layer size control strategies such as sizer-timer or adder-timer are ideal because they maintain constant mean concentrations whilst minimising concentration noise; (iv) accurate concentration homeostasis requires a fine tuning of dosage compensation, replication timing, and size-dependent gene expression; (v) deviations from perfect concentration homeostasis show up as deviations of the concentration distribution from a gamma distribution. Some of these predictions are confirmed using data for E. coli, fission yeast, and budding yeast.

Model reduction for the Chemical Master Equation: An information-theoretic approach.

The complexity of mathematical models in biology has rendered model reduction an essential tool in the quantitative biologist's toolkit. For stochastic reaction networks described using the Chemical Master Equation, commonly used methods include time-scale separation, Linear Mapping Approximation, and state-space lumping. Despite the success of these techniques, they appear to be rather disparate, and at present, no general-purpose approach to model reduction for stochastic reaction networks is known. In this paper, we show that most common model reduction approaches for the Chemical Master Equation can be seen as minimizing a well-known information-theoretic quantity between the full model and its reduction, the Kullback-Leibler divergence defined on the space of trajectories. This allows us to recast the task of model reduction as a variational problem that can be tackled using standard numerical optimization approaches. In addition, we derive general expressions for propensities of a reduced system that generalize those found using classical methods. We show that the Kullback-Leibler divergence is a useful metric to assess model discrepancy and to compare different model reduction techniques using three examples from the literature: an autoregulatory feedback loop, the Michaelis-Menten enzyme system, and a genetic oscillator.

Exact solution of a three-stage model of stochastic gene expression including cell-cycle dynamics.

The classical three-stage model of stochastic gene expression predicts the statistics of single cell mRNA and protein number fluctuations as a function of the rates of promoter switching, transcription, translation, degradation and dilution. While this model is easily simulated, its analytical solution remains an unsolved problem. Here we modify this model to explicitly include cell-cycle dynamics and then derive an exact solution for the time-dependent joint distribution of mRNA and protein numbers. We show large differences between this model and the classical model which captures cell-cycle effects implicitly via effective first-order dilution reactions. In particular we find that the Fano factor of protein numbers calculated from a population snapshot measurement are underestimated by the classical model whereas the correlation between mRNA and protein can be either over- or underestimated, depending on the timescales of mRNA degradation and promoter switching relative to the mean cell-cycle duration time.

Analytical distributions for detailed models of stochastic gene expression in eukaryotic cells.

The stochasticity of gene expression presents significant challenges to the modeling of genetic networks. A two-state model describing promoter switching, transcription, and messenger RNA (mRNA) decay is the standard model of stochastic mRNA dynamics in eukaryotic cells. Here, we extend this model to include mRNA maturation, cell division, gene replication, dosage compensation, and growth-dependent transcription. We derive expressions for the time-dependent distributions of nascent mRNA and mature mRNA numbers, provided two assumptions hold: 1) nascent mRNA dynamics are much faster than those of mature mRNA; and 2) gene-inactivation events occur far more frequently than gene-activation events. We confirm that thousands of eukaryotic genes satisfy these assumptions by using data from yeast, mouse, and human cells. We use the expressions to perform a sensitivity analysis of the coefficient of variation of mRNA fluctuations averaged over the cell cycle, for a large number of genes in mouse embryonic stem cells, identifying degradation and gene-activation rates as the most sensitive parameters. Furthermore, it is shown that, despite the model's complexity, the time-dependent distributions predicted by our model are generally well approximated by the negative binomial distribution. Finally, we extend our model to include translation, protein decay, and auto-regulatory feedback, and derive expressions for the approximate time-dependent protein-number distributions, assuming slow protein decay. Our expressions enable us to study how complex biological processes contribute to the fluctuations of gene products in eukaryotic cells, as well as allowing a detailed quantitative comparison with experimental data via maximum-likelihood methods.

MomentClosure.jl: automated moment closure approximations in Julia.

SUMMARY: MomentClosure.jl is a Julia package providing automated derivation of the time-evolution equations of the moments of molecule numbers for virtually any chemical reaction network using a wide range of moment closure approximations. It extends the capabilities of modelling stochastic biochemical systems in Julia and can be particularly useful when exact analytic solutions of the chemical master equation are unavailable and when Monte Carlo simulations are computationally expensive. AVAILABILITY AND IMPLEMENTATION: MomentClosure.jl is freely accessible under the MIT licence. Source code and documentation are available at https://github.com/augustinas1/MomentClosure.jl.

Stochastic Modeling of Autoregulatory Genetic Feedback Loops: A Review and Comparative Study.

Autoregulatory feedback loops are one of the most common network motifs. A wide variety of stochastic models have been constructed to understand how the fluctuations in protein numbers in these loops are influenced by the kinetic parameters of the main biochemical steps. These models differ according to 1) which subcellular processes are explicitly modeled, 2) the modeling methodology employed (discrete, continuous, or hybrid), and 3) whether they can be analytically solved for the steady-state distribution of protein numbers. We discuss the assumptions and properties of the main models in the literature, summarize our current understanding of the relationship between them, and highlight some of the insights gained through modeling.

A Stochastic Model of Gene Expression with Polymerase Recruitment and Pause Release.

Transcriptional bursting is a major source of noise in gene expression. The telegraph model of gene expression, whereby transcription switches between on and off states, is the dominant model for bursting. Recently, it was shown that the telegraph model cannot explain a number of experimental observations from perturbation data. Here, we study an alternative model that is consistent with the data and which explicitly describes RNA polymerase recruitment and polymerase pause release, two steps necessary for messenger RNA (mRNA) production. We derive the exact steady-state distribution of mRNA numbers and an approximate steady-state distribution of protein numbers, which are given by generalized hypergeometric functions. The theory is used to calculate the relative sensitivity of the coefficient of variation of mRNA fluctuations for thousands of genes in mouse fibroblasts. This indicates that the size of fluctuations is mostly sensitive to the rate of burst initiation and the mRNA degradation rate. Furthermore, we show that 1) the time-dependent distribution of mRNA numbers is accurately approximated by a modified telegraph model with a Michaelis-Menten like dependence of the effective transcription rate on RNA polymerase abundance, and 2) the model predicts that if the polymerase recruitment rate is comparable or less than the pause release rate, then upon gene replication, the mean number of RNA per cell remains approximately constant. This gene dosage compensation property has been experimentally observed and cannot be explained by the telegraph model with constant rates.

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.

Dynamical phase diagram of an auto-regulating gene in fast switching conditions.

While the steady-state behavior of stochastic gene expression with auto-regulation has been extensively studied, its time-dependent behavior has received much less attention. Here, under the assumption of fast promoter switching, we derive and solve a reduced chemical master equation for an auto-regulatory gene circuit with translational bursting and cooperative protein-gene interactions. The analytical expression for the time-dependent probability distribution of protein numbers enables a fast exploration of large swaths of the parameter space. For a unimodal initial distribution, we identify three distinct types of stochastic dynamics: (i) the protein distribution remains unimodal at all times; (ii) the protein distribution becomes bimodal at intermediate times and then reverts back to being unimodal at long times (transient bimodality); and (iii) the protein distribution switches to being bimodal at long times. For each of these, the deterministic model predicts either monostable or bistable behavior, and hence, there exist six dynamical phases in total. We investigate the relationship of the six phases to the transcription rates, the protein binding and unbinding rates, the mean protein burst size, the degree of cooperativity, the relaxation time to the steady state, the protein mean, and the type of feedback loop (positive or negative). We show that transient bimodality is a noise-induced phenomenon that occurs when the protein expression is sufficiently bursty, and we use a theory to estimate the observation time window when it is manifested.

Small protein number effects in stochastic models of autoregulated bursty gene expression.

A stochastic model of autoregulated bursty gene expression by Kumar et al. [Phys. Rev. Lett. 113, 268105 (2014)] has been exactly solved in steady-state conditions under the implicit assumption that protein numbers are sufficiently large such that fluctuations in protein numbers due to reversible protein-promoter binding can be ignored. Here, we derive an alternative model that takes into account these fluctuations and, hence, can be used to study low protein number effects. The exact steady-state protein number distribution is derived as a sum of Gaussian hypergeometric functions. We use the theory to study how promoter switching rates and the type of feedback influence the size of protein noise and noise-induced bistability. Furthermore, we show that our model predictions for the protein number distribution are significantly different from those of Kumar et al. when the protein mean is small, gene switching is fast, and protein binding to the gene is faster than the reverse unbinding reaction.

Stochastic time-dependent enzyme kinetics: Closed-form solution and transient bimodality.

We derive an approximate closed-form solution to the chemical master equation describing the Michaelis-Menten reaction mechanism of enzyme action. In particular, assuming that the probability of a complex dissociating into an enzyme and substrate is significantly larger than the probability of a product formation event, we obtain expressions for the time-dependent marginal probability distributions of the number of substrate and enzyme molecules. For delta function initial conditions, we show that the substrate distribution is either unimodal at all times or else becomes bimodal at intermediate times. This transient bimodality, which has no deterministic counterpart, manifests when the initial number of substrate molecules is much larger than the total number of enzyme molecules and if the frequency of enzyme-substrate binding events is large enough. Furthermore, we show that our closed-form solution is different from the solution of the chemical master equation reduced by means of the widely used discrete stochastic Michaelis-Menten approximation, where the propensity for substrate decay has a hyperbolic dependence on the number of substrate molecules. The differences arise because the latter does not take into account enzyme number fluctuations, while our approach includes them. We confirm by means of a stochastic simulation of all the elementary reaction steps in the Michaelis-Menten mechanism that our closed-form solution is accurate over a larger region of parameter space than that obtained using the discrete stochastic Michaelis-Menten approximation.

Revisiting the Reduction of Stochastic Models of Genetic Feedback Loops with Fast Promoter Switching.

Propensity functions of the Hill type are commonly used to model transcriptional regulation in stochastic models of gene expression. This leads to an effective reduced master equation for the mRNA and protein dynamics only. Based on deterministic considerations, it is often stated or tacitly assumed that such models are valid in the limit of rapid promoter switching. Here, starting from the chemical master equation describing promoter-protein interactions, mRNA transcription, protein translation, and decay, we prove that in the limit of fast promoter switching, the distribution of protein numbers is different than that given by standard stochastic models with Hill-type propensities. We show the differences are pronounced whenever the protein-DNA binding rate is much larger than the unbinding rate, a special case of fast promoter switching. Furthermore, we show using both theory and simulations that use of the standard stochastic models leads to drastically incorrect predictions for the switching properties of positive feedback loops and that these differences decrease with increasing mean protein burst size. Our results confirm that commonly used stochastic models of gene regulatory networks are only accurate in a subset of the parameter space consistent with rapid promoter switching.

Spatial Stochastic Intracellular Kinetics: A Review of Modelling Approaches.

Models of chemical kinetics that incorporate both stochasticity and diffusion are an increasingly common tool for studying biology. The variety of competing models is vast, but two stand out by virtue of their popularity: the reaction-diffusion master equation and Brownian dynamics. In this review, we critically address a number of open questions surrounding these models: How can they be justified physically? How do they relate to each other? How do they fit into the wider landscape of chemical models, ranging from the rate equations to molecular dynamics? This review assumes no prior knowledge of modelling chemical kinetics and should be accessible to a wide range of readers.

Insight into nuclear body formation of phytochromes through stochastic modelling and experiment.

Spatial relocalization of proteins is crucial for the correct functioning of living cells. An interesting example of spatial ordering is the light-induced clustering of plant photoreceptor proteins. Upon irradiation by white or red light, the red light-active phytochrome, phytochrome B, enters the nucleus and accumulates in large nuclear bodies (NBs). The underlying physical process of nuclear body formation remains unclear, but phytochrome B is thought to coagulate via a simple protein-protein binding process. We measure, for the first time, the distribution of the number of phytochrome B-containing NBs as well as their volume distribution. We show that the experimental data cannot be explained by a stochastic model of nuclear body formation via simple protein-protein binding processes using physically meaningful parameter values. Rather modelling suggests that the data is consistent with a two step process: a fast nucleation step leading to macroparticles followed by a subsequent slow step in which the macroparticles bind to form the nuclear body. An alternative explanation for the observed nuclear body distribution is that the phytochromes bind to a so far unknown molecular structure. We believe it is likely this result holds more generally for other nuclear body-forming plant photoreceptors and proteins.