Intrinsic noise in stochastic models of gene expression with molecular memory and bursting

0.43 MB·English

by  Tao Jia

Intrinsic noise in stochastic models of gene expression with molecular memory and
bursting
Tao Jia∗ and Rahul V. Kulkarni†
Department of Physics,
Virginia Polytechnic Institute and State University,
Blacksburg, VA 24061
(Dated: January 12, 2011)
Regulation of intrinsic noise in gene expression is essential for many cellular functions. Corre-
spondingly, there is considerable interest in understanding how different molecular mechanisms of
geneexpressionimpactvariationsinproteinlevelsacrossapopulationofcells. Inthiswork,weana-
lyzeastochasticmodelofburstygeneexpressionwhichconsidersgeneralwaiting-timedistributions
governing arrival and decay of proteins. By mapping the system to models analyzed in queueing
1 theory, we derive analytical expressions for the noise in steady-state protein distributions. The
1 derived results extend previous work by including the effects of arbitrary probability distributions
0 representingtheeffectsofmolecularmemoryandbursting. Theanalyticalexpressionsobtainedpro-
2 videinsightintotheroleoftranscriptional,post-transcriptionalandpost-translationalmechanisms
n in controlling the noise in gene expression.
a
J PACSnumbers: 87.10.Mn,82.39.Rt,02.50.-r,87.17.Aa

1
1
Regulation of gene expression is at the core of cellu-

] lar adaptation and response to changing environments.
N
Given that the underlying processes are intrinsically
M stochastic, cellular regulation must be designed to con-
trolvariability(noise)ingeneexpression[1]. Whilenoise
.
o reduction is essential in many cases, regulatory mecha-
i
b nisms can also exploit the intrinsic stochasticity to in-
- crease noise and generate phenotypic heterogeneity in a FIG. 1. Reaction scheme for the underlying gene expression
q
clonal population of cells [2]. Quantifying the contribu- model. Production of mRNAs occurs in bursts (character-
[
tionsofdifferentsourcesofintrinsicnoiseusingstochastic ized by random variable mb with arbitrary distribution) and
1 modelsofgeneexpression[3–5]isthusanimportantstep each mRNA gives rise to a burst of proteins (characterized
v by random variable p with arbitrary distribution) before it
towards understanding cellular processes and variations b
5 decays (with lifetimeτ ). The waiting-time distributionsfor
in cell populations. m
7 burst arrival and decay of proteins are characterized by the
1 Several recent studies have focused on quantifying functions f(t) and h(t) respectively.
2
noise in gene expression. Experiments have shown that
.
1 protein production often occurs in ‘bursts’ [6, 7] and
0
single-molecule measurements have also provided evi-
1 spectively [12]). The underlying reaction scheme for the
dencefortranscriptionalbursting,i.e.productionofmR-
1 models analyzed in this work is shown in Fig. 1. Pro-
: NAs in bursts [8–10]. The analysis and interpretation of
v duction of mRNAs occurs in independent bursts and the
suchexperimentalstudieshasbeenaidedbythedevelop-
i time interval between the arrival of consecutive mRNA
X ment of coarse-grained stochastic models of gene expres-
bursts is characterized by random variable T with corre-
r sion. The simplest of these considers the basic processes
a sponding probability density function (p.d.f) f(t). The
(transcription, translation and degradation) as elemen-
number of mRNAs produced in a single transcriptional
taryPoissonprocesses[11]withexponentialwaiting-time
burst is characterized by the random variable m . Each
distributions. However, since these processes are known b
mRNA independently gives rise to a random number of
to involve multiple biochemical steps, the corresponding
proteins (characterized by random variable p ) before it
waiting-time distributions can be more general than the b
is degraded. For the basic models of translation, p fol-
‘memoryless’exponentialdistribution[12]. Animportant b
lowsthegeometricdistribution[6,7,13]. However,more
questionthenarises: howdogeneexpressionmechanisms
general schemes of gene expression (e.g. involving post-
involvingmolecularmemoryeffectsinfluencethenoisein
transcriptional regulation [14]) can give rise to protein
protein distributions?
burst distributions that deviate significantly from a geo-
Motivatedbytheprecedingobservations,weintroduce metricdistribution. Proteinsaredegradedindependently
a model including general waiting-time distributions for and the waiting-time distribution for protein decay is
processes governing the arrival of bursts and the decay characterized by the p.d.f h(t).
ofproteins(termed‘gestation’and‘senescence’effectsre- In the limit that the mRNA lifetime (τ ) is much
m
2
shorter than the protein lifetime (τp), i.e. ττmp (cid:28) 1, the systemisGIX/M/∞whereM indicatesthattheprocess
evolution of cellular protein concentrations can be mod- ofcustomerdeparture, whichistheanalogofproteinde-
eled by processes governing arrival and decay of proteins cay, is Markovian. A(z) corresponds to the generating
alone [13, 15]. Unless otherwise stated, the analysis in function of burst size distribution (determined by ran-
this paper will focus on this ‘burst’ limit, in which pro- domvariablesmb andpb inFig. 1)andN(t)denotesthe
teins are considered to arrive in independent instanta- number of proteins in the cell at time t. The previous
neous bursts arising from the underlying mRNA burst. analysis [17] has derived expressions for the steady-state
In this limit, we have shown in recent work [16] that the mean and variance corresponding to N = limt→∞N(t)
processes involved in gene expression can be mapped on for the GIX/M/∞ queue as [18]:
to models analyzed in queueing theory. In this mapping,
1
individualproteinsaretheanalogsofcustomersinqueue- E[N]= µ (cid:104)T(cid:105)A1
ing models. The bursty synthesis of proteins then corre- p
f (µ ) A
sponds to the arrival of customers in ‘batches’, whereas Var[N]=E[N](1+ L p A −E[N]+ 2 ),(1)
the protein decay-time distribution is the analog of the 1−fL(µp) 1 2A1
service-time distribution for each customer. Given that
where (cid:104)T(cid:105) is the mean of p.d.f f(t) and f (s) is the
L
degradation of each protein is independent of others in
Laplace transform of f(t).
the system, the process maps on to queueing systems
To translate the result Eq.(1) into an expression for
with infinite servers. Correspondingly, the gene expres-
the noise in protein distributions, we derive expressions
sion model in Fig. 1 maps on to what is known as a
forA andA intermsofvariablescharacterizingmRNA
GIX/G/∞systeminthequeueingliterature. Inthisno- 1 2
and protein burst distributions. In general, each mRNA
tation, the symbol G refers to the general waiting-time
will produce a random number of proteins (p ) and fur-
distribution and IX indicates that the customers arrive b
thermore the number of mRNAs in the burst is also a
in batches of random size X, where X is drawn indepen-
randomvariable(m ). Thenumberofproteinsproduced
b
dently each time from an arbitrary distribution.
in a single burst is thus a compound random variable.
The GIX/G/∞ system has been analyzed in previ-
Correspondingly,usingstandardresultsfromprobability
ous work in queueing theory [17]. In the following, we
theory [19], we derive the following equations for burst
briefly review the notation and relevant results from the
size parameters (A and A ) in terms of m and p :
1 2 b b
queueing theory analysis. As in Fig. 1, f(t) and h(t)
denote the p.d.f. for the arrival time and service time A =(cid:104)m (cid:105)(cid:104)p (cid:105)
1 b b
respectively, with F(t) and H(t) as the corresponding A =(cid:104)m (cid:105)(σ2 −(cid:104)p (cid:105))+(σ2 +(cid:104)m (cid:105)2)(cid:104)p (cid:105)2, (2)
cumulative density functions (c.d.f). The distribution of 2 b pb b mb b b
batch size X has the corresponding generating function where the symbols (cid:104)..(cid:105) and σ represent the mean and
A(z), defined as A(z) =(cid:80)∞ P(X =i)zi. The kth fac- standard deviation respectively.
i=1
torial moment of batch size X, denoted by A , is given UsingEq.(2),incombinationwithidentificationofthe
k
by A = (dkA(z)/dzk)| . The number of customers randomvariableN withthecorrespondingvariablechar-
k z=1
in service at time t is denoted by N(t) and analytical acterizingthenumberofproteins(ps), weobtainthefol-
expressions have been derived for the rth binomial mo- lowing expressions for the mean and coefficient of vari-
mentB (t)ofN(t)[17]. Theseresultscanbeusedtode- ance (noise) of the steady-state protein distribution:
r
riveexpressionsforallthemomentsofN(t),forexample τ
E[N(t)]=B1(t)andVar[N(t)]=2B2(t)+B1(t)−B12(t). (cid:104)ps(cid:105)= (cid:104)Tp(cid:105)(cid:104)mb(cid:105)(cid:104)pb(cid:105)
In the following, we will focus on two general subcat-
σ2 1 (cid:104)T(cid:105) (cid:16)
egories of the GIX/G/∞ system for which closed-form ps = + × K +σ2 /(cid:104)m (cid:105)2
analytical expressions can be derived for the mean and (cid:104)ps(cid:105)2 (cid:104)ps(cid:105) 2τp g mb b
varianceofsteady-stateproteindistributions. Thesecor- σ2 /(cid:104)p (cid:105)2−1/(cid:104)p (cid:105)(cid:17)
+ pb b b , (3)
respond to two cases: A) arbitrary distributions for ges- (cid:104)m (cid:105)
b
tation and bursting with a Poisson process governing
where
protein degradation and B) arbitrary distributions for
bursting and senescence with a Poisson process govern- (cid:16) f (µ ) 1 (cid:17)
ing burst arrival. Kg =2 1−Lf (pµ ) − µ (cid:104)T(cid:105) +1, (4)
L p p
Consider first case A, for which arbitrary gestation
and bursting effects are included. In this case, the ran- is denoted as the gestation factor.
dom variable T characterizing the time interval between Different contributions to the noise in protein dis-
burstsisdrawnfromanarbitraryp.d.f. f(t). Theprotein tributions are highlighted in Eq.(3): gestation effects,
decay-timedistributionh(t)istakentobeanexponential mRNAtranscriptionalbursting, andtranslationalburst-
function with h(t) = µpe−µpt and the mean protein life- ing from a single mRNA, which correspond to the terms
time is given by τ =1/µ . The corresponding queueing K ,σ2 /(cid:104)m (cid:105)2 andσ2 /(cid:104)p (cid:105)2,respectively. Thefirsttwo
p p g mb b pb b
3
[3,20]. Usingtheapproximationthatthetime-averaging
factor is the same for general gestation and bursting dis-
tributions, we obtain
σ2 1 (cid:104)T(cid:105) (cid:16)
ps ≈ + × K +σ2 /(cid:104)m (cid:105)2
(cid:104)p (cid:105)2 (cid:104)p (cid:105) 2τ g mb b
s s p
+ σp2b/(cid:104)pb(cid:105)2−1/(cid:104)pb(cid:105)(cid:17)× τp , (5)
(cid:104)m (cid:105) τ +τ
b m p
FIG. 2. The noise vs µp(cid:104)T(cid:105) from analytical expressions and It is instructive to compare Eq.(5) with the result
stochasticsimulations. A)Thetimeintervalbetweenconsec- derived in previous work [12] which assumes the ba-
utiveburstsisfixedandonly1mRNAisproducedeachburst. sic protein production reaction scheme such that σ2 =
The protein production is under post-transcriptional regula- pb
(cid:104)p (cid:105)2+(cid:104)p (cid:105). Considering this specific case, we note that
tion [14] such that σ2 =0.67(cid:104)p (cid:105)2+(cid:104)p (cid:105) and τ /τ ≈0.02. b b
pb b b m p Eq.(5) is identical to the previous result [12] apart from
B)ThetimeintervalbetweenburstsisdrawnfromaGamma
distribution and the number of mRNAs created in one burst the terms corresponding to the gestation factor Kg. The
isdrawnfromaPoissondistribution. Thenumberofproteins connectiontothepreviousresultcanbeseenbyexpand-
createdbyeachmRNAfollowsageometricdistribution. The ing the Laplace transform, f (µ ), in terms of moments
L p
parameters are τm/τp = 0.2, (cid:104)mb(cid:105) = 10, σm2b/(cid:104)mb(cid:105)2 = 0.1 of T. By assuming µp(cid:104)T(cid:105) is small and (cid:104)Tn(cid:105) scales as
and σT2/(cid:104)T(cid:105)2 = 0.2. While Eq.(5) agrees with simulations, thenth powerof(cid:104)T(cid:105)orless, Kg canbeapproximatedby
theresultfromRef. [12]islessaccuratewhenµ (cid:104)T(cid:105)islarge.
p K ≈σ2/(cid:104)T(cid:105)2 which corresponds to the previous result.
g T
Since the parameter 1/(µ (cid:104)T(cid:105)) measures the mean num-
p
ber of bursts occurring during the protein lifetime, this
indicatesthatthepreviousresult[12]isvalidforthecase
terms can be modified by transcriptional regulation and offrequentburstingduringaproteinlifetime,andbreaks
the last term can be tuned by post-transcriptional reg- down when bursts occur over larger time intervals (Fig.
ulation. It is noteworthy that each source contributes 2B).
additively to the overall noise in the steady-state distri- We now consider case B, which corresponds to arbi-
bution. Moreover,whilethenoiseduetogestationeffects trary distributions for bursting and senescence effects
is independent of the degree of transcriptional bursting, along with exponential waiting-time distributions for
the noise contribution from translational bursting is ef- burst arrival. For this case, we take the waiting-time for
fectively reduced by transcriptional bursting. protein degradation to be drawn from an arbitrary dis-
While Eq.(3) is valid for general gestation effects, it is tribution characterized by p.d.f h(t) and c.d.f H(t). The
ofinteresttoconsiderspecificexamples. Weconsiderthe waiting-time between consecutive bursts is characterized
casesuchthatthereisaconstantdelaybetweenarrivalof by an exponential distribution with f(t) = λe−λt. The
consecutivemRNAbursts,i.e.thewaiting-timedistribu- correspondingsystem,followingthemappingtoqueueing
tion is f(t)=δ(t−T ). In this case, the gestation factor theory, is the MX/G/∞ queue. The steady-state mean
d
isgivenbyKg =2e−µpTd/(1−e−µpTd)−2/µpTd+1. The and variance of N for this queue has been obtained in
corresponding expression for the noise in protein distri- previous work [17]:
butions Eq.(3), considering a general case which also in-
cludes the effects of post-transcriptional regulation [14], (cid:90) ∞
E[N]=λA [1−H(t)]dt
isinexcellentagreementwithresultsfromstochasticsim- 1
0
ulations (Fig. 2A). It is noteworthy that K can be non- (cid:90) ∞
g
vanishing even though the time interval between consec- Var[N]=E[N]+λA2 [1−H(t)]2dt. (6)
utiveburstsisfixed(i.e.σ2 =0). Incontrasttoprevious 0
T
work [12], which suggests that the contribution of gesta-
By taking Eq.(2) and the relation (cid:104)T(cid:105) = 1/λ into ac-
tioneffectstothenoisevanisheswhenσ2 =0,ourresult
T count, the mean and the noise for arbitrary senescence
showsthatK canbetunedfrom0to1asµ T isvaried.
g p d and bursting distribution can be derived as:
While the results derived above are valid in the limit
τcmase(cid:28)(i.τep.,waintheoxuatctinevxopkriensgsitohnefcoorntdhietinoonisτemin(cid:28)thτepgaenndefroarl (cid:104)ps(cid:105)= (cid:104)AT1(cid:105)(cid:90) ∞[1−H(t)]dt= (cid:104)τTp(cid:105)(cid:104)mb(cid:105)(cid:104)pb(cid:105)
0
general gestation and bursting distributions) is difficult
σ2 1 (cid:104)T(cid:105) (cid:16)
to obtain. However, a useful approximation can be ob- ps = + × 1+σ2 /(cid:104)m (cid:105)2
tainedbynotingthat,forthebasicgeneexpressionmod- (cid:104)ps(cid:105)2 (cid:104)ps(cid:105) 2τp mb b
els, the exact result is obtained by scaling the terms in σ2 /(cid:104)p (cid:105)2−1/(cid:104)p (cid:105)(cid:17)
thebracketinEq.(3)withatime-averagingfactor τmτ+pτp + pb b(cid:104)mb(cid:105) b ×Ks, (7)
4
where tween models of stochastic gene expression and queue-
ingsystemswhichhaspotentialapplicationsforresearch
2(cid:82)∞[1−H(t)]2dt 2(cid:82)∞H(t)[1−H(t)]dt
K = 0 =2− 0 , in both fields. The extensive analytical approaches and
s τp τp toolsdevelopedinqueueingtheorycannowbeemployed
(8) to analyze stochastic processes in gene expression. It is
is denoted as the senescence factor. also anticipated that future analysis of regulatory mech-
ItisnoteworthyEq.(7)andEq.(3)havemultipleterms anismsforgeneexpressionwillleadtonewproblemsand
in common. The terms characterizing the noise from challenges for queueing theory.
transcriptional and translational bursting remain un-
changed. However, unlike the gestation factor that con-
The authors acknowledge funding support from NSF
tributes to the total noise additively, the senescence fac-
(PHY-0957430) and from ICTAS, Virginia Tech.
tor serves as a scaling factor for the total noise. While
there is no obvious upper limit on the value of K , the
g
upperboundforK is2asisevidentfromEq.(8). Ingen-
s
eral, as the distribution h(t) grows more sharply peaked,
theK valueincreases. Whenh(t)becomesadeltafunc-
s ∗ [email protected]
tion, Ks reaches its maximum value. † [email protected]
The general results derived in this work will serve as [1] M. Kaern, T. C. Elston, W. J. Blake, and J. J. Collins,
usefulinputsfortheanalysisandinterpretationofdiverse Nat Rev Genet 6, 451 (2005).
experimental studies of gene expression. Some examples [2] A. Raj and A. van Oudenaarden, Cell 135, 216 (2008).
are: 1)Recentexperimentsonsingle-cellstudiesofHIV-1 [3] J. M. Paulsson, Phys Of Life Rev 2, 157 (2005).
[4] S. Azaele, J. R. Banavar, and A. Maritan, Phys. Rev. E
viralinfectionshavefocusedonthefrequencyanddegree
80 (2009).
oftranscriptionalbursting[21]. Forsuchstudies, thede-
[5] B.Munsky,B.Trinh,andM.Khammash,Mol.Sys.Biol.
rivedresultscanbeusedtorelatemeasurementsofinter-
5 (2009).
arrivalwaiting-timedistributionsandburstdistributions [6] L. Cai, N. Friedman, and X. S. Xie, Nature 440, 358
to the noise in protein distributions. 2) Experimental (2006).
data and computational models of the cell-cycle in yeast [7] J. Yu, J. Xiao, X. Ren, K. Lao, and X. S. Xie, Science
indicate that modeling the basic processes of gene ex- 311, 1600 (2006).
[8] I. Golding, J. Paulsson, S. M. Zawilski, and E. C. Cox,
pression as Poisson processes gives rise to unrealistically
Cell 123, 1025 (2005).
large noise in protein distributions [22], thereby suggest-
[9] A. Raj, C. S. Peskin, D. Tranchina, D. Y. Vargas, and
ing that regulatory schemes which change distributions
S. Tyagi, PLoS Biol 4, e309 (2006).
toreducethenoiseareemployedbythecell. Theanalyt- [10] J.Chubb,T.Trcek,S.Shenoy,andR.Singer,Curr.Biol.
ical expressions derived highlight different contributions 16, 1018 (2006).
to noise and can thus provide insight into how different [11] M. Thattai and A. van Oudenaarden, Proc Natl Acad
regulatory schemes can lead to noise reduction. 3) More Sci U S A 98, 8614 (2001).
[12] J.M.PedrazaandJ.Paulsson,Science319,339(2008).
generally, the results derived can be used in the analysis
[13] N. Friedman, L. Cai, and X. S. Xie, Phys Rev Lett 97,
of inverse problems, i.e. using experimental measure-
168302 (2006).
ments of intrinsic noise to determine parameters of the
[14] T. Jia and R. Kulkarni, Phys. Rev. Lett. 105, 018101
underlyingkineticmodels. Suchefforts,inturn,canlead (2010).
to further insights into cellular factors that impact gene [15] V.ShahrezaeiandP.S.Swain,ProcNatlAcadSciUSA
regulation, based on experimnetal observations of noise 105, 17256 (2008).
in gene expression. [16] V. Elgart, T. Jia, and R. V. Kulkarni, Phys. Rev. E 82,
021901 (2010).
In summary, we have analyzed the noise in protein
[17] L.Liu,B.R.K.Kashyap,andJ.G.C.Templeton,Jour.
distributions for general stochastic models of gene ex-
Appl. Prob. 27, 671 (1990).
pression. The present work extends previous analysis by
[18] The result given in Ref. [17] has a minor error which is
deriving analytical results for the noise in protein distri- corrected here.
butions for arbitrary gestation, senescence and bursting [19] S. M. Ross, Introduction to Probability Models, Ninth
mechanisms. The expressions obtained provide insight Edition (Academic Press, Inc., 2006).
into how different sources contribute to the noise in pro- [20] A. Bar-Even, J. Paulsson, N. Maheshri, M. Carmi,
E. O’Shea, Y. Pilpel, and N. Barkai, Nat Genet 38, 636
tein levels which can lead to phenotypic heterogeneity in
(2006).
isogenic populations. The results derived will thus serve
[21] R. Skupsky, J. C. Burnett, J. E. Foley, D. V. Schaffer,
as useful inputs for the analysis and interpretation of
andA.P.Arkin,PLoSComputBiol6,e1000952(2010).
experiments probing stochastic gene expression and its [22] S. Kar, W. T. Baumann, M. R. Paul, and J. J. Tyson,
phenotypic consequences. At a broader level, this work Proceedings of the National Academy of Sciences 106,
demonstrates the benefits of developing a mapping be- 6471 (2009).