Optimal regulation of protein degradation to
schedule cellular
events with precision
Abstract
An important occurrence in many cellular contexts is the crossing of a prescribed threshold by a regulatory protein. The timing of such events is stochastic as a consequence of the innate randomness in gene expression. A question of interest is to understand how gene expression is regulated to achieve precision in event timing. To address this, we model event timing using the firstpassage time framework  a mathematical tool to analyze the time when a stochastic process first crosses a specific threshold. The protein evolution is described via a simple stochastic model of gene expression. Moreover, we consider the feedback regulation of protein degradation to be a possible noise control mechanism employed to achieve the precision. Exact analytical formulas are developed for the distribution and moments of the firstpassage time. Using these expressions, we investigate for the optimal feedback strategy such that noise (coefficient of variation squared) in event timing is minimized around a given fixed mean time. Our results show that the minimum noise is achieved when the protein degradation rate is zero for all protein levels. Lastly, the implications of this finding are discussed.
I Introduction
The process of gene expression is central to a cell’s function: the genetic information is used to make mRNA molecules which are further translated into proteins, the molecular machines to carry out different tasks. An important step in execution of many cellular functions is when a protein attains an effective level [1]. For example, in gene regulatory pathways, the protein expressed by one gene regulates expression of a downstream gene once its level crosses a certain threshold [2, 3, 4, 5]. Other examples include cellfate decisions which are made when specific regulatory proteins cross their respective critical levels [5, 6], and temporal order of gene activation which involves different thresholds of a given regulatory protein [7, 8, 9, 10].
The probabilistic nature of biochemical reactions and the fact that copy number of the constituents is small make expression of a gene a stochastic process [11, 12, 13, 14, 15, 16]. As a result, a population of cells induced at the same time sees difference in the times at which a critical protein level is attained in individual cells. While the stochasticity in timing does offer an advantage in terms of phenotype variability before a cell commits to an irreversible cell fate [17], in some cases precision in timing is paramount [10]. Although how cells achieve precision in timing is not understood very well, we expect that the gene expression might be regulated for this purpose. How some commonly found regulation motifs affect the noise in gene expression has been examined in several studies [18, 9, 19, 20, 21, 22, 23, 24].
In this paper, the regulation mechanism we analyze is degradation control of the protein via a feedback mechanism [19, 24]. Specifically, we investigate what form of feedback control of protein degradation results in minimum noise in event timing (quantified as the coefficient of variation squared, ), given a fixed mean time. The event time is modeled as a firstpassage time (), i.e., the first time at which a prescribed protein level is achieved. Previously, we have studied a similar question in [25, 26] assuming the regulation strategy to be selfregulation of transcription rate. In [25], a expression of a gene in bursts was considered along with the assumption that protein degradation rate is zero. Further, in [26] we revisited the same question assuming a constant degradation of the protein but simplified the production in bursts as birthdeath process. In another work [27], we considered both constant protein production in bursts and constant decay and developed the formulas for for the model to predict how changing various model parameters affect the statistics. This paper extends the calculations when the degradation rate is an arbitrary function of protein count. Employing exact solutions and semianalytical approach, we show that noise in is minimized when the protein does not degrade at all.
Remainder of the paper is organized as follows. In section II, we formulate the stochastic gene expression model wherein a protein is produced in bursts and its degradation is a function of the protein level. In the next section, for this gene expression model is determined. Section IV deals with determining the moments of , particularly the expressions of first two moments which are required to compute the noise (). Section V discusses the optimal regulation strategy that minimizes noise in . The results are discussed next in section VI. Some calculations to support the arguments in the main paper are provided in the appendices.
Ii Stochastic Gene Expression Model
Let denote the protein level at time . We consider expression of a gene as depicted in Fig. 1. We assume that the promoter is always active and consider the four fundamental processes, namely, transcription (making of mRNAs from gene), translation (production of protein from a mRNA), mRNA degradation and protein degradation. Also, the feedback regulation of protein degradation is formulated by assuming that if the protein level , then the degradation rate of one protein molecule is given by where represents an arbitrary function of .
The mRNA’s dynamics can be ignored in this model by assuming that its halflife is considerably smaller than that of the protein  an assumption which is true in most cases [28, 29, 30, 31]. This enables us to reduce the model to what is called a bursty birthdeath process wherein each transcription event creates a mRNA molecule which degrades immediately after synthesizing a burst of protein molecules [32, 30, 29]. The burst size is assumed to follow a geometric distribution which is consistent with previous theoretical and experimental studies [28, 33, 32, 34, 35, 36, 29, 30, 31]. Mathematically, the probabilities of having an arrival of a burst or a protein degradation in a small time interval are given by
(1a)  
(1b) 
Here is the transcription rate, denotes the geometrically distributed burst size while represents probability. The probability mass function of is given by:
(2) 
The mean burst size (average number of protein molecules produced in one mRNA lifetime), , can be expressed as:
(3) 
where and respectively denote the translation rate and the mRNA degradation rate. We are interested in determining the firsttime at which the aforementioned model reaches a fix threshold . In the next section, we present these calculations.
Iii First–passage time calculations
The firstpassage time for the random process describing protein level to cross a threshold is defined as:
(4) 
To compute , we employ same approach we used in our previous work [26, 27]. This involves considering a particle hopping on an integer lattice. A site on the lattice represents the values protein level can take. The forward jumps correspond to a birth in random bursts (protein production) while a backward jump corresponds to death (protein degradation). The sites with values greater than or equal to are considered to be absorbing, i.e., the process stops once the particle reaches one of these sites. For such a process, the probability that the threshold is crossed for the first time in a infinitesimal small time interval can be computed as
(5) 
Denoting the probability density function of by , we can write
(6a)  
(6b)  
(6c) 
Here we have made use of which comes from the fact that distribution of is geometric with parameter (see equation (2)). Equation (6c) can also be written as multiplication of two vectors and as
(7) 
where the row vector and the column vector are given as follows
(8)  
(9) 
In the expression for we have used the notation which stands for . To determine , we write the chemical master equation (or the forward Kolmogorov equation) for the bursty birthdeath process with absorbing state[37, 26, 27]:
(10a)  
(10b)  
(10c) 
These equations can be expressed as where elements of matrix are given by
(11) 
The above system of differential equations has the following solution
(12) 
where is the initial probability distribution (as at ). Using the expression of , the distribution of FPT can be written as
(13) 
Equation (13) gives the probability density function of in terms of known matrices , and . This expression can be generalized or simplified for a variety of cases such as a feedback regulation of the transcription rate, a different distribution of the burst size, a birthdeath process, etc. For each of these cases, the matrix and will change. Furthermore, can also be generalized to other distributions if the initial protein count is not assumed to be zero. In the next section, we use the expression in (13) to determine expressions of the first two moments of .
Iv Moments of FPT
In this section, we first give the expression of a general moment of in terms of the known matrices , and . Using this result, we further compute the exact formulas for the first two moments of .
The moment of can be computed as follows:
(14a)  
(14b) 
The matrix is a fullrank (invertible) Hurwitz matrix (see Appendices A, B). These properties can be used to get the following expression for a order moment [27]
(15) 
Using equation (15), we can get formulas for mean and second order moment of .
Iva Mean
The expression of mean is given by
(16) 
As discussed in Appendix B, we have where and are respectively given by equations (34) and (31). Therefore, we have
(17) 
The reason for not writing one of the in terms of and is that has a much simpler expression as described in Appendix C. We further observe that since has zero elements excepts for the first one, is just the first column of . Therefore is given by
(18a)  
(18b) 
By equation (39d), we have . Therefore the mean is negative of summation of the elements of the vector and is given by
(19) 
IvB Second order moment
The calculations of the second order moment are similar to those of the mean . For this reason, we skip the detailed steps and only provide the final formula.
(20a)  
(20b)  
the terms denoted by are given by  
(20c) 
Having developed expressions of the moments of , we next investigate the optimal feedback regulation of the protein degradation rate that minimizes noise in .
V Optimal regulation of protein degradation
We begin with considering an open loop control wherein the protein degradation rate is constant, i.e., . Using the formulas for the moments of that we developed in the previous section, we plot the of as the degradation rate is varied while keeping the mean constant with a simultaneous variation in the transcription rate .
As shown in Fig. 2, the of increases as is increased. That is, if there is no feedback regulation of the protein degradation, the best strategy is to have a zero degradation of the protein. It remains to be seen whether or not we can do better by having a feedback control of the protein degradation rate.
We explore this by first analyzing it for a birthdeath process which is a simplification of the bursty birthdeath process when the mean burst size [26].
Va Optimal feedback regulation of protein degradation for a birthdeath process
To analyze the birthdeath process, we use the formulas of moments developed in[38]. We would like to point out that in the limit or by carrying out calculation in the same fashion as this paper, we can derive moments (see Appendix D). The reason of going ahead with the formulas given in [38] is that they are easier to analyze in the present context. The moments are given by
(21)  
(22) 
We claim that the coefficient of variance squared (, defined as variance over mean squared) is lower bounded for a birthdeath process. The result is presented as a theorem.
Theorem 1
For a birthdeath process with a constant birth rate and protein level dependent death rate , the of is lower bounded by where is the threshold. The lower bound is achieved only when the degradation rates are zero.
We first show the equality. Let for . In this case, we have and . Therefore, . Next, we consider the case when . Using the CauchySchwarz inequality, we have
(23) 
The equality above holds only if
(24) 
which is true only in the limit when . Therefore, when , we have
(25a)  
(25b)  
(25c)  
(25d) 
This concludes the proof.
The above result shows that no matter what feedback strategy in degradation is employed, the in timing will always be greater than which is achieved when the degradation rates for all protein levels is zero. Next, we explore the optimal feedback for the bursty birthdeath process.
VB Optimal feedback regulation of protein degradation for a bursty birthdeath process
Intuitively, the result we obtained in the previous section for the birthdeath process should not change for the bursty birthdeath process. This is because the difference between both these processes is in the birth events and the feedback in death should affect them in similar fashion. We show numerical results to substantiate this argument. We assume that as the protein level is increased, it represses its degradation activity via a negative feedback given by
(26) 
Here represents the Hill coefficient, is the maximum degradation rate, is protein level and is referred to as the feedback strength. In Fig. 3, we show of as is varied for given value . The mean is kept constant by changing the transcription rate. It can be seen that as the feedback becomes stronger, the decreases. That is, faster the protein degradation rate decreases, lower the gets.
The global minimum in is the case when the protein does not degrade. For nonzero values of , the approaches this minimum value for very strong negative feedbacks (). We have not shown results for the positive feedback to the degradation rate but have checked that the gets worse as the feedback strength is increased  a result that is expected.
Vi Discussion
In this work, we characterized the time taken by a protein to reach a certain threshold using the wellknown firstpassage time framework. Further, we investigated if precision in can be achieved by regulating the degradation rate of the protein while keeping the mean fixed. We showed that this is not possible and the best strategy rather is to have no degradation of the protein.
This finding complements the results of our previous works in [25, 26] wherein we considered the regulation of transcription rate. We showed that when protein does not degrade, the best regulation strategy of the transcription rate is an open loop system [25]. In contrast, if a constant rate of protein decay is considered then the optimal transcription rates seem to be a mixture of both positive and negative feedbacks [26]. Combining these two results, it can be suggested that the best strategy for achieving minimum noise in is to have no degradation of the protein and have open loop (no feedback) production.
Questions similar to the one addressed in this work can also be asked about other commonly found gene regulation motifs such as feedforward loops [9]. Further, this work is confined to analyzing a reduced model of gene expression, ignoring the effects of promoter switching [39], and the parameter regimes when mRNA and protein halflives are comparable. In future, we would aim to investigate along these lines.
a Proof that is a Hurwitz matrix
To prove that the matrix defined in equation (11) is a Hurwitz matrix, we show that it satisfies the following conditions stated in [40]:

The diagonal elements for ,

.
Note that . Therefore, the first requirement is fulfilled. Also,
(27a)  
(27b) 
That is, satisfies the second condition as well. Hence, is Hurwitz.
B Calculations to determine
We use the same procedure as we did in our previous work [27]. We decompose as where consists of only the nondegradation rates terms. Specifically, and are respectively given by
(28) 
(29) 
Note that inverse of can be written as
(30) 
The term is easy to determine and is given by
(31) 
Further, we also need inverse of in order to determine . To this end, let us compute :
(32a)  
(32b)  
(32c) 
The matrix is a bidiagonal matrix with its diagonal elements for . The super diagonal elements are given by for . Using the result for inverse of a bidiagonal matrix derived in [41], we can write the element of as follows
(33) 
In matrix form, can be written as
(34) 
Thus, we have determined the matrices and which can be used to compute .
C Expression of
In [27], where both production and degradation of protein do not have feedback regulation, we show that simplifies to a row vector with all of its elements being . It turns out that even in the case when the protein degradation has a feedback, the same result is true. The calculations are as follows. Consider two matrices and such that where is a matrix
(35) 
while is a matrix
(36) 
We evoke the matrix inversion lemma to write as
(37a)  
(37b)  
This implies  
(37c) 
Let us compute :
(38a)  
(38b) 
Therefore, we can conclude that is in fact equal to which can be calculated by multiplying and as shown below.
(39a)  
(39b)  
(39c)  
(39d) 
D First two moments of for a birthdeath process
The matrixbased approach considered in this paper can be used to compute the formulas for moments of the firstpassage time. These are given by:
(40a)  
(40b)  
(40c) 
Though we have skipped the detailed steps here, these results are equivalent to the results derived in [26] and [38]. More specifically, if , we get the results of [26]. The results in [38] are more general and simplify to our results if the birth rate is made constant while the death rate is taken as instead of .
References
 [1] H. H. McAdams and A. Arkin, “Stochastic mechanisms in gene expression,” Proceedings of the National Academy of Sciences, vol. 94, pp. 814–819, 1997.
 [2] T. Lamb, “Gain and kinetics of activation in the Gprotein cascade of phototransduction,” Proceedings of the National Academy of Sciences, vol. 93, no. 2, pp. 566–570, 1996.
 [3] M. C. Gustin, J. Albertyn, M. Alexander, and K. Davenport, “Map kinase pathways in the yeast saccharomyces cerevisiae,” Microbiology and Molecular biology reviews, vol. 62, no. 4, pp. 1264–1300, 1998.
 [4] E. H. Davidson, J. P. Rast, P. Oliveri, A. Ransick, C. Calestani, C.H. Yuh, T. Minokawa, G. Amore, V. Hinman, C. ArenasMena, et al., “A genomic regulatory network for development,” science, vol. 295, no. 5560, pp. 1669–1678, 2002.
 [5] J. E. Ferrell and E. M. Machleder, “The biochemical basis of an allornone cell fate switch in xenopus oocytes,” Science, vol. 280, no. 5365, pp. 895–898, 1998.
 [6] P. Ingham and J. Smith, “Crossing the threshold,” Current Biology, vol. 2, no. 9, pp. 465 – 467, 1992.
 [7] A. Zaslaver, A. E. Mayo, R. Rosenberg, P. Bashkin, H. Sberro, M. Tsalyuk, M. G. Surette, and U. Alon, “Justintime transcription program in metabolic pathways,” Nature genetics, vol. 36, no. 5, pp. 486–491, 2004.
 [8] M. Fujita, J. E. GonzálezPastor, and R. Losick, “Highand lowthreshold genes in the spo0a regulon of bacillus subtilis,” Journal of bacteriology, vol. 187, no. 4, pp. 1357–1368, 2005.
 [9] U. Alon, An introduction to systems biology: design principles of biological circuits. CRC press, 2006.
 [10] A. Amir, O. Kobiler, A. Rokney, A. B. Oppenheim, and J. Stavans, “Noise in timing and precision of gene activities in a genetic cascade,” Molecular Systems Biology, vol. 3, 2007.
 [11] 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.
 [12] J. M. Raser and E. K. O’Shea, “Noise in gene expression: origins, consequences, and control,” Science, vol. 309, pp. 2010–2013, 2005.
 [13] J. Yu, J. Xiao, X. Ren, K. Lao, and X. S. Xie, “Probing gene expression in live cells, one protein molecule at a time,” Science, vol. 311, pp. 1600–1603, 2006.
 [14] A. Raj and A. van Oudenaarden, “Nature, nurture, or chance: stochastic gene expression and its consequences,” Cell, vol. 135, pp. 216–226, 2008.
 [15] M. Kærn, T. C. Elston, W. J. Blake, and J. J. Collins, “Stochasticity in gene expression: from theories to phenotypes,” Nature Reviews Genetics, vol. 6, pp. 451–464, 2005.
 [16] A. Singh and M. Soltani, “Quantifying intrinsic and extrinsic variability in stochastic gene expression models,” PloS One, vol. 8, p. e84301, 2013.
 [17] E. Yurkovsky and I. Nachman, “Event timing at the single–cell level,” Briefings in Functional Genomics, vol. 12, pp. 90–98, 2013.
 [18] A. Becskei and L. Serrano, “Engineering stability in gene networks by autoregulation,” Nature, vol. 405, pp. 590–593, 2000.
 [19] H. ElSamad and M. Khammash, “Regulated degradation is a mechanism for suppressing stochastic fluctuations in gene regulatory networks,” Biophysical journal, vol. 90, pp. 3749–3761, 2006.
 [20] U. Alon, “Network motifs: theory and experimental approaches,” Nature Reviews Genetics, vol. 8, pp. 450–461, 2007.
 [21] A. Singh and J. P. Hespanha, “Evolution of gene autoregulation in the presence of noise,” Systems Biology, IET, vol. 3, pp. 368–378, 2009.
 [22] Y. Tao, X. Zheng, and Y. Sun, “Effect of feedback regulation on stochastic gene expression,” Journal of Theoretical Biology, vol. 247, pp. 827–836, 2007.
 [23] A. Singh and J. P. Hespanha, “Optimal feedback strength for noise suppression in autoregulatory gene networks,” Biophysical Journal, vol. 96, pp. 4013–4023, 2009.
 [24] G. Chalancon, C. N. Ravarani, S. Balaji, A. MartinezArias, L. Aravind, R. Jothi, and M. M. Babu, “Interplay between gene expression noise and regulatory network architecture,” Trends in genetics, vol. 28, pp. 221–232, 2012.
 [25] K. R. Ghusinga and A. Singh, “Firstpassage time calculations for a gene expression model,” in IEEE 53rd Annual Conference on Decision and Control (CDC), 2014, pp. 3047–3052.
 [26] K. R. Ghusinga, P.W. Fok, and A. Singh, “Optimal autoregulation to minimize firstpassage time variability in protein level,” in American Control Conference (ACC), 2015, pp. 4411–4416.
 [27] K. R. Ghusinga and A. Singh, “Theoretical predictions on the firstpassage time for a gene expression model,” in IEEE 54th Annual Conference on Decision and Control (CDC) (accepted), 2015.
 [28] J. Paulsson, “Models of stochastic gene expression,” Physics of Life Reviews, vol. 2, pp. 157–175, 2005.
 [29] A. Singh and P. Bokes, “Consequences of mrna transport on stochastic variability in protein levels,” Biophysical journal, vol. 103, pp. 1087–1096, 2012.
 [30] A. Singh, B. S. Razooky, R. D. Dar, and L. S. Weinberger, “Dynamics of protein noise can distinguish between alternate sources of geneexpression variability,” Molecular systems biology, vol. 8, p. 607, 2012.
 [31] A. Singh, B. Razooky, C. D. Cox, M. L. Simpson, and L. S. Weinberger, “Transcriptional bursting from the hiv1 promoter is a significant source of stochastic noise in hiv1 gene expression,” Biophysical journal, vol. 98, pp. L32–L34, 2010.
 [32] V. Shahrezaei and P. S. Swain, “Analytical distributions for stochastic gene expression,” Proceedings of the National Academy of Sciences, vol. 105, pp. 17 256–17 261, 2008.
 [33] N. Friedman, S. Vardi, M. Ronen, U. Alon, and J. Stavans, “Precise temporal modulation in the response of the sos dna repair network in individual bacteria,” PLoS biology, vol. 3, no. 7, p. e238, 2005.
 [34] 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.
 [35] 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.
 [36] V. Elgart, T. Jia, A. T. Fenley, and R. Kulkarni, “Connecting protein and mrna burst distributions for stochastic models of gene expression,” Physical biology, vol. 8, p. 046001, 2011.
 [37] D. A. McQuarrie, “Stochastic approach to chemical kinetics,” Journal of applied probability, vol. 4, pp. 413–478, 1967.
 [38] O. Jouini and Y. Dallery, “Moments of first passage times in general birth–death processes,” Mathematical Methods of Operations Research, vol. 68, pp. 49–76, 2008.
 [39] A. Singh, C. A. Vargas, and R. Karmakar, “Stochastic analysis and inference of a twostate genetic promoter model,” American Control Conference (ACC), pp. 4563–4568, 2013.
 [40] X. Liao, L. Wang, and P. Yu, Stability of Dynamical Systems. Elsevier, 2007, vol. 5.
 [41] G. Chatterjee, “Negative integral powers of a bidiagonal matrix,” Mathematics of Computation, vol. 28, pp. 713–714, 1974.