Gene expression in individual cells is highly variable and sporadic, often resulting in the synthesis of mRNAs and proteins in bursts. Such bursting has important consequences for cell-fate decisions in diverse processes ranging from HIV-1 viral infections to stem-cell differentiation. It is generally assumed that bursts are geometrically distributed and that they arrive according to a Poisson process. On the other hand, recent single-cell experiments provide evidence for complex burst arrival processes, highlighting the need for analysis of more general stochastic models. To address this issue, we invoke a mapping between general stochastic models of gene expression and systems studied in queueing theory to derive exact analytical expressions for the moments associated with mRNA/protein steady-state distributions. These results are then used to derive noise signatures, i. e. explicit conditions based entirely on experimentally measurable quantities, that determine if the burst distributions deviate from the geometric distribution or if burst arrival deviates from a Poisson process. For non-Poisson arrivals, we develop approaches for accurate estimation of burst parameters. The proposed approaches can lead to new insights into transcriptional bursting based on measurements of steady-state mRNA/protein distributions.
The cellular response to fluctuating environments requires adjustments to cellular phenotypes driven by underlying changes in gene expression. Given the inherent stochasticity of cellular reactions, biological circuits controlling gene expression have to operate in the presence of significant noise [1–15]. While noise reduction and filtering is essential for several cellular processes [16], cells can also amplify and utilize intrinsic noise to generate phenotypic diversity that enables survival under stressful conditions [17]. Recent studies have demonstrated the importance of such bet-hedging survival strategies in diverse processes ranging from viral infections to bacterial competence [17]. Quantifying the kinetic mechanisms of gene expression that drive variations in a population of cells will thus contribute towards a fundamental understanding of cellular functions with important applications to human health. Recent experiments focusing on gene expression at the single-cell level have revealed striking differences from the corresponding population-averaged behavior. In particular, it has been demonstrated that transcription in single cells is sporadic, with mRNA synthesis often occurring in bursts followed by variable periods of inactivity [7,18–28]. Such transcriptional bursting can give rise to high variability in gene expression products and to phenotypic variations in a population of genetically identical cells [29–32]. Furthermore, dynamical parameters that characterize transcriptional bursting of key genes can significantly influence cell-fate decisions in diverse processes ranging from HIV-1 viral infections to stem-cell differentiation [17]. Correspondingly, there is significant interest in developing approaches for quantifying parameters related to transcriptional bursting such as frequency and mean burst size. In recent years, multiple studies have provided evidence for bursty synthesis of mRNAs [20–25,33,34] and proteins [35,36]. Experimental approaches in such studies include both steady-state measurements and time-dependent measurements of the mean and variance of gene expression products at the single-cell level. While obtaining time-lapse measurements of bursts at the single-cell level can be challenging, steady-state measurements at the single-cell level are now carried out routinely. It would thus be desirable to develop approaches for making inferences about burst parameters in gene expression using steady-state measurements at the single-cell level. As noted in [37], steady-state measurements of the mean and variance alone cannot be used for estimating burst parameters for general models of gene expression, e. g. when burst arrival is governed by complex promoter-based regulation [38]. Additional insights into processes leading to transcriptional bursting can potentially be obtained using measurements of higher moments. However, analytical results for higher moments of steady-state mRNA and protein distributions in general models of expression have not been obtained so far. The derivation of the corresponding analytical expressions will elucidate how measurement of higher moments can potentially lead to quantification of burst parameters. To address these issues, it is essential to develop and analyze a general class of stochastic models of gene expression. A simple stochastic model that is widely used in analyzing bursting in gene expression is the random telegraph model that takes into account the switching of promoter between transcriptionally active (ON) and inactive (OFF) states [39–41]. This model has been used as the basis for several studies focusing on inferring gene expression parameters based on observations of the mean and variance of mRNA/protein distributions [13,27,42]. In this model, in the limit that we have transcriptional bursting, the arrival of bursts is a Poisson process. Correspondingly, the waiting-time distribution between arrival of mRNA bursts is assumed to be exponential. In general, this assumption is not valid as there are multiple kinetic steps involved in promoter activation [37,43,44]. Recent experiments on mammalian genes [7,45,46] have demonstrated that the waiting-time for arrival of bursts does not have an exponential distribution. In view of these experimental observations, it is natural to ask: Using steady-state measurements, can we infer if the burst arrival process is not a Poisson process? If so, how can we estimate the corresponding burst parameters? Furthermore, in estimating burst size it is commonly assumed that mRNA/protein bursts are geometrically distributed. This assumption, which has been validated by experimental observations for some genes, is derived from the corresponding distribution of bursts in the random telegraph model. However, given the complexity and diversity of gene expression mechanisms, it is possible that several promoters involve multiple rate-limiting steps in the transition from the ON state to the OFF state. In such cases, the transcriptional burst size distribution will not be a geometric distribution. This observation leads to the following question: Can we use steady-state measurements of moments to determine if the burst distribution deviates from a geometric distribution? The aim of this paper is to address the above questions by considering models with general arrival processes for mRNA creation. The paper is organized as follows. First, we introduce a class of gene expression models with general arrival processes leading to mRNA/protein bursts with arbitrary burst distribution. Then we review the mapping from gene expression models to systems studied in queuing theory [43,47,48] and use this mapping to derive steady-state moments for mRNA/protein distributions. The analytical expressions obtained for the steady-state moments are used to develop approaches for estimating burst parameters for general arrival processes. Finally, we use the results obtained to derive conditions relating experimentally measurable quantities that determine if the arrival of mRNA bursts deviates from a Poisson process and if the distribution of mRNA bursts deviates from a geometric distribution.
We consider a general model of gene expression [43] as outlined in Fig 1. In the model, mRNAs are produced in bursts, with f (t) representing a general arrival time distribution for mRNA bursts. The mRNA burst distribution can be arbitrary. Each mRNA then produces proteins with rate kp, and finally, both mRNAs and proteins decay with rates μm and μp, respectively. Note that the model also allows for post-transcriptional regulation since the protein burst distribution from each mRNA can be arbitrary; the only assumption is that each mRNA produces proteins independently. In the limit μp ≪ μm, we can use the bursty synthesis approximation [40] for analyzing protein dynamics. This approximation consists of two steps: 1) obtaining the distribution of proteins produced from each mRNA and 2) assuming that the proteins are produced in instantaneous bursts. The corresponding distribution for the number of proteins created is referred to as the protein burst distribution. A detailed justification of the validity of this approximation has been provided in previous work [40,49]. Let a p (z) = ∑ n = 0 ∞ z n p (n) denote the generating function of the protein burst distribution p (n) produced by a single mRNA, and let A p (z) = ∑ n = 0 ∞ z n P (n) denote the generating function of the protein burst distribution P (n) produced by all the mRNAs in a burst. If we denote by Am (z) the generating function of the mRNA burst distribution, then we have the following relation between the generating functions A p (z) = A m [ a p (z) ]. (1) The above relation follows from the observation that the number of proteins produced in a burst is a compound random variable: the sum of m independent identical random variables, each of which corresponds to the number of proteins produced from a single mRNA in the burst and m itself is a random variable denoting the number of mRNAs produced in the burst. While the analytical results that we derive are valid for general mRNA and protein burst distributions, we will primarily focus on a specific class of burst distributions. Simple kinetic models and the results from multiple experiments indicate that mRNA burst distributions are geometric [35]. Similarly, the burst distribution of proteins produced from a single mRNA is a geometric distribution with mean ⟨pb⟩ = kp/μm. For a geometric distribution with mean ⟨pb⟩, the generating function is given by ap (z) =1[ 1+〈pb〉 (1−z) ]. If we condition the geometric distribution on the production of at least 1 mRNA, then the generating function for the corresponding conditional geometric distribution is given by Am (z) =z[ 1+〈mb〉 (1−z) ] with (1 + ⟨mb⟩) as the mean mRNA burst size. Note that in the limit ⟨mb⟩ → 0, this distribution reduces to exactly 1 mRNA produced in each burst. Thus the conditional geometric distribution provides a unified representation of both Poisson arrival process for mRNAs (⟨mb⟩ → 0) and processes leading to transcriptional bursting (⟨mb⟩ > 0). Consider now the protein burst distribution produced by an underlying conditional geometric mRNA burst distribution with mean (1 + ⟨mb⟩). Using Eq (1), we see that the corresponding generating function of the protein burst distribution is given by A p (z) = 1 1 + ⟨ m b ⟩ ⟨ p b ⟩ (1 - z). This is the generating function for a geometric distribution with mean b = (1 + ⟨mb⟩) ⟨pb⟩), where ⟨pb⟩ = kp/μp represents the mean protein burst size from a single mRNA. Single-cell experiments have demonstrated that the protein burst mean b can be directly measured in some cases [35]. However, if the protein production rate kp is not known, the preceding analysis implies that measurements of protein burst distributions (which determine b) cannot be used to determine the degree of transcriptional bursting (1 + ⟨mb⟩). Since the mean transcriptional burst size is an important parameter characterizing bursting, it is of interest to develop approaches for estimating it based on available experiments. Previous work [50] has argued that the mean transcriptional burst size cannot be determined using measurements of protein burst distributions alone or by using only protein steady-state distributions. It was suggested that combining such measurements can potentially provide a way of determining the mean transcriptional burst size. To explore this possibility, it is necessary to derive analytical results connecting moments of burst and steady-state distributions for general kinetic schemes. To obtain steady-state moments for the model outlined in Fig 1, we invoke the mapping of this gene expression model to systems studied in queueing theory [43,48,51,52]. Broadly speaking, queueing theory is the mathematical theory of waiting lines formed by customers who, arriving according to some random protocol, stay in the system until they receive service from a group of servers. Such queues are typically characterized by specifying a) the stochastic process governing the arrival of customers, b) distribution of number of customers in each arrival, c) the stochastic process governing departure of customers, and d) the number of servers. When the gene expression model in Fig 1 is expressed in the language of queueing theory, individual mRNAs/proteins are the analogs of customers in queueing models. The production of mRNAs/proteins in bursts corresponds to the arrival of customers in batches. Just as the customers leave the queue after receiving service, mRNAs/proteins exit the system upon degradation. Thus the waiting-time distribution for mRNA/protein decay is the analog of service time distribution for customers in queueing models. For the model in Fig 1, their decay time distribution is the exponential distribution. Also, since mRNAs/proteins are degraded independently of each other, the corresponding number of servers in queueing models is ∞ (which ensures that presence of a customer in the system does not affect the service time of any other customer in the system). Based on the above mapping, the queueing system corresponding to the model outlined in Fig 1 is the GIX/M/∞ system [43,48]. In this model, the symbol G refers to a general waiting-time distribution for the arrival process, IX denotes customers arriving in batches of independently distributed random sizes X, M stands for Markovian (i. e. exponential) service-time distribution for customers and ‘∞’ stands for infinite servers. For the GIX/M/∞ model, exact results for iteratively obtaining the moments of the steady-state distribution of the number of customers have been derived [48]. Using these results, explicit expressions for the first four moments of the steady-state distribution are provided in the Supplementary S1 Text. Applying the mapping discussed above, these results can be translated into exact expressions for the moments of mRNA/protein steady-state distributions, as discussed below. Let us first examine the expressions for steady-state means of mRNAs, ⟨ms⟩, and proteins, ⟨ps⟩, which are given by ⟨ m s ⟩ = k b μ m ⟨ m b ⟩, ⟨ p s ⟩ = k b μ p b, (2) where kb stands for the mean arrival rate of mRNA bursts and b = ⟨mb⟩⟨pb⟩ is the mean of the protein burst distribution (including contributions from all the mRNAs). Although Eq (2) has been derived by assuming that the arrival of mRNAs/proteins is a renewal process, it is valid for arbitrary arrival processes. This is because Eq (2) is a direct consequence of Little’s Law [47,53] which is valid for general arrival processes. The above equations, Eq (2), can be used to determine the mean transcriptional burst size, provided the protein burst distribution can be measured experimentally. To see this, we note that dividing the expressions for the mean mRNA and protein levels leads to b ⟨ m b ⟩ = μ p μ m ⟨ p s ⟩ ⟨ m s ⟩. (3) Since the steady-state means ⟨ms⟩ and ⟨ps⟩ as well as the degradation rates μm and μp are parameters that can be measured experimentally, the above equation implies that the ratio b/⟨mb⟩ can be determined experimentally. Given b/⟨mb⟩ = kp/μm, this implies that the mean protein production rate kp can also be determined experimentally. This is an important result since it provides an approach for determining the mean protein production rate kp that is valid for arbitrary arrival processes for mRNAs. Furthermore, the above equation implies that, if the mean of protein burst distribution b can be measured [36], then the mean transcriptional burst size ⟨mb⟩ can also be determined. Thus, if we have measurements for mean mRNA and protein numbers and also the mean of protein burst distribution, then these measurements can be used to determine the degree of transcriptional bursting ⟨mb⟩ as well as the parameters ⟨pb⟩ and kp. It is noteworthy that this procedure for estimating the burst parameters is valid for arbitrary stochastic processes corresponding to mRNA transcription. We next turn to expressions for higher moments of mRNA and protein steady-state distributions. The noise in mRNA steady-state distributions is given by σ m s 2 ⟨ m s ⟩ 2 = 1 ⟨ m s ⟩ + μ m k b + μ m 2 k b [ K g (μ m) - 1 + σ m b 2 ⟨ m b ⟩ 2 - (1 + 1 ⟨ m b ⟩) ], (4) where σ m b 2 is the variance of mRNA burst distribution and Kg (μm) is the so-called gestation factor, K g (μ m) = 1 + 2 [ f L (μ m) 1 - f L (μ m) - k b μ m ], (5) with fL (s) denoting the Laplace transform of arrival time distribution of mRNA bursts. The function Kg (μm) encodes information about the arrival process. Specifically, we note that for Poisson arrivals, we have Kg (μm) = 1. For proteins (in the burst limit μm ≫ μp), we obtain [43] σps2〈ps〉2=1〈ps〉+μpkb+μp2kb[ Kg (μp) −1+σmb2〈mb〉2 − (1−1〈mb〉) + (σpb2〈pb〉2− (1+1〈pb〉) ) 1〈mb〉 ] (6) where Kg (μp) is given by Eq (5) and σ p b 2 is the variance of protein burst distribution produced by a single mRNA. The expression for protein noise is composed of the noise term for the basic two-stage model of gene expression [40] and additive noise contributions due to: a) deviations from exponential waiting-time distribution for the arrival process, b) deviations from conditional geometric distributions for mRNA burst distributions and c) deviations from geometric distributions for protein burst distributions. For both mRNAs and proteins, the noise in steady-state distributions depends on all the moments of the burst arrival time distribution through the term Kg. Therefore, arrival processes corresponding to different kinetic schemes for transcription will make different contributions to the overall noise, even if they have identical means and variances for the the burst arrival time distribution. We note from Eq (4) that, for Poisson arrivals, i. e. Kg = 1, and geometrically distributed burst, i. e. σ m b 2 = 〈 m b 〉 (〈 m b 〉 - 1), the equations for the noise and mean have only two unknown burst parameters, kb and ⟨mb⟩. In this case, experimental measurements of the first two moments of the steady-state distribution are sufficient to estimate the burst parameters, as has been done in multiple studies. However, when the arrival process is non-Poisson or if the burst distribution deviates from a geometric distribution, measurements of the first two steady-state moments are not sufficient for estimating the burst parameters. This observation motivates the need for analytical expressions for the higher moments which we turn to next. We now derive analytical expressions for the third moment, specifically the skewness parameter. For mRNAs, the exact expression for skewness γms is given by γ m s σ m s 3 m s = 1 + ⟨ m s ⟩ ⟨ m b ⟩ K 1 (μ m) + 2 ⟨ m b ⟩ 2 K 2 (μ m, ⟨ m b ⟩) + (σ m b 2 + ⟨ m b ⟩ 2 - ⟨ m b ⟩) K 3 (μ m, ⟨ m b ⟩) + ⟨ m b (m b - 1) (m b - 2) ⟩ 3 ⟨ m b ⟩, (7) where we have defined K 1 (μ m) = K g (2 μ m) - K g (μ m), K 2 (μ m, ⟨ m b ⟩) = K g (μ m) - 1 4 (3 ⟨ m b ⟩ + K g (2 μ m) - 1), K 3 (μ m, ⟨ m b ⟩) = 3 2 ⟨ m b ⟩ + K g (μ m) + K g (2 μ m) 2 - 1. (8) For proteins, we obtain in the burst-limit (μm ≫ μp), γ p s σ p s 3 p s = 1 + (A 1 p) 2 [ ⟨ p s ⟩ b K 1 (μ p) + 2 K 2 (μ p, A 1 p) ] + A 2 p K 3 (μ p, A 1 p) + A 3 p 3 A 1 p, (9) where the functions K 1, K 2, K 2 are given by Eq (8), A k p is defined by A k p = d k A p (z) / d z k | z = 1 and, using Eq (1) [54], we obtain the parameters A k p as: A1p=〈mb〉〈pb〉, A2p=〈mb〉 (σpb2−〈pb〉) + (σmb2+〈mb〉2) 〈pb〉2, A3p=〈pb〉3〈mb (mb−1) (mb−2) 〉+3〈mb (mb−1) 〉〈pb〉 〈pb (pb−1) 〉+〈mb〉 〈pb (pb−1) (pb−2) 〉. (10) Similarly, expressions for higher order moments of protein and mRNA steady-state distributions can be obtained iteratively. The corresponding expressions for the kurtosis are provided in the S1 Text. The analytical results derived above for proteins are exact in the burst limit, which assumes that proteins are produced instantaneously from all the mRNAs in a burst. Going beyond the burst limit (i. e. not limited to μm ≫ μp), exact results for the higher moments of the protein steady-state distribution will, in general, depend on the details of the kinetic scheme for gene expression. However, we can derive approximate analytical expressions for general schemes by requiring that: a) the results reduce to the exact results in the burst limit and b) they match the exact results for the two-stage model of gene expression. For the two-stage model, exact results for the first four moments have been derived by Bokes et. al [55]. Comparing these exact results with our results derived in the burst limit, we observe that results of [55] can be reproduced by a suitable scaling of the burst-size parameters A k p. For example, the exact expression for the noise is obtained by the following scaling [43]: (σ p s 2 ⟨ p s ⟩ 2 - 1 ⟨ p s ⟩) → (σ p s 2 ⟨ p s ⟩ 2 - 1 ⟨ p s ⟩) 1 1 + μ p μ m. (11) Similarly, for the expression for skewness, the parameters A 2 p and A 3 p are scaled as: A 2 p → A 2 p 1 1 + μ p μ m and A 3 p → A 3 p 1 (1 + μ p μ m) (1 + 2 μ p μ m). (12) As shown in Fig 2 (for the random telegraph model) the analytical expressions using this approach are in good agreement with results from simulations [56]. It is noteworthy that the results derived are valid for a general class of kinetic schemes of gene expression. For a specific kinetic scheme, we can determine the corresponding waiting-time distribution for the arrival process and the burst distributions for mRNA and proteins. Substituting these results in the equations derived leads to the corresponding expressions for the moments of the steady-state distribution. The results obtained can thus provide insight into how specific kinetic schemes of gene expression (e. g. combining promoter-based regulation and post-transcriptional regulation) can be used to impact the noise and higher moments of steady-state distributions. The results derived for the steady-state moments indicate that, if the burst arrival process is not a Poisson process, then it is no longer accurate to estimate burst parameters based on measurements of mean and variance only, as has been done in previous studies [13]. In the following, we present approaches for estimating burst parameters in the general case. We begin by considering the general kinetic scheme shown in Fig 3. This form for the kinetic scheme is supported by recent experiments in mammalian cells which suggest the presence of multiple rate-limiting steps between transition of the promoter from OFF to ON state [45,57]. However, as observed in these experiments, a promoter in the ON state switches to the OFF state by a single rate-limiting step. We model the promoter switching from OFF to ON state by a general waiting-time distribution, g (t). The switching rate from ON to OFF state is given by β. The analytical expressions derived for the steady-state moments for mRNAs and proteins can also be used to make inferences about the burst arrival process based on steady-state measurements. Since multiple studies assume that the burst arrival process is characterized by an exponential waiting-time distribution, it would be useful to determine if this assumption is invalid using measurements of steady-state distributions. As shown below, we can obtain conditions for the same using the results derived for higher moments. In the following, we will first focus on the cases that the mRNA burst distribution is conditional geometric and the protein burst distribution is geometric, which is consistent with multiple experimental observations. As discussed, choosing the conditional geometric distribution for mRNAs allows us to consider both single mRNA arrivals and geometric mRNA bursts in one framework. Since experiments can provide measurements of both mRNA and protein steady-state distributions, it is useful to have conditions for the arrival process using either mRNA data or protein data or both mRNA and protein data. Based on these three possibilities, we present three different conditions in the following. As discussed in the previous section, it is widely assumed that the mRNA burst distribution can be represented by a conditional geometric distribution (i. e. including both single mRNA arrivals and geometrically distributed burst arrivals). While this assumption is consistent with multiple experimental observations, for general kinetic schemes the possibility of non-geometric mRNA burst distributions has to be considered. Thus, it is of interest to examine if the results obtained can be used to determine if the mRNA burst distribution deviates from a conditional geometric distribution. To address the possibility of non-geometric mRNA burst distributions, let us first consider that the random variable corresponding to the mRNA burst distribution (mb) has a conditional geometric distribution. That is, the probability that a burst produces n mRNA molecules is given by P (m b = n) = (1 - p) n - 1 p, (29) where 0 < p ≤ 1, and n = 1,2, 3…∞. This distribution leads to σ m b 2 = ⟨ m b ⟩ (⟨ m b ⟩ - 1). (30) Using Eqs (2) and (30) in Eq (4), and denoting F m = σ m s 2 / 〈 m s 〉 as the Fano factor of mRNA copy number, Eq (4) can be rewritten as F m (μ m) = ⟨ m b ⟩ 2 [ 1 + K g (μ m) ]. (31) Similarly, using the burst size distribution from Eq (29), the skewness in Eq (7) is given by γ m s σ m s 3 ⟨ m s ⟩ = 1 + ⟨ m s ⟩ ⟨ m b ⟩ K 1 (μ m) + 2 ⟨ m b ⟩ 2 K 2 (μ m, ⟨ m b ⟩) + 2 (⟨ m b ⟩ - 1) [ (1 + K 3 (μ m, ⟨ m b ⟩) ) ⟨ m b ⟩ - 1 ]. (32) We note that Eq (32) connects experimentally measurable moments of the steady-state distribution to the parameters Kg (μm), Kg (2μm) and ⟨mb⟩. Furthermore, note that Eq (31) can be recast as Kg (μm) = (2Fm (μm) /⟨mb⟩) −1. Now, considering a change in the degradation rate from μm to 2μm (keeping the mean burst size, ⟨mb⟩ invariant), we obtain K g (2 μ m) = (2 F m (2 μ m) / ⟨ m b ⟩) - 1. (33) Using the above in Eq (32), we get an expression connecting experimentally measurable quantities associated with moments of the mRNA steady-state distribution. The resulting expression is: G m ≡ γ m s σ m s 3 ⟨ m s ⟩ - 1 2 F m (2 μ m) (⟨ m s ⟩ - 1) + F m (1 - 2 ⟨ m s ⟩ + 2 F m (2 μ m) ) = 1. (34) We note that the above expression has been derived by making just one assumption, namely, the mRNA burst distribution is a conditional geometric distribution. The derived expression thus indicates that a combination of experimentally measurable quantities has to deviate from 1 if the mRNA burst distribution deviates from a conditional geometric distribution. Thus the analytical results derived provide a signature for deviation from conditional geometric mRNA bursts using measurements of the first three moments of the mRNA steady-state distribution. The main requirement for using the above relation is that measurements of mRNA steady-state distribution can be carried out at two different rates of the mRNAs μm and 2μm. Given that mRNA degradation rates can be tuned experimentally, a straightforward strategy to ensure that the degradation rate is tuned to twice the original value (2μm) is to compare the mean mRNA levels at μm and 2μm. Given these measurements, a value of G m ≠ 1 implies that bursts are not distributed geometrically. The strength of this result lies in the fact that it holds for general arrival processes for mRNA bursts with arbitrary waiting-time distributions. Let us consider a specific simple model to illustrate the condition derived above. First, let the arrival process for mRNA bursts be a Poisson process. For this, arrival time distributions of mRNA bursts in the time domain, t, and in the Laplace domain, s, are given by f (t) = k b e - k b t, f L (s) = k b / (k b + s), (35) where kb is the rate of arrival of mRNA bursts. For the mRNA burst distribution, let us assume that it is given by the negative binomial distribution, i. e. P (m b = n) = (n + r - 1)! n! (r - 1)! p n (1 - p) r, (36) where 0 < p ≤ 1, r ≥ 1, and n = 0,1, 2,3…∞. For r = 1, the above reduces to the geometric distribution and therefore we expect Gm = 1 in this limit. Using the expressions for the moments derived in the previous section, we obtain an explicit expression for Gm (Supplementary S3 Text): G m = 1 3 (- p + 1 p r + 1 + 4 p (r - 1) + 2 + 2). (37) Notice that for the geometric bursts (r = 1) we get G m = 1, as expected. However, for non-geometric bursts, deviations of G m from 1 are observed (also see Fig S3–1 in Supplementary S3 Text). Two additional examples of microscopic models for non-geometric bursts (the two state random telegraph model and a model with three promoter states where mRNAs are produced from two states) are discussed in the. The preceding analysis can be extended to protein steady-state distributions to derive a similar condition for deviations from geometric burst distributions in terms of steady state moments associated with proteins (see Supplementary S4 Text). In this paper we study stochastic gene expression models with a general renewal-type arrival process for mRNAs. By mapping such a generic model of gene expression to systems studied in queueing theory, we derive analytical expressions for the moments for mRNA and protein steady-state distributions. While the focus of this work is on using approaches drawn from queueing theory, it is noteworthy that the kinetic scheme defined in Fig 1 can also be analyzed using the general theory of branching processes with immigration [58]. In future work, it would be of interest to explore potential connections between complementary approaches to such models based on branching processes and queueing theory. While previous studies [37,43] have focused on protein noise, in the present work we derive analytic expressions for higher order moments of both mRNA and protein steady-state distributions. For arbitrary kinetic schemes, the results obtained determine how the moments of steady-state distributions depend on model parameters. They elucidate how different sources (promoter-based regulation, transcriptional bursting, post-transcriptional regulation) combine to determine the overall noise and higher moments. Furthermore, the results derived show how parameters of interest (such as mean protein production rate kp) can be estimated for general models (i. e. without making any assumptions about specific features of the models). The expressions derived for the moments can also be used to infer if the arrival process for mRNAs is non-Poisson or if the mRNA burst distribution deviates from the geometric distribution. Correspondingly, we obtain analytic conditions that provide signatures for non-Poisson arrivals of mRNA bursts and for non-geometric mRNA burst distributions. These conditions involve relations between combinations of of experimentally measurable quantities and can thus be tested by using measurements of either mRNA steady-state distributions or protein steady-state distributions or both. Apart from obtaining insights into the statistics of the arrival process, we can use the results derived for steady-state moments for accurately estimating burst parameters using an iterative approach. Notably, the results and the approaches developed in this work are valid for general models of gene expression i. e. , given the general assumptions made, they do not depend on the specifics of the kinetic schemes. It is important to note that the burst parameter estimation approaches presented in this paper rely on the accurate measurements of higher order moments, such as skewness or kurtosis. This, in turn requires that we have relatively large sample sizes. For example, simulations of two state random telegraph model (see Supplementary S5 Text) indicate that for the standard error in skewness to be below 10%, the sample size should be ∼ 1000. Current experimental limitations on measurements of mRNA distributions (e. g. using RNA FISH) do not allow for such large sample sizes and thus do not lead to accurate computation of skewness or kurtosis. While accurate measurements of higher moments are not readily available in the existing data, it is hoped that our results will provide motivation for carrying out the corresponding experiments in future. The combination of these experimental results with our theoretical approaches can be used in obtaining accurate representations of the arrival process and burst parameters for a wide range of cellular systems.
One of the fundamental problems in biology is understanding how phenotypic variations arise among individuals in a population. Recent research has shown that phenotypic variations can arise due to probabilistic cell-fate decisions driven by inherent randomness (noise) in the process of gene expression. One of the manifestations of such stochasticity in gene expression is the production of mRNAs and proteins in bursts. Bursting in gene expression is known to impact cell-fate in diverse systems ranging from latency in HIV-1 viral infections to cellular differentiation. Recent single-cell experiments provide evidence for complex arrival processes leading to bursting, however an analytical framework connecting such burst arrival processes with the corresponding higher moments of mRNA/protein distributions is currently lacking. We address this issue by invoking a mapping between general models of gene expression and systems studied in queueing theory. The framework developed and the results derived lead to new approaches for testing commonly used assumptions in modeling gene expression and for accurate estimation of burst parameters. Notably, the phenomenon of stochastic bursting has been observed in a wide range of disciplines ranging from neuroscience and finance to cell biology. The approaches developed and results obtained in this work will thus contribute towards quantitative characterization of burst processes in diverse systems of current interest.
Transcriptional Bursting in Gene Expression: Analytical Results for General Stochastic Models