Pathway dynamics can delineate the sources of transcriptional noise in gene expression

  1. Lucy Ham  Is a corresponding author
  2. Marcel Jackson
  3. Michael PH Stumpf  Is a corresponding author
  1. School of BioSciences, University of Melbourne, Australia
  2. Department of Mathematics and Statistics, La Trobe University, Australia
  3. School of Mathematics and Statistics, University of Melbourne, Australia

Abstract

Single-cell expression profiling opens up new vistas on cellular processes. Extensive cell-to-cell variability at the transcriptomic and proteomic level has been one of the stand-out observations. Because most experimental analyses are destructive we only have access to snapshot data of cellular states. This loss of temporal information presents significant challenges for inferring dynamics, as well as causes of cell-to-cell variability. In particular, we typically cannot separate dynamic variability from within cells (‘intrinsic noise’) from variability across the population (‘extrinsic noise’). Here, we make this non-identifiability mathematically precise, allowing us to identify new experimental set-ups that can assist in resolving this non-identifiability. We show that multiple generic reporters from the same biochemical pathways (e.g. mRNA and protein) can infer magnitudes of intrinsic and extrinsic transcriptional noise, identifying sources of heterogeneity. Stochastic simulations support our theory, and demonstrate that ‘pathway-reporters’ compare favourably to the well-known, but often difficult to implement, dual-reporter method.

eLife digest

In biology, seemingly random variation within or between cells can have significant effects on a number of cellular processes, like how cells divide and develop. For example, how often a gene is switched on, or ‘expressed’, can randomly fluctuate over time. This ‘noise’ may lead to a cell having slightly more of a particular molecule, causing it to behave differently to other cells in the population. However, it is currently unclear how this random variation is created and controlled in cells, and what effect this has on biological systems as a whole.

When a gene is expressed, its sequence typically gets copied in to a molecule called mRNA, which is then processed and used to build the protein encoded by the gene. By measuring the levels of mRNA molecules in individual cells, researchers have been able to investigate how gene expression varies within populations. These experiments are carried out on dead cells at a single point in time, and mathematical models are then applied to detect noise in the molecular data.

This approach, however, precludes how noise changes over time, making it difficult to determine the source of cell-to-cell variability. In particular, whether the variation detected is the result of genuine random molecular changes (intrinsic noise), or external factors – such as temperature and pH – fluctuating in the cells environment (extrinsic noise).

Here, Ham et al. have built on previous mathematical models to identify a new approach for investigating the source of molecular noise. They found that for any given gene it is impossible to understand what causes its activity levels to vary just from data on its mRNA levels. Instead, information on other molecules that are affected by expression of the gene (termed ‘pathway reporters’) can provide a clearer picture of whether molecular variability is the result of intrinsic or extrinsic noise.

The mathematical models developed by Ham et al. reveal what can and cannot be learned about noise from gene expression data. Furthermore, pathway-reporters are easier to measure experimentally than other reporters that are typically used to study the origins and effects of cell-to-cell variability. These findings could help researchers design single cell experiments that are better for studying noise, leading to a deeper understanding of how different types of variation impact cell biology.

Introduction

Noise is a fundamental aspect of every cellular process Shahrezaei and Swain, 2008b. Frequently it is even of functional importance, for example in driving cell-fate transitions. Sometimes it can afford evolutionary advantages, for example, in the context of bet-hedging strategies. Sometimes, it can be a nuisance, for example, when it makes cellular signal processing more difficult. But noise is nearly ubiquitous at the molecular scale, and its presence has profoundly shaped cellular life. Analysing and understanding the sources of noise, how it is propagated, amplified or attenuated, and how it can be controlled, has therefore become a cornerstone of modern molecular cell biology.

Noise arising in gene expression has arguably attracted most of the attention so far (but see e.g. Filippi et al., 2016 and Jetka et al., 2018 for the analysis of noise at the signalling level). Generally speaking, gene expression noise is separable into two sources of variability, as pioneered by Swain et al., 2002. Intrinsic noise is generated by the dynamics of the gene expression process itself. The process, however, is often influenced by other external factors, such as the availability of promoters and of RNA polymerase, the influence of long noncoding RNA as a transcriptional regulator Goodrich and Kugel, 2006, as well as differences in the cellular environment. Such sources of variability contribute extrinsic noise, and reflect the variation in gene expression and transcription activity across the cell population. As such, understanding extrinsic noise lies at the heart of understanding cell-population heterogeneity.

So far, identifying the sources of gene expression noise from transcriptomic measurements alone has proven difficult Paulsson, 2004; Pedraza and Paulsson, 2008. The fundamental hindrance lies in the fact that single-cell RNA sequencing, which provides most of the available data, is destructive, so that datasets reflect samples from across a population, rather than samples taken repeatedly from the same cell. As temporal information is lost in such measurements Komorowski et al., 2011, it may be impossible to distinguish temporal variability within individual cells (e.g. burstiness), from ensemble variability across the population (i.e. extrinsic noise). A number of numerical and experimental studies have suggested this confounding effect Jones and Elf, 2018; Jones et al., 2014; Zopf et al., 2013, showing that systems with intrinsic noise alone exhibit behaviour that is indistinguishable from systems with both extrinsic and intrinsic noise. This is examined more formally in Ham et al., 2020a, where we show that the moment scaling behaviour and transcript distribution may be indistinguishable from situations with purely intrinsic noise. The limitations in inferring dynamics from population data are becoming increasingly evident, and a number of studies that seek to address some of these problems have emerged Skinner et al., 2016; Gorin and Pachter, 2020a.

Here we provide adetailed analysis of the extent to which sources of variability are identifiable from population single-cell omics data. We are able to prove rigorously that it is in general impossible to identify the sources of variability, and consequently, the underlying transcription dynamics, from observed transcript abundance distributions alone. Systems with intrinsic noise alone can always present identically to similar systems with extrinsic noise. For practical purposes, the effect does not appear to depend on the precise choice of distribution, but holds more generally. Moreover, we demonstrate that extrinsic noise invariably distorts the apparent degree of burstiness of the underlying system: data which seems ‘bursty’ is not necessarily generated by a bursty process, if there is cell-to-cell variability across the population. This is a stronger non-identifiability result than has previously been obtained Ingram et al., 2008; Shahrezaei and Swain, 2008a; Khammash, 2009; Munsky et al., 2009; Komorowski et al., 2011, and has important ramifications for our analysis of experimental data: we cannot assess causes of transcriptional variability, if we do not simultaneously assess cell-to-cell variability in the transcriptional machinery. Our results highlight (in fact prove mathematically) the requirement for additional information, beyond the observed copy number distribution, in order to constrain the space of possible dynamics that could give rise to the same distribution.

This seemingly intractable problem can at least partially be resolved with a brilliantly simple approach: the dual-reporter method Swain et al., 2002. In this approach, noise can be separated into extrinsic and intrinsic components, by observing correlations between the expression of two independent, but identically distributed fluorescent reporter genes. Dual-reporter assays have been employed experimentally to study the noise contributed by both global and gene-specific effects Elowitz et al., 2002; Raser and O'Shea, 2005; Raj et al., 2006. A particular challenge, however, is that dual reporters are rarely identically regulated Raj et al., 2006; Quarton et al., 2020, and are not straightforward to set up in every system. More recently, it has been shown that the independence of dual reporters cannot be guaranteed in some systems Naik et al., 2021. As a result, there have been efforts in developing alternative methods for decomposing noise Quarton et al., 2020; Singh, 2014; Lin and Amir, 2021. Here we develop a widely applicable generalisation (and simplification) of the original dual-reporter approach Swain et al., 2002. We demonstrate that non-identical and not-necessarily independent reporters can provide an analogous noise decomposition. Our result shows that measurements taken from the same biochemical pathway (e.g. mRNA and protein) can replace dual reporters, enabling the noise decomposition to be obtained from a single gene. This completely circumvents the requirement of conditionally independent and identically regulated reporter genes. The results obtained from our ‘pathway-reporter’ method are also borne out by stochastic simulations, and compare favourably with the dual-reporter method. In the case of constitutive expression, the results obtained from our decomposition are identical to those obtained from dual reporters. For bursty systems, we show that our approach provides a satisfactorily close approximation, except in extreme cases. Our methodology is verified mathematically for the most common models of gene transcription, but holds independently of the specific nature of the gene expression dynamics, as we verify in silico across a range of more detailed models.

Materials and methods

A simple model for stochastic mRNA dynamics is the Telegraph model: a two-state model describing promoter switching, transcription, and mRNA degradation. In this model, all parameters are fixed, and gene expression variability arises due to the inherent stochasticity of the transcription process. As discussed above, this process will often be influenced by extrinsic sources of variability, and so modifications to account for this additional source of variability are required.

The telegraph model

Request a detailed protocol

The Telegraph model was first introduced in Ko, 1991, and has since then been widely employed in the literature to model bursty gene expression in eukaryotic cells Bahar Halpern et al., 2015; So et al., 2011; Suter et al., 2011; Larsson et al., 2019. In this model, the gene switches probabilistically between an active state and an inactive state, at rates λ (on-rate) and μ (off-rate), respectively. In the active state, mRNAs are synthesised according to a Poisson process with rate K, while in the inactive state, transcription does either not occur, or possibly occurs at some lower Poisson rate, K0K. Degradation of mRNA molecules occurs independently of the gene promoter state at rate δ. Figure 1A shows a schematic of the Telegraph model. Throughout the discussion here, we will rescale all parameters of the Telegraph model by the mRNA degradation rate, so that δ=1. The steady-state distribution for the mRNA copy number can be explicitly calculated as Peccoud and Ycart, 1995

(1) p~T(n;θ)=Knλ(n)n!(μ+λ)(n)F11(λ+n,λ+μ+n,-K).

Here, θ denotes the parameter vector (μ,λ,K,δ), the function F11 is the confluent hypergeometric function Abramowitz and Stegun, 1965, and, for real number x and positive integer n, the notation x(n) abbreviates the rising factorial of x (also known as the Pochhammer function). Throughout, we refer to the probability mass function p~T(n;θ) as the Telegraph distribution with parameters θ.

Modeling the effects of both intrinsic and extrinsic noise.

(A) A schematic of the Telegraph process, with nodes A (active) and I (inactive) representing the state of the gene. Transitions between the states A and I occur stochastically at rates μ and λ, respectively. The parameter K is the mRNA transcription rate, and δ is the degradation rate. (B) The compound model incorporates extrinsic noise by assuming that parameters θ of the Telegraph model vary across an ensemble of cells, according to some probability distribution f(θ;η). (C) Variation in the parameters across the cell population leads to greater variability in the mRNA copy number distribution.

Constitutive gene expression is a limiting case of the Telegraph model, which arises when the off-rate μ is 0, so that the gene remains permanently in the active state. In this case, the Telegraph distribution simplifies to a Poisson distribution with rate K; the distribution Pois(K).

At the opposite extreme is instantaneously bursty gene expression in which mRNA are produced in very short bursts with potentially prolonged periods of inactivity in between. This mode of gene expression has been frequently reported experimentally, particularly in mammalian genes Raj et al., 2006; Bahar Halpern et al., 2015; Suter et al., 2011; Larsson et al., 2019. Transcriptional bursting may be treated as a limit of the Telegraph model, where the off-rate, μ, tends to infinity, while the on-rate λ remains fixed. In this limit, it can be shown Jones and Elf, 2018; Ham et al., 2020a that the Telegraph distribution converges to the negative binomial distribution NegBin(λ,Kμ+K).

The compound distribution

Request a detailed protocol

We can account for random cell-to-cell variation across a population by way of a compound distribution Ham et al., 2020b 

(2) q~(n;η)=p~(n;θ)f(θ;η)𝑑θ,

where p~(n;θ) is the stationary probability distribution of a system with fixed parameters θ and f(θ;η) denotes the multivariate distribution for θ with hyperparameters η. Often we will take p~(n;θ) to be the stationary probability distribution of the Telegraph model ((1)), and refer to (2) as the compound Telegraph distribution. Sometimes p~(n;θ) will be the Poisson distribution or the negative binomial distribution, depending on the underlying mode of gene activity. Figure 1B gives a pictorial representation of the compound distribution.

The compound distribution is valid in the case of ensemble heterogeneity, that is, when parameter values differ between individual cells according to the distribution f(θ;η), but remain constant over time Filippi et al., 2016. This model is also a valid approximation for individual cells with dynamic parameters, provided these change sufficiently slower than the transcriptional dynamics Lenive et al., 2016. In general, the compound Telegraph distribution q~(n;η) will be more dispersed than a Telegraph distribution to account for the uncertainty in the parameters; see Figure 1C. Such dispersion is widely observed experimentally, and as demonstrated in Ham et al., 2020a, reflects the presence of extrinsic noise.

In the context of gene expression, it has been shown experimentally that some of the primary sources of extrinsic noise have an autocorrelation time comparable to the cell cycle Rosenfeld et al., 2005. It is these slow changes in variability that justify the assumptions of the compound model. Typical sources of extrinsic variability for each parameter in the Telegraph model are summarised in Table 1 of Ham et al., 2020a. A further significant source of heterogeneity arises from the differences in cell-cycle phases across the population Skinner et al., 2016. Such effects have been shown to obscure the precise underlying transcriptional dynamics Zopf et al., 2013; Huh and Paulsson, 2011, and impede the inference of transcriptional parameters from experimental data Beentjes et al., 2020. The compound model we consider here is able to capture this variability, provided that the parameter change within each cell phase is relatively slow, and any dynamic parameter changes during the transition between phases can be considered as ephemeral. A more explicit treatment of the mechanisms and changes during the cell cycle is challenging to study analytically, and theoretical modelling is only in its infancy Cao and Grima, 2020; Beentjes et al., 2020; Perez-Carrasco et al., 2020. Later, we verify our proposed noise decomposition on detailed models of gene transcription, incorporating salient features of the cell-cycle, such as gene replication, dosage compensation, binomial partitioning of products due to cell division, and cell-cycle length variability.

Table 1
Summary of the non-identifiability results.

in lines 1, 3, and 5 are our contributions, while the remaining representations (lines 2 and 4) are known and can be obtained as special cases of our results. Note that here we use Tele(λ,μ,K) to denote a Telegraph distribution with parameters λ,μ,K. In lines 3 and 4, the parameter β>0 can be chosen freely and determines the mean burst intensity in the resulting compound system. In line 5, the parameters b,θ>0 are again mean burst intensities, and b can be chosen freely in the determination of the distribution of θ.

Copy no. distribution q~(n;η)Underlying distribution p~(n;θ)Noise distribution f(θ;η)
Tele(λ,μ,K)Tele(λ,μ,K)KBetaK(λ+μ,μμ)
Tele(λ,μ,K)Pois(K)KBetaK(λ,μ)
NegBin(λ,ββ+1)Tele(λ,μ,K)KGamma(λ+μ,β)
NegBin(λ,ββ+1)Pois(K)KGamma(λ,β)
NegBin(λ,bb+1)NegBin(λ,θθ+1)θBetaPrimeb(λλ,λ)

Results

Identifiability considerations

Decoupling the effects of extrinsic noise from experimental measurements has been notoriously challenging. In the context of (2), the distribution f(θ;η) reflects the population heterogeneity, but experimental data provides samples only of q~(n;η). How much can we deduce of the underlying dynamics (that is, p~(n;θ)), and the population heterogeneity (f(θ;η)), from measurements of transcripts from across the cell population (q~(n;η))?

Of course, even though we may be able to infer the parameter η from experimental data, the expression p~(n;θ) is really a family of distributions, parameterised by θ. This presents two fundamental challenges. The first is the possibility that there are different families of distributions p~(n;θ) that can yield the same compound distribution, q~(n;η), but which are generated by different mechanisms, that is noise distributions, f(θ;η). The second, perhaps more subtle, challenge is that, even for a fixed family of distributions p~(n;θ) it may be possible that different choices of the noise distribution f(θ;η) could still yield the same compound distribution q~(n;η). In these situations, we cannot hope to infer a unique solution for the noise distribution. This belongs to the important class of identifiability problems Villaverde, Villaverde and Banga, 2017; it has important ramifications for the interpretability of parameter estimates obtained from experimental data Ingram et al., 2008. Indeed, if two or more model parameterisations are observationally equivalent (in this case, in the form of the transcript abundance distribution q~(n;η)), then not only does this cast doubt upon the suitability of the model to represent (and subsequently predict) the system, but it also obstructs our ability to infer mechanistic insight from experimental data.

An example of the first identifiability problem arises from a well-known example of a compound distribution, (2): when f(θ;η) is a gamma distribution and p~(n;θ) is a Poisson distribution, corresponding to constitutive gene expression, the resulting compound distribution q~(n;η) is a negative binomial distribution Greenwood and Yule, 1920. But this is the same distribution as that arising from instantaneously bursty gene expression Ham et al., 2020a. Such identifiability instances may be circumvented if there is confidence in the basic mode of gene activity, that is, if there is reason to believe that a gene is not constitutively active, for example. We find, however, that there are numerous instances that can present insurmountable identifiability problems.

Bursty gene expression

We first observe that any Telegraph distribution with fixed parameters can be identically obtained from a Telegraph distribution with parameter variation. As shown in the supplementary material (Appendix Pathway dynamics can delineate the sources of transcriptional noise in gene expression), any Telegraph distribution p~T(n;λ,μ,K) can be written as, 

(3) p~T(n;λ,μ,K)=0Kp~T(n;λ,μ,t)fK(t;λ+μ,μ-μ)𝑑t,

where μ<μ and fK(t;λ+μ,μ-μ) is the probability density function for a scaled beta distribution BetaK(λ+μ,μμ) with support [0,K]. Thus, any Telegraph distribution can be obtained by varying the transcription rate parameter on a narrower Telegraph distribution (i.e. with a smaller off-rate) according to a scaled beta distribution. Figure 2A (top panel) compares the representation obtained in (3) with the corresponding fixed-parameter Telegraph distribution for two different sets of parameters. When μ=0 the representation given in (3) simplifies to the well-known Poisson representation of the Telegraph distribution in terms of the scaled beta distribution Srividya et al., 2009.

Accuracy of our integral representations for the Telegraph and negative binomial distribution.

(A) For each of the results in (3 - 5), we compare the (fixed-parameter) Telegraph and negative binomial distributions with their respective compound representations for two different sets of parameter values. The top panel (pink) shows comparisons for (3), with parameter values (left) λ=2, μ=12, K=100, μ=3, and KBetaK(5,9), and (right) λ=1, μ=20, K=100, μ=2 and KBetaK(3,18). The middle panel (green) gives comparisons for (4), with parameter values (left) λ=10, β=2, μ=2 and KGamma(12,2) and (right) λ=1, β=1, μ=2 and KGamma(3,1). The bottom panel (coral) gives comparisons for (5). The parameter values (left) are λ=10, λ=15 and c=2 and (right) are λ=2, λ=5 and c=3. (B) The top figure compares a Telegraph(2,4,60) distribution with samples from a compound Telegraph distribution with normal noise Norm(37,10) on the transcription rate parameter. The middle figure compares a NegBin(5,0.5) with samples from a compound Telegraph distribution with normal noise Norm(5.5,2.3) on the transcription rate parameter. The bottom figure compares a NegBin(5,1) distribution with samples from a compound negative binomial distribution with normal noise Norm(2.3,0.6) on the burst intensity parameter.

Instantaneously bursty gene expression

The previous result extends to instantaneously bursty systems. The copy number distribution of an instantaneously bursty system can be obtained from both bursty and instantaneously bursty dynamics, provided that there is appropriate parameter variation. The supplementary material contains the relevant derivations; refer to Appendix Pathway dynamics can delineate the sources of transcriptional noise in gene expression. In the following, we let p~NB(n;r,β) denote the probability mass function of a NegBin(r,ββ+1) distribution, where β>0. Then for any negative binomial distribution of the form NegBin(λ,ββ+1) we have,

(4) p~NB(n;λ,β)=0p~T(n;λ,μ,t)f(t;λ+μ,β)𝑑t,

where f(t;λ+μ,β) is the probability density function of a Gamma(λ+μ,β) distribution. This result generalises the aforementioned well-known representation of the negative binomial distribution Greenwood and Yule, 1920, which corresponds to the particular case of μ=0. In Figure 2A (middle panel), we compare the representation obtained in (4) with the corresponding fixed-parameter negative binomial distribution for two different sets of parameters.

We also obtain the following representation for a negative binomial distribution in terms of a scaled beta prime distribution,

(5) p~NB(n;λ,b)=1bp~NB(n;λ,θ)fb(bθ-1;λ-λ,λ)𝑑θ,

where fb(bθ-1;λ-λ,λ) is the probability mass function of a scaled beta prime BetaPrimeb(λλ,λ) distribution, where b>0 and λ>λ. This equivalently corresponds to scaled beta noise Betab(λλ,λ) on the inverse of the expected burst intensity. Thus, the distribution of any instantaneously bursty system with mean burst intensity b can be obtained from one with greater burst frequency, by varying the mean burst intensity θ according to a shifted beta prime distribution. Figure 2A (bottom panel), compares the representation obtained in (5) with the associated fixed-parameter negative binomial distribution for two different sets of parameters.

An exception: constitutive expression

It has long been known Feller, 1943 that a compound Poisson distribution uniquely determines the compounding distribution. In the context of (2), this means the full extrinsic noise distribution f(θ,η) is identifiable from q~(n;η). As we will demonstrate in related work (Ham et al., in preparation) , it is therefore possible to extract the extrinsic noise distribution, f(θ,η), from transcript copy number measurements.

Implications for parameter inference

Estimates of kinetic parameters from experimental data suggest that gene expression is often either bursty or instantaneously bursty (i.e., μλ). In turn, the assumption that gene-inactivation events occur far more frequently than gene-activation events is often used to derive other models of stochastic mRNA dynamics Jia and Grima, 2020; Cao and Grima, 2020; Beentjes et al., 2020. The representations given in (3 – 5), however, show that both estimating parameters and the underlying dynamics from the form of the copy number distribution alone can be misleading. Noise on the transcription rate will invariably produce copy number data that is suggestive of a more bursty model. To illustrate this, consider an example in which the underlying process is a (mildly) bursty Telegraph system with distribution p~T(n;2,3,K). Now assume that noise on the transcription rate parameter K follows the scaled Beta distribution on the interval [0,100] with α=λ+μ=2+3=5 and β=μ-μ=12-3=9. The shape of this noise distribution closely resembles a slightly skewed Gaussian distribution, with the majority of transcription rates between around 10 and 60. This noise on the transcription rate K within the Telegraph system p~T(n;2,3,K) will present identically to the significantly burstier system p~T(n;2,12,100).

It is of practical importance to recognise that, while the non-identifiability results (summarised in Table 1) are dependent on specific noise distributions, for practical purposes any similar distribution will produce a similar effect. To demonstrate this, we replace the various noise distributions required for the representations in (3 – 5), with suitable normal distributions truncated at 0. In each case, we sample 1000 data points from the corresponding compound distribution, and compared this with the associated fixed-parameter copy number distribution. The results are shown in Figure 2B. The truncated normal distribution is not chosen on the basis of biological relevance, but rather to demonstrate that even a symmetric noise distribution (except for truncation at 0) produces qualitatively similar results to the distributions used in the precise non-identifiability results. In every case, the effect of varying the transcription rate or burst intensity parameter according to a unimodal noise distribution is to produce copy number distributions that are generally consistent with systems that appear burstier.

Finally, we note that our results explain previous empirical observations that static measurements of mRNA are not always sufficient to infer the underlying dynamics of gene activity. Skinner et al., 2016 address some of these limitations by quantifying both nascent and mature mRNA in individual cells, as well incorporating cell-cycle effects into their analysis of two mammalian genes. A more developed treatment of model identifiability is given in Gorin and Pachter, 2020a, where it shown how a stochastic model incorporating the downstream processing of mRNA can be used to distinguish a particular instance of non-identifiability. More specifically, the authors consider the non-identifiability problem noted in Ham et al., 2020a, arising from the Gamma-Poisson compound representation of the negative binomial distribution; a particular case of (4) above. Despite identical distributions at the nascent level, the marginal distributions of the processed (mature) mRNA are found to be substantially different. It is likely that a similar analysis will be valuable in the context of the other identifiability problems we have given in Table 1. In the next section, we take a more general approach to resolving non-identifibiality and exploit the properties of complex gene expression dynamics to determine not only the presence of extrinsic noise, but also estimate its magnitude.

Resolving non-identifiability

The results of the previous section show that additional information, beyond the observed copy number distribution, is required to constrain the space of possible dynamics that could give rise to the same distribution. One way to narrow this space of possibilities, is to determine the intrinsic and extrinsic contributions to the total variation in the system.

The dual-reporter method

The total gene expression noise, as measured by the squared coefficient of variation η2, can be decomposed exactly into a sum of intrinsic and extrinsic noise contributions Swain et al., 2002. The decomposition applies to dynamic noise Hilfinger and Paulsson, 2011, and generalises to higher moments in Hilfinger et al., 2012. Sets of dual reporters at multiple levels of the transcriptional pathway has been shown to achieve a finer breakdown of noise into subcategories Bowsher and Swain, 2012. As shown in Hilfinger and Paulsson, 2011, the noise decomposition is equivalent to the normalised Law of Total Variance (Ross, 2014). Indeed, if X is the random variable denoting the number of molecules of a certain species (e.g. mRNA or protein) in a given cell, then we can decompose the total noise by conditioning X on the state 𝐙=(Z1,,Zn) of the environmental variables Z1,…,Zn,

(6) ηX2=E(Var(X;Z))E(X)E(Y)+Var(E(X;Z))E(X)E(Y)ηint2+ηext2.

It has been shown Swain et al., 2002; Hilfinger and Paulsson, 2011 that if X1 and X2 are random variables denoting the expression levels of independent (conditional on 𝐙) and identically distributed gene reporters, then the extrinsic noise contribution ηext2 in (6) can be identified by the normalised covariance between X1 and X2,

(7) ηext2=Cov(X1,X2)E(X1)E(X2).

Decomposing noise with non-identical reporters

The dual-reporter method requires distinguishable measurements of transcripts or proteins from two conditionally independent and identically distributed reporter genes integrated into the same cell. In practice, however, dual reporters rarely have identical dynamics, which is widely considered to be a significant challenge to interpreting experimental results Quarton et al., 2020. We show that, under certain conditions, the decomposition in (6) can alternatively be obtained from non-identically distributed and not-necessarily independent reporters.

Our result relies on the observation that the covariance of any two variables can be decomposed into the expectation of a conditional covariance and the covariance of two conditional expectations (the Law of Total Covariance Ross, 2014). If X and Y denote, for example, the numbers of molecules of two chemical species (eg. mRNA and protein) in a given cell, then the covariance of X and Y can be decomposed by conditioning on the state 𝐙=(Z1,,Zn) of the environmental variables Z1,…,Zn,

(8) Cov(X,Y)=E(Cov(X,Y;Z))intrinsic+Cov(E(X;Z),E(Y;Z))extrinsic.

We will see that in many cases of interest the random variable E(X;Z) (as a function of 𝐙) splits across common variables with E(Y;Z). By this we mean that E(X;Z)=f(ZX)hX(Z) and E(Y;Z)=g(ZY)hY(Z), where ZX are the variables of 𝐙 that appear in E(X;Z) but not in E(Y;Z), and dually, 𝐙Y are those in E(Y;Z) that are not in E(X;Z). The variables 𝐙 are those variables from 𝐙 not in ZXZY. In these cases, the component of Cov(X,Y) that is contributed by the variation in 𝐙 (the extrinsic component) may be written as the covariance of the functions hX(𝐙) and hY(𝐙). Conveniently, in the cases of interest here, the two functions hX and hY coincide, and this is the form we use in the following decomposition principle. The supplementary material (Appendix A) contains the proof of this result.

The noise decomposition principle (NDP)

Assume that there are measurable functions f, g, and h such that E(X;Z) and E(Y;Z) split across common variables by way of E(X;Z)=f(ZX)h(Z) and E(Y;Z)=g(ZY)h(Z). Then, provided that the variables Z1,,Zm are mutually independent, the normalised covariance of E(X;Z) and E(Y;Z) will identify the total noise on h(𝐙) (i.e. ηh(𝐙)2).

As we show in the next section, there are many situations where the random variable E(X;Z) is precisely the common part of E(Y;Z) and E(X;Z) (i.e., h(Z)=E(X;Z)), and the normalised intrinsic contribution to the covariance is either zero or negligible. In these cases, the normalised covariance of X and Y will identify precisely the extrinsic noise contribution ηext2 to the total noise ηX2. To see this, consider the situation where E(Y;Z)=f(ZY)E(X;Z). Then provided f(ZY) and E(X;𝐙) are independent random variables, the extrinsic contribution to the covariance of X and Y is given by,

(9) Cov(E(X;Z),E(Y;Z))=Cov(E(X;Z),f(ZY)E(X;Z))=E(f(ZY))Cov(E(X;Z),E(X;ZY))=E(f(ZY))Var(E(X;Z)).

If the normalised intrinsic contribution to the covariance is either zero or is negligible, it follows from (8) that

(10) Cov(X,Y)E(X)E(Y)Var(E(X;Z))E(X)2=ηext2.

Thus, under certain conditions, measuring the covariance between two non-identically distributed and not-necessarily independent reporters can replace dual reporters.

The pathway-reporter method

We show that for some reporters X and Y belonging to the same biochemical pathway, the covariance of X and Y continues to identify the extrinsic, and subsequently intrinsic, noise contributions to the total noise. The basis of the pathway-reporter method depends on the emergent covariances between the various species (e.g. nascent/mature mRNA and protein) in the gene expression pathway. Qualitatively, this effect can be seen in Figure 3, which compares simulated sample distributions of a simple four-stage model of gene transcription (refer to the model 𝐌4 below) in the case of moderate extrinsic noise to the case with no extrinsic noise. The plots are the bivariate distributions for nascent-mature, nascent-protein, and mature-protein levels, respectively. This will be made more precise below, where we find that it is possible to extract quantitative information about the extrinsic noise distribution itself from this data.

A comparison of joint distributions in the case of moderate extrinsic noise and no extrinsic noise.

The plots are generated from a three-stage model of gene transcription, incorporating the production of nascent mRNA, mature mRNA and protein. Details of the model can be found in Figure 4 (model 𝐌4) and the associated text. The top panel shows nascent-mature, nascent-protein and mature-protein joint distributions in the case of extrinsic noise, while the bottom panel displays the corresponding plots in the case of no extrinsic noise. Extrinsic noise produces a visibly more correlated joint distribution, which forms the basis of the pathway-reporter method.

Stochastic models of gene expression.

(A) The model 𝐌1 is the simplest model of mRNA maturation. Here, nascent (unspliced) mRNA are shown in red/blue wavy lines; the blue segments represent introns and the red segments represent the exons. Nascent mRNA are synthesised at the rate KN, and spliced into mature mRNA (blue wavy lines) at rate KM. Degradation of the mature mRNA occurs at rate δM. The model 𝐌2 is the well-known two-stage model of gene expression. The model 𝐌3 is the extension of the two-stage model to include promoter switching. The nodes A (active) and I (inactive) represent the state of the gene, with transitions between states occurring at rates λ and μ. The remaining parameters are the same as those in the model 𝐌2. The model 𝐌4 extends the model 𝐌3 by incorporating mRNA maturation. Here, KN is the transcription rate parameter, and KM is the maturation rate. All other parameters are the same as in 𝐌3. (B) Time series simulation of the copy number and activity state of a gene modelled by 𝐌4. For ease of visualisation, the parameters were artificially chosen as λ=2, μ=2.5, KN=40, KM=4, Kp=4 and δp=1, with all parameters scaled relative to δm=1. (C) As λ approaches 0, we see a higher correlation in the copy numbers of nascent mRNA, mature mRNA and protein. Again, the parameters are artificially chosen to be λ=0.1, μ=2.5, KN=80, KM=4, Kp=4 and δp=1, with all parameters scaled relative to δm=1.

Throughout this section, we assume that extrinsic noise sources act independently i.e., the environmental variables Z1,,Zn of 𝐙 are mutually independent. Additionally, our modelling focuses only on a single gene copy, although the same analysis applies to multiple but indistinguishable gene copies; we refer to the supplementary material (Appendix B) for more details.

Measuring noise from a constitutive gene

We consider first the simplest case where the underlying process is constitutive. We begin with a stochastic model of mRNA maturation, which we denote by 𝐌1; Figure 4A (top left) gives a schematic of the constitutive model. In this model, the gene continuously produces nascent mRNA according to a Poisson process at constant rate KN, which are subsequently spliced into mature mRNA according at rate KM. Degradation of mature mRNA occurs as a first-order Poisson process with rate δM. The model 𝐌1, together with its extensions, has been considered in a number of recent studies Gorin and Pachter, 2020a; Cao and Grima, 2020; La Manno et al., 2018; Gorin and Pachter, 2020b; Bergen2020, and has a known solution for the stationary joint probability distribution Jahnke and Huisinga, 2007 given by,

(11) p~(n,m;θ)=e-KnKM(KNKM)nn!e-Knδ(KNδM)mm!,

where n is the number of nascent mRNA, m is the number of mature mRNA, and the parameter θ=(KN,KM,δM). We use XN to denote the number of nascent mRNA, XM the number of mature mRNA produced from the same constitutive gene, and 𝐙=(KN,KM,δM). To simplify notation, we abbreviate the variables in 𝐙XN as 𝐙N, and similarly for 𝐙XM. It follows immediately from (11) that XN and XM are independent conditional on 𝐙, and so the intrinsic contribution to the covariance of XN and XM (the first term of (8)) is 0. It is also easy to see from (11) that E(XN;𝐙)=f(𝐙N)KN and E(XM;𝐙)=g(𝐙M)KN, where f(𝐙N)=1KM and g(𝐙M)=1δM. Since the extrinsic noise sources are assumed to act independently, it follows that the Noise Decomposition Principle (NDP) of the previous section holds. We then have that Cov(XN,XM)=ηKN2, where ηKN2 is the total noise on the transcription rate parameter KN. Thus, measuring Cov(XN,XM) can replace dual reporters to decompose the gene expression noise at the transcriptional level.

To support our mathematical results, we simulate the model 𝐌1 subject to parameter variation using the stochastic simulation algorithm (SSA). Table 2 compares the extrinsic noise contributions found from various simulations with the corresponding theoretical values. In each simulation, the degradation rate δm is fixed at 1, with the other parameters scales accordingly. The maturation rate KM is sampled according to a Gamma(8,0.0125) distribution, which has coefficient of variation 0.125. We consider different noise distributions on KN, producing a range of noise strengths. Our theory predicts that pathway-reporters will identify the total noise on KN. Overall, we observe an excellent agreement between the results obtained by the pathway-reporter method, the dual-reporter method and the theoretical noise. There is consistently slightly more variation in the pathway-reporter results compared with the dual-reporter results.

Table 2
A comparison of the pathway-reporter method and the dual-reporter method for constitutive expression under the model 𝐌1.

Here, PR (NM) gives the results of the nascent and mature pathway reporters, while DR (Mat) gives the results of dual reporters calculated from the mature mRNA. We considered noise on both the transcription rate (KN) and the maturation rate (KM). The decay rate is fixed at one, with the other parameters scaled accordingly. In each case, the maturation rate KM is varied according to a Gamma(8,1.25) distribution, which has coefficient of variation 0.125. The values given are the average of 100 simulations, each calculated from 500 copy number samples, and the errors are ± one standard deviation. Our theory predicts that pathway-reporters will identify the noise on the nascent transcription rate KN (ηext2). The noise distribution parameters are chosen to produce an average nascent mRNA copy number of approximately five and an average mature mRNA copy number of approximately 50.

TheorySimulation
(r)1-2 ηext2Noise (KN)Pr (NM)DR (Mat)
0.00KN=500.00±0.010.00±0.00
0.10Beta133.3˙(6,10.5)0.10±0.010.10±0.01
0.20Gamma(5,10)0.20±0.020.20±0.01
0.50Beta300(1.5,7.5)0.49±0.040.50±0.03
Table 2—source data 1

This is an Excel spreadsheet containing the data used to produce the final values in Table 2.

https://cdn.elifesciences.org/articles/69324/elife-69324-table2-data1-v2.xlsx

To explore the pathway-reporter method more widely, we consider 60 different parameter combinations to produce a range of mean copy numbers consistent with those reported experimentally. We also consider different noise distributions taken from the scaled Beta distribution family in order to produce a range of noise strengths; refer to Supplementary file 1. The pathway-reporter method performs favourably compared to the dual-reporter method calculated from mature mRNA, and consistently outperforms the dual-reporter method on nascent mRNA.

Next, we consider the simplest stochastic model of gene expression that includes both mRNA and protein dynamics: the well-known ‘two-stage model’ of gene expression, which, together with its three-stage extension to include promoter switching has been widely studied Raser and O'Shea, 2005; Raj et al., 2006; Thattai and van Oudenaarden, 2001; Friedman et al., 2006; Singh and Hespanha, 2007; Shahrezaei and Swain, 2008a; Bokes et al., 2012; Molina et al., 2013. We denote this model by 𝐌2; see Figure 4A (top right) for a schematic of this model. In this model, mRNA are synthesised according a Poisson process at rate Km, which are then later translated into protein at rate Kp. Degradation of mRNA and protein occur as first-order Poisson processes with rates δm and δp, respectively. If Xm denotes the number of mRNA, Xp the number of proteins produced from the same constitutive gene, and if 𝐙=(Km,Kp,δm,δp), then the stationary means and covariance are given by Thattai and van Oudenaarden, 2001; Singh and Hespanha, 2007:

(12) E(Xm;Z)=Kmδm,E(Xp;Z)=KpδpKmδmandCov(Xm,Xp;Z)=KmKpδm(δm+δp).

It is easily verified that E(Xp;Z)=f(Zp)E(Xm;Z), where f(𝐙p)=Kpδp. Thus, it follows from the NDP that the normalised contribution of Cov(Xm,Xp) contributed by 𝐙 will identify the extrinsic noise contribution to the total noise on Xm. Now, if we assume that δm is fixed across the cell-population, and all parameters are scaled so that δm=1, we have the following expression for the intrinsic contribution to the covariance of Xm and Xp; refer to the supplementary material (Appendix B) for details.

(13) E(Cov(Xm,Xp;Z))E(Xm)E(Xp)=αE(Km),where α=E(1/(δp+1))E(1/δp).

Since mRNA tends to be less stable than protein, we have δp<1, and often δp1Bernstein et al., 2002; Schwanhäusser et al., 2011. So, we can expect α1. Further, for many genes we can expect the number of mRNA per cell (Km) to be in the order of tens, so 1/E(Km)<1. It follows that E(Cov(Xm,Xp;Z)) 1, so that Cov(Xm,Xp) will closely approximate the extrinsic noise at the transcriptional level.

We test our theory using stochastic simulations of the model 𝐌2 subject to parameter variation. Table 3 gives a comparison of the results of the mRNA-protein reporters and dual reporters. In each case, we varied Kp according to a Gamma(5,0.4) distribution and δp according to a Gamma(8,0.125) distribution; the corresponding noise strengths are 0.20 and 0.125, respectively. We consider different noise distributions on Km, which produce a range of noise strengths, and the noise distribution parameters are selected to produce a mean mRNA of approximately 50 and a mean number of approximately 1000 proteins in each simulation. As our theory predicts, the mRNA-protein reporters identify the extrinsic noise contribution to the total noise on Xm. Again, we can see from Table 3 that there is excellent agreement between the results of the pathway reporters and the dual reporters, with slightly more variation in the pathway-reporter results. A larger exploration of the parameter space reveals similar results; these are provided in Supplementary file 1. Thus, despite mRNA and protein numbers not being strictly independent, they can, for practical purposes, replace dual reporters to decompose the noise at the transcriptional level.

Table 3
A comparison of the pathway-reporter method and the dual-reporter method for constitutive expression under the model 𝐌2.

Here PR (MP) gives the results of the mRNA-protien pathway reporters, while DR (Mat) gives the results of dual reporters calculated from the mature mRNA. We considered noise on the transcription rate (Km), the protein synthesis rate (Kp), and the protein decay rate (δp). The mRNA decay rate is fixed at one. In each case, we varied Kp according to a Gamma(5,0.4) distribution and δp according to a Gamma(8,0.125) distribution; the corresponding noise strengths are 0.20 and 0.125, respectively. We considered different noise distributions on Km, which produce a range of noise strengths. The noise distribution parameters are selected to produce a mean mRNA of approximately 50 and a mean number of approximately 1000 proteins in each simulation. The values given are the average of 100 simulations, each calculated from 500 copy number samples, and the errors are ± one standard deviation. As our theory predicts, the mRNA-protein reporters identify the noise on the transcription rate parameter Km (ηext2).

TheorySimulation
(r)1-2 ηext2Noise (Km)Pr (MP)DR (Mat)
0.00Km=500.00±0.010.00 ± 0.00
0.10Beta133.3˙(6,10.5)0.10±0.010.10 ± 0.01
0.20Gamma(5,10)0.20±0.020.20±0.01
0.50Beta300(1.5,7.5)0.51±0.040.50±0.03
Table 3—source data 1

This is an Excel spreadsheet containing the data used to produce the final values in Table 3.

https://cdn.elifesciences.org/articles/69324/elife-69324-table3-data1-v2.xlsx

We note that both the pathway-reporter (nascent-mature or mature-protein) and dual-reporter methods show slower convergence when copy numbers are low. Pathway reporters usually show fractionally slower convergence and fractionally more variation than a dual reporter, as suggested by the standard deviations in Tables 2 and 3. A more detailed exploration of convergence is given in the supplementary material (Appendix C).

Measuring noise from a facultative (bursty) gene

The most common mode of gene expression that is reported experimentally is burstiness Jones and Elf, 2018; Raj et al., 2006; Bahar Halpern et al., 2015; Suter et al., 2011; Larsson et al., 2019; Golding et al., 2005, in which mRNA are produced in short bursts with periods of inactivity in between. One example is gene regulation via repression, which naturally leads to periods of gene inactivity. Here, we consider a four-stage model of bursty gene expression, which incorporates both promoter switching and mRNA maturation; we denote this model by 𝐌4; see Figure 4A (bottom left). This model has recently been considered in Cao and Grima, 2020, where the marginal distributions are solved in some limiting cases. In this model, the gene switches probabilistically between an active state (A) and an inactive state (I), at rates λ (on-rate) and μ (off-rate), respectively. In the active state, nascent mRNA is synthesised according to a Poisson process at rate KN, while in the inactive state transcription does not occur. Nascent mRNA is spliced into mature mRNA at rate KM; this in turn is later translated into protein, with rate KP. Degradation of mRNA and protein occur independently of the promoter state at rates δM and δP, respectively.

For this model, we have three natural candidates for pathway reporters: (a) nascent and mature mRNA (b) mature mRNA and protein, and (c) nascent mRNA and protein reporters. Below we show that nascent mRNA–protein reporters yield consistently good estimates of the extrinsic noise contribution ηext2 to the total gene expression noise, while nascent–mature and mature RNA–protein reporters are reliable in some restricted cases. We begin by showing that each of the reporter pairs (a), (b), and (c) satisfy the Noise Decomposition Principle. We then demonstrate computationally, that despite a lack of independence between these reporter pairs, the pathway-reporter method can still be used to decompose the total gene expression noise at the transcriptional level. Throughout, we let XN denote the number of nascent mRNA, we let XM denote the number of mature mRNA, and let XP denote the number of proteins produced from the same gene. We also let 𝐙={λ,μ,KN,KM,KP,δM,δP}.

Following Raj et al., 2006, we assume that the transcription rate KN is large relative to the other parameters. We further assume that the maturation rate KM is large (i.e. KM>δM), which is supported by experiments Cao and Grima, 2020. Then, using the results of Raj et al., 2006 and arguments similar to those in Cao and Grima, 2020, it can be shown that the stationary averages for the nascent mRNA, mature mRNA and protein levels are given by,

(14) E(XN;Z)=KNKMλ(λ+μ), E(XM;Z)=KNδMλ(λ+μ)andE(XP;Z)=KPδPKNδMλ(λ+μ),

respectively. The supplementary material (Appendix B) provides more detail on how these expressions can be obtained.

We consider first the nascent-mature pathway reporters, case (a). From (14), it is easily seen that E(XN;Z)=f(ZN)KNλ(λ+μ) and E(XM;Z)=g(ZM)KNλ(λ+μ), where f(ZN)=1KM and g(𝐙M)=1δM. So the NDP holds, and the normalised covariance of E(XN;Z) and E(XM;Z) will identify the extrinsic noise on the transcription component KNλ(λ+μ). Consider now the mature-protein reporters, case (b). Again, we can see from (14) that E(XM;Z)=f(ZP)E(X;Z), where f(ZP)=KPδP. Thus, the NDP holds, and so the normalised covariance of E(XM;Z) and E(XP;Z) will identify the total noise on E(XM;Z) (the extrinsic noise on XM). For the nascent-protein reporters, case (c), it is easy to see that E(XN;Z)=f(ZN)KNλ(λ+μ), where f(ZN)=1KM, and E(XP;Z)=g(ZP)KNλ(λ+μ), where g(ZP)=KPδMδP. Thus, again the NDP holds, and the normalised covariance of E(XN;Z) and E(XP;Z) will identify the noise on the transcriptional component KNλ(λ+μ).

In order for the pathway-reporter method to provide a close approximation to the extrinsic noise in cases (a), (b), and (c), we require that the normalised intrinsic contribution to the covariance is either zero or negligible. This condition will hold provided there is sufficiently small correlation between the reporter pairs. In the case of (prokaryotic) mRNA and protein, this lack of correlation has been been verified experimentally in E. coli (Taniguchi et al., 2010). More generally, it is possible to provide an intuition about the conditions under which the lack of correlation might hold. The time series of copy numbers for each of nascent mRNA, mature mRNA and protein broadly follow each other, each with delay from its predecessor (Figure 4B). Parameter values that reduce this delay will tend to increase correlation, and thereby increase the normalised intrinsic contribution to the covariance. The primary example of this effect is seen when δp approaches, or even exceeds δm (or for nascent-mature reporters, when δM approaches the maturation rate). A further contributor to high correlation between mRNA and protein, is when the system undergoes long timescale changes. In this situation, the copy numbers tend to drop to very low values for extended periods. The primary parameter influencing this type of behaviour is the active-rate λ, specifically, when λ tends to 0 (Figure 4C). An illustrative example of this can be seen by considering a Telegraph system in the limit of slow switching, which produces a copy number distribution that converges to a scaled Poisson-Bernoulli compound distribution: even without any extrinsic noise, the pathway reporter method will identify the η2 value of the corresponding scaled Bernoulli distribution.

Figure 5 with 1 supplement see all
Heatmaps for the intrinsic contribution to the covariance.

These heatmaps estimate the level of overshoot in the pathway-reporter approach for the nascent-protein and mature-protein reporters; blue regions show an overshoot of less than 0.05.

Here, the intrinsic contribution is calculated using stochastic simulations of the model 𝐌4. For the mature-protein and nascent-protein reporters, we consider three different values of the parameter μ, specifically μ=2, μ=10, and μ=20. In all cases, the parameter δp and the on-rate λ are varied between 0.01 and 0.5, and 0.5 and 5, respectively. The parameters of the model 𝐌4 are scaled so that δM=1. The maturation rate is fixed at 20, with the parameters KN and KP chosen to produce a mean protein level of 1000, a mean nascent mRNA level of 5 and a mean mature mRNA level of 50. Each individual pixel is generated from a sample of size 3000, although there is still some instability in the convergence for the nascent-protein reporter, particularly as the overshoot estimation starts to increase, and particularly as μ is larger. To produce more accurate values, the case of μ=2 was averaged over two full experiments while μ=20 was averaged over three. This was also done for the mature-protein reporter, however for these images there was almost no visible difference between the various runs of the experiment and their averages. Each of the three μ values takes approximately 7–10 hr of computation, depending on lead in time before sampling within a simulation. Figure 5—figure supplement 1 gives a heatmap for the overshoot in the pathway-reporter approach for nascent-mature pathway reporters.

An extensive computational exploration of the parameter space (Supplementary file 2) supports our intuition, though the strength of the effect varies across the three different reporter pairs. This is further corroborated by the heatmaps shown in Figure 5, which for three fixed values of μ and a broad spectrum of δp and λ values, give the intrinsic term ηint2 in (8), for fixed 𝐙. Thus, the heatmaps provide an estimate for the overshoot error in the pathway-reporter approach. Note that blue pixels represent an overshoot estimate of less than 0.05.

For nascent-protein reporters, the normalised intrinsic contribution to the covariance is satisfactorily small (less than 0.05) except for (a) high values of δp in unison with (b) low values of λ (less than 1, although lower values are acceptable if δp is small). The assumption (a) that δP<δM is known to be true for a large number of genes, and is justified by the difference in the mRNA and protein lifetimes. While there is of course variation across genes and organisms, values of δP0.5δM and even δP0.2δM seem reasonable for the majority of genes. In E. coli Taniguchi et al., 2010 and yeast Belle et al., 2006, for example, mRNA are typically degraded within a few minutes, while most proteins have lifetimes at the level of 10 s of minutes to hours. For mammalian genes Schwanhäusser et al., 2011, it has been reported that the median mRNA decay rate δM is (approximately) five times larger than the median protein decay rate δP, determined from 4200 genes. Assumption (b) requires that the gene is sufficiently active. In a recent paper by Larsson et al., 2019, the promoter-switching rates λ, μ, and the transcription rate K of the Telegraph model are estimated from single-cell data for over 7000 genes in mouse and human fibroblasts. Of those genes with mean mRNA levels greater than 5, we found that over 90% have a value for λ of at least 0.5, and over 65% have a value for λ greater than 1. In Cao and Grima, 2020, a comprehensive list of genes (ranging from yeast cells through to human cells) with experimentally inferred parameter values are sourced from across the literature (see Table 1 in Cao and Grima, 2020). After scaling the parameter values of the 26 genes reported there, we find that around 88% have a value for λ of at least 0.5, and approximately 58% have a value for λ greater than 1. Thus, the heatmaps given in Figure 5 (top panel) suggest that nascent-protein reporters will provide a satisfactory estimate of the extrinsic noise level for a substantial fraction of genes.

The mature nRNA–protein reporters are less reliable, with the requirement of slow protein decay and higher on-rate being more pronounced than for the nascent mRNA–protein reporters; this is evident from Figure 5 (second panel). The performance of the nascent–mature reporter is of course independent of δp, but is only viable in the case of a large on-rate (see Figure 5—figure supplement 1).

We test our approach for each of the pathway reporter pairs (a), (b), and (c) against dual reporters using stochastic simulations. Table 4 shows the results from six simulations across a spectrum of behaviours from moderately slow switching, to fast switching as well as a range of levels of burstiness. For each of the parameters μ,λ,KP,δP, we selected a scaled Beta(5,6) distribution, with squared coefficient of variation η2=0.1; the scaling is chosen in each case to achieve a mean value equal to the parameter value given in Table 4. It is routinely verified that scaling these distributions does not change the value of η2. The parameter KN is given the noise distribution Beta(3,6), which has a slightly higher coefficient of variation η2=0.2. In order to achieve direct benchmarking against the dual reporters, the parameter KM is fixed; we select KM=10. This is because the nascent-protein pathway reporter estimates noise on the value of KNλλ+μ, while the mature mRNA dual-reporter measures noise on KNKMλλ+μ, and these coincide only when KM is fixed. The mean values of KN and KP are chosen to achieve approximate average nascent mRNA levels, mature mRNA levels and protein levels at 5, 50, and 1000 respectively, given the chosen values of λ,μ,δP.

Table 4
A comparison of the pathway-reporter method and dual-reporter method for bursty expression.

Here PR (NP) gives the results of the nascent and protein pathway reporters, PR (MP) gives the results of the mRNA and protein reporters, while DR (Mat) gives the results of the dual reporters calculated from the mature mRNA. We consider noise on all of the parameters except for δM and KM; see discussion in main text. The values given are the average of 100 simulations, each calculated from 500 copy number samples, and the errors are ± one standard deviation. Our theory predicts that pathway-reporters will identify the noise at both the promoter level (λ,μ) and transcriptional level (KN); the total extrinsic noise in each case is given by ηext2. As before, the noise distribution parameters are chosen to produce an average nascent mRNA copy number of 5 and an average mature mRNA copy number of 50, and an average number of 1000 proteins.

MeanSimulation
(r)1-5 λμKNKPδPPr (MP)Pr (NP)DR (Mat)
0.5115020.10.46±0.060.38±0.070.32±0.07
1215020.10.39±0.050.34±0.070.32±0.05
120105020.10.66±0.150.52±0.220.50±0.15
2210060.30.35±0.040.29±0.050.27±0.03
22055060.30.61±0.090.47±0.150.47±0.09
101010060.30.29±0.030.27±0.040.27±0.02
Table 4—source data 1

This is an Excel spreadsheet containing the data used to produce the final values in Table 4.

https://cdn.elifesciences.org/articles/69324/elife-69324-table4-data1-v2.xlsx

The results for the nascent mRNA–protein reporters, case (c), given in Table 4 show comparable performance to dual reporters, with only modest overshoot; even in the worst performing case of λ=0.5, μ=1 the result of the pathway reporters is within one standard deviation, in a very tight distribution. The error heatmaps of Figure 5 provide an accurate estimate of the overshoot in the nascent-protein results in Table 4. As an example, the first row is most closely matched by the heatmap at top left of Figure 5, which at λ=0.5 and δP=0.1 is suggestive of an error around the boundary between blue and red (around 0.06). The same accuracy is obtained for the other rows. As predicted, the mature-protein reporters show significantly more overshoot, especially with the less active genes. Improved accuracy can again be obtained by subtracting the estimated overshoot given in the error heatmaps from the obtained value. Thus for example, the error heatmap for μ=2 (Figure 5 lower left) gives an error approximately 0.07 for λ=1,δP=0.1, which agrees very closely to the actual overshoot of 0.07 shown in the corresponding row of Table 4. An overshoot of approximately 0.06 is suggested by the heatmap for μ=2, when λ=2,δP=0.3, which leads to a correction from 0.35 in Table 4 to a value of 0.29. This is quite consistent with the dual reporter benchmark of 0.27. As expected (based on Figure 5—figure supplement 1), nascent-mature reporters do not perform well on bursty systems except for high λ and so the values are not included in Table 4; only in the case of λ=μ=10 does the result begin to approach the dual reporter value, returning 0.32±0.03.

Generality of the pathway reporter method

To test the robustness of our pathway reporter approach, we validate our theoretical results on various other gene expression dynamics. (1): We begin by considering a more detailed model of the mRNA maturation process, where the nascent mRNA maturate after a fixed amount of time. The assumption of a fixed maturation time is sometimes taken to approximate the combined effect of intermediate maturation processes such as initiation, elongation, and release Xu et al., 2016. More specifically, we consider the model 𝐌4 above (Figure 4D) in the case of constitutive expression (λ=1, μ=0), and replacing the first-order reaction KM by a fixed interval of time. Additionally, we explore maturation times sampled from Erlang distributions, to account for the fact that maturation can involve several shorter stochastic processes. We find that the extrinsic noise contribution obtained using the nascent and mature mRNA reporters match closely to the true (dual reporter) values across a range of maturation times; refer to the supplementary material for details.

(2): Next we consider an extension of a model of transcriptional bursting given in Bartman et al., 2019; Cao et al., 2020. The model we consider is the same as in Bartman et al., 2019; Cao et al., 2020, however, is extended to include protein synthesis and degradation. This model captures the transcriptional process at a finer level of detail, and is argued in Cao et al., 2020 to be the model most closely matching experimental observations. In this more nuanced ‘multiscale’ model, the gene stochastically switches between three states: two active states S10 and S11, and one inactive state S0. The activation of the gene occurs in two steps, initially by the binding of transcriptional factors (transition from S0 to S10 at rate λ1, and reversible at rate μ1), and then as a secondary step, by the binding and pause of the mRNA polymerase (transition from S10 to S11 at rate λ2). Transitions from S11 to S0 also occur at rate μ1, due to detachment of both the transcriptional factors and polymerase. Transcription of nascent mRNA (at rate KN) occurs only in state S11 and results in immediate transition to state S10. Nascent mRNA maturate at rate KM, and are subsequently translated into protein at rate Kp. Degradation of mRNA and protein occur with rates δm and δp, respectively. All reactions are considered to be first-order with exponentially distributed waiting times between successive reactions. A schematic of the model is given in Figure 6 (inner rectangle). In this case, we are again able to prove that the Noise Decomposition Principle holds for all reporter pairs taken from this pathway using existing formulæ derived in Cao et al., 2020 for the marginal means. For details refer to the supplementary material (Appendix Convergence of Pathway and Dual Reporters).

Multiscale model of transcriptional bursting with additional features of the cell cycle.

In this model, the gene stochastically switches between three states: two active states, S10 and S11, and one inactive state S0. Gene activation occurs in two steps, initially by the binding of transcription factors (at rate λ1, reversible at rate μ1), and then as a secondary step by the binding and pause of the mRNA polymerase (at rate λ2). Transitions from S11 to S0 also occur at rate μ1, due to detachment of both the transcriptional factors and polymerase. Transcription of nascent mRNA (at rate KN) occurs only in state S11 and results in immediate transition to state S10. Nascent mRNA mature at rate KM, and are subsequently translated into protein at rate Kp. Degradation of mRNA and protein occur with rates δm and δp, respectively. We verify our pathway reporter method on three variations of the multiscale model. First, we assume all reactions are first-order Poisson processes (Case (2) in the main text). We then incorporate further details of the mRNA maturation process, where nascent mRNA occurs after a fixed amount of time (Case (3)). Finally, we incorporate features of the cell-cycle such as gene replication, dosage compensation, cell division, and cell-cycle length variability, as well as incorporating more realistic Erlang distributed maturation times (Case (4)).

(3): We combine models (1) and (2) above, incorporating the fixed time maturation of (1) with the multiscale approach of (2).

(4): The cell cycle is a major contributing factor to extrinsic noise, introducing population heterogeneity (as cells are at different stages of the cell cycle), as well as internal dynamics to parameter values. Here we incorporate the salient features of the cell cycle into model (3), which is measurable as extrinsic noise by our methodology. Specifically, we model the effects of (i) gene replication, (ii) dosage compensation, (iii) binomial partitioning of products due to cell division, as well as (iv) cell-cycle length variability. Refer to Figure 6 for a schematic. This model is an extension of that solved in Cao and Grima, 2020. Our model further incorporates the multi-scale dynamics of model (3) and the Erlang-distributed maturation times of model (1). As far as we are aware, this model has not been explored even by stochastic simulations before. A detailed description of how the above cell-cycle effects are captured in our model is given as follows. (i) Replication results in a doubling of the gene from one to two at the replication time, tr. This replication occurs over a period which is much shorter than the length of the cell cycle, and we follow the assumptions in Cao and Grima, 2020 by considering it to occur instantaneously. (ii) Dosage compensation changes the rate at which the gene switches from inactive to active (λ1) upon replication at time tr. Again following Cao and Grima, 2020, the activation rate after replication is 70% of the rate prior to replication. (iii) Binomial partitioning of nascent mRNA, mature mRNA and protein at cell division. We assume that nascent mRNA, mature mRNA, and protein segregate into the two daughter cells, with independent probability 1/2. We follow just one of the daughter cells, chosen at random with equal probability. (iv) Cell-cycle length variability. The time between successive cell divisions is stochastic, and is assumed to be sampled from an Erlang distribution. This has been shown to be consistent with experimental data Cao and Grima, 2020. The time to replication, and subsequently to cell division, are both chosen from an Erlang distribution with shape parameter 12, which produces a total cell cycle length distributed according to Erlang(24,λ); this matches the Erlang(24,λ) cell cycle length in Cao and Grima, 2020, where replication is at the exact halfway point. We similarly choose a rate parameter λ=λ1 chosen at a commensurate proportion to our mRNA decay rate of δ=1.

In each of the cases (1 – 4) above, we find that the correlations between reporter pairs is negligible, and the predicted contribution of extrinsic noise matches those obtained from the dual reporter method across a range of parameter combinations. Details of the simulation methods and results can be found in the supplementary material (Appendix Convergence of Pathway and Dual Reporters). In summary, the results show that our pathway reporter approach is remarkably insensitive to the specific dynamics of mRNA and protein synthesis. In particular, the correlations between reporter pairs do not strongly depend on the details of the gene expression model used.

Discussion

Despite the proliferation of experimental methods for single-cell profiling, the ability to extract transcriptional dynamics from measured distributions of mRNA copy numbers is limited. In particular, the multiple factors that contribute to mRNA heterogeneity can confound the measured distribution, which hinders analysis. Theoretical innovations that allow us to quantify and help in identifying the causes of these observable effects are therefore of great importance. In this work, we have demonstrated, through a series of mathematical results, that it is impossible to delineate the relative sources of heterogeneity from the measured transcript abundance distribution alone: multiple possible dynamics can give rise to the same distribution. Our approach involves establishing integral representations for distributions that are commonly encountered in single-cell data analysis, such as the negative binomial distribution and the stationary probability distribution of the Telegraph model. We show that a number of well-known representations can be obtained from our results. A particular feature of our non-identifiability results is that population heterogeneity inflates the apparent burstiness of the system. It is therefore necessary to collect further information, beyond measurements of the transcripts alone, in order to constrain the number of possible theoretical models of gene activity that could represent the system. In particular, additional work may be required to determine the true level of burstiness of the underlying system.

We have developed a theoretical framework for estimating levels of extrinsic noise, which can assist in resolving non-identifiability problems. The dual reporter method of Swain et al., 2002 already provides one such approach; but it is experimentally challenging to set up in many systems, and requires strictly identical and conditionally independent pairs of gene reporters. Our Noise Decomposition Principle directly generalises the theoretical underpinnings of the dual reporter method and related approaches Bowsher and Swain, 2012; Jetka et al., 2018; here we have used it to identify a practical approach—the pathway-reporter method—for obtaining an effective and experimentally more easily implementable noise decomposition. Our approach allows us to use measurements of two different species from the transcriptional pathway of a single gene copy instead of having to set up a more cumbersome dual reporter assay. The accuracy of the pathway-reporter method is provably identical for constitutive gene expression, and in the case of nascent-mature mRNA reporters, the measurements are readily obtainable from current single-cell data Shah et al., 2018; La Manno et al., 2018; Skinner et al., 2016. For bursty systems, the method in general provides only an approximation of the extrinsic noise. We are, however, able to demonstrate computationally, that one of the proposed pathway reporters provides a satisfactory estimate of the extrinsic noise for most genes. The other pathway reporters also provide viable estimates of the extrinsic noise in some cases. We further validate our theoretical framework on synthetic data for genes with various underlying gene expression dynamics. Our results show that the pathway reporter method is independent of the specific dynamics of mRNA and protein synthesis, and therefore should be applicable to a broad range of experimental scenarios.

Despite the generality of our theoretical contribution, our pathway-reporter approach has some caveats. In particular, the approach relies on the assumption that extrinsic noise sources act independently. Experimentally, however, these may be correlated. For example, it has been suggested Hilfinger et al., 2016; Cole and Luthey-Schulten, 2017 that the transcription and translation rates in E. coli anticorrelate. Additional work is required to determine degree to which the independence of noise sources is a reasonable assumption.

Recent developments in single-cell profiling now allow simultaneous measurements of transcripts and proteins in thousands of single cells Stoeckius et al., 2017; Peterson et al., 2017; Reimegård et al., 2019. As discussed in Gorin and Pachter, 2020a, experimental improvements that would additionally allow measurements of nascent transcripts, coupled with theoretical developments such as those presented here, will allow for identification of noise sources on a genome-wide scale. Our work reveals that extrinsic noise distorts the multivariate copy number distribution of the different species in the gene expression pathway. We have exploited this to yield reliable estimates of noise strength, which we are confident will assist in setting better practices for model fitting and inference in the analysis of single-cell data. A more nuanced analysis of this multivariate distribution may offer even further insight into model and noise identifiability.

Appendix 1

Derivations of the non-identifiability results

In the main text we provide a number of examples of when the compound distribution does not have a unique representation. We here provide the full details of the derivations of these results.

Bursty expression: telegraph representation

Consider first a Telegraph distribution p~T(n;λ,μ,K) and the probability density function of the scaled beta distribution BetaK(λ+μ,μμ), given by

fK(t;λ+μ,μ-μ)=Γ(λ+μ)Γ(λ+μ)Γ(μ-μ)K1-λ-μtλ+μ-1(K-t)μ-μ-1.

Note that this distribution has support [0,K]. In the main text, we claim that the Telegraph distribution p~T(n;λ,μ,K) can be obtained from compounding a Telegraph distribution by a scaled beta distribution with pdf f(t;λ+μ,μ-μ). In other words, that p~T(n;λ,μ,K) can be written as:

(15) p~T(n;λ,μ,K)=0Kp~T(n;λ,μ,t)fK(t;λ+μ,μμ)dt.

This is Equation (3) in the main text. Starting from the right hand side of (15), we have

(16) 0Kp~T(n;λ,μ,t)fK(t;λ+μ,μμ)dt=1n!Γ(λ+μ)Γ(λ+μ)Γ(μμ)Γ(λ+n)Γ(λ)Γ(λ+μ)Γ(λ+μ+n)0K1F1(λ+n,λ+μ+n,t)(K)1λμtλ+μ+n1(Kt)μμ1dt.

Substituting t=KT and simplifying, the right hand side of (16) becomes

1n!Γ(λ+μ)Γ(μμ)Γ(λ+n)Γ(λ)(K)nΓ(λ+μ+n)011F1(λ+n,λ+μ+n,KT)Tλ+μ+n1(1T)μμ1dT.

Now, using the integral representation Olver et al., 2010 [13.4.2] of the confluent hypergeometric function (with a=λ+n, b=λ+μ+1 and c=λ+μ+1), this becomes

1n!Γ(λ+μ)Γ(μμ)Γ(λ+n)Γ(λ)(K)nΓ(λ+μ+n)Γ(λ+μ+n)Γ(μμ)Γ(μ+λ+n)1F1(λ+n,λ+μ+n,K)=Γ(λ+μ)Γ(λ+μ+n)Γ(λ+n)Γ(λ)(K)nn!1F1(λ+n,λ+μ+n,K)=p~T(n;λ,μ,K),

which is the left hand side of (15). Hence, we have that (15) holds.

Instantaneously bursty expression: negative binomial representations

We consider first the representation given in Equation (4) in the main text. Recall that we let p~NB(n;r,β) denote the probability mass function of a NegBin(r,ββ+1) distribution, where β>0. In the main text, we claim that for any negative binomial distribution of the form NegBin(λ,ββ+1) (where β>0) we have, 

(17) p~NB(n;λ,β)=0p~T(n;λ,μ,t)f(t;λ+μ,β)𝑑t,

where f(t;λ+μ,β) is the probability mass function of a Gamma(λ+μ,β) distribution. Beginning with the right hand side of (17) we have,

(18) 0p~T(n;λ,μ,t)f(t;λ+μ,β)dt=1n!Γ(λ+n)Γ(λ)Γ(λ+μ)Γ(λ+μ+n)β(λ+μ)Γ(λ+μ)01F1(λ+n,λ+μ+n,t)tλ+μ+n1eβtdt.

We now apply the following identity given in Saad and Hall, 2003 [Appendix 1] with a=λ+n, b=λ+μ+n, k=-1, d=λ+μ+n and h=β.

(19) 01F1(a,b,kt)td1eht=Γ(d)hd(1kh)a.

So the left hand side of (18) becomes

1n!Γ(λ+n)Γ(λ)Γ(λ+μ)Γ(λ+μ+n)β(λ+μ)Γ(λ+μ)Γ(λ+μ+n)βλ+μ+n(1+1β)(λ+n)=1n!Γ(λ+n)Γ(λ)βλ(β+1)λ+n=p~NB(n;λ,β),

which is the right hand side of (17). Hence, we have that (17) holds. We now consider the representation given in Equation (5) of the main text. Here we claim that any negative binomial distribution of the form NegBin(λ,b1+b) (where b>0) can be written as,

(20) p~NB(n;λ,b)=1bp~NB(n;λ,θ)fb(bθ1;λλ,λ)dθ,

 where fb(bθ-1;λ-λ,λ) is the probability mass function of a scaled beta prime BetaPrimeb(λλ,λ) distribution, where b>0 and λ>λ. Starting from the right hand side of (20), we have

1bp~NB(n;λ,θ)fb(bθ1;λλ,λ)dθ=1bΓ(λ+n)Γ(n+1)Γ(λ)(1θθ+1)n(θθ+1)λbΓ(λ)Γ(λθ)Γ(θ)(bθ1)λλ1(bθ)λ=b1λΓ(λ+n)Γ(n+1)Γ(λ)Γ(λλ)1b(bθ1)λλ1(1+θ)λ+n.

Now substituting H=bθ-1 and simplifying, we obtain

(21) bnΓ(λ+n)Γ(n+1)Γ(λ)Γ(λλ)0Hλλ1(H+b+1)λ+ndH,

and letting y=b+1 and simplifying we have that

(22) 0Hλλ1(H+b+1)λ+ndH=1(b+1)λ+n0yλλ1(y+1)λ+ndy
(23) =1(b+1)λ+nΓ(λλ)Γ(λ+n)Γ(λ+n).

Here, we used the fact that the integral in the variable y is the probability density function of a BetaPrime(λλ,λ+n) distribution. Thus, Equation (21) simplifies to

(24) Γ(λ+n)Γ(n+1)Γ(λ)bn(b+1)λ+n=p~NB(λ,b),

 which is the right hand side of (20). Thus, we have that (20) holds.

Appendix 2

Proof of the Noise Decomposition Principle

Here we provide the proof of the Noise Decomposition Principle given in the main text. For convenience we restate the principle below.

The Noise Decomposition Principle (NDP)

Assume that there are measurable functions f, g and h such that E(X;Z) and E(Y;Z) split across common variables by way of E(X;Z)=f(ZX)h(Z) and E(Y;Z)=g(ZY)h(Z). Then provided that the variables Z1,,Zm are mutually independent, the normalised covariance of E(X;Z) and E(Y;Z) will identify the total noise on h(Z) (i.e., ηh(Z)2).

Consider first the covariance of E(X;Z) and E(Y;𝐙). We have

(25) Cov(E(X;Z),E(Y;Z))=Cov(f(ZX)h(Z),g(ZY)h(Z))=E(f(ZX)g(ZY)(h(Z))2)E(f(ZX)h(Z))E(g(ZY)h(Z))=E(f(ZX))E(g(ZY))[E((h(Z))2)(E(h(Z)))2].

Note that here we used the fact that the variables in Z=(Z1,,Zn) are mutually independent. Now using the fact that E(X)=E(E(X;Z)) and E(Y)=E(E(Y;Z)) (the Law of Total Expectation), and then normalising we obtain

(26) Cov(E(X;Z),E(Y;Z))E(X)E(Y)=E(f(ZX))E(g(ZY))[E((h(Z))2)(E(h(Z)))2]E(f(ZX))E(g(ZY))(E((h(Z)))2=Var(h(Z))(E(h(Z)))2=ηh(Z)2.

Hence, under certain conditions the normalised covariance of E(X;Z) and E(Y;Z) will identify the total noise on h(𝐙).

Appendix 3

Pathway reporters

Justification for simulating only one copy of the gene

Our simulations and theory have been based over reporters from a single gene copy, whereas in practice there may be multiple copies of the gene that contribute to the overall mRNA. If there are mechanisms in place to distinguish the mRNA or protein of one gene copy from another then the theory and analysis we have developed in main paper holds without change. In the case where it is not possible to distinguish mRNA and protein (respectively) from the multiple gene copies, then we now observe that the general theory continues to hold, provided there is independence between the gene copies; this assumption has been verified experimentally Skinner et al., 2016. Here, we demonstrate this in the case of two gene copies, though the general case for more than two genes is essentially identical but is notationally cumbersome. We are considering a situation where the variable X in the Noise Decomposition Principle is the sum of two independent variables X1,X2 and Y is the sum of two independent variables Y1,Y2. We assume common dependence on the environmental variables 𝐙 so that E(X1;𝐙)=E(X2;𝐙), E(Y1;𝐙)=E(Y2;𝐙). Using these equalities and the independence of X1,X2 and Y1,Y2 in X=X1+X2, Y=Y1+Y2, we find the numerator of

Cov(E(X;Z),E(Y;Z))E(X)E(Y)

is simply 4Cov(E(X1;Z),E(Y1;Z)), while the denominator is 4E(X1)E(Y1). Thus the noise decomposition coincides with that for the single copy X1,Y1. Further work may be required to consider systems where there is independence between, or there is significant deviations in the gene copies.

Upper bound for the intrinsic contribution to the covariance: constitutive expression

In the main text, we claim that the error in the pathway-reporter approach in the case of mRNA-protein reporters is negligible (i.e. the error is 1); refer to Equation (13) in the main text. We here provide full details of the derivation of this expression. First let Xm and Xp be the number of mRNA and protein produced from the same constitutive gene modelled by the 'two-stage' model, 𝐌2 (see Figure 4A (top right) of the main text). Also, let Z={Km,δm,Kp,δp}.

We restate here the expression for the intrinsic contribution to the covariance of Xm and Xp given as Equation (13) of the main text.

(27) ECov(Xm,Xp;Z))E(Xm)E(Xp)=αE(Km),where α=E(1/(δp+1))E(1/δp).

We require the following expressions for the stationary mean mRNA level and protein level of the two-stage model Thattai and van Oudenaarden, 2001; Singh and Hespanha, 2007.

(28) E(Xm;Z)=Kmδm,E(Xp;Z)=KpδpKmδmandCov(Xm,Xp;Z)=KmKpδm(δm+δp).

Assuming that δm is fixed across the cell-population, and all parameters are scaled so that δm=1, it follows from (28), that

(29) E(Cov(Xm,Xp;Z))=E(1δp+1)E(Km)E(Kp).

Using the Total Law of Expectation, we also have,

(30) E(Xm)=E(Km)andE(Xp)=E(1δp)E(Km)E(Kp).

Thus, it follows that

(31) E(Cov(Xm,Xp;Z))E(Xm)E(Xp)=αE(Km),where α=E(1/(δp+1))E(1/δp).

Determination of the marginal means for model 𝐌4

In order to establish the Noise Decomposition principle in the case of bursty gene expression, we rely on expressions for the marginal stationary means of the model 𝐌4; see Figure 4A (bottom left) and the associated caption for more details. To derive the marginal means for nascent and mature mRNA and for protein (Equation (14) of the main text), we first observe that the nascent mRNA population may treated identically to that of mRNA in general (that is, no distinction between nascent and mature), as in Peccoud and Ycart, 1995, except that mRNA decay is replaced by the sum of decay and maturation. As in the work of Cao and Grima, 2020, the assumption of fast maturation allows us to ignore decay completely in the nascent phase, so that the marginal distribution is identical to that of Peccoud and Ycart, 1995, except with decay replaced by maturation. This leads to a marginal nascent mRNA mean of

E(XN;Z)=KNδMλ(λ+μ)

The marginal means for mature mRNA and protein are derived in Raj et al., 2006 under the assumption that the transcription rate parameter KN is large relative to the other parameters. The expressions are given by

E(XM;Z)=KNδMλ(λ+μ) and E(XP;Z)=KPδPKNδMλ(λ+μ).

Formally, the marginal means in Raj et al., 2006 are for the three-stage model 𝐌3, which ignores the downstream processing of mRNA, such as splicing. The assumption of fast maturation however, justifies the treatment of the nascent phase of mRNA as an ephemeral step within the Poissonian modelling of mRNA transcription.

There are a number of possibly compounding assumptions on the parameters here, but simulations show that there is a lot of tolerance, with even only moderate maturation and transcription still returning sample means consistent with the formulas.

Measuring noise from an instantaneously bursty gene

Here we use the simple stochastic model of Singh and Bokes, 2012 that includes both instantaneous transcriptional bursting and mRNA maturation. This model can be obtained as a limiting case of the model 𝐌4 (Figure 4A (top left) in the main text), where the off-rate μ has tended toward infinity, while the on-rate λ remains fixed. Under this condition, transcription is rare enough to be considered instantaneous, leading to ‘bursts’ of transcriptional activity. The intensity of the bursts M is known to be geometrically distributed with mean burst size B=KNμ.

If we let XN and XM denote the number of nascent and mature mRNA, respectively, then according to Singh and Bokes, 2012, the steady-state marginal means and covariance are given by

(32) E(XN;Z)=λBKM,E(XM;Z)=λBδMandCov(XN,XM)=λB2KM+δM.

It is easily seen from (32) that E(XN;Z)=f(ZN)Bλ and E(XM;Z)=g(ZM)Bλ, where f(ZN)=1KM and g(ZM)=1δM. Thus, the Noise Decomposition Principle holds. Using arguments similar to those given in the above section, it is straightforward to derive the following expression for the error in the pathway reporter approach:

(33) E(Cov(XN,XM;Z))E(XN)E(XM)=αE(λ)(ηB2+1),where α=E(1/(KM+1))E(1/KM).

We can expect that α1, so that unless the burst frequency is substantial, that is λ1, the overshoot in the the nascent-mature pathway-reporter approach is significant.

Appendix 4

Convergence of pathway and dual reporters

Appendix 4—figure 1
Comparison of convergence of η2 estimates for low and high mRNA levels by way of nascent-mature reporters and mature-mature reporters.

Low level corresponds to mean nascent mRNA level of 0.5, and mean mature mRNA level of 5. High level corresponds to mean nascent mRNA level of 5, and mean mature mRNA level of 50. In both cases, the simulated genes are constitutive and noise is on all parameters except for δM=1. The green line gives the squared coefficient of variation for KN, set to 0.2, which is the value the various reporters are expected to estimate. (A) Convergence of the η2 estimate over the first 2000 samples in the low- and high-output genes. (B) Convergence of the η2 estimate over 100,000 samples in the low-output gene only.

Convergence of both pathway reporters and dual reporters in low output genes is slower than for high output genes, and this affect is most pronounced for reporters taken from nascent mRNA. Fig. Convergence of Pathway and Dual Reporters compares convergence in some low output genes and high output genes for nascent-mature pathway reporters and the mature-mature dual reporter, in the case of constitutive expression. From Equation (11) of the main text and the NDP, these should measure precisely the extrinsic noise on the transcription rate parameter KN. In the low output gene, both the dual reporter and pathway reporters are yet to show accurate measurement of the overall extrinsic noise after 1000 samples (Figure Convergence of Pathway and Dual ReportersA), although they are providing rough estimates after as little as a few hundred samples. Pathway reporters for nascent-mature typically exhibit slightly slower convergence than mature-mature dual reporter; the difference is marginal, but can be seen for both the high output gene (compare the values after 500 samples) and the low output gene (compare the values after 1000 samples). In this figure, all parameters (except the reference parameter δM) were given scaled Beta distribution noise. The noise on KN is a Beta(3,6) distribution scaled to achieve E(KN)=5 in the low output case and E(KN)=50 in the high output case. The squared coefficient of variation is 0.2, and is shown in green in the figure.

Appendix 4—figure 2
Comparison of convergence for low and high mRNA levels by way of mature-protein and mature-mature reporters.

Low level corresponds to mean nascent mRNA level of 0.5, and mean mature mRNA level of 5. High level corresponds to mean nascent mRNA level of 5, and mean mature mRNA level of 50. In both cases the simulated genes are constitutive and noise is on all parameters except for δM=1. The noise on KN has squared coefficient of variation equal to 0.2, which is shown as the red horizontal line. Our theory shows that mature-protein reporters will return an overshoot that is negligible in the high-output gene (the blue horizontal line), but larger in the low output gene (light blue horizontal line); these values are calculated in the text. (A) Comparison of convergence for low and high mRNA levels over the first 2000 samples. (B) Convergence of the η2 estimate over 20,000 samples in the case of the low-output gene only. Two examples of each are given, to show the variation in behaviour.

Convergence for the mature-protein reporters is considered in Figure Convergence of Pathway and Dual Reporters. Here, Equation (13) of the main text shows that we should expect an overshoot in comparison to the mature-mature dual reporter. We have again compared a low output gene with a high output gene, and with all parameters (except δM) experiencing noise. The transcription is again given scaled Beta(3,6) distribution, having squared coefficient of variation 0.2, with scaling to achieve E(KN)=5 in low output and E(KN)=50 in the high output. The noise distribution for δP is a Beta(8,6) distribution, scaled to achieve a mean value of 0.2. Computational sampling from 107 samples finds the value of E(1/(δp+1))E(1/δp) as approximately 0.1573. The high output gene then expects an overshoot of 0.003146 from the mature-protein reporters, while the low output gene expects an overshoot of approximately 0.03146. Figure Convergence of Pathway and Dual ReportersA shows comparison of convergence over the first 2000 samples, with the theoretical values accommodating the calculated overshoot. Figure Convergence of Pathway and Dual ReportersB shows the same over 20,000 samples; this time two low-output two high-output genes are considered, so that the variation around the expected long term value can be seen. It is evident that reasonable estimates are given after a relatively modest number of samples, but there is a very long delay to a highly accurate convergence for the pathway reporters in the low output gene.

Appendix 4—figure 3
Convergence for reporter pairs for low gene activity.

In each case, the mean nascent mRNA level is 0.5, the mean mature mRNA level of 5, and the mean protein level is 500. The simulated genes are constitutive and noise is on all parameters except for δM=1. Each graph shows the convergence of 600 individual reporter simulations, for each combination of reporters from nascent mRNA, mature mRNA and protein. Each reporter simulation is from 10,000 samples, with the reporter estimates calculated at intervals of 100. The noise on KN has squared coefficient of variation equal to 0.2, which should be identified by both the nascent-mature reporter and mature-mature dual reporter. As in Figure Convergence of Pathway and Dual Reporters, the mature-protein reporter should converge to an estimate of approximately 0.2315. Nascent-nascent and protein-protein reporters identify combined noise on more parameters, so do not converge to 0.2. The lower graph shows each of nascent-mature, mature-mature, mature-protein and nascent-protein in the same plot for direct comparison.

Appendix 5

Generality of the pathway reporter method: model and simulation details

We verify our pathway reporter approach on synthetic data for genes with various underlying gene expression dynamics. In the main text, we consider four different scenarios, beyond the standard four-stage model of gene transcription. Below we give details of these models, and their method of simulation.

(1): First, we consider a more detailed model of the mRNA maturation process, where the nascent mRNA maturate after a more finely controlled time interval. Specifically we adjust the constitutive model of gene transcription (𝐌4 of the main text with λ=1 and μ=0) to include a fixed-duration of maturation. We also explore time-intervals sampled from Erlang distributions, to account for the fact that maturation can involve several shorter stochastic processes. The process is simulated using an adaptation of SSA, whereby the time of maturation of an individual nascent mRNA is calculated (by sample from the chosen distribution) at the point it is transcribed, and this timepoint is queued in a register of maturation times; this can be trivially adapted to other distributions; generalised-Erlang for example. At each iteration of the simulation, determination of the next event involves taking the smallest of the various exponentially sampled Markov processes (on, off, mRNA transcription, etc) and the next maturation event (the minimum of the maturation register). In the case that maturation is the next event, the corresponding maturation time is removed from the maturation register.

(2): We next consider a model of transcriptional bursting given in Bartman et al., 2019; Cao et al., 2020. In this ‘multiscale’ model, gene activation occurs in two steps, initially by the binding of transcription factors, and then as a secondary step by the binding and pause of the mRNA polymerase. For details of the model refer to Figure 6 and the surrounding paragraphs in the main text. We here show that our Noise Decomposition Principle holds for all reporters taken from this gene pathway, and verify computationally that correlations between some of the reporter pairs are negligible; the process is simulated by the SSA. In the following, we assume all reactions are first-order, characterised by exponential waiting times between successive reactions. Let XN denote the number of nascent mRNA, let XM denote the number of mature mRNA, and let XP denote the number of proteins produced from the same gene. Also let 𝐙={λ1,μ1,λ2,KN,KM,KP,δM,δP}. We assume that the maturation rate KM is large (i.e. KM>δM), which is supported by experiments Cao and Grima, 2020. Then using the results of Cao et al., 2020, the stationary averages for the nascent mRNA, mature mRNA and protein levels are given by,

(34) E(XN;Z)=KNKMλ1λ2γ1γ2, E(XM;Z)=KNδMλ1λ2γ1γ2 and E(XP;Z)=KPδPKNδMλ1λ2γ1γ2,

 respectively. Here γ1=λ1+μ1 and γ2=KN+λ2+μ1. We begin by considering the nascent-mature pathway reporters. From (34), it is easily seen that E(XN;Z)=f(ZN)KNλ1λ2γ1γ2 and E(XM;Z)=g(ZM)KNλ1λ2γ1γ2, where f(ZN)=1KM and g(𝐙M)=1δM. So the NDP holds, and the normalised covariance of E(XN;Z) and E(XM;Z) will identify the extrinsic noise on the transcriptional component KNλ1λ2γ1γ2. For the mature-protein reporters, we can again see from (34) that E(XM;Z)=f(ZP)E(X;Z), where f(ZP)=KPδP. Thus the NDP holds, and so the normalised covariance of E(XM;Z) and E(XP;Z) will identify the strength of extrinsic noise in the mRNA copy number distribution. For the nascent-protein reporters, it is easy to see that E(XN;Z)=f(ZN)KNλ1λ2γ1γ2, where f(𝐙N)=1KM, and E(XP;Z)=g(ZP)KNλ1λ2γ1γ2, where g(𝐙P)=KPδMδP. Thus, again the NDP holds, and the normalised covariance of E(XN;Z) and E(XP;Z) will identify the noise on the transcriptional component KNλ1λ2γ1γ2.

(3): We combine models (1) and (2) above, incorporating the fixed time maturation of (1) with the fine-level transcriptional approach of (2). Again, we employ the adaptation of the SSA described in (1).

(4): We incorporate the salient features of the cell-cycle into model (3), as well as introducing Erlang-distributed maturation times. Specifically, the model accounts for the effects of gene replication, dosage compensation, binomial partitioning of products due to cell division, as well as cell-cycle length variability. Within each phase of the cell cycle, the handling of maturation times for a given gene copy is exactly as in Case (1) above, but within the overall system described in Case (3). To capture the cell cycle, the SSA is now further modified to handle the two key sporadic changes: gene replication and cell division. After gene replication, we perform two independent modified SSA simulations, with common parameter values. The final copy numbers are obtained as a sum of those from each gene copy. At cell division, we follow just one daughter cell, selecting the inherited copy numbers by way of binomial partitioning from the mother cell values at the point of division (including the maturation register). The length of these cell-cycle phases is sampled from the chosen Erlang distribution, on commencement of the simulation of the phase.

The results of pathway reporters for each of the Cases (1 – 4) above are given in Appendix 5—tables 18 below. In all cases, the values given are the average of 20 simulations, each calculated from 500 copy number samples, and the errors are ± one standard deviation. For each model, we consider two parameter sets, each without extrinsic noise (Subcase A) and with extrinsic noise (Subcase B). For Case (1), the model parameters are chosen to produce an average nascent mRNA copy number of 10, and an average number of 2000 proteins; the average mature mRNA copy number varies according to the maturation time. We find that all of the pathway reporters accurately predict the true extrinsic noise levels (as given by dual reporters) across a range of maturation times; refer to Appendix 5—table 1, 2. For Cases (2 – 4), the model parameters are chosen to produce an average nascent mRNA copy number of 5 and an average mature mRNA copy number of 100, and an average number of 2000 proteins. We verify computationally that the mature-protein and nascent-protein reporter pairs provide accurate predictions of extrinsic noise contributions across a range of noise levels. For these results refer to Appendix 5—tables 38. We mention that the results of Case (4)B (Appendix 5—tables 8) may suggest a slight undershoot in the pathway reporter values in comparison to dual reporters. The difference is, however, within one standard deviation, and are calculated from a relatively small sample size due to the complexity of the simulation and corresponding run time. Investigating the presence and possible causes of this marginal effect for this model and other complicated models may be of interest in further work.

Appendix 5—table 1
A comparison of the pathway-reporter method and dual-reporter method for constitutive gene expression with Erlang-distributed maturation times (Case (1)A).

Here, PR (NP) gives the results of the nascent and protein pathway reporters, PR (MP) gives the results of the mRNA and protein reporters, while DR (Mat) gives the results of the dual reporters calculated from the mature mRNA. The maturation time Tmat is chosen to be Erlang distributed with mean length 1/30,0.05, and 0.1, respectively. We consider the rate parameters for the remaining exponentially distributed times to be constant, so that there is no extrinsic noise. The pathway-reporters correctly identify the zero extrinsic noise contribution.

ParametersSimulation
(r)1-6 KNTmat(mean)KPδPPr(NM)Pr (MP)Pr (NP)DR (Mat)
3000.03˙0.06˙0.10.00±0.0010.00±0.00010.00±0.00030.00±0.0001
2000.0510.10.00±0.0010.00±0.00010.00±0.00040.00±0.0001
1000.120.10.00±0.0010.00±0.00020.00±0.00010.00±0.0004
Appendix 5—table 2
A comparison of the pathway-reporter method and dual-reporter method for constitutive gene expression and fixed maturation time (Case (1)B).

For each of the parameters KP,δP we selected a scaled Beta(5,6) distribution, with squared coefficient of variation η2=0.1; the scaling is chosen in each case to achieve a mean value equal to the parameter value. The parameter KN is given the noise distribution Beta(3,6), which has a slightly higher coefficient of variation η2=0.2. In order to benchmark against dual reporters, the maturation time was fixed in each case. The extrinsic noise contribution predicted by the pathway-reporters matches well with the dual reporter values.

MeanSimulation
(r)1-4 TmatKNKPδPPr(NM)Pr (MP)Pr (NP)DR (Mat)
0.0520010.10.20±0.010.20±0.020.20±0.030.20±0.01
0.110020.10.20±0.010.21±0.030.21±0.030.20±0.01
Appendix 5—table 3
A comparison of the pathway-reporter method and dual-reporter method for the multiscale model (Case (2)A above).

We consider fixed parameters values (that is, no extrinsic noise). As our theory predicts, the pathway-reporters correctly identify zero extrinsic noise.

ParametersSimulation
(r)1-6 λ1μ1λ2=KNKMKPδPPr (MP)Pr (NP)DR (Mat)
224002020.10.02±0.0030.00±0.010.00±0.002
42012102020.10.02±0.0030.01±0.010.00±0.01
Appendix 5—table 4
A comparison of the pathway-reporter method and dual-reporter method for the multiscale model (Case 2.B).

For each of the parameters λ1,μ1,KP,δP we selected a scaled Beta(5,6) distribution, with squared coefficient of variation η2=0.1; the scaling is chosen in each case to achieve a mean value equal to the parameter value. The parameter λ2=KN is given the noise distribution Beta(3,6), which has a slightly higher coefficient of variation η2=0.2. In order to benchmark against dual reporters, the maturation rate was fixed at 20. As our theory suggests, the extrinsic noise contribution predicted by the pathway reporters matches well with the dual-reporter values.

MeanSimulation
(r)1-5 λ1μ1λ2=KNKPδPPr (MP)Pr (NP)DR (Mat)
2240020.10.18±0.030.16±0.040.16±0.02
420121020.10.29±0.040.27±0.070.28±0.03
Appendix 5—table 5
A comparison of the pathway-reporter method and dual-reporter method for the multiscale model with a fixed duration of maturation (Case (3)A).

Here the time to maturation Tmat is chosen to be consistent with the mean of the stochastic maturation time used in our other models (where the maturation time is exponentially distributed). We consider all rate parameters to be constant, that is, there is no extrinsic noise. Pathway-reporters correctly identify the zero extrinsic noise contribution.

ParametersSimulation
(r)1-6 λ1μ1λ2=KNTmatKPδPPr (MP)Pr (NP)DR (Mat)
224000.0520.10.02±0.0040.00±0.010.00±0.01
42012100.0520.10.02±0.0030.00±0.010.00±0.01
Appendix 5—table 6
A comparison of the pathway-reporter method and dual-reporter method for the multiscale model with a fixed duration of maturation (Case (3)B).

Here the maturation time, Tmat, is set to 0.05. For each of the parameters λ1,μ1,KP,δP, we selected a scaled Beta(5,6) distribution, with squared coefficient of variation η2=0.1; the scaling is chosen in each case to achieve a mean value equal to the parameter value. The parameter λ2=KN is given the noise distribution Beta(3,6), which has a slightly higher coefficient of variation η2=0.2. The extrinsic noise values given by pathway reporters match well with those obtained by dual reporters.

MeanSimulation
(r)1-5 λ1μ1λ2=KNKPδPPr (MP)Pr (NP)DR (Mat)
2240020.10.18±0.020.17±0.040.15±0.02
420121020.10.30±0.050.30±0.120.28±0.05
Appendix 5—table 7
A comparison of the pathway-reporter method and dual-reporter method for the multiscale model with Erlang-distributed maturation times and cell-cycle effects (Case (4)A).

Here, the time to maturation, Tmat, is chosen to be consistent with the mean of the stochastic maturation time used in our other models (where the maturation time is exponentially distributed). Specifically, we choose TmatErlang(3,60), with mean length 3/60=0.05, matching our earlier benchmarking using exponentially-distributed maturation time, with mean length 0.05. We consider the rate parameters for the remaining exponentially distributed times to be constant, that is, there is no extrinsic noise beyond that contributed by the cell-cycle effects.

ParametersSimulation
(r)1-5 λ1μ1λ2=KNKPδPPr (MP)Pr (NP)DR (Mat)
2240020.10.04±0.010.02±0.0010.04±0.01
420121020.10.02±0.0030.01±0.010.02±0.01
Appendix 5—table 8
A comparison of the pathway-reporter method and dual-reporter method for the multiscale model with Erlang-distributed maturation times and cell-cycle effects (Case (4)B).

The Erlang distributed maturation time is chosen as in Appendix 5—tables 7. For each of the parameters λ1,μ1,KP,δP, we selected a scaled Beta(5,6) distribution, with squared coefficient of variation η2=0.1; the scaling is chosen in each case to achieve a mean value equal to the parameter value. The parameter λ2=KN is given the noise distribution Beta(3,6), which has a slightly higher coefficient of variation η2=0.2.

MeansSimulation
(r)1-6 λ1μ1λ2=KNTmatKPδPPr (MP)Pr (NP)DR (Mat)
224000.0520.10.21±0.040.18±0.040.23±0.02
42012100.0520.10.34±0.050.30±0.120.36±0.02

Data availability

All methods and simulation results are shared via a github site https://github.com/leham/PathwayReporters (copy archived at https://archive.softwareheritage.org/swh:1:rev:269e0fffe4fc716db6991ccf78ad2191e509c2e1). There is no original data associated with this manuscript.

References

  1. Book
    1. Abramowitz M
    2. Stegun IA
    (1965)
    Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables
    Massachusetts, United States: Courier Corporation.
  2. Book
    1. Olver FW
    2. Lozier DW
    3. Boisvert RF
    4. Clark CW
    (2010)
    NIST Handbook of Mathematical Functions
    Cambridge, United Kingdom: Cambridge University Press.
  3. Book
    1. Ross S
    (2014)
    A First Course in Probability
    London, United Kingdom: Pearson.

Decision letter

  1. Ramon Grima
    Reviewing Editor
  2. Patricia J Wittkopp
    Senior Editor; University of Michigan, United States
  3. Ramon Grima
    Reviewer

In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.

Acceptance summary:

The revisions have substantially improved the manuscript. In particular the authors have clarified the connection of the present work to previous work and shown that their method's accuracy in identifying noise sources does not depend on the details of the gene expression model used. The novel method is based on multiple generic reporters from the same biochemical pathways. Since now it is possible to do simultaneous measurements of transcripts and proteins in single cells, this method offers a viable alternative to the standard dual reporter method, as a means to infer the magnitudes of intrinsic and extrinsic transcriptional noise.

Decision letter after peer review:

[Editors’ note: the authors submitted for reconsideration following the decision after peer review. What follows is the decision letter after the first round of review.]

Thank you for submitting the paper "Pathway dynamics can delineate the sources of transcriptional noise in gene expression" for consideration at eLife. Your article has been reviewed by 3 reviewers, one of whom is a member of our Board of Reviewing Editors, and the evaluation has been overseen by a Senior Editor. Although the work is of interest, we are not convinced that the findings presented have the potential significance that we require for publication in eLife.

Specifically, some of the main critical comments overlapping between one or more reviews are: the need for experimental confirmation of the approach, major issues with the model assumptions and potential overlap with previously published works. Several recommendations have been made by the referees and we hope you find these useful. However it appears that the changes needed are considerable and the paper would need to be rewritten, hence the decision to reject the paper in its current form.

Reviewer #1:

Ham et al. present a manuscript that attempts to answer two main questions:

(i) Is it possible to identify the relative sources of population heterogeneity from

the measured transcript abundance distribution alone?

(ii) Can one develop a method that can reliably estimate the strength of intrinsic and extrinsic noise and that does not require identical and independent pairs of

gene reporters?

In the context of this paper, extrinsic noise is assumed to be variations in parameters across cells; these variations are static and do not depend on time. The authors then proceed to answer question (i) by writing the observed distribution of transcripts as a compound distribution given by Equation 2 which takes as input the distribution of the telegraph model (that is analytically known in closed-form) – the underlying dynamics – and the distribution of parameters across a cellular population – the population heterogeneity. As summarized in Table I, they find that there are various different choices of these 2 distributions that can lead to the same compound (observed) distribution, thus highlighting the issue of non-identifiability if one only has available transcription abundance distributions. I note that the issue of non-identifiability is already known and described in the literature, however the authors identify some cases that were previously unknown. This framework is useful and the results are interesting; however the present formalism does not take into account some of the most important and ubiquitous sources of noise, namely those due to cell division (binomial partitioning of products upon division) and the cell-cycle (replication and cell-cycle length variability) because the authors use the conventional telegraph model that does not account for these phenomena – see Ref 38 where it is shown that these sources of noise can easily mimic that due to transcriptional bursting.

In the second part of the paper, the authors address question (ii). The dual reporter method of Swain et al. (Ref 4) is the standard method to decompose noise into its intrinsic and extrinsic contributions by using two independent reporter genes integrated into the same cell. The assumption that the reporters have identical dynamics is a strong one and presents difficulties in interpreting the results. Ham et al. present an alternative method which does not need dual reporters, rather makes use of 2 pieces of data (for e.g. mRNA and protein) belonging to the same biochemical pathway. The result is based on the decomposition of the covariance of any two variables according to the law of total covariance. This "pathway-reporter method" is shown to be accurate provided there is a small correlation between the reporter pairs – this is a strong assumption, in my view, and it is difficult to make a strong case for it generally. Also when carrying out verification of the method using synthetic data, the correlations between reporter pairs may strongly depend on the details of the model used. For e.g. nascent mRNA maturation to mRNA is in the present paper modeled via a one-step first-order process but more biologically faithful models (see for example PMID: 27667861 and 33976195) model this via a reaction step with a fixed time delay (modelling elongation + termination). The correlations calculated using the latter are likely stronger than the former because the former assumes exponentially distributed times. My main concern however is that the method presented is similar to another one described in the paper PMID: 22529351 which also uses the law of total covariance to understand how the variation in a species is determined by other species in the same pathway. There is also no comparison of the present method with another common one whereby the extrinsic noise is estimated as the part of the expression for the coefficient of variation that is not dependent on the mean molecule numbers; see for e.g. Ref 38 and Ref. 55 (Figure 2B).

Reviewer #1:

I found this paper very interesting to read and I think there is ample scope for the development of new methods in this field. Hence I am generally supportive of this paper.

My main concerns that I would like the authors to address, in order of importance, are:

(i) Novelty – is their method different than PMID: 22529351? They seem to be based on similar mathematical formalisms.

(ii) Nascent mRNA is more faithfully modeled using a delay step rather than a first-order step. Using the latter is approximative at best but will certainly decrease the correlations and maybe making your method perform pathway-reporter better than it would otherwise. The issue might be that the inclusion of such a delay step will lead to a non-Markovian model in which case I am not sure that the same mathematical steps follow

(iii) Cell-cycle and cell division effects are strong sources of noise and it would be ideal that these are at least discussed extensively; even better if they can be integrated into a modified telegraph model and then one can ask the same questions about what one can identify from a compound distribution based on Equation 2.

(iv) Discussion of their method vs the other method of obtaining the extrinsic noise from plots of the coefficient of variation squared versus molecule numbers.

(iv) The manuscript should be checked carefully as I found a number of grammatical mistakes scattered throughout.

Reviewer #2:

The paper deals, theoretically, with the question of decomposing noise into intrinsic and extrinsic components. This has previously been done primarily using the dual reporter method. Here, the authors suggest alternative methods that could in principle bypass the need for using two reporters.

I do not find the results very relevant from the biological perspective. The authors only verify their methods on stochastic simulations, where the model assumptions are specified by hand and thus fit perfectly into their theoretical framework. But how the proposed protocols would perform on experimental data remains unclear. The authors also seem unaware of recent work where essentially the same problem is addressed (Lin and Amir, PRL 2021) – inferring extrinsic vs. intrinsic noise from data – and a method not relying on dual reporters is tested on two existing experimental datasets. Related to the above, it was unclear to me precisely what sort of data would be needed in practice to implement their pathway-reporter method – and whether such data are currently available.

Another example of the potential discrepancy between mathematics and working with actual data is the authors' discussion of the identifiability problem. It seems to me that in reality distributions are always subject to noise (e.g. due to sampling errors) and it is unclear to me whether the sort of rigorous analysis they perform in the paper (that shows whether or not identifiability is possible or not in principle, assuming a perfectly measured distribution) is relevant at all to real data. (e.g.. even if two distributions are mathematically different, would it be possible to distinguish them in practice?)

It seems to me that the paper in its current format is a much better fit for specialized journals in biophysics/mathematics such as the biophysical journal or physical review.

Reviewer #3:

It is now well-known that single-cell expression data exhibits significant cell-to-cell heterogeneity, which can stem from both "extrinsic" factors (e.g. enzyme concentrations, energy content, cellular environment etc.) or "intrinsic" noise (e.g. firing of reactions). Ham et al. rigorously show how extrinsic factors can cause non-identifiability of a model from single time-point gene-expression data, and also lead to incorrect conclusions about the underlying process if they are excluded from the model. Then they propose a novel pathway reporter (PR) approach for quantifying the contributions of the extrinsic and the intrinsic factors to the overall cell-to-cell heterogeneity. These results are illustrated with examples based on the telegraph gene-expression model which includes both nascent and mature mRNA. While the developed results are interesting and mathematically accurate, they come with many caveats, and they do not present a significant advance over existing works on this topic. My reasons for this assessment are mentioned below:

1. Non-identifiability results: The observations regarding non-identifiability of the dynamic parameters from the compound distribution model are not surprising and even though specific examples in Table 1 are nicely worked out, I do not see how they contribute substantially to the existing knowledge on the subject. Essentially non-identifiability will hold for nearly all choices of tilde{p} and f, and identifiability, if it exists, is an exception rather than the norm. This is consistent with what the authors observed in the simulation study outlined at the end of page 10.

2. The Pathway Reporter (PR) scheme: The PR scheme is presented as an alternative to the dual reporter (DR) scheme. However PR also suffers from issues similar to the DR scheme:

– DR assumes conditional independence of reporters. Likewise, PR also requires that such an independence (at least at the level of covariance) in order to argue that the intrinsic noise contribution is zero in the reporter covariance. The paper establishes this approximately for certain parameter regimes by using analytical expressions for the steady-state covariances. However generally such neat expressions are not readily available. It is unclear how one can check that this expected covariance is small. Of course if a fully identified model is available, then one can check it with simulations, but then one can simply compute the extrinsic noise directly.

– Another critical assumption of the PR scheme is that the conditional expectation for the reporters given the extrinsic factors Z, should nicely separate into a function of the common parameters and a function of independent parameters. Even though this works for the specific linear gene-expression model considered in the paper, it is unclear if this would work more generally. Moreover having such a splitting is (up to an additional scaling factor) is not much more general than having identical conditional expectations, as needed by the DR approach.

3. Applicability to experimental studies: As the authors mention in the paper, it is difficult to experimentally construct reporters satisfying the conditions of the DR scheme. However the paper does not discuss how for unknown experimental systems one can identify/construct reporters that satisfy the conditions of the PR scheme. In light of my previous comment, this seems to be as challenging as the DR approach, if one does not have a well-characterized mathematical model.

4. Lack of a meaningful connection: There is very little connection between the two parts of the paper (non-identifiability and PR scheme). A more substantial connection is expected, as in the abstract it is said that: "Here we mathematically formalize this non-identifiability; but we also use this to identify how new experimental set-ups coupled to statistical noise decomposition can resolve this non-identifiability."

An example showing how noise decomposition helped turn a non-identifiable system into an identifiable one is lacking. Moreover, it is unclear how quantification of the extrinsic noise enables resolution of the challenging non-identifiability issue in the compound model.

5. All the results are presented for a specific gene-expression model and its derivatives. As the paper claims to provide a general approach for analyzing noise sources, more examples need to be provided, for a clearer evaluation of the PR method and comparison with the dual reporter method.

6. The paper only considers single time-point snapshot data measured at steady-state. What happens if multiple time-points, including the transience, is also taken into account. Certainly model identifiability would improve and the PR scheme should still work. The authors should consider such examples in the paper.

7. The description of the dual reporter method is misleading – at several places the paper says that the dual reporter method requires "independent gene-reporters". However only conditional independence is required by the method. If there is full independence, then the covariance (extrinsic noise) would always be zero.

8. On page 21, first paragraph it is mentioned that "For nascent-protein reporters, the normalized intrinsic contribution to the covariance is satisfactorily small (less than 0.05) for (a) high values of δp in unison with (b) low values of λ (less than 1, though lower values are acceptable if δp is small)". However based on the discussion prior/post this statement and also Figure 5, it seems that λ values should be high (not low!) for the normalized intrinsic contribution to the covariance to be small.

https://doi.org/10.7554/eLife.69324.sa1

Author response

[Editors’ note: The authors appealed the original decision. What follows is the authors’ response to the first round of review.]

Reviewer #1:

[…] I found this paper very interesting to read and I think there is ample scope for the development of new methods in this field. Hence I am generally supportive of this paper.

My main concerns that I would like the authors to address, in order of importance, are:

(i) Novelty – is their method different than PMID: 22529351? They seem to be based on similar mathematical formalisms.

These are based on similar mathematical formalisms, but beyond this there is no appreciable overlap. To clarify, the primary similarity is that both papers involve an analysis of the covariance between reporters. In our submission, this analysis is between two reporters in the same pathway. In the article PMID: 2259351, this is between reporters in different genes. While it is true that PMID: 2259351 considers dual reporters (not pathway reporters) from multiple levels of a system, we emphasise that this is not the same as an analysis of covariance between single reporters from multiple levels of a system, as we do. There is certainly relevance of the PMID: 2259351 method in our broader discussion of noise decomposition, and we have now included this in the revised manuscript.

(ii) Nascent mRNA is more faithfully modeled using a delay step rather than a first-order step. Using the latter is approximative at best but will certainly decrease the correlations and maybe making your method perform pathway-reporter better than it would otherwise. The issue might be that the inclusion of such a delay step will lead to a non-Markovian model in which case I am not sure that the same mathematical steps follow

While we understand that this is not the primary concern of the referee, we agree that it is a point deserving consideration. We would first like to address the validity of the model for nascent mRNA maturation used in the submitted manuscript. This model is standard in the literature and the article PMID: 27667861 states that the model we use is a “reasonable interpretation of the release kinetics" of the maturation process. PMID: 27667861 also notes that the release kinetics are currently not well-understood but that the nascent mRNA distribution is insensitive to whether the maturation time is exponentially distributed or is modelled as a deterministic time delay, provided that the transcription rate is significantly greater than the switching (on/off rates) rates; this assumption is generally the case experimentally (refer to Section 14.3 and Figure S3 of the Supplementary Material). Our expectation then is that the nascent/mature-mRNA pathway reporter scheme will produce similar results irrespective of the details of the model used.

In order to substantiate this intuition, we have now tested the performance of our method for an alternative model of maturation, where the nascent mRNA maturate after a fixed amount of time. We also explore maturation times sampled from Erlang distributions to account for the fact that maturation can involve several shorter stochastic processes. As expected, we find that pathway reporters continue to correctly identify the extrinsic noise contribution across a range of maturation times; a discussion of these results have now been added to the revised manuscript (see the new section titled “Generality of the Pathway Reporter method"). We mention, in addition, that the issue does not arise for the mature-mRNA/protein (or the nascent-mRNA/protein) pathway reporter method in our article, which for at least some reasonable parameter ranges, is sufficient to provide good estimates of extrinsic noise.

Furthermore, we have now additionally considered various other gene expression dynamics, including a nuanced model of transcriptional bursting (capturing polymerase recruitment and pause release), as well as models that account for cell-cycle effects such as gene replication, dosage compensation, binomial partitioning due to cell division, and cell-cycle length variability. In all cases, the pathway reporter method is able to accurately identify noise sources, showing that our approach does not depend on the precise dynamics of mRNA and protein synthesis. In particular, we find that the correlations between pathway reporter pairs do not depend on the details of the gene expression model used. A new section that encompasses these results, as well as addressing the generality of the pathway reporter method has been added to the revised manuscript; please refer to the new section titled “Generality of the Pathway Reporter method" and the supplementary material. Together our results continue to suggest that our method will serve as a useful tool for gene expression analysis.

(iii) Cell-cycle and cell division effects are strong sources of noise and it would be ideal that these are at least discussed extensively; even better if they can be integrated into a modified telegraph model and then one can ask the same questions about what one can identify from a compound distribution based on Equation 2.

Cell-cycle and cell division/replication effects do indeed have a significant impact on gene expression variability. Cell-cycle effects arise predominantly from differences in the transcription rate across the cell population, and can be modelled already by Equation 2 under reasonable assumptions, without modification to the underlying telegraph model. A more explicit treatment of the mechanisms and changes along the cell cycle are challenging to study analytically (see Reference 38 for example), and theoretical models to account for such details are only starting to emerge. However, we agree that an interesting future direction would be to address how taking into account such details in the context of Equation 2 may also affect the inference of biochemical parameters. At this stage, we feel that a more detailed discussion of these effects should be included, and we have now amended the manuscript to reflect this; please refer to the final paragraph of the section titled “The Compound Distribution".

(iv) Discussion of their method vs the other method of obtaining the extrinsic noise from plots of the coefficient of variation squared versus molecule numbers.

The method of obtaining extrinsic noise from such plots is mentioned already in the introduction of the submitted manuscript. Recent work of two of the authors proves that this method is subject to non-identifiability issues (at least from measurements of transcript abundance); see Ham et al., Phys. Rev. Lett, 2020 (Reference 12 of the submitted manuscript). There it is demonstrated that the coefficient of variation squared as a function of copy numbers cannot distinguish the presence of extrinsic noise.

(iv) The manuscript should be checked carefully as I found a number of grammatical mistakes scattered throughout.

The manuscript has been thoroughly checked, and all grammatical mistakes that we can find have been amended.

Reviewer #2:

The paper deals, theoretically, with the question of decomposing noise into intrinsic and extrinsic components. This has previously been done primarily using the dual reporter method. Here, the authors suggest alternative methods that could in principle bypass the need for using two reporters.

I do not find the results very relevant from the biological perspective. The authors only verify their methods on stochastic simulations, where the model assumptions are specified by hand and thus fit perfectly into their theoretical framework. But how the proposed protocols would perform on experimental data remains unclear.

This is a very valid observation in general concerning the translation from mathematical models to practice. The article considers already some existing databases of transcriptomic data, finding that around 90% of the genes lie within the parameter range for reasonable estimates via our method (and certainly identification of the presence of extrinsic noise), with 70% in the range for accurate measurement. We also give discussion of newly emerging experimental approaches required for the most widely applicable of our pathway combinations; these are in the final discussion of the article.

While for some models we have been able to mathematically prove that our method applies, we have now considered the robustness of our pathway reporter method on more nuanced models of gene transcription. More specifically, we consider models incorporating complex maturation processes (at the response of Referee #1), nuanced models of transcriptional bursting (capturing polymerase recruitment and pause release), as well as incorporating the effects of the cell cycle such as gene replication, dosage compensation, binomial partitioning of products due to cell division and cell-cycle length variability. In all cases, the pathway reporter method is able to accurately identify noise sources, showing that our approach does not depend on the precise dynamics of mRNA and protein synthesis. We hope this will help in assuring the referee that the method does not significantly depend on the precise details of the gene expression dynamics, nor on the availability of analytical formulas. A new section that encompasses these results has been added to the revised manuscript; please refer to the section titled “Generality of the Pathway Reporter Method" in the revised version.

The authors also seem unaware of recent work where essentially the same problem is addressed (Lin and Amir, PRL 2021) – inferring extrinsic vs. intrinsic noise from data – and a method not relying on dual reporters is tested on two existing experimental datasets. Related to the above, it was unclear to me precisely what sort of data would be needed in practice to implement their pathway-reporter method – and whether such data are currently available.

We thank the reviewer for drawing our attention to this work. The article by Lin and Amir was yet to appear at the time of submission, which partly explains why the work slipped the attention of the authors. Lin and Amir’s work does offer an alternative approach to decomposing noise, which we now reference in the revised manuscript. The generality of their approach appears comparable to ours, and is demonstrated on what is now a similar range of model variations. The applicability also appears comparable: their approach relies upon lineage (single-cell trajectory) measurements of protein concentration over many generations, and such high throughput measurements have only recently become available. Some of our reporter pairs also require experimental data that is only recently emerging, while others are obtainable from nascent/mature-mRNA measurements, which have been routinely achieved for some time.

Another example of the potential discrepancy between mathematics and working with actual data is the authors' discussion of the identifiability problem. It seems to me that in reality distributions are always subject to noise (e.g. due to sampling errors) and it is unclear to me whether the sort of rigorous analysis they perform in the paper (that shows whether or not identifiability is possible or not in principle, assuming a perfectly measured distribution) is relevant at all to real data. (e.g.. even if two distributions are mathematically different, would it be possible to distinguish them in practice?)

We respectfully disagree that mathematically different distributions in general will not be distinguishable in practice: the calibration of models to data is one of the central pillars of the scientific method, and in gene expression in particular, has lead to significant improvements in the understanding of the underlying mechanisms. The non-identifiability results presented in the manuscript show that the concern expressed by the referee does indeed hold in the case of the analysis of transcript abundance alone. It is a particularly strong instance, as it shows that totally different processes can present identically in experimental data. This is proved using precise distributions, but the manuscript verifies that the effect is robust, and that “for practical purposes any similar distribution will produce a similar effect”; see Figure 2B for example. A direct consequence of our results is that extrinsic noise invariably leads to data appearing more bursty than the true underlying system; we prove this mathematically for the most popular model underlying current gene expression studies, and demonstrate this through simulations for other instances. Discussion of this important observation can be found in the introductory section as well as in the section just mentioned. To ensure that these points are made as clearly as possible, we have now amended some of the phrasing in the manuscript.

It seems to me that the paper in its current format is a much better fit for specialized journals in biophysics/mathematics such as the biophysical journal or physical review.

Reviewer #3:

It is now well-known that single-cell expression data exhibits significant cell-to-cell heterogeneity, which can stem from both "extrinsic" factors (e.g. enzyme concentrations, energy content, cellular environment etc.) or "intrinsic" noise (e.g. firing of reactions). Ham et al. rigorously show how extrinsic factors can cause non-identifiability of a model from single time-point gene-expression data, and also lead to incorrect conclusions about the underlying process if they are excluded from the model. Then they propose a novel pathway reporter (PR) approach for quantifying the contributions of the extrinsic and the intrinsic factors to the overall cell-to-cell heterogeneity. These results are illustrated with examples based on the telegraph gene-expression model which includes both nascent and mature mRNA. While the developed results are interesting and mathematically accurate, they come with many caveats, and they do not present a significant advance over existing works on this topic. My reasons for this assessment are mentioned below:

1. Non-identifiability results: The observations regarding non-identifiability of the dynamic parameters from the compound distribution model are not surprising and even though specific examples in Table 1 are nicely worked out, I do not see how they contribute substantially to the existing knowledge on the subject. Essentially non-identifiability will hold for nearly all choices of tilde{p} and f, and identifiability, if it exists, is an exception rather than the norm. This is consistent with what the authors observed in the simulation study outlined at the end of page 10.

We agree that our non-identifiability results may not be surprising to some audiences (non-identifiability issues have already been reported empirically), and we are open about this in the manuscript. This does not, however, undermine the value in proving that they are true and uncovering previously unknown, and more encompassing, non-identifiability cases. One of the results, for example, has recently formed the basis of the research paper doi.org/10.1101/2020.09.25.312868.

Our non-identifiability results are in many ways more of an entry point to our main noise decomposition methods. However, there is a very important consequence of this preliminary section: (as noted in the comments to Referee #2) the results demonstrate that extrinsic noise invariably leads to data appearing more bursty than the actual underlying system, and that this holds for all parameter regimes of the standard models. We feel this is a particularly important message to be highlighted, especially to an interdisciplinary community where the analysis of transcriptional processes increasingly relies on calibrating models to experimental data.

2. The Pathway Reporter (PR) scheme: The PR scheme is presented as an alternative to the dual reporter (DR) scheme. However PR also suffers from issues similar to the DR scheme:

– DR assumes conditional independence of reporters. Likewise, PR also requires that such an independence (at least at the level of covariance) in order to argue that the intrinsic noise contribution is zero in the reporter covariance. The paper establishes this approximately for certain parameter regimes by using analytical expressions for the steady-state covariances. However generally such neat expressions are not readily available. It is unclear how one can check that this expected covariance is small. Of course if a fully identified model is available, then one can check it with simulations, but then one can simply compute the extrinsic noise directly.

While we used the specific expressions to prove approximate independence (conditional on the cell environment), the approximate independence can also be verified across a range of models, where analytical expressions are not available. Our computational efforts already provided this corroboration to some extent, and we have now explored a broader range of models, with the same result. Please refer to our response of point (5) below for more details of these newly added models.

– Another critical assumption of the PR scheme is that the conditional expectation for the reporters given the extrinsic factors Z, should nicely separate into a function of the common parameters and a function of independent parameters. Even though this works for the specific linear gene-expression model considered in the paper, it is unclear if this would work more generally. Moreover having such a splitting is (up to an additional scaling factor) is not much more general than having identical conditional expectations, as needed by the DR approach.

The Noise Decomposition Principle presented in the paper is a convenient and powerful framework that appears to be general enough to capture a flexible range of reporter pairs to estimate extrinsic noise. The assumptions are not restrictive. We have focused on the linear gene expression model because it is broadly applicable, and we can use the available analytical formulas to provide a mathematical verification of the method. As in the previous response, however, the computational simulations also provide a demonstration of what our method is capable of measuring. We have now tested the theoretical framework in extensive simulations for a range of other, more nuanced, gene expression dynamics, to explore applicability in the absence of analytical solutions. For these results, refer to the new section in the manuscript titled "Generality of the Pathway Reporter Method". The results show that our approach is remarkably insensitive to the specific dynamics of mRNA and protein synthesis, and therefore should be broadly applicable to a range of experimental scenarios.

3. Applicability to experimental studies: As the authors mention in the paper, it is difficult to experimentally construct reporters satisfying the conditions of the DR scheme. However the paper does not discuss how for unknown experimental systems one can identify/construct reporters that satisfy the conditions of the PR scheme. In light of my previous comment, this seems to be as challenging as the DR approach, if one does not have a well-characterized mathematical model.

As discussed, the general approach is verified computationally across a range of models, including those where we cannot rely on precise formulas. Additionally, the models we consider are widely used for identifying modes of gene activity from experimental data. In such situations, our assumptions have been implicitly invoked already, and the consequences of our observations would appear to hold the same validity as other consequences frequently being drawn from experimental data.

4. Lack of a meaningful connection: There is very little connection between the two parts of the paper (non-identifiability and PR scheme). A more substantial connection is expected, as in the abstract it is said that: "Here we mathematically formalize this non-identifiability; but we also use this to identify how new experimental set-ups coupled to statistical noise decomposition can resolve this non-identifiability."

An example showing how noise decomposition helped turn a non-identifiable system into an identifiable one is lacking. Moreover, it is unclear how quantification of the extrinsic noise enables resolution of the challenging non-identifiability issue in the compound model.

We initially show that mRNA measurements alone cannot be used to identify even the presence of extrinsic noise; we then proceed to show that with pathway reporters it can; and moreover, we can measure its strength. We feel that this is a very clear connection. The pathway reporter scheme resolves non-identifiability to the extent that the presence of extrinsic noise can be identified and measured, whereas from copy number data alone it cannot, as we demonstrate in the first part of the manuscript. There is no claim in the article that we can identify the precise noise distribution in a compound model, and we are not sure this can be done in general, though further consideration of that claim will be an interesting topic for future work. It is, we believe, important to be able to demonstrate convincingly and unequivocally the presence of extrinsic noise from experimental data.

5. All the results are presented for a specific gene-expression model and its derivatives. As the paper claims to provide a general approach for analyzing noise sources, more examples need to be provided, for a clearer evaluation of the PR method and comparison with the dual reporter method.

We thank the reviewer for this suggestion, as we believe the addition of a more robust evaluation of the pathway reporter method has substantially improved the paper. We have now considered various other gene expression dynamics, including detailed models of the maturation process (at the suggestion of Referee #1), a more nuanced model of transcriptional bursting (capturing polymerase recruitment and pause release), as well as accounting for the cell-cycle effects such as gene replication, dosage compensation, binomial partitioning of products due to cell division, and cell-cycle length variability. In all cases, the pathway reporter method is able to accurately identify noise sources, showing that our approach does not depend on the precise dynamics of mRNA and protein synthesis. In particular, we find that the correlations between pathway reporter pairs do not depend on the details of the gene expression model used. A new section that encompasses these results, as well as addressing the generality of the pathway reporter method has been added to the revised manuscript. Together our results continue to suggest that our method will serve as a useful tool for gene expression analysis.

6. The paper only considers single time-point snapshot data measured at steady-state. What happens if multiple time-points, including the transience, is also taken into account. Certainly model identifiability would improve and the PR scheme should still work. The authors should consider such examples in the paper.

The extension to multiple time points is certainly of interest to us. In previous and related work, we have considered time course data for cells and the identifiability of model characteristics (PMID:21551095) in contrast to snapshot data. So this is definitely an avenue for further exploration that we intend to follow. But in the submitted manuscript, the focus is on the development of an alternative noise decomposition for snapshot data as such data sets are becoming readily available. The analysis of the transient behaviour of dissipative dynamical systems in the presence of noise requires a much more detailed analysis than we attempt here and is clearly outside the scope of what can be dealt with here.

7. The description of the dual reporter method is misleading – at several places the paper says that the dual reporter method requires "independent gene-reporters". However only conditional independence is required by the method. If there is full independence, then the covariance (extrinsic noise) would always be zero.

We have now amended the manuscript to make this implicit understanding explicit.

8. On page 21, first paragraph it is mentioned that "For nascent-protein reporters, the normalized intrinsic contribution to the covariance is satisfactorily small (less than 0.05) for (a) high values of δp in unison with (b) low values of λ (less than 1, though lower values are acceptable if δp is small)". However based on the discussion prior/post this statement and also Figure 5, it seems that λ values should be high (not low!) for the normalized intrinsic contribution to the covariance to be small.

We thank the reviewer for drawing our attention to this. It appears that this error was introduced at some stage in preparation (it is expressed exactly as the referee suggests in the original BioaRχiv version of the manuscript, for example). We regret this potential cause for confusion, and have carefully checked and corrected all typographical mistakes that we can find in the paper.

https://doi.org/10.7554/eLife.69324.sa2

Article and author information

Author details

  1. Lucy Ham

    School of BioSciences, University of Melbourne, Melbourne, Australia
    Contribution
    Conceptualization, Resources, Software, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing - original draft, Writing - review and editing
    For correspondence
    lucy.ham@unimelb.edu.au
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-5177-4434
  2. Marcel Jackson

    Department of Mathematics and Statistics, La Trobe University, Melbourne, Australia
    Contribution
    Conceptualization, Software, Formal analysis, Validation, Investigation, Methodology, Writing - original draft, Writing - review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-8149-1141
  3. Michael PH Stumpf

    School of Mathematics and Statistics, University of Melbourne, Melbourne, Australia
    Contribution
    Conceptualization, Supervision, Funding acquisition, Methodology, Writing - original draft, Project administration, Writing - review and editing
    For correspondence
    mstumpf@unimelb.edu.au
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-3577-1222

Funding

University of Melbourne (DRM)

  • Lucy Ham
  • Michael PH Stumpf

Volkswagen Foundation (93 062)

  • Michael PH Stumpf

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

The authors gratefully acknowledge Rowan D Brackston for helpful discussions in the early stages of this research. We thank Lior Pachter and Gennady Gorin for fruitful discussions on noise identifiability. We also wish to thank Arjun Raj for providing valuable feedback on this work. LH. and MPHS. were supported by the University of Melbourne Driving Research Momentum initiative.

Senior Editor

  1. Patricia J Wittkopp, University of Michigan, United States

Reviewing Editor

  1. Ramon Grima

Reviewer

  1. Ramon Grima

Publication history

  1. Preprint posted: September 30, 2020 (view preprint)
  2. Received: April 12, 2021
  3. Accepted: October 11, 2021
  4. Accepted Manuscript published: October 12, 2021 (version 1)
  5. Version of Record published: November 22, 2021 (version 2)

Copyright

© 2021, Ham et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 2,246
    Page views
  • 320
    Downloads
  • 5
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Lucy Ham
  2. Marcel Jackson
  3. Michael PH Stumpf
(2021)
Pathway dynamics can delineate the sources of transcriptional noise in gene expression
eLife 10:e69324.
https://doi.org/10.7554/eLife.69324

Further reading

    1. Chromosomes and Gene Expression
    Desislava P Staneva, Stefan Bresson ... Robin C Allshire
    Research Article

    Kinetoplastids are a highly divergent lineage of eukaryotes with unusual mechanisms for regulating gene expression. We previously surveyed 65 putative chromatin factors in the kinetoplastid Trypanosoma brucei. Our analyses revealed that the predicted histone methyltransferase SET27 and the Chromodomain protein CRD1 are tightly concentrated at RNAPII transcription start regions (TSRs). Here we report that SET27 and CRD1, together with four previously uncharacterized constituents, form the SET27 promoter-associated regulatory complex (SPARC), which is specifically enriched at TSRs. SET27 loss leads to aberrant RNAPII recruitment to promoter sites, accumulation of polyadenylated transcripts upstream of normal transcription start sites, and conversion of some normally unidirectional promoters to bidirectional promoters. Transcriptome analysis in the absence of SET27 revealed upregulated mRNA expression in the vicinity of SPARC peaks within the main body of chromosomes in addition to derepression of genes encoding variant surface glycoproteins (VSGs) located in subtelomeric regions. These analyses uncover a novel chromatin-associated complex required to establish accurate promoter position and directionality.

    1. Chromosomes and Gene Expression
    2. Genetics and Genomics
    Meng Huang, Minjie Hong ... Xuezhu Feng
    Research Article Updated

    Histone methylation plays crucial roles in the development, gene regulation, and maintenance of stem cell pluripotency in mammals. Recent work shows that histone methylation is associated with aging, yet the underlying mechanism remains unclear. In this work, we identified a class of putative histone 3 lysine 9 mono/dimethyltransferase genes (met-2, set-6, set-19, set-20, set-21, set-32, and set-33), mutations in which induce synergistic lifespan extension in the long-lived DAF-2 (insulin growth factor 1 [IGF-1] receptor) mutant in Caenorhabditis elegans. These putative histone methyltransferase plus daf-2 double mutants not only exhibited an average lifespan nearly three times that of wild-type animals and a maximal lifespan of approximately 100 days, but also significantly increased resistance to oxidative and heat stress. Synergistic lifespan extension depends on the transcription factor DAF-16 (FOXO). mRNA-seq experiments revealed that the mRNA levels of DAF-16 Class I genes, which are activated by DAF-16, were further elevated in the daf-2;set double mutants. Among these genes, tts-1, F35E8.7, ins-35, nhr-62, sod-3, asm-2, and Y39G8B.7 are required for the lifespan extension of the daf-2;set-21 double mutant. In addition, treating daf-2 animals with the H3K9me1/2 methyltransferase G9a inhibitor also extends lifespan and increases stress resistance. Therefore, investigation of DAF-2 and H3K9me1/2 deficiency-mediated synergistic longevity will contribute to a better understanding of the molecular mechanisms of aging and therapeutic applications.