Effect of gene-expression bursts on stochastic timing of cellular events
EEffect of gene-expression bursts on stochastic timing of cellular events
Khem Raj Ghusinga
Department of Electrical and Computer EngineeringUniversity of Delaware, Newark, DE, USA { [email protected] } Abhyudai Singh ∗ Department of Electrical and Computer Engineering,Department of Mathematical Sciences,Department of Biomedical Engineering,University of Delaware, Newark, DE, USA { [email protected] } Abstract
Gene expression is inherently a noisy process which manifests as cell-to-cell variability in time evo-lution of proteins. Consequently, events that trigger at critical threshold levels of regulatory proteinsexhibit stochasticity in their timing. An important contributor to the noise in gene expression is transla-tion bursts which correspond to randomness in number of proteins produced in a single mRNA lifetime.Modeling timing of an event as a first-passage time (FPT) problem, we explore the effect of burst sizedistribution on event timing. Towards this end, the probability density function of FPT is computedfor a gene expression model with burst size drawn from a generic non-negative distribution. Analyticalformulas for FPT moments are provided in terms of known vectors and inverse of a matrix. The effect ofburst size distribution is investigated by looking at how the feedback regulation strategy that minimizesnoise in timing around a given time deviates from the case when burst is deterministic. Interestingly,results show that the feedback strategy for deterministic burst case is quite robust to change in burstsize distribution, and deviations from it are confined to about 20% of the optimal value. These findingsfacilitate an improved understanding of noise regulation in event timing.
Gene expression, the process by which a gene is transcribed to mRNAs and each mRNA is subsequentlytranslated in to proteins, plays a central role in determining cellular behavior. As the biochemical reactionsare innately probabilistic, and species such as gene, mRNA, etc. often occur in low copy numbers at single ∗ Corresponding author a r X i v : . [ q - b i o . S C ] S e p ell level, expression of a gene is a stochastic process [1–5]. Consequently, even if an isogenic population isinduced at the same time, two cells might have different evolution of protein level over time.An important point in triggering of several cellular events is attainment of a threshold level of theprotein being expressed. For example, an environmental cue or internal signal usually induces expressionof a regulatory protein which subsequently activates an appropriate cellular response. The activationtakes place when a certain threshold level of the regulatory protein is achieved [6, 7]. Other examples ofsuch events include cell-fate decisions [8–17], temporal program of gene activation, etc. [18, 19]. Becauseexpression of a gene is a stochastic process, a target protein level might be reached at different times inindividual cells of an isogenic clonal population induced at the same time.The timing of such threshold crossing events can be mathematically formulated using a first-passagetime (FPT) process. First-passage time is defined as the first time at which a random walker reacheda certain critical level, and has been used in several fields to study threshold crossing phenomena [20–32]. Particularly, in the context of stochastic gene expression models, previous works have developedexact analytical formulas for FPT moments when gene expresses in translation bursts [32–35]. In thesemodels, one usually assumes a geometrically distributed translation burst distribution which arises fromthe assumption that both mRNA translation and mRNA degradation are one step processes occurring atexponentially distributed times [36].Here, we relax these assumptions and extend the FPT calculations to include any arbitrary burst sizedistribution. In addition to modeling the mRNA translation and degradation as multistep processes [37],the arbitrary burst size distribution allows one to potentially model post transcriptional regulation aswell [38]. We consider a minimal gene expression model that consists of protein arrival in bursts, andits degradation. The arrival rates of bursts is assumed to be dependent on the protein level therebyimplementing a feedback transcriptional control. A first-passage problem for this model is formulated, andthe probability density function of FPT is determined exactly. Furthermore, an exact analytical formulafor the FPT moments is also provided in terms of product of known vectors and inverse of a matrix.These formulas are written as simple series summations for a simple case when the protein of interestis assumed to be stable. Considering different possible distributions of the burst size, we investigate theoptimal feedback strategy that might minimize noise in timing of an event around a fixed time. We showthat the effect of bursting can deviate the optimal feedback strategy from a no feedback that results from abirth-death model of gene expression. However, these deviations are within 20% of the optimal no feedbackstrategy.Remainder of the paper is organized as follows. In section II, we formulate a gene expression model thatdescribes production of a protein in bursts and its degradation. In the next section, FPT computationsfor this gene expression model are performed. Section IV deals with determining the moments of FPT,particularly the expressions of first two moments for various distributions of burst size. The effect ofdifferent burst size distributions on deviations from the optimal feedback strategy are examined in sectionV. Finally, the results and potential directions of research are discussed in section VI.2 Model Description
We consider expression of a protein from its gene that is induced at t = 0 as shown in Fig. 1. The modelconsists of four fundamental components of the gene expression process, namely, transcription (productionof mRNAs from gene), translation (production of proteins from a mRNA), mRNA degradation, and proteindegradation. We assume the transcription rate to be an arbitrary function of the protein level whichcorresponds to feedback regulation. In terms of notations, we denote the transcription rate when proteinlevel x ( t ) = i by k i , and the protein degradation rate by γ . Typically the mRNA half-life is much smallerthan the protein half-life [36, 39–41], and this time-scale separation can be exploited to ignore the mRNAdynamics. As a result the model reduces to a bursty birth-death process in which each transcription eventcreates a mRNA molecule and it degrades immediately after synthesizing a burst of protein molecules. Gene ∅ ∅ Protein mRNA 𝑘 𝑖 γ Figure 1: Model schematic showing a gene transcribing to mRNAs which are translated into proteins. Afeedback regulation is implemented by assuming the transcription rate to be function of the protein level(the transcription rate is k i when protein level is i ) . The mRNA translation and degradation processesare assumed to be general, giving rise to some arbitrary non-negative distribution of burst size. Finally,each protein molecule is assumed to degrade with a rate γ .In essence, the probabilities of occurrence of an arrival event and a degradation event in an infinitesimaltime interval ( t, t + dt ) are given by P { x ( t + dt ) = i + B | x ( t ) = i } = k i dt, (1a) P { x ( t + dt ) = i − | x ( t ) = i } = iγdt, (1b)where B represents the burst size. We consider that B follows an arbitrary non-negative distribution as P ( B = i ) = β i , i ∈ { , , , . . . , ∞} , (2a) ∞ (cid:88) i =0 β i = 1 . (2b)The arbitrary distribution of B allows us to relax assumptions on the translation, and degradation stepof mRNA being one step processes [37] and also incorporate any post transcriptional regulation [38]. For3xample, consider one-step degradation of mRNA, i.e., exponentially distributed mRNA lifetime. Denotingthe average mRNA lifetime by 1 / ˜ γ and the translation rate of a protein from a mRNA by ˜ k , the burst sizedistribution can be computed as P ( B = i ) = (cid:90) ∞ ˜ γe − ˜ γs (˜ ks ) i i ! e − ˜ ks ds (3a)= (cid:32) ˜ k ˜ k + ˜ γ (cid:33) i (cid:18) ˜ γ ˜ k + ˜ γ (cid:19) (3b)which is a geometric distribution [36]. If instead of an one-step mRNA degradation, a multistep degradationis considered that corresponds to an Erlang distributed mRNA lifetime, then using a similar integral as(3) results in a negative binomial distribution. In the extreme case when the mRNA lifetime is consideredto be deterministic, the integral results in a Poisson distributed burst size.Another advantage of considering a general burst size distribution is that it also accounts for staticextrinsic noise. As an example, we can consider that the burst size follows a geometric distribution givenby (3), and the average burst size is affected by some enzyme Z that is drawn from a distribution, i.e.,its levels do not fluctuate over the time-scale of the event. The factor Z here might represent cell-to-cellvariability in some factor that affects the translation machinery. The resulting burst size distribution canbe written as P ( B = i ) = (cid:88) z ∈ support of Z (cid:32) ˜ kz ˜ kz + ˜ γ (cid:33) i (cid:18) ˜ γ ˜ kz + ˜ γ (cid:19) P ( Z = z ) , (4)which can be well described by (2). In the next section, we compute the first-passage time distribution ofthe bursty birth-death process. The first-passage time (FPT) is defined as the first-time at which a stochastic process x ( t ) crosses athreshold X . Mathematically, we are interested in determining the probability distribution function (pdf)of the following random variable T = inf { t ≥ x ( t ) ≥ X } . (5)As done in our previous work [34], the first-passage time can be computed by constructing an equivalentbursty birth-death process wherein all states greater than or equal to X are absorbing. The protein countevolves as per probabilities of occurrences given in (1) until one of the absorbing state is achieved. Thedifference between this equivalent formulation and the original bursty birth-death process is that in theoriginal formulation, the protein count can return back to the states less than X . However, in termsof first-passage times both processes are equivalent. The first-passage time probability density can be4omputed as f T ( t ) = X − (cid:88) i =0 k i P ( B ≥ X − i ) p i ( t ) , (6)where p i ( t ) = P ( x ( t ) = i ). Intuitively, this formula can be interpreted as follows: the process crosses thethreshold X for the first time at time t + dt if the protein count was equal to i at time t and a burst ofsize greater than or equal to X − i occurred in the next infinitesimal time interval ( t, t + dt ).For sake of convenience, we can express (6) as f τ ( t ) = U (cid:62) P ( t ) , (7a)where P ( t ) = (cid:104) p ( t ) p ( t ) p ( t ) . . . p X − ( t ) (cid:105) (cid:62) , (7b)and U is given as U = (cid:104) k (cid:16) − (cid:80) X − j =0 β j (cid:17) . . . k X − (1 − β ) (cid:105) , (7c)where we have used P ( B ≥ i ) = 1 − i − (cid:88) j =0 β j (7d)which is obtained from (2).The time evolution of P ( t ) can be obtained from the forward Kolmogorov equation (also called thechemical master equation) for the equivalent bursty birth-death process. This can be written as˙ p ( t ) = − k (1 − β ) p ( t ) + γp ( t ) , (8a)˙ p i ( t ) = − ( k i (1 − β ) + iγ ) p i ( t ) + ( i + 1) γp i +1 ( t )+ i − (cid:88) n =0 k n β i − n p n ( t ) , ≤ i ≤ X − , (8b)˙ p X − ( t ) = − ( k X − (1 − β ) + ( X − γ ) p X − ( t )+ X − (cid:88) n =0 k n β X − n − p n ( t ) . (8c)These equations can be compactly written in form of a linear matrix differential equation˙ P = AP , (9a)5here a i,j , the element in i th row and j th column of the Hessenberg matrix A , is given as a i,j = , j > i + 1 , ( i − γ, j = i + 1 , − k i − (1 − β ) − ( i − γ, j = i,k j − β i − , j < i. (9b)The solution (9a) is given by the following P ( t ) = exp( A t ) P (0) , (10)where P (0) is the vector consisting of probabilities of the initial protein count. In this work, we consider x (0) = 0 which leads to P (0) = (cid:104) . . . (cid:105) (cid:62) , (11)although any other value of x (0) can be taken as long as it is less than X . One can also choose a probabilitydistribution for x (0).Using (10) in (7a) yields the following for the first-passage time probability density function f T ( t ) = U (cid:62) exp( A t ) P (0) . (12)One can numerically compute the pdf in (12) for given model parameters. In Figure 2, we plot it forseveral different distributions of the burst size. The production and degradation rates, the event threshold X , and the average burst size are taken to be constant across these distributions. It can be seen that thedeterministic burst size distribution has a tighter distribution than others – a feature that is expected.Interestingly, the qualitative shape of the FPT distribution does not appear to vary between different burstsize distributions. Next, we discuss how the moments of FPT can be computed using the pdf in (12). One is often interested in estimating a few low order moments, particularly the average time to an event andits second order moment. Using the probability density function in (12), the moments of the first-passagetime can be calculated as (cid:104) T m (cid:105) = (cid:90) ∞ t m U (cid:62) exp( A t ) P (0) dt (13a)= U (cid:62) (cid:18)(cid:90) ∞ t m exp( A t ) (cid:19) P (0) dt. (13b)The above integral can be explicitly written in terms of inverse of A (cid:104) T m (cid:105) = m !( − m +1 U (cid:62) (cid:0) A − (cid:1) m +1 P (0) . (14)6 eterministic GeometricPoissonNegativebinomial Bimodal Time (t) P r obab ili t y den s i t y f T ( t ) Figure 2: First-passage time distributions for different translation burst distributions show qualitativelysimilar shapes. The mean burst size of 1 molecule, transcription rate k i = 2 min − , degradation rate γ =0 .
01 min − and the event threshold X = 20 molecules are taken to be same across different distributions.Whereas the Geometric, Poisson, and the Deterministic distributions are uniquely determined by the meanburst size, the Negative binomial and Bimodal distributions have additional parameters. For the Negativebinomial distribution the shape parameter is taken to be 5. For the Bimodal distribution, the probabilityof having bursts of size 1 and 2 are respectively taken to be 0 .
50 and 0 .
25 whereas probability of havingburst of size 0 is taken to be 0 . A is Hurwitz stable. As shown in appendix 6, the matrix A in (9b) is indeed Hurwitzstable regardless of the distribution of the burst size.To be able to find a certain moment of FPT, one is required to compute A − . Recall that A is alower Hessenberg matrix, i.e., all elements above its first super-diagonal are zero. There are recursionformulas available to invert Hessenberg matrices which can be used to compute (14), e.g., see [42, 43]. Inour experience, we have seen that computation of A − is well-behaved if one writes A = A + A γ , (15)and computes A − as A − = (cid:0) I + A − A γ (cid:1) − A − . (16)Here A is the lower triangular matrix which is equal to A when γ = 0 and A γ contains the γ terms. Insome special cases such as when the burst size follows a geometric distribution, an exact expression of A − can be obtained and the first two FPT moments can be further simplified to series summations (see eqs.(16)-(17) in [34]). Generally speaking, the resulting formulas obtained from the above expression are quiteconvoluted, and it is hard to get any physical insights from them.It turns out that for a simpler case when γ = 0 (consequently A γ = 0), the resulting formulas can beexpressed in rather elegant forms as presented below. Theorem 1
Consider the FPT probability density function in (12) . When the protein of interest does not egrade, i.e., γ = 0 , the first two moments of FPT are given by (cid:104) T (cid:105) = − X − (cid:88) i =0 α i k i (17a) (cid:10) T (cid:11) = 2 X − (cid:88) i =0 α i k i τ i , τ i = X − (cid:88) j = i α j k j , (17b) where the coefficients α i depend upon the burst size distribution in (3) and are recursively computed as α = − − β , (17c) α i = α β i + α β i − + . . . + α i − β − β , i = 1 , , . . . , X − . (17d) Proof:
When γ = 0, the matrix A is equal to A and given by A = − k (1 − β ) 0 · · · k β − k (1 − β ) · · · k β X − k β X − · · · − k X − (1 − β ) . (18)As A is a lower triangular matrix, its inverse is also lower triangular. It can be shown that A − hasfollowing form A − = α k · · · α k α k · · · α X − k X − α X − k X − · · · α k X − α k X − , (19)where the coefficients α , α , . . . α X − are determined by the burst size distribution as given in (17c)-(17d)(also see Remark 1).Recall from (14) that the mean FPT is given by (cid:104) T (cid:105) = U (cid:62) A − A − P (0) . (20)Since P (0) = (cid:104) . . . (cid:105) (cid:62) , A − P (0) is just the first column of A − . Furthermore, as shown in appendix6 we have that U (cid:62) A − = − (cid:104) . . . (cid:105) . Therefore, (cid:104) T (cid:105) is equal to the negative sum of first columnelements of A − , resulting in (cid:104) T (cid:105) = − X − (cid:88) i =0 α i k i . (21)8n the same manner, we can compute the second order moment as (cid:10) T (cid:11) = − U (cid:62) A − A − A − P (0) (22)= 2 (cid:62) α k · · · α k α k · · · α X − k X − α X − k X − · · · α k X − α k X − α k α k ... α X − k X − . (23)Defining τ i = (cid:80) X − j = i α j k j , the above expression simplifies to the following (cid:10) T (cid:11) = 2 X (cid:88) i =1 α i k i τ i . (24) Remark:
The coefficients α , α , . . . greatly simplify for some distributions of the burst size B . Forexample, if the burst size is assumed to be drawn from a deterministic distribution such that β i = 1when i = 1 and zero otherwise, then the coefficients α i = − i = 0 , , . . . , X −
1. On the otherhand, if the burst distribution is assumed to follow a geometric distribution with parameter µ such that β i = (1 − µ ) i µ , then one obtains α = − − µ , α i = − µ − µ for i ≥
1. Similarly, for a Poisson distributedburst with parameter λ , i.e., β i = λ i i ! e − λ , the resulting coefficients are α = 1 , α i = e λ e λ − λ i i ! ( e λ − i i − (cid:88) m =0 a [ i, m ] e mλ , with a [ r, s ] represents an Eulerian number a [ r, s ] = s +1 (cid:88) i =0 ( − i (cid:18) r + 1 i (cid:19) ( s + 1 − i ) r . So far we have determined the expressions for FPT moments for a general distribution of the burstsize. Using these, we investigate whether for a given burst size distribution there is a specific feedbackregulation mechanism that can schedule cellular events with precision.
As event timing exhibits cell-to-cell variability arising from stochastic nature of gene expression, a problemof interest is to investigate the optimal feedback mechanism that can attenuate variability. Such mecha-nisms could be employed by cells in cases where precision in timing is important. As an simple example, ifone ignores the protein degradation and assumes the burst size to be deterministic, then the FPT moments9re given by (cid:104) T (cid:105) = X − (cid:88) i =0 k i , (cid:10) T (cid:11) − (cid:104) T (cid:105) = X − (cid:88) i =0 k i . (25)For this model, if one minimizes the variance in timing around a fixed mean (cid:104) T (cid:105) = c , the optimal tran-scription rates are given by k i = Xc , i = { , , . . . , X − } . (26)Importantly, these transcription rates are equal to each other which represents a no feedback regulation [44].Furthermore, considering geometrically distributed burst size, and solving for optimal transcription ratesyields k = 2(2 − µ ) + ( X − µ (2 − µ )(1 − µ ) c , (27a) k i = 2(2 − µ ) + ( X − µ (1 − µ ) c , i = { , , . . . , X − } , (27b)where µ is the parameter of geometric distribution [35]. Interestingly, here except for the first transcriptionrate (when the protein level is zero), the other transcription rates are equal, and this optimal feedbackstrategy is quite similar to a no feedback mechanism.Motivated from these findings, we ask how the optimal feedback strategy deviates from a no feedbackcase for other burst size distributions. To this end, we consider various burst size distributions and numer-ically find the corresponding optimal feedback strategies. Our results show that even though the optimalfeedback strategy shows deviations from a no feedback strategy, these deviations are, however, within 20%of the transcription rate for no feedback strategy (Fig. 3). Similar to the geometric distributed burstsize, the first transcription rate is seen to be significantly different than others even for other distributions.We have also checked the optimal feedback strategies for other values of mean burst size and they arequalitatively similar to the ones shown in Fig. 3. Important cellular events are typically governed by accumulation of a regulatory protein up to a criticalthreshold [6, 7, 9–19]. As the expression of the protein is a stochastic process, the resulting timing ofevents is stochastic as well. Such events can be studied as first-passage time problems. To this end, weconsidered a stochastic gene expression model that includes translational bursting in protein production andits degradation. One of the typical key modeling assumptions takes the translation and mRNA degradationevents as one step process. The model considered here relaxes this key assumption and considers anarbitrary non-negative discrete distribution of the burst size. We carried out the FPT calculations for thismodel and found tractable forms of the FPT probability density, and moments. For a special case whenthe protein of interest does not decay, we found elegant series summations that describe the FPT moments.Our results show that even though the FPT distribution is affected by the underlying burst size dis-10 eometricNegativebinomialDeterministicPoissonBimodal
Protein level T r an sc r i p t i on r a t e Figure 3: Optimal transcription rates for different burst size distributions. The transcription rates arecomputed by solving an optimization problem that constraints the mean FPT to be 10 minutes and theevent threshold to be 10 molecules. The mean burst size for each distribution is assumed to be 1 molecule.While the three distributions (Deterministic, Poisson, and Geometric) are uniquely determined by themean burst size, the other three require additional parameters. The shape parameter for the Negativebinomial distribution is taken to be 5. For the Bimodal distribution, the probability of having bursts ofsize 1 and 2 are respectively taken to be 0 .
50 and 0 .
25 whereas probability of having burst of size 0 is takento be 0 .
25. For the Zipf’s law, the exponent is taken to be 10, and the number of elements is taken to be30. 11ribution, its shape does not change much with the burst distribution (Fig. 2). Furthermore, for a simplecase when protein does not degrade, introduction of burst does quantitatively change the optimal feedbackstrategy that would minimize the noise in timing around a fixed given time. However, qualitatively, thefeedback strategy remains close to a no feedback strategy which is optimal for the case when the burstis deterministic (Fig. 3). These observations suggest the robustness of optimal feedback regulation withrespect to burst size distribution.There are several possible directions of future work. For example, one step could be to take proteindegradation into account and perform a systematic analysis of different feedback strategies for various burstsize distributions. Furthermore, the production of mRNAs (transcription) and degradation of proteinsmight also be considered general processes as translation and mRNA degradation.
Appendix
Proof that A is a Hurwitz matrix
An important requirement for the formula in (14) is that inverse of the matrix A should exist, and thatmatrix A itself must be a Hurwitz matrix. To show that the matrix is Hurwitz, we prove that it fulfillsthe following two requirements [45]:1. The diagonal elements a ii < i = 1 , , · · · , X ,2. max ≤ j ≤ X X (cid:88) i =1 j (cid:54) = i (cid:12)(cid:12)(cid:12)(cid:12) a ij a jj (cid:12)(cid:12)(cid:12)(cid:12) < A are a ii = − k i − (1 − β ) − ( i − γ i − <
0, the first requirement above issatisfied. To check the second requirement for each of j = 1 , , · · · , X , note that X (cid:88) i =1 j (cid:54) = i (cid:12)(cid:12)(cid:12)(cid:12) a ij a jj (cid:12)(cid:12)(cid:12)(cid:12) = ( j − γk j − (1 − β ) + ( j − γ + k j − (cid:80) Xi = j +1 β i − j k j − (1 − β ) + ( j − γ ≤ k j − (1 − β ) + ( j − γk j − (1 − β ) + ( j − γ = 1 . (28)Thus A is a Hurwitz matrix. Expression of U (cid:62) A − Using (7c) and (19), the expression of U (cid:62) A − is U (cid:62) A − = k P ( B ≥ X )... k X − P ( B ≥ (cid:62) α k · · · α k α k · · · α X − k X − α X − k X − · · · α k X − α k X − , (29) P ( B ≥ i ) = 1 − (cid:80) i − j =0 β j for i = 1 , , . . . , X . The above expression reduces to U (cid:62) A − = (cid:80) X − i =0 α i (cid:16) − (cid:80) X − i − j =0 β j (cid:17) ... α (1 − β − β ) + α (1 − β ) α (1 − β ) (cid:62) . (30)From the relationship between the coefficients α i and β i in (17c)-(17d), it can be established that each ofthe elements of the above vector is equal to − References [1] W. J. Blake, M. Kærn, C. R. Cantor, and J. J. Collins, “Noise in eukaryotic gene expression,”
Nature ,vol. 422, pp. 633–637, 2003.[2] J. M. Raser and E. K. O’Shea, “Noise in gene expression: origins, consequences, and control,”
Science ,vol. 309, pp. 2010–2013, 2005.[3] A. Raj and A. van Oudenaarden, “Nature, nurture, or chance: stochastic gene expression and itsconsequences,”
Cell , vol. 135, pp. 216–226, 2008.[4] M. Kærn, T. C. Elston, W. J. Blake, and J. J. Collins, “Stochasticity in gene expression: from theoriesto phenotypes,”
Nature Reviews Genetics , vol. 6, pp. 451–464, 2005.[5] A. Singh and M. Soltani, “Quantifying intrinsic and extrinsic variability in stochastic gene expressionmodels,”
PloS One , vol. 8, p. e84301, 2013.[6] H. H. McAdams and A. Arkin, “Stochastic mechanisms in gene expression,”
Proceedings of the Na-tional Academy of Sciences , vol. 94, pp. 814–819, 1997.[7] A. Amir, O. Kobiler, A. Rokney, A. B. Oppenheim, and J. Stavans, “Noise in timing and precision ofgene activities in a genetic cascade,”
Molecular Systems Biology , vol. 3, 2007.[8] J. J. Dennehy and N. Wang, “Factors influencing lysis time stochasticity in bacteriophage λ ,” BMCMicrobiology , vol. 11, p. 174, 2011.[9] K. C. Chen, L. Calzone, A. Csikasz-Nagy, F. R. Cross, B. Novak, and J. J. Tyson, “Integrative analysisof cell cycle control in budding yeast,”
Molecular Biology of the Cell , vol. 15, pp. 3841–3862, 2004.[10] J. M. Bean, E. D. Siggia, and F. R. Cross, “Coherence and timing of cell cycle start examined atsingle-cell resolution,”
Molecular Cell , vol. 21, pp. 3–14, 2006.[11] D. L. Satinover, D. L. Brautigan, and P. T. Stukenberg, “Aurora-A kinase and inhibitor-2 regulatethe cyclin threshold for mitotic entry in xenopus early embryonic cell cycles,”
Cell Cycle , vol. 5,pp. 2268–2274, 2006. 1312] X. Liu, X. Wang, X. Yang, S. Liu, L. Jiang, Y. Qu, L. Hu, Q. Ouyang, and C. Tang, “Reliable cellcycle commitment in budding yeast is ensured by signal integration,” eLife , vol. 4, p. e03977, 2015.[13] S. L. Spencer, S. Gaudet, J. G. Albeck, J. M. Burke, and P. K. Sorger, “Non-genetic origins ofcell-to-cell variability in TRAIL-induced apoptosis,”
Nature , vol. 459, pp. 428–432, 2009.[14] J. Roux, M. Hafner, S. Bandara, J. J. Sims, H. Hudson, D. Chai, and P. K. Sorger, “Fractional killingarises from cell-to-cell variability in overcoming a caspase activity threshold,”
Molecular SystemsBiology , vol. 11, p. 803, 2015.[15] M. Kracikova, G. Akiri, A. George, R. Sachidanandam, and S. Aaronson, “A threshold mechanismmediates p53 cell fate decision between growth arrest and apoptosis,”
Cell Death & Differentiation ,vol. 20, pp. 576–588, 2013.[16] K. Carniol, P. Eichenberger, and R. Losick, “A threshold mechanism governing activation of thedevelopmental regulatory protein σ f in bacillus subtilis,” Journal of Biological Chemistry , vol. 279,pp. 14860–14870, 2004.[17] P. J. Piggot and D. W. Hilbert, “Sporulation of bacillus subtilis,”
Current Opinion in Microbiology ,vol. 7, pp. 579–586, 2004.[18] I. Nachman, A. Regev, and S. Ramanathan, “Dissecting timing variability in yeast meiosis,”
Cell ,vol. 131, pp. 544–556, 2007.[19] J. M. Pedraza and J. Paulsson, “Random timing in signaling cascades,”
Molecular Systems Biology ,vol. 3, 2007.[20] D. Middleton, A. Veitch, and R. Nisbet, “The effect of an upper limit to population size on persistencetime,”
Theoretical Population Biology , vol. 48, pp. 277–305, 1995.[21] J. Grasman and R. HilleRisLambers, “On local extinction in a metapopulation,”
Ecological Modelling ,vol. 103, pp. 71–80, 1997.[22] P. Fauchald and T. Tveraa, “Using first-passage time in the analysis of area-restricted search andhabitat selection,”
Ecology , vol. 84, pp. 282–288, 2003.[23] O. Ovaskainen and B. Meerson, “Stochastic models of population extinction,”
Trends in Ecology &Evolution , vol. 25, pp. 643–652, 2010.[24] G. H. Weiss and M. Dishon, “On the asymptotic behavior of the stochastic and deterministic modelsof an epidemic,”
Mathematical Biosciences , vol. 11, pp. 261–265, 1971.[25] A. L. Lloyd and R. M. May, “How viruses spread among computers and people,”
Science , vol. 292,pp. 1316–1317, 2001.[26] D. Volovik and S. Redner, “First-passage properties of bursty random walks,”
Journal of StatisticalMechanics: Theory and Experiment , vol. 2010, p. P06018, 2010.1427] D. A. Charlebois, N. Abdennur, and M. Kaern, “Gene expression noise facilitates adaptation and drugresistance independently of mutation,”
Physical Review Letters , vol. 107, p. 218101, 2011.[28] T. Chou and Y. Wang, “Fixation times in differentiation and evolution in the presence of bottlenecks,deserts, and oases,”
Journal of Theoretical Biology , vol. 372, pp. 65–73, 2015.[29] G. Bel, B. Munsky, and I. Nemenman, “The simplicity of completion time distributions for commoncomplex biochemical processes,”
Physical biology , vol. 7, p. 016003, 2010.[30] S. Iyer-Biswas and Z. Anton, “First passage processes in cellular biology,” arXiv preprint , 2015.[31] W. Dai, A. M. Sengupta, and R. M. Levy, “First passage times, lifetimes, and relaxation times ofunfolded proteins,”
Physical Review Letters , vol. 115, p. 048101, 2015.[32] K. R. Ghusinga and A. Singh, “First-passage time calculations for a gene expression model,” in
IEEE53rd Annual Conference on Decision and Control (CDC) , pp. 3047–3052, 2014.[33] A. Singh and J. J. Dennehy, “Stochastic holin expression can account for lysis time variation in thebacteriophage λ ,” Journal of The Royal Society Interface , vol. 11, p. 20140140, 2014.[34] K. R. Ghusinga and A. Singh, “Theoretical predictions on the first-passage time for a gene expressionmodel,” in
IEEE 54th Annual Conference on Decision and Control (CDC) , pp. 3864–3869, 2015.[35] K. R. Ghusinga, J. J. Dennehy, and A. Singh, “Controlling noise in the timing of intracellular events:A first-passage time approach,” bioRxiv , p. 056945, 2016.[36] V. Shahrezaei and P. S. Swain, “Analytical distributions for stochastic gene expression,”
Proceedingsof the National Academy of Sciences , vol. 105, pp. 17256–17261, 2008.[37] H. Kuwahara, S. T. Arold, and X. Gao, “Beyond initiation-limited translational bursting: the effectsof burst size distributions on the stability of gene expression,”
Integr. Biol. , vol. 7, pp. 1622–1632,2015.[38] N. Kumar, A. Singh, and R. V. Kulkarni, “Transcriptional bursting in gene expression: Analyticalresults for general stochastic models,”
PLoS Comput Biol , vol. 11, pp. 1–22, 2015.[39] J. Paulsson, “Models of stochastic gene expression,”
Physics of Life Reviews , vol. 2, pp. 157–175,2005.[40] O. G. Berg, “A model for the statistical fluctuations of protein numbers in a microbial population,”
Journal of Theoretical Biology , vol. 71, pp. 587–603, 1978.[41] D. R. Rigney, “Stochastic model of constitutive protein levels in growing and dividing bacterial cells,”
Journal of Theoretical Biology , vol. 76, pp. 453–480, 1979.[42] X. Zhong, “On inverse and generalized inverses of hessenberg matrices,”
Linear Algebra and its Ap-plications , vol. 101, pp. 167 – 180, 1988. 1543] Y. Ikebe, “On inverses of hessenberg matrices,”
Linear Algebra and its Applications , vol. 24, pp. 93 –97, 1979.[44] K. R. Ghusinga, P.-W. Fok, and A. Singh, “Optimal auto-regulation to minimize first-passage timevariability in protein level,” in
American Control Conference (ACC) , pp. 4411–4416, 2015.[45] X. Liao, L. Wang, and P. Yu,