Time-resolved analysis of transcription kinetics in single live mammalian cells

In eukaryotic cells, RNA polymerase II synthesizes mRNA in three stages, initiation, elongation, and termination, and numerous factors determine how quickly a gene is transcribed to produce mRNA molecules through these steps. However, there are few techniques available to measure the rate of each step in living cells, which prevents a better understanding of transcriptional regulation. Here, we present a quantitative analysis method to extract kinetic rates of transcription from time-lapse imaging data of fluorescently labeled mRNA in live cells. Using embryonic fibroblasts cultured from two knock-in mouse models, we monitored transcription of β-actin and Arc mRNA labeled with MS2 and PP7 stem–loop systems, respectively. After inhibiting transcription initiation, we measured the elongation rate and the termination time by fitting the time trace of transcription intensity with a mathematical model function. We validated our results by comparing them with those from an autocorrelation analysis and stochastic simulations. This live-cell transcription analysis method will be useful for studying the regulation of elongation and termination steps, providing insight into the diverse mechanisms of transcriptional processes.


Introduction
Complex regulation of eukaryotic gene expression starts with transcription, through which information in DNA is copied into mRNA. At the site of transcription, RNA polymerase II (Pol II) reads a DNA code and synthesizes mRNA, and three distinct steps can be considered to comprise the transcription procedure: initiation, elongation, and termination. Initiation is the first step in which Pol II binds to a promoter sequence at the 5ʹ end of a gene [1]. In the elongation step, Pol II reads DNA, synthesizes mRNA, and moves from the 5ʹ end to the 3ʹ end [2]. Termination is the last step in which mRNA and Pol II are released from the 3ʹ end of the gene [3,4]. Thus, three kinetic rates corresponding to each step can be defined: the initiation rate c (the number of Pol II binding per unit time), the elongation rate k (the number of base pairs transcribed per unit time), and the termination time T t (the time for mRNA separated from DNA after reaching the transcription end site). Despite many efforts to understand how transcription is regulated, limitations and challenges in quantitative measurement of these three kinetic rates still remain.
In recent years, various RNA sequencing methods have been developed to study elongation and pausing of Pol II at the promoter proximal region, intron-exon junctions, and nucleosomes [5][6][7]. These methods revealed not only the nascent RNA density in the steady state but also the elongation rate in time-resolved experiments. For such timeresolved measurements, transcription was induced or inhibited to track progression of Pol II after drug treatment, whereby the progression speed of the sequence read density directly reflected the elongation rate [8][9][10][11]. Nevertheless, population studies cannot resolve the transcriptional processes occurring at a single gene locus and lack the temporal resolution to investigate fast dynamics in real time.
Live-cell imaging has been applied to study transcriptional dynamics in single cells by labeling mRNA, Pol II, or transcription factors with fluorescent proteins [12,13]. Fluorescence recovery after photobleaching (FRAP) technique has been used to measure transcriptional kinetic rates by fitting recovery curves with various model functions: half recovery time [14], linear and exponential fitting [15][16][17], and solutions of rate equations [18][19][20][21]. However, there is concern about the accuracy of kinetic modeling using FRAP because recovery curves can be affected by photobleaching, diffusion, and binding of fluorescent particles; additionally, different kinetic models can fit the same experimental data almost equally well [22]. More recently, fluctuation correlation analysis of fluorescently labeled mRNA was used to measure the transcription initiation rate and dwell time of a transcript of a yeast gene [23]. In this analysis, the dwell time consists of the elongation and termination times, which cannot be separated in the steady state. Nevertheless, the authors were able to calculate the elongation rate and the termination time by comparing results from two different constructs bearing a PP7 RNA stem-loop cassette in either the 5ʹ untranslated region (UTR) or 3ʹ UTR [23]. Likewise, another study measured the elongation rate of a reporter gene in living Drosophila embryos by comparing the results from two different positions of an MS2 stem-loop cassette [24]. After that, Coulon et al. developed a dual-color fluctuation correlation analysis for human β-globin reporter mRNA labeled with a PP7 cassette in the second intron and an MS2 cassette in the 3ʹ UTR [25]. Using two autocorrelation functions and a single crosscorrelation function, the authors measured the kinetics of transcription and splicing of single mRNAs. Liu et al. also used a two-color MS2/PP7 reporter to infer the initiation rate, elongation rate and termination time of the hunchback gene in developing fruit fly embryos [26]. However, it would be challenging to insert~1.5-kblong MS2 and PP7 cassettes into two different locations in an endogenous gene, which limits the general application of the technique in mammalian cells and tissues.
In this study, we developed a time-resolved transcriptional analysis method to measure the elongation rate and termination time of an endogenous mRNA labeled with one stem-loop cassette. Previously, two knock-in (KI) mouse models, Actb-MBS [27] and Arc-PBS [28], were generated to label β-actin (Actb) and Arc mRNA, respectively. In the former, 24 repeats of the MS2 binding site (MBS) were inserted in the 3′ UTR of the Actb gene [27] ( Figure 1A); in the latter, 24 repeats of the PP7 binding site (PBS) were knocked in the 3′ UTR of the Arc gene [28] ( Figure 1B). Using mouse embryonic fibroblasts (MEFs) derived from these two mouse models, we investigated the transcriptional dynamics of the Actb and Arc genes. First, by applying our time-resolved analysis method, we measured the elongation rate and the termination time of each transcription allele. Then we measured the initiation rate and the total dwell time using the autocorrelation analysis method [23] and compared the total dwell time with the time-resolved analysis results. We found a highly heterogeneous distribution of the kinetic parameters in each locus of both Actb and Arc genes. Yet, there was a significant difference between Actb and Arc regarding the average transcriptional kinetic rates. This approach allows quantitative measurement of all three kinetic rates of single-color labeled mRNA transcription in live cells, providing opportunities to study diverse transcriptional regulation at the single-allele level with high spatiotemporal resolution.

Time-resolved transcription analysis model
We developed a time-resolved transcription analysis model to determine the transcription elongation rate k and the termination time T t . The dynamics of translation have previously been analyzed using the SunTag system [29] with a mathematical model describing the decreasing fluorescence intensity of a translation site after inhibiting translation initiation [30]. Because RNA stem-loop systems are analogous to the SunTag system, we adapted this model for application to transcriptional dynamics. While translation termination time was not considered in the previous model, we included transcription termination time in our model and derived an analytical function for the decreasing fluorescence intensity after inhibition of transcription initiation.
We assumed that Pol IIs are distributed uniformly on the gene and move toward the 3ʹ end at a constant speed k. Although the elongation process can be complex with pausing, traffic jams, backtracking, etc.
[31], we sought to measure the average elongation rate k over the entire length of a gene. For simplicity, we modeled transcription termination as a firstorder irreversible process with a rate constant a. Hence, the number of Pol IIs residing at the 3ʹ end of the gene at time t, N(t) can be calculated as follows.
The transcription termination time was defined as T t 1/a. Suppose that there is no Pol II at the 3ʹ end and that the first Pol II Frontiers in Physics frontiersin.org 02 just reaches the end at time t 0. With the transcription initiation rate c, the number of Pol II reaching the 3ʹ end of the gene during the first time interval Δt is cΔt. After another time interval Δt, the number of Pol II accumulated at the 3ʹ end becomes If we assume that Pol II and nascent RNAs fall off at the same time from the 3ʹ end of the gene, Eq. 2 also describes the number of nascent RNAs accumulated at the 3ʹ end at time 2Δt. After an infinite number of Δt intervals, the number of nascent RNAs at the 3ʹ end, N(∞, Δt) becomes In the limit that the time interval Δt goes to zero, This result shows that cT t nascent RNAs accumulate at the 3ʹ end in the steady state. Now let us suppose transcription initiation is blocked at time t 0 while in the steady state. After time t, the number of nascent RNAs in the 3ʹ end becomes Based on the previous study of translational dynamics [30], we set the nascent RNA density as a function of the position on the gene and the time after inhibitor addition. Following the formula for the ribosome probability density in the previous Frontiers in Physics frontiersin.org study [30], we first modeled the nascent RNA density function R(x, t) using a Heaviside step function: Because k is the elongation rate of Pol II, c/k indicates the number of nascent RNAs per unit length in the Pol II occupied region. The density function in Eq. 6 propagates to the 3′ end direction over time. In this model, however, it is not considered that the nascent RNAs are accumulated at the 3′ end (x L). We expressed the nascent RNA accumulation using a Dirac delta function: In Eq. 7, the summation of two Heaviside functions H(L/k − t) and H(t − L/k)e −(t−L/k)/Tt is multiplied by the Dirac delta function. Before the last Pol II enters the 3′ end (t < L/k), the amplitude of Dirac delta function (the number of accumulated nascent RNAs) is steady. After the last Pol II enters the 3′ end (t ≥ L/k), the amplitude decreases exponentially as shown in Eq. 5. Thus, the nascent RNA density at position x and time t, R(x, t) can be described by adding Eq 7 to Eq 6: This density function is depicted in Figure 2A. The first term in Eq. 8 is the nascent RNA density without the termination time effect. In the second term, the amplitude of the Dirac delta function means the number of accumulated nascent RNAs at position L.
When Pol II starts to transcribe the stem-loop region [X 1 , X 2 ], the fluorescence intensity of a single nascent RNA I(x) increases as more stem-loops are synthesized. In the post stem-loop region [X 2 , L], I(x) remains constant ( Figures 1A,B). Then, the intensity with the Pol II position x becomes I mRNA is the fluorescence intensity of a single mRNA that is fully transcribed (all 24 stem-loops). Finally, the fluorescence intensity of a transcription site at time t, F(t) can be derived The fluorescence intensity at time t becomes In the Eq. 11, the exponential term in the fourth interval is due to the termination time effect. A normalized fluorescence intensity This is the model function of time-resolved transcription analysis that can be used for fitting experimental data. With the three fitting parameters, F(0), k and T t , the elongation and terminationkineticscanbedistinguishedina single-cellexperiment.

Autocorrelation analysis model
Larson et al. introduced a mathematical model to analyze transcriptional dynamics using an autocorrelation function of a transcription site intensity trace [23]. We rederived the autocorrelation function in a more concise manner and applied the model to Actb-MBS ( Figure 1A) and Arc-PBS ( Figure 1B) mRNAs. The autocorrelation function is defined as where δI(t) I(t) − 〈I(t)〉 and I(t) is the intensity of a transcription site. The bracket denotes the average intensity over time. The experimental data can be fit with the following mathematical model: where N and M are the length of the stem-loop [X 1 , X 2 ] and the post stem-loop [X 2 , L] region scaled with the stem-loop size l s . k s is the effective transition rate of Pol II. T d is the total dwell time, which includes the elongation time through the fluorescent region (the stem-loop and the post stem-loop region) and the termination time. I n is the fluorescence intensity of a pre-mRNA after Pol II transcribes the n th stem-loop from the 5ʹ end of the stem-loop region. This pre-mRNA intensity is described as follows: P(n, t + τ|m, t) in Eq. 14 is the conditional probability that a Pol II is at position n at time (t + τ) given that it has been at position m at time t.
Larson et al. approximated the autocorrelation function in Eq. 14 to consider two special cases in which the position of the stem-loop region was at either the 5ʹ or 3ʹ end of an mRNA. Such approximation cannot be applied for Actb-MBS and Arc-PBS mRNAs; the full autocorrelation function in Eq. 14 needs to be used for a general position of the stem-loop region. The autocorrelation function was calculated based on three assumptions: 1) initiation occurs randomly and evenly; 2) Pol II elongates only in one direction from the 5ʹ end to the 3ʹ end, and there is no reverse translocation; and 3) pre-mRNA is released immediately at the 3ʹ end of the gene.
With the assumption that Pol II can progress only in one direction, the probability density of Pol II dwell time t n in an arbitrary position n can be calculated for the irreversible process.
The time spent for the transition to the next position is exponentially distributed. If t i is the time spent for the transition from the i th stem-loop to the (i + 1) th stem-loop, the conditional probability P(n, t + τ|m, t) becomes while satisfying constraints given as The (n − m) dimensional integral denoted as d n−m t is dt m /dt n−1 . Each exponential term represents the probability density that Pol II spends time t i at position i. Because Pol II is in the n th stem-loop at time (t + τ), if Pol II spends time τ′ from the m th stem-loop to the (n − 1) th stem-loop, the possible dwell time in the n th stem-loop is from τ − τ′ to ∞. Then, the conditional probability is described as By integrating with dt n first, The last integral term is the (n − m) dimensional volume in the domain denoted in Eq. 19. This domain is a (n − m)-dimensional hyper pyramid in which (n − m) segments are perpendicular to each other, and their lengths are τ. The last integral in Eq. 21 is τ n−m /(n − m)!. Finally, the conditional probability becomes This is the Poisson distribution; k s τ is the expected number of stem-loops transcribed during time τ. The autocorrelation function can be calculated by substituting Eq. 22 into Eq. 14. The initiation rate c and the effective transition rate k s can be obtained by fitting the autocorrelation curve of the experimental data with Eq. 14. Using the k s value and Eq. 15, we can calculate the total dwell time T d , which includes the elongation time through the fluorescent region (the stem-loop and the post stem-loop region) and the termination time T t . Although we assumed in our derivation that pre-mRNA is released immediately at the 3ʹ end of the gene, the total dwell time calculated from k s effectively includes the termination time.

Fluorescence microscopy
All images were taken using a wide-field fluorescence microscopy system based on an IX-83 inverted microscope (Olympus) equipped with a UAPON 150×/1.45 NA oil objective (Olympus), iXon Life 888 electron-multiplying charge-coupled device (EMCCD) camera (Andor), SOLA SE u-niR light engine (Lumencor), Chamlide TC top-stage incubator system (Live Cell Instrument), and ET-EGFP filter set (Chroma, ET470/40x, T495lpxr, ET525/50 m). Live-cell imaging was performed at 37°C, and time-lapse z-stack images were taken at an interval of every 8 s over a period of 67 min. Z-stacks were imaged at an interval of 0.5 μm for a total height of 6 μm. The imaging was started right after the inhibitor addition for the time-resolved measurement and serum addition for the auto-correlation analysis (Supplementary Figure S1).

Image analysis
All z-stack images were maximum projected and bleachcorrected via an exponential bleach correction fit using Fiji [34]. Fluorescence intensity time traces of TSs were generated by using HybTrack [35]. The background level was determined by the average intensity of the period in which TS was not automatically detected by HybTrack. For each time trace, transcription 'ON' states were distinguished from 'OFF' states by two criteria: 1) the TS intensity was higher than 1.5 times the background level, and 2) the durations of both ON and OFF states were equal to or longer than 1 min.
For autocorrelation analysis, we calculated the autocorrelation function of the time traces in the ON-state; time traces longer than the average ON-time were used for analysis. The autocorrelation function was calculated by the multitau algorithm [36]. This algorithm reduces noise in the correlation function at a long lag time. To measure the initiation rate and the total dwell time, we fitted an autocorrelation Frontiers in Physics frontiersin.org function of a single ON-time trace with Eq. 14. We used weighted least squares fitting with weights equal to the inverse standard error of the mean. The fitting parameters were the initiation rate c and the effective transition rate of Pol II k s . For time-resolved measurements after inhibiting transcription initiation, we analyzed time traces showing a plateau followed by a decreasing signal. In our analysis, we assumed that the gene is fully and evenly occupied by Pol IIs before inhibition of transcription initiation; thus, the fluorescence intensity of TS cannot increase after inhibition. Nonetheless, some data exhibited increasing fluorescence intensity after treatment with triptolide; those time traces were considered not to satisfy the condition of our kinetic model. The selected intensity time traces were fit with the model function with a nonlinear least squares method. The fitting parameters were the elongation rate k, the termination time T t , and the normalization factor F(0).

Time-resolved transcription analysis
To label Actb and Arc mRNA with green fluorescent proteins (GFPs), we used the MS2-and PP7-GFP systems, respectively. The MS2-GFP system utilizes highly specific binding between the MS2 bacteriophage capsid protein (MCP) and the MBS RNA stem-loop [37]. Constitutive expression of MCP fused with GFP (MCP-GFP) labeled all endogenous β-actin mRNA in mouse embryonic fibroblasts (MEFs) cultured from Actb-MBS KI mice in which 24 repeats of the MBS were inserted in the 3′ UTR of the β-actin gene [27] (Figure 1A). In a similar manner, endogenous Arc mRNA was visualized by expressing PP7 capsid protein (PCP) fused with GFP (PCP-GFP) in MEFs derived from Arc-PBS KI mice bearing 24 repeats of PBS in the Arc gene 3′ UTR [28] ( Figure 1B). In Figures 1A,B, the positions of the start and end of the stem-loop cassettes are denoted as X 1 and X 2 , and the intensity change of a precursor mRNA (pre-mRNA) is shown as a function of the Pol II position x. Detailed information of the stem-loop position is described in Table 1. The MS2 or PP7 capsid protein fused with GFP (CP-GFP) complexes carry a nuclear localization sequence (NLS), which targets the protein to the nucleus. When pre-mRNAs are transcribed, multiple CP-GFPs bind to the RNA binding sites (Figure 1C), and the transcription site appears as a bright spot in fluorescence images ( Figure 1D, Supplementary Videos S1, S2).
To induce transcription, we starved homozygous Actb-MBS and Arc-PBS MEFs overnight and added serum-containing medium. A few minutes after serum induction, up to four transcription sites appeared in nuclei ( Figure 1D). We performed time-resolved experiments by using a small molecule transcription inhibitor, triptolide (Trp) (MW 360.6). It has been shown that triptolide inhibits transcription initiation by preventing transcription bubble formation without significantly affecting the elongation rate [38,39]. After inhibiting transcription initiation, the region occupied by Pol IIs decreased over time in the direction of elongation ( Figure 2A). Then, the normalized intensity of a transcription site declined in 4 steps ( Figure 2B). First, the intensity remained constant in 0 ≤ t < X 1 /k after inhibitor addition. As the last Pol II entered the stem-loop region, the intensity decreased quadratically in the interval of X 1 /k ≤ t < X 2 /k. After that, the intensity decreased linearly in the time interval of X 2 /k ≤ t < L/k as the last Pol II entered the post stem-loop region. After time t L/k, there were no elongating Pol IIs, and the intensity decreased exponentially due to transcription termination and release of Pol II and mRNA. The model function behaviors for the different values of the termination time T t and the elongation rate k are shown in Figures 2C,D, respectively. The increase in the termination time T t resulted in a longer decay time while the time length of the plateau remained the same ( Figure 2C). The length of the plateau X 1 /k is solely determined by the elongation rate k ( Figure 2D), which denotes the speed of Pol II in the pre stem-loop region. Therefore, the elongation rate k and termination time T t can be acquired by fitting the normalized transcription site intensity of the experimental data with the model function in Eq. 12.
Immediately after adding the initiation inhibitor, we imaged transcription sites in living cells by using fluorescence microscopy. Representative data of the normalized fluorescence intensities of transcription sites in Actb-MBS and Arc-PP7 MEFs are shown in Figures 3A,B, respectively. By fitting the data with Eq. 12, we obtained an average elongation rate of Actb transcription to be 10.4 ± 2.9 nt/s, which was similar to the elongation rate measured by fluorescence in situ hybridization (FISH) in normal rat kidney (NRK) cells [40]. On the other hand, the average elongation rate of Arc was significantly different from that of Actb; the average elongation rate of Arc transcription was 29.0 ± 17.0 nt/s, almost 3-fold faster than that of Actb  Figure 3C). Termination times also differed significantly. As illustrated in Figure 3D, the average termination time of Actb mRNA (233 ± 59 s) was 2.7-fold longer than that of Arc mRNA (86 ± 27 s).
To validate our measurements, we used another transcription inhibitor, flavopiridol (FP) (MW 438.30), which is known to inhibit cyclin-dependent kinase 9 (Cdk9) in the positive transcription elongation factor b (P-TEFb) complex [41]. This P-TEFb inhibitor FP prohibits Pol II release from the promoter proximal pausing to the productive elongation phase, whereas Trp inhibits new Pol II initiation [8]. Although the mechanisms of inhibition by Trp and FP are different, they result in a similar physical state that inhibits transcription initiation. Using FP with the same procedure as in the Trp experiment, we obtained elongation rate and termination time consistent with the Trp results. We measured the elongation rates of Actb and Arc transcription from the FP experiment as 13.3 ± 5.6 and 29.0 ± 15.3 nt/s, respectively, similar to those from the Trp experiment ( Figures 4A,B). The termination times obtained from the FP experiments were 223 ± 59 and 60 ± 25 s for Actb and Arc transcription, respectively. These values were also consistent with the Trp results ( Figures 4C,D). The consistency between the two independent experiments using Trp and FP supports that our time-resolved analysis approach can measure transcriptional elongation and termination kinetics without distinct toxicity effect by the drugs. We determined the average elongation rate and termination time of each gene by pooling the singleallele data set from the Trp and FP experiments (Table 1).

Autocorrelation analysis
Next, we examined whether the time-resolved results are valid by using the autocorrelation analysis described in the section of mathematical models. For the autocorrelation analysis, we recorded the intensity of transcription sites, which showed bursting transcriptional activity without the initiation inhibitors after serum induction ( Figures 5A,B). Both the Actb and Arc genes exhibited stochastic switching between active (ON) and inactive (OFF) states after serum induction. The average duration of the ON-state was 15.4 ± 6.6 min (means ± SD) for Actb, which was significantly longer Frontiers in Physics frontiersin.org than the 8.6 ± 3.4 min observed for Arc (p < 0.01, t test; Figure 5C). The duration of the OFF state was 6.5 ± 3.4 min for Actb, shorter than the 9.4 ± 4.0 min for Arc (p < 0.01, t test; Figure 5D). These results demonstrate that Actb exhibits transcriptional bursts more frequently in longer time periods than Arc does. The initiation rate and total dwell time were measured from the autocorrelation function of transcription intensity traces during ON-times. Because the autocorrelation analysis model was developed for a steady-state condition, we calculated the autocorrelation function from the selected ON-state traces of which durations were longer than the average ON-time of each gene. Then, the selected ON-time duration was 28.9 ± 9.9 min for Actb (183 traces) and 15.5 ± 7.3 min for Arc (82 traces) on average. Before the autocorrelation analysis, we characterized the general behavior of the model function shown in Eq. 14. The autocorrelation amplitude G(0) in Eq. 14 decreases as either the initiation rate c or the total dwell time T d increases ( Figures 6A,B). The decay time of the autocorrelation function is independent of c ( Figure 6A) but increases as T d increases ( Figure 6B). Representative autocorrelation curves of Actb and Arc transcription are depicted in Figures 6C,D, respectively. We plotted c and T d values obtained from individual ON-traces in Figures 6E,F, respectively. The average initiation rate and total dwell time were c 0.12 ± 0.10 s −1 and T d 348 ± 104 s for Actb and c 0.21 ± 0.17 s −1 and T d 179 ± 38 s for Arc. These average total dwell times T d were about 20% of the average of the subjected ON-time duration.
Next, we compared the total dwell times from the timeresolved measurements with those from the autocorrelation analysis. Using the results of the time-resolved measurements, the total dwell time can be calculated as.
With this, the total dwell times from the time-resolved measurements were calculated as T TR d 357 ± 77 s for Actb Frontiers in Physics frontiersin.org and 179 ± 66 s for Arc transcription. The two different dwell times, T d and T TR d , obtained from the autocorrelation and the time-resolved analyses, respectively, were quite similar for each gene. We note that the total dwell time was defined as the nascent RNA dwell time in the fluorescence region (stemloop and post stem-loop regions) in the autocorrelation analysis. In the time-resolved analysis, however, the total dwell time was calculated using the elongation rate k which Frontiers in Physics frontiersin.org was measured in the region before the stem-loops. Therefore, we compared T d and T TR d under the assumption of a constant Pol II elongation rate k over the whole gene. Although we made such a simple assumption, T d and T TR d matched reasonably well.

Error estimation by simulations
Using simulations, we estimated the errors in the timeresolved measurements for each gene structure of Actb and Arc labeled with the stem-loops. The input parameters of the Frontiers in Physics frontiersin.org simulations were initiation rate c, elongation rate k, and termination time T t . For our analytical model, we made a simple assumption of uniform Pol II distribution with a constant initiation rate c. However, it is known that the initiation time interval between two Pol IIs is dependent on the promoter dynamics [31,42]. For a Poissonian promoter, the polymerase initiation time interval follows an exponential distribution. Therefore, we initiated Pol IIs with an exponential time interval T int in computer simulations: The mean time interval is 1/c as in our time-resolved analytical model.
The elongation process follows the Poisson distribution as shown in Eq. 22. Then, the probability of x nucleotides elongated in time Δt is In the simulation, the time Δt was 8 s, which was the same as the interval of our time-lapse imaging data. Eq. 25 shows that the number of nucleotides elongated by Pol II during Δt follows a Poisson distribution with a mean of kΔt. As Pol II transcribes a stem loop, the intensity of the transcription site increases by a unit brightness. After summation of all pre-mRNA intensities, the simulation generated the intensity time trace of a transcription site. For the simulation of the time-resolved analysis, Pol IIs were distributed at random positions on the gene considering the Poissonian initiation and elongation, and cT t Pol IIs were placed at position L at time t 0. The termination time was given as an exponential random value, and there was no newly initiated Pol II after t 0. Figure 7 shows the average normalized simulation intensity of Actb 1) and Arc 2) transcription sites (blue dots) when transcription initiation is inhibited at time t 0. These intensity time traces were generated by repeating the simulation 100 times with input parameters from the experimental results (Actb: c 0.12 s −1 , k 12.0 nt/s, and T t 228 s. Arc: c 0.21 s −1 , k 29.0 nt/s, and T t 179 s). We fitted the single intensity time traces with the model function from Eq. 12. The red lines in Figures 7A,B show the model function with the average fitted values. The average elongation rate and termination time obtained by fitting the simulation curves were k 12.4 ± 1.9 nt/s and T t 226 ± 45 s for Actb; k 30.6 ± 13.2 nt/s and T t 174 ± 28 s for Arc. These fitting results are quite similar to their ground truth values used as the input parameters in the simulations. Figure 8 shows the errors in the elongation rate k and the termination time T t obtained by fitting the simulation results corresponding to various input parameters. The input initiation rate c varied from 0.01 s −1 to 1 s −1 in the x-axis. The input elongation rate k varied from 10 0.5 to 10 2 nt/s in the y-axis. The input termination time T t is shown on top in Figure 8. The error in the output value was defined as |input − output|/input. The color maps in Figure 8 display the errors in the output values of elongation rate k and total dwell time T t . The mean errors are generally larger for shorter termination times T t and lower initiation rates c. Because the assumption of uniformly distributed Pol IIs is not valid for the low c limits, the timeresolved analysis is not appropriate in this realm.
The areas inside the magenta lines in Figure 8 indicate the regime in which the mean errors are under 10%. The valid region of the time-resolved model is on the lower and right sides of the color map. The time-resolved analysis of our Actb and Arc transcription data yielded the elongation rate k and the termination time T t within the valid region of less than 10% error. These results support the overall validity of the elongation rate and termination time measured by the time-resolved method.

Discussion
In this study, we investigated the transcription kinetics of the Actb and Arc genes by using two analysis methods of live-cell imaging. Because transcription initiation, elongation, and termination dynamic factors cannot be measured separately in the steady state, we used Trp and FP to inhibit transcription initiation and performed time-resolved measurements. For simplicity, we assumed that Pol IIs accumulate at the end of the gene for transcription termination and described it as a Dirac delta function in the Pol II density distribution R(x, t). By fitting our data with this mathematical model, we obtained the elongation rate and the termination time separately and confirmed that the total dwell times of Actb and Arc nascent RNAs were consistent with those from the autocorrelation measurement. To identify the regime where the model function fittings are valid, we checked the error behavior of the time-resolved model with various input parameters. The elongation rates and the termination times of Actb and Arc measured in this study were well within the valid region of our analysis method.
We obtained significantly different elongation rates and termination times for Actb and Arc transcription. The elongation process can include different events such as pausing, traffic jams, backtracking, and so forth [31]. Therefore, the elongation rate obtained in this study represents the overall average elongation rate affected by those complex processes. The Actb and Arc elongation rates were 12.0 ± 4.6 bp/s and 29.0 ± 15.5, respectively. So, it appears that Pol II elongation was interrupted more often on the Actb gene than on the Arc. There are many factors that affect the elongation rate of Pol II, such as exon density, nucleosome structure, and CpG content [8,11]. Previous studies have revealed that Pol II pauses or slows down at intron-exon junctions [8,9,11] and that exon density has the greatest effect on the elongation rate among the factors mentioned above [8]. The gene lengths of Actb and Arc are similar (Actb: 3,640 bp, Arc: 3,490 bp), but Actb has 6 exons and Arc has 3. Thus, the exon density of Actb (1.6 exons/kb) is about two times higher than that of Frontiers in Physics frontiersin.org  Frontiers in Physics frontiersin.org 13 Arc (0.9 exons/kb), which could be the reason why the elongation rate of Actb is much slower than that of Arc.
Furthermore, the termination time of Actb (228 ± 59 s) was 3.5-fold longer than that of Arc (65 ± 25 s). There are several different types of transcription termination in mammalian cells [4], and the termination of human Actb transcription occurs by Pol II pause-dependent termination type. This Pol II pausing occurs at the G-rich pause sequence downstream of the poly(A) signal (PAS) [43]. In the human Actb gene, nascent transcripts hybridize to the antisense DNA strand and displace the sense DNA strand, forming a structure called an R-loop [44]. The R-loop is formed prevalently over G-rich pause sites behind the elongation complex and induces Pol II pausing. Senataxin, a DNA-RNA helicase, resolves R-loops to allow the 5ʹ-3ʹ exonuclease Xrn2 to degrade the nascent RNA fragment from the poly(A) cleavage site [44]. Similar to human Actb, mouse Actb also has a G-rich region at 450 nucleotides downstream of the PAS; hence, transcription of mouse Actb is also likely to be terminated by the R-loop-dependent pausing mechanism. We determined the termination time, which was defined as the time from reaching the transcription end site (TES) to the release of mRNA, to be approximately 230 s for mouse Actb. This data provides us with information on the timescale of transcription termination by R-loop-dependent pausing. In contrast, the mouse Arc gene does not have a G-rich region downstream of the PAS and exhibits a much shorter transcription termination time of~65 s. Therefore, we presume that Arc transcription termination occurs via a different mechanism. The Arc gene has a viral origin and contains sequences that are related to retrotransposon Gag genes [45,46]. It is also interesting that the transcription termination time of Arc is similar to that of an HIV-1 RNA reporter gene (63.5 s) in U2OS cells [15]. Our new analysis technique to measure termination time can be useful to unravel diverse transcription termination processes.
We note that our time-resolved model has some limitations due to the simple assumption of uniform transcription initiation and elongation rates. In fact, the initiation and elongation processes are more complex, depending on the promoter dynamics [31,42]. For instance, the inter-polymerase distance distribution is exponential for a Poissonian promoter and is much more complicated for a bursty promoter. Although Actb and Arc genes exhibit ON and OFF transitions as shown in Figure 5, we investigated the transcriptional activity only during the ON states. Because the average duration of the ON states was sufficiently long, on the order of 10 min for both genes, we were able to analyze the transcription dynamics within a single ON state. The average number of Pol IIs on a gene can be estimated as c(L/k) + cT t . Using the data summarized in Table 1, we estimated the average number of Pol IIs on Actb and Arc genes to be 77 and 50, respectively. These numbers are sufficiently large to approximate the Poissonian initiation as a uniform procedure. Therefore, our simple assumption of the uniform Pol II distribution did not result in significant errors for these genes as demonstrated in Figure 8, where we compared our results from the time-resolved analysis with simulation results. However, our simple analysis method may not be generally applicable for a Poissonian promoter with a smaller number of Pol IIs or a bursty promoter with a much faster switching behavior.
The transcription initiation rate was not determined by our timeresolved model because we normalized the decreasing intensity by the plateau level F(0). If the intensity of a single mRNA I mRNA can be measured, we can calculate the initiation rate c as c F(0) I mRNA 2L−X1−X2 2k Because of the high GFP background in the nucleus, it was difficult to measure the intensities of single mRNA and transcription sites simultaneously in our live-cell experiments. We expect that background-free RNA imaging techniques such as the MS2-PP7 hybrid stem-loop system combined with split-GFP approaches [47,48] will allow for measuring both single mRNA and transcription site intensities and extracting all three transcription kinetic rates, c, k, and T t , simultaneously in a single cell.
In addition, we expect this analysis method to be used in conjunction with RNA-targeting CRISPR-Cas systems [49][50][51], through which we could analyze the transcription kinetics of any endogenous RNA without requiring genetic manipulation. Our timeresolved model function was calculated by integrating the product of the Pol II density function and the fluorescence intensity of a nascent RNA. Such model function for a different RNA tagging system can be easily derived following the same steps described herein. This analysis technique will expand our capability to investigate the transcription kinetics of endogenous genes in a highly quantitative manner.

Author contributions
HP designed research; HC developed the mathematical model and performed all experiments, analysis, and simulations; HC and BL developed software; and HC and HP wrote the paper.

Funding
This work was supported by the Creative-Pioneering Researchers Program through Seoul National University, Howard Frontiers in Physics frontiersin.org