# Dynamical Regulation of mRNA Distribution by Cross-Talking Signaling Pathways

Authors: Chunjuan Zhu; Guosheng Han; Feng Jiao
Journal: Complexity
 (2020)

---

## 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.


---

## Body

## 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 Pmt, the probability that there are exactly m mRNA molecules of the gene of our concern at time t in one cell  [3, 4]. The distribution profile of Pmt contains a panoramic information for distinct cellular phenotypes  [5, 6]. The commonly observed modes are the decaying distribution that Pmt decreases in m, the unimodal distribution that Pmt peaks uniquely at some m>0, and the bimodal distribution that Pmt takes exactly two peaks with the first one at m=0 and the other one at some m>0. 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 ofPmt 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:(1)gene off⇄γλgene on⟶vmRNA⟶δ∅,the dwell times in the off and on states are independent and exponentially distributed with the activation rate λ>0 and inactivation rate γ>0, 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 v>0 and degradation rate δ>0, 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 Pmt 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, 12–14]. 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, 14–16].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  [21–23]. 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:(2)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λ1>0 for the first pathway and the strength λ2>0 for the other, satisfying(3)0<λ1≤λ2<∞.We rename the gene-off state asO1 if it is transformed to the gene-on state by the first pathway, and as O2 otherwise. Then, the pathway selection probabilities, denoted by(4)q1=ProbO=O1,q2=ProbO=O2,satisfy(5)0<q1,q2<1,q1+q2=1.Also, the gene inactivation from theon state to the off state and the birth and the death of mRNA molecules are separately controlled by the first-order kinetic rates γ, v, and δ, as modeled in the two-state model (1).Whent⟶∞, the exact form for the stationary mass function limt⟶∞Pmt 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 q1=1 or λ1=λ2. 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 t is finite, however, the analytical formulas of Pmt have remained elusive due to its intrinsic complexity. Only a few papers considered the case of a single pathway and expressed Pmt 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 Pmt 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 calculatePmt 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 Pmt 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 Pmt are in general formidable and pose a major obstacle in the study.We have successfully derived many exact forms ofPmt 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γ=2δ, q1>q2, and λ2=λ1+δ/q1−q2, then the probability mass function Pmt in the cross-talking pathway model can be expressed as(6)P0t=q1e−λ1t+q2e−λ2t+λ1λ2λ2−λ1∫0te−λ1s−e−λ2seve−δt−e−δs/δds+q1λ1+q2λ2λ2−λ1∫0tδ−λ1e−λ1s−δ−λ2e−λ2seδs−teve−δt−e−δs/δds,and for m=1,2,…,(7)Pmt=vmq1λ1+q2λ2δmm!λ2−λ1λ1λ2q1λ1+q2λ2∫0te−λ1s−e−λ2se−δs−e−δtmeve−δt−e−δs/δds+∫0tδ−λ1e−λ1s−δ−λ2e−λ2seδs−te−δs−e−δtmeve−δt−e−δs/δds,or alternatively, for m=0,1,…,(8)Pmt=vmδmm!λ2−λ1∫0tq1λ1+q2λ2e−λ1s−e−λ2se−δt−λ2e0−λ1+δs+λ1e−λ2+δs×mδ+ve−δt−e−δse−δs−e−δtm−1×eve−δt−e−δs/δds+vm1−e−δtmδmm!ev1−e−δt/δ.Proof.
LetPm,jt,j=e,1,2, be the respective probabilities of m mRNA molecules at time t that the gene is residing in the on state and the off state with the jth pathway being selected. Then, the total probability mass function is(9)Pmt=Pm,1t+Pm,2t+Pm,et.
Following 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:(10)Pm,1′t=q1γPm,et−mδ+λ1Pm,1t+m+1δPm+1,1t,(11)Pm,2′t=q2γPm,et−mδ+λ2Pm,2t+m+1δPm+1,2t,(12)Pm,e′t=λ1Pm,1t+λ2Pm,2t−v+mδ+γPm,et+m+1δPm+1,et+vPm−1,et,where, by convention, we set P−1,et=0. Without loss of generality, we assume that the gene is inactive and the number of transcripts is zero at t=0 with the initial conditions:(13)P0,10=q1,P0,20=q2,P0,e0=Pm,i0=0,fori=1,2,e,m≥1.
Following the standard procedure  [6, 12], we introduce the probability generating functions(14)Viz,t=∑m=0∞z+1mPm,it,i=1,2,e,to transform the initial value problem for the master equations (10)–(12) into the problem for the system of first-order partial differential equations:(15)∂V1∂tz,t=−λ1V1z,t+q1γVez,t−δz∂V1∂zz,t,(16)∂V2∂tz,t=−λ2V2z,t+q2γVez,t−δz∂V2∂zz,t,(17)∂Ve∂tz,t=λ1V1z,t+λ2V2z,t−γVez,t+vzVez,t−δz∂Ve∂zz,t,(18)V1z,0=q1,V2z,0=q2,Vez,0=0.
If the initial value problem (15)–(18) can be solved analytically, then we may obtain Pmt by adding the three solutions together and then applying the conversion formula: (19)Pmt=1m!∂zVz,t∂zmz=−1,whereVz,t=V1z,t+V2z,t+Vez,t,m≥0.
We continue from there and transform (15)–(17) into ordinary differential equations through the method of characteristics [12, 28, 30]. Let z0 be a parameter. Let(20)uit=Viz0eδt,t,fori=1,2,e,ut=u1t+u2t+uet,be the restriction of V1, V2, Ve, and V on the characteristic curve z=z0eδt. Then,(21)u1′t=−λ1u1t+q1γuet,u2′t=−λ2u2t+q2γuet,(22)ue′t=λ1u1t+λ2u2t−γuet+vz0eδtuet,(23)u10=q1,u20=q2,ue0=0.
By introducing nondimensionalized system parameters,(24)λ¯1=λ1δ,λ¯2=λ2δ,γ¯=γδ,v¯=vδ,and from the further transformation,(25)wix=uit,fori=1,2,e,wx=ut,wherex=v¯z0eδt,we find ui′t=δxwi′x, and transform (21)–(23) into(26)xw1′x=−λ¯1w1x+q1γ¯wex,xw2′x=−λ¯2w2x+q2γ¯wex,(27)xwe′x=λ¯1w1x+λ¯2w2x−γ¯wex+xwex,(28)w1v¯z0=q1,w2v¯z0=q2,wev¯z0=0.
It is interesting to note that equations (26) and (27) are indeed identical with equations (17)–(19) in  [25]. Therefore, by introducing two real numbers α¯>β¯>0 determined by algebra equations  [24],(29)α¯+β¯=λ¯1+λ¯2+γ¯,α¯β¯=λ¯1λ¯2+q1λ¯2γ¯+q2λ¯1γ¯,the same calculation verifying (21) in  [25] shows that wx is the unique solution of the initial value problem:(30)ℒα¯,β¯λ¯1,λ¯2w=0,wv¯z0=1,w′v¯z0=0,w″v¯z0=q1λ¯1+q2λ¯2v¯z0,for the linear operator ℒ given by(31)ℒc,da,bf=x2f‴x+x1−x+c+df″x+cd−xa+b+1f′x−abfx,with real constants a, b, c, and d and a smooth function f=fx of x.
To proceed, we define by(32)yx=x2w″x+xλ¯1+λ¯2+1w′x+λ¯1λ¯2wx,and utilize notations (24) and (29) and the assumption of theorem to verify(33)β¯=λ¯1+1,α¯=λ¯2+1.
We then reformulateℒα¯,β¯λ¯1,λ¯2w=0 of (30) in the form(34)x2w‴x+x1+α¯+β¯w″x+α¯β¯w′x=x2w″x+xλ¯1+λ¯2+1w′x+λ¯1λ¯2wx.
The substitution of (32) and (33) into the above equation gives y′x=yx, which readily converts (30) to the simple problem:(35)y′x=yx,yv¯z0=λ¯1λ¯2+v¯z0q1λ¯1+q2λ¯2.
The unique solution of (35) is y=yv¯z0expx−v¯z0. We replace its left side by (32) and the right side by the initial condition in (35). This gives an inhomogeneous equation:(36)x2w″x+xλ¯1+λ¯2+1w′x+λ¯1λ¯2wx=λ¯1λ¯2+v¯z0q1λ¯1+q2λ¯2ex−v¯z0,subject to the initial condition given in (30)(37)wv¯z0=1,w′v¯z0=0.
Taking into account that the homogeneous counterpart of (36) is Euler’s equation which possesses two independent solutions x−λ¯1 and x−λ¯2, 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:(38)w¯x=λ¯1λ¯2+v¯z0q1λ¯1+q2λ¯2λ¯2−λ¯1∫v¯z0xsλ¯1−1xλ¯1−sλ¯2−1xλ¯2es−v¯z0ds.
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:(39)wx=C1x−λ¯1+C2x−λ¯2+w¯x,C1=λ2v¯z0λ¯1λ2−λ1,C2=−λ1v¯z0λ¯2λ2−λ1.
Now, the transformationsut=Vz0eδt,t defined in (20) and ut=wv¯z0eδt defined in (25) convert (39) back to(40)Vz0eδt,t=λ2e−λ1tλ2−λ1−λ1e−λ2tλ2−λ1+λ¯1λ¯2λ¯2−λ¯1∫v¯z0v¯z0eδtsλ¯1−1e−λ1tv¯z0λ¯1−sλ¯2−1e−λ2tv¯z0λ¯2es−v¯z0ds+v¯z0q1λ¯1+q2λ¯2λ¯2−λ¯1∫v¯z0v¯z0eδtsλ¯1−1e−λ1tv¯z0λ¯1−sλ¯2−1e−λ2tv¯z0λ¯2es−v¯z0ds.
By making the change of variables=v¯z0expδτ in the integral, we rewrite V as(41)Vz0eδt,t=λ2e−λ1tλ2−λ1−λ1e−λ2tλ2−λ1+λ1λ2λ2−λ1∫0teλ1τ−t−eλ2τ−tev¯z0eδτ−1dτ+vz0q1λ1+q2λ2λ2−λ1∫0teλ1τ−t−eλ2τ−tev¯z0eδτ−1dτ.
In turn, by using conversionz0=ze−δt and then letting s=t−τ in the integral, we obtain(42)Vz,t=λ2e−λ1tλ2−λ1−λ1e−λ2tλ2−λ1+λ1λ2λ2−λ1∫0te−λ1s−e−λ2sev¯ze−δs−e−δtds+vze−δtq1λ1+q2λ2λ2−λ1∫0te−λ1s−e−λ2sev¯ze−δs−e−δtds.
We 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:(43)vze−δt∫0te−λ1s−e−λ2sev¯ze−δs−e−δtds=e−δt∫0teδ−λ2s−eδ−λ1sdev¯ze−δs−e−δt=e−λ2t−e−λ1t−∫0tδ−λ2e−λ2s−δ−λ1e−λ1seδs−tev¯ze−δs−e−δtds,which helps us rewrite the generating function Vz,t as(44)Vz,t=q1e−λ1t+q2e−λ2t+λ1λ2λ2−λ1∫0te−λ1s−e−λ2sev¯ze−δs−e−δtds+q1λ1+q2λ2λ2−λ1∫0tδ−λ1e−λ1s−δ−λ2e−λ2seδs−tev¯ze−δs−e−δtds.
To finally verify (6) and (7), we extract the probability mass function Pmt from (44) through the conversion formula (19). Since the first two terms in (44) are eliminated by differentiation with respect to z, P0t in (6) needs to be derived from (44) directly as P0t=V−1,t. For m≥1,(45)∂mVz,t∂zm=v¯mq1λ1+q2λ2λ2−λ1λ1λ2q1λ1+q2λ2∫0te−λ1s−e−λ2se−δs−e−δtmev¯ze−δs−e−δtds+∫0tδ−λ1e−λ1s−δ−λ2e−λ2seδs−te−δs−e−δtmev¯ze−δs−e−δtds.
By substitutingv¯=v/δ and z=−1 into this expression and dividing it by m!, we obtain (7).
Derivation of (8): by changing the first integral in (42) in the form(46)λ1λ2∫0te−λ1s−e−λ2sev¯ze−δs−e−δtds=∫0tev¯ze−δs−e−δtdλ1e−λ2s−λ2e−λ1s=λ1e−λ2t−λ2e−λ1t+λ2−λ1ev¯z1−e−δt−vz∫0tλ2e−λ1+δs−λ1e−λ2+δsev¯ze−δs−e−δtds,we rewrite (42) as(47)Vz,t=ev¯z1−e−δt−vzλ2−λ1∫0tλ2e−λ1+δs−λ1e−λ2+δsev¯ze−δs−e−δtds+vze−δtq1λ1+q2λ2λ2−λ1∫0te−λ1s−e−λ2sev¯ze−δs−e−δtds.
With the help of the elementary identity(48)∂mzeaz∂zm=m+azam−1eaz,we find(49)∂mVz,t∂zm=v¯m1−e−δtmev¯z1−e−δt+v¯mλ2−λ1∫0tq1λ1+q2λ2e−λ1s−e−λ2se−δt−λ2e−λ1+δs+λ1e−λ2+δs×mδ+vze−δs−e−δte−δs−e−δtm−1ev¯ze−δs−e−δtds.
Then, (8) follows from (19) immediately.
Next, we present a result for the caseγ=0, 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 (γ≪1) under the highest induction level  [15, 16, 32].Theorem 2.
Ifγ=0, then the probability mass function Pmt takes the form(50)Pmt=q1e−λ1t+q2e−λ2t+∫0tq1λ1eλ1s−t+q2λ2eλ2s−tev¯e−δs−1ds,form=0,v¯mm!∫0tq1λ1eλ1s−t+q2λ2eλ2s−t1−e−δsmev¯e−δs−1ds,form=1,2,…,or alternatively, for m=0,1,…,(51)Pmt=v¯mm!1−e−δtmev¯e−δt−1−e−δt∫0tq1eδ−λ1t−s+q2eδ−λ2t−s×mδ−v1−e−δs1−e−δsm−1ev¯e−δs−1ds.Proof.
As in the proof of Theorem1, our basic strategy is to solve the initial value problem (30). If γ=0, then α=λ2 and β=λ1 by definitions (24) and (29). Combining this with (31), we see that ℒα¯,β¯λ¯1,λ¯2w=0 reduces to(52)x2w‴x−w″x+xλ¯1+λ¯2+1w″x−w′x+λ¯1λ¯2w′x−wx=0.
According to the initial conditions in (30), we denote by(53)y¯x=w′x−wx,withy¯′v¯z0=−1,y¯″v¯z0=q1λ¯1+q2λ¯2v¯z0.
Then, (52) is Euler’s equation of y¯x with two independent solutions y¯1x=x−λ¯1 and y¯2x=x−λ¯2. Hence, y¯x can be expressed as a linear combination of y¯1x and y¯2x. By determining the coefficients through initial conditions of y¯x, it readily follows that(54)y¯x=−q1v¯z0λ¯1x−λ¯1−q2v¯z0λ¯2x−λ¯2.
Together with (53), we find that wx can be obtained by solving the simple problem:(55)w′x−wx=−q1v¯z0λ¯1x−λ¯1−q2v¯z0λ¯2x−λ¯2,wv¯z0=1.
By a routine calculation of the first-order linear ordinary differential equation  [31], we obtain(56)wx=ex−v¯z0−∫v¯z0xq1v¯z0λ¯1s−λ¯1+q2v¯z0λ¯2s−λ¯2ex−sds.
By replacingx with v¯z and z0 with ze−δt, we convert wx to the generating function Vz,t:(57)Vz,t=ev¯z0eδt−1−∫v¯z0v¯z0eδtq1v¯z0λ¯1s−λ¯1+q2v¯z0λ¯2s−λ¯2ev¯z0eδt−sds.
Substitutings=v¯z0expδτ in the integral simplifies the expression to(58)Vz,t=ev¯z0eδt−1−vz0∫0tq1e−λ1τ+q2e−λ2τeδτev¯z0eδt−eδτdτ.
After substitutingz0=ze−δt and then s=t−τ, we finally obtain(59)Vz,t=ev¯z1−e−δt−vze−δt∫0tq1eδ−λ1t−s+q2eδ−λ2t−sev¯z1−e−δsds.
Derivation of (50): through integration by parts, we express the last integral of (59) in an equivalent form:(60)vze−δt∫0tq1eδ−λ1t−s+q2eδ−λ2t−sev¯z1−e−δsds=e−δt∫0tq1eδ−λ1t−s+q2eδ−λ2t−seδsdev¯z1−e−δs=ev¯z1−e−δt−q1e−λ1t−q2e−λ2t−∫0tq1λ1eλ1s−t+q2λ2eλ2s−tev¯z1−e−δsds.
The substitution of this into (59) helps us rewrite the generating function Vz,t as(61)Vz,t=q1e−λ1t+q2e−λ2t+∫0tq1λ1eλ1s−t+q2λ2eλ2s−tev¯z1−e−δsds.
By the conversion formula (19), P0t in (50) is derived from (61) directly as P0t=V−1,t. For m≥1,(62)∂mVz,t∂zm=v¯m∫0tq1λ1eλ1s−t+q2λ2eλ2s−t1−e−δsmev¯z1−e−δsds,which leads to (50) by substituting v¯=v/δ and z=−1 into this expression and dividing it by m!.
Derivation of (51): with the help of (59) and the elementary identity (48), we find(63)∂mVz,t∂zm=v¯m1−e−δtmev¯z1−e−δt−v¯me−δt∫0tq1eδ−λ1t−s+q2eδ−λ2t−s×mδ+vz1−e−δs1−e−δsm−1ev¯z1−e−δsds.
Then, (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 theSTL1 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 inS. cerevisiae range around a median of 11 min  [33], so we fix δ=ln2/11≈0.063min−1. We set the transcription rate as its upper bound v/δm=10  [34], and thus v=0.63min−1. Also, we set a relative small strength rate λ1=0.1min−1 for the basal pathway and a relative large strength rate λ2=3min−1 for the signaling pathway. Such parameter set generates a prolonged region of the gene-off time (denoted by Toff):(64)Toff=q1λ1+q2λ2∈1λ2,1λ1,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 q1≈q2. We set q1=0.51 and q2=0.49 that generates transcriptional efficiency1/Toff=0.19min−1. 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 q2>q1. We then set q1=0.3 and q2=0.7 that generates 1/Toff=0.3min−1. Also, it was reported that the gene-on state could be significantly stabilized in response to increasing stress level, and thus γ≪1  [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 q1=1 so that the gene is activated by a single pathway with strength rate λ=λ1. 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].Figure 1
Cross-talking pathways generate typical dynamical transitions from mRNA decaying distribution to unimodal distribution undergoing the significant intermediate bimodal stage for stress genes in yeast withv=0.63,δ=0.063,λ1=0.1,andλ2=3 (min−1). (a) Dynamical distribution of the stress gene under mild stress level equipped with q1=0.51,q2=0.49, and γ=0.126min−1. (b) Distribution dynamics of the stress gene under high stress level with q1=0.3,q2=0.7, and γ=0min−1. In (a) and (b), the solid circles in lines represent analytical solution, whereas the hollow circles represent the numerical result obtained by Gillespie’s algorithm. The embedded distribution dynamics are generated by the single pathway in which durations of bimodal distribution are much shorter than that generated by cross-talking pathways.
(a)(b)


## 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  [1–3]. 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 Tat 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λ1 and selection probability q2, or alternately by an inducible signaling pathway with strength λ2 and selection probability q2=1−q1. The probability q1 (or q2) 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 λ2 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(65)λ1≪λ2.On the assumption of Theorem1, (65) leads to q1∼q2, 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 106 folds for the Lac promoter in E. coli  [32] and possible 105 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.



---

*Source: 6402703-2020-06-17.xml*