Stochastic Model Formulation

Consider a gene that is switched on at time t=0 and begins to express a timekeeper protein. The intracellular event of interest is triggered once the protein reaches a critical level in the cell. We describe a minimal model of gene expression that assumes translation in bursts and incorporates feedback regulation by considering the transcription rate as a function of the protein level (Fig. 1). More precisely, if x(t)∈{0,1,…} denotes the level of a protein in a single cell at time t, then the gene is transcribed at a Poisson rate ki when x(t)=i. Any arbitrary form of feedback can be realized by appropriately defining ki as a function of i. For example, increasing (decreasing) ki’s correspond to a positive (negative) feedback loop in protein production, and a fixed transcription rate implies no feedback. The translation burst approximation is based on assuming short-lived mRNAs, that is, each mRNA degrades instantaneously after producing a burst of B protein molecules. In agreement with experimental and theoretical studies (38⇓–40), B is assumed to follow a geometric distributionℙ(B=j)=bj(b+1)j+1, b∈(0,∞),j={0,1,2,…},[1]where b denotes the mean protein burst size and the symbol ℙ is the notion for probability. Finally, each protein molecule degrades with a constant rate γ. The time evolution of x(t) is described through the following probabilities of occurrences of burst and decay events in the next infinitesimal time (t,t+dt]:ℙ(x(t+dt)=i+B|x(t)=i)=kidt,[2a]ℙ(x(t+dt)=i−1|x(t)=i)=iγdt.[2b]Note that in this representation of gene expression as a bursty birth–death process, the mRNA transcription rate ki is the burst arrival rate, whereas the rate at which proteins are translated from an mRNA determines the mean protein burst size b. Next, we formulate event timing through the FPT framework.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 1.

Modeling event timing as an FPT problem. (Left) Schematic depicting a gene transcribing mRNAs, which are further translated into proteins. The transcription rate is assumed to be regulated by the protein level, creating a feedback loop. (Right) The timing of an intracellular event is formulated as the FPT for the protein level to reach a critical threshold. Sample trajectories for the protein level over time obtained via Monte Carlo simulations are shown, and they cross the threshold at different times due to stochasticity in gene expression. The histogram of the event timing is shown on the top.

Computing Event Timing Distribution

The time to an event is the FPT for x(t) to reach a threshold X starting from a zero initial condition x(0)=0 (Fig. 1). It is mathematically described by the following random variable:T:=min{t:x(t)≥X|x(0)=0},[3]and can be interpreted as the time taken by a random walker to first reach a defined point. For the bursty birth–death process in Eq. 2, the probability density function (pdf) of the FPT is given byfT(t)=∑i=0X−1ki(bb+1)X−ipi(t),[4]where pi(t)=ℙ(x(t)=i) (SI Appendix, section S1).

The FPT pdf in Eq. 4 can be compactly written as product of two vectors:fT(t)=𝐔⊤𝐏(t),[5a]𝐔=[k0(bb+1)Xk1(bb+1)X−1⋯kX−1bb+1]⊤,[5b]𝐏(t)=[p0(t)p1(t)⋯pX−1(t)]⊤.[5c]Here, the dynamics of 𝐏(t) can be written as a linear system𝐏˙=𝐀𝐏[6]derived from the chemical master equation corresponding to the bursty birth–death process (SI Appendix, section S1) (23, 41). It turns out that in this case the matrix 𝐀 is a Hessenberg matrix whose ith row and jth column element aij is given byaij={0,j>i+1(i−1)γ,j=i+1−ki−1bb+1−(i−1)γ,j=iki−1bi−j(b+1)i−j+1,j<i[7]i,j∈{1,…,X}. Solving Eq. 6 and using Eq. 5a yields the following pdf for the FPT:fT(t)=𝐔⊤𝐏(t)=𝐔⊤exp(𝐀t)𝐏(0),[8]where 𝐏(0)=[10⋯0]⊤ is vector of probabilities at t=0 that follows from x(0)=0. Although this pdf provides complete characterization of the event timing, we are particularly interested in the lower-order statistical moments of FPT. Next, we exploit the structure of matrix 𝐀 to obtain analytical formulas for the first- and second-order moments of FPT.

Moments of the FPT

From Eq. 8, the mth-order uncentered moment of the FPT is given by⟨Tm⟩=𝐔⊤(∫0∞tmexp(𝐀t)dt)𝐏(0)[9a]=(−1)m+1m!𝐔⊤(𝐀−1)m+1𝐏(0).[9b]Here, in computing the above integral, we used the fact that the matrix 𝐀 is full-rank with negative eigenvalues (SI Appendix, section S2). One can also find explicit formulas for the first two moments as series summations in terms of event threshold X, mean burst size b, protein decay rate γ, and transcription rates ki, 0≤i≤X−1 (SI Appendix, section S3).

Analysis of the FPT moments in some limiting cases gives important insights (SI Appendix, section S4). For the simplest case of a stable protein (γ=0) and a constant transcription rate (no feedback; ki=k), the moment expressions simplify to⟨T⟩=1k(Xb+1)≈Xbk,CVT2=b2+X+2bX(b+X)2≈1+2bX,[10]where CVT2 represents the noise in FPT as quantified by its coefficient of variation squared (variance/mean²;⟨T2⟩/⟨T⟩2−1). The approximate formulas in Eq. 10 are valid for a high event threshold compared with the mean protein burst size (X/b≫1). The mean FPT formula can be interpreted as the time taken to reach X with an accumulation rate bk. Further, X/b represents the average number of burst events required for the protein level to cross the threshold, and increasing X/b leads to noise reduction through more efficient averaging of the bursty process. One can also gain important insights, such as, the noise in FPT is invariant of the transcription rate k. Therefore, ⟨T⟩ and CVT2 can be independently tuned—increasing the event threshold and/or reducing the burst size will lower the noise level. Once CVT2 is sufficiently reduced, k can be altered to obtain a desired mean event timing.

Interesting features of the FPT statistics are revealed when γ≠0 (unstable protein) is considered with a constant transcription rate ki=k. In this case, the expressions of FPT moments are quite involved (SI Appendix, section S4), and we investigate the effect of various parameters numerically. It turns out that changing the event threshold leads to a U-shaped profile for CVT2, where noise in event timing first decreases and then increases (Fig. 2). Intuitively, when γ≠0, the protein level approaches a steady-state level xss=kb/γ. When the event threshold X is sufficiently below xss, CVT2 reduces with increasing X, similar to the γ=0 case. As X approaches xss, the protein trajectories start saturating and crossing the threshold becomes a noise-driven event. This results in an increase in CVT2 and ultimately leads to CVT2→1 as X≫xss. Recall that the coefficient of variation of an exponentially distributed random variable is exactly equal to one. Thus, when X is much larger than xss, the timing process becomes memoryless, yielding exponentially distributed FPTs. The minimum value of CVT2 is achieved at an intermediate threshold level X≈xss/2 for a birth–death process (b→0), and the dip in the U-shape shifts to the right as the protein expression becomes more bursty (Fig. 2). Another interesting point to note is that whereas increasing b increases noise in event timing when X≪xss, it has a contrasting effect when X≫xss, where increasing b can sometimes reduce CVT2. Next, we explore how feedback regulation of the transcription rate affects noise in timing, for a given X and b.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

Relative positions of the event threshold and the steady-state protein level determine noise in timing. Noise in the FPT (CVT2) is plotted as a function of the event threshold (X) for different mean protein burst sizes b and decay rates γ. For γ≠0, the noise is high at small values of X and decreases with increasing X, which is similar to the CVT2 versus X trend when γ=0 (dashed lines). After attaining a minima at an intermediate value of X, CVT2 increases with further increase in X. The minimum value of CVT2 is achieved at X≈xss/2 for a birth–death process, and this optimal point shifts to the right as b is increased while keeping xss fixed by commensurate change in k. The steady-state protein level and the decay rate are taken as xss=300 molecules, and γ=0.01 per min.

Optimal Feedback Strategy

Having derived the FPT moments, we investigate optimal forms of transcriptional feedback that schedule an event at a given time with the lowest CVT2. Because ⟨T⟩ is assumed to be fixed, minimizing CVT2 is equivalent to minimizing ⟨T2⟩. Thus, the problem mathematically corresponds to a constraint optimization problem: Find transcription rates k0,k1,⋯,kX−1 that minimize ⟨T2⟩ for a fixed ⟨T⟩. We first consider a stable protein whose half-life is much longer than the event timescale, and, hence, degradation can be ignored (γ=0).

Optimal Feedback for a Stable Protein.

When the protein of interest does not decay (γ=0), the expressions for the FPT moments take much simpler forms:⟨T⟩=1k0+1b∑i=0X−11ki,[11a]⟨T2⟩=2b2(τ0bk0+∑i=0X−1τiki),τi:=bki+∑j=iX−11kj.[11b]Note that, in Eq. 11a, the contribution of k0 (transcription rate when there is no protein) differs from the other transcription rates ki, i∈{1,2,⋯,X−1}. For instance, when the event threshold is large compared with the mean burst size (X≫b), then the term 1/k0 can be ignored and ⟨T⟩≈∑i=0X−11/bki. In contrast, if the burst size is large (b≫X) then ⟨T⟩≈1/k0, because a single burst event starting from zero protein molecules is sufficient for threshold crossing. A similar observation for different contributions of k0 can be made about Eq. 11b.

It turns out that, for these simplified formulas, the problem of minimizing ⟨T2⟩ given ⟨T⟩ can be solved analytically using the method of Lagrange multipliers (SI Appendix, section S5). The optimal transcription rates are given byk0=1+b1+2b2b+Xb⟨T⟩, ki=1+2b1+bk0, 1≤i≤X−1,[12]and all rates are equal to each other except for k0. Intuitively, the difference for k0 comes from the fact that it contributes differently to the FPT moments compared with other rates. Note that for a small mean burst size (b≪1), k0=ki, whereas k0=ki/2 for a sufficiently large b. Despite this slight deviation in k0, for the purposes of practical implementation, the optimal feedback strategy in this case is to have a constant transcription rate (i.e., no feedback in protein expression).

We tested the above result for a more complex stochastic gene expression model that explicitly includes mRNA dynamics via Monte Carlo simulations (Fig. 3). For ease of implementation, the feedbacks are assumed to be linear:ki=c1+c2i, i∈{0,1,…},[13]where c2=0 represents no feedback, and c2>0(c2<0) denotes a positive (negative) feedback. In agreement with Eq. 12, a no-feedback strategy outperforms negative/positive feedbacks in terms of minimizing noise in FPT around a given mean event time. The qualitative shape of trajectories in Fig. 3 is determined by the feedback strategy used, with no feedback resulting in linear time evolution of protein levels. This provides an intriguing geometric interpretation of our results—an approximate linear path from zero protein molecules at t=0 to X molecules at time ⟨T⟩ provides the highest precision in timing. Next, we discuss the optimal feedback strategy when protein degradation is taken into consideration.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 3.

For a stable protein, no feedback provides the lowest noise in event timing for a fixed mean FPT. Protein trajectories obtained using the stochastic simulation algorithm for a stochastic gene expression model with positive feedback (Left), no feedback (Middle), and negative feedback (Right) (42). The threshold for event timing is assumed to be 500 protein molecules and feedback is implemented by assuming a linear form of transcription rates as ki=c1+c2i. The value of c2 is taken as 0.05min−1 for positive feedback and −0.05min−1 for negative feedback. For each feedback, the parameter c1 is taken such that the mean FPT is kept constant (40 min). The mRNA half-life is assumed to be 2.7min, and proteins are translated from mRNAs at a rate 0.5min−1, which corresponds to a mean burst size of b=2. Histograms on the top represent distribution of FPT from 10,000 Monte Carlo simulations.

Optimal Feedback for an Unstable Protein.

Now consider the scenario where protein degradation cannot be ignored over the event timescale (γ≠0). Unfortunately, the expressions of the FPT moments are too convoluted for the optimization problem to be solved analytically (SI Appendix, section S3, Eq. S3.10), and the effect of different feedbacks is investigated numerically.

We implement the feedbacks using physiologically relevant Hill functions, where the transcription rates for a negative feedback mechanism take the following form:ki=kmax1+(ci)H, i∈{0,1,…}.[14]Here H denotes the Hill coefficient, kmax corresponds to the maximum transcription rate, and c characterizes the negative feedback strength, with c=0 representing no feedback (29, 43). Similarly, a positive feedback is assumed to take the following form:ki=kmax(r+(1−r)(ci)H1+(ci)H)=kmaxr+(ci)H1+(ci)H.[15]Note that an additional parameter r∈(0,1), referred to as the basal strength, is introduced in Eq. 15. This is to ensure that the transcription rate in protein absence k0=kmaxr>0, and this is necessary to prevent protein levels from getting stuck at zero molecules.

To find the optimal feedback mechanism, our strategy is as follows: For given r and H, choose a certain feedback strength c in Eq. 14/Eq. 15, appropriately tune kmax for the desired mean event timing, and explore the corresponding noise in FPT as measured by its coefficient of variation squared CVT2. Counterintuitively, results show that for a given value of γ, a negative feedback loop in gene expression has the highest CVT2, and its performance deteriorates with increasing feedback strength (Fig. 4, Top). In contrast, CVT2 first decreases with increasing strength of the positive feedback and then increases after an optimal feedback strength is crossed (Fig. 4, Top). Thus, when the protein is not stable, precision in timing is attained by having a positive feedback in protein synthesis with an intermediate strength.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 4.

For an unstable protein, positive feedback provides the lowest noise in event timing for a fixed mean FPT. (Top) Noise in timing (CVT2) as a function of the feedback strength c for different control strategies. The value of kmax is changed in Eq. 14/Eq. 15 so as to keep ⟨T⟩=40mins fixed. The performance of the negative feedback worsens with increasing feedback strength. In contrast, positive feedback with an optimal value of c provides the highest precision in event timing. Other parameters used are γ=0.05min−1, X=500 molecules, H=1, b=2, and for positive feedback r=0.05. (Bottom) The minimum value of CVT2 obtained via positive feedback increases monotonically with the protein degradation rate. A smaller basal promoter strength r=0.01 in Eq. 15 gives better noise suppression than a larger value r=0.05. For comparison purposes, CVT2 obtained without any feedback (c=0), and a linear feedback with c1 and c2 in Eq. 13 chosen to minimize CVT2 for a given ⟨T⟩=40mins are also shown. The parameter values used are X=500 molecules, H=1, and b=2.

We next explore how the minimal achievable noise in event timing, for a fixed ⟨T⟩, varies with the protein decay rate γ. Our analysis shows that for a given basal strength r, the minimum CVT2 obtained via positive feedback increases monotonically with γ, and CVT2→1 as γ becomes large (Fig. 4, Bottom). A couple of interesting observations can be made from Fig. 4, Bottom: (i) The difference in CVT2 for optimal feedback and no feedback is indistinguishable when the protein is stable (γ=0) or highly unstable (γ→∞); (ii) for a range of intermediate protein half-lives the optimal feedback strategy provides better reduction of CVT2, as compared with no feedback regulation (which also corresponds to minimum CVT2 obtained via a negative feedback); (iii) lowering the basal strength r results in better performance in terms of noise suppression; and (iv) a linear feedback based on Eq. 13 outperforms feedbacks based on Hill functions and provides significantly lower levels of CVT2 for high protein decay rates. It is also worth pointing out that the qualitative shape of curves in Fig. 4, Bottom does not change for different values of event threshold X or mean burst size b (SI Appendix, section S6).

Why is positive feedback the optimal control strategy for ensuring precision in event timing? One way to understand this result is to consider the linear feedback form Eq. 13, in which case the mean protein levels evolve according to the following ordinary differential equation:dx(t)dt=b(c1+c2x)−γx, x(0)=0.[16]

Recall the geometric argument presented in Fig. 3, where an approximately linear path for the protein to reach the prescribed threshold in a given time provides the highest precision in event timing. Whereas no feedback (c2=0) and negative feedback (c2<0) in Eq. 16 will create nonlinear protein trajectories, choosing a positive value c2≈γ/b results in linear x(t), and hence minimal noise in event timing. Indeed, our detailed stochastic analysis shows that the optimal feedback strength that minimizes CVT2 in the stochastic model is qualitatively similar to c2≈γ/b (SI Appendix, section S6).