Abstract

Gene transcription is a random process in single cells manifested by the observed distribution of mRNA copy numbers in homogeneous cell populations. A central question is to understand how mRNA distribution is modulated under environmental changes. In this work, we initiate a theoretical study on mRNA distribution dynamics for the stochastic transcription model that involves cross-talking signaling pathways to direct gene activation in response to external signals. We first express the distribution in mathematical dynamical formulas under both moderate and high transcriptional upregulations. In each scenario, our further numerical examples display an observed dynamical transition type among three distribution modes for stress genes in yeast. In particular, the intermediate bimodal stage sustains within a certain length of early time and lasts much longer than that generated by the single pathway. This shows the general and robust bimodal transcription regulated by the cross-talk of signaling pathways.

1. Introduction

Recent single-cell measurements have generated massive data on the histogram of mRNA copy numbers for the target gene in homogeneous cell populations [1, 2]. This provides a good approximation for data fitting by the probability mass function , the probability that there are exactly mRNA molecules of the gene of our concern at time in one cell  [3, 4]. The distribution profile of contains a panoramic information for distinct cellular phenotypes  [5, 6]. The commonly observed modes are the decaying distribution that decreases in , the unimodal distribution that peaks uniquely at some , and the bimodal distribution that takes exactly two peaks with the first one at and the other one at some . A decaying or unimodal distribution suggests a phenotypic homogeneity, while a bimodal distribution supports a binary process that steers cells into subpopulations with distinct cell identities  [1, 7, 8].

It has been a central question to understand how varying external signals influence the profile of for inducible genes  [1, 2]. As suggested by the prevailing two-state model  [1, 9], it may involve intrinsic mechanisms that regulate a random switching between gene-off (inactive) and gene-on (active) states  [6, 10, 11]. As shown in the following equation:the dwell times in the off and on states are independent and exponentially distributed with the activation rate and inactivation rate , respectively. In an active gene, RNA polymerase can bind to the promoter and traverses the template DNA strand to direct mRNA synthesis. The waiting time for the birth of a new mRNA molecule and its lifetime before death are independent and exponentially distributed with the synthesis rate and degradation rate , respectively.

The two-state model (1) has emerged as a standard framework for quantifying the deviation of the mRNA level in individual cells from the mean, through calculating the noise (variance over mean squared), fano factor (variance over mean), and probability distribution of random mRNA copy numbers  [1, 2, 8]. These indexes have been expressed as analytical formulas for the two-state model or its more elaborated extension regarding a larger degree of biological realism such as chromatin remodeling, mRNA maturation, and cell division  [6, 1214]. Together with experimental data, those formulas can be reversed to the variation of system parameters and, in turn, uncover a large spectrum of regulation modes that cells utilize under different cellular environments  [2, 4, 1416].

On the contrary, the two-state model (1) implicitly assumes that the gene activation from the off state to the on state is directed by a single signaling pathway. Such assumption may not be justified for a large class of inducible genes which maintain a low transcription level induced by the spontaneous basal pathway under normal cellular growth conditions [17, 18], while upregulating the transcription through specific signaling pathways in response to acute external changes  [19, 20]. Also, the activation of other genes that are involved in stem cell renewal, development, and immunity is often mediated by two signal transduction pathways  [2123]. In all these cases, the downstream specific transcription factors (TFs) in one pathway compete with the other TFs in another pathway for binding at their shared promoter sites to direct the formation of the intermediate complex  [19, 24]. Therefore, the two competitive cross-talking pathways could be generated to activate the target gene.

By integrating competitive cross-talking pathways into the gene activation process  [24, 25], the two-state model (1) can be generalized as the following equation:

The transcription of the gene can be either induced by the first pathway, or alternatively, the second pathway. The sojourn times in the two pathways are independent and exponentially distributed, with the induction strength for the first pathway and the strength for the other, satisfying

We rename the gene-off state as if it is transformed to the gene-on state by the first pathway, and as otherwise. Then, the pathway selection probabilities, denoted bysatisfy

Also, the gene inactivation from the state to the state and the birth and the death of mRNA molecules are separately controlled by the first-order kinetic rates , , and , as modeled in the two-state model (1).

When , the exact form for the stationary mass function was stated in [25]. The subsequent numerical examples showed that cross-talking pathways are more likely to generate bimodal distribution compared to the single pathway with or . This naturally gives rise to an interesting question of whether the cross-talking regulation still maintains the robust bimodal distribution as time develops. When time is finite, however, the analytical formulas of have remained elusive due to its intrinsic complexity. Only a few papers considered the case of a single pathway and expressed or the corresponding generating function in the form of integrals  [12], infinite series  [26, 27], and hypergeometric functions  [6, 28] under certain parameter regions. We shall derive generated by cross-talking pathways in simple mathematical formulas in Section 2, and then discuss their dynamical profiles and implications in Section 3. The main conclusion and its discussion are given in the last section.

2. Analytical Formulas for Dynamical Distribution

In this section, we have endeavored to calculate by solving the system of master equations. The standard approach is to transform it into a system of first-order partial differential equations through the method of generating functions  [6, 12]. The analytical form of has been found to be the solutions for several special types of third-order ordinary differential equations. Solving their corresponding initial value problems and then extracting are in general formidable and pose a major obstacle in the study.

We have successfully derived many exact forms of within a large range of system parameters. In particular, our first result assumes that the weaker signaling pathway is more frequently selected. In this case, the upregulation of the target gene may not be efficient due to the exposure of cells to mild external signals or the intrinsic mechanism that avoids overexuberant transcription of the target gene  [17, 20, 29].

Theorem 1. If , , and , then the probability mass function in the cross-talking pathway model can be expressed asand for ,or alternatively, for ,

Proof. Let , be the respective probabilities of mRNA molecules at time that the gene is residing in the on state and the off state with the th pathway being selected. Then, the total probability mass function isFollowing the standard procedure in the theory of stochastic processes [13, 24], we find that the three partial mass functions in the cross-talking pathway model satisfy the following system of master equations:where, by convention, we set . Without loss of generality, we assume that the gene is inactive and the number of transcripts is zero at with the initial conditions:Following the standard procedure  [6, 12], we introduce the probability generating functionsto transform the initial value problem for the master equations (10)–(12) into the problem for the system of first-order partial differential equations:If the initial value problem (15)–(18) can be solved analytically, then we may obtain by adding the three solutions together and then applying the conversion formula: We continue from there and transform (15)–(17) into ordinary differential equations through the method of characteristics [12, 28, 30]. Let be a parameter. Letbe the restriction of , , , and on the characteristic curve . Then,By introducing nondimensionalized system parameters,and from the further transformation,we find , and transform (21)–(23) intoIt is interesting to note that equations (26) and (27) are indeed identical with equations (17)–(19) in  [25]. Therefore, by introducing two real numbers determined by algebra equations  [24],the same calculation verifying (21) in  [25] shows that is the unique solution of the initial value problem:for the linear operator given bywith real constants , , , and and a smooth function of .
To proceed, we define byand utilize notations (24) and (29) and the assumption of theorem to verifyWe then reformulate of (30) in the formThe substitution of (32) and (33) into the above equation gives , which readily converts (30) to the simple problem:The unique solution of (35) is . We replace its left side by (32) and the right side by the initial condition in (35). This gives an inhomogeneous equation:subject to the initial condition given in (30)Taking into account that the homogeneous counterpart of (36) is Euler’s equation which possesses two independent solutions and , their linear combination, adding a particular solution of (36), constitutes the general solution of (36). Following the standard method of undetermined coefficients in the theory of ordinary differential equation  [31], we find that one particular solution takes the following form:By fixing the coefficients in the linear combination through the initial condition (37), we obtain the unique solution of (36) and (37) in the following form:Now, the transformations defined in (20) and defined in (25) convert (39) back toBy making the change of variable in the integral, we rewrite asIn turn, by using conversion and then letting in the integral, we obtainWe divide the rest of the proof into two parts. We derive (6) and (7) in the first part and (8) in the second part.
Derivation of (6) and (7): we first note that the second integral in (42) can be expressed in an equivalent form through integration by parts:which helps us rewrite the generating function asTo finally verify (6) and (7), we extract the probability mass function from (44) through the conversion formula (19). Since the first two terms in (44) are eliminated by differentiation with respect to , in (6) needs to be derived from (44) directly as . For ,By substituting and into this expression and dividing it by , we obtain (7).
Derivation of (8): by changing the first integral in (42) in the formwe rewrite (42) asWith the help of the elementary identitywe findThen, (8) follows from (19) immediately.
Next, we present a result for the case , for which gene activation is irreversible  [26]. It may give an approximation to the case when the gene-on state is regulated to be very stable () under the highest induction level  [15, 16, 32].

Theorem 2. If , then the probability mass function takes the formor alternatively, for ,

Proof. As in the proof of Theorem 1, our basic strategy is to solve the initial value problem (30). If , then and by definitions (24) and (29). Combining this with (31), we see that reduces toAccording to the initial conditions in (30), we denote byThen, (52) is Euler’s equation of with two independent solutions and . Hence, can be expressed as a linear combination of and . By determining the coefficients through initial conditions of , it readily follows thatTogether with (53), we find that can be obtained by solving the simple problem:By a routine calculation of the first-order linear ordinary differential equation  [31], we obtainBy replacing with and with , we convert to the generating function :Substituting in the integral simplifies the expression toAfter substituting and then , we finally obtainDerivation of (50): through integration by parts, we express the last integral of (59) in an equivalent form:The substitution of this into (59) helps us rewrite the generating function asBy the conversion formula (19), in (50) is derived from (61) directly as . For ,which leads to (50) by substituting and into this expression and dividing it by .
Derivation of (51): with the help of (59) and the elementary identity (48), we findThen, (51) follows from the conversion formula (19) immediately. The proof is completed.

3. Dynamical Bimodal Distribution for Inducible Genes

We elucidate the regulation of cross-talking pathways on dynamical mRNA distribution for inducible genes under acute stresses, which have to fulfill two requirements. First, a default or spontaneous mechanism maintains transcription at the basal level when cells are under normal cellular environments to avoid the overexuberant expression  [17, 29]. Second, the exposure of cells to the extracellular stresses results in the transient activation of dedicated intracellular signaling pathways, which enables cells to adapt immediately to the new environment for maximal cell survival  [19, 20]. For instance, the osmotic shock induces the MAPK HOG pathway in S. cerevisiae through the phosphorylation of MAPK Hog1 that leads to the rapid translocation of Hog1 into the nucleus to launch a transcriptional program  [19]. Therefore, the cross-talk between a weak basal pathway and an inducible signaling pathway could be generated to initiate the stress gene transcription.

We shall discuss the observed dynamical transition among three distribution modes for the STL1 gene in S. cerevisiae in response to acute osmotic stress  [19]. Under mild or high NaCl concentrations, no expression was detected in most cells at the initial time referred to as the decaying distribution and full expression was detected at a large time scale referred to as the unimodal distribution. The scenario becomes more fascinated within the intermediate time as both nonexpressing and fully expressing cell subpopulations were presented, characterized as the bimodal distribution. It is believed that such bimodal transcription output is controlled by both the retention time and concentration of Hog1 in the nucleus and results in a larger phenotypic variability for cell survival in harsh environments [19].

We demonstrate through numerical examples that cross-talking pathways can generate dynamical bimodality of osmoresponsive gene transcription under mild or high stress level. Most mRNA half-lives in S. cerevisiae range around a median of 11 min  [33], so we fix . We set the transcription rate as its upper bound   [34], and thus . Also, we set a relative small strength rate for the basal pathway and a relative large strength rate for the signaling pathway. Such parameter set generates a prolonged region of the gene-off time (denoted by ):that covers the lag time of 1 to 3 minutes for the activation of most typical osmoresponsive genes, such as STL1, GRE2, or GPD1  [20].

The mild stress level results in a moderate increasing nuclear accumulation of activated Hog1  [19]. This leads to the well-matched competition between two pathways with their selected probabilities . We set and that generates transcriptional efficiency. On the other hand, the high stress level induces saturated nuclear Hog1 concentration  [19] and thus may direct gene activation more frequently through the signaling pathway with . We then set and that generates . Also, it was reported that the gene-on state could be significantly stabilized in response to increasing stress level, and thus   [2, 15, 16, 32].

We next take advantage of exact forms (8) and (51) to quantify dynamical mRNA distribution under mild and high stress levels, respectively. In both scenarios, we observe obvious transitions from the decaying distribution to unimodal distribution as time develops, going through a bimodal stage within a 30 min length of early time (Figure 1). To test to what extent the cross-talking pathways modulate the bimodality, we reset so that the gene is activated by a single pathway with strength rate . In turn, we let be equal to transcriptional efficiencies for both cases of mild and high stress levels. Interestingly, it shows that the intermediate bimodality becomes much more fragile that lasts less than 5 min, see Figure 1 (inset). This confirms that the cross-talking regulation indeed prolongs duration of mRNA bimodal distribution, reinforcing the general feature of bimodal transcription for stress genes  [19, 35].

4. Conclusion and Discussion

Recent single-cell studies have generated massive data on the distribution of mRNA copy numbers for many genes under various experimental conditions  [13]. The most commonly observed modes for the distribution are the decaying distribution, unimodal distribution, and bimodal distribution. It is reported that the bimodal distribution could help differentiate an isogenic cell population into two dynamically stable groups with distinct phenotypes. For instance, the bimodal Nanog expression plays a key role in mediating cell differentiation in embryonic stem cells  [7]; the transcriptional bimodality of generates one HIV highly expressed subpopulation and the other with HIV latency  [36]; the bimodal transcription pattern of stress genes in yeast increases phenotypic variability in response to unpredictable environmental changes  [19]. In the classical two-state model (1), the bimodal distribution has been suggested to originate from the random switch between gene-off and gene-on states as bimodality would disappear when the gene is always active  [7, 28].

One of our major concerns is how the gene activation process modulates the bimodal distribution for inducible genes in face of external cues. We introduce a cross-talking pathway model (2) by integrating two competitive signaling pathways into the two-state model  [24, 25]. Comparing to other statistical concepts such as the noise or fano factor that can be computed relatively easily, exploring transparent analytical formulas for mRNA distribution remains a challenging task. For instance, a traditional approach for calculating mRNA distribution is to express the corresponding generating function in the form of hypergeometric functions  [6, 28]. However, this does not permit an easy tractable way to further transform the generating function to mRNA distribution when the time and system parameters are finite. We endeavor to express time-dependent mRNA distribution in simple mathematical formulas within a large range of finite parameters for the cross-talking pathway model. The formulas are obtained by solving the initial value problems of third-order ordinary differential equations, which are derived from the master equations of the cross-talking pathway model.

We consider stress-responsive genes under acute stresses as the implement of dynamical analytical formulas for mRNA distribution. The stress gene can be activated either by a weak spontaneous basal pathway with strength and selection probability , or alternately by an inducible signaling pathway with strength and selection probability . The probability (or ) may be governed by the number and availability of activated transcription factors (TFs) near the binding sites  [19, 20] and thus is related with the stress level; the inducible strength can be governed by the binding strength between the TFs and the DNA sites [11, 20], and thus may be tightly related with the type of external signals. We suggest that the signaling pathway is efficiently induced with

On the assumption of Theorem 1, (65) leads to , which shows a mild stress level that may induce a similar concentration of activated TFs in each pathway. The assumption of Theorem 2 presents an active system in which the gene inactivation process is blocked. It suggests a very high stress level that may induce a significant decrease of the gene inactivation rate   [2, 16, 32]. For instance, compared with normal cellular conditions, the inductions with the highest strength lead to the decrease of for over folds for the Lac promoter in E. coli  [32] and possible fold change for the POL1 promoter in yeast  [34].

Under both mild and high stress levels, we utilize the parameter rates for stress genes in yeast to understand their dynamical distributions. We reveal a typical dynamical transition from the decaying distribution to unimodal distribution as time develops, going through an intermediate bimodal stage with over 30 minutes long. This result strongly suggests the general and robust bimodal transcription of stress genes regulated by the cross-talking pathways in response to acute stresses  [19, 35]. First, noting that the system takes about 120 minutes to reach the steady state, bimodal distribution maybe clearly visible in the course of real-time imaging due to its relatively large time window of 30 minutes. Second, under the same durations for the gene-off and gene-on states, as well as the same mRNA synthesis and degradation rates, the two-state model (a single pathway) can only generate a much fragile intermediate bimodal stage that lasts less than 5 minutes.

The general intermediate bimodal stage is consistent with recent observations that the cross-talking pathway model (2) is more likely to generate high transcription noise and bimodal distribution at the steady state compared with the two-state model (1) [24, 25]. It suggests that natural selection may favor two or more parallel signaling pathways for gene activation to enhance expression variability which may provide a clear benefit in the face of an environmental stress  [19, 37]. On the other hand, the cross-talking pathway model (2) assumes the constant system rates and cannot take into account some inducible genes and immune genes regulated by the oscillating signal  [38], temporal feedback  [39], or different cell cycle phases  [14]. Interestingly, the dynamical bimodality can be easily generated through the other very different mechanisms such as the system involves positive or negative genetic feedback loops  [39]. The extensions of our model that include those salient biological features may require assuming time-varying parameters with biological realism accordingly. Future work will focus on how those extended models modulate the dynamical bimodal or even multimodal transcription distributions.

Data Availability

The data used to support findings of this study are included within the article: (1) the degradation rate of the mRNA molecule in S. cerevisiae equals ln(2)/11 ≈ 0.063 (1/min), which is calculated through the 11 min average half-life obtained in reference [33]. (2) The mRNA synthesis rate in S. cerevisiae equals 0.063 ∗ 10=0.63 (1/min), which is calculated according to the mRNA degradation rate 0.063 (1/min) estimated in (1) and the upper bound of its transcription rate 10 (per mRNA lifetime) given in reference [34]. (3) The 1–3 min gene-off duration for most typical osmoresponsive genes in S. cerevisiae is reported in reference [20].

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This work was supported by the Natural Science Foundation of China grants (nos. 11631005 and 11871174), Program for Changjiang Scholars and Innovative Research Team in University (no. IRT_16R16), Key Project of Hunan Educational Committee (no. 19A497), and China Postdoctoral Science Foundation (Grant no. 2019M660196).