
==== Front
Research (Wash D C)
Research (Wash D C)
RESEARCH
Research
2639-5274
AAAS

10.34133/research.0253
0253
Research Article
Regulatory Mechanisms for Transcriptional Bursting Revealed by an Event-Based Model
Wu Renjie 1
Zhou Bangyan 1
Wang Wei 1 2 *
https://orcid.org/0000-0001-9560-8518
Liu Feng 1 2 *
1 National Laboratory of Solid State Microstructures, Department of Physics, and Collaborative Innovation Center of Advanced Microstructures, Nanjing University, Nanjing 210093, P. R. China.
2 Institute for Brain Sciences, Nanjing University, Nanjing 210093, P. R. China.
* Address correspondence to: wangwei@nju.edu.cn (W.W.); fliu@nju.edu.cn (F.L.)
24 10 2023
2023
6 025325 5 2023
01 10 2023
Copyright © 2023 Renjie Wu et al.
2023
Renjie Wu et al.
https://creativecommons.org/licenses/by/4.0/ Exclusive licensee Science and Technology Review Publishing House. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution License 4.0 (CC BY 4.0).

Gene transcription often occurs in discrete bursts, and it can be difficult to deduce the underlying regulatory mechanisms for transcriptional bursting with limited experimental data. Here, we categorize numerous states of single eukaryotic genes and identify 6 essential transcriptional events, each comprising a series of state transitions; transcriptional bursting is characterized as a sequence of 4 events, capable of being organized in various configurations, in addition to the beginning and ending events. By associating transcriptional kinetics with mean durations and recurrence probabilities of the events, we unravel how transcriptional bursting is modulated by various regulators including transcription factors. Through analytical derivation and numerical simulation, this study reveals key state transitions contributing to transcriptional sensitivity and specificity, typical characteristics of burst profiles, global constraints on intrinsic transcriptional noise, major regulatory modes in individual genes and across the genome, and requirements for fast gene induction upon stimulation. It is illustrated how biochemical reactions on different time scales are modulated to separately shape the durations and ordering of the events. Our results suggest that transcriptional patterns are essentially controlled by a shared set of transcriptional events occurring under specific promoter architectures and regulatory modes, the number of which is actually limited.
==== Body
pmcIntroduction

Gene transcription is a complicated dynamic process integrating regulatory signals with genetic information. mRNAs can be produced at constant rates or in short intervals followed by long periods of inactivity, termed transcriptional bursting [1–3]. Although burst traces are qualitatively similar across organisms, the duration of the burst cycle varies remarkably, ranging from minutes to hours [4]. Does this great variability reflect diverse molecular mechanisms or stem from a shared set of highly adjustable molecular events? Similarly, molecular processes involved in transcription typically last from milliseconds to hours [3–5]. How can they be affected by regulators to underlie transcriptional bursting? How can bursting kinetics be modulated across such a wide range of time scales? These topics are fundamental to comprehending the regulation of gene transcription.

Live-cell fluorescence measurements have yielded large amounts of single-cell data [3,6,7]. Durations of transcriptionally active and inactive states and burst size were extracted to characterize transcriptional bursting and expression noise. To theoretically interpret them, phenomenological models have been built upon independent effective processes [8–14], enabling the association between regulatory signals and experimental observables. By incorporating more transcriptional details such as multistep initiation and multistep degradation of mRNA, newly developed models gain more insights into regulatory mechanisms [15–22]. To unravel gene regulation and promoter configurations through the lens of input–output relationships, however, precise knowledge of biochemical reaction mechanisms is essential; yet, such knowledge is often unavailable and must be inferred on the basis of assumptions. This could result in discrepancy between model predictions and experimental data, ignorance of global constraint on transcriptional noise, and inaccurate mapping between rate-limiting steps and effective processes under limiting conditions when deducing the underlying regulatory mechanism [23,24]. A different modeling framework is required.

On the basis of a set of common biochemical reactions essential to mRNA synthesis, a generic model has been developed, which is compatible with those simplified models and fills in their gaps [25–27]. Nevertheless, there always exist additional gene-specific reactions; some common reactions actually comprise multiple steps, and the transition process often occurs through distinct paths, leaving a memory in the transcriptional path [9,28]. While these details could be incorporated to construct more intricate models, we may develop a simpler and more computationally efficient method using queuing theory, which is widely leveraged to understand the behavior of systems with queues. Indeed, queuing models [28–30] were proposed to analyze empirical data on mRNA copy numbers and waiting-time distributions of biochemical reactions while neglecting the underlying molecular processes. To elucidate the molecular mechanism of gene regulation, we may focus on essential transcriptional events, with each comprising multiple molecular processes, and simplify the transcription into a sequence of events.

Gene transcription is modulated by various regulatory factors, such as transcription factors (TFs) that bind DNA to activate transcription and modifiers and remodelers that induce changes to histone markers and nucleosome arrays to overcome structural barriers [31]. They may exert a wide influence on transcription by altering the duration and directionality of molecular processes [32–34]. To unravel the roles for regulators in gene expression, it is essential to differentiate their impacts on rate-limiting steps and probability flows of event transitions in transcription cycle.

Here, we first propose a network of first-order Markovian biochemical reactions. Because transcriptional bursting is composed of inactive and active phases involving 3 major stages, i.e., chromatin opening for promoter access, assembly of the preinitiation complex (PIC), and mRNA production, all gene states can be divided into 5 sets to separate the stages, with functional state transitions classified into 9 categories. A series of state transitions is combined into an event, and a 6-event model is developed for transcriptional bursting. This model is generic and requires no free parameters. Transcriptional commonality is reflected in the sequential order of the 4 events, while gene specificity is embodied in the duration distributions and ordering of the events. The kinetic features of transcription, such as the mean mRNA number, burst size, burst frequency, and duration of the (in)active phase, are interconnected and functions of 6 model parameters depending on the regulator concentration (i.e., 4 average event durations and 2 repetition probabilities).

On the basis of this event model, we clarify typical features of and intrinsic constraints on transcriptional bursting, illuminate how the regulatory mechanism can be inferred by fitting experimental data at the single-gene and genome levels, and explore the transient transcriptional dynamics upon stimulation. Both numerical and analytical results are presented to provide an integrated view of how reactions on different time scales are modulated to underpin specific burst profiles. The pivotal regulatory modes and their functional implications are elaborated in detail.

Results

Establishment of the event model

Although gene transcription exhibits heterogeneity, its common features allow for the development of a generic model for transcriptional dynamics. To reveal the role for regulatory factors in transcriptional modulation, we need establish a model based on transcriptional events, which is both accurate and simple enough for theoretical analysis.

We first proposed a network model of inducible eukaryotic transcription (Fig. 1A). The transcription begins with activators binding to the promoter/enhancer [35], which promotes the recruitment of cofactors, chromatin remodeling [34,36] and histone modifications [31,37]. General TFs (GTFs; including TFIIA, TFIIB, TFIID, TFIIE, TFIIF, and TFIIH) and RNA polymerase II (Pol II) are then recruited to the core promoter, forming the PIC. The enhancer-bound activator and the PIC are connected by the mediator complex [38,39] via enhancer-promoter proximity. The PIC turns into the open complex after the DNA template strand is positioned into the active center cleft of Pol II. With the preparations completed, Pol II gets away into elongation to synthesize nucleotide chains [40], while the remaining scaffold complex (composed of activators, TFIIA, TFIID, TFIIE, TFIIH, and mediator) on the promoter facilitates transcriptional reinitiation. During early transcript synthesis, the early elongation complex (EEC) retains a measurable tendency to undergo backtracking until Pol II pauses about 30-base-pair downstream of the transcription start site (i.e., promoter-proximal pausing) [41] and then is released for productive elongation, restricting the frequency of transcriptional initiation [42].

Fig. 1. Model of the transcriptional cycle and transcriptional regulation. (A) Schematic of biochemical reactions underlying the gene states (top) and state transitions (middle; large circles denote gene states, while small circles denote omitted states). According to whether the core promoter is open or activated, all states are divided into 5 subsets (C, F, A, I, and E; bottom). (B) Flow-process diagram of the event-based model simplified from (A). There are 6 non-Markovian events (rectangles) and 2 checkpoints (diamonds). “Y” and “N” separately represent whether the criterion is satisfied or not; ri and rj are random numbers from the uniform distribution on the unit interval. Transcription begins with Ebegin, goes through ES2, EP, ES1, and ES3, and ends with Eend. The ordering of the events is determined by p1 and p2, which are decided by the event flows (J1+, J1−, J2+, and J2−). mRNAs are synthesized in the active phase (A, red area) comprising a series of EP, while no mRNA is produced in the inactive phase (I, gray area) consisting of sequences of ES1, ES2, and ES3. (C) Gene transcription is simulated based on the distribution functions fX (X = P, S1, S2, S3), p1 and p2, which depend on the regulator occupancy rate (O). “B” denotes regulators bound and “U” regulator unbound. The left panel schematically shows fX for different O values. When O rises from 0 to 1, (fXU, p1U, p2U) turns into (fXB, p1B, p2B). (D) Eight basic regulatory modes in which regulator binding can promote transcription by decreasing the mean event duration τX and p2 or by increasing p1. Only 1 of the 6 bursting parameters is adjustable in each mode. The first “M” in the leftmost column stands for modulation, while “RF” stands for regulatory factor. The black arrows signify that RF binding accelerates the events or increases the event flows (J1+ or J2−), whereas the black lines with a flat head signify that RF binding decreases the event flows (J1− or J2+).

Each node in the network represents a gene state, and the transition between adjacent states is a Markov process (Fig. 1A). The same reactions may occur at diverse sites around the promoter, and some reactions are not exclusive with each other, leading to various transition paths between 2 states. All states are divided into 5 sets (Table 1), separately responsible for maintaining chromatin inaccessibility (C), chromatin remodeling and histone modification when the core promoter is a nucleosome-free region (F), PIC assembly (A), initiation of mRNA synthesis by Pol II and promoter escape (I), and promoter-proximal elongation and pausing of Pol II (E). When a Pol II enters productive elongation (P), the gene state returns to set I, and another Pol II is recruited for transcription reinitiation. Each state can transition to a state in the same or adjacent set.

Table 1. Gene state sets.

State set	Description of state sets	
C	States in which the access of GTFs to the core promoter is inhibited by occluded canonical nucleosomes or other macromolecules.	
F	States in which the core promoter is a nucleosome free region and ready for GTF recruitment.	
A	States in which the core promoter is occupied by GTFs but the PIC has not yet formed completely.	
I	States in which the scaffold complex is at the promoter without the EEC formed.	
E	States in which Pol II undergoes promoter-proximal elongation until it enters productive elongation.	

This network model is too complicated to be analyzed easily. To simplify it, we exploited queueing theory [28,43] to turn the network model into an event-based model. First, 9 types of functional transitions are specified: CF, FC, FA, AF, AI, IA, IE, EI, and EIP, each comprising a series of reactions with the same function (Table 2). XX′ (X, X′ ∈ C, F, A, I, or E) represents the process where the initial state stems from set X and all state transitions occur in X until the last state belongs to X′. EIP involves Pol II released into productive elongation, such that a nascent full-length transcript can be produced, whereas EI does not. The time scales of (IE, EI, EIP), (FA, AI, IA, AF), and (CF, FC) are in the order of seconds to minutes, minutes, and minutes to hours, respectively.

Table 2. Reactions in functional transitions.

Transition	Reactions during state transitions	
CF	Reactions that make the core promoter nucleosome-free to promote the recruitment of GTFs, e.g., removal of nucleosomes at the promoter, altering the accessibility of DNA on the surface of nucleosomes by remodelers, replacement with certain histone variants, or destabilization of internucleosome contacts by acetylation of lysine residues. Reactions that change the microenvironment when the promoter is occluded by canonical nucleosomes or other macromolecules different from GTFs, e.g., formation or disruption of topologically associating domain boundaries and changes in distribution modes and interactions of chromatin marks [36,75].	
FC	Reactions that change the microenvironment when the promoter is a nuclosome-free region [75].	The core promoter is wrapped around histone octamers to form canonical nucleosomes or occupied by other macromolecules [36].	
FA	Binding of GTFs, especially TFIID [75–77].	
AF	Reactions that change the microenvironment when the promoter is occupied by GTFs; enhancer-promoter proximity (hubs or looping) [75].	Dissociation of GTFs bound to the core promoter [75–77].	
AI	Formation of the PIC scaffold [78–80].	
IA	Reactions that change the microenvironment when the PIC scaffold is formed without the EEC on the DNA template; enhancer-promoter proximity [75,81].	Dissociation of the PIC scaffold components [78–80].	
IE	Stabilizing the initially transcribing complex to form the EEC via promoter escape, e.g., the B-finger of TFIIB inserts into the polymerase active site to complete escape commitment [81].	
EI	Reactions that change the microenvironment when the EEC begins elongation; enhancer-promoter proximity [75,81].	Transcript slippage and backtracking [81].	
EIP	Transcript synthesis; recovery of the elongation competency of arrested Pol II via TFIIS; release of paused Pol II to enter productive elongation [81].	

Second, we defined 6 events, a detailed description of which is presented in Fig. S1. EP refers to a series of IE and EI transitions plus the ending transition EIP, resulting in synthesis of mRNA transcripts. ES1 comprises a series of IE and EI plus IP, responsible for disassembly of the scaffold complex. ES2 consists of a series of AF and FA plus AI, enabling assembly of the PIC. ES3 comprises a series of AF, FC, CF, and FA, maintaining the inaccessibility of GTFs to the core promoter and delaying the entry into ES2. As the beginning event of the entire transcriptional process, Ebegin comprises a series of CF and FC plus FA, contributing to chromatin opening; as the ending event, Eend is composed of a series of AF and FA plus FC, resulting in termination of transcription. Notably, each event is not fixed in composition but amenable to modulation; Ebegin and Eend are governed by the initial and final states in set C, respectively. A single transcriptional burst undergoes 2 stages: The inactive phase refers to no mRNA production, including ES1, ES2, and ES3, while the active phase corresponds to mRNA synthesis, comprising a sequence of EP.

Last, we developed an event-based model for transcription (Fig. 1B). Ebegin first appears, and then ES2 follows, or ES3 occurs once or repetitively until ES2 arises. Subsequently, different paths lead to the (repetitive) occurrence of EP, with a burst of mRNAs generated. After bursting, ES1 emerges, then ES2 appears with or without ES3 preceding, and a new burst may ensue. Notably, there may exist a refractory period in ES1 where immediate reactivation of transcription is prohibited. The above process repeats until ES1 and Eend appear successively, with the gene returning to the silent state. Thus, the transcription is characterized by a sequence of ordered events. The process of mRNA splicing, nuclear export, and degradation of mature mRNA is simplified as an effective Poisson process with rate constant 1/τm. τm is set to 5 min unless specified otherwise. The biochemical master equation describing mRNA production is presented in Text S1.

Given the stochasticity in transcription, we focused on the distribution function fi of the duration of event i (i = P, S1, S2, S3) in steady state with the mean τi (Text S2 and Figs. S2 and S3). The probability of repeated occurrence of EP (ES3) is p1 (p2) (Fig. S4). In general, fP, fS1, fS2, and fS3 are primarily controlled by Pol II-dependent reactions, disassembly of the scaffold complex, GTF-dependent reactions, and nucleosome-associated reactions near the core promoter, respectively. fi is an exponential or a unimodal function in most cases [11,44]. τi is mainly determined by rate-limiting steps, whose durations in a burst are the main constituent of the burst period. p1 is elevated via accelerating such reactions as recruitment, phosphorylation, and release of Pol II, which promote the transition of A → I → E → P, while p2 is boosted by accelerating reactions that promote I → A → F → C.

Once fi, p1 and p2 were given or derived from experimental data, numerical simulation was performed to depict mRNA production using the Gillespie algorithm (Text S3). The default setting of fi, p1, and p2 is denoted as (fiU, p1U, p2U), corresponding to the situation without regulators binding to cognate sites, and is determined by the gene specificity [45] and local environment; it changes into (fiB, p1B, p2B) with bound regulators (Fig. 1C). At various regulator concentrations, its rate of occupancy at the regulatory site varies between 0 and 1, and the resulting (fi, p1, p2) lies between (fiU, p1U, p2U) and (fiB, p1B, p2B). That is, the regulation of transcription is mapped to changes in fi, p1, and p2. We found that the main conclusions do not rely heavily on the concrete form of fi provided that it is exponential or unimodal (Figs. S5 and S6). All the results presented here are based on the exponential distributions. Since an exponential distribution is determined by its mean value, the roles for regulators in gene transcription are classified on the basis of their influence on τi, p1, and p2. There exist 255 (28 − 1) regulatory modes, and Fig. 1D illuminates 8 basic modes, under each of which regulators affect τP, τS1, τS2, τS3, p1 (through J1+ or J1−), or p2 (through J2+ or J2−) alone.

Compared with previous models where independent quantities characterized experimental observables, the current model focuses on a sequence of ordered events and is more concerned with the impact of biochemical reactions on the event durations and transitions. Thus, the quantities such as the transcription rate constant (1/τP), the duration of the active phase [τA = τP/(1 − p1)] and of the inactive phase (τI = τS/p1) are all correlated (see below), which will have important implications as shown later.

Regulatory modes promoting transcriptional sensitivity and specificity

Gene transcription is mediated by a multitude of regulatory factors, such as TFs and cofactors. Without loss of generality, they are assumed to function by binding their cognate sites to affect biochemical reactions. For the general multibinding site scenario, the occupancy rate of regulators can be expressed as RnHRnH+KdnH, where [R] denotes the regulator concentration, Kd is a constant, and Hill coefficient nH could be a noninteger. As shown in Texts S4 and S5, if [R′] and K′are separately substituted for RnH and KdnH, either parts of original analytical expressions are identical to those in the case of nH = 1, or relationships between specific quantities can be preserved. Consequently, the occupancy rate can simply be written asRKd+R, corresponding to the simplest single-binding site case, for the purpose of simplifying the exploration of gene regulation.

As mentioned above, bursting parameters can be modulated to various extents under 255 regulatory modes. In mode X, the mean transcription rate 𝜐 in steady state is approximated asυ≈υ0βXn+RnΩXn+Rn,(1)

where ΩX refers to [R] at which υ is halfway between its maximum and minimum, βX is associated with the basal transcription without bound regulators, and n is a fitting number (see Text S6 for details). Kd/ΩX reflects the sensitivity of transcription to changes in [R] (Fig. 2A); high sensitivity allows for efficient regulation of transcription [46]. Fυ = υ([R] = ∞)/υ([R] = 0) = (Ωx/βx)n measures how adjustable the mean transcription rate is.

Fig. 2. Sensitivity to changes in regulator concentration ([R]) and resistance against crosstalk. (A) Schematic of the mean regulator occupancy rate (top) and steady-state transcription rate (bottom) versus the normalized concentration of regulators. Kd denotes the equilibrium dissociation constant for regulator binding, and Kd/ΩX reflects the transcriptional sensitivity. (B) Schematic of the mean transcription rate (top) and specificity S (bottom) versus [R]/Kd. The equilibrium dissociation constant is Kd and 100Kd, respectively, for cognate and nontarget binding. (C) Regulation of EIP, IE, AI, FA, or CF promotes the gene’s sensitivity to the regulator when its binding facilitates state transitions. (D) Sm versus ΩX/Kd under different regulatory modes. The default parameters are τPU = 1 min, τS1U = 5 min, τS2U =10 min, τS3U = 20 min, p2U = 0.5, and p1U = 0.5, leading to τAU < τIU (top); τPU = 10 min, τS1U = 2 min, τS2U =5 min, τS3U = 10 min, p2U = 0.5, and p1U = 0.9, leading to τAU > τIU (bottom). τPB, τS1B, τS2B, τS3B, p1B, and p2B are decided by ΩX/Kd and regulatory modes (see Table 3).

Regulators actually bind to both cognate and nontarget sites, and the latter may also induce mRNA production. To quantify the difference between the 2 scenarios, we introduced the specificity S defined as the ratio of the average expression resulting from specific binding of regulators to that resulting from nonspecific binding. S is a unimodal function of [R] for nonzero β (Fig. 2B); its maximum Sm approximately equals Fυ, given that cognate binding has a much higher affinity than nontarget binding (see Text S6 for details). Collectively, Kd/ΩX and Sm provide crucial information about transcriptional regulation.

Among 8 basic regulatory modes, ΩX < Kd is realizable for X = MEE/EA/FI/CC/FA/IE, whereas ΩX > Kd for X = MIA/FC (Fig. 1D and Table 3). In the former, the transcription is activated by reducing τP, τS1, τS2, τS3, p2 (via J2−) and increasing p1 (via J1+), which are achieved mainly by accelerating EIP and IE, IA and EI, AI and FA, FA and CF, FA and AI, and IE and EIP, respectively. In the latter, transcription is activated by increasing p1 (via J1−) and reducing p2 (via J2+), which can be acquired by slowing down IA and EI and AF and FC, respectively. Therefore, with ΩX < Kd, the transcriptional rate reaches its maximum before the regulator occupancy saturates, while the regulator binding promotes transcription by speeding up some state transitions, which involve the recruitment of TFs, enzymes, and Pol II and are part of EIP, IE, AI, FA, or CF. By contrast, with ΩX > Kd, the regulator binding contributes to transcription by slowing down some state transitions, which are engaged in EI, IA, AF, or FC. Notably, most reactions involved in EI, IA, AF, and FC are irrespective of regulator binding, like disintegration of the enhancer-promoter looping and scaffold complex. Together, the regulation of EIP, IE, AI, FA, or CF may be engaged more frequently to gain high transcriptional sensitivity (Fig. 2C).

Table 3. Gene regulatory function.

Mode	Changed parameter	Parameter change	υ max	β	Ω	
MEE	τ P	τPB = εEEτPU εEE < 1	11−p11εEEτAU+τI	ε EE K d	τAU+τIτAU+τIεEEKd<Kd	
MEA	τ S1	τS1B = εEAτS1U εEA < 1	11−p11τA+εEA1p1τS1U+τS2+τI2	ε EA K d	εEAτI1U+εEAτA+τI2εEAτI1U+τA+τI2Kd<Kd	
MFI	τ S2	τS2B = εFIτS2U εFI < 1	11−p11τA+εFI1p1τS1+τS2U+τI2	ε FI K d	εFIτI1U+εFIτA+τI2εFIτI1U+τA+τI2Kd<Kd	
MCC	τ S3	τS3B = εCCτS3U εCC < 1	11−p11τA+τI1+εCCτI2U	ε CC K d	εCCτA+τI1+εCCτI2UτA+τI1+εCCτI2UKd<Kd	
MIE	p1 via J1+	p1B = εIEp1U εIE > 1	11−p1U1τAU+1−εIEp1UεIE1−p1UτIU	1−εIEp1UεIE1−p1UKd	τAU+τIUεIE1−p1U1−εIEp1UτAU+τIUKd<Kd	
MIA	p1 via J1−	p1B = εIAp1U εIA > 1	11−p1U1τAU+1−εIAp1UεIA1−p1UτIU	K d	τAU+τIUτAU+1−εIAp1UεIA1−p1UτIUKd>Kd	
MFC	p2 via J2+	p2B = εFCp2U εFC <1	11−p11τA+τI1+εFC1−p2U1−εFCp2UτI2U	K d	τA+τI1+τI2UτA+τI1+εFC1−p2U1−εFCp2UτI2UKd>Kd	
MFA	p2 via J2−	p2B = εFAp2U εFA <1	11−p11τA+τI1+εFA1−p2U1−εFAp2UτI2U	εFA1−p2U1−εFAp2UKd	τA+τI1+τI2UτA+τI11−εFAp2UεFA1−p2U+τI2UKd<Kd	
Note: υ = υmax(β + [R])/(Ω + [R]), τA = τP/(1 − p1),τI = τI1 + τI2, τI1 = (τS1 + τS2)/p1, τI2 = p2τS3/[p1(1 − p2)]. “B” and “U” denote the cases with and without bound regulators, respectively. No “U” emerges if regulator binding has no influence. εX (X = EE, EA, FI, CC, IE, IA, FC, FA) is the ratio of the corresponding parameter with bound regulators to that without bound regulators.

The type of regulatory mode determines the trend in Sm changing with ΩX. For ΩX < Kd, Sm drops with increasing ΩX/Kd (Fig. 2D); Sm varies sharply around ΩX/Kd = 1 under MIE (MEE) for τAU < τIU (τAU > τIU). The maximum of Sm is achieved via MIE for τAU < τIU and via MEE otherwise, suggesting that regulating EIP or IE is a desirable mode. For ΩX > Kd, Sm rises linearly with increasing ΩX/Kd (on a log–log scale); obtaining a large Sm through MIA/FC requires that the binding sites should always be occupied, which is biologically implausible.

According to the event model, we directly infer that the regulation of EIP, IE, AI, FC, or CF can enhance the transcriptional sensitivity while promoting mRNA synthesis. The mode modulating EIP and/or IE has a marked role in enlarging the transcriptional adjustability and reducing the crosstalk effects resulting from noncognate regulator binding.

Major regulatory modes underlying the observed transcriptional noise

Transcription is a complex stochastic process; for simplicity, we only explored intrinsic noise in transcription, neglecting variations in parameters between cells or across time (i.e., extrinsic noise). The relationship between the mean (<m>) and variance (σm2) of mRNA copy number or Fano factor (F = σm2/<m>) is a key representation of transcriptional regulation. <m>, σm2, and F in steady state take the following forms (see Text S7 for details):m=τmτP+1−p1p1τS,(2)

σm2≈m+m2τIτm+τAτmτSτm−τAτIτm2τIτm+τAτm+τAτIτm2,(3)

F≈1+bc+ddc+db−1/b−cdc+d+cd,(4)

where τm is the lifetime of mRNA, τA (τI) is the mean duration of the active (inactive) phase, and τS is the mean interval between ES1 and ES2: τA = τP/(1 − p1), τI = τS/p1, and τS = τS1 + τS2 + p2τS3/(1 − p2). The Fano factor is governed by the mean burst size b [the number of mRNAs per burst; b = 1/(1 − p1)], c = τA /τm, and d = τI /τm. A small burst size, a plateau-like bursting profile, or a short inactive phase can contribute to reducing F.

Figure 3A displays how F varies with <m> over a wide range of parameter values. For each curve, only one of the 6 parameters is altered, while the others are fixed at default values; we also performed simulations by directly setting τS rather than τS1, τS2, τS3, and p2 and assumed τS to obey an exponential distribution. When y is varied, the corresponding F is denoted as Fy. As Fτs differs slightly from Fx (x = τS1, τS2, τS3, or p2) (Fig. 3A, inset), only the curve for Fτs is presented. With τP regulated alone, the mean height of burst profiles {h ≈ τm(1 − e−c)(1 − e−d)/[τP(1 − e−c−d)] for τP < τm and h ≈ 1 otherwise}, rather than burst size, varies prominently. At small τP, F rises sharply with increasing <m> due to little degradation of mRNA within a short active period. At large τP, the transcription rate is small such that a nascent mRNA has a high probability of being degraded before new mRNA production, and the active phase is much longer than the inactive phase, i.e., the transcription is nearly a Poisson process (F = 1). With τS modulated alone, b and τA show little change, whereas τI varies markedly, such that F drops toward one with increasing <m>. F varies slightly at small <m> but drops markedly at large <m> due to a continuous reduction indFdm≈1−p12τm3τP31+τmτPm−1−p1τm2τP2−τmτP2−1.(5)

Fig. 3. Relationships between the mean mRNA number <m> and Fano factor F under various regulatory modes. Time is in units of minutes. (A and B) Solid lines are from analytic expressions (Eqs. G21 to G23 in Text S7), while symbols label numerical results. (A) F versus <m> for different parameter values. The dashed line represents F = 1 (Poisson noise). The default parameters without bound regulators are p1U = 0.95 and τPU = 0.165, and τSU is assumed to be exponentially distributed with an average of 100 min (cyan dot). On each curve, only one parameter (p1, τP, or τS) varies, while the others are fixed at default values. The top panel depicts the sketch of transcriptional bursts (red bars represent Pol II entering productive elongation). The right bottom panel shows the corresponding burst height h versus m, while the right top one shows Fτs versus FX (X = τS1, τS2, τS3, p2) with τS1U = 10, τS2U = 30, τS3U = 60, and p2U = 0.5. (B) F versus <m> in different regulatory modes. Regulators affect the transcription by changing the following parameters, while the others remain fixed (p1U = 0.95, τPU = 0.165, τSU = 100): τP (red circle; τPB = 0.01), τS (black square; τSB = 1), p1 (blue diamond; p1B = 0.999), τP and p1 (pink rightward triangle; τPB = 0.01 and p1B = 0.999), p1 and τS (green leftward triangle; p1B = 0.999 and τSB = 1), τP and τS (yellow upward triangle; τPB = 0.01and τSB = 1); and τP, p1, and τS (cyan downward triangle; τPB = 0.01, p1B = 0.999, and τSB = 1). Each data point represents an average over 10 thousand trajectories. Here, the regulator binding promotes transcription. Two right panels show the cases where both p1 and τS are regulated with τP fixed (p1U = 0.7, τSU = 500, p1B = 0.999, τSB = 1, τP = 0.165; top), or p1, τS, and τP are regulated (τPB = 0.01, p1B = 0.999, τSB = 1, τPU = 0.2, p1U = 0.5, τSU = 100; bottom). (C and D) Experimental data on gap genes in early Drosophila embryos [47] (C) and on CDKN1A (measured by p21-MS2 signals) in MCF7 cells [49] (D). The red curve is a fit to all data with the event model. a.u., arbitrary units.

With p1 altered alone, b rises with increasing <m>, and h rises toward saturation; thus, F first rises because of an increase in b and then drops because of saturation in h and plateau-like bursting. Strikingly, all numerical results agree perfectly with the analytical expressions.

The dependence of F on <m> is more complicated when 2 or more parameters vary concurrently. Notably, the curves for FτP and FτS constitute the boundaries, and the others lie in between and exhibit unique features (Fig. 3B and Fig. S7). Providing that the regulator binding promotes transcription, F rises markedly at small <m> together with a notable change in h or varies slightly when h changes little. If F falls at intermediate or large <m>, then at least τS or p1 is regulated, and <m> tends toward τm/τP when F approaches 1; specifically, F drops rapidly when τP is sufficiently large or changes little. Together, the regulation of τP or p1 is indispensable for a rise in F, and its rapid rise is accompanied by sharp bursting; the regulation of τS or p1 is essential to a drop in F, and its fast fall is associated with plateau-like bursting. Similarly, the numerical and analytic results match exactly.

Indeed, the fitting curves to experimental data resemble the typical curves above. Activated by the Bicoid (Bcd) protein, major gap genes in early Drosophila embryo are expressed at distinct levels along the anterior–posterior axis of its body due to the Bcd concentration gradient [47], sharing a similar trend in F versus <m> (i.e., F first rises and then drops toward 1 with increasing <m>) (Fig. 3C; see also the right top panel in Fig. 3B). All the data points can be fitted by a single curve with the event model (see parameter inference in Text S4). For each gene, F is greater than 3.8 at small <m> and falls rapidly after the maximum, indicating that τP is nearly consistent across body locations. With τP fixed, there should exist the upper bounds, <m> < τm/τP and F < 1 + τm/τP (see Text S7 for details). As these genes display similar limits in <m> and F and have similar τm values [48], τP (and the frequency of transcription initiation) should also be comparable among them, which implies that variations in Bcd concentration have a minor impact on the major rate-limiting steps in Pol II recruitment and promoter-proximal pausing. The increase in Bcd concentration may primarily facilitate the activation of Pol II and its entry into productive elongation.

Another example involves the expression of CDKN1A activated by the tumor suppressor p53 in MCF7 cell lines [49]. F varies slightly at small <m>, rises to a peak, and then drops gradually (Fig. 3D); to fit the data with the event model, p1, τS, and τP vary concurrently with small p1U and τP (see also the right bottom panel in Fig. 3B). At small p1U, the gene state tends to leave set E via EI rather than EIP or leave set I via IA rather than IE, associated with slow reactions in recruiting, stabilizing, or releasing Pol II. At small τP, fast reactions are required before Pol II enters productive elongation. Thus, p53 may prefer to regulate CDKN1A expression via speeding up EIP through cellular signaling.

Further analysis of other experimental data reveals the diversity in transcriptional regulation. When the ctgf expression is mediated by transforming growth factor-β (TGF-β), the transcription rate rises with increasing the TGF-β concentration, and τA is nearly fixed [50], which requires that τP is inversely proportional to b. This further implies that the likelihood of Pol II entering productive elongation is prominently higher than its probability of escaping from the core promoter, i.e., the modulation lies in EIP rather than IE. When glucagon-like peptide 1 (GLP-1)/Notch induces sygl-1 expression in Caenorhabditis elegans, the burst size is nearly irrelevant to the concentration of GLP-1/Notch ligand [51], i.e., τP varies markedly but b (p1) is fixed, implying that the early elongation of Pol II until promoter-proximal pausing can be modulated. Obviously, different regulations of τP and p1 are gene-specific, reflecting which steps are rate-limiting and which reactions are being adjusted in EP.

Together, the event model enables the deduction of transcriptional modulatory mechanisms. The results above suggest that the changing trend of burst profiles is influenced by the adjustable range of Pol II-dependent reaction rates (e.g., recruitment and promoter-proximal elongation of Pol II). In the 3 limiting cases where the transcription occurs with fixed τA, b, or τP, the burst height, burst duration, and area under the curve are modulated, respectively.

Global constraints on transcriptional bursting

The results above unravel an important regulatory mode involving the concurrent modulation of p1 and τS. Under this mode, the mean mRNA number (<m> = υτm) and burst size b rise toward saturation with increasing [R], while the burst frequency f first rises to a peak and then drops toward saturation, corresponding to τA < τI and τA > τI, respectively (Fig. 4A). Ωm, Ωb, and Ωf refer to the operating points where m, b, and f reach half of their maximum values, respectively. Because Ωm < Kd, Ωb ≈ Kd when EIP, IE, AI, FA, or CF is modulated, and <m> = bfτm, Ωf < Ωm < Ωb generally holds true, and thus f, <m>, and b reach their respective maximum values successively (see Text S4 and Fig. S8 for details). Furthermore, f is sensitive to changes in <m> at low expression, while b is insensitive, and vice versa at high expression (Fig. 4B, top).

Fig. 4. Modulation of transcriptional bursting. Time is in units of minutes. (A) Mean mRNA count <m>, burst size b, and burst frequency f, normalized by their respective maximum values, versus the normalized regulator concentration [R]. The regulator binding affects p1 and τS; τP = 0.2, p1B = 0.995, p1U = 0.5, τSU = 250, τSB = 5 and n = 1. (B) Burst size and frequency versus the mean mRNA number (top); Sb and Sf versus <m> (bottom). The dashed line marks the case of τI = τA. The black arrows point to <m> with Sf = Sb. The parameters are the same as in (A). The curves in (A) and (B) are plotted on the basis of analytic expressions. (C) A fit to the experimental data (symbols) in [47] using the event model with τP = 0.139, p1U = 0.401, p1B = 0.995, τSU = 498, and τSB = 0.694 (solid line). (D) Experimental data on the mean expression level and Fano factor in [52]. Each point is from an individual gene with the shadow showing the range. The red fitting curve is obtained by Gaussian regression. The black curve shows the lower bound on Fano factor at large <m>. The inset shows ηm2=σm2/<m>2 versus <m>, with the red dashed line representing the mean. (E) The red lines show Fano factor versus <m>, and the black line shows the lower limit of Fano factor at large <m> with τS ≥ 1. The default parameters are the same as in (A). (F) Fano factor versus <m> with the event model or telegraph model (inset). Ten thousand sets of bursting parameters are used in each panel, where the duration of a burst varies between 1 and 1,000 min, the transcription rate is from 0.1 to 100/min and τA < τI. The refractory period is set to 1 (blue), 5 (red), or 10 min (yellow).

We also used Sb=dlnbdlnm and Sf=dlnfdlnm (Sb + Sf = 1; Sb and Sf > 0 for τA < τI), i.e., their sensitivity to changes in mRNA number, to dissect the contribution of size and frequency modulation to mRNA production. Sf and Sb drop and rise respectively with increasing <m> (Fig. 4B, bottom), indicating that the modulation of burst frequency and size predominates in distinct regimes and a combination of both mechanisms emerges over intermediate ranges. A similar relationship between f(b) and <m> was observed across gap genes in early Drosophila embryos [47] (Fig. 4C), where their expression levels are graded along the anterior–posterior axis, responsible for body segmentation.

With high-throughput experimental data available, we calculated the mean mRNA numbers and Fano factors for 9196 genes in fibroblasts [52]; to exclude the influence of extrinsic noise, we preprocessed the data using the method in [53]. Each mRNA count is normalized by dividing it by the total count in each sample and scaled by the median count across all samples. Each data point in Fig. 4D denotes <m> and F for one gene, and the red curve is a fitted trendline for all data. There exist lower bounds on Fano factor and burst size at each expression level; meanwhile, ηm2=σm2/<m>2 fluctuates around the mean for large <m> (Fig. 4D, inset), while the total noise intensity in the original data from [52] remains nearly constant at high expression levels [54]. The lower bound on F is nearly irrelevant to <m> at low expression but rises prominently with increasing <m> at high expression. That is, higher expression is accompanied by increased burstiness. This phenomenon has been widely observed in mammalian cells [44,55–58], fly embryos [59–61], and even Escherichia coli [62,63], suggesting that transcriptional regulation may obey a universal principle across the genome, i.e., high-level expression is associated with the bursting with large size.

This global constraint results from the correlations among bursting parameters. Given that m=τmτP+1−p1p1τS, τp is expressed as a function of p1, τS, and <m>, and, thus, F = 1 + 1−p1mmτS−p1τmτSp12τm2+1−p1mτS2. For fixed <m>, the minimum of F, Fmin, is determined by the ranges of τS and p1; if <m> > τm/(2τS), Fmin equals 1 + 1−p10mmτS0−p10τmτS0p102τm2+1−p10mτS02, where τS0 is the minimum of τS and p10 is the maximum of p1. Fmin is greater than one and rises with increasing <m> at high expression, as illustrated by the black curve in Fig. 4E, which is similar in trend to the black one in Fig. 4D (to justify that this global restriction is intrinsic to the transcriptional machinery, few constraints are imposed on parameter values; otherwise, if τI > τA were guaranteed, the 2 curves would match). A cluster of red curves is plotted in Fig. 4E, each beginning with the same parameters and ending with a point on the black line. Each curve represents the F-<m> curve for some gene, similar to that shown in Fig. 3B.

To recapitulate the lower bound on F at high expression levels in Fig. 4D, we systematically altered the values of p1, τS, and τP and calculated the mean mRNA number and Fano factor in the case of τA < τI. The scatter plot of F versus m has similarities to that in Fig. 4D (Fig. 4F). The global constraint becomes more striking when there exists a longer refractory period in the inactive phase or a higher transcription rate, such as in the case of large τS and small τP, leading to sharp bursting peaks and enhancing the discernibility of burst profiles. By contrast, high-level expression is not necessarily associated with large Fano factor in the telegraph model and other phenomenological models (inset of Fig. 4F and Fig. S9), which arises from the independence of bursting parameters. Collectively, the global constraint reflects an inherent feature of transcriptional bursting.

Transient dynamics reflect the underlying regulatory mode

After showing the typical features of transcriptional bursting in steady state, we turned to explore the transcriptional response to a step rise in regulator concentration. A new steady state will be reached in diverse manners because of different regulatory modes and composition of the inactive phase. With the initial state completely silent (i.e., p1 → 0, p2 → 1, or τi → ∞ and no mRNA production), the mean probability PA(t) of gene activation and mean mRNA number are determined as follows (see Text S8 for details):PAt≈fIinitt∗δt∑i=0∞∗fAt∗fIti∗e−tτA<m>t≈1τPe−tτm∗fI_initt∗δt∑i=0∞∗fAt∗fIti∗e−tτA(6)

where fA(tA), fI(tI), and fI_init(tI_init) separately denote the distribution functions of the duration of the active (inactive) phase in the new steady state and the time taken to enter the active phase (Fig. 5A, top), with fAtA=1τAe−tAτA and fI(tI) ≈ Γ(tI, αI, θI) (Gamma distribution), while * represents convolution. Depending on the initial condition, fI_init can assume different forms (Fig. S10). ffirst = fA * fI_init represents the distribution function of the first burst duration tfirst, which is the time from stimulus onset to completion of mRNA synthesis in the first active phase, and can be expressed as Γ(tfirst, αf, θf). fT = fA * fI denotes the distribution function of burst duration in steady state (tss). The mean and coefficient of variability (CV) of tfirst are τfirst and CVfirst, respectively. The peak and steady-state values of PA are separately PAP and PAss = τA/T (T is the mean burst period) (Fig. 5A, middle); <m>p and <m>SS = bτm/T denote the mean peak and steady-state values of mRNA count, respectively (Fig. 5A, bottom). CVZ (Z = P, S1, S2, S3, S) refers to the CV for τZ in steady state.

Fig. 5. Response patterns upon stimulation. Time is in units of minutes. (A) Schematic of time courses of [R], the mRNA number (black; red bars represent Pol II entering productive elongation), the probability of gene activation, and the mean mRNA count (from top to bottom). The panel in the second row also marks the durations of the active phase (tA), inactive phase (tI), first inactive phase (tI_init), first burst (tfirst), and a burst in steady state (tss). (B) Time courses of PA with τA = 10, αI = 2, and θI = 10. The α and θ of ffirst are 3 and 5 (τfirst = 15, CVfirst2 = 1/3), 3 and 10 (τfirst = 30, CVfirst2 = 1/3), 3 and 15 (τfirst = 45, CVfirst2 = 1/3), 5 and 6 (τfirst = 30, CVfirst2 = 1/5), and 2 and 15 (τfirst = 30, CVfirst2 = 1/2), respectively. The data points are from simulation, while the solid curves are from analytical expression (Eq. H2 in Text S8). (C) Time courses of the mean mRNA number. fA and fI are the same as in (B); the α and θ of ffirst are 3 and 5 (τfirst = 15, CVfirst2 = 1/3) for τm = 1, 10, or 20, and PA overshoots. In the inset, the α and θ of ffirst are 2 and 15 (τfirst = 30, CVfirst2 = 1/2) for τm = 1, 10, or 20, and PA does not overshoot. <m>(t) is normalized by <m>ss. The data points are from simulation, while the solid curves are from Eq. H5. (D) Color-coded HP controlled by τfirst/T and CVfirst. The black points signify that PA does not overshoot with hPT < 1, while the red curve corresponds to hPT = 1. Small values of τfirst/T and CVfirst lead to large HP. (E) Color-coded Hm controlled by τfirst/T, CVfirst, and T/τm. The black points denote that m does not overshoot with hmT < 1, while the red surface represents hmT = 1. (F) Color-coded tre controlled by τfirst/T, CVfirst, and T/τm. In (D) to (F), the initial state is completely silent, and τA ∈ (0.1,10), τI ∈ (0.1,100), CVP = 1, CVi ∈ (0.1,1) (i = S1, S2, S3), and τI_init = xτI with x ∈ (0.1, 10).

Figure 5B shows the typical dynamics of PA(t): PA either overshoots PAss or monotonically rises toward PAss. Without PA overshooting, <m>p cannot exceed <m>ss, and its dynamics also markedly depend on τm (Fig. 5C). To determine the condition for overshoot, we introduced HP = (PAP − PAss)/PAss and Hm = (<m>p − <m>ss)/<m>ss. PA (<m>) must overshoot for hPT > 1 (hmT > 1), where hx=αxδx+τfirstαx−1αx−1Γαxe−αx−1 (x = P or m) with αP=1CVfirst2, αm=τm+τfirst2τm2+τfirst2CVfirst2, δP = 0, and δm = τm. HP is governed by τfirst/T and CVfirst; the area below the red curve of hp = 1/T corresponds to Hp > 0 (Fig. 5D). Hm also relies on τm/T, and the phase plane of hm = 1/T separates Hm = 0 from Hm > 0 (Fig. 5E). Strikingly, overshoot appears at small values of τfirst/T, CVfirst, and τm/T.

On the other hand, the time taken to reach half of PAss (<m>ss) is denoted as tGA(tre); tre is affected by mRNA production and degradation. If the rising phase of <m>(t) is governed by mRNA production over the time scale of τm due to strong correlation of bursts among individual simulations, tre can be less than τmln2; otherwise, the dynamics of <m> are controlled by mRNA degradation, leading to tre > τmln2. Thus, decreasing τfirst/T and CVfirst contributes to reducing tre; changing T/τm has a dual effect, and exclusively increasing it leads to a first drop and then rise in tre (Fig. 5F).

In general, τfirst/T and CVfirst are modulated by the initial state distribution and constrained by the promoter architecture. When the gene state immediately enters the active phase upon stimulus onset due to regulation of Pol II-dependent reactions, τfirst/T and CVfirst are separately close to τA/T and the CV of the active phase duration (CA, which is often ~1), whereas τfirst/T and CVfirst separately approach 1 and the CV of tss (i.e., CVss) when the state just leaves the active phase. If there are multiple rate-limiting steps in burst cycle before stimulus onset [11,44], the initial state is mostly distributed around those states, usually leading to τA/T < τfirst/T < 1 and small CVfirst, whereas τfirst/T ~ 1 and CVfirst ~ 1 only if there is one rate-limiting step. The number of rate-limiting steps is governed by the promoter architecture [11,64,65] and is associated with the magnitude of p2. For large p2, the sum of all durations of ES3 within a burst cycle often accounts for the majority of the burst period, and its distribution is approximately exponential, implying that large p2 is often associated with only one rate-limiting step (i.e., nucleosome clearance) (see Text S2). By contrast, formation of the scaffold complex and pause of Pol II are the major rate-limiting steps for small p2. For instance, the nucleosome occupancy rate is rather low at target genes of white collar complex in Neurospora crassa, which is linked with small p2 and relatively large τS1 + τS2, leading to the observable overshoot [46]. Collectively, analyzing the transient response also provides insight into the modulatory mechanism and promoter architecture.

Notably, only a finite number of scenarios emerge from the comparison of transient transcriptional responses across different regulatory modes. When the gene is activated from the silent state to the same steady state (with identical fA and fI), for example, 4 scenarios appear upon examining the ordering of tre values across 6 basic regulatory modes (Fig. 6A). These relationships primarily depend on the composition of the inactive phase, i.e., the values of p2, τS2/τS1, and τS3/τS1 (Fig. 6B). Nevertheless, the fastest response always occurs via MEE, as paused Pol II is poised to enter productive elongation upon an excitatory signal [66]; the second fastest response takes place via MFI, as the transcription machinery is ready for assembly on the promoter (Figs. S11 and S12). Under combinatory modes, a faster response can be induced to reach the same <m>ss when τP is regulated (Fig. 6C). Therefore, accelerating the reactions involving Pol II can elicit fast transcription from a completely silent state.

Fig. 6. Transcriptional responses in different regulatory modes. Time is in units of minutes. (A) Transcriptional responses under different regulatory modes and initial conditions. There emerge 4 scenarios in terms of tre: MEE < MFI < MFA, MCC < MIE < MEA for τS1 = 15, τS2 = 2, and τS3 = 4 after the stimulation; MEE < MFI < MIE < MFA, MCC < MEA for τS1 = 10, τS2 = 20, and τS3 = 6; MEE < MFI < MEA < MIE < MFA, MCC for τS1 = 3, τS2 = 10, and τS3 = 15; and MEE < MFI < MIE < MEA < MFA, MCC for τS1 = 5, τS2 = 20, and τS3 = 8. The other parameters are fixed: τP = 1, τm = 10, p1 = 0.9, and p2 = 0.5. <m> is normalized by its steady-state value. (B) Two-parameter phase diagram of x = τS2/τS1 and y = τS3/τS1 at p2 = 0.5. The black curves constitute the boundaries. The blue region with y < 1/(x+1), red region with 1/(x+1) < y <1, purple region with 1<y<1+1+4x1−p2/p2/2, and yellow region with y>1+1+4x1−p2/p2/2 correspond to the 4 scenarios in (A). (C) <m> and tre for different regulatory modes and bursting parameters. The color represents the bursting parameters regulated. Default parameters are τPU = 1, τSU = 10, p1U = 0.5, τm = 5, CVP2 = 1, and CVS2 = 1. Under a regulatory mode, the regulated parameters are decided by the regulator concentration, default parameters, and the corresponding parameters with bound regulators (τPB = 0.1, τSB = 1, and p1B = 0.95).

Moreover, it may be beneficial for an inducible gene to respond fast to stimuli and generate a broad range of outputs. The rapid response largely results from the overshoot caused by regulating τP via MEE, while regulating τP in the case of τA > τI or regulating p1 and τS in the case of τA < τI elicits a broad range in <m>. Thus, integrating MEE with MX (X = IE, EA, FI, CC, or FA) allows for a fast response and highly tunable expression (Fig. 6C). In this sense, regulating EIP and IE is conducive to inducible genes [67].

Discussion

The event model showcases several strengths when compared with previous phenomenological models. First, biochemical reactions involved in transcription are accurately mapped to the events in our model, and parameter values can be directly inferred from experimental data. Second, the event model retains the correlations among the bursting variables [transcription rate constant, burst size, and duration of the (in)active phase]. This naturally imposes constraints on burst profiles, guaranteeing bursty transcription in highly expressed genes [52,68,69] and enabling the modulation of burst frequency and burst size to separately predominate at low and high expression levels [47,57,70]. Third, 9 kinds of functional state transitions are distinctly separated in time scales from seconds to hours, and the durations and occurrence probabilities of events are independently modulated. These enable transcriptional bursting to manifest across a broad spectrum of temporal scales, exhibiting diverse dynamic patterns. With these features, the event model could have wide applications.

In principle, we can deduce the changes in (τP, τS, p1) through the correlations among observables and further reveal the promoter architecture and molecular regulatory mechanisms. For example, EP consists of sequential stages (initiation, promoter clearance, elongation, and promoter-proximal pausing). The tendency of Pol II to continually synthesize mRNAs and the duration of promoter-proximal pause determine p1 and τP, respectively. For a group of gap genes in early Drosophila embryos [47], τP is nearly fixed, but p1 is regulated, suggesting that regulators strongly affect the stability of the transcription complex and weakly affect the motion of Pol II along DNA templates, e.g., modulating the phosphorylation of Pol II C-terminal domain. In the regulation of ctgf expression by TGF-β1 [50], the duration of the active phase varies only slightly, while the transcription rate constant rises with increasing the TGF-β1 concentration, suggesting that TGF-β1 affects the ratio of forward to backward rates of Pol II movement along DNA templates during the stalling and reverting of the EEC. In the regulation of sygl-1 expression in C. elegans, the burst size remains relatively constant, while τA experiences changes [51], suggesting that GLP-1/Notch influences the duration of a highly irreversible stage, e.g., early elongation of Pol II until promoter-proximal pausing. Together, although the major stages in the event model are similar across genes, the reversibility and rate limits of those stages govern the regulatory modes, reflected in the changes to the durations and ordering of the events.

Cells are exposed to various fluctuations, and some have to rapidly respond and adapt to new environments via transcription [67]. For inducible genes with rapid transcriptional activation [66], it is desirable to trigger highly tunable expression and allow transcription overshoot to lift the limitation of mRNA lifetime on response speed. This can be achieved when τP and/or p1 (with τA > τI) or τS and/or p1 (with τA < τI) are regulated. Meanwhile, transcriptional overshoot can be evoked when Pol II is poised for transcription, or there exist multiple rate-limiting steps, e.g., a refractory period and PIC recruitment during the inactive phase. The genes are indeed prone to overshoot and respond fast in the presence of a refractory period [46], helpful for the rapid and accurate transmission of information. On the contrary, eukaryotic genes with TATA box often burst with a large size and seldom overshoot [11], conducive to filtering environmental fluctuations, and the TATA box seems to be correlated with lack of a refractory period.

The event model engages various types of regulatory factors; nevertheless, only the concentration of one regulator was altered, while the concentrations of others were kept constant in simulations. When transcriptional modulation involves coordinating multiple regulators and the relationships between their concentrations are known, it is possible to elucidate the combinatorial control of bursting parameters. We could determine the dependence of bursting parameters on the concentration of one of the involved regulators and generate <m>-F curves under diverse conditions. Note that the conclusions drawn in Figs. 3A, 4D to F, and 5 rely solely on the values of τP, τS1, τS2, τS3, p1, and p2 and are not influenced by the single-regulator assumption. In Figs. 3C and D and 4A to C, the experimental data were collected when only one kind of regulator was monitored. If the concentrations of multiple regulators change simultaneously, the trends depicted in Fig. 3B may undergo changes, although the <m>-FτP and <m>-FτS curves still constitute the boundaries.

Transcriptional noise can be divided into intrinsic and extrinsic noise [71]. Here, the intrinsic noise is mainly manifested in the stochasticity in event orders, event durations, and mRNA degradation. Treating mRNA degradation as a one-step process is a simplification; however, it was reported that, in a broad range of parameter space, models assuming either multistep or one-step degradation of mRNA yield indistinguishable mRNA count distributions [16]. The current work did not explore extrinsic noise, which is influenced by various factors including cell volume, gene copies, DNA replication, and cell cycle progression [72,73]. For example, it has been shown that the presence of 2 uncorrelated alleles halves the noise intensity [56]. It was reported that extrinsic noise typically leads to a more dispersed mRNA count distribution, such as an increased probability of low numbers or a longer tail, and that extrinsic noise contributes substantially to the total noise at high expression levels, leading to a higher noise plateau [11]. It was also demonstrated that parameter variability among cells should be taken into account to fully interpret the experimental data [54]. Furthermore, when analyzing data obtained through single-cell measurements, technological noise is also unavoidable [53]. All these suggest that the event model has to be extended to incorporate more sources of noise, and it would be interesting to dissect the contributions of intrinsic and extrinsic noise.

The current work explored the generic case where the dwell times of regulators are much shorter than the time taken to switch between the inactive and active phases, regulator binding promotes transcription with monotonic gene-regulatory functions, and only one gene locus is involved [5,34,35,74]. More advanced models could be developed to probe the gene-specific transcriptional activity. If the concentration of regulators varies periodically (e.g., the pulsing of p53 and nuclear factor κB levels), how transcriptional bursting is modulated to reliably code signals could be the focus of further study.

In summary, the event model, built upon a small set of adjustable functional events, offers a clear and straightforward description of transcriptional progression. This approach facilitates the investigation of both shared and specific mechanisms for gene regulation. Modulation can occur at any stage of a burst cycle, but the underlying regulatory modes are limited in number because of intrinsic correlations among transcriptional events. Our results suggest that large Fano factor is required for transcriptional bursting at high expression levels and the transitions EIP and IE, involved in the recruitment, pause and release, or elongation of Pol II, are pivotal regulatory targets for enhancing the transcriptional sensitivity, anti-crosstalk, adaptability, and responsiveness. Combining experimental data with the event model allows for deducing the molecular processes that predominantly govern transcriptional activity.

Methods

The development of the event model is an important part of this work, and the model is the basis of the whole study. Thus, we provided the details on building the model in the “Establishment of the event model” section. Texts S1 to S8 further present the derivation of formula and explanations. We made a detailed reference to the Supplementary Materials in the text.

Acknowledgments

We thank Y. Wang for comments on an early version of this paper and 2 anonymous reviewers for constructive comments and suggestions. The numerical calculations in this paper have been performed on the computing facilities at the High Performance Computing Center (HPCC) of Nanjing University.

Funding: This work was supported by the National Natural Science Foundation of China (11874209 and 12090052).

Author contributions: F.L., R.W., and W.W. conceived and designed the project. R.W. performed the study. F.L., R.W., and W.W. wrote the paper. All authors analyzed the data and read/edited the paper.

Competing interests: The authors declare that they have no competing interests.

Data Availability

All data needed to evaluate the conclusions in the paper are present in the paper. Custom codes will be available upon request to F.L.

Supplementary Materials

Supplementary 1 Text S1. Establishment of the event model.

Text S2. Structure of the event model.

Text S3. Simulation of the event model and phenomenological models.

Text S4. Regulation of transcriptional bursting.

Text S5. Comparison between the single-binding and multibinding sites cases.

Text S6. Gene regulatory function.

Text S7. Mean and variance of mRNA numbers.

Text S8. Transient response to stimulation.

Fig. S1. Events in transcriptional bursting.

Fig. S2. Duration distribution of the process that consists of 2 sequential processes.

Fig. S3. Duration distribution of a mixed process.

Fig. S4. Probability of repeating events EP and ES3.

Fig. S5. Relationships between the mean (<m>) and variance (σm2) of mRNA numbers for different distributions of fP and fS.

Fig. S6. Relationships between the mean (<m>) and relative noise strength (𝜂m2 = 𝜎m2/<m>2) of mRNA numbers for different fP and fS.

Fig. S7. Relationships between the mean (<m>) and variance (σm2) of mRNA numbers under different regulatory modes.

Fig. S8. Transcriptional kinetics under different combination modes.

Fig. S9. Difference between the telegraph and event models under different regulatory modes.

Fig. S10. Response curves for different distribution functions of the time taken to activate the gene from the initial state.

Fig. S11. Regulatory range (υmax) and response speed for basic regulatory modes.

Fig. S12. Transcriptional sensitivity under different regulatory modes.
==== Refs
References

1. Neuert G, Munsky B, Tan RZ, Teytelman L, Khammash M, van Oudenaarden A. Systematic identification of signal-activated stochastic gene regulation. Science. 2013;339 (6119 ):584–587.23372015
2. Ezer D, Moignard V, Göttgens B, Adryan B. Determining physical mechanisms of gene expression regulation from single cell gene expression data. PLOS Comput Biol. 2016;12 (8 ): e1005072.27551778
3. Lis JT. A 50 year history of technologies that drove discovery in eukaryotic transcription regulation. Nat Struct Mol Biol. 2019;26 (9 ):777–782.31439942
4. Coulon A, Chow CC, Singer RH, Larson DR. Eukaryotic transcriptional dynamics: From single molecules to cell populations. Nat Rev Genet. 2013;14 (8 ):572–584.23835438
5. Lu F, Lionnet T. Transcription factor dynamics. Cold Spring Harb Perspect Biol. 2021;13 (11 ): a040949.34001530
6. Wissink EM, Vihervaara A, Tippens ND, Lis JT. Nascent RNA analyses: Tracking transcription and its regulation. Nat Rev Genet. 2019;20 (12 ):705–723.31399713
7. Chen X, Zhang D, Su N, Bao B, Xie X, Zuo F, Yang L, Wang H, Jiang L, Lin Q, et al. Visualizing RNA dynamics in live cells with bright and stable fluorescent RNAs. Nat Biotechnol. 2019;37 (11 ):1287–1293.31548726
8. Peccoud J, Ycart B. Markovian modeling of gene-product synthesis. Theor Popul Biol. 1995;48 (2 ):222–234.
9. Pedraza JM, Paulsson J. Effects of molecular memory and bursting on fluctuations in gene expression. Science. 2008;319 (5861 ):339–343.18202292
10. Zhang J, Chen L, Zhou T. Analytical distribution and tunability of noise in a model of promoter progress. Biophys J. 2012;102 (6 ):1247–1257.22455907
11. Zoller B, Nicolas D, Molina N, Naef F. Structure of silent transcription intervals and noise characteristics of mammalian genes. Mol Syst Biol. 2015;11 (7 ):823.26215071
12. Corrigan AM, Tunnacliffe E, Cannon D, Chubb JR. A continuum model of transcriptional bursting. eLife. 2016;5 : e13501.
13. Tantale K, Mueller F, Kozulic-Pirher A, Lesne A, Victor J-M, Robert M-C, Capozi S, Chouaib R, Bäcker V, Mateos-Langerak J, et al. A single-molecule view of transcription reveals convoys of RNA polymerases and multi-scale bursting. Nat Commun. 2016;7 :12248.27461529
14. Zhang J, Zhou T. Promoter-mediated transcriptional dynamics. Biophys J. 2014;106 (2 ):479–488.24461023
15. Szavits-Nossan J, Grima R. Steady-state distributions of nascent RNA for general initiation mechanisms. Phys Rev Res. 2023;5 (1): 013064.
16. Braichenko S, Holehouse J, Grima R. Distinguishing between models of mammalian gene expression: Telegraph-like models versus mechanistic models. J R Soc Interface. 2021;18 (183 ):20210510.34610262
17. Karmakar R. Control of noise in gene expression by transcriptional reinitiation. J Stat Mech. 2020;2020 : 063402.
18. Karmakar R, Das AK. Effect of transcription reinitiation in stochastic gene expression. J Stat Mech. 2021;2021 : 033502.
19. Cao Z, Filatova T, Oyarzún DA, Grima R. A stochastic model of gene expression with polymerase recruitment and pause release. Biophys J. 2020;119 (5 ):1002–1014.32814062
20. Muthukrishnan AB, Kandhavelu M, Lloyd-Price J, Kudasov F, Chowdhury S, Yli-Harja O, Ribeiro AS. Dynamics of transcription driven by the tetA promoter, one event at a time, in live Escherichia coli cells. Nucleic Acids Res. 2012;40 (17 ):8472–8483.22730294
21. Lloyd-Price J, Startceva S, Kandavalli V, Chandraseelan JG, Goncalves N, Oliveira SMD, Häkkinen A, Ribeiro AS. Dissecting the stochastic transcription initiation process in live Escherichia coli. DNA Res. 2016;23 (3 ):203–214.27026687
22. Weidemann DE, Holehouse J, Singh A, Grima R, Hauf S. The minimal intrinsic stochasticity of constitutively expressed eukaryotic genes is sub-Poissonian. Sci Adv. 2023;9 (32 ):eadh5138.37556551
23. Hansen AS, Zechner C. Promoters adopt distinct dynamic manifestations depending on transcription factor context. Mol Syst Biol. 2021;17 (2 ): e9821.33595925
24. Scholes C, DePace AH, Sánchez Á. Combinatorial gene regulation through kinetic control of the transcription cycle. Cell Syst. 2017;4 (1 ):97–108.e9.28041762
25. Wang Y, Liu F, Wang W. Dynamic mechanism for the transcription apparatus orchestrating reliable responses to activators. Sci Rep. 2012;2 :422.22639730
26. Wang Y, Ni T, Wang W, Liu F. Gene transcription in bursting: A unified model for realizing accuracy and stochasticity. Biol Rev. 2019;94 (1 ):248–258.30024089
27. Wang Y, Liu F, Wang W. Kinetics of transcription initiation directed by multiple cis-regulatory elements on the glnAp2 promoter. Nucleic Acids Res. 2016;44 (22 ):10530–10538.27899598
28. Jia T, Kulkarni RV. Intrinsic noise in stochastic models of gene expression with molecular memory and bursting. Phys Rev Lett. 2011;106 (5 ): 058102.21405439
29. Zhang Z, Deng Q, Wang Z, Chen Y, Zhou T. Exact results for queuing models of stochastic transcription with memory and crosstalk. Phys Rev E. 2021;103 (6 ): 062414.34271765
30. Schwabe A, Rybakova KN, Bruggeman FJ. Transcription stochasticity of complex gene regulation models. Biophys J. 2012;103 (6):1152–1161.22995487
31. Nicolas D, Zoller B, Suter DM, Naef F. Modulation of transcriptional burst frequency by histone acetylation. Proc Natl Acad Sci USA. 2018;115 (27 ):7153–7158.29915087
32. Donovan BT, Huynh A, Ball DA, Patel HP, Poirier MG, Larson DR, Ferguson ML, Lenstra TL. Live-cell imaging reveals the interplay between transcription factors, nucleosomes, and bursting. EMBO J. 2019;38 (12): e100809.31101674
33. Zabidi MA, Stark A. Regulatory enhancer-core-promoter communication via transcription factors and cofactors. Trends Genet. 2016;32 (12 ):801–814.27816209
34. Voss TC, Hager GL. Dynamic regulation of transcriptional states by chromatin and transcription factors. Nat Rev Genet. 2014;15 (2 ):69–81.24342920
35. Suter DM. Transcription factors and DNA play hide and seek. Trends Cell Biol. 2020;30 (6 ):491–500.32413318
36. Brown CR, Mao CH, Falkovskaia E, Jurica MS, Boeger H. Linking stochastic fluctuations in chromatin structure and gene expression. PLOS Biol. 2013;11 (8 ): e1001621.23940458
37. Muramoto T, Müller I, Thomas G, Melvin A, Chubb JR. Methylation of H3K4 is required for inheritance of active transcriptional states. Curr Biol. 2010;20 (5 ):397–406.20188556
38. Casamassimi A, Napoli C. Mediator complexes and eukaryotic transcription regulation: An overview. Biochimie. 2007;89 (12 ):1439–1446.17870225
39. Malik S, Roeder RG. Dynamic regulation of pol II transcription by the mammalian mediator complex. Trends Biochem Sci. 2005;30 (5 ):256–263.15896744
40. Jonkers I, Kwak H, Lis JT. Genome-wide dynamics of pol II elongation and its interplay with promoter proximal pausing, chromatin, and exons. eLife. 2014;3 : e02407.24843027
41. Tome JM, Tippens ND, Lis JT. Single-molecule nascent RNA sequencing identifies regulatory domain architecture at promoters and enhancers. Nat Genet. 2018;50 (11):1533–1541.30349116
42. Shao W, Zeitlinger J. Paused RNA polymerase II inhibits new transcriptional initiation. Nat Genet. 2017;49 (7):1045–1051.28504701
43. Kumar N, Singh A, Kulkarni RV. Transcriptional bursting in gene expression: Analytical results for general stochastic models. PLOS Comput Biol. 2015;11 (10 ): e1004292.26474290
44. Suter DM, Molina N, Gatfield D, Schneider K, Schibler U, Naef F. Mammalian genes are transcribed with widely different bursting kinetics. Science. 2011;332 (6028 ):472–474.21415320
45. Hendy O, Campbell J Jr, Weissman JD, Larson DR, Singer DS. Differential context-specific impact of individual core promoter elements on transcriptional dynamics. Mol Biol Cell. 2017;28 (23 ):3360–3370.28931597
46. Li C, Cesbron F, Oehler M, Brunner M, Hӧfer T. Frequency modulation of transcriptional bursting enables sensitive and rapid gene regulation. Cell Syst. 2018;6 (4 ):409–423.e11.29454937
47. Zoller B, Little SC, Gregor T. Diverse spatial expression patterns emerge from unified kinetics of transcriptional bursting. Cell. 2018;175 (3):835–847.30340044
48. Garcia HG, Tikhonov M, Lin A, Gregor T. Quantitative imaging of transcription in living Drosophila embryos links polymerase activity to patterning. Curr Biol. 2013;23 (21 ):2140–2145.24139738
49. Hafner A, Reyes J, Stewart-Ornstein J, Tsabar M, Jambhekar A, Lahav G. Quantifying the central dogma in the p53 pathway in live single cells. Cell Syst. 2020;10 (6 ):495–505.32533938
50. Molina N, Suter DM, Cannavo R, Zoller B, Gotic I, Naef F. Stimulus-induced modulation of transcriptional bursting in a single mammalian gene. Proc Natl Acad Sci USA. 2013;110 (51 ):20563–20568.24297917
51. Lee CH, Shin H, Kimble J. Dynamics of Notch-dependent transcriptional bursting in its native context. Dev Cell. 2019;50 (4 ):426–435.e4.31378588
52. Larsson AJM, Johnsson P, Hagemann-Jensen M, Hartmanis L, Faridani OR, Reinius B, Segerstolpe Å, Rivera CM, Ren B, Sandberg R. Genomic encoding of transcriptional burst kinetics. Nature. 2019;565 (7738):251–254.30602787
53. Grün D, Kester L, van Oudenaarden A. Validation of noise models for single-cell transcriptomics. Nat Methods. 2014;11 (6):637–640.24747814
54. Grima R, Esmenjaud P-M. Systematic biases in transcriptional parameters inferred from single-cell snapshot data. bioRxiv. 2023;2023.06.19.545536.
55. Yunger S, Rosenfeld L, Garini Y, Shav-Tal Y. Single-allele analysis of transcription kinetics in living mammalian cells. Nat Methods. 2010;7 (8):631–633.20639867
56. Raj A, Peskin CS, Tranchina D, Vargas DY, Tyagi S. Stochastic mRNA synthesis in mammalian cells. PLOS Biol. 2006;4 (10 ): e309.17048983
57. Dar RD, Razooky BS, Singh A, Trimeloni TV, McCollum JM, Cox CD, Simpson ML, Weinberger LS. Transcriptional burst frequency and burst size are equally modulated across the human genome. Proc Natl Acad Sci USA. 2012;109 (43 ):17454–17459.23064634
58. Itzkovitz S, Lyubimova A, Blat IC, Maynard M, van Es J, Lees J, Jacks T, Clevers H, van Oudenaarden A. Single-molecule transcript counting of stem-cell markers in the mouse intestine. Nat Cell Biol. 2012;14 (1 ):106–114.
59. Little SC, Tikhonov M, Gregor T. Precise developmental gene expression arises from globally stochastic transcriptional activity. Cell. 2013;154 (4 ):789–800.23953111
60. Pare A, Lemons D, Kosman D, Beaver W, Freund Y, McGinnis W. Visualization of individual Scr mRNAs during Drosophila embryogenesis yields evidence for transcriptional bursting. Curr Biol. 2009;19 (23 ):2037–2042.19931455
61. Boettiger AN, Levine M. Rapid transcription fosters coordinate snail expression in the Drosophila embryo. Cell Rep. 2013;3 (1 ):8–15.23352665
62. Engl C, Jovanovic G, Brackston RD, Kotta-Loizou I, Buck M. The route to transcription initiation determines the mode of transcriptional bursting in E. coli. Nat Commun. 2020;11 (1 ):2422.32415118
63. Zaslaver A, Bren A, Ronen M, Itzkovitz S, Kikoin I, Shavit S, Liebermeister W, Surette MG, Alon U. A comprehensive library of fluorescent transcriptional reporters for Escherichia coli. Nat Methods. 2006;3 (8 ):623–628.16862137
64. Valen E, Sandelin A. Genomic and chromatin signals underlying transcription start-site selection. Trends Genet. 2011;27 (11 ):475–485.21924514
65. Lenhard B, Sandelin A, Carninci P. Metazoan promoters: Emerging characteristics and insights into transcriptional regulation. Nat Rev Genet. 2012;13 (4):233–245.22392219
66. Muse GW, Gilchrist DA, Nechaev S, Shah R, Parker JS, Grissom SF, Zeitlinger J, Adelman K. RNA polymerase is poised for activation across the genome. Nat Genet. 2007;39 (12 ):1507–1511.17994021
67. Weake VM, Workman JL. Inducible gene expression: Diverse regulatory mechanisms. Nat Rev Genet. 2010;11 (6 ):426–437.20421872
68. Sanchez A, Golding I. Genetic determinants and cellular constraints in noisy gene expression. Science. 2013;342 (6163 ):1188–1193.24311680
69. Battich N, Stoeger T, Pelkmans L. Control of transcript variability in single mammalian cells. Cell. 2015;163 (7 ):1596–1610.26687353
70. Lammers NC, Galstyan V, Reimer A, Medin SA, Wiggins CH, Garcia HG. Multimodal transcriptional control of pattern formation in embryonic development. Proc Natl Acad Sci USA. 2020;117 (2 ):836–847.31882445
71. Swain PS, Elowitz MB, Siggia ED. Intrinsic and extrinsic contributions to stochasticity in gene expression. Proc Natl Acad Sci USA. 2002;99 (20 ):12795–12800.12237400
72. Cao Z, Grima R. Analytical distributions for detailed models of stochastic gene expression in eukaryotic cells. Proc Natl Acad Sci USA. 2020;117 (9 ):4682–4692.32071224
73. Fu X, Patel HP, Coppola S, Xu L, Cao Z, Lenstra TL, Grima R. Quantifying how post-transcriptional noise and gene copy number variation bias transcriptional parameter inference from mRNA distributions. eLife. 2022;11 : e82493.36250630
74. Zhang Z, English BP, Grimm JB, Kazane SA, Hu W, Tsai A, Inouye C, You C, Piehler J, Schultz PG, et al. Rapid dynamics of general transcription factor TFIIB binding during preinitiation complex assembly revealed by single-molecule analysis. Genes Dev. 2016;30 (18 ):2106–2118.27798851
75. Roeder RG. 50+ years of eukaryotic transcription: An expanding universe of factors and mechanisms. Nat Struct Mol Biol. 2019;26 (9 ):783–791.31439941
76. Yakovchuk P, Gilman B, Goodrich JA, Kugel JF. RNA polymerase II and TAFs undergo a slow isomerization after the polymerase is recruited to promoter-bound TFIID. J Mol Biol. 2010;397 (1 ):57–68.20083121
77. Warfield L, Ramachandran S, Baptista T, Devys D, Tora L, Hahn S. Transcription of nearly all yeast RNA polymerase II-transcribed genes is dependent on transcription factor TFIID. Mol Cell. 2017;68 (1 ):118–129.e5.28918900
78. Blair RH, Goodrich JA, Kugel JF. Single-molecule fluorescence resonance energy transfer shows uniformity in TATA binding protein-induced DNA bending and heterogeneity in bending kinetics. Biochemistry. 2012;51 (38 ):7444–7455.22934924
79. Bartman CR, Hsu SC, Hsiung CCS, Raj A, Blobel GA. Enhancer regulation of transcriptional bursting parameters revealed by forced chromatin looping. Mol Cell. 2016;62 (2 ):237–247.27067601
80. Chen H, Levo M, Barinov L, Fujioka M, Jaynes JB, Gregor T. Dynamic interplay between enhancer-promoter topology and gene activity. Nat Genet. 2018;50 (9 ):1296–1303.30038397
81. Saunders A, Core LJ, Lis JT. Breaking barriers to transcription elongation. Nat Rev Mol Cell Biol. 2006;7 (8 ):557–567.16936696
