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

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


Introduction
Transcription, the production of RNA from a gene, is an inherently stochastic process. Specifically, the interval of time between two successive transcription events is a random variable whose statistics depend on multiple single-molecule events behind transcription (Sanchez and Golding 2013). When the distribution of this random variable is exponential, we say that expression is constitutive; in that case, the number of transcripts produced in a certain interval of time follows a Poisson distribution. On the other hand, when the distribution of times between two successive transcripts is nonexponential, then the number of transcripts is non-Poissonian. A special case of such non-constitutive behaviour is bursty expression, whereby transcripts are produced in short bursts that are separated by long silent intervals (Suter et al. 2011;Halpern et al. 2015). In yeast, genes whose expression is constitutive include MDN1, KAP104, and DOA1, whereas PDR5 is an example of a gene whose expression is bursty (Zenklusen et al. 2008).
For two decades, mathematical models of gene expression have been developed to predict the distribution of RNA abundance. By matching the theoretical distribution with experimental measurements from microscopy-based methods (Raj et al. 2008), one hopes to obtain insight into the underlying kinetics of transcription and to estimate transcriptional parameters. The standard model of gene expression which has been used for these analyses is the telegraph model (Peccoud and Ycart 1995), whereby a gene can be in two states. Transcription occurs in one of the states, whereupon RNA degrades; first-order kinetics is assumed for all processes. While the distribution obtained from the telegraph model can typically fit cellular RNA abundance data, there are innate difficulties with the interpretation of that fit: fluctuations in cellular RNA numbers and, hence, the shape of the experimental RNA distribution do not only reflect transcription, but also many processes downstream thereof, such as splicing, RNA degradation, and partitioning during cell division.
To counteract these difficulties, in the past few years, mathematical models (Choubey et al. 2015;Choubey 2018;Heng et al. 2016;Cao and Grima 2020) have been developed to predict the statistics of nascent RNA, i.e. of RNA in the process of being synthesised by the RNA polymerase molecule (RNAP), which can be visualised and quantified due to recent advances in fluorescence microscopy (Lenstra et al. 2016;Skinner et al. 2016;Larson et al. 2011;Antoine et al. 2014;Brouwer and Lenstra 2019). In contrast to cellular RNA, the statistics of nascent RNA is a direct reflection of the transcription process; hence, these models can potentially give more insight than the simpler, but cruder telegraph model. Choubey and collaborators (Choubey et al. 2015;Choubey 2018) have developed a stochastic model with the following properties: (i) a gene can be in two states (active or inactive); (ii) from the active state, transcription initiation occurs in two sequential steps: the pre-initiation complex is formed, after which the RNA polymerase escapes the promoter; (iii) once on the gene, the polymerase moves from one base pair to the next (with some probability) until the end of the gene is reached, when transcription is terminated and polymerase detaches. Queuing theory is used to derive analytical expressions for the transient and steady-state means and variances of numbers of RNAP that are attached to the gene in the long-gene limit when the elongation time is practically deterministic. Heng et al. (2016) have considered a coarse-grained version of that model, whereby the movement of RNAP from one base pair to the next is not explicitly modelled, obtaining an analytical expression for the total RNAP distribution in steady-state conditions. More recently, Cao and Grima (2020) have studied a model of eukaryotic gene expression that yields approximate time-dependent distributions of both nascent and cellular RNA abundance as a function of the parameters controlling gene switching, DNA duplication, partitioning at cell division, gene dosage compensation, and RNA degradation; in their coarse-grained model, the movement of RNAP is not explicitly modelled, while the elongation time is assumed to be exponentially distributed, which simplifies the requisite analysis.
The complexity of nascent RNA models has thus far not allowed the same detailed level of analysis as has been possible with the much simpler telegraph model. A few shortcomings of current models can be summarised as follows: (i) distributions of nascent RNA have been derived from models that do not explicitly model the movement of RNAP along a gene (Heng et al. 2016;Cao and Grima 2020), resulting in a disconnect between theoretical description and the microscopic processes underlying transcription; (ii) while the analysis of single-cell sequencing data and electron micrograph data yields the positions of individual polymerases along the gene, allowing for the calculation of statistics (means and variances) of the numbers of RNAP on gene segments that are obtained after binning, detailed models of RNAP elongation (Choubey et al. 2015;Choubey 2018) provide analytical results only for total RNAP on a gene and hence cannot be used to understand gene segment data; (iii) analytical calculations of the statistics of nascent RNA ignore important details of the transcription process such as pausing, traffic jams, backtracking, and premature termination, some of which have to date been explored via stochastic simulation (Klumpp and Hwa 2008;Rajala et al. 2010;Choubey et al. 2015;Rodriguez et al. 2019;Md Zulfikar et al. 2020).
In this paper, we overcome some of the aforementioned shortcomings of analytically tractable models for the transcription process. In Sect. 2, we study a stochastic model for promoter switching and the stochastic movement of RNAP along a gene, allowing for premature termination. We derive exact closed-form expressions for the first and second moments (means and variances) of local RNAP fluctuations on gene segments of arbitrary length, which allows us to study how these statistics vary along a gene as a function of transcriptional parameters; we also obtain expressions for the mean and variance of the total RNAP on the gene which generalise previous work by Choubey et al. (2015). In Sect. 3, we investigate approximations for the distributions of total RNAP and mature RNA, showing in particular that Negative Binomial distributions can provide an accurate approximation in certain biologically meaningful limits. In Sect. 4, we illustrate the difference between the statistics of local and total RNAP fluctuations and those of light fluorescence due to tagged nascent RNA. In Sect. 5, we extend our model to include pausing by deriving approximate expressions for the mean, variance, and distribution of observables. We conclude with a discussion of our results in Sect. 6.

Detailed Stochastic Model of Transcription: Set-up and Analysis
In this section, we specify the stochastic model studied here; then, we derive closedform expressions for the moments of mature RNA and of local and total RNAP fluctuations in various parameter regimes.

Set-up of Model
We consider a stochastic model of transcription that includes the processes of initiation, elongation, and termination, as illustrated in Fig. 1. For simplicity, we divide the gene into L segments; the RNAP on gene segment i is then denoted by P i . The promoter can be either in the inactive state (G off ) or the active state (G on ), switching from the inactive state to the active one with rate s u and from the active state to the inactive one with rate s b . When the promoter is active, initiation commences via the binding of an RNAP with rate r , denoted by P 1 . Subsequently, the RNAP either moves from a gene segment to the neighbouring segment with rate k, or it prematurely detaches with rate d. Note that here we have made two assumptions: (i) the movement of RNAP is unidirectional, away from the promoter site and hence left to right, with no pausing or backtracking allowed; (ii) the detachment and elongation rates are independent of the position of RNAP on the gene. Each RNAP has associated with it a nascent RNA tail that grows longer as the RNAP transcribes more of the gene. When the RNAP reaches the last gene segment, termination occurs, i.e. the RNAP-nascent RNA complex gets dissociated from the gene leading to a mature RNA (M) which degrades with rate d m . Note that for simplicity, we have not considered excluded-volume interaction between adjacent RNAPs here; hence, we make the implicit assumption of low 'traffic', which is plausible when the initiation rate is sufficiently low. (We test the validity of this assumption through simulations below.) Note that, while the choice of L is arbitrary, it should be kept in mind that L needs to be sufficiently large for the dynamics to be described at a fine spatial resolution. However, L also has to be small enough for the length of each gene segment to be much  Figure Online) Model of transcription. a The gene is arbitrarily divided into L segments, with RNAP (blue) on gene segment i denoted by P i . The promoter switches from the active state G on to the inactive state G off with rate s b , while the reverse switching occurs with rate s u . When the promoter is active, initiation of RNAP occurs with rate r . Initiation is followed by elongation, which is modelled as RNAP 'hopping' from gene segment i to the neighbouring segment i + 1 with rate k, i.e. as the transformation of species P i to P i+1 . RNAP prematurely detaches from the gene with rate d. A nascent RNA tail (red), attached to the RNAP, grows as elongation proceeds. Termination is modelled by the change of P L with rate k to mature RNA (M), which subsequently degrades with rate d m . In b, we show the probability distribution P(T ) of the total elongation time T -the time between initiation and termination-as predicted by the stochastic simulation algorithm (SSA; histogram) and our theory (Erlang distribution with shape parameter L and rate k + d; solid line). The parameter values used are L = 50, k = 10/min, and d = 1.5/min. In c, we show the dependence of the mean of the distribution P(T ) on the RNAP detachment rate (d), as predicted by SSA (dots) and our theory ( T = L/(k + d); solid line). The relevant parameter values are L = 50 and k = 10/min (Color figure online) larger than the footprint of an RNAP; the latter is needed to ensure the validity of the low-traffic assumption. The elongation time which is the total time T from initiation to termination, that is, conditioning on those realisations for which the RNAP does not prematurely detach, is Erlang distributed with mean L/(k + d) and coefficient of variation 1/ √ L; see 'Appendix A' for a derivation and Fig. 1b, c for verification through stochastic simulation (SSA).
Note that the total number of RNAPs transcribing the gene is equal to the number of nascent RNA molecules present, irrespective of their lengths; to shed light on the fluctuations of nascent RNA, in this section we therefore focus on the calculation of statistics of local and total RNAP fluctuations. We define the vector of molecule numbers m = (n 0 , n 1 , . . . , n L , n), and we write n 0 , n i (i = 1, 2, . . . , L), and n for the average numbers of molecules of active gene, RNAP, and mature RNA, respectively. The above model can then be conveniently described by L + 2 species interacting via a set of 2L + 4 reactions with the following rate functions: Note that G off is not an independent species; the reason is that the binary state of the gene implies a conservation law, with the sum of the numbers of G on and G off equalling 1. Hence, the number of independent species in the model is L+2. The rate functions f j are the averaged propensities from the underlying chemical master equation (CME); note that, because our reaction network is composed of first-order reactions, these rate functions also equal the reaction rates in the corresponding deterministic rate equations. The description of our model is completed by the (L + 2) × (2L + 4)dimensional stoichiometric matrix S; the element S i j of S gives the net change in the number of molecules of the ith species when the jth reaction occurs. Given the ordering of species and reactions as described in the tables above, it follows that the matrix S has the simple form where i = 2, . . . , L + 1.

Closed-Form Expressions for Moments of Mature RNA and Local RNAP
In this subsection, we outline the derivation of the steady-state means and variances of local RNAP fluctuations (on each gene segment), as well as of mature RNA. Our results are summarised in the following two propositions. respectively.
Proposition 1 can be proved in a straightforward fashion, as follows. Using the underlying CME, one can show from the corresponding moment equations (Warren et al. 2006) that the time evolution of the vector m of mean molecule numbers in a system of zeroth-order or first-order reactions, i.e. with propensities that are linear in the number of molecules, is given by the time derivative d m /dt = S · f . Given the form of the stoichiometric matrix S and of the rate functions f j , as described in Sect. 2.1, it follows that the mean numbers of all species in steady state can be obtained by solving the following system of L + 2 algebraic equations: These equations can easily be solved simultaneously to yield the steady-state value of m , as given in Eq. (2). Proposition 2 Let τ p = 1/(d +k), τ g = 1/(s u +s b ), and τ m = 1/d m be the timescales of fluctuations of RNAP, gene, and mature RNA, respectively, and define the three new parameters Furthermore, let β = s b /s u denote the ratio of gene inactivation and activation rates. Then, the variances and covariances of molecule number fluctuations of active gene, RNAP, and mature RNA are given by and where i, j = 1, . . . , L. Here, δ i j is the Kronecker delta; moreover, where 2 F 1 denotes the generalised hypergeometric function of the second kind (Digital Library of Mathematical Functions 2020a), which is defined as Given the structure of the matrices J and D above, the Lyapunov Eq. (5) can be solved explicitly for the covariance matrix C whose elements are given by Eq. (4). The solution by induction is involved and can be found in 'Appendix B', which proves Proposition 2.

Simplification in Bursty and Constitutive Limits
Bursty limit: We now consider a particular parameter regime-the limit of large initiation rate r and large gene inactivation rate s b such that b = r /s b is constant. Since the fraction of time spent in the active state is η, it follows that the gene is mostly in the inactive state in that limit. During the short periods of time when it transitions to the active state, a burst of initiation events occur; in particular, a mean number b of RNAPs bind to the promoter during activation. Hence, such genes are often termed bursty, since transcription proceeds via sporadic bursts of activity and b is called the mean transcriptional burst size. For r and s b large with b constant, the expressions for the first two moments of RNAP at every gene segment and of mature RNA from Eqs. (2) and (4), respectively, simplify to here, the subscript b denotes the moments in the bursty limit. Moreover, υ k = s u /k, υ m = s u /d m , and h i j = f i j | α→0 denotes the simplified function f i j in the limit of α−→ 0, which is achieved when s b → ∞. We note that the above expressions for the functions h i j are derived from the expressions for f i j that are given in Eq. (B.33), rather than from those in Eq. (4d). The reason is that, in the bursty limit, we have that 1 2α → ∞, in which case the identity in Eq. (B.36) does not hold. The bursty limit in Eq. (B.33) is simply taken by collecting terms that are not dependent on α, since α −→ 0 in that limit.
To test the accuracy of our theory, in Fig. 2 we compare our analytical expressions for the mean of local RNAP numbers, as well as for various measures of local RNAP fluctuations-the coefficient of variation CV, the Fano factor FF, and the Pearson correlation coefficient CC-with those calculated from stochastic simulation using Gillespie's algorithm (SSA) (Gillespie 1977). Simulations are performed for two different scenarios: (i) without volume exclusion, where the footprint of RNAPs is not taken into account; and (ii) with volume exclusion, where RNAPs are treated as solid objects with a footprint of 35 bp, which is the value reported in Md Zulfikar et al. (2020). For our simulations in Fig. 2, we use parameter values characteristic for the gene PDR5 of length 3070 bp, as reported in Zenklusen et al. (2008). Our choice of L = 30 implies that the length of each gene segment is about 100 bp and, hence, that at most 3 RNAPs can fit in each segment when volume exclusion is taken into account. In this case, Gillespie's algorithm is modified such that the initiation and RNAP 'hopping' rates are proportional to the available volume in the gene segment which the RNAP is moving to. That is achieved by rescaling the transcription initiation rate as r → r (1 − n 1 /3) and the RNAP hopping rate from the ith to the (i + 1)th gene segment as k → k(1 − n i+1 /3). Since we use parameters measured for a gene that demonstrates bursty expression (PDR5) (Zenklusen et al. 2008), we test the accuracy of both the exact theory from Eqs. (2) and (4) and the approximate expressions given in Eq. (8).
The perfect agreement between our exact theory (solid lines) and simulation without volume exclusion (dots) provides a numerical validation of that theory. Our approximate theory (dashed lines) also yields a reasonably good approximation; the mismatch can be decreased if the degree of burstiness is increased, i.e. by increasing the parameters r and s b relative to the other rates in the model. We also note that the theory is in good agreement with simulation with volume exclusion (open circles), which shows that the 'low traffic' assumption upon which our theory is based is valid.
The following interesting observations can be made from these figures: (i) if the rate of premature detachment is greater than zero, then the mean of local RNAP decreases monotonically with the distance i from the promoter according to a power law, whereas that mean is constant along the gene if there is no premature detachment, as expected; (ii) the size of RNAP fluctuations, as measured by CV, decreases with i for small premature detachment rates, but increases with i for sufficiently large values of the detachment rate; (iii) the Fano factor approaches 1-the value of FF for a Poissonian distribution-as i increases, which is due to the dispersal of the burst as stochastic elongation proceeds; (iv) the correlation coefficient between the local RNAP on two neighbouring gene segments decreases monotonically with i, which is exacerbated by premature detachment and is a direct result of the stochasticity inherent in the elongation process.
The observation in (iii) can be explained in detail as follows. When the detachment rate is zero, a burst of RNAPs rapidly bind to the promoter, leading to large fluctuations near that site; however, thereafter each RNAP moves distinctly from all others due to stochastic elongation. Hence, the burst is gradually dispersed as elongation proceeds, which implies a decrease in the variance of fluctuations with increasing i. When the detachment rate is nonzero, then the same effect is at play; however, the increase in the variance of fluctuations along the gene is now counteracted by the decrease in mean RNAP numbers, which leads to two types of behaviour: for small i, CV decreases with i, since the variance dominates over the mean, while for large i, the opposite occurs and CV increases with i.

Constitutive limit:
The other common parameter regime is that of constitutive gene expression, where the gene spends most of its time in the active state and transcription is continuous, which corresponds to the limit of very small s b . In that limit, the expressions from Eqs. (2) and (4) simplify to while the covariances Cov(n i , n j ) c and Cov(n i , n) c between the species are zero; here, the subscript c denotes the constitutive limit. This drastic simplification reflects the fact that, in the constitutive limit, the distributions of mature RNA and local RNAP are Poissonian: as the regulatory network is effectively given by ∅ → P 1 → P 2 → ... → P L → M → ∅ then, the result follows directly from the exact solution provided in Jahnke and Huisinga (2007).
To further test the accuracy of our theory, in Fig. 3 we compare our analytical expressions for the mean of local RNAP numbers, as well as for various measures of local RNAP fluctuations, with those calculated from stochastic simulation using Gillespie's algorithm, where we use parameters measured for a gene that demonstrates constitutive expression (DOA1) (Zenklusen et al. 2008). As before, we test the accuracy of both the exact theory given by Eqs.
(2) and (4) and the approximate expressions from Eq. (9). Unsurprisingly, we observe agreement between exact theory (solid lines) and simulation (dots); the mismatch between our approximate theory and simulation is due to the fact that the gene does not spend 100% of its time in the active state-the true constitutive limit-but, rather, s u /(s u + s b ) ≈ 85%. The local mean RNAP number decreases with distance from the promoter, as was the case for bursty expression in the previous subsubsection, which is to be expected. The various measures which depend on the second moments are, however, considerably different: CV increases monotonically with i, independently of the rate of premature detachment, while FF and CC are very close to 1 and zero, respectively; moreover, the latter two measures practically show very little variation along the gene. The lack of transcriptional bursting explains all these effects in a straightforward fashion.
Finally, we remark that the accuracy of our expressions for the mean and variance of mature RNA, as given in Eq.  Table 2 of Zenklusen et al. (2008). The number of gene segments is arbitrarily chosen to be L = 30. The total elongation time T = 4.5 min is also reported for PDR5, described as the synthesis time and denoted by τ in Zenklusen et al. (2008). The elongation rate by definition takes the value of the ratio k = L/ T − d ≈ L/ T , since d k. The detachment rate d is arbitrarily chosen to be d = 0.01/min (red lines and dots) or d = 0.2/min (black lines and dots). Note that, for the SSA, moments are calculated from one long trajectory with a few million time points, sampled at unit intervals (Color figure online) b for parameters typical of the bursty PDR5 gene. The meaning of the dependence of descriptive statistics on L is discussed in the next section.

Closed-Form Expressions for Moments of Total RNAP
While local RNAP fluctuations are measurable in experiment, as discussed in the Introduction, measurements of total RNAP on a gene are typically reported. Hence, in this section, we briefly discuss descriptive statistics of total RNAP fluctuations.
Recalling that n i is the number of RNAP molecules on the ith gene segment, the total number of RNAPs on the gene-arbitrarily divided into L segments-is given by n tot = L i=1 n i . Given Eq.
(2) and (4), the steady-state mean n tot = L i=1 n i (2) and (4); solid lines), the approximate theory in the constitutive limit (Eq. (9); dashed lines), and simulation via Gillespie's stochastic simulation algorithm (SSA; dots), respectively. The parameters are fixed to s u = 0.7/min, s b = 0.12/min and r = 0.14/min, which are characteristic of the DOA1 gene in yeast, as reported in Supplemental Table 2 of Zenklusen et al. (2008). The number of gene segments is arbitrarily chosen to be L = 30. The total elongation time T = 2.9 min is also reported for DOA1, described as the synthesis time and denoted by τ in Zenklusen et al. (2008). The elongation rate by definition takes the value of the ratio k = L/ T − d ≈ L/ T , since d k. The detachment rate d is arbitrarily chosen to be d = 0.01/min (red lines and dots) or d = 0.2/min (black lines and dots). Note that, for the SSA, moments are calculated from one long trajectory with a few billion time points, sampled at unit intervals (Color figure online) and the steady-state variance Var(n tot ) = L i, j=1 Cov(n i , n j ) of the total RNAP distribution are given by For a detailed derivation of the variance in Eq. (10), we refer to 'Appendix C'. These expressions for the mean and variance of the total RNAP distribution simplify in the bursty and constitutive limits, as can be seen in 'Appendix D'. The accuracy of Eq. (10) is tested by comparing against stochastic simulation with SSA in Fig. 4c, d. Both mean and variance are seen to increase monotonically with the number of (2) and (4); solid lines) and SSA (dots). In c, d, we show the dependence of the moments of total RNAP on L, as predicted by our exact theory (Eq. (10); solid lines) and SSA (dots). The parameters s u , s b , r , and T are characteristic of the PDR5 gene and are the same as in Fig. 2. The premature detachment rate is chosen to be d = 0.01/min; the elongation rate is then given by k ≈ L/ T . The degradation rate of mature RNA is d m = 0.04/min, which is chosen such that the mean mature RNA is roughly consistent with that reported in Fig. 6(b) of Zenklusen et al. (2008). Note that, for the SSA, moments are calculated from one long trajectory with a few billion time points, sampled at unit intervals gene segments L, as we keep the mean elongation time constant; the mean shows very little dependence on L, while the dependence of the variance is more pronounced. We recall that, while the parameter L is arbitrary in principle, it actually determines the size of fluctuations in the elongation time. Since that time is the sum of L independent exponential variables with mean 1/(k + d) each, it follows that the distribution of the elongation time T is Erlang with mean T = L/(k + d) and coefficient of variation squared equal to 1/L. Hence, the larger L is, the narrower is the distribution of T and the more deterministic is elongation itself. Thus, Fig. 4c, d predicts that the mean and variance of total RNAP increase rapidly with decreasing fluctuations in the elongation time T . It hence follows that models in which the elongation rate is assumed to be exponentially distributed (Cao and Grima 2020), which correspond to the case where L = 1 in our model, underestimate the size of nascent RNA fluctuations.

Special Case of Deterministic Elongation
Next, we derive expressions for the descriptive statistics of total RNAP and mature RNA in the limit of large L taken at constant mean elongation time, which corresponds to deterministic elongation. As is shown in Fig. 4, these statistics converge quickly to the ones obtained in the large-L limit; hence, the resulting limiting expressions are likely to be useful across a variety of genes.

Moments of total RNAP distribution:
We define the non-dimensional parameters δ g = τ g /τ d , T g = T /τ g , and T d = T /τ d , which correspond to the ratio of the gene timescale and the polymerase detachment timescale, the ratio of the mean elongation time and the gene timescale, and the ratio of the mean elongation time and the polymerase detachment timescale, respectively; here, τ d = 1/d, as before.
Substituting k → L/ T − d into Eq. (10) and taking the limit of deterministic elongation, i.e. letting L → ∞ at constant T , we obtain the following expressions for the mean, variance, and CV 2 of total RNAP: Here, the subscript ∞ denotes the limit of L → ∞. A detailed derivation of the variance in Eq. (11) can be found in Lemma C.1 of 'Appendix C'. In the special case when RNAP does not prematurely detach from the gene, i.e. for d = 0, the expressions in Eq. (11) simplify to n tot (∞;0) = ηr T , where the subscript (∞; 0) denotes the limit of (L, d) → (∞, 0). The expressions in Eq. (12) have been previously reported in Choubey et al. (2015), where they were derived using queuing theory. Hence, our expressions in Eq. (11) constitute a generalisation of known results, by further taking into account premature detachment of RNAP from the gene. Equation (12) shows that the coefficient of variation squared of total RNAP, denoted by CV 2 (∞;0) , can be written as the sum of two terms: (i) the inverse of the mean which is expected if the distribution of total RNAP is Poissonian, and (ii) a term that increases with increasing β and decreasing T g . Hence, the latter term provides a measure for the deviation of the total RNAP distribution from a Poissonian. In particular, it shows that the deviation is significant in genes for which (i) the fraction of time spent in the inactive state is large (large β), and (ii) the elongation time is much shorter than the switching time between the active and inactive states (small T g ).
Moments of mature RNA distribution: Similarly, in the limit of deterministic elongation, it is straightforward to show that the expressions for the mean and variance of the distribution of mature RNA given by Eqs. (2) and (4) reduce to These expressions can be further simplified in the special case of no premature detachment to read Note that the mean and variance are precisely the same as would be obtained from the telegraph model, for which the corresponding Fano factor in the bursty limit is given by Eq. (16) below. Hence, we anticipate that, in the limit of no premature detachment and deterministic elongation, the distribution of mature RNA from our transcription model is the same as the distribution obtained from the coarser telegraph model. A formal proof of that claim will be given in Sect. 3.

Relationship between Fano factors of total RNAP and mature RNA:
Specifying to the case of no premature detachment, it is interesting to note that in the bursty limit, i.e. for r , s b → ∞ at constant mean burst size b = r /s b in Eq. (12), the Fano factor of total RNAP is given by see also Eq. (D.3) in 'Appendix D'. Here, the subscript n denotes nascent RNA (total RNAP). Eq. (15) is in contrast to the Fano factor of mature RNA in the same bursty limit: see Eq. (D.8) in 'Appendix D', where the subscript m denotes mature RNA. (Note that FF m (b;∞;0) also equals the Fano factor of the telegraph model in the same bursty limit (Raj et al. 2006).) Hence, by comparing Eqs. (15) and (16), we can deduce the following for bursty expression: (i) if the telegraph model is used to estimate the mean transcriptional burst size from total RNAP data where the elongation time is deterministic, then the mean burst size will be overestimated by a factor of two-in other words, the implicit assumption that the elongation time is exponentially distributed is inadequate; (ii) fluctuations in total RNAP (nascent RNA) deviate more from Poisson statistics, for which the Fano factor equals one, than fluctuations in mature RNA. More generally, if we do not enforce the bursty limit, then we find the following relationship between the Fano factors of total RNAP and mature RNA, which are calculated from Eqs. (12) and (14), respectively: Here, Hence, the Fano factor of nascent RNA is larger than that of mature RNA if and only if the above (approximate) condition is satisfied. In the bursty limit, T g → ∞ due to s b → ∞ which, together with T m > 0, implies that Eq. (19) holds; the condition is also satisfied if promoter switching is very fast compared to elongation. By contrast, if T m < 1 and T g < 1, then it is possible to have the opposite scenario where the Fano factor of mature RNA is larger than that of nascent RNA, which occurs, for example, if promoter switching and mature RNA decay are very slow compared to elongation.

Sensitivity of coefficient of variation of total RNAP and mature RNA:
Since we have found explicit expressions for the first two moments of the distributions of total RNAP and of mature RNA, we can now estimate the sensitivity of the noise in each of those to small perturbations in the transcriptional parameters. Specifically, we calculate the logarithmic sensitivity (LS), which is also known as the relativity sensitivity, of the coefficient of variation (CV) to a parameter s, which is defined as s = (s/CV)(∂CV/∂s). (That definition implies that a 1% change in the value of the parameter s results in a change of s % in CV.) In Table 1b, we report the logarithmic sensitivity of the coefficient of variation of total RNAP fluctuations, which is obtained from Eq. (12), to perturbations in the parameters s u , s b , r , and T . Similarly, in Table 1c, we report the logarithmic sensitivity of the coefficient of variation of mature RNA fluctuations from Eq. (14) to perturbations in the parameters s u , s b , r , and d m . In both cases, these sensitivities are calculated for parameter values estimated for five genes in yeast, as reported in Zenklusen et al. (2008); see Table 1a.
The following observations can be made regarding the sensitivity of the noise in total RNAP fluctuations: (i) for the two genes PDR5 and POL1 which spend most of their time in the inactive state due to s b s u , CV is most sensitive to changes in the parameters s u and T ; (ii) for the genes DOA1, MDN1, and KAP104 which spend most of their time in the active state due to s u s b , CV is most sensitive to changes in the parameters r and T ; (iii) the size of mature RNA fluctuations is found to be most sensitive to perturbations in s u and d m for PDR5 and POL1, and to perturbations in r and d m for the other three genes. We furthermore note that for both total RNAP and mature RNA, r is the least sensitive parameter for the genes which are mostly inactive, whereas it is among the most sensitive parameters for genes that are mostly active.  Zenklusen et al. (2008). The degradation rate d m of mature mRNA is estimated from the reported mean number of mature RNA, the parameters s u , s b , r , and Eq. (14) for the mean. (b) Logarithmic sensitivity of CV of total RNAP fluctuations. (c) Logarithmic sensitivity of CV of mature mRNA fluctuations. The most sensitive parameter and the next most sensitive one are marked in dark bold and italic, respectively  such as to infer parameters from data using a Bayesian approach. Moreover, to our knowledge, no exact solutions are known for the distribution of mature RNA in our model. In this section, we aim to devise a simple approximation for the distribution of total RNAP numbers in terms of the Negative Binomial (NB) distribution; these simple distributions have shown great flexibility in describing complex gene expression models with a large number of parameters (Cao and Grima 2020). Finally, by means of singular perturbation theory, we will obtain the distribution of mature RNA under the assumption that RNA polymerase elongation is faster than degradation of mature RNA.

Approximation of Total RNAP Distribution
We approximate the distribution of total RNAP transcribing the gene via a Negative Binomial distribution, as follows. The mean and variance of the Negative Binomial distribution NB(q, p) are given by pq/(1 − p) and pq/(1 − p) 2 , respectively. By assuming that these are equal to the exact mean and variance, respectively, of the total RNAP distribution, see Eq. (10), we obtain effective values for the parameters p and q: In Fig. 6, we show a comparison between the distributions of total RNAP obtained from SSA (dots) and the Negative Binomial approximation in Eq. (20) (solid lines). Our results are presented for two different values of the number of gene segments: L = 1 (exponentially distributed elongation time; left column) and L = 50 (quasideterministic elongation time; right column). Additionally, we rescale our gene inactivation rate as s b → s b , and we present results for three different values of the parameter : 10 −3 , the constitutive limit of the gene being mostly in the active state (top row); 10 −1 , where the gene spends almost equal amounts of time in the active and inactive states, with s b ≈ s u (middle row); and 1, the bursty limit, where the gene spends most of its time in the inactive state (bottom row).
We can make several observations, as follows. For both L = 1 and L = 50, the Negative Binomial approximation performs well for bursting and constitutive expression (top and bottom rows), whereas it is appreciably poor when expression is in between those two limits (middle row). Intuitively, this observation can be explained via the following reasoning. In the limits of the gene being mostly in the active state (constitutive expression) or the inactive state (bursty expression), the distribution of total RNAP is necessarily unimodal. However, when the gene spends a considerable amount of time in each state, the distribution is the sum of two conditional distributions which can manifest either as bimodality or as a wide unimodal distribution, neither of which can be captured by a Negative Binomial distribution. Assuming bursty expression, the Negative Binomial distribution is a more accurate approximation to the distribution obtained from SSA for L = 1 than it is for L = 50; the reason is that L = 1 corresponds to the telegraph model (Raj et al. 2006), in which case it can be proven analytically that the distribution reduces to a Negative Binomial in the limit of bursty expression. For constitutive expression, the Negative Binomial approximation is equally good for L = 1 and L = 50, as the distribution is necessarily Poissonian then and as it is well known that a Negative Binomial distribution can approximate a Poissonian to a high degree of accuracy. In summary, our results hence indicate that Eq. (20) yields a good approximation for the total RNAP distribution of bursty and constitutively expressed genes.
We also note from Fig. 6 that the comparison between the SSA distributions for L = 1 and L = 50, with equal mean elongation times, highlights the importance of modelling elongation with the correct distribution of elongation times for genes that are non-constitutive, i.e. for = 10 −1 or = 1. In particular, if the elongation time is quasi-deterministic (L = 50), there appears to be a significant increase in the probability of observing zero total RNAP transcribing the gene compared to models with an exponentially distributed elongation time (L = 1).  Steady-state distribution of total RNAP and its approximation by a Negative Binomial distribution. We compare the approximation from Eq. (20) (blue lines) with the distribution of total RNAP obtained from stochastic simulation (SSA; red dots). With the exception of s b , the parameters are for the PDR5 gene in yeast and are hence the same as in Fig. 2, with d = 0.01/min. Results are presented for two different values of L, corresponding to an exponentially distributed elongation time (L = 1) and a quasideterministic elongation time (L = 50); k is rescaled such that the two have the same mean elongation time. Additionally, we rescale the gene inactivation rate via s b → s b , where = 10 −3 , 10 −1 , 1, corresponding to constitutive, general, and bursty expression, respectively. (Here, general expression is neither clearly constitutive nor bursty, since the gene spends roughly equal amounts of time in the inactive and active states.) Note that = 1 results in a distribution of nascent RNA that is consistent with that measured for PDR5; the experimental data from Fig. 6 Zenklusen et al. (2008) are plotted for comparison. The Negative Binomial approximation is found to be accurate in the limits of constitutive and bursty expression (top and bottom rows), independently of L

Approximation of Mature RNA Distribution
Next, we apply singular perturbation theory to formally derive the distribution of mature RNA when the elongation rate is much larger than the degradation rate of mature RNA.
We start by defining P j ( n; t) ( j = 0, 1) as the probability of the state n = (n 1 , . . . , n L , n) at time t while the gene is either active (0) or inactive (1). Note that n i is the number of RNAPs on gene segment i for i = 1, . . . , L, while n is the number of mature RNAs. The time evolution of the probabilities P j ( n; t) can be described by a system of coupled CMEs: where . . , n i +c, . . . , n L , n), with c ∈ Z, denotes the standard step operator. We assume that the elongation rate k is faster than the degradation rate d m of mature RNA, i.e. that k/d m 1. Since k = L/ T − d, it follows that in the limit of deterministic elongation (k → ∞), i.e. for L → ∞ at constant mean elongation time T , the condition k/d m 1 is naturally satisfied. In order to find an analytical expression for the propagator probabilities P( n; t) which satisfies the system of CMEs in Eq. (21), we define the probability-generating function as F = j F j , with F j ( z; t) = ∞ n= 0 P j ( n; t) z n ; here, z = (z 1 , . . . , z L , z) is a vector of variables corresponding to the state n. Given the equations for P j ( n; t) from Eq. (21), we obtain the following systems of PDEs for the corresponding generating functions F j ( z; t): where is a differential operator acting on the generating functions F 0 and F 1 . Eq. (22) represents a system of coupled, linear, first-order partial differential equations (PDEs). Now, we introduce the new variables u i = z i − 1 (i = 1, . . . , L) and u = z − 1 to rewrite Eq. (22) as here, the operator in Eq. (23) now takes the form In order to find an analytical solution to Eq. (24), we rescale all rates and the time variable by the decay rate of mature RNA; then, we apply the method of characteristics, with s being the characteristic variable. The first characteristic equation gives d m (dt/ds) = 1, with solution s ≡ t = d m t; hence, we can use the variable t as the independent variable and thus convert the system of PDEs in Eq. (24) into a characteristic system of ordinary differential equations (ODEs), where the overdot denotes differentiation with respect to t . The existence of an integral-form solution to Eq. (26) follows from the fact that the reaction scheme in Fig. 1 contains first-order reactions only. Under the assumption that k d m , we define ε = d m /k; then, we apply Geometric Singular Perturbation Theory (GSPT) (Fenichel 1979;Jones 1995), with 0 < ε 1 as the (small) singular perturbation parameter. We hence separate the system in Eq. (26) into fast and slow dynamics, which will allow us to find an asymptotic approximation for F 0 and F 1 in steady state. A brief introduction to GSPT can be found in 'Appendix E'. Given the above definition of ε, Eqs. (26a) and (26b), the governing equations for u i in the 'slow system', become where u i (i, . . . , L) are the fast variables and u, F 0 , and F 1 are the slow ones. Setting ε = 0 in Eq. (27), we can express the variables u i as u i = μ·u i+1 , with μ = k/(k +d) for i = 1, . . . , L. Finally, we write the variable u 1 as u 1 = μ L ·u. Next, given Eq. (26c), we apply the chain rule, with dt ≡ du · u, to rewrite Eqs. (26d) and (26e) as where the prime now denotes differentiation with respect to u. Solving Eq. (28a) for F 1 and substituting the result into Eq. (28b), we obtain the second-order ODE for F 0 (u). Eq. (29) is a confluent hypergeometric differential equation (Kummer's equation) (Digital Library of Mathematical Functions 2020b) which admits the solution where 1 F 1 denotes the confluent hypergeometric function; here, we consider only one of two independent fundamental solutions of Kummer's differential equation, as we are seeking a solution in steady state where the variable u is bounded. The constant C in Eq. (30) is a constant of integration that is determined from the normalisation condition on the full generating function: Making use of Eq. (31) and applying the normalisation condition F| u=0 = 1, we find that the generating function in steady state reads The probability distribution P(n) of mature RNA can be found from the formula which yields the analytical expression where (·) n is the Pochhammer symbol, as before. Note that the mean and variance of mature mRNA, as calculated from the distribution in Eq. (33), agree exactly with Eqs. (2c) and (4f) in the limit of fast elongation (k → ∞). Note also that the solution in Eq. (33) depends on the parameter μ L , which represents the survival probability of an RNAP molecule, i.e. the probability that RNAP will not prematurely detach from the gene. Finally, we take the limit of deterministic elongation, letting L → ∞ at constant T , which leads to  33) is derived under the assumption that k d m . Note that the distribution is practically independent of L, since Eq. (33) depends on L only through μ L , which for small premature detachment rates d implies μ L ≈ 1 for any L (Color figure online) Note that in the limit of no premature detachment (d = 0), Eq. (34) is precisely equal to the distribution of mature RNA predicted by the telegraph model, which is in wide use in the literature (Raj et al. 2006). Hence, our perturbative approach can be seen as a means to formally derive the conventional telegraph model of gene expression starting from a more fundamental and microscopic model. In Fig. 7, we verify our analytical solution with stochastic simulation for two different genes in yeast. We also note that, for nonzero premature detachment rates (d = 0), Eq. (34) is the steady-state solution predicted by the telegraph model, with parameter r renormalised to r e −d T ; that is to be expected, as the latter is the rate at which RNAPs undergo termination, leading to mature RNAs.

Statistics of Fluorescent Nascent RNA Signal
Thus far, we have determined the statistics of the total number of RNAP transcribing the given gene; these are also the statistics of the number of nascent RNA molecules. However, in experiments using single-molecule fluorescence in situ hybridisation [smFISH (Heng et al. 2016)], molecule numbers of nascent RNA cannot be directly determined. Rather, the experimentally measured RNA 'abundance' is the fluorescent signal emitted by oligonucleotide probes bound to the RNA. Since the length of the nascent RNA grows as RNAP moves away from the promoter, it follows that we must account for the increase in the fluorescent signal as elongation proceeds.
In this section, we take into account these experimental details to obtain closedform expressions for the mean and variance of the fluorescent signal of local and total nascent RNA. We assume that the signal from nascent RNA on the ith gene segment is given by r i = (ν/L)in i for i = 1, . . . , L, where ν is some experimental constant; the value of the parameter (ν/L)i is increasing with i, which models the fact that the fluorescent signal becomes stronger as RNAP moves along the gene. The formula for the mean fluorescent signal at gene segment i is then given by r i = (ν/L)i n i , where n i follows from Eq. (2b); the covariance of two fluorescent signals along the gene, r i and r j (i, j = 1, . . . , L), is given by Cov(r i , r j ) = (ν/L) 2 i jCov(n i , n j ), where Cov(n i , n j ) is obtained from Eq. (4d). In Fig. 8a, b, we plot the mean and Fano factor of the local signal as a function of the gene segment i; note the contrast between the statistics of the fluorescent signal and the corresponding statistics of local RNAP-which is the statistics of nascent RNA-shown in Fig. 2a, c. Similarly, denoting by r tot = L i=1 r i the total fluorescent signal across the gene, we find the following expressions for the steady-state mean r tot = L i=1 r i and the steady-state variance Var(r tot ) = L i, j=1 Cov(r i , r j ): For a detailed derivation of the variance in Eq. (35), see Eq. (F.1) in 'Appendix F'; see also 'Appendix G' for the corresponding expressions in the bursty, constitutive, and deterministic elongation limits. In Fig. 8c, d, we show the mean and Fano factor of the total signal as a function of the number of gene segments (L); as above, we note the contrasting difference between the statistics of the fluorescent signal and the corresponding statistics of total RNAP-which is the statistics of total nascent RNA-shown in Fig. 4c, d. Hence, the calculation of the statistics of the number of nascent RNAs from the raw signal intensity presents a challenge and has to be approached carefully. The expressions presented above allow for the inference of transcriptional parameters from the first two moments of the fluorescent signal by means of moment-based inference techniques (Zechner et al. 2012). Quantitative information about nascent RNA can also be obtained from electron micrograph images (El Hage et al. 2010), which avoids the challenges presented by smFISH.

Model Extension with Pausing of RNAP
Thus far, we have studied a model where RNAPs do not pause as they move along the gene. A natural extension is provided by a modified model in which RNAPs pause along the gene at random sites and elongation is characterised by three processes: forward hopping, pausing, and unpausing of RNAP. The motivation for studying this extended model, which has recently been considered via stochastic simulation in Md Zulfikar et al. (2020), is that experiments have revealed that RNAP exhibits  Fig. 2, as do the rates of elongation and RNAP detachment. The value of the parameter ν is arbitrarily chosen to be ν = 10 pauses of varying duration, typically on the timescale of few seconds (Forde et al. 2002;Adelman et al. 2002).

Closed-Form Expressions for Moments of Local RNAP Fluctuations
We extend the model described in Fig. 1 by assuming that the RNAP on gene segment i can switch between a non-paused (actively moving) state P i and a paused stateP i . The actively moving state P i switches toP i with rate r p , while the reverse reaction occurs with rate r a . Premature detachment from the actively moving RNAP occurs with rate d a , whereas it occurs with rate d p from the paused RNAP. The resulting extended model is illustrated in Fig. 9a. In 'Appendix A', we derive the mean and variance of the corresponding elongation time, which is not Erlang distributed now, as was the case for the model without pausing. Furthermore we find two interesting properties of the coefficient of variation CV 2 T of the elongation time: (i) in the limit of large L at constant mean elongation time, CV 2 T does not tend to zero, which implies (a) (b) Fig. 9 Model of transcription that includes RNAP pausing. In a, we extend the model in Fig. 1 so that it takes into account pausing of RNAP at random segments on the gene. Pausing on gene segment i is modelled by the transition from the active state P i to the paused stateP i with rate r p , while the reverse ('unpausing') transition occurs with rate r a . Premature termination of RNAP occurs with rate d a from the actively moving state, and with rate d p from the paused state. In b, we show the dependence of the coefficient of variation squared (CV 2 T ) of the elongation time distribution on the pausing rate (r p ), as predicted from SSA (dots) and theory (Eq. (A.7); solid lines). Results are shown for two different parameter regimes: D 0 ≡ {d a = 0/min = d p } (no premature polymerase detachment) and D 1 ≡ {d a = 0.05/min = d p } (premature polymerase detachment). The remaining parameters are fixed to L = 50, k = 10/min, and r a = 0.1/min that elongation is not deterministic; (ii) for small rates of premature detachment, CV 2 T is at its maximum when r p ≈ r a , i.e. when RNAP spends roughly half of its time in the paused state. See 'Appendix A' for details and Fig. 9b for a confirmation through stochastic simulation.

Proposition 3 Let the number of RNAP molecules in the active state P i be denoted by n a
i , let the number of molecules in the paused stateP i be n p i , and let the number of molecules of mature RNA be denoted by n. Let σ = r p /r a be the ratio of the pausing and activation rates, let π r a = r a /(r a + d p ) be the probability of RNAP switching to the actively moving state from the paused state, and let π d p = d p /(r a + d p ) be the probability of premature RNAP detachment from the paused state. Furthermore, define the new parametersμ = k/(k + d a + r p π d p ) and λ = σ π r a .
Then, it follows that the steady-state mean number of RNAP molecules in the active and paused states on gene segment i (i = 1, . . . L) is given by n a i = ηρ kμ i and n p i = n a i λ.
Hence, the total mean number of RNAP molecules on each gene segment i reads The proof of Proposition 3 can be found in 'Appendix H'. Note that in the limit of no pausing, i.e. for r p = 0, Eq. (37) reduces to the expression for the mean of RNAP reported in Eq. (2b).
Proposition 4 Let τ r a = 1/r a be the timescale of RNAP activation from the paused state, let τ d p = 1/d p be the timescale of premature termination of paused RNAP, let τ p = 1/(k + d a ) be the typical time that an actively moving RNAP spends on a gene segment, and let τ pp = 1/(r a + d p ) be the typical time spent in the paused state. Furthermore, define the new parameters λ r p = π r p /(1 − π r p ), where π r p = (a) (b) (c) (d) Fig. 10 Dependence of the steady-state probability distributions of total RNAP and mature RNA on the RNAP pausing rate r p for two different genes in yeast. In a, b, we compare the distribution P(n tot ) of the total number of RNAP molecules, as predicted by our model (solid lines), with that obtained from SSA (dots) for yeast genes PDR5 and DOA1, respectively. The model prediction involves fitting a Negative Binomial distribution with a mean and variance given by the closed-form expressions in Eqs. (41) and (42). In c, d, we compare the distribution P(n) of mature RNA, as obtained from singular perturbation theory (Eq. (43); solid lines) with the SSA (dots) for yeast genes PDR5 and DOA1, respectively. Note that for both genes, we keep all parameters fixed (including the elongation rate k) while varying the pausing rate r p to simulate an experiment where the pausing rate can be perturbed directly. The parameters for each gene can be found in Table 1a; we furthermore used L = 50 and fixed k to L/ T , where T is the mean elongation time measured experimentally and reported in Table 1a. Note that the actual mean elongation time is not fixed, as it depends on the pausing rate (r p ) via Eq. (40). The remaining parameters are fixed to r a = 0.1/min, d a = 0.01/min, and d p = 0.03/min. The value of d a is taken from Table 1 in Rajala et al. (2010),where it is reported as the premature termination rate of polymerase in E. coli; the value of d p was chosen to be larger than that of d a to simulate a scenario where premature detachment is enhanced in the paused state. Note that our theory is less accurate for PDR5 than it is for DOA1, as all parameters are very small compared to the elongation rate in the latter case, hence satisfying better the assumptions behind the theory (Color figure online) is the probability of the actively moving RNAP switching to the paused state, as well as ω r a = π r a τ g π r a τ r a + τ g ,α = τ g + λ r p π d p τ g τ g + τ p + λ r p τ g (1 − ω r a ) , and ω = τ g τ pp + τ g .
Assume that the elongation rate is faster than the rates of RNAP pausing, activation, and premature termination, i.e. that k r a , r p , d a , d p . Then, it follows that to leading order in 1/k, asymptotic expressions for the variances and covariances of molecule number fluctuations of active and paused RNAP are given by Cov(n a i , n a j ) = δ i j n a i + n a i n a j αβ · g aa i j , where g aa i j = g aa (i, j) + g aa ( j, i), Cov(n a i , n here, i, j = 1, 2, . . . , L and These results are proved in full in 'Appendix H'. From 'Appendix A', we also have that the mean elongation time in the pausing model is given by Solving Eq. (40) for the elongation rate k, we find that in the limit of L → ∞ taken at constant mean elongation time, k tends to infinity and hence is much larger than r a , r p , d a , and d p , which implies that the results of Proposition 4 hold naturally in that limit.

Approximate Distributions of Total RNAP and Mature RNA
Negative Binomial approximation of total RNAP distribution: We define the total number of RNAP molecules as n tot = L i=1 n i . It then immediately follows from Eq. (37) that the mean of the total RNAP distribution in the pausing model is given by It can also be shown that the variance of total RNAP fluctuations reads see 'Appendix H'. Next, we approximate the distribution of total RNAP by a Negative Binomial distribution whose mean and variance match those just derived, i.e. we consider Eq. (20) with the mean and variance of the total RNAP distribution given by Eqs. (41) and (42) now, respectively. The resulting approximate Negative Binomial distribution is compared with the distribution obtained from SSA in Fig. 10a, b for two different yeast genes, PDR5 and DOA1. The results verify that our approximation is accurate provided the elongation rate k is significantly larger than the other parameters, as assumed in Proposition 4. Perturbative approximation of mature RNA distribution: We can apply singular perturbation theory to formally derive the distribution of mature RNA, assuming that k/d m 1 and r a /d m 1. Following the derivation in Sect. 3.2, we find the following analytical expression for the steady-state probability distribution of mature RNA: see 'Appendix I' for details. Note that the solution in Eq. (43) is dependent on the parameterμ L , which gives the probability that an RNAP molecule does not prematurely detach before termination; see 'Appendix A'. Also, note that in the limit of zero premature termination, i.e. for d a = 0 = d p , Eq. (43) is identical to the distribution of mature RNA predicted by the telegraph model. Finally, by solving Eq. (40) for k, then substituting the resulting expression into Eq. (43) and taking the long-gene limit of L → ∞ at constant T , we obtain that the probability distribution of mature RNA has the same functional form as in Eq. (43), albeit with Note that Eqs. (43) and (44) equal the steady-state solution predicted by the telegraph model, with the initiation rate r renormalised to rμ L or r e −ψ T , respectively. In Fig. 10c, d, we verify the accuracy of our analytical solution using stochastic simulation for two different genes in yeast. Note that a change in the pausing rate r p has relatively little effect on the distribution of mature RNA, as compared to the effect on the distribution of total RNAP; cf. panels (a) and (b) of Fig. 10 in comparison with panels (c) and (d), respectively.

Summary and Conclusion
In this paper, we have analysed a detailed stochastic model of transcription. Our model extends previous analytical work (Choubey et al. 2015;Heng et al. 2016) by (i) taking into account salient processes, such as premature detachment and pausing of RNAP, that were previously not considered analytically; (ii) deriving explicit expressions for the mean and variance of RNAP numbers (nascent RNA) on gene segments as well as on the entire gene; (iii) deriving explicit expressions for the mean and variance of the fluorescent nascent RNA signal obtained from smFISH and identifying differences between the statistics thereof and those of direct measurements of nascent RNA; and (iv) finding approximate distributions of total nascent RNA fluctuations on a gene, The cartoon represents our model in various limits: no pausing (r p = 0), pausing (r p = 0), stochastic elongation (T Erlang distributed), deterministic elongation (T fixed), bursty limit (r , s b → ∞), and premature RNAP detachment (d, d a , d p = 0). We summarise our analytical expressions for the approximate distributions and moments of total RNAP and mature RNA Table 3 Definition of parameters and functions Fraction of time the gene spends in the active state ρ k = r /k Mean number of bound RNAPs in the time 1/k ρ = r /d m Mean number of bound RNAPs in the time Local RNAP survival probability (no-pausing case) Timescale of fluctuations of gene Ratio of the pausing and activation rates π r a = r a /(r a + d p ) Probability of RNAP activation Table 3 continued Probability of premature RNAP detachment from the paused state λ = σ π r a Probability of RNAP pausing from active statẽ μ = k/(k + d a + r p π d p ) Local RNAP survival probability (in pausing case) τ r a = 1/r a Timescale of RNAP activation from the paused state Typical time that an actively moving RNAP spends on a gene segment τ pp = 1/(r a + d p ) Typical time spent in the paused state λ r p = π r p /(1 − π r p ) Ratio of active RNAP timescale over RNAP pausing timescale π r p = r p /(r p + k + d a ) Probability of the actively moving RNAP switching to the paused state ω r a = π r a τ g /(π r a τ r a + τ g ) Non-dimensional parameter Non-dimensional parameter.
without assuming slow promoter switching. A number of interesting observations from our work include the following: (i) When the premature detachment rate of RNAP is nonzero and gene expression is bursty, the coefficient of variation of local RNAP fluctuations can either decrease or increase with distance from the promoter. By contrast, when expression is constitutive, the coefficient of variation increases monotonically with distance from the promoter. Other statistical measures such as the mean, Fano factor, and correlation coefficient of local RNAP numbers decrease monotonically with distance from the promoter. (ii) In the limits of bursty expression, deterministic elongation, and no premature detachment or pausing, the Fano factor of total nascent RNA equals 1+2b, whereas that of mature RNA is 1 + b, where b denotes the mean burst size. An implication is that the telegraph model will result in an overestimate of the mean burst size from nascent RNA data by a factor of 2. Another implication is that deviations from Poisson fluctuations are more apparent in data for nascent RNA than they are for mature RNA. One can further state the following relationship: the Fano factor of nascent RNA equals twice the Fano factor of mature RNA, minus 1. If expression is non-bursty, then the Fano factor of nascent RNA can be larger or smaller than that of mature RNA, as determined by the condition in Eq. (19). (iii) For genes characterised by bursty expression, the sensitivity of the noise in total RNAP fluctuations is highest to perturbations in the gene activation rate and the mean elongation time; for constitutive genes, the most sensitive parameters are the initiation rate and the mean elongation time.
(iv) A Negative Binomial distribution, parameterised with the expressions for the mean and variance of total nascent RNA derived here, provides a good approximation to the true distribution of total nascent RNA fluctuations on a gene when expression is either bursty or constitutive; the approximation is not accurate when the gene spends roughly equal amounts of time in the active and inactive states. We show that the distribution of nascent RNA is highly sensitive to the distribution of elongation times. In particular, if the elongation time is assumed to be exponentially distributed, as is implicitly assumed by telegraph models of nascent RNA, then the probability of observing zero RNA is much lower than if the elongation time is assumed to be fixed. (v) Using geometric singular perturbation theory (GSPT), we have rigorously proven that, in the limit of deterministic elongation (or fast elongation), no pausing and premature detachment, the steady-state distribution of mature RNA in our model is identical to that in the telegraph model (Raj et al. 2006). Consideration of pausing and premature detachment leads to a distribution that can also be obtained from a telegraph model with appropriately renormalised parameters.
A summary of the main theoretical results can be found in Table 2, with all requisite parameters and functions defined in Table 3. The main limiting assumption of our theoretical approach is that the initiation rate is slow enough such that RNAP molecules do not frequently collide with each other while moving along the gene. Hence, the expressions we have derived are reasonable for all but the strongest promoters which are characterised by very fast initiation rates. We anticipate that approximate closedform expressions for the corresponding moments can also be derived when volume exclusion between RNAPs is taken into account by a modification of methods previously devised to understand molecular movement and kinetics in crowded conditions (Cianci et al. 2016;Smith et al. 2017). It is also possible to extend our model by including translation of mature RNA to protein; one can then again apply GSPT to derive distributions for protein numbers in the limit of RNA decaying much faster than protein; however, given item (v) above, we anticipate that the resulting protein distribution will be very similar to those derived from models that do not explicitly take into account nascent RNA (Shahrezaei and Swain 2008;Popović et al. 2016). Further research is required to develop simple approximations of the nascent RNA distribution that are accurate independently of the ratio of gene switching rates. Finally, given the strong recent interest in the development of statistical inference techniques in molecular biology (Gorin et al. 2020;Zechner et al. 2012;Kaan Öcal et al. 2019), we expect that our closed-form expressions for the moments and distributions of nascent and mature RNA will be useful for developing computationally efficient and accurate methods for estimating transcriptional parameters. material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

A. Distribution of Elongation Time
In this section, we answer the following question: what is the distribution of the elongation time, i.e. the time between initiation and termination? In other words, with reference to Fig. 9-which includes the non-pausing model in Fig. 1 as a special casewe want to find the distribution of the time at which RNAP leaves gene segment L (termination) if it was in the active state on gene segment 1 at time t = 0 (initiation).
Let z i (t) be the probability of an RNAP to be on gene segment i in the active state at time t, letz i (t) be the probability of the RNAP to be on gene segment i in the paused state at time t, and let z * i (t) be the probability of the RNAP moving to gene segment i + 1 at time t; note that z * L (t) is the probability of the RNAP falling off the gene and forming a mature RNA, since for i = L, gene segment L + 1 does not exist. Then, it follows from the reaction scheme illustrated in Fig. 9 that the master equations describing the Markovian dynamics on gene segment i are given by Now, we use these equations to find the distribution of the time when RNAP jumps to gene segment i + 1, given that it is on gene segment i in the active state at t = 0, i.e. that z i (0) = 1 andz i (0) = 0. Taking the Laplace transform of Eqs. (A.1a) and (A.1b), we find Solving these equations simultaneously, we obtain Let w(t)dt be the probability that the RNAP moves from segment i to i + 1 in the time interval (t, t + dt). Then, it follows from Eq. (A.1c) that w(t) = ∂ t z * i (t) = kz i (t). Integrating w(t) over all times gives us the probability that the RNAP ultimately moves to the next segment i + 1, Note thatŵ(0) is identical to the parameterμ, as defined in Proposition 3. Let y(t)dt be the probability that the RNAP moves from gene segment i to segment i + 1 in the time interval (t, t + dt), conditioned on those realisations that lead to an RNAP moving to the next gene segment i + 1. (In other words, we exclude those realisations that lead to premature detachment.) Then, it follows by the definition of conditional probabilities that y(t) = w(t)/ŵ(0), which implieŝ . (A.5) It follows that the mean t and variance Var(t) of the time t it takes RNAP to move to the next gene segment are given by respectively. Since RNAP can only move forwards in our model (irreversible motion), it follows that the time it takes an RNAP to move from the ith to the (i + 1)th gene segment is independent of the time taken to move from another, jth segment to the ( j + 1)th segment. Hence, the time required for an RNAP to move across the entire gene from the first to the Lth segment, i.e. the 'elongation' time T from initiation to termination, is a sum of L independent and identical random variables. Thus, we can immediately state that the mean elongation time is T = L t , whereas the variance of the elongation time is Var(T ) = LVar(t). The coefficient of variation squared takes the form From Eq. (A.7), it can be shown that for small premature detachment rates, the coefficient of variation of the elongation time is maximised when r p ≈ r a . Taking the limit of infinitely many gene segments at constant mean elongation time, i.e. solving for k from the expression for the mean elongation time in Eq. (A.6), substituting into Eq. (A.7), and taking the limit of L → ∞, we obtain For the non-pausing model shown in Fig. 1, the above results simplify considerably due to r p = 0 = d p and d a = d; in that case, the inverse Laplace transform of Eq. (A.5) implies that y(t) is an exponential distribution with parameter k + d. Hence, the total time it takes an RNAP to move across the entire gene is the sum of L independent and identically distributed exponential random variables, i.e. an Erlang distribution with shape parameter L and rate k + d, which implies that the mean elongation time is L/(k + d), with coefficient of variation 1/ √ L. It can be seen from Eq. (A.8) that deterministic elongation can only be observed when there is no pausing, i.e. when r p = 0.

B. Solution of Lyapunov Equation
Proof of Proposition 2 We start by defining the symmetric functions f i j = f ji for i, j = 1, . . . , L as where the non-dimensional parameters α, γ , and θ are defined in Proposition 2. The elements of the Lyapunov equation given by Eq. (5) can be written explicitly as a set of simultaneous equations: Now, we substitute the elements of the Jacobian matrix J and the diffusion matrix D from Eqs. (6) and (7), respectively, into the above system of algebraic equations, which we then solve to find the elements of the covariance matrix C. Note that, for the following mathematical derivation, we take into account the expressions for the steady-state mean numbers of species given in Eq.
From Eq. (B.2g), we have that, for j = 4, . . . , L + 1, The proof of Eq. (B.8) is given in Lemma B.1. The above expression for C 2 j can be further simplified to For the proof of the last equality in Eq. (B.9), see Lemma B.2. From Eq. (B.2h), we have that where f 1M is defined in Eq. (B.1). Eqs. (B.2i) through (B.2k) yield the system which can be rewritten more compactly as C i j = δ i j n i−1 + n i−1 n j−1 αβ · f i−1, j−1 for i, j = 3, . . . , L + 1, (B.12) where δ i j is the Kronecker delta. A detailed derivation is given in Lemma B.3. From Eq. (B.2l), we have that for i = 3, . . . , L + 1, where f M M = f L M is defined in Eq. (B.1). Summarising the above results, we conclude that the solution for the symmetric covariance matrix C is given by the system in Eq. (4), where we have that Cov(n i , n j ) = C i+1, j+1 , Cov(n i , n) = C i+1,L+2 for i, j = 0, . . . , L, and Var(n, n) = C L+2,L+2 .
Here, the functions f i j are defined as in Eq. (B.1). Now, the recurrence relation f i j = ( f i−1, j + f i, j−1 )/2 in Eq. (B.1) can be solved for i, j = 1, 2, . . . , L via the method of generating functions, which gives the following analytical expression: see Lemma B.5 for a detailed derivation. Additionally, we can easily prove that the function f i M in Eq. (B.1) can be rewritten as as stated in Eq. (B.8).
Proof The identity in Eq. (B.17) will be proved by induction: one can easily show that it holds for j = 4. Now, we assume that Eq. (B.17) is true for some j ≥ 5; hence, for j + 1, we have as claimed. Hence, the equality in Eq. (B.19) is true for all j = 3, . . . , L.
The proof is by induction: for i = 1, the identity is obvious. We now suppose that Eq. (B.24) is true for some i ≥ 2; hence, for i + 1, we have Proof In order to solve the recurrence relation for the function f i j , we take into account the initial conditions f 00 = 1 and f 0 j = f j0 = α j−1 . Then, we define a generating function g(x, y) via where the last term can be rewritten as (B.28) Hence, Eq. (B.27) becomes which is equivalent to Taking into account the initial conditions, we find that which we substitute into Eq. (B.29) to obtain Making use of the well-known symmetric, bivariate generating function of the binomial coefficients we can rewrite Eq. (B.31) as Rearranging sums in the above expression, we find Hence, we obtain the following exact expression for the function f i j , The expression in Eq. (B.33) can be simplified further due to its symmetry with respect to the indices i and j: we write The function f (i, j) can be further simplified as next, we use the identity where 2 F 1 is again the generalised hypergeometric function of the second kind (Digital Library of Mathematical Functions 2020a). Note that the above identity can be used only when |x| < 1, as the hypergeometric function 2 F 1 is not defined otherwise.
Hence, Eq. (B.35) becomes Given the expression for f (i, j) in Eq. (B.37), one can find the corresponding expression for f ( j, i) by exchanging the indexes i ↔ j.

C. Variance of Total RNAP Distribution
In this section, we derive the exact expression for the variance of the total RNAP distribution, as stated in Eq. (10), which is given by the sum over the covariances Cov(x i , x j ) (i, j = 1, . . . , L), as defined in Eq. (4d). Hence, we have where the function f i j is given in Eq. (10). The first term in Eq. (C.1) equals n tot , the mean of the total RNAP distribution, as stated in Eq. (10); substituting in the expressions for the means n i from Eq. (2b), as well, we obtain Lemma C.1 In the limit of deterministic elongation, i.e. for L → ∞, the expression for Var(n tot ) in Eq. (10) simplifies to which can be further simplified to the expression in Eq. (11).
Proof In order to find the limit of L → ∞ in Eq. (10) (or Eq. (C.2)), we have to evaluate the term L i, j=1 μ i+ j · f i j in that limit. For the following derivation, we consider the function Substituting k → L/ T − d in Eq. (C.4) and taking the limit of L → ∞, we have that G 1 −→ L→∞ 0; hence, Var(n tot ) evaluates to in that limit, which yields the expression in Eq. (C.3), as can easily be verified with the computer algebra package Mathematica. Hence, in the limit of deterministic elongation, the expression for the variance of the RNAP distribution in Eq. (10) reduces to the one in Eq. (11), as claimed.

D. Moments of Total RNAP and Mature RNA in Bursty and Constitutive Limits
Moments of total RNAP in the bursty limit: In the bursty limit, the expressions for the mean and variance of the total RNAP distribution given in Eq. (10) simplify to If, furthermore, we take the limit of deterministic elongation, with L → ∞ at constant T , Eq. (D.1) simplifies to where the subscript (b; ∞) denotes the bursty limit with infinite L. In the limit of zero RNAP detachment, Eq. (D.2) further simplifies to where the subscript (b; ∞; 0) denotes the bursty limit, with L → ∞ and d → 0. Moments of total RNAP in the constitutive limit: In the constitutive limit, Eq. (10) simplifies to If, furthermore, we take the limit of deterministic elongation, i.e. L → ∞ at constant Moments of mature RNA distribution in the bursty limit: In that limit, the closedform expressions in Eq. (8) are given by

E. Introduction to Geometric Singular Perturbation Theory (GSPT)
We consider a system of first-order autonomous ordinary differential equations in the general ('standard') form εẋ = f(x, y, ε), (E.1) y = g(x, y, ε), where (x, y) ∈ R m × R l , with m, l ∈ N. Here, 0 < ε 1 is a (real) singular perturbation parameter, and the overdot denotes differentiation with respect to the 'slow' time t. (Correspondingly, Eq. (E.1) is referred to as the 'slow' system.) The variable x is referred to as the 'fast variable', while y is the 'slow variable'. For simplicity, the functions f : R m × R l × R + → R m and g : R m × R l × R + → R l are assumed to be C ∞ -smooth in all their arguments. In the context of our analysis of the characteristic system in Eq. (26), we have the 'slow system' By comparing the system of equations in Eq. (E.3) with the general form in Eq. (E.1), we see that u i (i = 1, . . . , L) are the fast variables, while u, F 0 , and F 1 are slow. Correspondingly, we have m = L and l = 3 in the above notation, which implies . Now, we introduce a new 'fast' time τ = t/ε, which we substitute into Eq. (E.1) to find the 'fast system' corresponding to Eq. (E.1); here, the prime denotes the derivative with respect to τ . Hence, rewriting Eq. (E.3) in the fast formulation, we find For positive ε, the systems in Eqs. (E.1) and (E.4)-and, correspondingly, the systems in Eqs. (E.3) and (E.5)-are equivalent; however, in the singular limit of ε → 0, we obtain two different systems: setting ε = 0 in Eq. (E.1), we have the 'reduced problem' while we obtain the 'layer problem' for ε = 0 in Eq. (E.4). The 'reduced problem' for the system in Eq. (E.3) implies that the flow of (u, F 0 , F 1 ) is constrained to lie on the (l = 3)-dimensional 'critical manifold' S 0 that is defined by f = 0: where u L+1 ≡ u and (F 0 , F 1 ) are assumed to vary in an appropriately chosen subset of R 2 . From the 'layer problem' of the system in Eq. (E.3), we conclude that y = (u, F 0 , F 1 ) is a parameter which parameterises the (m = L)-dimensional flow of u i = f i (i = 1, . . . , L), the equilibria of which are located on S 0 .
The Jacobian matrix D x f(x, y, 0) of the 'layer problem' corresponding to Eq. (E.5) about S 0 has the eigenvalues Since our definition of the generating function F(z, τ ) in Sect. 3.2 assumed z ∈ [−1, 1], we may restrict to u ∈ [−2, 0] which, by Eq. (E.9), implies that λ i > 0. Hence, the critical manifold S 0 is 'normally hyperbolic'-and, in fact, normally repellingwith an (m + l = L + 3)-dimensional unstable manifold W u (S 0 ). The geometric singular perturbation theory due to Fenichel (1979) thus implies that S 0 will persist, for ε positive and sufficiently small, as a slow manifold' S ε that is (locally) invariant, smooth, and O(ε)-close to S 0 . (As the unstable manifold W u (S 0 ) equals the entire phase space of Eq. (E.3), it trivially persists as the unstable manifold W u (S ε ) for S ε .) In particular, as S 0 is repelling in forward time, it follows that the inverse characteristic transformation corresponding to Eq. (26) is well defined in backward time; details can be found in Veerman et al. (2018), Popović et al. (2016).

F. Variance of Fluctuating Total Fluorescent Signal
By definition, the variance of the total fluorescent signal is given by the sum over all elements Cov(r i , r j ) for i, j = 1, . . . , L, where r i = (ν/L)in i ; the corresponding definitions can be found in Sect. 4 of the main text. Hence, we have that Substituting the expressions for the means n i from Eq. (2b) into Eq. (F.1), we obtain which is the expression stated in Eq. (35).

G. Moments of Fluctuations in Total Fluorescent Signal in Various Limits
Deterministic elongation Substituting k → L/ T −d and taking the long-gene limit of L → ∞ in Eq. (35), we obtain the simplified expressions Var(r tot ) ∞ = r tot ∞ · F 0 + r tot 2 ∞ · βδ g F 1 + F 2 + F 3 2(δ g − 1) 2 (δ g + 1) 2 [1 − (1 + T d )e −T d ] 2 , Here, we note that we are only interested in expressions for the covariances of fluctuations in active and paused RNAP, but not of mature RNA fluctuations; hence, we require closed-form expressions for the elements C i j with i, j = 2L + 2, which we derive by following the same procedure as in 'Appendix B'. Now, we recall that β = s b /s u is the ratio of gene deactivation and activation rates, while τ p = 1/(k + d a ) is the typical time that an actively moving RNAP spends on a gene segment. Additionally, let τ r a = 1/r a be the timescale of RNAP activation from the paused state, let τ d p = 1/d p be the timescale of premature termination of paused RNAP, and let τ pp = 1/(r a + d p ) be the typical time spent in the paused state. Finally, we define the following new parameters: λ r p = π r p /(1 − π r p ), where π r p = r p /(r p + k + d a ) is the probability of actively moving RNAP switching to the paused state, as well as ω r a = π r a τ g π r a τ r a + τ g ,α = τ g + λ r p π d p τ g τ g + τ p + λ r p τ g (1 − ω r a ) , and ω = τ g τ pp + τ g ; (H.6) g aa (i, j) =α from Eq. (H.7); the corresponding solution is then given by g ap i j = ωα j−1 . Finally, the solution of the recurrence relation in Eq. (H.10c) for g pp i j is given by g pp i j = ω(α j−1 +α i−1 )/2. In sum, the leading-order asymptotics (in 1/k) of the covariances between the various RNAP species for k large is hence given by Eq. where the expressions for the corresponding covariances are given in Eq. (39). In order to simplify the above expression, we consider each term on the right-hand side in Eq. (H.12) separately, as follows:

I. Approximation of Mature RNA Distribution in Extended Model
Similarly to Sect. 3.2, we apply geometric singular perturbation theory (GSPT) to formally derive the distribution of mature RNA for the extended pausing model. As was done there, we define P j ( n; t) ( j = 0, 1) as the probability of the state n = (n a 1 , . . . , n a L , n p 1 , . . . , n p L , n) at time t while the gene is either active (0) or inactive (1); then, the time evolution of these probabilities can be described by a system of coupled CMEs: (E n p i − 1)n p i P 0 + d m (E n − 1)n P 0 , − 1)n a i P 1 + k(E n a L E −1 n − 1)n a L P 1 (E n p i − 1)n p i P 1 + d m (E n − 1)n P 1 . (I.1) In order to find analytical expressions for the propagator probabilities P( n; t) which satisfy the system of CMEs in Eq. (I.1), we define the probability-generating functions F j ( z; t), where z = (z a 1 , . . . , z a L , z p 1 , . . . , z p L , z) is a vector of variables corresponding to the state n. Given the equations for P j ( n; t) from Eq. (I.1), we obtain the following system of PDEs for the corresponding generating functions F j ( z; t): L[F 0 ] = s u F 1 − s b F 0 + r (z a 1 − 1)F 0 , L[F 1 ] = s b F 0 − s u F 1 ; (I.2) here, is a differential operator acting on the functions F 0 and F 1 . Eq. (I.2) represents a system of coupled, linear, first-order PDEs. Now, we introduce new variables u a i = z a i − 1, u p i = z p i − 1, and u = z − 1; we also rescale all rates and the time variable with the degradation rate d m of mature RNA. Next, we apply the method of characteristics, with s being the characteristic variable. The first characteristic equation will give us d m (dt/ds) = 1, with solution s ≡ d m t; hence, we can use the variable t = d m t as the independent characteristic variable and thus convert the system of PDEs in Eq. (I.2) into a characteristic system of ODEs: where the overdot denotes differentiation with respect to t. Here, we assume that k/d m 1 and r a /d m 1; hence, we define ε = d m /k as the singular perturbation parameter, and we write d m /r a = εδ, where δ = k/r a = O(1) by assumption. Since 0 < ε 1 is small, we can apply GSPT in order to separate the system in Eq. (I.4) into fast and slow dynamics, which will allow us to find an asymptotic approximation for F 0 and F 1 in steady state. With the above definitions, the governing equations for It follows that u a i and u p i (i = 1, . . . , L) are the fast variables in our system, while u, F 0 , and F 1 are the slow ones; see 'Appendix E'. Setting ε = 0 and solving the system in Eq. (I.5), we find u a 1 =μ L · u, whereμ = k/(k + d a + r p π d p ) has previously been defined in Proposition 3. Now, given Eq. (I.4d), we apply the chain rule, dt ≡ du · u, to rewrite Eqs. (I.4e) and (I.4f) as: where the prime now denotes differentiation with respect to u. The system in Eq. (I.6) is the same as that in Eq. (28), with the substitution μ →μ; hence, following the same derivation as in Sect. 3.2, we conclude that the steady-state analytical expression for the probability distribution of mature RNA is given by