Regulation increases temporal precision

To investigate the effects of regulation on temporal precision, we consider the timing variance as a function of the parameters k and K, or μ and K. The special case of no regulation corresponds to the limits k → ∞ and K → 0 in the case of activation, or μ → ∞ and K → ∞ in the case of repression. In this case, the production of X occurs at the constant rate α. Reaching the threshold requires x_* sequential events, each of which occurs in a time that is exponentially distributed with mean 1/α. The total completion time for such a process is given by a gamma distribution with mean and variance [19]. Ensuring that requires α = x_*/t_*, for which the variance satisfies . This expression gives the timing variance for the unregulated process.

In Fig 2 we plot the scaled variance as a function of the parameters k and K, or μ and K, for cooperativity H = 3 (color maps). In the case of activation (Fig 2A), the variance decreases with increasing k and K. This means that the temporal precision is highest for an activator that accumulates quickly and requires a high abundance to produce X. In the case of repression (Fig 2B), the variance has a global minimum as a function of μ and K. This means that the temporal precision is highest for a repressor with a particular well-defined degradation rate and abundance threshold. Importantly, we see that for both activation and repression, the scaled variance can be less than one, meaning that regulation allows improvement of the temporal precision beyond that of the unregulated process. We have checked that this result holds for H ≥ 1.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Optimal regulatory strategies.

Timing variance as a function of the regulatory parameters reveals (A) a trajectory along which the variance decreases in the case of the activator and (B) a global minimum in the case of the repressor. White dashed line in A and white dot in B show the analytic approximations in Eqs 9 and 11, respectively. Parameters are N = 15 in B, and x_* = 15 and H = 3 in both.

https://doi.org/10.1371/journal.pcbi.1006201.g002

Optimal regulation balances regulator and target noise

To understand the dependencies in Fig 2, we develop analytic approximations. First, we assume that H → ∞, such that the regulation functions in Eqs 1 and 2 become threshold functions. In this limit, the production rate of X is zero if a < K or r > K, and α otherwise. The deterministic dynamics of X become piecewise-linear, (5) where t₀ is determined by either or according to Eqs 3 and 4. Then, to set α, we use the condition , which results in α = x_*/(t_* − t₀).

Lastly, we approximate the variance in the first-passage time using the variance in the molecule number and the time derivative of the mean dynamics [13]. Specifically, the timing variance of X arises from two sources: (i) uncertainty in the time when the regulator crosses its threshold K, which determines when the production of the target X begins, and (ii) uncertainty in the time when x crosses its threshold x_*, given that production begins at a particular time. The first source is regulator noise, and the second source is target noise. We estimate these timing variances from the associated molecule number variances, propagated via the time derivatives, (6) where y ∈ {a, r} denotes the regulator molecule number.

For the activator, which undergoes a pure production process with rate k, the molecule number obeys a Poisson distribution with mean kt. Therefore, the molecule number variance at time t₀ is . For the repressor, which undergoes a pure degradation process with rate μ starting from N molecules, the molecule number obeys a binomial distribution with number of trials N and success probability e^−μt. Therefore, the molecule number variance at time t₀ is . For the target molecule, which undergoes a pure production process with rate α starting at time t₀, the molecule number obeys a Poisson distribution with mean α(t − t₀). Therefore, the molecule number variance at time t_* is . Inserting these expressions into Eq 6, along with the derivatives calculated from Eqs 3–5 and the appropriate expressions for α and t₀, we obtain (7) (8) As a function of kt_* and K, the global minimum of Eq 7 occurs as kt_* → ∞ and K → ∞. The path of descent toward this minimum is given by differentiating with respect to K at fixed kt_* and setting the result to zero, which yields the curve (9) along which the variance satisfies (10) where the first case comes from the fact that K must be nonnegative. In contrast, the global minimum of Eq 8 occurs at finite μt_* and K: differentiating with respect to each and setting the results to zero gives the values (11a) (11b) (12) where we have assumed that K/N ≪ 1 (see Materials and methods), which is justified post-hoc by Eq 11a.

These analytic approximations are compared with the numerical results for the activator in Fig 2A (white dashed line, Eq 9) and for the repressor in Fig 2B (white circle, Eq 11). In Fig 2A we see that the global minimum indeed occurs as kt_* → ∞ and K → ∞, and the predicted curve agrees well with the observed descent. In Fig 2B we see that the predicted global minimum lies very close to the observed global minimum. We have also checked along specific slices in the (K, kt_*) or (K, μt_*) plane and found that the analytic approximations generally differ from the numerical results by about 10% or less, despite the fact that the approximations take H → ∞ whereas the numerics in Fig 2 use H = 3.

The success of the approximations means that Eq 6 describes the key mechanism leading to the optimal temporal precision. Eq 6 demonstrates that the optimal regulatory strategy arises from a tradeoff between minimizing regulator and target noise. On the one hand, minimizing only the regulator noise would require that the regulator cross its threshold K with a steep slope and therefore at an early time, meaning that the target molecule would be effectively unregulated and would increase linearly over time. On the other hand, minimizing only the target noise would require that the regulator cross its threshold only shortly before the target time t_*, such that the target molecule would cross its threshold x_* with a steep slope , leading to a highly nonlinear increase of the target molecule over time. In actuality, the optimal strategy is somewhere in between, with the regulator crossing its threshold at some intermediate time t₀, and the target molecule exhibiting moderately nonlinear dynamics as in Fig 1C and 1D.

Eqs 10 and 12 demonstrate that the timing variance is small for large kt_*/x_* in the case of activation, and small for large N/x_* in the case of repression. This makes intuitive sense because each of these quantities scales with the number of regulator molecules: k is the production rate of activator molecules, while N is the initial number of repressor molecules. To make this intuition quantitative, we define a cost as the time-averaged number of regulator molecules, (13) (14) where the second steps follow from Eqs 3 and 4. We see that, indeed, 〈a〉 scales with k, and 〈r〉 scales with N. Thus, Eqs 10 and 12 demonstrate that increased temporal precision comes at a cost, in terms of the number of regulator molecules that must be produced.

Model predictions are consistent with neuroblast migration data

We test our model predictions using data from neuroblast cells in C. elegans larvae [5]. During C. elegans development, particular neuroblast cells migrate from the posterior to the anterior of the larva. It has been shown that the migration terminates not at a particular position, but rather after a particular amount of time, and that the termination time is controlled by a temporal increase in the expression of the mig-1 gene [5]. Since mig-1 is known to be subject to regulation [21], we investigate the extent to which the dynamics of mig-1 can be explained by the predictions of our model.

Fig 3A shows the number x of mig-1 mRNA molecules per cell as a function of time t, obtained by single-molecule fluorescent in situ hybridization (from [5]). We analyze these data in the following way (see Materials and methods for details). First, noting that the dynamics are nonlinear, we quantify the linearity using the area under the curve, normalized by that for a perfectly linear trajectory x_*t_*/2, (15) By this definition, ρ = 1 for perfectly linear dynamics, and ρ → 0 for maximally nonlinear dynamics (a sharp rise at t_*). Then, we estimate x_*, t_*, and the timing variance from the data. Specifically, migration is known to terminate between particular reference cells in the larva [5], which gives an estimated range for the termination time t_*. This range is shown in magenta in Fig 3A and corresponds to a threshold within the approximate range 10 ≤ x_* ≤ 25. Therefore, we divide the x axis into bins of size Δx, choose bin midpoints x_* within this range, and for each choice compute the mean t_* and the variance of the data in that bin. Fig 3B shows the average and standard deviation of results for different values of x_* and Δx (blue circle).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Model predictions agree with neuroblast migration data.

(A) Number of mig-1 mRNA molecules per cell as a function of time t, obtained by single-molecule fluorescent in situ hybridization, from [5]. Magenta shows approximate range of times when cell migration terminates. Black lines show mean (dashed) and standard deviation σ_d of cell division times (black points). (B) Timing variance vs. linearity of x(t), both for experimental data in A (blue circle) and our model (curves, Eqs 16 and 17). Data analyzed using ranges of threshold 10 ≤ x_* ≤ 25 and bin size 3 ≤ Δx ≤ 12; error bars show standard deviations of these results. We see that for sufficiently large cost 〈a〉/x_* or 〈r〉/x_*, model predictions agree with experimental data point.

https://doi.org/10.1371/journal.pcbi.1006201.g003

The experimental data point in Fig 3B exhibits two clear features: (i) the dynamics are nonlinear (ρ is significantly below 1), and (ii) the timing variance is low ( is significantly below 1). Neither feature can be explained by a model in which the production of x is unregulated, since that would correspond to a linear increase of molecule number over time (ρ = 1) and a timing variance that satisfies (square in Fig 3B). Furthermore, since auto-regulation has been shown to generically increase timing variance beyond the unregulated case [12], it is unlikely that feature (ii) can be accounted for by a model with auto-regulation alone. Can these data be accounted for by our model with regulation?

To address this question we calculate ρ and from our model. For simplicity, we focus on the analytic approximations in Eqs 7 and 8, since they have been validated in Fig 2. In these approximations, since is piecewise-linear (Eq 5), calculating ρ via Eq 15 is straightforward: ρ = 1 − t₀/t_*, where t₀ is once again determined by either or according to Eqs 3 and 4. For a given ρ and cost 〈a〉/x_* or 〈r〉/x_*, we calculate the minimum timing variance . For the activator, we use the expression for ρ along with Eq 13 to write Eq 7 in terms of ρ and 〈a〉/x_*, (16) For the repressor, we use the expression for ρ along with Eq 14 to write Eq 8 in terms of ρ and 〈r〉/x_*, and then minimize over N (see Materials and methods) to obtain (17) Eqs 16 and 17 are shown in Fig 3B (green solid and red dashed curves, respectively), and we see the same qualitative features for both cases: all curves satisfy at ρ = 1, as expected; and as ρ decreases, each curve exhibits a minimum whose depth and location depend on cost. Specifically, as cost increases (lighter shades of green or red), the variance decreases, as expected. Importantly, we see that at a cost on the order of 〈a〉/x_* = 〈r〉/x_* ∼ 10, the model becomes consistent with the experimental data: both the low timing variance and the low degree of linearity predicted by either the activator or repressor case agree quantitatively with the experiment. This suggests that either an accumulating activator or a degrading repressor is sufficient to account for the temporal precision observed in mig-1-controlled neuroblast migration.

Results are robust to additional complexities including cell division

Our minimal model neglects common features of gene expression such as bursts in molecule production [22] and additional sources of noise. Therefore we test the robustness of our findings to these effects here. First, we test the robustness of the results to the presence of bursts by replacing the Poisson process governing the activator production with a bursty production process. Specifically, we assume that each production event increases the activator molecule count by an integer in [1, ∞) drawn from a geometric distribution with mean b [23, 24]. The limiting case of b = 1 recovers the original Poisson process. The results are shown in Fig 4A for b = 1, 3, and 5 (green solid, cyan dashed, and magenta dashed curves). We see that bursts in the activator increase the timing variance of the target molecule, as expected, but that there remain parameters for which the variance is less than that for the unregulated case, (dashed black line). This result shows that even with bursts, regulation by an accumulating activator is beneficial for timing precision.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Results are robust to additional complexities including cell division.

(A, B) Green solid curves show slices from Fig 2 with K = 10 while black dashed line shows unregulated limit . We see that regulation can reduce timing variance even with bursts in activator production of mean size b (A, cyan and magenta dashed), initial Poisson noise in repressor number (B, green dashed), or steady state k/μ in regulator dynamics (blue) unless it approaches regulation threshold K (red). (C) Mean dynamics of activator model (solid) and repressor model (dashed) in which cell division occurs at time on average. Abrupt reductions in molecule numbers are smoothed by noise in t_d and by binomial partitioning of molecules. (D) Timing variance approaches that with no division (dashed) within experimental division region (gray). In A and B, parameters are as in Fig 2. In C and D, parameters are x* = 15, 〈a〉/x_* = 〈r〉/x_* = 10, and H = 3, with kt_*, μt_*, and K set to optimal values (Fig 2) and and σ_d set to experimental values. In all cases, α is set to ensure that mean threshold crossing time equals t_*.

https://doi.org/10.1371/journal.pcbi.1006201.g004

We also recognize that whereas the activator can be assumed to start with exactly zero molecules, it is more realistic for the repressor to start with an initial number of molecules that has its own variability. We incorporate this additional variability into the model by performing stochastic simulations [25] of the reactions in Fig 1A and drawing the initial repressor molecule number from a Poisson distribution across simulations. The result is shown by the green dashed curve in Fig 4B. We see that the additional variability gives rise to an increase in the timing variance of the target molecule, as expected (compare with the green solid curve). However, for most of the range of degradation rates, including the optimal degradation rate, the variance remains less than that of the unregulated case, (dashed black line). This result indicates that the benefit of repression is robust to this additional source of noise.

Then, we test the robustness of the results to our assumptions that the activator undergoes pure production and the repressor undergoes pure degradation. Specifically, we introduce a degradation rate μ for the activator, and a production rate k for the repressor, such that either regulator reaches a steady state of k/μ. The blue curves in Fig 4A and 4B show the case where the increasing activator’s steady state k/μ is twice its regulation threshold K, or the decreasing repressor’s steady state k/μ is half its regulation threshold K, respectively. In both cases, we see that the timing variance of the target molecule increases because the regulator slows down on the approach to its regulation threshold. Nonetheless, we see that it is still possible for the variance to be lower than that of the unregulated case. The red curves show the case where the regulator’s steady state is equal to its regulation threshold. Here we are approaching the regime in which threshold crossing is an exponentially rare event. As a result, the variance further increases, to the point where it is above that of the unregulated case for the full range of parameters shown. These results demonstrate that the benefit of regulation is robust to more complex regulator dynamics, but only if the regulator still crosses its regulation threshold at a reasonable mean velocity.

Finally, we test the robustness of the results to a feature exhibited by the experimental mig-1 data: near the end of migration, cell division occurs (Fig 3A, black data). One daughter cell continues migrating (Fig 3A, dark blue data), while the other undergoes programmed cell death [5]. To investigate the effects of cell division, we perform stochastic simulations, and at a given time t_d we assume that the cell volume V is reduced by a factor of two. For each simulation, t_d is drawn from a Gaussian distribution with mean and variance determined by the data (Fig 3A, black). At t_d, we reduce the molecule numbers of both the regulator and the target molecule assuming symmetric partitioning, such that the molecule number after division is drawn from a binomial distribution with total number of trials equal to the molecule number before division and success probability equal to one half. We also reduce the molecule number threshold K by a factor of two because it is proportional to the cell volume via K = K_dV, where K_d is the dissociation constant.

Fig 4C shows the dynamics of the mean molecule numbers of the activator (green solid) and its target (blue solid), or the repressor (red dashed) and its target (blue dashed). We see that the activator, repressor, and target drop in molecule number at division but that the abruptness of the drop is smoothed by the variability in the division time. The smoothing is more pronounced in the cases of the repressor and the target because the molecule numbers of these species are smaller at division. Thus, for either the activator or repressor mechanism, we see that the experimentally observed variability in division time is sufficient to smooth out the dynamics of the target molecule number, consistent with the experimental data in Fig 3A.

Additionally, we see in Fig 4D that the timing variance of the target molecule in both the activator and repressor cases is similar to that without division in the region of the experimental division time. This further indicates that either model remains sufficient to account for the low experimental timing variance, even with the experimentally observed cell division. Taken together, the results in Fig 4C and 4D show that the key results of the model are robust to the effects of cell division.

Computation of the first-passage time statistics

We compute the first-passage time statistics and numerically from the master equation following [12], generalized to a two-species system. Specifically, the probability F(t) that the molecule number crosses the threshold x_* at time t is equal to the probability P_{y, x*−1}(t) that there are y regulator molecules (where y ∈ {a, r}) and x_* − 1 target molecules, and that a production reaction occurs with rate f_±(y) to bring the target molecule number up to x_*. Since this event can occur for any regulator molecule number y, we sum over all y, (18) where Y = {a_max, N}. The repressor has a maximum number of molecules N, whereas the activator number is unbounded, and therefore we introduce the numerical cutoff . The probability P_yx evolves in time according to the master equation corresponding to the reactions in Fig 1A, (19a) (19b) The moments of Eq 18 are (20) where and . Therefore computing and requires solving for the dynamics of P_yx.

Because Eq 19 is linear in P_yx, it is straightforward to solve by matrix inversion. We rewrite P_yx as a vector by concatenating its columns, , such that Eq 19 becomes , where (21) Here, for i, j ∈ {0, …, Y}, the x_* − 1 subdiagonal blocks are the diagonal matrix , and the x_* diagonal blocks are the subdiagonal matrix or the superdiagonal matrix for the activator or repressor case, respectively. The term prevents activator production beyond a_max molecules. The final M⁽¹⁾ matrix in Eq 21 contains f_± production terms that are not balanced by equal and opposite terms anywhere in their columns. These terms correspond to the transition from x_* − 1 to x_* target molecules, for which there is no reverse transition. Therefore, the state with x_* target molecules (and any number of regulator molecules) is an absorbing state that is outside the state space of [12]. Consequently, probability leaks over time, and . Equivalently, the eigenvalues of M are negative.

The solution of the dynamics above Eq 21 is . Therefore, Eq 20 becomes , where is a row vector of length x_*(Y + 1) consisting of [f_±(0), …, f_±(Y)] preceded by zeros. We solve this equation via integration by parts [12], noting that the boundary terms vanish because the eigenvalues of M are negative, to obtain (22) We see that computing and requires inverting M, which we do numerically in Matlab. We initialize as P_ax(0) = δ_a0 δ_x0 or P_rx(0) = δ_rN δ_x0 for the activator or repressor case, respectively.

When including cell division, we compute and from 50,000 stochastic simulations [25].

Deterministic dynamics

The dynamics of the mean regulator and target molecule numbers are obtained by calculating the first moments of Eq 19, and , where y ∈ {a, r}. For the regulator we obtain or in the activator or repressor case, respectively, which are solved by Eqs 3 and 4. For the target molecule we obtain , which is not solvable because f_± is nonlinear (i.e., the moments do not close). A deterministic analysis conventionally assumes , for which the equation becomes solvable by separation of variables. For example, in the case of H = 1, using Eqs 1–4, we obtain (23) Eq 23 is plotted in Fig 1C and 1D.

When including cell division, we compute the mean dynamics from the simulation trajectories (Fig 4C).

Details of the analytic approximations

To find the global minimum of Eq 8, we differentiate with respect to kt_* and K and set the results to zero, giving two equations. kt_* can be eliminated, leaving one equation for K, (24) This equation is transcendental. However, in the limit K ≪ N, we neglect the last term, which gives Eq 11.

To derive Eq 17, we use (25) where the second step follows from according to Eq 4; and, from Eq 14, (26) where the second step assumes that the repressor is fast-decaying, μt_* ≫ 1. We use Eqs 26 and 25 to eliminate μt_* and K from Eq 8 in favor of ρ and 〈r〉, (27) For nonlinear dynamics (ρ < 1) we may safely neglect the −1 in Eq 27. Then, differentiating Eq 27 with respect to N and setting the result to zero, we obtain N = 3〈r〉/(1 − ρ). Inserting this result into Eq 27 produces Eq 17.

Analysis of the experimental data

To estimate the time at which migration terminates in Fig 3A, we refer to [5]. There, the position at which neuroblast migration terminates is measured with respect to seam cells V1 to V6 in the larva (see Fig. 4D in [5]). In particular, in wild type larvae, migration terminates between V2 and the midpoint of V2 and V1. This range corresponds to the magenta region in Fig 3A (see Fig. 4B, upper left panel, in [5]). Under the assumptions of constant migration speed and equal distance between seam cells, the horizontal axis in Fig 3A represents time.

To compute ρ for the experimental data in Fig 3A according to Eq 15 we use a trapezoidal sum. For the choices of x_* and t_* described in the text, this produces the ρ values in Fig 3B.