\documentclass{aims}
\usepackage[pagewise]{lineno}
\usepackage{bm}% bold math
\usepackage{dcolumn}
\usepackage{graphicx}
\usepackage[version=4]{mhchem}

\usepackage{txfonts,cite}
\def\typeofarticle{Research article}
\def\currentvolume{21}
\def\currentissue{6}
\def\currentyear{2024}
\def\currentmonth{}
\def\ppages{6225--6262.}
\def\DOI{10.3934/mbe.2024273}
\def\Received{19 November 2023}
\def\Revised{01 March 2024}
\def\Accepted{24 April 2024}
\def\Published{13 June 2024}

\allowdisplaybreaks

\DeclareMathOperator{\var}{var}
\DeclareMathOperator{\erf}{erf}

\makeatletter
\DeclareRobustCommand{\verybig}{\bBigg@{4.5}}
\makeatother

\begin{document}

\title{Analytic delay distributions for a family of gene transcription models}

\author{S. Hossein Hosseini and Marc R. Roussel\corrauth}
% \shortauthors is used in copyright information in the end of the paper
\shortauthors{the Author(s)}
\address{Alberta RNA Research and Training Institute, Department of Chemistry and Biochemistry, University of Lethbridge,
Lethbridge, AB T1K 3M4, Canada }

\corraddr{Email: roussel@uleth.ca; Tel: +1 403 329 2326.}



\begin{abstract}
Models intended to describe the time evolution of a gene network must somehow include transcription, the DNA-templated synthesis of RNA, and translation, the RNA-templated synthesis of proteins. In eukaryotes, the DNA template for transcription can be very long, often consisting of tens of thousands of nucleotides, and lengthy pauses may punctuate this process. Accordingly,  transcription  can last for many minutes, in some cases hours. There is a long history of introducing delays in gene expression models to take the transcription and translation times into account. Here we study a family of detailed transcription models that includes initiation, elongation, and termination reactions. We establish a framework for computing the distribution of transcription times, and work out these distributions for some typical cases. For elongation, a  fixed delay is a good model provided elongation is fast compared to initiation and termination, and there are no sites where long pauses occur. The initiation and termination phases of the model then generate a nontrivial delay distribution, and elongation shifts this distribution by an amount corresponding to the elongation delay. When initiation and termination are relatively fast, the distribution of elongation times can be approximated by a Gaussian. A convolution of this Gaussian with the initiation and termination time distributions gives another analytic approximation to the transcription time distribution. If there are long pauses during elongation, because of the modularity of the family of models considered, the elongation phase can be partitioned into reactions generating a simple delay (elongation through regions where there are no long pauses), and reactions whose distribution of waiting times must be considered explicitly (initiation, termination, and motion through regions where long pauses are likely). In these cases, the distribution of transcription times  again involves a nontrivial part and a shift due to fast elongation processes.
\end{abstract}

%\pacs{87.14.G-,87.10.Mn,82.39.Pj}
\keywords{{gene transcription; delay distribution; transcriptional pausing}}

\setcounter{page}{6225}
\maketitle

\section{Introduction}
 Transcription is a complex process for the biosynthesis of RNA  from a DNA template. Roughly, this process can be divided into three phases: initiation, elongation and termination \cite{Larson-etal2011}.
During initiation, a set of protein factors bind the DNA template upstream of the start site and facilitate the binding and positioning of RNA polymerase \cite{Buratowski-etal1989,Dvir2001209}. The first few nucleotides are added to the nascent RNA, and the polymerase  transitions into productive elongation thereafter. Elongation consists of a sequence of rapid nucleotide addition steps \cite{Darzacq-etal2007}. While each individual step is fast (over $70\,\mathrm{nt\,s^{-1}}$ \cite{Darzacq-etal2007}; nt = nucleotide), genes can be very long. The average human gene, which contains many introns that must be spliced out to form the final transcript, has a length of 28\,000\,nt \cite{Venter-etal2001}, and much longer genes are known \cite{tennyson1995human}, so the elongation time can be significant. Finally, transcription must be terminated. In eukaryotes, redundant mechanisms cooperate to ensure termination \cite{LB2015}.

Transcription is only the first step in gene expression. In eukaryotes, the primary transcript is subjected to processing~\cite{BS2019}, which includes splicing, 5' capping, and 3' polyadenylation; the completed messenger RNA (mRNA) is packaged with proteins into a messenger ribonucleoprotein  (mRNP) complex that is required for export from the nucleus~\cite{Daneholt2001}; it must make its way to the nuclear pore and be exported to the cytoplasm~\cite{SS2013}; it may need to be shuttled to a specific part of the cell~\cite{CDD2020}; and finally it is translated to a protein. Since splicing \cite{Schmidt-etal2011}, 5' capping \cite{Cramer-etal2001} and mRNP assembly \cite{BDS2012} are generally co-transcriptional, a minimal model of gene expression in eukaryotes might consider just transcription, nuclear export, and translation.  All three of these major processes may be separately regulated, so it is not, in general, possible to collapse them into a single protein-synthesis step as is often done in gene expression models.
If we want to model a gene network however, we almost certainly do not want
to model transcription, nuclear export and translation in detail. Rather, we
would prefer to treat each of these processes as effective reactions; one
reaction for transcription, one for export, and one for translation. If we
intend to model the temporal evolution of a gene network, it will be
necessary to consider the time required for each of these processes. Human
genes are very large with, as noted above, an average size of 28\,kb
(kilobases)~\cite{Venter-etal2001}.\footnote{For  practical purposes,
nucleotides and bases are interchangeable in the discussion of  rates of
templated polynucleotide synthesis (transcription, DNA replication, etc.)
and related statistics. Typically, nucleotides are used to count the units
being added to the polynucleotide (RNA, DNA), while bases (usually kb) are
used to measure the template. Since one nucleotide is added for each
complementary base in the template, the rate of synthesis can equivalently
be expressed in nt or in bases per unit time. However, `knt' is, to our
knowledge, never used, so when SI prefixes are to be used, distances along
the template are always measured in bases, hence rates measured in
$\mathrm{nt\,s^{-1}}$ or $\mathrm{kb\,min^{-1}}$.} In the absence of
complicating factors (pauses, co-transcriptional splicing), elongation proceeds at $4.3\,\mathrm{kb\,min^{-1}}$ (approximately $70\,\mathrm{nt\,s^{-1}}$)~\cite{Darzacq-etal2007}, so the transcription elongation phase for a typical human gene would take approximately 6.5\,min. Some genes are larger, and pausing for splicing can add substantially to the transcription time~\cite{tennyson1995human}. Moreover, some events in transcription initiation can be similarly slow~\cite{CP1997,KG2000}. Nuclear export can take 10 to 20 minutes~\cite{KKS2015}. Because the coding sequence is much shorter than the original gene---in humans, the coding sequence is only about 4\% of the gene, the rest being intronic sequences removed by splicing~\cite{Venter-etal2001}---translation will usually be much faster. Accordingly, transcription, including pauses associated with splicing, and nuclear export will tend to dominate the gene expression time.  One way to avoid detailed models of these processes while still considering the expression time of a gene is to incorporate delays \cite{BPM82,BS1987,BM1988,Lewis2003,Monk2003}, which may be fixed or distributed \cite{CLW2022,JP2015,RW2007}, into a model containing effective reactions. But this of course means that we need to know something about the distributions of these delays.

In this paper, we will compute the transcriptional delay distribution, the transcriptional delay being defined here as the delay before the appearance of the product RNA \emph{given that} the gene is active at time zero, in a class of models of varying degrees of complexity. Although the transcription model studied here is similar to prokaryotic models presented elsewhere \cite{roussel2013distribution,roussel2006stochastic}, we focus  on a model and on parameters appropriate to eukaryotes~\cite{Vashishtha}. We further narrow the focus to protein-coding genes transcribed by RNA polymerase II (RNAPII). In eukaryotes, most protein-coding genes are rarely transcribed so that it will be rare for two polymerases to interact \cite{PCP2010}. Even if we take into account that RNAs are often produced in bursts \cite{MCG2012,RPT+2006,Suter-etal2011}, polymerases tend to accelerate as they travel down the gene~\cite{JKL2014,PBS+2018}, which should minimize polymerase-polymerase interactions, at least for polymerases that have translocated sufficiently far along the gene. Accordingly, we study a single-polymerase model.

There have been a number of  models of the kinetics of transcription. Several models have focused on the mechanics and statistical thermodynamics of elongation \cite{bai2007mechanochemical,bai2004sequence,julicher1998motion,wang1998force,yager1991thermodynamic}, and some of these models have yielded insights into the factors responsible for the transcriptional delay.
In prior work, we proposed a single-gene prokaryotic transcription model  from which we obtained analytical and numerical results for  the distribution of the transcriptional delay~\cite{roussel2006stochastic}. We also studied the case where many polymerases transcribe the same gene and obtained stochastic simulation results \cite{roussel2013distribution,roussel2006stochastic,Vashishtha}. A similar approach was followed by von Hippel's group, obtaining time-dependent solutions from which a delay distribution could be inferred~\cite{greive2011development}. Filatova and coworkers have obtained distributions of elongation times as well as means and variances of various stages of transcription for a somewhat simpler transcription model than is considered here, but also considering random activation and deactivation of the gene~\cite{FPG2021}. Boettiger and coworkers also obtained distributions of transcription times focusing on initiation and treating elongation and termination as a single Poisson process~\cite{BRE2011}. Xu and colleagues have examined the dwell time in paused states~\cite{XKK2013}. There have been a number of other modeling studies of transcription of varying complexity and with a variety of objectives \cite{Hamano2016,Klumpp18112008,ribeiro2009delayed,SN2010} (e.g.).  Of course, there is a formal similarity between transcription and translation, which are both templated processes, so we find similar work also in translation modeling. Garai and coworkers have computed both dwell times at a codon, and distributions of translation times \cite{garai2009stochastic,GCR2009}. Perhaps most closely related to the present contribution is the work of Mier-y-Ter{\'a}n-Romero and coworkers, who showed that for a long mRNA, an advection equation is obtained in the continuum limit such that a fixed elongation delay emerges for the process of translation \cite{TSH2010}.

Single-molecule optical trapping experiments show that the movement of RNA
polymerase along the DNA strand is not continuous and is punctuated many
times by pauses \cite{CW2011,Neuman2003}. Mechanistically, there are two
main classes of pauses: (1) backtracking pauses, in which RNAPII slides
backward reversibly  along both the DNA and the RNA, and
 (2) non-backtracking pauses, in which conformational changes in the RNA polymerase active site stop the nucleotide addition cycle \cite{landick2006regulatory,Neuman2003}. Backtracking pauses can be affected by the presence of a trailing RNA polymerase that can restrict how far an RNA polymerase can backtrack \cite{epshtein2003transcription}, leading to density dependence of the transcription rate \cite{klumpp2011pausing}. Voliotis {et al.}\ \cite{voliotis2008fluctuations} studied a model of backtracking transcriptional pauses, for which they obtained a distribution of transcription times. They found that backtracking leads to a heavy-tailed distribution, a conclusion also reached by Klumpp~\cite{klumpp2011pausing}.
 
Pauses can occur in many parts of a gene, and for many different reasons. In eukaryotes, the polymerase almost always pauses just downstream of the start site in an event known as promoter-proximal pausing \cite{LG2011,Nechaev201134,Rahl2010432}. These pauses are thought to have a regulatory role. Pausing can be associated with features of the sequence \cite{bai2004sequence,FXF2021,ZBI+2012}. Splicing can also cause pausing \cite{AIB+2010}. Moreover, pausing  facilitates termination of transcription \cite{GWP2006}.
 
We treat transcription both with and without pauses, and we consider both short ubiquitous pauses and specific pause sites where lengthy pauses occur.
In Section~\ref{sec:Eukaryoticmodel}, we  present a simplified version of  Vashishtha's  model of eukaryotic transcription \cite{Vashishtha} that includes assembly of the pre-initiation complex, initiation, elongation and termination. Vashishtha additionally considered abortive initiation, and  distinguished early elongation (a promoter escape phase with the possibility of promoter-proximal pausing) from productive elongation. With the exception of abortive initiation, the neglected kinetic processes are captured in later sections of this work, albeit in somewhat different form. In Section~\ref{sec:kineticmodel}, we write the  chemical master equations corresponding to each phase of transcription. For both initiation and termination phases, since there are few ordinary differential equations (ODEs), the time-dependent probability   for being in each state of the model can easily be  obtained. For the elongation phase, we again write down the chemical master equation. We find that there are two distinct cases. When the elongation phase is short in duration compared to initiation or termination, which will be the case for many genes for which initiation is slow, the polymerase's motion can be treated as  advective, and the elongation phase can then be treated as introducing a fixed delay between initiation and termination. As noted earlier, the advective treatment developed here is similar in spirit, if not in detail, to the work of Mier-y-Ter{\'a}n-Romero and coworkers on translation~\cite{TSH2010}.  On the other hand, when the duration of the elongation phase becomes significant compared to the slower of initiation or termination, given that most genes are relatively long, an application of the central-limit theorem predicts a Gaussian distribution of elongation times.
In both cases, it is  possible to obtain an analytic expression for the distribution of the total transcription time. 
Because of the modularity of the class of models considered~\cite{greive2011development} such that the total delay is just a sum of delays at each nucleotide, we can treat long pause sites separately, which we do in Section~\ref{sec:pause_in_one_nucleotide}. For simplicity, we go back to the case where elongation is fast, and we consider a single strong pause site, obtaining the distribution of the total transcriptional delay in this case. In all our calculations to this point, we assume that binding of the RNA polymerase to the promoter is irreversible. In Section~\ref{sec:revRNAPbinding}, we relax this assumption and discover that the distribution of transcription times can be  sensitive to the reversibility of RNAP binding when the transition to elongation is slow. Finally, in Section~\ref{sec:Discussion}, we discuss a number of possible extensions and refinements, notably the treatment of pause exit due to the action of a release factor, and the possibility of using our results as a basis for simulating genes expressed in bursts.

This work is situated in a larger literature on the use of delays in chemical and biochemical models. We can think of the states appearing between  initiation and the completion of transcription as intermediates in a synthesis pathway.
The problem of replacing chains repeated steps generating a sequence of
intermediates by distributed delays has  been intensely studied, notably by MacDonald \cite{MacDonald1976,MacDonald1977a,MacDonald89},  and by Barrio and coworkers \cite{BLM2013,LBM2014}. Epstein has additionally considered the issue of ``bottleneck intermediates''~\cite{Epstein90}, of which strong pause sites would be an example. We have already mentioned the contribution of Mier{-}y{-}Ter{\'a}n{-}Romero and coworkers~\cite{TSH2010}, which provided a justification for the use of delays in modeling translation. The distinction made in this paper between cases where elongation is fast and those where elongation is relatively slow adds to our understanding of the proper use of fixed and distributed delays in gene expression models. Moreover, the availability of analytic expressions for the distribution of transcription times in a variety of scenarios will be useful for delay-stochastic simulations of gene expression networks \cite{Bratsun-etal2005,roussel2006validation}.

\section{Eukaryotic transcription model}
\label{sec:Eukaryoticmodel}

We start with a simplified version of Vashishtha's kinetic model for the eukaryotic transcription process~\cite{Vashishtha}, removing abortive initiation, early elongation/promoter escape and promoter-proximal pausing. Promoter-proximal pausing \cite{Jennings2013,KFC2013} can, at least to a first approximation, be modeled as a distinct, long pause, as we will do in section~\ref{sec:pause_in_one_nucleotide}. Similarly, the early elongation phase, during which the polymerase adds nucleotides more slowly than it does in the productive elongations phase \cite{HBG+2006,SHS+2014}, could be treated analogously to productive elongation, but with different kinetic parameters. Abortive initiation~\cite{SJG+2018} represents an alternative pathway that can be included in a transcription model~\cite{Vashishtha}. The effect of abortive initiation would be to reduce the frequency of initiation. Since we focus here on events that lead to the production of a protein, we leave this out. We similarly leave the other refinements mentioned in this paragraph to later work. Models of a similar degree of complexity are sufficient to analyze some experimental datasets~\cite{LHE2021}. Moreover, note that the simplest model considered in this contribution is similar to previously studied models of bacterial transcription \cite{roussel2013distribution,roussel2006stochastic}. With some  reinterpretation of symbols, the results obtained here could thus be applied to bacterial models, and perhaps also to structurally similar translation models \cite{Kremling2007,TSH2010,MP2013}.

Before describing the transcription model, it is useful to think about the context in which the delay distributions derived in later sections are likely to be used. In delay-stochastic simulations \cite{Bratsun-etal2005,roussel2006validation}, we might model transcription as the simple delayed mass-action \cite{Roussel96b} process
\begin{equation}
\ce{pro($t$) + RNAPII($t$)
->[$k_\mathrm{ini}$] pro($t$) + RNAPII($t+\tau_\mathrm{tr}$) + RNA($t+\tau_\mathrm{tr}$)},
\label{DStranscription}\end{equation}
where $\mathrm{pro}(t)$ is understood to be a promoter that is able to initiate transcription, i.e.,\ one that is in the `on' state.
The activity of a gene is regulated by transcription factors acting either as repressors or as activators. These interactions would be modeled by separate mass-action processes, e.g.,\
\begin{subequations}
\begin{gather}
\ce{pro + R <=> pro.R},\\
\ce{R + L <=> R.L}
\end{gather}
\end{subequations}
for a repressor R that is inactivated by binding to a ligand~L. Thus, regulatory processes that turn a gene on and off would be modeled by their own reactions and would not be involved in determining the distribution of the delay~$\tau_\mathrm{tr}$, which is the focus of this contribution. In other words, we study the `pure' transcriptional process on the assumption that binding or conformational transition events that turn a gene on or off are handled separately.

Having narrowed our focus as noted above, the transcription model is divided into the usual three phases:  initiation, elongation, and termination. This section is correspondingly organized.

\subsection{Initiation}\label{sec:init}
The initiation phase starts with binding of initiation factors, including the TATA-binding protein (TBP),  to the TATA box of the promoter region (pro) of a gene [reaction (\ref{bind})], following which RNAPII can bind to the promoter, forming the pre-initiation complex (PIC)   [reaction~(\ref{PIC})] \cite{Buratowski-etal1989}:
\begin{equation}\label{bind}
\ce{TBP + pro  ->[$k_\mathrm{bind}^{(2)}$] TBP.pro},
\end{equation}
\begin{equation}\label{PIC} 
\ce{TBP.pro + RNAPII ->[$k_\mathrm{PIC}^{(2)}$] PIC}.
\end{equation}
The rate constants $k_\mathrm{bind}^{(2)}$ and $k_\mathrm{PIC}^{(2)}$ are second-order rate constants; the superscripts appear for emphasis. The reverse of reaction~\eqref{bind} is neglected because TBP dissociation from DNA is extremely slow, with a half-life of several minutes~\cite{CP1995,DJS+2005}. TBP probably dissociates rapidly following transcription initiation~\cite{SKM+2008}, but we need not consider the recycling of this factor given that we focus here on a single polymerase. Dissociation of RNAPII from the promoter is, on the other hand, a frequent event~\cite{Darzacq-etal2007}, so $k_\mathrm{PIC}^{(2)}$ is best thought of as an effective rate constant that takes into account the repeated binding and dissociation of the polymerase prior to forming a stable PIC, at least for now. We consider the effect of reversible binding of the polymerase explicitly in Section~\ref{sec:revRNAPbinding}.

After forming a stable PIC, the  RNA polymerase moves to the first nucleotide of the DNA template. Reaction (\ref{init}) represents the translocation of the active site of the RNA polymerase to the first  nucleotide (nucleotide 1) if this nucleotide is unoccupied (U state), thus converting $\mathrm{U_{1}}$ into an occupied  nucleotide (O state):
\begin{equation}\label{init}
\ce{PIC + U_{1} ->[\ce{$k_\mathrm{init}$}] O_{1}}.
\end{equation}
Since we only consider a single polymerase here, the forward site will always be unoccupied. We maintain this notation for consistency with our earlier work \cite{roussel2013distribution,roussel2006stochastic,Vashishtha}. Note that the engagement of the polymerase represented by reaction~\eqref{init} is likely a complex process. However, reaction~\eqref{init} will be a reasonable approximation if this process has a single rate-limiting step.

The minimal initiation model considered here captures the slowest steps in a single round of initiation at promoters that include a TATA box. Initiation models can be much more complicated~\cite{BRE2011}. We could, for example, allow for reinitiation by having either the promoter or \ce{TBP.pro} become available again after promoter clearance (i.e.,\ after the polymerase has traveled sufficiently far downstream).\footnote{\emph{In vivo}, \ce{TBP.pro} is rapidly removed after initiation by Mot1~\cite{SKM+2008}. Including a reaction to represent the action of Mot1 following clearance of the promoter would likely lead to transcriptional bursts, a topic to which we will return in the discussion.} This is easily done in stochastic simulations~\cite{Vashishtha}. However, here we want to focus on analytic distributions obtained for a single polymerase. Some analytic progress can be made for a model with multiple polymerases~\cite{Hosseini2016}. Some of these results however depend on the single-polymerase distributions, so the work presented here is an important starting point. Moreover, not all genes display transcriptional bursting~\cite{LHE2021}, so it is worth considering  initiation of a single polymerase.

From here and until we reach the termination phase, transcription consists of cycling through binding and recognition of a nucleoside triphosphate complementary to the template nucleotide, translocation to the next template position, and formation of a phosphodiester bond, which adds the selected nucleotide to the growing RNA chain. The early steps are slower than the corresponding steps in productive elongation~\cite{HBG+2006} but, for simplicity, we ignore these kinetic differences in this contribution. Thus, the binding of the first nucleotide is treated as indistinguishable from binding of any other nucleotide and is included in the elongation part of the model.

\subsection{Elongation}
Although conceptually there are at least three distinct processes in elongation (binding of a complementary nucleoside triphosphate, translocation, catalysis of phosphodiester bond formation), the simplest models of transcription describe the process in terms of two states \cite{roussel2006stochastic,woo2006analytical}, roughly corresponding to the pre- and post-translocated states. More \cite{bai2007mechanochemical,bai2004sequence,LZK+2012} or fewer \cite{greive2011development} states can be used to model elongation.
In our models \cite{roussel2013distribution,roussel2006stochastic,Vashishtha}, the two states involved in productive elongation are called occupied ($\mathrm{O}$) and activated~($\mathrm{A}$) (Figure~\ref{fig:elongation_states}), corresponding respectively to the pre-  and post-translocated states.

 \begin{figure}[H]\centering
  \includegraphics[width=0.8\textwidth]{elongation_scheme.pdf}
 \caption{Schematic of the elongation phase showing the occupied (O), activated (A) and paused (P) states at three consecutive nucleotides.}
\label{fig:elongation_states}
\end{figure}
 
The process of activation and translocation of RNA polymerase repeats itself as RNAPII translocates nucleotide by nucleotide along the gene.  These repeated reactions are responsible for productive elongation [reactions (\ref{ae}) and (\ref{add})]. 
\begin{eqnarray}\label{ae}
&\ce{O_{$i$} ->[$k_{\mathrm{ae}}$] A_{$i$}},
&i = 1,2,\ldots,N;\\
\label{add}
&\ce{A_{$i$} + U_{$i+1$} ->[$k_{\mathrm{add}}$] O_{$i+1$} + U_{$i$}}, 
 &i = 1,2,\ldots,N-1.
\end{eqnarray}
Here, $N$ is the length of the template, i.e.,\ the transcribed sequence.

During elongation, the polymerase may pause for any of a number of reasons: There are so-called ``ubiquitous'' pauses that may occur in any part of the gene \cite{CW2011}, promoter-proximal pauses associated with a change in state of the polymerase as it transitions from initiation to elongation \cite{CC2011}, pauses associated with splicing \cite{BS2015}, and pauses associated with termination \cite{GWP2006}. Pausing models can be more or less complex and incorporate backtracking \cite{voliotis2008fluctuations} or arrest \cite{IKL+2013}. Here we opt for a simple model of a reversible pause. We neglect backtracking, which is likely appropriate for short, ubiquitous pauses, but may be less so for long pauses. It appears that pausing proceeds from the pre-translocated (O) state \cite{IKL+2013,Roberts2014}, so we write
\begin{equation}
\ce{O_{$i$} <=>[$k_{p,i}$][$k_{r,i}$] P_{$i$}},\qquad i = 1,2,\ldots,N.
\end{equation}
We allow here for the possibility that the rate constants for pausing, $k_{p,i}$, and for the release from the pause, $k_{r,i}$, might be template-dependent. On the other hand, we need make no such allowance for the rate constants associated with productive elongation, $k_\mathrm{ae}$ and $k_\mathrm{add}$, which are believed to be constant over the length of a gene and across genes within a particular organism \cite{SP2009}. Figure~\ref{fig:elongation_states} shows the off-pathway P state along with the O and A states.

\subsection{Termination}
According to the allosteric
model of termination, transcription of the poly(A) signal at the end of the gene can cause a termination-inducing conformational change  
in the elongation complex \cite{rosonina2006terminating}. After the last reaction of elongation [activation at nucleotide $N$, reaction (\ref{ae})], the elongation complex is converted into a termination complex (TC) [reaction (\ref{tc})]. Finally, the termination complex dissociates and releases the RNA [reaction (\ref{term})].  
\begin{equation}\label{tc}
\ce{A_{$N$} ->[$k_{\mathrm{TC}}$] TC},
\end{equation}
\begin{equation}\label{term}
\ce{TC ->[$k_{\mathrm{term}}$] RNA + RNAPII + U_{$N$}}.
\end{equation}

The above reactions assume 100\% efficiency of allosteric termination. In reality, there is at least one, and possibly two failsafe mechanisms backing up allosteric termination \cite{LB2015}. We assume a single mechanism here for simplicity. In terms of the completion of the transcript, the torpedo model will behave similarly to allosteric termination since its first step is to cleave the transcript at the poly(A) signaling sequence, releasing the former from the transcriptional machinery \cite{LB2015}.

\section[Probability distributions for a single RNA Polymerase]{Time-dependent probability distribution for transcription by RNAPII}
\sectionmark{Probability distribution for a single RNA Polymerase}
\label{sec:kineticmodel}

\subsection{Initiation}
Figure~\ref{figsolplot} shows a schematic of the model. RNA polymerases move one nucleotide at a time along the template strand~\cite{abbondanzieri2005direct}, so each of the numbered sites in the diagram represents one nucleotide. The polymerase initiates at the promoter, represented by a single site labeled~`0' in the diagram, and terminates at the end of the transcribed sequence, represented here as the $N$'th nucleotide.  By solving the chemical
master equation (CME) for each state in the three phases of transcription, we can obtain time-dependent probability distributions for the progress of transcription, from which we can recover the  transcriptional delay distribution.

\begin{figure}[H]\centering
  \includegraphics[width=0.8\textwidth]{rnap2.pdf}
 \caption{A schematic representation of RNA polymerase on the DNA template strand.}
\label{figsolplot}
\end{figure}

In principle, the propensities of reactions (\ref{bind}) and (\ref{PIC}) depend on the populations of available TBP and RNAPII, respectively. Here, we assume that there are constant  pools of these two species in the nucleus so that we can treat reactions (\ref{bind}) and~(\ref{PIC}) as pseudo-first order reactions, with rate constants $k_\mathrm{bind}=k_\mathrm{bind}^{(2)}[\mathrm{TBP}]$, and $k_\mathrm{PIC}=k_\mathrm{PIC}^{(2)}[\mathrm{RNAPII}]$, respectively.
The governing master equations  are
\begin{subequations}
    \begin{align}\label{MA-pro}
\frac{dP_{\mathrm{pro}} (t)}{dt}&=-k_{\mathrm{bind}} P_{\mathrm{pro}}(t),\\
\label{MA-TBP.pro}
\frac{dP_{\mathrm{TBP \cdot pro}} (t)}{dt}&=k_{\mathrm{bind}} P_{\mathrm{pro}} (t)-k_{\mathrm{PIC}} P_{\mathrm{TBP \cdot pro}}(t),\\
\label{MA-PIC}
\frac{dP_{\mathrm{PIC}} (t)}{dt}&=k_{\mathrm{PIC}} P_{\mathrm{TBP \cdot pro}} (t)-k_{\mathrm{init}} P_{\mathrm{PIC}}(t),
\end{align}
\end{subequations}
with initial condition $P_{\mathrm{pro}}(0)=1$, $P_{\mathrm{TBP \cdot pro}}(0)= P_{\mathrm{PIC}}(0)=0$. In these equations, $P_{\mathrm{pro}}(t)$ is the probability that the promoter has never been bound up to time~$t$, i.e.,\ it is the survival probability of the empty promoter. Since we treat the single-polymerase case, we do not need to consider promoter clearance~\cite{Hsu2002}. $P_{\mathrm{TBP\cdot pro}}(t)$ is the probability that, at time~$t$, TBP is bound to the promoter. Similarly, $P_{\mathrm{PIC}}(t)$ is the probability that a complete PIC is assembled at the promoter at time~$t$. The initial condition corresponds to a gene that has been turned on at $t=0$, prior to which TBP would not have been able to bind the promoter. It is possible to consider other initial conditions. For example, for some modes of regulation, it might be possible for TBP to bind the promoter while the polymerase is unable to bind until the gene is turned on. In this case, we would choose an initial condition such that $P_{\mathrm{pro}}(0)+P_{\mathrm{TBP \cdot pro}}(0)=1$ and $P_{\mathrm{pro}}(0)/P_{\mathrm{TBP \cdot pro}}(0)=K_d^{\mathrm{(eff)}}$, where $K_d^{\mathrm{(eff)}}$ is the effective dissociation constant for the $\mathrm{TBP \cdot pro}$ complex at cellular TBP concentrations.

The time-dependent rate at which probability exits  initiation  and enters  elongation  is, from reaction~(\ref{init}),
\begin{equation}
v_\mathrm{init}(t) \equiv \rho_\mathrm{init}(t) = k_\mathrm{init}P_\mathrm{PIC}(t),
\label{rhoinit_def}\end{equation}
where $\rho_\mathrm{init}(t)$ is the distribution of initiation times. The equality of $v_\mathrm{init}$ and $\rho_\mathrm{init}$ follows because the integral of the rate at which initiation is completed is the cumulative distribution of initiation times.

The linear ordinary differential equations (\ref{MA-pro}) to (\ref{MA-PIC}) can be solved to find the probability distribution of being in each state. Equation~\eqref{rhoinit_def} can then be used to obtain the distribution of initiation times. The latter is (assuming all of the rate constants are distinct)
\begin{equation}\label{ppic2}\small
 \rho_\mathrm{init}(t) = \frac{k_\mathrm{bind}\,k_\mathrm{PIC}\,k_\mathrm{init}}{(k_\mathrm{PIC}
-k_\mathrm{bind})  (k_\mathrm{bind}-k_\mathrm{init})  (k_\mathrm{init}-k_\mathrm{PIC}) }
   \bigg[   (k_\mathrm{PIC}-k_\mathrm{init}) e^{-k_\mathrm{bind}t}+ (k_\mathrm{init} -k_\mathrm{bind})e^{-k_\mathrm{PIC}t}
+ (k_\mathrm{bind}-k_\mathrm{PIC})e^{-k_\mathrm{init}t}\bigg].
\end{equation}
Figure \ref{rho3} shows $\rho_\mathrm{init}(t)$. The parameters used in this study were estimated by Vashishtha \cite{Vashishtha} and are given in Table~\ref{tab:param}. They are likely representative of genes expressed at lower rates, although Vashishtha showed that they were consistent with the average behavior of  transcription in mammalian cells~\cite{KSC2002}. Briefly, the rate of TBP binding to the TATA box is assumed not to be limited by reaction~(\ref{bind}), but by dissociation of TBP dimers, so the effective rate constant for this non-elementary step is the rate constant for dimer dissociation, $k_\mathrm{bind}=0.0016\,\mathrm{s^{-1}}$ \cite{CP1997}. The rate constant for PIC assembly, $k_\mathrm{PIC}=0.0029\,\mathrm{s^{-1}}$, was taken from an \emph{in vitro} study \cite{FKG2001}. Finally, based on experimental data for the progress of the polymerase through the first few nucleotides~\cite{FKG2001}, Vashishtha inferred a lower bound for $k_\mathrm{init}$ of $0.3\,\mathrm{s^{-1}}$ \cite{Vashishtha}. The value of $k_\mathrm{init}$ eventually adopted is somewhat arbitrary, but in any event, it is clear that binding of TBP and PIC assembly are the rate limiting steps. Although the parameters chosen for this study may represent a gene at which initiation is particularly inefficient, typical eukaryotic initiation times are relatively long \cite{FKG2001,Larson-etal2011,SJG+2018}, resulting in well-spaced polymerases \cite{JIM1998}, which supports our decision to focus on a single-polymerase model.
\begin{table}
    \centering
    \caption{Default values of the rate constants used in this study.}
    \begin{tabular}{lD{.}{.}{4}}
    \hline
    \textbf{Rate constant}   & \multicolumn{1}{c}{\textbf{Value/$\mathrm{s^{-1}}$}}\\
    \hline
        $k_\mathrm{bind}$ & 0.0016 \\
        $k_\mathrm{PIC}$ & 0.0029 \\
         $k_\mathrm{init}$ & 0.6 \\
    $k_\mathrm{ae}$ & 144 \\
    $k_\mathrm{add}$ & 144 \\
    $k_{p,i}$ & 100 \\
    $k_{r,i}$ & 3.6 \\
    $k_\mathrm{TC}$ & 0.0032 \\
    $k_\mathrm{term}$ & 0.003\\
    \hline
    \end{tabular}
    \label{tab:param}
\end{table}
\begin{figure}\centering
  \includegraphics{rhoinit}
 \caption{$\rho_\mathrm{init}(t)$ from Eq~\eqref{ppic2} with the parameters of Table~\ref{tab:param}.}
\label{rho3}
\end{figure}

\subsection{Elongation} \label{secelongation}
\subsubsection{Continuum approximation}\label{continuumsection}
The distribution of transcription times is an $N$-fold convolution of passage times through each nucleotide. For very simple models, this convolution is tractable, and is the basis of the well-known ``linear chain trick'' \cite{HK2019,MacDonald1976,MacDonald1977a,MacDonald89}. This problem is less straightforward if we take into account deviations from the sequential addition of nucleotides such as pausing. For typical protein-coding genes in eukaryotes, $N$ might range from several hundred nucleotides to  tens of thousands \cite{Venter-etal2001}. The large value of $N$ for typical genes and the high rates of elongation~\cite{Darzacq-etal2007} afford an opportunity to make a continuum approximation. The methods applied here are somewhat analogous to those used to obtain a wave equation in a chain of coupled oscillators~\cite{Goldstein-wave-eq} and have a long history in mathematical biology \cite{Othmer1983,RR2004}. Note also the close analogy to the derivation of delays in age-structured models of populations~\cite{SVN89}.

Denote by $P_{i,\sigma}$ the probability that the polymerase has its active site at nucleotide $i$ and that it is in state $\sigma$, where $\sigma$ can take on the values O, P or A. The master equations governing the time evolution of the elongation phase are 
\begin{subequations}\label{CMEonesite}\begin{align}
\frac{dP_{i,\mathrm{O}}(t)}{dt}&=k_{\mathrm{add}} P_{i-1,\mathrm{A}}(t)-k_{\mathrm{ae}} P_{i,\mathrm{O}}(t)
-k_{p,i}P_{i,\mathrm{O}}(t) + k_{r,i}P_{i,\mathrm{P}}(t),\label{MA-iO}\\
\frac{dP_{i,\mathrm{A}}(t)}{dt}&=k_{\mathrm{ae}} P_{i,\mathrm{O}}(t)-k_{\mathrm{add}} P_{i,\mathrm{A}}(t),\label{MA-iA}\\
\frac{dP_{i,\mathrm{P}}(t)}{dt} &=k_{p,i}P_{i,\mathrm{O}}(t) - k_{r,i}P_{i,\mathrm{P}}(t),\label{MA-iP}
\end{align}\end{subequations}
where $i=2,3,...,N-1$. (Equation~\eqref{MA-iO} is slightly modified for $i=1$ since the O state is reached from the PIC state according to reaction~\eqref{init}. This does not alter the argument below in any substantial way.) 

Because the states are mutually exclusive, the probability that the polymerase is at nucleotide $i$ in any state, $P_i$, is the sum of the probabilities $P_{i,\sigma}$, that is
\begin{equation}\label{p-i}
P_{i}(t) = P_{i,\mathrm{O}}(t) + P_{i,\mathrm{A}}(t) + P_{i,\mathrm{P}}(t).
\end{equation}
Using Eqs\  (\ref{CMEonesite})--(\ref{p-i}), we obtain
\begin{equation}\label{MA-i}
\frac{dP_{i}}{dt}=k_\mathrm{add} \left(P_{i-1,\mathrm{A}}- P_{i,\mathrm{A}}\right).
\end{equation}

We want to treat the motion on the time scale of elongation. We therefore introduce the rescaled time
\begin{equation}
\theta = t/\tau_e,
\label{thetadef}
\end{equation}
where $\tau_e$ is the mean total elongation time. This quantity will be computed below. Introducing this scaling in Eq~(\ref{MA-i}), we get
\begin{equation}
\frac{dP_i}{d\theta} = k_\mathrm{add}\tau_e\left(P_{i-1,\mathrm{A}}- P_{i,\mathrm{A}}\right).
\label{dPidtheta}
\end{equation}

For long genes, an increment of one nucleotide is a small change in position measured relative to the length of the gene. Therefore define
\begin{align}
x &= i/N,\\
\delta &= 1/N.\label{deltadef}
\end{align}
In a continuous representation, the $P_{i,\mathrm{A}}(\theta)$, which are supported on the natural numbers, will become continuous functions of $x$, $\rho_\mathrm{A}(x,\theta)$, with the correspondence  $P_{i,\mathrm{A}}(\theta)=\delta\rho_\mathrm{A}(i\delta,\theta)$. Similarly, we define the probability density for all states combined by $P_i(\theta)=\delta\rho_e(i\delta,\theta)$. Because $\delta$ is small, from the definition of the derivative, we have
\begin{equation}
\begin{split}
P_{i-1,\mathrm{A}}(\theta)- P_{i,\mathrm{A}}(\theta)&=\delta\left(\rho_A(x-\delta,\theta)-\rho_A(x,\theta)\right)\\
&\approx -\delta^2\frac{\partial \rho_\mathrm{A}(x,\theta)}{\partial x}.
\end{split}
\label{dPAdx}\end{equation}
This will hold provided the $P_{i,\mathrm{A}}$ vary sufficiently slowly. Otherwise, the estimate of the derivative on which Eq~(\ref{dPAdx}) is based is meaningless. This approach will therefore fail if the gene contains a specific site at which probability  accumulates for a time during elongation due to significant slowing of elongation over a short stretch of the sequence. Take the extreme case where the residence time at one particular nucleotide, $i$, is several orders of magnitude larger than at adjacent nucleotides: then there will be a period of time of $O(\tau_i)$ over which $P_i(t)\sim O(1)$ with $P_j(t)\approx 0$ at adjacent nucleotides and where the derivative estimate of Eq~\eqref{dPAdx} is therefore invalid. Slowing could be due to any of the elongation parameters, but will most commonly   occur because of a strong pause site. We treat the latter in Section~\ref{sec:pause_in_one_nucleotide}. On the other hand, slowly varying parameters, including those associated with pausing, do not in principle pose any particular difficulty for this approach. We note in passing that the treatment of translation by Mier-y-Ter{\'a}n-Romero and coworkers also required sufficiently slowly varying parameters. 

Substituting (\ref{dPAdx}) and the definitions of the densities into (\ref{dPidtheta}), we get the PDE
\begin{equation}
\frac{\partial \rho_e}{\partial\theta} =-k_\mathrm{add}\tau_e\delta\frac{\partial \rho_\mathrm{A}}{\partial x}.\label{dPidtheta_vs_gradA}
\end{equation}
Equation (\ref{dPidtheta_vs_gradA}) is not closed as it relates the time derivative of $\rho_e$ to the gradient of $\rho_A$. If we examine the state of a specified nucleotide from an ensemble of similarly prepared systems at an arbitrary moment in time,\footnote{Because the $\tau_i$ are short, Eq~\eqref{PiA_related_to_Pi} should hold true provided observations are not made at the earliest time that some molecules in the ensemble reach nucleotide~$i$. In other words, \eqref{PiA_related_to_Pi} should hold except in a brief window of duration $O(\tau_i)$ when $P_i(t)$ is beginning to rise.} the conditional probability of catching it in a particular state given that the polymerase has reached this nucleotide is proportional to its mean dwell time in that state. In particular,
\begin{equation}
P_{i,\mathrm{A}} = P_iP\left(\mathrm{A}|i\right) = P_i\frac{\tau_{i,\mathrm{A}}}{\tau_i},
\label{PiA_related_to_Pi}\end{equation}
where $\tau_{i,\mathrm{A}}=k_\mathrm{add}^{-1}$ is the mean dwell time in the $\mathrm{A}$ state at nucleotide~$i$, and $\tau_i$ is the total dwell time at nucleotide~$i$. Converting to densities, Eq~(\ref{dPidtheta_vs_gradA}) becomes
\begin{equation}
\frac{\partial \rho_e}{\partial\theta} =-\delta\tau_e\frac{\partial\phantom{x}}{\partial x}\left(\frac{\rho_e}{\tau(x)}\right),\label{elongPDE}
\end{equation}
where $\tau(i\delta)=\tau_i$, and $\tau(x)$ interpolates between these values for $x\ne i\delta$. A similar equation has been obtained by Mier-y-Ter{\'a}n-Romero and coworkers for translation by the ribosome~\cite{TSH2010}. Following  these authors, provided the velocity is sufficiently slowly varying, we obtain a mean elongation delay of, converting time back to dimensional units,
\begin{equation}
    \tau_e = N\int_0^1\tau(x)\,dx.
\end{equation}
We see from this equation that $\tau(x)$ has the interpretation of the time required for the polymerase to move the distance of one nucleotide ($\delta$~in the $x$~scaling), which is consistent with the original description of $\tau(x)$ as interpolating the $\tau_i$.

We will focus for now on the case in which the kinetic parameters are identical at every nucleotide. Then, $\tau(x)=\tau_i$ is a constant, and 
\begin{equation}
\tau_e=(N-1)\tau_i+k_\mathrm{ae}^{-1}\approx N\tau_i
\end{equation}
for large $N$.
Taking into account~(\ref{deltadef}), (\ref{elongPDE}) reduces to
\begin{equation}
\frac{\partial \rho_e}{\partial\theta} = -\frac{\partial \rho_e}{\partial x},
\label{elongPDEsimple}\end{equation}
i.e.,\ a simple advective equation with unit velocity.
(An alternative derivation of this equation is to be found in~\cite{Hosseini2016}.)
At the left boundary ($x=0$), probability flows in at a rate $k_\mathrm{init}P_\mathrm{PIC}(t)$ [reaction~(\ref{init})] or, in the scaling~(\ref{thetadef}), $\tau_ek_\mathrm{init}P_\mathrm{PIC}(\tau_e\theta)$. The appropriate boundary condition for~(\ref{elongPDEsimple}) is therefore
\begin{equation}
\rho_e(0,\theta) = \tau_ek_\mathrm{init}P_\mathrm{PIC}(\tau_e\theta).
\label{bcrhoe}\end{equation}
(Technically, since the nucleotides are numbered from 1, the boundary condition should be applied at $x=\delta$ rather than $x=0$. However, since $\delta$ is assumed small, there is no significant error in the boundary condition given above.) Since we deal with just one polymerase, and since this polymerase cannot start transcription before $t=0$, we also have the initial condition $\rho_e(x,0)=0$.

The PDE~(\ref{elongPDEsimple}) with the boundary condition~(\ref{bcrhoe}) simply transports probability through the elongation compartment at unit velocity. Accordingly, the solution is
\begin{equation}
\rho_e(x,\theta) = \tau_ek_\mathrm{init}P_\mathrm{PIC}\left(\tau_e(\theta-x)\right)H(\theta-x),
\end{equation}
where $H(\cdot)$ is the Heaviside function. The flux exiting  elongation, at $x=1$, transformed back to the original time scale, is therefore
\begin{equation}
J_e(t) = k_\mathrm{init}P_\mathrm{PIC}(t-\tau_e)H(t-\tau_e).
\label{elongationflux}\end{equation}
This solution holds for any fixed kinetic parameters, with or without pausing. We conclude that for transcription, as was the case for translation \cite{TSH2010}, elongation normally acts as a simple delay. A similar result was obtained elsewhere using a Laplace transform technique  \cite{roussel2013distribution}.

Figure \ref{position} shows the position of the polymerase against time obtained from a stochastic simulation for a single polymerase in the absence of pausing, while Figure~\ref{position_w_pauses} shows representative trajectories from similar simulations  with frequent short pauses. These simulations confirm that  for long genes without long pause sites, the elongation process is essentially advective, with  RNAPII moving at a roughly constant  velocity. From an ensemble of 40\,000 realizations, we found a mean and standard deviation of the elongation rate of $72.0\pm 1.8\,\mathrm{nt\,s^{-1}}$. The situation with frequent short pauses is not so different. Over 40\,000 realizations, we found that the velocities during the elongation phase varied from 3.7 to $6.5\,\mathrm{nt\,s^{-1}}$, with a mean of $4.8\,\mathrm{nt\,s^{-1}}$. While the observed range of elongation rates is relatively large, the variability in the elongation times is tiny compared to the variability in the initiation and termination times, as is evident from Figure~\ref{position_w_pauses}. Interestingly, the relatively constant time required to get through elongation has been observed experimentally \cite{Larson-etal2011}. Bel and coworkers have also provided theoretical arguments demonstrating that the distribution of a process consisting of sequential steps, even if the individual steps have some complexity, would be expected to converge to a narrow distribution when the number of steps is large~\cite{BMN2010}.
\begin{figure} \centering
 \includegraphics[width=0.43\textwidth]{p_t.pdf}
 \caption[Position of the RNA polymerase vs time]{Position of the RNA polymerase vs. time obtained from the stochastic
simulation of a single polymerase in the absence of pausing.  Parameter values are as in Table~\ref{tab:param} except $k_{p,i}=k_{r,i}=0$ with  $N=1000$.
\label{position}}
\end{figure}
\begin{figure} \centering
 \includegraphics[scale=0.6]{fig5}
 \caption{Several realizations of the RNA polymerase position vs. time obtained from  stochastic
simulations of single polymerases with ubiquitous pausing.   The parameter values are  as in Table~\ref{tab:param} with $N=1000$.
\label{position_w_pauses}}
\end{figure}

\subsubsection{The elongation delay}
Equation (\ref{elongationflux}) shows that, subject to the underlying assumption that elongation is a fast process, elongation only contributes a delay. Mier-y-Ter{\'a}n-Romero and coworkers have shown that the same holds true for slowly varying delays~\cite{TSH2010}. In this model, which leaves out backtracking that would require hydrolysis of the 3' end of the transcript, addition of nucleotides is strictly sequential. Accordingly,
\begin{equation}
\tau_e = \sum_{i=1}^{N-1}\tau_i + k_\mathrm{ae}^{-1}.
\label{taue_def}\end{equation}
(Note that reaction~(\ref{ae}) includes $N$ steps, but reaction~(\ref{add}) only includes $N-1$ steps, where $N$ is the length of the transcribed sequence.)
Similarly, the transition from the pre- to the post-translocated state is irreversible, so
\begin{equation}
\tau_i = \tau_{i,\mathrm{O}\vee\mathrm{P}} + \tau_{i,\mathrm{A}},
\label{tauisum}\end{equation}
where $\tau_{i,\mathrm{O}\vee\mathrm{P}}$ is the mean dwell time in the aggregated pre-translocation O and P states. As noted previously, $\tau_{i,\mathrm{A}}=k_\mathrm{add}^{-1}$. We now turn to the calculation of $\tau_{i,\mathrm{O}\vee\mathrm{P}}$.

We focus on the steps
\begin{equation}
\ce{P_{$i$} <=>[$k_{r,i}$][$k_{p,i}$] O_{$i$} ->[$k_\mathrm{ae}$]}\end{equation}
which define the escape process from the $\mathrm{O}\vee\mathrm{P}$ aggregated state. The master equation for this escape process is
\begin{subequations}
\begin{align}
\frac{dP_{i,\mathrm{P}}}{dt} &= k_{p,i}P_{i,\mathrm{O}} - k_{r,i}P_{i,\mathrm{P}},\\
\frac{dP_{i,\mathrm{O}}}{dt} &= -k_{p,i}P_{i,\mathrm{O}} + k_{r,i}P_{i,\mathrm{P}} - k_\mathrm{ae}P_{i,\mathrm{O}},
\end{align}\end{subequations}
with initial conditions $P_{i,\mathrm{O}}(0) = 1$, $P_{i,\mathrm{P}}(0)=0$. The probability density for the dwell time in the aggregated O and P states at nucleotide $i$, i.e.,\ for making a transition from state O to state A, is $\rho_{\mathrm{OA},i}=k_\mathrm{ae}P_{i,\mathrm{O}}(t)$. After some tedious algebra, we find
\begin{subequations}
    \begin{equation}
\rho_{\mathrm{OA},i}(t) = \frac{k_\mathrm{ae}}{2\sqrt{D_i}}\Bigl[e^{-(S_i+\sqrt{D_i})t/2}\left(\sqrt{D_i}+k_\mathrm{ae}+k_{p,i}-k_{r,i}\right)
+e^{-(S_i-\sqrt{D_i})t/2}\left(\sqrt{D_i}-k_\mathrm{ae}-k_{p,i}+k_{r,i}\right)\Bigr],
\label{rhoOAi}\end{equation}
where
\begin{align}
S_i &= k_\mathrm{ae}+k_{p,i}+k_{r,i},\\
D_i &= S_i^2-4k_\mathrm{ae}k_{r,i}.
\end{align}\end{subequations}

As usual, we compute the mean time spent in the O and P states by
\begin{equation}
\tau_{i,\mathrm{O}\vee\mathrm{P}} = \int_0^\infty t\rho_{\mathrm{OA},i}dt = \frac{k_{p,i}+k_{r,i}}{k_\mathrm{ae}\,k_{r,i}}.
\label{tauiOorP}\end{equation}
Equation (\ref{tauisum}) now implies the dwell time at a particular nucleotide:
\begin{equation}
\tau_i = \frac{k_{p,i}+k_{r,i}}{k_\mathrm{ae}k_{r,i}} + \frac{1}{k_\mathrm{add}} = \frac{1}{k_\mathrm{ae}}\left(1+\frac{k_{p,i}}{k_{r,i}}\right) + \frac{1}{k_\mathrm{add}}.
\label{tauival}\end{equation}

Note that Eqs~(\ref{rhoOAi}) and~(\ref{tauival}) were derived without any approximations. Thus, they apply to any nucleotide, including those where long pauses occur, discussed in Section~\ref{sec:pause_in_one_nucleotide}.

It is instructive to consider a few special limits of Eq~(\ref{tauival}). First, if there is no pausing, $k_{p,i}=0$, and we get
\begin{equation}
\left.\tau_i\right|_{k_{p,i}=0}=\frac{1}{k_\mathrm{ae}}+\frac{1}{k_\mathrm{add}}\end{equation}
which is, as expected, the sum of the mean dwell times in the O and A states. On the other hand, for very strong pause sites with $k_{p,i}/k_{r,i}\gg \max\left(1,k_\mathrm{ae}/k_\mathrm{add}\right)$, we get
\begin{equation}
\tau_i \approx \frac{k_{p,i}}{k_{r,i}}\frac{1}{k_\mathrm{ae}}.\end{equation}
In this case, the dwell time in the O and P states dominates, and is computed from the pause-free O-state dwell time stretched by the proportion of the time spent in the paused state.

$\tau_i$ is the mean time required to get from nucleotide $i$ to nucleotide $i+1$. If the kinetic parameters are identical at each nucleotide, then $\tau_i^{-1}$ is  the mean velocity of the polymerase. For the parameters of Figure~\ref{position}, we can calculate a velocity $v=72\,\mathrm{nt\,s^{-1}}$, which is in perfect agreement with the mean velocity obtained from the numerical experiments mentioned earlier. For the case with pausing (Figure~\ref{position_w_pauses}), Eq~(\ref{tauival}) leads to a velocity of $v=4.8\,\mathrm{nt\,s^{-1}}$, which again agrees with the mean over many realizations of the stochastic simulations.

It will also be useful to consider the variance in the total elongation time $\tau_e$.  According to the central limit theorem \cite{Billingsley1995}, assuming $N$ is large, the distribution of $\tau_e$ approaches a Gaussian with variance
\begin{equation}
    \var(\tau_e) = \sum_{i=1}^N\var(\tau_i),
\label{eq:vartaue}\end{equation}
making the obvious large-$N$ approximations. The $\tau_i$ need not be identically distributed, but there should not be a small number of nucleotides where passage is particularly slow. This case will be considered in Section~\ref{sec:pause_in_one_nucleotide}.
The variance of $\tau_i$ is easily calculated:
 \begin{equation}
     \var(\tau_i) = \var(\tau_{i,\mathrm{O}\vee\mathrm{P}}) + \var(\tau_{i,\mathrm{A}})
 \end{equation}
and $\var(\tau_{i,\mathrm{A}})=k_\mathrm{add}^{-2}$. To calculate $\var(\tau_{i,\mathrm{O}\vee\mathrm{P}})$, we first compute $\mathrm{E}(\tau_{i,\mathrm{O}\vee\mathrm{P}}^2)$, where $\mathrm{E}(\cdot)$ denotes the expectation value:
\begin{equation}
\mathrm{E}(\tau_{i,\mathrm{O}\vee\mathrm{P}}^2) = \int_0^\infty t^2\rho_{\mathrm{OA},i}dt = \frac{2(k_{p,i}+k_{r,i})^2 + 2k_\mathrm{ae}k_{p,i}}{k_\mathrm{ae}^2k_{r,i}^2}.
\end{equation}
Then,
\begin{equation}
     \var(\tau_{i,\mathrm{O}\vee\mathrm{P}}) = \mathrm{E}(\tau_{i,\mathrm{O}\vee\mathrm{P}}^2)-\tau_{i,\mathrm{O}\vee\mathrm{P}}^2 = \frac{(k_{p,i}+k_{r,i})^2 + 2k_\mathrm{ae}k_{p,i}}{k_\mathrm{ae}^2k_{r,i}^2}
\end{equation}
so that
\begin{equation}
     \var(\tau_i)=\frac{(k_{p,i}+k_{r,i})^2 + 2k_\mathrm{ae}k_{p,i}}{k_\mathrm{ae}^2k_{r,i}^2}+\frac{1}{k_\mathrm{add}^2}.
\end{equation}

In the special case where all sites are identical,
the distribution of $\tau_e$ therefore converges to a Gaussian with mean given by Eqs\ \eqref{taue_def} and~\eqref{tauival} and variance
\begin{equation}
    \var(\tau_e) = N\var(\tau_i) = N\left(\frac{(k_{p,i}+k_{r,i})^2 + 2k_\mathrm{ae}k_{p,i}}{k_\mathrm{ae}^2k_{r,i}^2}+\frac{1}{k_\mathrm{add}^2}\right).
\end{equation}
If we let $\sigma^2_e\equiv\var(\tau_e)$, then the distribution of elongation times can be written
\begin{equation}
    \rho_e(t) = \frac{e^{-(t-\tau_e)^2/2\sigma^2_e}}{\sqrt{2\pi\sigma^2_e}}.
\label{eq:Gaussian_rhoe}\end{equation}
Several comments about this equation are in order:
\begin{enumerate}
    \item Provided we compute $\tau_e$ from Eq~\eqref{taue_def} and $\sigma^2_e$ from Eq~\eqref{eq:vartaue}, this equation is applicable to any situation where there are no sites where transit is unusually slow, i.e.,\ no long pausing sites. It is not necessary for the distributions of the dwell sites to be identical at each nucleotide.
    \item In cases like those of Figures\ \ref{position} and~\ref{position_w_pauses} where the distribution of elongation times is narrow relative to the distribution of at least one of elongation or termination, $\rho_e(t)$ can be replaced by a Dirac $\delta$~distribution. Then elongation contributes a simple fixed delay to the overall distribution of transcription times. This is the situation to which Eq~(\ref{elongationflux}) applies.
    \item Technically, the expression for $\rho_e(t)$ is supported on $t\in(-\infty,\infty)$. Clearly, the negative values of $t$ are meaningless. Thus, Eq~\eqref{eq:Gaussian_rhoe} will only apply when $\tau_e\gg\sigma_e$ such that $\rho_e(0)$ is extremely small and the $t<0$ tail of the distribution contributes negligibly to the total probability. In these cases, we can treat $\rho_e(t)$ as if it is  supported on the physically meaningful interval $t\in[0,\infty)$ with negligible error. Since $\sigma_e$ grows as $\sqrt{N}$ and $\tau_e$ grows as $N$, this will always be the case for sufficiently large $N$.
\end{enumerate}

% Finally, the coefficient of variation for the residence time at a single nucleotide, which measures the relative width of a distribution, is given by the ratio of the standard deviation to the mean:
% \begin{equation}
%     \mathrm{CV}(\tau_i)=\frac{\sqrt{\var(\tau_i)}}{\tau_i} = \frac{\sqrt{k_\mathrm{add}^2\left[(k_{p,i}+k_{r,i})^2 + 2k_\mathrm{ae}k_{p,i}\right] + k_\mathrm{ae}^2k_{r,i}^2}}{k_\mathrm{add}(k_{p,i}+k_{r,i})+k_\mathrm{ae}k_{r,i}}.
% \end{equation}

% For the parameters of Fig.~\ref{position}, the coefficient of variation of the  dwell time works out to 0.7. Fluctuations away from the mean dwell time are relatively small so that a constant velocity is an excellent approximation in the elongation phase. For the parameters of Fig.~\ref{position_w_pauses}, the CV is 1.9. The velocity is still roughly constant, but occasional longer pauses are visible and the resulting stochasticity accounts for the difference between the calculated mean velocity and the velocity calculated by regression from the stochastic simulation.

\subsection{Termination}
Assume that we are in the normal situation where the flux out of  elongation  and into  termination  is given by Eq\ (\ref{elongationflux}).
Thus, the following ODEs govern the evolution during  termination:
\begin{subequations}\label{termME}
    \begin{align}
\frac{dP_{N,\mathrm{A}}(t)}{dt}&=J_e(t) - k_\mathrm{TC} P_{N,\mathrm{A}}(t),\\
\frac{dP_\mathrm{TC}(t)}{dt}&=k_\mathrm{TC} P_{N,\mathrm{A}}(t)-k_\mathrm{term} P_\mathrm{TC}(t),\\
\frac{dP_\mathrm{T}(t)}{dt}&=k_\mathrm{term} P_\mathrm{TC}(t),
\end{align}\end{subequations}
where $P_\mathrm{TC}(t)$ is the probability that there is a termination complex (TC) at time $t$, and  $P_\mathrm{T}(t)$ represents the probability that the transcript (T) has been  released, with $J_e(t)$ given by Eq~(\ref{elongationflux}). 

\subsection{Distributions of transcription times}
We can solve the ODEs~\eqref{termME} with the initial conditions of $P_{N,\mathrm{A}}(0)= P_\mathrm{TC}(0)=P_\mathrm{T}(0)=0$ to derive an expression for probability density of the entire transcription process, $\rho_{\mathrm{tr}}(t)$ ($=dP_\mathrm{T}/dt=k_\mathrm{term} P_\mathrm{TC}(t)$):\footnote{There is at least one other way to arrive at Eq~\eqref{rhoT1}: Below, we will derive the distribution of termination times, Eq~\eqref{rho-T}. We can compute the convolution of the distributions of initiation times, Eq~\eqref{ppic2}, with the distributions of termination times and with a $\delta$~distribution representing the distribution of elongation times.}
\begin{equation}\label{rhoT1}\begin{split}
 \rho_{\mathrm{tr}}(t) =
H(t-\tau_{e}) \, &k_\mathrm{bind}\,k_\mathrm{PIC}\,k_\mathrm{init}\,k_\mathrm{TC}\,k_\mathrm{term}\\
\times\Bigl\{&\bigl[\left( k_{\mathrm{bind}}-k_{\mathrm{PIC}} \right)\left( k_{\mathrm{bind}}-k_{\mathrm{init}} \right)  \left(k_{\mathrm{bind}}- k_{\mathrm{TC}} \right)  \left( k_{\mathrm{bind}}-k_{\mathrm{term}} \right)   \bigr]^{-1}\ e ^{-k_{\mathrm{bind}} \, (t-\tau_{e}) }\\
&+\bigl[\left( k_{\mathrm{PIC}}-k_{\mathrm{bind}} \right) \left( k_{\mathrm{PIC}}-k_{\mathrm{init}} \right)  \left( k_{\mathrm{PIC}}- k_{\mathrm{TC}} \right)  \left( k_{\mathrm{PIC}}-k_{\mathrm{term}} \right)   \bigr]^{-1}\ e ^ {-k_{\mathrm{PIC}} \,  (t-\tau_{e})}\\
&  +\bigl[ \left(k_{\mathrm{init}}- k_{\mathrm{bind}} \right) \left(k_{\mathrm{init}}- k_{\mathrm{PIC}}  \right) \left(k_{\mathrm{init}}- k_{\mathrm{TC}} \right) \left(k_{\mathrm{init}}- k_{\mathrm{term}}  \right)   \bigr]^{-1}\ e ^{-k_{\mathrm{init}}\, (t-\tau_{e})}\\
 &+\bigl[\left( k_{\mathrm{TC}}-k_{\mathrm{bind}} \right) \left( k_{\mathrm{TC}}-k_{\mathrm{PIC}} \right)  \left( k_{\mathrm{TC}}-k_{\mathrm{init}} \right)     \left( k_{\mathrm{TC}}-k_{\mathrm{term}} \right)\bigr]^{-1}\ e ^{-k_{\mathrm{TC}}\,  (t-\tau_{e}) }\\
 &+\bigl[\left(k_{\mathrm{term}}- k_{\mathrm{bind}} \right) \left(k_{\mathrm{term}}- k_{\mathrm{PIC}} \right)\left(k_{\mathrm{term}}- k_{\mathrm{init}} \right)\left(k_{\mathrm{term}}- k_{\mathrm{TC}} \right)      \bigr]^{-1} \ e ^ {-k_{\mathrm{term}}\, (t-\tau_{e}) }\Bigr\}.
 \end{split}\end{equation}

Figures \ref{compare2} and \ref{ubiquitous_fig} show comparisons between the analytical solution and stochastic simulations of $\rho_{\mathrm{tr}}(t)$, showing excellent agreement. Note the effect of ubiquitous pausing in Figure~\ref{ubiquitous_fig}: the shape of the distribution is not changed, but the distribution is shifted to the right with increasing pausing probability. In other words, ubiquitous pauses only contribute to the constant elongation delay under normal conditions.
   \begin{figure}\centering
  \includegraphics[width=0.45\textwidth]{transcription-density2.pdf}
 \caption[A comparison between stochastic simulation  and analytical solutions for $\rho_T(t)$]{A comparison between stochastic simulation  (circles) and analytical (red line) solutions for $\rho_{\mathrm{tr}}(t)$ in the absence of pausing. Parameter values are as in Table~\ref{tab:param} except $k_{p,i}=k_{r,i}=0$ with  $N=1000$. The numerical distribution was computed from 40\,000 realizations of the transcription process.}
\label{compare2}
\end{figure}
 \begin{figure}\centering
  \includegraphics[width=0.45\textwidth]{transcription-density-ubiq5.pdf}
 \caption{Probability density for transcription time with ubiquitous pausing. We compare stochastic simulation results with analytical results. Parameter values are as in Table~\ref{tab:param} with  $N=1000$ except for the $k_{p,i}$ values. The pausing rate constant, $k_{p,i}$, is  $7.6\,\mathrm{s}^{-1}$ (red curve with triangles), $100\,\mathrm{s}^{-1}$ (blue curve with circles) and $500\,\mathrm{s}^{-1}$ (black curve with squares). The numerical distribution was computed from 40\,000 realizations of the transcription process.}
\label{ubiquitous_fig}
\end{figure}

An alternative approach to termination, which will shortly prove useful, is  to solve directly for the distribution of termination times. We do this by solving the equations
\begin{subequations}
    \begin{align}
\frac{dP_{N,\mathrm{A}}(t)}{dt}&=-k_\mathrm{TC} P_{N,\mathrm{A}}(t),\\
\frac{dP_\mathrm{TC}(t)}{dt}&=k_\mathrm{TC} P_{N,\mathrm{A}}(t)-k_\mathrm{term} P_\mathrm{TC}(t),
\end{align}\end{subequations}
subject to the initial conditions $P_{N,\mathrm{A}}(0)=1$, $P_\mathrm{TC}(0)=0$.
Some elementary mathematics shows the distribution of termination times to be
\begin{equation}\label{rho-T}
\rho_\mathrm{T}(t)=k_\mathrm{term}P_\mathrm{TC}=
\frac{k_{\mathrm{term}}\,k_{\mathrm{TC}}}{k_{\mathrm{TC}}-k_{\mathrm{
term}}} \,\, \left( { e ^{-k_{\mathrm{term}}t}} -{ e ^{-k_{\mathrm{TC}}
t}} \right).
\end{equation}

Equation~\eqref{rhoT1} will provide an accurate approximation to  the  distribution of transcription times provided elongation is must faster than initiation and termination combined. Otherwise, we must consider the Gaussian distribution of elongation times. There are various ways in which we can imagine globally slowing elongation. Numerical experiments reveal that they all have the same effect since the distribution of elongation times is Gaussian, regardless of the details of the kinetics at any particular nucleotide, again with the proviso that there there are not narrow regions of a few nucleotides where transit is unusually slow. Consider in particular a situation in which the supply of nucleotides is strongly limited. In this model, this corresponds to decreasing the value of~$k_\mathrm{ae}$.
Because of the role of nucleotides in the Brownian ratchet mechanism of RNA polymerases \cite{BarNahum-etal2005,LZK+2012}, low occupancy of the nucleotide addition site would be expected to cause increased backtracking and pausing as well \cite{IKL+2013,Shaevitz-etal2003}, further reducing the polymerase's average elongation rate.

We can derive an equation for the distribution of transcription times of a slow polymerase by taking a convolution of the initiation, termination and elongation time distributions under the assumption that Eq~\eqref{eq:Gaussian_rhoe} accurately describes the distribution of elongation times. This is most conveniently done by first taking the convolution of the initiation and termination time distributions, which gives an equation similar to Eq~\eqref{rhoT1} but with $\tau_e=0$. We then take the convolution of this initiation-termination distribution, $\rho_\mathrm{IT}$, with the elongation time distribution as follows:
\begin{equation}
    \rho_\mathrm{tr}(t) = \int_{\tau=0}^t\rho_\mathrm{IT}(t-\tau)\rho_e(\tau)\,d\tau.
\end{equation}
The result of this calculation is
\begin{equation}\begin{split}
\rho_\mathrm{tr}(t) = \frac{1}{2}k_\mathrm{bind}&\,k_\mathrm{PIC}\,k_\mathrm{init}\,k_\mathrm{TC}\,k_\mathrm{term}\\
\times\verybig\{&\vcenter{\hbox{$\displaystyle\frac{e^{-k_\mathrm{bind}(t-\tau_e-k_\mathrm{bind}\sigma^2_e/2)}\left[\erf\left(\frac{k_\mathrm{bind}\sigma^2_e+\tau_e}{\sqrt{2\sigma^2_e}}\right)-\erf\left(\frac{k_\mathrm{bind}\sigma^2_e-t+\tau_e}{\sqrt{2\sigma^2_e}}\right)\right]}{(k_\mathrm{bind}-k_\mathrm{PIC})(k_\mathrm{bind}-k_\mathrm{init})(k_\mathrm{bind}-k_\mathrm{TC})(k_\mathrm{bind}-k_\mathrm{term})}$}}\\
&+\frac{e^{-k_\mathrm{PIC}(t-\tau_e-k_\mathrm{PIC}\sigma^2_e/2)}
\left[\erf\left(\frac{k_\mathrm{PIC}\sigma^2_e+\tau_e}{\sqrt{2\sigma^2_e}}\right)-\erf\left(\frac{k_\mathrm{PIC}\sigma^2_e-t+\tau_e}{\sqrt{2\sigma^2_e}}\right)\right]}{(k_\mathrm{PIC}-k_\mathrm{bind})(k_\mathrm{PIC}-k_\mathrm{init})(k_\mathrm{PIC}-k_\mathrm{TC})(k_\mathrm{PIC}-k_\mathrm{term})}\\
&+\frac{e^{-k_\mathrm{init}(t-\tau_e-k_\mathrm{init}\sigma^2_e/2)}\left[\erf\left(\frac{k_\mathrm{init}\sigma^2_e+\tau_e}{\sqrt{2\sigma^2_e}}\right)-\erf\left(\frac{k_\mathrm{init}\sigma^2_e-t+\tau_e}{\sqrt{2\sigma^2_e}}\right)\right]}{(k_\mathrm{init}-k_\mathrm{bind})(k_\mathrm{init}-k_\mathrm{PIC})(k_\mathrm{init}-k_\mathrm{TC})(k_\mathrm{init}-k_\mathrm{term})}\\
&+\frac{e^{-k_\mathrm{TC}(t-\tau_e-k_\mathrm{TC}\sigma^2_e/2)}\left[\erf\left(\frac{k_\mathrm{TC}\sigma^2_e+\tau_e}{\sqrt{2\sigma^2_e}}\right)-\erf\left(\frac{k_\mathrm{TC}\sigma^2_e-t+\tau_e}{\sqrt{2\sigma^2_e}}\right)\right]}{(k_\mathrm{TC}-k_\mathrm{bind})(k_\mathrm{TC}-k_\mathrm{PIC})(k_\mathrm{TC}-k_\mathrm{init})(k_\mathrm{TC}-k_\mathrm{term})}\\
&+\vcenter{\hbox{$\displaystyle\frac{e^{-k_\mathrm{term}(t-\tau_e-k_\mathrm{term}\sigma^2_e/2)}\left[\erf\left(\frac{k_\mathrm{term}\sigma^2_e+\tau_e}{\sqrt{2\sigma^2_e}}\right)-\erf\left(\frac{k_\mathrm{term}\sigma^2_e-t+\tau_e}{\sqrt{2\sigma^2_e}}\right)\right]}{(k_\mathrm{term}-k_\mathrm{bind})(k_\mathrm{term}-k_\mathrm{PIC})(k_\mathrm{term}-k_\mathrm{init})(k_\mathrm{term}-k_\mathrm{TC})}$}}\verybig\}
\end{split}\label{rhoquasiexact}\end{equation}

We now consider a polymerase whose rate constant for the O to A transition (i.e.,\ binding of the appropriate nucleotide) has been reduced by a factor of 100 while the pausing rate constant has been set to $10\,\mathrm{s}^{-1}$. The results of stochastic simulations are plotted along with  Eqs\ \eqref{rhoT1} and~\eqref{rhoquasiexact} in Figure~\ref{fig:slow_polym}(a). Because, at these high relative rates of pausing, much of the elongation time is spent in pauses, this reduction in $k_\mathrm{ae}$ and $k_{p,i}$ from the default values has the effect of increasing $\tau_i$  by a factor of  more than~10, from 0.21 to $2.63\,\mathrm{s}$. Even at this extremely slow rate of translation (about $0.4\,\mathrm{nt\,s^{-1}}$), the advective approximation to elongation still holds and either expression for the distribution of elongation times fits the simulation data equally well. It is only at even lower elongation rates that the more complex Eq~\eqref{rhoquasiexact} is necessary. For example, Figure~\ref{fig:slow_polym}(b) shows the distribution at $k_{p,i}=200\,\mathrm{s}^{-1}$ with the same value of $k_\mathrm{ae}$ as in panel~(a). At these parameters, $\tau_i$ has increased to 39\,s, corresponding to a glacial rate of transcription. The distribution of transcription times is now strongly affected by the Gaussian distribution of elongation times. Equation~\eqref{rhoT1} no longer applies, and Eq~\eqref{rhoquasiexact}  now gives the correct distribution.
\begin{figure}
    \centering
    \includegraphics[width=0.5\linewidth]{slow_polym_mod_pausing}\includegraphics[width=0.5\linewidth]{slow_polym_lots_pausing}
    \caption{Transcription time distributions for polymerases that are transcribing slowly. In both panels, $k_\mathrm{ae}=1.44\,\mathrm{s}^{-1}$ and $N=1000$. In panel~(a), $k_{p,i}=10\,\mathrm{s^{-1}}$ while in panel~(b), $k_{p,i}=200\,\mathrm{s^{-1}}$. All other parameters are as in Table~\ref{tab:param}. Circles show simulation results, the dotted blue lines show the transcription time distribution in the advective approximation, Eq~\eqref{rhoT1}, while the solid red line is the quasi-exact Eq~\eqref{rhoquasiexact}.
    \label{fig:slow_polym}}
\end{figure}


While Figure~\ref{fig:slow_polym} was produced at parameters that generate unrealistically slow elongation rates, the same issue would arise if both initiation and termination were fast processes, as they will sometimes be. Figure~\ref{fig:fast_init_and_term} shows the simulated distribution along with the approximations \eqref{rhoT1} and~\eqref{rhoquasiexact} for such a case. We see again that the advective approximation is a poor fit while Eq~\eqref{rhoquasiexact} provides an excellent approximation to the distribution of transcription times. The reason is exactly the same as it was for the parameters of Figure~\ref{fig:slow_polym}(b): when initiation or termination are slow compared to the duration of the elongation process, then Eq~\eqref{rhoT1} applies; otherwise, we must consider the distribution of elongation times as well as the distributions of initiation and termination times, and Eq~\eqref{rhoquasiexact} must be used, as is the case in Figures \ref{fig:slow_polym}(b) and~\ref{fig:fast_init_and_term}.
\begin{figure}
    \centering
    \includegraphics[scale=0.7]{fast_init_and_term}
    \caption{Transcription time distributions for the case where both
    initiation and termination are fast processes, with parameters
    $k_\mathrm{bind}=1$, $k_\mathrm{PIC}=0.2$, $k_\mathrm{init}=0.6$, $k_\mathrm{ae}=144$, $k_\mathrm{add}=144$, $k_{p,i}=20$, $k_{r,i}=3.6$, $k_\mathrm{TC}=0.3$ and $k_\mathrm{term}=0.4\,\mathrm{s}^{-1}$ for a 10\,000\,nt gene. Circles show the distribution computed from 40\,000 stochastic simulations, the dotted blue line is Eq~\eqref{rhoT1}, and the solid red line is Eq~\eqref{rhoquasiexact}.
    \label{fig:fast_init_and_term}}
\end{figure}

We emphasize that the simulation data are always well fit by Eq~\eqref{rhoquasiexact}, while Eq~\eqref{rhoT1} is only valid if elongation is fast compared to initiation and termination. Thus, Eq~\eqref{rhoquasiexact} is more general. The only difficulty with using Eq~\eqref{rhoquasiexact} is that its evaluation generally requires higher precision than standard double-precision floating-point arithmetic.
We used Maple to evaluate this expression, setting the number of digits of precision to values up to 35 (adjusted upward when anomalies were noted), roughly corresponding to quadruple-precision floating-point numbers.

\section{Effect of a long pause along a gene}
\label{sec:pause_in_one_nucleotide}
The case of a site at which a lengthy pause may occur has to be treated separately, for the reasons outlined in Section~\ref{continuumsection}. We take advantage here of the modular nature \cite{greive2011development} of our model. Specifically, we can break up the elongation region into segments where only short ubiquitous pauses occur, dealt with using the theory of Section \ref{secelongation}, and segments where long pause sites are found, as illustrated in Figure~\ref{longpause}. At any given long pause site, the distribution of the dwell time in the combined O and P states is given by Eq~(\ref{rhoOAi}). We can take a convolution of this dwell time distribution with the other relevant distributions to obtain the  overall transcriptional delay distribution. For simplicity, we will focus in this section on the case where elongation is advective, although it should equally be possible to carry out these calculations for the case of elongation times with a Gaussian distribution.
\begin{figure}
\centering\includegraphics[scale=0.7]{long_pause_regions}
\caption{Schematic illustration of a long pause region, here consisting of two nucleotides. Outside the region where long pauses may occur, elongation is treated as contributing a delay. In the region of the long pause, the dwell time distribution is explicitly computed.
\label{longpause}}
\end{figure}

The convolution of probability densities for initiation and termination times ($\rho_\mathrm{I,T}$) is calculated by
\begin{equation}\label{conv1}
\rho_\mathrm{I,T} (\tau)= \int_{0}^\tau \rho_\mathrm{PIC}(\tau-t)\,\rho_\mathrm{T}(t)\,dt.
\end{equation}
This equation gives us an alternative path to Eq~(\ref{rhoT1}), by simply shifting the distribution $\rho_\mathrm{I,T}(t)$ by the elongation delay,~$\tau_e$.

For a long pause at a single nucleotide at $i=i_p$, the case on which we focus here,  we can obtain the convolution of $\rho_\mathrm{I,T}(t)$ and $\rho_{OA,i_p}(t)$, Eqs\ \eqref{conv1} and~\eqref{rhoOAi}, respectively:
\begin{equation}
\rho_\mathrm{I,T,P} (\tau)= \int_{0}^\tau \rho_{OA,i_p}(\tau-t)\,\rho_\mathrm{I,T}(t)\,dt.
\end{equation}
In the event that long pauses can occur at multiple nucleotides, this last step is iterated, i.e.,\ we take the $q$-fold convolution with $\rho_{OA,i_p}$, where $q$ is the number of nucleotides at which long pauses may occur and $i_p$ represents each of the positions where a long pause occurs, potentially with different kinetic parameters at each of these sites.
Once the convolution has been calculated, the final result can be obtained by shifting the distribution by $\tau_e'$, the  delay due to all processes in elongation not involving the $\mathrm{O}\vee\mathrm{P}$ aggregated state at the long pause site. The distribution for the transcription time is therefore
\begin{equation}
\rho_\mathrm{tr}(t) = H(t-\tau_e')\,\rho_{I,T,P}(t-\tau_e').
\end{equation}
The expression that results from this calculation is awkward and is not reproduced here, but is available in a Maple worksheet provided in the Supplemental Materials.

Let $\mathcal{P}$ be the set of nucleotides at which a long pause may occur,
and $\mathcal{R}$ be the set of nucleotides where ``regular'' elongation
occurs. $\mathcal{R}\cup\mathcal{P}$ is the entire transcribed sequence,
i.e.,\ $\mathcal{R}$ is the complement of $\mathcal{P}$ in the transcribed
sequence. The set $\mathcal{P}$ consists of $q$~nucleotides.
$\tau_e'$ includes $N-1$ extension steps [reaction~(\ref{add})], and $N-q$ steps with short pauses for which (\ref{tauiOorP}) applies. Thus, $\tau_e'$ is, in general,
\begin{equation}
\tau_e' = \frac{N-1}{k_\mathrm{add}} + \sum_{i\in\mathcal{R}}\frac{k_{p,i}+k_{r,i}}{k_\mathrm{ae}\,k_{r,i}}.
\label{taueprime}\end{equation}
For the specific case where long pauses only occur at one nucleotide, this reduces to $\tau_e' = (N-1)\tau_i$, with $\tau_i$ given by Eq~(\ref{tauival}).

Figure~\ref{strong} shows the probability density of the transcription time for three different values of pausing rate constant in a case with a very long pause (mean pause time of $500\,\mathrm{s}$). Pauses of this duration might be associated with splicing, which is quite a slow process \cite{SP2009}. 

 \begin{figure}\centering
  \includegraphics[width=0.45\textwidth]{strong-kpause}
 \caption{Effect of a strong pause in the gene on the probability density of transcription comparing stochastic simulation results (symbols) and the analytical solution (lines). Parameter values are as in Table~\ref{tab:param} with  $N=1000$ except that $k_{p,i}=k_{r,i}=0$ everywhere except at nucleotide $i_p$. At this nucleotide, $k_{r,i_p}=0.002\,\mathrm{s}^{-1}$ with each curve representing a  different value of $k_{p,i_p}$:  $k_{p,i_p}=1000\,\mathrm{s}^{-1}$ for the black line and squares, $500\,\mathrm{s}^{-1}$ for the blue line and circles, and $200\,\mathrm{s}^{-1}$ for the red line and triangles.  Each numerical distribution was computed from $10^5$ realizations of the transcription process.}
\label{strong}
\end{figure}

\section{Reversible RNAPII binding}\label{sec:revRNAPbinding}
It is possible to imagine many variations on the model presented in section~\ref{sec:Eukaryoticmodel}. Because RNA polymerase binding is likely reversible on time scales relevant to transcription initiation~\cite{Darzacq-etal2007}, we consider a model in which reaction~\eqref{PIC} is replaced by
\begin{equation} 
\ce{TBP.pro + RNAPII <=>[$k_\mathrm{PIC}^{(2)}$][$k_{-\mathrm{PIC}}$] PIC}.
\end{equation}
The master equation for initiation becomes
\begin{subequations}
\begin{align}
    \frac{dP_\mathrm{pro}}{dt} &= -k_\mathrm{bind}P_\mathrm{pro},\\
    \frac{dP_\mathrm{TBP\cdot pro}}{dt} &= k_\mathrm{bind}P_\mathrm{pro} - k_\mathrm{PIC} P_\mathrm{TBP\cdot pro} + k_{-\mathrm{PIC}}P_\mathrm{PIC},\\
    \frac{dP_\mathrm{PIC}}{dt} &= k_\mathrm{PIC} P_\mathrm{TBP\cdot pro} - k_{-\mathrm{PIC}}P_\mathrm{PIC} - k_\mathrm{init}P_\mathrm{PIC}
\end{align}
\end{subequations}
in the same pseudo-first-order regime as was previously assumed to treat RNAPII binding to the promoter.
Also as before, we use the initial condition $P_{\mathrm{pro}}(0)=1$, $P_{\mathrm{TBP \cdot pro}}(0)= P_{\mathrm{PIC}}(0)=0$. Define
\begin{subequations}
    \begin{align}
        s &= \sqrt{k_\mathrm{init}^2 + 2k_\mathrm{init}\left(k_{-\mathrm{PIC}}-k_\mathrm{PIC}\right) + \left(k_\mathrm{PIC}+k_{-\mathrm{PIC}}\right)^2},\\
        \lambda_1 &= \frac{1}{2}\left(k_\mathrm{PIC}+k_{-\mathrm{PIC}}+k_\mathrm{init}+s\right),\\
        \lambda_2 &= \frac{1}{2}\left(k_\mathrm{PIC}+k_{-\mathrm{PIC}}+k_\mathrm{init}-s\right).
    \end{align}
\end{subequations}
The resulting distribution of initiation times can then be written
\begin{equation}
    \rho_\mathrm{init}(t) = \frac{k_\mathrm{bind}\,k_\mathrm{PIC}\,k_\mathrm{init}}{s(\lambda_1-k_\mathrm{bind})(\lambda_2-k_\mathrm{bind})}\left[e^{-\lambda_2t}(k_\mathrm{bind}-\lambda_1) - e^{-\lambda_1t}(k_\mathrm{bind}-\lambda_2) + se^{-k_\mathrm{bind}t}\right].
\label{rhoinit_w_rev_RNAP_binding}\end{equation}
It is easy to verify that~\eqref{rhoinit_w_rev_RNAP_binding} reduces to~\eqref{ppic2} when $k_{-\mathrm{PIC}}=0$.

For the parameters used in our previous figures, the distribution of initiation times is not very sensitive to the value of $k_{-\mathrm{PIC}}$. Figure~\ref{fig:vary_kmPIC} shows the distribution of initiation times at different values of $k_{-\mathrm{PIC}}$ for a set of parameters inspired by the data of Darzacq and coworkers~\cite{Darzacq-etal2007}. The value of $k_\mathrm{init}$ suggested by their data is much smaller than the one used by Vashishtha~\cite{Vashishtha}. As a result, dissociation of the polymerase and entry into the elongation phase compete directly, which explains why $\rho_\mathrm{init}$ is much more sensitive to $k_{-\mathrm{PIC}}$ in this regime.
\begin{figure}
    \centering
    \includegraphics[width=0.7\textwidth]{vary_kmPIC}
    \caption{Distribution of initiation times for $k_\mathrm{bind}=1.6\times 10^{-3}\,\mathrm{s^{-1}}$, $k_\mathrm{PIC}=0.022\,\mathrm{s^{-1}}$, $k_\mathrm{init}=1.59\times 10^{-3}\,\mathrm{s^{-1}}$, and $k_{-\mathrm{PIC}}=0$ (red solid line), 0.03 (blue dotted line) or $0.15\,\mathrm{s^{-1}}$ (black dashed line).
    \label{fig:vary_kmPIC}}
\end{figure}

In Section~\ref{sec:init}, we suggested that we could treat binding of RNA polymerase as an irreversible process with an effective rate constant $k_\mathrm{PIC}$, i.e.,\ that we could use a reduced value of this rate constant instead of taking reversible binding into account. To test this idea, we calculated the average initiation time in the model with reversible RNAPII binding, and then chose $k_\mathrm{PIC}$ in the irreversible model to match the average initiation time. We carried out this calculation for the parameters of Figure~\ref{fig:vary_kmPIC} where the distribution is more sensitive to the value of $k_{-\mathrm{PIC}}$. The result of this calculation is shown in Figure~\ref{fig:deflate}. Differences in the shapes of the two distributions can clearly be seen. These differences persist to surprisingly small values of  $k_{-\mathrm{PIC}}$, but eventually vanish (results not shown), as they must. Thus, there appears to be a wide range of parameters where reversibility has a noted effect on the distribution of initiation times. Whether these differences are sufficient to be distinguishable based on noisy experimental data is an open question.
\begin{figure}
    \centering
    \includegraphics[width=0.5\textwidth]{deflate}
    \caption{Distribution of initiation times for the model with reversible RNAPII binding (dashed black curve) and for the model with irreversible RNAPII binding (solid red curve). The value of $k_\mathrm{PIC}$ was adjusted in the latter model to match the average initiation time of the former. For the model with reversible binding (dashed black curve),  $k_\mathrm{PIC}=0.022\,\mathrm{s^{-1}}$ and $k_{-\mathrm{PIC}}=0.15\,\mathrm{s^{-1}}$. For the model with irreversible binding (solid red curve), $k_\mathrm{PIC}=2.31\times 10^{-4}\,\mathrm{s^{-1}}$.
    Common parameters: $k_\mathrm{bind}=1.6\times 10^{-3}\,\mathrm{s^{-1}}$, $k_\mathrm{init}=1.59\times 10^{-3}\,\mathrm{s^{-1}}$.
    \label{fig:deflate}}
\end{figure}

\section{Conclusions}\label{sec:Discussion}
In  this work, a simplified version of the transcription model  proposed by Vashishtha \cite{Vashishtha} for eukaryotic cells was utilized to describe the kinetics of a single RNAPII during transcription. The master equations for  initiation and termination  were solved exactly. As in the work of Mier{-}y{-}Ter{\'a}n{-}Romero et al.~\cite{TSH2010}, for the case where elongation is fast, the motion during elongation  led to an advection equation so that the elongation phase has the effect of a simple delay connecting  initiation and termination.
There are other ways of showing that simple elongation will often be adequately  modeled by a fixed delay \cite{roussel2013distribution}. Ultimately, all of these methods depend, directly or indirectly, on the linear growth of both the mean elongation time and the variance with the number of nucleotides such that the width of the distribution becomes relatively narrower as $N$ grows~\cite{BMN2010}. Contrastingly, in the case where elongation is relatively slow, we found that the distribution of transcription times was poorly reproduced unless we took into account the Gaussian distribution of elongation times. It then becomes an interesting question whether the resulting distribution of transcription times is dynamically important or if, given the relatively narrow and roughly symmetric distribution typically obtained in these cases, a simple fixed delay would yield the same phenomenology in most gene expression models. These questions can be explored using delay-stochastic simulations \cite{GB2000,Bratsun-etal2005}, where it is easy to vary the delay distribution~\cite{roussel2006validation} and thus to determine the sensitivity of the dynamics to these distributions. The few studies available so far suggest that the delay distribution is sometimes important whether the dynamics are treated stochastically~\cite{BCH2013} or deterministically \cite{RW2007,JP2015}, but sometimes not~\cite{FDT2016}.

In addition to delay distributions for homogeneous or near-homogeneous
elongation processes, we considered the possibility of long pause sites,
which are extreme cases of systems with variable transcription rate. Elsewhere, one of us (SHH) found that
 variable velocity along the template  can lead to some non-trivial phenomenology~\cite{Hosseini2016}.
Interestingly, elongation in regions without long pauses just contributes a fixed delay according to Eq~\eqref{taueprime}, regardless of where the pause occurs relative to the rapid-elongation regions. This result is reminiscent of Epstein's identification of ``bottleneck intermediates'' in his treatment of delays arising from  consecutive first-order reactions~\cite{Epstein90}.

The delay distributions presented here are likely to be particularly useful for delay-stochastic simulations \cite{Bratsun-etal2005,LGR2012,StochPy,RL2007a,roussel2006validation}. The delay-stochastic simulation algorithm (DSSA) is an extension of the Gillespie algorithm \cite{gillespie1976general} that allows for delayed product formation. Thus, transcription can, in the simplest case, be written as the single delayed reaction~\eqref{DStranscription}.
The transcriptional delay may be fixed or  distributed. In the latter case, a new value of the delay is generated from the appropriate distribution each time this reaction ``fires'' \cite{roussel2006validation}. Using the DSSA therefore requires that we know the distribution of transcription times, $\rho_\mathrm{tr}(t)$, if only to determine whether or not a fixed delay is appropriate. The distributions derived here can be compared to facilitate these decisions once the parameters of a model have been selected. When Eqs\ \eqref{rhoT1} and~\eqref{rhoquasiexact} agree, elongation contributes a fixed delay; otherwise, it will likely be necessary to consider a distribution of elongation delays.

While we have focused on the total transcription time, including the events
of initiation, there are a few other ways our results could be used. We
could explicitly model the initiation steps, and use a (distributed) delay
to represent only elongation and termination. This would be the simplest
approach to modeling the expression of a gene whose transcription is
regulated at initiation, and again, modularity plays a role in allowing us
to use parts of the calculations presented here in more complex situations.
For a gene regulated at a pause site \cite{sims2004elongation}, elongation
could be broken up into pieces, very much as we did to obtain the
distributions shown in Figure~\ref{strong}, treating regulation at the pause
site explicitly. For example, if a release factor (RF) such as P-TEFb
\cite{Jennings2013} is required to release a paused complex (PC) from an
obligatory pause, we could write, neglecting promoter clearance,\footnote{If
we wanted to deal with promoter clearance, we would replace
$\mathrm{pro}(t)$ on the right-hand side of Eq~\eqref{formPC} by $\mathrm{pro}(t+\tau_c)$, where $\tau_c$ is the time required for RNAPII to clear the promoter and therefore make the latter available for binding another molecule of polymerase~\cite{TR2020}.}
\begin{subequations}\label{model_w_PC}
\begin{gather}
\ce{pro($t$) + RNAPII($t$) ->[$k_\mathrm{ini}$] pro($t$) + PC($t+\tau_p$)},\label{formPC}\\
\ce{PC($t$) + RF($t$) ->[$k_\mathrm{rel}$]  RNAPII($t+\tau_f$) + RNA($t+\tau_f$) + RF($t+\tau_\mathrm{RFC}$)}\label{release_from_pause}
\end{gather}
\end{subequations}
where $\tau_p$ is the time required to reach the pause site, $\tau_f$ is the time required to finish transcription after release from the pause and $\tau_\mathrm{RFC}$ is the lifetime of the polymerase--release-factor complex.
If there are no other sites in the gene where long pauses occur and if the sequences preceding and following the pause site are sufficiently long, $\tau_p$ and $\tau_f$ could be treated as fixed delays as seen in Section~\ref{secelongation}. The interesting dynamics would then lie in the production and/or destruction of RF, which would itself be described by appropriate reactions, including possibly delayed reactions representing its transcription and translation.

In the case of a gene expressed in bursts, the results obtained here could be used to describe the transcription time distribution of the first polymerase in a burst, on the assumption that interactions between polymerases do not affect the behavior of the pioneer polymerase.\footnote{It is likely that polymerases do affect each other's motions \cite{epshtein2003cooperation,epshtein2003transcription,JKL2014}. However, the nature of these interactions is, to our knowledge, not well characterized at this time, which makes them difficult to incorporate into models. Treating the first polymerase in a burst as traveling freely along the gene may be the best we can do right now, give or take the possibility of considering the prevention of backtracking by trailing polymerases. Moreover, the approximation that the pioneer polymerase travels freely is likely to run into trouble in the ``Kurtz'' limit \cite{JWY2019} where frequent bursts mean that the leading polymerase in a burst may well interact with the trailing polymerase of the previous burst. Thus, Eq~\eqref{burst_react} is likely more appropriate in Jia et al.'s ``L{\'e}vy'' limit where bursts are well separated in time.} We could write a delayed mass-action reaction for a burst as follows:
\begin{equation}
    \mathrm{pro}(t) \xrightarrow{k_\mathrm{ini}} \mathrm{pro}_i(t + \tau_b) + \mathrm{RNA}(t+\tau_\mathrm{tr}) + \sum_{i=2}^b \mathrm{RNA}(t+\tau_i),
\label{burst_react}\end{equation}
where $\mathrm{pro}_i$ is an inactive form of the promoter. This reaction assumes that the nuclear concentration of RNAPII is roughly constant in time so that this concentration can be included in the pseudo-first-order rate constant $k_\mathrm{ini}$.
The delay $\tau_b$ is the duration of a burst generating $b$ RNA molecules.
This delay is required to ensure that a  burst does not start at a
particular copy of a gene before the previous burst has ended. The delay
$\tau_\mathrm{tr}$ is the transcription time of the pioneer polymerase
which, as we argued above, could be randomly generated according to a distribution from a single-polymerase model. The $\tau_i$ are the completion times of the subsequent transcripts measured from the beginning of the burst. These completion times could be computed sequentially as follows:
\begin{equation}
    \tau_i = \tau_\mathrm{tr} + \sum_{j=1}^{i-1}\Delta\tau_i,
\end{equation}
where $\Delta\tau_i=\tau_{i+1}-\tau_i$ (identifying
$\tau_1\equiv\tau_\mathrm{tr}$) is the time between completions of
subsequent transcripts. Although we have not tried this, the distribution of
this quantity should be derivable from our model in the low-traffic limit by
considering the time required for promoter clearance, initiation of a
subsequent transcript (possibly considering retention of TBP for a time
after the pioneer polymerase has initiated transcription) and differences in
termination times between consecutive polymerases. Alternatively,
Szavits-Nossan and Grima have recently obtained some results on the
distribution of $\Delta\tau_i$ \cite{SG2023}. Both the distributions of
burst duration $\tau_b$ and of burst size $b$  could be generated from a
model \cite{BKW+2013,DB2009}, or experimental distributions could be used
\cite{CFX2006,LJH2019,Suter-etal2011,WAR2021}.\footnote{Note however that
the value of $\tau_b$ is implied by the values of the random variables $b$ and $\Delta\tau_i$, so experimental distributions of the burst time could be used to validate the models for the distributions of $b$ and/or $\Delta\tau_i$.} The model would be completed by including one or more reactions that return the promoter to the active state. For telegrapher's noise, this would be the simple Poisson process $\mathrm{pro}_i\rightarrow \mathrm{pro}$. More complicated reactivation mechanisms can of course be contemplated \cite{JL2023}. Given the distributions described above, genes displaying bursting could be studied using delay-stochastic simulations, with the single-polymerase transcription time distribution as a key element.

As we have emphasized, the construction of the delay-stochastic models described above is facilitated by the modularity of the transcription model studied here. As noted by Greive and coworkers for a related model \cite{greive2011development}, modularity means that we can mix and match different descriptions of initiation, elongation, pausing and termination to suit any given gene. The mathematical analysis of these models is similarly modular: we can take the pieces derived here and combine them as needed by convolution.

As explained in Section \ref{sec:Eukaryoticmodel}, we have substantially simplified Vashishtha's model. Studying the effects of various omitted details on transcription dynamics is of course a potential avenue for extending this work.
The experimental observation of, as Larson and coworkers  put it, ``deterministic'' kinetics during the elongation phase of transcription \cite{Larson-etal2011} suggests that a more refined model of productive elongation would not behave differently than the simple model studied here. This is a consequence of the law of large numbers (LLN)~\cite{LMB1978LLN} and of the large sizes of many eukaryotic genes: If the kinetics are identical at every nucleotide, then $\tau_e = \sum_i\tau_i=N\bar{\tau_i}$, where $\bar{\tau_i}$ is the sample mean of the residence times at each nucleotide. The LLN says that $\bar{\tau_i}$ converges to the mean of the underlying distribution. Thus, $\tau_e$ converges to a fixed value for long genes. If the residence times are not identically distributed, then there are generalizations of the law of large numbers leading to the same conclusion provided the $\tau_i$ satisfy some additional (physically reasonable) conditions~\cite{Janisch2021}: $\tau_e$ converges to a constant for sufficiently large $N$. Pausing is a different matter, as Voliotis and coworkers have shown \cite{voliotis2008fluctuations}: backtracked pausing can contribute a heavy tail to the distribution of transcription times. To model genes that are prone to backtracking pauses in specific regions of the sequence, it would be possible to combine the Voliotis \emph{et al.}\ pausing module with our model to obtain the corresponding distribution of elongation times.

Shortly after initiation, RNAPII tends to pause 20--40\,nt downstream of the start site  in an event known as promoter-proximal pausing~\cite{Li-etal1996}. A promoter-proximal pause that is not subject to specific regulation is relatively easy to model since it can be dealt with using the splitting trick illustrated in Figure~\ref{longpause}. The modularity of the model in this case means that the exact location of the pause site is irrelevant. The situation is different if the polymerase is released from its promoter-proximal pause by a specific release factor with its own dynamics. Consider again model~(\ref{model_w_PC}) assuming that the concentration of the release factor cannot be treated as being constant in time. This will be the case if (e.g.) a signaling cascade modulates the activity of the release factor. In a dynamic model of such a system, we would have to explicitly include reaction~\eqref{release_from_pause}  so that the time dependence of the release from the pause could be properly captured. Thus, we would need to know the distribution of~$\tau_p$ since this interval of time necessarily precedes the action of the release factor.  In particular, note that the concentration of the release factor could change during the time interval~$\tau_p$. However, we probably cannot treat the initial elongation leading to the pause using the theory of Section~\ref{secelongation} because only a short run of nucleotides precedes promoter-proximal pausing, thus violating the central assumption in the development of the advection equation, namely that $N$ is large. We could solve the first-passage time for reaching the pause site, but even for the moderate number of nucleotides separating the start and promoter-proximal pause sites, this is cumbersome unless the kinetics are particularly simple. An additional difficulty is that the promoter-proximal pause probably occurs in a region covering several nucleotides, and not at a specific nucleotide in any given gene. Fortunately, compared to typical initiation times (Figure~\ref{rho3}), the contribution of the elongation time to $\tau_p$ may well be negligible despite the slower elongation rate in this region of a gene~\cite{FKG2001}. Even if the time required to reach the promoter-proximal pausing site is not negligible, it is likely possible to develop useful approximations to the distribution of~$\tau_p$.

Termination is another area ripe for investigation in terms of its effect on transcription dynamics as it is quite slow in eukaryotes \cite{Darzacq-etal2007}. Indeed, we have initiated such studies in our laboratory \cite{Murphy2017}. It is also interesting in that multiple independent mechanisms cooperate to ensure termination \cite{LB2015}, introducing an additional element of stochasticity right at the end of the process of transcription.

Even within the context of the model studied here, we have not exhausted the possibilities for further study. In particular, while we have developed the relevant equation of motion, we have not studied in much detail the effects of variable kinetic parameters on the density of polymerases. At higher initiation rates, this could have an effect on polymerase-polymerase interactions. Some preliminary calculations on a system in which the polymerase slows down dramatically in one region of the gene have been carried out \cite{Hosseini2016}, but this hardly scratches the surface of what could be done.

The delays that appear in models of gene expression are measurable quantities. There has been recent progress towards inferring these delays from experimental data \cite{CCC2020,HCC2023}, with some caveats: depending on what is measured, not all kinetic parameters, indeed not all delays, are identifiable. This creates an opportunity for a two-way discussion with experiment: Clearly, these experimental measurements can provide data to be used in modeling. In addition, the current contribution and similar efforts in other groups can provide parameterizable forms for the delay distributions that can be used to fit experimental results, and thereby to obtain estimates of some kinetic parameters that determine these distributions. 

As noted in the introduction, the eventual goal of this research is to provide delay distributions that can be used in models of gene expression. In addition to the distribution of transcription times which is the topic of this contribution, we will also need distributions for the nuclear export time, on which some progress has been made \cite{HS2014,RT2012,Tang2010} but which remains very much an open problem, and for translation \cite{garai2009stochastic,GCR2009,TSH2010}. We can look forward to soon having all the pieces necessary either for the direct inclusion in a model of these three major steps in eukaryotic gene expression, or at the very least to enable rational decisions about whether they need to be explicitly included in a given model.

%%\section*{CRediT authorship contribution statement}
%\textbf{S. Hossein Hosseini:} Formal analysis, Investigation, Methodology, Validation, Visualization, Writing --- original draft.
%\textbf{Marc R. Roussel:} Conceptualization, Funding acquisition, Formal analysis, Methodology, Resources, Writing --- review \& editing.

\section*{Use of AI tools declaration}
The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.
\section*{Acknowledgments}
This work was funded by the Natural Sciences and Engineering Research Council (NSERC) of Canada.

\section*{Conflict of interest}
The authors declare there is no conflict of interest.

\providecommand{\href}[2]{#2}
\providecommand{\arxiv}[1]{\href{http://arxiv.org/abs/#1}{arXiv:#1}}
\providecommand{\url}[1]{\texttt{#1}}
\providecommand{\urlprefix}{URL }
\begin{thebibliography}{100}
\setlength{\baselineskip}{14pt} 

\bibitem{Larson-etal2011}
 D.~R. Larson, D.~Zenklusen, B.~Wu, J.~A. Chao, R.~H. Singer, Real-time observation of transcription initiation and elongation on  an endogeneous yeast gene, \emph{Science}, \textbf{332} (2011), 475--478. \doilink{https://doi.org/10.1126/science.1202142}


\bibitem{Buratowski-etal1989}
 S.~Buratowski, S.~Hahn, L.~Guarente, P.~A. Sharp, Five intermediate complexes in transcription initiation by {RNA}  polymerase {II}, \emph{Cell}, \textbf{56} (1989), 549--561. \doilink{https://doi.org/10.1016/0092-8674(89)90578-3}


\bibitem{Dvir2001209}
 A.~Dvir, J.~W. Conaway, R.~C. Conaway, Mechanism of transcription initiation and promoter escape by {RNA} polymerase {II}, \emph{Curr. Opin. Genet. Dev.}, \textbf{11} (2001), 209--214. \doilink{https://doi.org/10.1016/S0959-437X(00)00181-7}



\bibitem{Darzacq-etal2007}
 X.~Darzacq, Y.~Shav-Tal, V.~de~Turris, Y.~Brody, S.~M. Shenoy, R.~D.
 Phair, R.~H. Singer, \emph{In vivo} dynamics of {RNA} polymerase {II} transcription, \emph{Nat. Struct. Mol. Biol.}, \textbf{14} (2007), 796--806.  \doilink{https://doi.org/10.1038/nsmb1280}



\bibitem{Venter-etal2001}
 J.~C. Venter, M.~D. Adams, E.~W. Myers, P.~W. Li, R.~J. Mural, G.~G.
 Sutton, et al., The sequence of the human genome, \emph{Science},
 \textbf{291} (2001), 1304--1351.
 \doilink{https://doi.org/10.1126/science.1058040}

\bibitem{tennyson1995human}
 C.~N. Tennyson, H.~J. Klamut, R.~G. Worton, The human dystrophin gene requires 16 hours to be transcribed and is  cotranscriptionally spliced, \emph{Nat. Genet.}, \textbf{9} (1995), 184--190.
 \doilink{https://doi.org/10.1038/ng0295-184}


\bibitem{LB2015}
 J.-F. Lemay, F.~Bachand, Fail-safe transcription termination: Because one is never enough, \emph{RNA Biol.}, \textbf{12} (2015), 927--932.  \doilink{https://doi.org/10.1080/15476286.2015.1073433}




\bibitem{BS2019}
 R.~Ben-Yishay, Y.~Shav-Tal, The dynamic lifecycle of {mRNA} in the nucleus, \emph{Curr. Opin. Cell Biol.}, \textbf{58} (2019), 69--75.   \doilink{https://doi.org/10.1016/j.ceb.2019.02.007}

\bibitem{Daneholt2001}
 B.~Daneholt, Assembly and transport of a premessenger {RNP} particle, \emph{Proc. Natl. Acad. Sci. U.S.A.}, \textbf{98} (2001), 7012--7017. \doilink{https://doi.org/10.1073/pnas.111145498}

\bibitem{SS2013}
 J.~Sheinberger, Y.~Shav-Tal, The dynamic pathway of nuclear {RNA} in eukaryotes, \emph{Nucleus}, \textbf{4} (2013), 195--205.   \doilink{https://doi.org/10.4161/nucl.24434}

\bibitem{CDD2020}
 A.~Chaudhuri, S.~Das, B.~Das, Localization elements and zip codes in the intracellular transport
  and localization of messenger {RNAs} in \emph{Saccharomyces cerevisiae},
  \emph{WIREs RNA}, \textbf{11} (2020), e1591. \doilink{https://doi.org/10.1002/wrna.1591}

\bibitem{Schmidt-etal2011}
 U.~Schmidt, E.~Basyuk, M.-C. Robert, M.~Yoshida, J.-P. Villemin,  D.~Auboeuf, et al.,
 Real-time imaging of cotranscriptional splicing reveals a kinetic  model that reduces noise: Implications for alternative splicing regulation, \emph{J. Cell Biol.}, \textbf{193} (2011), 819--829. \doilink{https://doi.org/10.1083/jcb.201009012}

\bibitem{Cramer-etal2001}
 P.~Cramer, A.~Srebrow, S.~Kadener, S.~Werbajh, M.~de~la Mata, G.~Melen, et al., Coordination between transcription and pre-{mRNA} processing, \emph{FEBS Lett.}, \textbf{498} (2001), 179--182. \doilink{https://doi.org/10.1016/S0014-5793(01)02485-1}

\bibitem{BDS2012}
 A.~Babour, C.~Dargemont, F.~Stutz, Ubiquitin and assembly of export competent {mRNP}, \emph{Biochim. Biophys. Acta}, \textbf{1819} (2012), 521--530. \doilink{https://doi.org/10.1016/j.bbagrm.2011.12.006}



\bibitem{CP1997}
 R.~A. Coleman, B.~F. Pugh, Slow dimer dissociation of the {TATA} binding protein dictates the  kinetics of {DNA} binding, \emph{Proc. Natl. Acad. Sci. U.S.A.}, \textbf{94} (1997), 7221--7226.
\doilink{https://doi.org/10.1073/pnas.94.14.7221}

\bibitem{KG2000}
 J.~F. Kugel, J.~A. Goodrich, A kinetic model for the early steps of {RNA} synthesis by human {RNA}  polymerase {II}, \emph{J. Biol. Chem.}, \textbf{275} (2000), 40483--40491. 
\doilink{https://doi.org/10.1074/jbc.M006401200}

\bibitem{KKS2015}
 A.~Kalo, I.~Kanter, A.~Shraga, J.~Sheinberger, H.~Tzemach, N.~Kinor, et al.,  Cellular levels of signaling factors are sensed by $\beta$-actin  alleles to modulate transcriptional pulse intensity, \emph{Cell Rep.}, \textbf{11} (2015), 419--432.  \doilink{https://doi.org/10.1016/j.celrep.2015.03.039}

\bibitem{BPM82}
 R.~D. Bliss, P.~R. Painter, A.~G. Marr, Role of feedback inhibition in stabilizing the classical operon, \emph{J. Theor. Biol.}, \textbf{97} (1982), 177--193. \doilink{https://doi.org/10.1016/0022-5193(82)90098-4}

\bibitem{BS1987}
 F.~Buchholtz, F.~W. Schneider, Computer simulation of {T3}/{T7} phage infection using lag times,
 \emph{Biophys. Chem.}, \textbf{26} (1987), 171--179.  \doilink{https://doi.org/10.1016/0301-4622(87)80020-0}

\bibitem{BM1988}
 S.~N. Busenberg, J.~M. Mahaffy, The effects of dimension and size for a compartmental model of
  repression, \emph{SIAM J. Appl. Math.}, \textbf{48} (1988), 882--903.  \doilink{https://doi.org/10.1137/0148049}

\bibitem{Lewis2003}
 J.~Lewis, Autoinhibition with transcriptional delay: A simple mechanism for the  zebrafish somitogenesis oscillator, \emph{Curr. Biol.}, \textbf{13} (2003), 1398--1408.  \doilink{https://doi.org/10.1016/S0960-9822(03)00534-7}

\bibitem{Monk2003}
 N.~A.~M. Monk, Oscillatory expression of {Hes1}, p53, and {NF}-$\kappa${B} driven by  transcriptional time delays, \emph{Curr. Biol.}, \textbf{13} (2003), 1409--1413. \doilink{ https://doi.org/10.1016/S0960-9822(03)00494-9}

\bibitem{CLW2022}
 L.-J. Chiu, M.-Y. Ling, E.-H. Wu, C.-X. You, S.-T. Lin, C.-C. Shu, The distributed delay rearranges the bimodal distribution at protein  level, \emph{J. Taiwan Inst. Chem. Eng.}, \textbf{137} (2022), 104436. \doilink{https://doi.org/10.1016/j.jtice.2022.104436}


\bibitem{JP2015}
 M.~Jansen, P.~Pfaffelhuber, Stochastic gene expression with delay, \emph{J. Theor. Biol.}, \textbf{364} (2015), 355--363.  \doilink{https://doi.org/10.1016/j.jtbi.2014.09.031}

\bibitem{RW2007}
 K.~Rateitschak, O.~Wolkenhauer, Intracellular delay limits cyclic changes in gene expression,
 \emph{Math. Biosci.}, \textbf{205} (2007), 163--179.  \doilink{https://doi.org/10.1016/j.mbs.2006.08.010}

\bibitem{roussel2013distribution}
 M.~R. Roussel, On the distribution of transcription times, \emph{BIOMATH}, \textbf{2} (2013), 1307247. \doilink{https://doi.org/10.11145/j.biomath.2013.07.247}

\bibitem{roussel2006stochastic}
 M.~R. Roussel, R.~Zhu, Stochastic kinetics description of a simple transcription model, \emph{Bull. Math. Biol.}, \textbf{68} (2006), 1681--1713.  \doilink{https://doi.org/10.1007/s11538-005-9048-6}

\bibitem{Vashishtha}
 S.~Vashishtha, \emph{Stochastic modeling of eukaryotic transcription at the single
  nucleotide level}, {M.Sc.}\ thesis, University of Lethbridge, 2011, \urlprefix\url{https://www.uleth.ca/dspace/handle/10133/3190}.

\bibitem{PCP2010}
 V.~Pelechano, S.~Ch{\'a}vez, J.~E. P{\'e}rez-Ort{\'\i}n, A complete set of nascent transcription rates for yeast genes, \emph{PLoS One}, \textbf{5} (2010), e15442. \doilink{https://doi.org/10.1371/journal.pone.0015442}

\bibitem{MCG2012}
 T.~Muramoto, D.~Cannon, M.~Gierli{\'n}ski, A.~Corrigan, G.~J. Barton, J.~R. Chubb,  Live imaging of nascent {RNA} dynamics reveals distinct types of  transcriptional pulse regulation, \emph{Proc. Natl. Acad. Sci. U.S.A.}, \textbf{109} (2012),  7350--7355. \doilink{https://doi.org/10.1073/pnas.1117603109}

\bibitem{RPT+2006}
 A.~Raj, C.~S. Peskin, D.~Tranchina, D.~Y. Vargas, S.~Tyagi, Stochastic {mRNA} synthesis in mammalian cells, \emph{PLoS Biol.}, \textbf{4} (2006), e309.  \doilink{https://doi.org/10.1371/journal.pbio.0040309}

\bibitem{Suter-etal2011}
 D.~M. Suter, N.~Molina, D.~Gatfield, K.~Schneider, U.~Schibler,  F.~Naef, Mammalian genes are transcribed with widely different bursting  kinetics, \emph{Science}, \textbf{332} (2011), 472--474. \doilink{https://doi.org/10.1126/science.1198817}

\bibitem{JKL2014}
 I.~Jonkers, H.~Kwak, J.~T. Lis, Genome-wide dynamics of {Pol II} elongation and its interplay with  promoter proximal pausing, chromatin, and exons, \emph{eLife}, \textbf{3} (2014), e02407.
\doilink{https://doi.org/10.7554/eLife.02407}

\bibitem{PBS+2018}
 P.~K. Parua, G.~T. Booth, M.~Sans{\'o}, B.~Benjamin, J.~C. Tanny,  J.~T.
 Lis, R.~P. Fisher,
 A {Cdk9}-{PP1} switch regulates the elongation-termination transition  of {RNA} polymerase {II}, \emph{Nature}, \textbf{558} (2018), 460--464.  \doilink{https://doi.org/10.1038/s41586-018-0214-z}


\bibitem{bai2007mechanochemical}
 L.~Bai, R.~M. Fulbright, M.~D. Wang, Mechanochemical kinetics of transcription elongation, \emph{Phys. Rev. Lett.}, \textbf{98} (2007), 068103. \doilink{https://doi.org/10.1103/PhysRevLett.98.068103}

\bibitem{bai2004sequence}
 L.~Bai, A.~Shundrovsky, M.~D. Wang, Sequence-dependent kinetic model for transcription elongation by  {RNA} polymerase, \emph{J. Mol. Biol.}, \textbf{344} (2004), 335--349. \doilink{https://doi.org/10.1016/j.jmb.2004.08.107}

\bibitem{julicher1998motion}
 F.~J{\"u}licher, R.~Bruinsma, Motion of {RNA} polymerase along {DNA}: a stochastic model, \emph{Biophys. J.}, \textbf{74} (1998), 1169--1185.  \doilink{https://doi.org/10.1016/S0006-3495(98)77833-6}

\bibitem{wang1998force}
 H.-Y. Wang, T.~Elston, A.~Mogilner, G.~Oster, Force generation in {RNA} polymerase, \emph{Biophys. J.}, \textbf{74} (1998), 1186--1202. \doilink{https://doi.org/10.1016/S0006-3495(98)77834-8}

\bibitem{yager1991thermodynamic}
 T.~D. Yager, P.~H. Von~Hippel, A thermodynamic analysis of {RNA} transcript elongation and  termination in \emph{Escherichia coli}, \emph{Biochemistry}, \textbf{30} (1991), 1097--1118.
\doilink{https://doi.org/10.1021/bi00218a032}

\bibitem{greive2011development}
 S.~J. Greive, J.~P. Goodarzi, S.~E. Weitzel, P.~H. von Hippel, Development of a ``modular'' scheme to describe the kinetics of  transcript elongation by {RNA} polymerase, \emph{Biophys. J.}, \textbf{101} (2011), 1155--1165. \doilink{https://doi.org/10.1016/j.bpj.2011.07.042}

\bibitem{FPG2021}
 T.~Filatova, N.~Popovic, R.~Grima, Statistics of nascent and mature RNA fluctuations in a stochastic  model of transcriptional initiation, elongation, pausing, and termination, \emph{Bull. Math. Biol.}, \textbf{83} (2021), 3. \doilink{ https://doi.org/10.1007/s11538-020-00827-7}

\bibitem{BRE2011}
 A.~N. Boettiger, P.~L. Ralph, S.~N. Evans, Transcriptional regulation: Effects of promoter proximal pausing on  speed, synchrony and reliability, \emph{PLoS Comput. Biol.}, \textbf{7} (2011), e1001136.  \doilink{ https://doi.org/10.1371/journal.pcbi.1001136}

\bibitem{XKK2013}
 X.~Xu, N.~Kumar, A.~Krishnan, R.~V. Kulkarni, Stochastic modeling of dwell-time distributions during  transcriptional pausing and initiation, in \emph{52nd IEEE Conference on Decision and Control}, 2013, 4068--4073.

\bibitem{Hamano2016}
 M.~Hamano, Stochastic transcription elongation via rule based modelling, \emph{Electron. Notes Theor. Comput. Sci.}, \textbf{326} (2016),  73--88. \doilink{https://doi.org/10.1016/j.entcs.2016.09.019}

\bibitem{Klumpp18112008}
 S.~Klumpp, T.~Hwa, Stochasticity and traffic jams in the transcription of
 ribosomal  {RNA}: Intriguing role of termination and antitermination,
 \emph{Proc. Natl. Acad. Sci. U.S.A.}, \textbf{105} (2008), 18159--18164.
 \doilink{https://doi.org/10.1073/pnas.0806084105}

\bibitem{ribeiro2009delayed}
 A.~S. Ribeiro, O.-P. Smolander, T.~Rajala, A.~H{\"a}kkinen,  O.~Yli-Harja, Delayed stochastic model of transcription at the single nucleotide  level, \emph{J. Comput. Biol.}, \textbf{16} (2009),  539--553.  \doilink{https://doi.org/10.1089/cmb.2008.0153}


\bibitem{SN2010}
 M.~J. Schilstra, C.~L. Nehaniv, Stochastic model of template-directed elongation processes in  biology, \emph{BioSystems}, \textbf{102} (2010), 55--60.   \doilink{https://doi.org/10.1016/j.biosystems.2010.07.006}


\bibitem{garai2009stochastic}
 A.~Garai, D.~Chowdhury, D.~Chowdhury, T.~V. Ramakrishnan, Stochastic kinetics of ribosomes: Single motor properties and collective behavior, \emph{Phys. Rev. E}, \textbf{80} (2009), 011908.  \doilink{https://doi.org/10.1103/PhysRevE.80.011908}

\bibitem{GCR2009}
 A.~Garai, D.~Chowdhury, T.~V. Ramakrishnan, Fluctuations in protein synthesis from a single {RNA} template:  Stochastic kinetics of ribosomes, \emph{Phys. Rev. E}, \textbf{79} (2009), 011916.  \doilink{https://doi.org/10.1103/PhysRevE.79.011916}
 

\bibitem{TSH2010}
 L.~Mier{-}y{-}Ter{\'a}n{-}Romero, M.~Silber, V.~Hatzimanikatis, The origins of time-delay in template biopolymerization processes, \emph{PLoS Comput. Biol.}, \textbf{6} (2010), e1000726. \doilink{https://doi.org/10.1371/journal.pcbi.1000726}

\bibitem{CW2011}
 L.~S. Churchman, J.~S. Weissman, Nascent transcript sequencing visualizes
 transcription at nucleotide  resolution, \emph{Nature}, \textbf{469}
 (2011), 368--373. \doilink{https://doi.org/10.1038/nature09652}

\bibitem{Neuman2003}
 K.~C. Neuman, E.~A. Abbondanzieri, R.~Landick, J.~Gelles, S.~M.  Block, Ubiquitous transcriptional pausing is independent of {RNA} polymerase  backtracking, \emph{Cell}, \textbf{115} (2003), 437 -- 447. \doilink{https://doi.org/10.1016/S0092-8674(03)00845-6}

\bibitem{landick2006regulatory}
 R.~Landick, The regulatory roles and mechanism of transcriptional pausing, \emph{Biochem. Soc. Trans.}, \textbf{34} (2006), 1062--1066. \doilink{https://doi.org/10.1042/BST0341062}

\bibitem{epshtein2003transcription}
 V.~Epshtein, F.~Toulm{\'e}, A.~R. Rahmouni, S.~Borukhov,  E.~Nudler, Transcription through the roadblocks: the role of {RNA} polymerase  cooperation, \emph{EMBO J.}, \textbf{22} (2003), 4719--4727. \doilink{https://doi.org/10.1093/emboj/cdg452}

\bibitem{klumpp2011pausing}
 S.~Klumpp, Pausing and backtracking in transcription under dense traffic  conditions, \emph{J. Stat. Phys.}, \textbf{142} (2011), 1252--1267.  \doilink{https://doi.org/10.1007/s10955-011-0120-3}


\bibitem{voliotis2008fluctuations}
 M.~Voliotis, N.~Cohen, C.~Molina-Par{\'\i}s, T.~B. Liverpool, Fluctuations, pauses, and backtracking in {DNA} transcription, \emph{Biophys. J.}, \textbf{94} (2008), 334--348.
 \doilink{https://doi.org/10.1529/biophysj.107.105767}


\bibitem{LG2011}
 J.~Li, D.~S. Gilmour, Promoter proximal pausing and the control of gene expression, \emph{Curr. Opin. Genet. Dev.}, \textbf{21} (2011), 231--235.
 \doilink{https://doi.org/10.1016/j.gde.2011.01.010}


\bibitem{Nechaev201134}
 S.~Nechaev, K.~Adelman, Pol {II} waiting in the starting gates: Regulating the transition  from transcription initiation into productive elongation, \emph{Biochim. Biophys. Acta}, \textbf{1809} (2011), 34--45.
\doilink{https://doi.org/10.1016/j.bbagrm.2010.11.001}

\bibitem{Rahl2010432}
 P.~B. Rahl, C.~Y. Lin, A.~C. Seila, R.~A. Flynn, S.~McCuine, C.~B.  Burge, et al.,
 c-{Myc} regulates transcriptional pause release, \emph{Cell}, \textbf{141} (2010), 432--445. \doilink{https://doi.org/10.1016/j.cell.2010.03.030}

\bibitem{FXF2021}
 P.~Feng, A.~Xiao, M.~Fang, F.~Wan, S.~Li, P.~Lang, et al., A machine learning-based framework for modeling transcription  elongation, \emph{Proc. Natl. Acad. Sci. U.S.A.}, \textbf{118} (2021),  e2007450118.  \doilink{https://doi.org/10.1073/pnas.2007450118}


\bibitem{ZBI+2012}
 B.~Zamft, L.~Bintu, T.~Ishibashi, C.~Bustamante, Nascent {RNA} structure modulates the transcriptional dynamics of  {RNA} polymerases, \emph{Proc. Natl. Acad. Sci. U.S.A.}, \textbf{109} (2012),  8948--8953. \doilink{https://doi.org/10.1073/pnas.1205063109}

\bibitem{AIB+2010}
 R.~D. Alexander, S.~A. Innocente, J.~D. Barrass, J.~D. Beggs, Splicing-dependent {RNA} polymerase pausing in yeast, \emph{Mol. Cell}, \textbf{40} (2010), 582--593. \doilink{https://doi.org/10.1016/j.molcel.2010.11.005}

\bibitem{GWP2006}
 N.~Gromak, S.~West, N.~J. Proudfoot, Pause sites promote transcriptional termination of mammalian {RNA}  polymerase {II}, \emph{Mol. Cell. Biol.}, \textbf{26} (2006), 3986--3996. \doilink{https://doi.org/10.1128/MCB.26.10.3986-3996.2006}

\bibitem{MacDonald1976}
 N.~MacDonald, Time delay in prey-predator models, \emph{Math. Biosci.}, \textbf{28} (1976), 321--330. \doilink{https://doi.org/10.1016/0025-5564(76)90130-9}

\bibitem{MacDonald1977a}
 N.~MacDonald, Time lag in a model of a biochemical reaction sequence with end  product inhibition, \emph{J. Theor. Biol.}, \textbf{67} (1977), 549--556. \doilink{https://doi.org/10.1016/0022-5193(77)90056-X}

\bibitem{MacDonald89}
 N.~MacDonald, \emph{Biological Delay Systems: Linear Stability Theory}, Cambridge, Cambridge, 1989.


\bibitem{BLM2013}
 M.~Barrio, A.~Leier, T.~T. Marquez-Lago, Reduction of chemical reaction networks through delay distributions, \emph{J. Chem. Phys.}, \textbf{138} (2013), 104114. \doilink{https://doi.org/10.1063/1.4793982}

\bibitem{LBM2014}
 A.~Leier, M.~Barrio,  T.~T. Marquez-Lago, Exact model reduction with delays: Closed-form distributions and  extensions to fully bi-directional monomolecular reactions, \emph{J. R. Soc. Interface}, \textbf{11} (2014), 20140108. \doilink{https://doi.org/10.1098/rsif.2014.0108}

\bibitem{Epstein90}
 I.~R. Epstein, Differential delay equations in chemical kinetics: Some simple linear  model systems, \emph{J. Chem. Phys.}, \textbf{92} (1990), 1702--1712.  \doilink{https://doi.org/10.1063/1.458052}


\bibitem{Bratsun-etal2005}
 D.~Bratsun, D.~Volfson, L.~S. Tsimring, J.~Hasty, Delay-induced stochastic oscillations in gene regulation, \emph{Proc. Natl. Acad. Sci. U.S.A.}, \textbf{102} (2005), 14593--14598.  \doilink{https://doi.org/10.1073/pnas.0503858102}


\bibitem{roussel2006validation}
 M.~R. Roussel, R.~Zhu, Validation of an algorithm for delay stochastic simulation of  transcription and translation in prokaryotic gene expression, \emph{Phys. Biol.}, \textbf{3} (2006), 274.  \doilink{https://doi.org/10.1088/1478-3975/3/4/005}

\bibitem{Jennings2013}
 B.~H. Jennings, Pausing for thought: Disrupting the early transcription elongation  checkpoint leads to developmental defects and tumourigenesis, \emph{BioEssays}, \textbf{35} (2013), 553--560.  \doilink{https://doi.org/10.1002/bies.201200179}

\bibitem{KFC2013}
 H.~Kwak, N.~J. Fuda, L.~J. Core, J.~T. Lis, Precise maps of {RNA} polymerase reveal how promoters direct  initiation and pausing, \emph{Science}, \textbf{339} (2013), 950--953. \doilink{https://doi.org/10.1126/science.1229386}

\bibitem{HBG+2006}
 A.~R. Hieb, S.~Baran, J.~A. Goodrich, J.~F. Kugel, An 8nt {RNA} triggers a rate-limiting shift of {RNA} polymerase {II}  complexes into elongation, \emph{EMBO J.}, \textbf{25} (2006), 3100--3109.
\doilink{https://doi.org/10.1038/sj.emboj.7601197}

\bibitem{SHS+2014}
 T.~J. Stasevich, Y.~Hayashi-Takanaka, Y.~Sato, K.~Maehara, Y.~Ohkawa, K.~Sakata-Sogawa, et al.,
 Regulation of {RNA} polymerase {II} activation by histone acetylation  in single living cells,
 \emph{Nature}, \textbf{516} (2014), 272--275.  \doilink{https://doi.org/10.1038/nature13714}

\bibitem{SJG+2018}
 B.~Steurer, R.~C. Janssens, B.~Geverts, M.~E. Geijer, F.~Wienholz,  A.~F. Theil, et al., Live-cell analysis of endogeneous {GFP-RPB1} uncovers rapid turnover  of initiating and promoter-paused {RNA} polymerase {II}, \emph{Proc. Natl. Acad. Sci. U.S.A.}, \textbf{115} (2018),  E4368--E4376.
\doilink{https://doi.org/10.1073/pnas.1717920115}

\bibitem{LHE2021}
 J.~Liu, D.~Hansen, E.~Eck, Y.~J. Kim, M.~Turner, S.~Alamos, H.~G. Garcia, Real-time single-cell characterization of the eukaryotic  transcription cycle reveals correlations between {RNA} initiation,  elongation, and cleavage, \emph{PLoS Comput. Biol.}, \textbf{17} (2021), e1008999.
\doilink{https://doi.org/10.1371/journal.pcbi.1008999}

\bibitem{Kremling2007}
 A.~Kremling, Comment on mathematical models which describe transcription and  calculate the relationship between mRNA and protein expression ratio, \emph{Biotech. Bioeng.}, \textbf{96} (2007), 815--819.  \doilink{https://doi.org/10.1002/bit.21065}


\bibitem{MP2013}
 N.~Mitarai, S.~Pedersen, Control of ribosome traffic by position-dependent choice of  synonymous codons, \emph{Phys. Biol.}, \textbf{10} (2013), 056011.  \doilink{https://doi.org/10.1088/1478-3975/10/5/056011}

\bibitem{Roussel96b}
 M.~R. Roussel, The use of delay differential equations in chemical kinetics, \emph{J. Phys. Chem.}, \textbf{100} (1996), 8323--8330.  \doilink{https://doi.org/10.1021/jp9600672}


\bibitem{CP1995}
 R.~A. Coleman, B.~F. Pugh, Evidence for functional binding and stable sliding of the {TATA}  binding protein on nonspecific {DNA}, \emph{J. Biol. Chem.}, \textbf{270} (1995), 13850--13859.
\doilink{https://doi.org/10.1074/jbc.270.23.13850}

\bibitem{DJS+2005}
 A.~Dasgupta, S.~A. Juedes, R.~O. Sprouse, D.~T. Auble, Mot1-mediated control of transcription complex assembly and activity, \emph{EMBO J.}, \textbf{24} (2005), 1717--1729. \doilink{https://doi.org/10.1038/sj.emboj.7600646}

\bibitem{SKM+2008}
 R.~O. Sprouse, T.~S. Karpova, F.~Mueller, A.~Dasgupta, J.~G. McNally, D.~T. Auble, Regulation of {TATA}-binding protein dynamics in living yeast cells, \emph{Proc. Natl. Acad. Sci. U.S.A.}, \textbf{105} (2008),  13304--13308. \doilink{https://doi.org/10.1073/pnas.0801901105}

\bibitem{Hosseini2016}
 S.~H. Hosseini, \emph{Analytic Solutions for Stochastic Models of Transcription},
 Master's thesis, University of Lethbridge, 2016, \urlprefix\url{https://www.uleth.ca/dspace/handle/10133/4791}.


\bibitem{woo2006analytical}
 H.-J. Woo, Analytical theory of the nonequilibrium spatial distribution of {RNA}  polymerase translocations, \emph{Phys. Rev. E}, \textbf{74} (2006), 011907.  \doilink{https://doi.org/10.1103/PhysRevE.74.011907}






\bibitem{LZK+2012}
 M.~H. Larson, J.~Zhou, C.~D. Kaplan, M.~Palangat, R.~D. Kornberg,
 R.~Landick, S.~M. Block,
 Trigger loop dynamics mediate the balance between the transcriptional  fidelity and speed of {RNA} polymerase {II}, \emph{Proc. Natl. Acad. Sci. U.S.A.}, \textbf{109} (2012),  6555--6560.  \doilink{https://doi.org/10.1073/pnas.1200939109}


\bibitem{CC2011}
 A.~C.~M. Cheung, P.~Cramer, Structural basis of {RNA} polymerase {II} backtracking, arrest and
  reactivation, \emph{Nature}, \textbf{471} (2011), 249--253.  \doilink{https://doi.org/10.1038/nature09785}


\bibitem{BS2015}
 G.~Brzy{\.z}ek, S.~{\'S}wie{\.z}ewski, Mutual interdependence of splicing and transcription elongation, \emph{Transcription}, \textbf{6} (2015), 37--39. \doilink{ https://doi.org/10.1080/21541264.2015.1040146}

\bibitem{IKL+2013}
 M.~Imashimizu, M.~L. Kireeva, L.~Lubkowska, D.~Gotte, A.~R. Parks, J.~N.
 Strathem, M.~Kashlev,
 Intrinsic translocation barrier as an initial step in pausing by {RNA} polymerase {II}, \emph{J. Mol. Biol.}, \textbf{425} (2013), 697--712. \doilink{https://doi.org/10.1016/j.jmb.2012.12.002}

\bibitem{Roberts2014}
 J.~W. Roberts, Molecular basis of transcription pausing, \emph{Science}, \textbf{344} (2014), 1226--1227.  \doilink{https://doi.org/10.1126/science.1255712}


\bibitem{SP2009}
 J.~Singh, R.~A. Padgett, Rates of \emph{in situ} transcription and splicing in large human  genes, \emph{Nat. Struct. Mol. Biol.}, \textbf{16} (2009), 1128--1133.  \doilink{https://doi.org/10.1038/nsmb.1666}




\bibitem{rosonina2006terminating}
 E.~Rosonina, S.~Kaneko, J.~L. Manley, Terminating the transcript: breaking up is hard to do, \emph{Genes Dev.}, \textbf{20} (2006), 1050--1056. \doilink{https://doi.org/10.1101/gad.1431606}

\bibitem{abbondanzieri2005direct}
 E.~A. Abbondanzieri, W.~J. Greenleaf, J.~W. Shaevitz, R.~Landick,  S.~M. Block, Direct observation of base-pair stepping by {RNA} polymerase, \emph{Nature}, \textbf{438} (2005), 460--465. \doilink{https://doi.org/10.1038/nature04268}

\bibitem{Hsu2002}
 L.~M. Hsu, Promoter clearance and escape in prokaryotes, \emph{Biochim. Biophys. Acta}, \textbf{1577} (2002), 191--207. \doilink{https://doi.org/10.1016/S0167-4781(02)00452-9}

\bibitem{KSC2002}
 H.~Kimura, K.~Sugaya, P.~R. Cook, The transcription cycle of {RNA} polymerase {II} in living cells, \emph{J. Cell Biol.}, \textbf{159} (2002), 777--782. \doilink{https://doi.org/10.1083/jcb.200206019}

\bibitem{FKG2001}
 H.~A. Ferguson, J.~F. Kugel, J.~A. Goodrich, Kinetic and mechanistic analysis of the {RNA} polymerase {II}  transcription reaction at the human interleukin-2 promoter, \emph{J. Mol. Biol.}, \textbf{314} (2001), 993--1006. \doilink{https://doi.org/10.1006/jmbi.2000.5215}

\bibitem{JIM1998}
 D.~A. Jackson, F.~J. Iborra, E.~M.~M. Manders, P.~R. Cook, Numbers and organization of {RNA} polymerases, nascent transcripts, and transcription units in {HeLa} nuclei, \emph{Mol. Biol. Cell}, \textbf{9} (1998), 1523--1536.   \doilink{https://doi.org/10.1091/mbc.9.6.1523}


\bibitem{HK2019}
 P.~J. Hurtado, A.~S. Kirosingh, Generalizations of the `linear chain trick': Incorporating more  flexible dwell time distributions into mean field {ODE} models, \emph{J. Math. Biol.}, \textbf{79} (2019), 1831--1883.  \doilink{https://doi.org/10.1007/s00285-019-01412-w}


\bibitem{Goldstein-wave-eq}
 H.~Golstein, \emph{Classical Mechanics}, chapter~12, Addison-Wesley, Reading, Massachusetts, 1980.


\bibitem{Othmer1983}
 H.~G. Othmer, A continuum model for coupled cells, \emph{J. Math. Biol.}, \textbf{17} (1983), 351--369. \doilink{https://doi.org/10.1007/BF00276521}

\bibitem{RR2004}
 C.~J. Roussel, M.~R. Roussel, Reaction-diffusion models of development with state-dependent  chemical diffusion coefficients, \emph{Prog. Biophys. Mol. Biol.}, \textbf{86} (2004), 113--160. \doilink{https://doi.org/10.1016/j.pbiomolbio.2004.03.001}

\bibitem{SVN89}
 D.~Sulsky, R.~R. Vance, W.~I. Newman, Time delays in age-structured populations, \emph{J. Theor. Biol.}, \textbf{141} (1989), 403--422. \doilink{https://doi.org/10.1016/S0022-5193(89)80122-5}

\bibitem{BMN2010}
 G.~Bel, B.~Munsky, I.~Nemenman, The simplicity of completion time distributions for common complex  biochemical processes, \emph{Phys. Biol.}, \textbf{7} (2010), 016003. \doilink{https://doi.org/10.1088/1478-3975/7/1/016003}

\bibitem{Billingsley1995}
 P.~Billingsley, \emph{Probability and Measure}, Wiley, New York, 1995.

\bibitem{BarNahum-etal2005}
 G.~Bar-Nahum, V.~Epshtein, A.~E. Ruckenstein, R.~Rafikov, A.~Mustaev, E.~Nudler,
 A ratchet mechanism of transcription elongation and its control, \emph{Cell}, \textbf{120} (2005), 183--193. \doilink{https://doi.org/10.1016/j.cell.2004.11.045}

\bibitem{Shaevitz-etal2003}
 J.~W. Shaevitz, E.~A. Abbondanzieri, R.~Landick, S.~M. Block, Backtracking by single {RNA} polymerase molecules observed at  near-base-pair resolution, \emph{Nature}, \textbf{426} (2003), 684--687. \doilink{https://doi.org/10.1038/nature02191}

\bibitem{GB2000}
 M.~A. Gibson, J.~Bruck, Efficient exact stochastic simulation of chemical systems with many  species and many channels, \emph{J. Phys. Chem. A}, \textbf{104} (2000), 1876--1889.
\doilink{https://doi.org/10.1021/jp993732q}

\bibitem{BCH2013}
 H.~T. Banks, J.~Catenacci, S.~Hu, A comparison of stochastic systems with different types of delays, \emph{Stoch. Anal. Appl.}, \textbf{31} (2013), 913--955. \doilink{https://doi.org/10.1080/07362994.2013.806217}

\bibitem{FDT2016}
 Y.-L. Feng, J.-M. Dong, X.-L. Tang, Non-{Markovian} effect on gene
 transcriptional systems, \emph{Chin. Phys. Lett.}, \textbf{33} (2016),
 108701.  \doilink{https://doi.org/10.1088/0256-307X/33/10/108701}

\bibitem{LGR2012}
 J.~Lloyd-Price, A.~Gupta, A.~S. Ribeiro, SGNS2: A compartmentalized stochastic chemical kinetics simulator for  dynamic cell populations, \emph{Bioinformatics}, \textbf{28} (2012), 3004--3005.
 \doilink{https://doi.org/10.1093/bioinformatics/bts556}


\bibitem{StochPy}
 T.~Maarleveld, \emph{StochPy User Guide, Release 2.3.0}, 2015,  \urlprefix\url{https://sourceforge.net/projects/stochpy/files/stochpy_userguide_2.3.pdf/download}.


\bibitem{RL2007a}
 A.~S. Ribeiro, J.~Lloyd-Price, {SGN} {Sim}, a stochastic genetic networks simulator, \emph{Bioinformatics}, \textbf{23} (2007), 777--779.  \doilink{https://doi.org/10.1093/bioinformatics/btm004}


\bibitem{gillespie1976general}
 D.~T. Gillespie, A general method for numerically simulating the stochastic time  evolution of coupled chemical reactions, \emph{J. Comput. Phys.}, \textbf{22} (1976), 403--434. \doilink{https://doi.org/10.1016/0021-9991(76)90041-3}

\bibitem{sims2004elongation}
 R.~J. Sims~III, R.~Belotserkovskaya, D.~Reinberg, Elongation by {RNA} polymerase {II}: the short and long of it, \emph{Genes Dev.}, \textbf{18} (2004), 2437--2468.  \doilink{https://doi.org/10.1101/gad.1235904}


\bibitem{TR2020}
 E.~A.~M. Trofimenkoff, M.~R. Roussel, Small binding-site clearance delays are not negligible in gene  expression modeling, \emph{Math. Biosci.}, \textbf{325} (2020), 108376.
\doilink{https://doi.org/10.1016/j.mbs.2020.108376}

\bibitem{epshtein2003cooperation}
 V.~Epshtein, E.~Nudler, Cooperation between {RNA} polymerase molecules in transcription  elongation, \emph{Science}, \textbf{300} (2003), 801--805.  \doilink{https://doi.org/10.1126/science.1083219}


\bibitem{JWY2019}
 C.~Jia, L.~Y. Wang, G.~G. Yin, M.~Q. Zhang, Single-cell stochastic gene
 expression kinetics with coupled  positive-plus-negative feedback,
 \emph{Phys. Rev. E}, \textbf{100} (2019), 052406. \doilink{https://doi.org/10.1103/PhysRevE.100.052406}


\bibitem{SG2023}
 J.~Szavits-Nossan, R.~Grima, Uncovering the effect of {RNA} polymerase steric interactions on gene expression noise: Analytical distributions of nascent and mature {RNA}  numbers, \emph{Phys. Rev. E}, \textbf{108} (2023), 034405.  \doilink{https://doi.org/10.1103/PhysRevE.108.034405}

\bibitem{BKW+2013}
 P.~Bokes, J.~R. King, A.~T.~A. Wood, M.~Loose, Transcriptional bursting diversifies the behaviour of a toggle  switch: Hybrid simulation of stochastic gene expression, \emph{Bull. Math. Biol.}, \textbf{75} (2013), 351--371.  \doilink{https://doi.org/10.1007/s11538-013-9811-z}


\bibitem{DB2009}
 M.~Dobrzy{\'n}ski, F.~J. Bruggeman, Elongation dynamics shape bursty transcription and translation, \emph{Proc. Natl. Acad. Sci. U.S.A.}, \textbf{106} (2009), 2583--2588.\doilink{https://doi.org/10.1073/pnas.0803507106}


\bibitem{CFX2006}
 L.~Cai, N.~Friedman, X.~S. Xie, Stochastic protein expression in individual cells at the single  molecule level, \emph{Nature}, \textbf{440} (2006), 358--362.  \doilink{https://doi.org/10.1038/nature04599}



\bibitem{LJH2019}
 A.~J.~M. Larsson, P.~Johnsson, M.~Hagemann-Jensen, L.~Hartmanis,  O.~R. Faridani, B.~Reinius, et al., Genomic encoding of transcriptional burst kinetics, \emph{Nature}, \textbf{565} (2019), 251--254. \doilink{https://doi.org/10.1038/s41586-018-0836-1}

\bibitem{WAR2021}
 Y.~Wan, D.~G. Anastasakis, J.~Rodriguez, M.~Palangat, P.~Gudla,  G.~Zaki, et al., Dynamic imaging of nascent {RNA} reveals general principles of  transcription dynamics and stochastic splice site selection, \emph{Cell}, \textbf{184} (2021), 2878--2895.  \doilink{https://doi.org/10.1016/j.cell.2021.04.012}

\bibitem{JL2023}
 C.~Jia, Y.~Li, Analytical time-dependent distributions for gene expression models  with complex promoter switching mechanisms, \emph{SIAM J. Appl. Math.}, \textbf{83} (2023), 1572--1602.
\doilink{https://doi.org/10.1137/22M147219X}

\bibitem{LMB1978LLN}
 B.~W. Lindgren, G.~W. McElrath, D.~A. Berry, \emph{Introduction to Probability and Statistics}, 154--156, 4th edition, Macmillan, New York, 1978.


\bibitem{Janisch2021}
 M.~Janisch, Kolmogorov's strong law of large numbers holds for pairwise  uncorrelated random variables, \emph{Theory Probab. Appl.}, \textbf{66} (2021), 263--275.  \doilink{https://doi.org/10.4213/tvp5459}


\bibitem{Li-etal1996}
 B.~Li, J.~A. Weber, Y.~Chen, A.~L. Greenleaf, D.~S. Gilmour, Analyses of promoter-proximal pausing by {RNA} polymerase {II} on the  \emph{hsp70} heat shock gene promoter in a \emph{Drosophila} nuclear extract, \emph{Mol. Cell. Biol.}, \textbf{16} (1996), 5433--5443.
\doilink{https://doi.org/10.1128/MCB.16.10.5433}

\bibitem{Murphy2017}
 R.-J. Murphy, \emph{Stochastic Modeling of the Torpedo Mechanism of Eukaryotic  Transcription Termination}, Master's thesis, University of Lethbridge, 2017, \urlprefix\url{https://www.uleth.ca/dspace/handle/10133/4906}.



\bibitem{CCC2020}
 B.~Choi, Y.-Y. Cheng, S.~Cinar, W.~Ott, M.~R. Bennett, K.~Josi\'c, J.~K.
 Kim, Bayesian inference of distributed time delay in transcriptional and  translational regulation, \emph{Bioinformatics}, \textbf{36} (2020), 586--593. \doilink{https://doi.org/10.1093/bioinformatics/btz574}

\bibitem{HCC2023}
 H.~Hong, M.~J. Cortez, Y.-Y. Cheng, H.~J. Kim, B.~Choi, K.~Josi{\'c}, J.~K.
 Kim, Inferring delays in partially observed gene regulation processes, \emph{Bioinformatics}, \textbf{39} (2023), btad670.  \doilink{https://doi.org/10.1093/bioinformatics/btad670}


\bibitem{HS2014}
 D.~Holcman, Z.~Schuss, The narrow escape problem, \emph{SIAM Rev.}, \textbf{56} (2014), 213--257.
\doilink{https://doi.org/10.1137/120898395}


\bibitem{RT2012}
 M.~R. Roussel, T.~Tang, Simulation of {mRNA} diffusion in the nuclear environment, \emph{IET Syst. Biol.}, \textbf{6} (2012), 125--133.  \doilink{https://doi.org/10.1049/iet-syb.2011.0032}


\bibitem{Tang2010}
 S.~Tang, \emph{Mathematical Modeling of Eukaryotic Gene Expression}, PhD thesis, University of Lethbridge, 2010, available at: \url{https://www.uleth.ca/dspace/handle/10133/2567}.


\end{thebibliography}


\end{document}
