Mean-Independent Noise Control of Cell Fates via Intermediate States.

Stochasticity affects accurate signal detection and robust generation of correct cell fates. Although many known regulatory mechanisms may reduce fluctuations in signals, most simultaneously influence their mean dynamics, leading to unfaithful cell fates. Through analysis and computation, we demonstrate that a reversible signaling mechanism acting through intermediate states can reduce noise while maintaining the mean. This mean-independent noise control (MINC) mechanism is investigated in the context of an intracellular binding protein that regulates retinoic acid (RA) signaling during zebrafish hindbrain development. By comparing our models with experimental data, we find that the MINC mechanism allows for sharp boundaries of gene expression without sacrificing boundary accuracy. In addition, this MINC mechanism can modulate noise to levels that we show are beneficial to spatial patterning through noise-induced cell fate switching. These results reveal a design principle that may be important for noise regulation in many systems that control cell fate determination.

Introduction

Stochasticity is prevalent in cell signaling networks, yet organisms manage to determine cell fates and coordinate their development in response to signals robustly in spite of this noise. Stochasticity has been directly observed in gene expression (Elowitz et al., 2002, Singh et al., 2012) and identified as a cause of cell-to-cell variation (Marinov et al., 2014, Raj and van Oudenaarden, 2008). The amount of stochasticity in gene expression can change the qualitative features of the system, such as introducing bistability in a deterministically monostable system, and vice versa (Kepler and Elston, 2001), or allowing switching between stable states (Paul, 2017, Koseska et al., 2009, Wu et al., 2013, Wu et al., 2017). Noise can have adverse effects by distorting downstream signals (Filippi et al., 2016) and by disrupting entrainment of biochemical oscillators (Gupta et al., 2016), suggesting the existence of noise control mechanisms to ensure robustness. However, these mechanisms are not well understood.

There is increasing evidence that such systems of noise control exist in signaling pathways. Differences in regulatory network architecture correlate with both the amplification and attenuation of expression noise (Chalancon et al., 2012, Paulsson, 2004), and increased complexity of a biochemical network has been shown to correlate with reduced noise (Cardelli et al., 2016). For one autoregulatory protein it has been suggested that negative feedback decreases system noise (Thattai and van Oudenaarden, 2001). In addition, increased growth rates in a single-celled organism (yeast) can increase noise in gene expression (Keren et al., 2015), and studies in Drosophila suggest that these principles apply to noise in spatial signals in a multicellular (albeit syncytial) context (Gregor et al., 2007). Multiple binding sites for the Bicoid morphogen in the Hunchback promoter buffer spatial noise in the Bicoid morphogen gradient (Holloway et al., 2011). In addition, specific noise levels are advantageous for establishing population heterogeneity (Kaern et al., 2005), and noise can facilitate sharp segmental boundaries of gene expression in the zebrafish hindbrain, but only within a limited range of levels of signaling noise (Zhang et al., 2012). It has become increasingly clear not only that noise needs to be attenuated but also that appropriate levels of noise are critical for accurate cellular responses to a signal.

One example of noise regulation occurs in the developing vertebrate hindbrain, which is patterned by a retinoic acid (RA) morphogen gradient. RA is a well-known signaling molecule (Niederreither and Dolle, 2008) that controls the formation of hindbrain segments (Sirbu et al., 2005), as well as the patterning and differentiation of many other cell types and tissues. Models of morphogen gradients suggest that they can cause all-or-none regulatory responses (Meinhardt, 2009, Wartlick et al., 2009), even in the presence of stochasticity, and consistent with this, RA signaling specifies segmental patterns of hindbrain gene expression despite substantial levels of noise (Schilling et al., 2012). We used fluorescence lifetime imaging microscopy (FLIM) to measure stochasticity in the RA gradient and directly demonstrated the existence of large fluctuations in the gradient (Sosnik et al., 2016). One of the four cellular RA-binding proteins (Crabp2a) in zebrafish is essential for robust patterning in the developing hindbrain (Cai et al., 2012). Experiments increasing or decreasing Crabp2a levels show that it acts to attenuate noise in RA while leaving the mean unchanged (Sosnik et al., 2016), pointing to the existence of a mechanism that allows for the desired noise levels to be achieved without shifting the gradient.

Using the RA signaling network in zebrafish hindbrain development as an example, we investigate design principles for controlling noise in the signal without affecting the mean levels. We find that a coupling of reversible reactions gives rise to mean-independent noise control (MINC), and this coupling naturally arises in complex systems through the presence of an intermediate state. In the example of spatial patterning of hindbrain segments, the MINC mechanism involving Crabp2a enables mean-independent sharpening of gene expression boundaries in the correct locations in response to RA. In addition, we find that the degree to which the mean and variance are coupled distinguishes between different noise origins, and when analyzed with the FLIM data obtained from zebrafish embryos, our results suggest that the dominant noise source is endogenous to the RA signaling pathway.

Results

Proportional-Reversibility Enables Mean-Independent Noise Control

First, we considered a simplification of the differential equation model (Schilling et al., 2012) to capture the essential qualitative features. In the RA signaling pathway, extracellular RA enters the cell and binds to an intermediate (CRABP), which shuttles it to the nucleus to bind to a receptor (RAR) and form a compound (RA-RAR), which binds to the DNA to both signal downstream targets and produce a protein (CYP), which in turn inactivates RA (depicted in Figure 1A). The basic interaction (denoted as the Simple Model [SM]) can be modeled as a two-state stochastic differential equation (SDE):where the deterministic portion of the equation is due to mass action laws and the additional term (σdW) describes stochasticity in the production and degradation of RA (schematized in Figure 1B). The steady-state mean values areand the steady-state variances are(see the Transparent Methods for details of the derivation). Assume that the rates for the reversible binding and unbinding of the morphogen to its receptor are coupled via a constant C:

Figure 1

The Proportional-Reversibility Strategy for Noise Attenuation

(A) A diagram of the retinoic acid (RA) signaling network. RA (shown in red) enters the cell to bind with the cellular retinoic acid-binding proteins (CRABP, shown in green), which shuttles RA into the nucleus. There the RA binds with its receptor (RAR, shown in purple) to produce Cyp26 (shown in pink). The Cyp26 proteins in turn deactivate the RA signaling proteins.

(B) A schematic depiction of the two-state simplified model (SM) of the RA network. The red node in the graph is the concentration of RA, and the blue node is RA bound to its receptor (RA-RAR). Birth and death of RA are allowed, along with reversible binding/unbinding with RAR.

(C) Representative time-series solutions to SM. These were obtained using the Euler-Maruyama method. (I) Shows the solution using parameters from Table S1 while (II) uses the same parameters with γ and δ decreased by a factor of 0.005. Both trajectories start from the expected value of 10. Notice that in (I) the blue line is mostly covered by the orange line.

(D) The covariances from experiments (I) and (II).

(E) A schematic depiction of the full RA model RM. The red node depicts the concentration of RA, which binds with the binding protein (BP, shown in teal) to form the RA-BP intermediate complex (shown in green). From there, the RA binds to its receptor RAR (shown in pink) to form RA-RAR (shown in purple). This causes the transcription of Cyp26a1, which degrades RA.

(F–H) Representative time-series solutions to RM. The parameters for the model were chosen randomly according to the method described in the Transparent Methods. The resulting stochastic differential equations were solved using the Euler-Maruyama method for 100 s to give the blue line for the wild-type value (F). The value γ was reduced by 90% and the simulation was re-solved to give the orange line (G). The value for γ was reset and the value for α was reduced by 90%, and the simulation was re-solved to give the green line (H).

Under this coupling assumption, the steady-state values become proportional,with the mean concentrations determined solely by the production rate β and the decay rate η of RA. However, under the same conditions,the variance directly depends on γ, the rate of conversion from RA to RA-RAR. Note that increasing γ attenuates noise in the RA concentration (derived in the Transparent Methods). Thus, with this coupling assumption, changing the reaction rates of RA binding and dissociation from its receptor has no effect on steady-state mean concentrations but has a direct and directional effect on their variances. Comparing the temporal dynamics of the system to a setup with reduced binding rates suggests that this mean-independent variance effect is due to changes in the binding rates (Figure 1C).

The mean amounts of the proteins are left unchanged by the coupled reaction rates since the amount of RA that leaves its unbound state increases by the same amount that leaves its bound state due to the coupling assumption. To intuit why the variance decreases, we computewhich is an increasing function in γ. Thus when the binding and unbinding rates are higher, the concentrations of RA and RA-RAR tend to be in sync (Figure 1D). Therefore if [RA] is above its steady state and the binding/unbinding rates are high, then it is highly likely that [RA−RAR] will also be above its steady-state levels. One can heuristically understand that when random fluctuations cause an excess of morphogen, the natural force toward steady state will drive the excess morphogen toward degradation. If [RA] is below its steady-state levels but [RA−RAR] is even lower, the system could push some of the morphogen to its receptor-bound state even if the total morphogen in the system is lower than its equilibrium value, thus decreasing the total pull toward steady state. Therefore the total pull toward equilibrium is greatest when the correlation is highest, which explains why the variance decreases and saturates to a constant as leads to near-perfect correlation. This intuition suggests that this feature is a property of the coupling of the reversibility. To see if this extends to more complex systems, we note that [RA−RAR] enhances the production of a protein Cyp26a1, which in turn degrades [RA] (White et al., 2007, White and Schilling, 2008). We show that when incorporating this nonlinearity into the previous model, the essential feature of mean-independence from the reversible reaction rates holds under the same coupling assumptions (see the Transparent Methods). In addition, the mean-independence, along with the noise attenuation and increase in covariance due to γ, holds for the general master equation formulation of the model (see the Transparent Methods).

The previous results rely on the assumption that the binding and unbinding rates are proportional. However, next we asked if this coupling naturally occurs in the presence of an intermediate state. The RA signaling network includes an intermediate state RA-BP where RA is bound to its binding protein Crabp2a (BP), which shuttles RA to the nucleus where it binds to its receptor to form RA-RAR (Cai et al., 2012). Adding this interaction to the previous model gives the Intermediate Model (IM). Notice that in this model the flux into [RA−RAR] at steady state is ν[RA−BP], whereas the internal flux back to [RA] is δ[RA−BP]. Therefore near steady state, [RA−BP] is a natural coupling constant between the influx of RA and the influx of RA-RAR, suggesting a generalization of the proportional-reversibility mechanism. The heuristic derivation is confirmed since one can show that a proportional coupling of the rate from [RA] to [RA−BP] and the rate from [RA−RAR] to [RA−BP] produces the same mean-independence and variance-dependence on the coupled reaction rates behavior as in the previous model and that E[RA−BP] is proportional to the proportionality constant (see the Transparent Methods for details). In addition, we can extend the signaling pathway by adding an intermediate unbinding in the nucleus step or model Crabp2a's shuttling of RA to Cyp26a1 as a separate pool of RA, and in this extended model the MINC property still holds (see the Transparent Methods).

We build upon this initial analysis to generate a more mechanistic model that fully incorporates the BP and receptor concentrations, which we refer to as the Retinoic Acid Model (RM) (Figure 1E). In this model, we derive the same natural coupling via E[RA−BP] and also show that E[RA−BP] is proportional to E[BP], which in turn are directly determined by the production and degradation rates of BP (see the Transparent Methods section for details of the derivation). Therefore the production and degradation rates of BP have no effect on the mean amounts of [RA] and [RA−RAR], whereas they directly affect the variance. Comparisons of [RA], [RA−BP], and [RA−RAR], respectively, for two separate trajectories of the system show that the means of [RA] and [RAR] are unchanged by perturbations in [BP] but that the variance changes (Figures 1F–1H). Perturbations in [Cyp] change the mean as well as the variance, but the latter to a much lesser extent.

Last, we note the experimental evidence that RA-RAR upregulates the production of Crabp2a but no other Crabps (Cai et al., 2012, Schilling et al., 2012). To simulate the effect of RA-RAR signaling on Crabp production, along with other possibly indirect downstream effects such as Eph and Ephrins on the RA signaling pathway (Wang et al., 2017), we introduce the Retinoic Acid Model with Feedback (RMF) by adding upregulation of BP by [RA−RAR] to RM. We derive that the mean-independent variance control via a coupling related to E[BP] still holds in the RMF model. This shows that by incorporating the dynamics of an intermediate state in the model we find a natural and controllable coupling of mean reaction rates that allows for the amount of intermediate to attenuate the system noise without changing the mean levels of the signal.

Active Intermediate States Naturally Lead to Proportional-Reversibility Control

The previous results were determined by linearizations of stochastic models around steady-state values. To establish the mean-independent noise attenuation property directly on nonlinear complex models, we performed simulations. First, we introduced a mean-variance dependence index (MVDI),to encapsulate the relation of mean and variance changes into a single non-dimensionalized value. If the variance changes but the mean does not, ζ is exactly 1. If the mean changes but the variance does not, then ζ is 0. Thus the MVDI ζ offers a measure of the degree of isolation of the mean relative to variance changes with respect to a perturbation. Four separate models of the RA signaling pathway (SM, IM, RM, and RMF) were simulated to determine the robustness of the mean-independent noise attenuation property to the network topology for RA alone (Figures 2A–2D) or bound to RAR (Figures 2E–2H). Each model was simulated using random reaction rates 100,000 times and ζ′s calculated from reductions in BP and Cyp (Table 1). Over 85% of the simulations resulted in a ζ > 0.9 with a decrease in BP, indicating that for most parameter sets the relative change in variance was much larger than the relative change in mean. In addition, when the same computational experiment was applied with multiplicative noise, over 90% of the simulations resulted in a ζ > 0.9 with BP knockdown (Figure S1). We note that from the experimental data (Sosnik et al., 2016) we estimate that the true ζ after depleting BP is ζ = 0.97. This shows that in almost all the trajectories in each model, the variance changes with only minor changes in the mean, confirming that proportional-reversibility and intermediate states naturally give rise to MINC. In addition, over 90% of simulations result in ζ < 0.6 with reductions of Cyp under both additive and multiplicative noise, showing that this property is specific to manipulation of BP levels.

Figure 2

Mean and Variance Knockdown Distributions

(A–H) Histograms depicting the ζ distribution due to additive noise. The models were solved using the Euler-Maruyama method on the time span of t = [0, 200] and re-solved with a 90% knockdown to the associated parameter. The value ζ was calculated from these two series, and this was repeated 100,000 times. More details outlining the experiments are found in the Transparent Methods. (A–D) ζ′s Calculated for [RA] on the respective models. (E–H) ζ′s Calculated for [RA−RAR] on the respective models.

Table 1

Numerical Knockdown Experiment Base Exponents

Parameter	baseExpp
β	2
a	0
b	1
α	2
u	1
ζ	0
η	1
r	1
All others	3

Each parameter was chosen by taking x uniformly from [-2,2] and calculating . On the left is the parameter and on the right is the corresponding baseExp value. All values not listed in the table had the value baseExp = 3.

Spatial MINC Enables Accurate Specification of Landmark Locations

We hypothesize that the mean-independent variance attenuation property of BP should locally smooth spatial gradients. To study the MINC property in spatial signaling, we first extended the RMF to a two-dimensional stochastic reaction-diffusion system with space-time white noise where the extracellular RA, RA, diffuses throughout space. This model, RMFS, is thus given by a stochastic partial differential equation (SPDE) (defined in the Transparent Methods). As counterparts to the mean and variance in space, we measured the mean and variance of the location where the gradient hits specific values. To determine if the same mean/variance relationship is associated with these properties of the morphogen gradient, we simulated BP depletion experiments on the RMFS model with random parameters (Figure S2 and detailed in the Transparent Methods). Gradients of [RA] and [RA−RAR] exhibit less noise in wild-type than after BP depletion (Figures 3A–3D). Given that the average size of cells in the zebrafish hindbrain is around 10 μm, our simulations suggest that changes in the amount of BP change the sharpness of the RA signaling boundary without changing its mean location by more than one cell diameter over 90% of the tested parameters (Figures 3E and 3F). In addition, the mean shift in the RA-RAR gradient with decreased Cyp shifts the threshold approximately 4 cell diameters, showing that spatial MINC is a special feature that does not extend to other interactions (Figure S3B).

Figure 3

Location-Independent Boundary Sharpening

(A–D) Representative solution of RMFS. Shown are the two-dimensional gradients of intracellular [RA] (plots A and C) and [RA−RAR] (plots B and D). Random parameters were chosen according to the Transparent Methods. The x axis runs from the anterior to the posterior of the zebrafish hindbrain. The simulations were run to determine the steady-state gradients and then were run for 100 s more to give a snapshot of the spatial stochasticity. (A and B) The wild-type results for a given parameter set, and (C and D) the results when simulated again but with the parameter for binding protein (BP) reduced by 90%. The color levels are fixed between the wild-type and BP-deficient plots to provide accurate comparisons.

(E) Scatterplot of the percent change in variance versus the change in the average boundary location (μm). 100 Parameter sets were chosen, and for each one the threshold concentration was taken to be 60% of the maximum value. The mean and the variance of the boundary location was calculated, the simulation was solved once more using the BP-deficient value, and the mean and variance of the boundary location were calculated using the same threshold.

(F) Histogram of the number of simulations with a given change in mean boundary location (μm).

The receptor-bound RA signal induces the expression of Hoxb1a and Krox20, which mutually repress each other (Zhang et al., 2012), which we include downstream in the RMFS model (depicted in Figure S3A). By preferentially upregulating Krox20, the boundary between rhombomeres 4 and 5 (r4/5) is established at a threshold where the RA gradient is sufficiently high enough for the initial Hoxb1a expression to be replaced by Krox20. Stochasticity in the Hox-Krox interaction allows for the initial r4/5 boundary to sharpen in the wild-type organism (Figures 4A–4D). In addition, reducing BP disrupts the sharpening of the r4/5 boundary without moving its location, whereas reducing Cyp causes a shift in the segmental boundary position. Over time the number of predominantly Hox-expressing cells displaced in the wild-type saturates to approximately 3, whereas with depletion of BP the number of displaced cells saturates to approximately 10 (Figure 4E). On average a stochastic trajectory of the wild-type has a maximal displacement of predominantly Hox-expressing cells by 1–2 cell diameters, whereas with reduced BP, there is a maximal displacement of predominantly Hox-expressing cells of around 4 cell diameters. Segmental sharpening in terms of a previously defined sharpening index (SI) from Zhang et al. (2012) reveals a similar disruption of the sharpening mechanism under this measurement (Figure S3C). Together, these results show that the loss of BP disrupts downstream segmental sharpening, consistent with the previous in vivo experimental observations in zebrafish (Sosnik et al., 2016).

Figure 4

Disruption of Downstream Boundary Sharpening Due to Perturbations in Binding Protein (BP) and Cyp26 Concentration

(A–D) Representations of the rhombomere 4/5 segmental boundary. Representation at (A) 10 hr post-fertilization (hpf), (B) 10.5 hpf, (C) 11 hpf, and (D) 12 hpf. The trajectories were calculated according to the RMFS model with the Hox-Krox extension (described in the Transparent Methods). The color corresponds to the concentration of Hox at a particular region in the anterior-posterior versus the left-right plane. The top panels correspond to the wild-type model, the middle panels are the models with reduced BP, and the bottom panels are the models with reduced Cyp26.

(E) Number of displaced cells over time. The y axis corresponds to the number of displaced cells, calculated as a predominantly Hox-expressing cell lying one cell length posterior to a predominantly Krox-expressing cell. The x axis shows time in terms of hpf. Each condition was run 10 times, and the results were averaged.

(F) Maximum displacement over time. The y axis corresponds to the maximum displacement, calculated as the maximum distance between a predominantly Hox-expressing cell that is posterior to a predominantly Krox-expressing cell. The axis is in cell diameters, which correspond to 10 μm. The x axis shows time in terms of hpf. Each condition was run 10 times, and the results were averaged.

Noise Levels are Regulators for Patterning and Indicative of the Sources of Stochasticity

To further investigate the significance of noise levels in RA signaling, we next investigated the requirements for an optimal range of noise levels in gene expressions and we explored the relationship between the noise level and the origins of the stochasticity. To establish a direct relationship between the noise levels controlled via BP and successful segmentation of the zebrafish hindbrain, we defined the effective noise in the [RA−RAR] signal as a measure of the variance of the gradient (described in the Transparent Methods). We simulated the full spatial model RMFS with varied [BP] and Hox-Krox signaling noise to determine if noise levels affected the ability to sharpen the r4/5 boundary. A successful sharpening event was defined as having less than or equal to 3 cells displaced by more than 1 cell diameter. All the successful sharpening events had Hox-Krox noise levels in the range 0.175–0.275 with an upper limit on the effective [RA−RAR] signaling noise of approximately 10−3 (Figure 5A). Similar qualitative results were obtained when successful sharpening events were defined instead in terms of a threshold on mean displacement and maximal displacement (Figure S4). This shows that an optimal range for the effective noise is required for segmental patterning to occur, indicating the necessity of noise control mechanisms for properly regulating downstream signals.

Figure 5

Noise Levels Distinguish Between Models and Developmental Phenotypes

(A) Scatterplot of the successful and unsuccessful sharpening events. BP production was taken as 5, 15, …, 105 and Hox-Krox regulatory noise was taken as 0.1, 0.125, …, 0.3 and each pairing in the grid was solved once. The effective RA noise was calculated according to the measure from the Transparent Methods, and a sharpening event was declared successful if after 2 hr less than 3 cells were displaced by more than 1 cell diameter.

(B) Depiction of the probability distribution for ζ according to the parameter search scheme from the parameter search scheme on RMF. The simulations which produced these distributions are discussed in the Transparent Methods. Shown are the kernel density estimates from the ζ values from the stochastic simulations. The different colored lines show the distributions for different noise types. The red circles depict experimental values for ζ computed pairwise.

ζ Calculations for changes in BP do not distinguish between common noise types such as multiplicative or additive noise, but the probability distribution of ζ with respect to changes in the amount of Cyp strongly depends on the choice of the noise term (compare Figure 2 with Figure S1). Therefore we investigated the validity of particular noise terms by comparing these probability distributions with experimental data. We examined data collected from free intracellular RA in the zebrafish hindbrain morphogen gradients with a morpholino knockdown of Cyp26a1 (Sosnik et al., 2016). The methodology for the analysis (see Transparent Methods) estimates the average experimental value as ζ ≈ 0.62.

To determine the likelihood that this system uses one type of noise versus another, we utilized the same parameter search scheme as previously described to calculate probability distribution for ζ with different noise types and depletion of Cyp on the model RMF (Figure 5B). An experimental ζ was calculated pairwise between each wild-type and Cyp-deficient embryo since there does not exist a canonical pairing, thus giving a distribution of 27 experimental ζ values, which center around ζ = 0.6 with a tail toward zero. ζ Was less than the mean experimental ζ for >90% of the parameters when the noise was additive for [RA] or multiplicative for [RA−RAR]. In contrast, the ζ distributions with multiplicative noise for [RA], [RA−BP], and [RA] all peak around the mean ζ value. Similar results were obtained for the cumulative distribution of ζ (Figure S5). This reliance of the noise on the levels of [RA]in itself (all these concentrations are directly linked) strongly suggests that the dominant form of noise in the zebrafish hindbrain signaling network is intrinsic to stochastic processes related to [RA]in and establishes that exogenous noise is not likely to fit the data.

Discussion

Stochasticity in gene regulatory networks (GRNs) naturally exists (Elowitz et al., 2002), and noise attenuation or control is needed to enable proper biological functions. For example, in the zebrafish hindbrain, noise regulation is required for subsequent boundary sharpening processes to occur properly (Zhang et al., 2012). Here we uncover an MINC mechanism that can tune the level of noise in the downstream components of a GRN, without affecting the mean of the signal. In the zebrafish hindbrain system, this MINC mechanism provides a way (through Crabp2a) to achieve the required noise levels in the RA morphogen without disturbing other aspects of its spatial gradient. Together, we directly link the preservation of a stochastic spatial phenotype to a noise control mechanism, demonstrating a potential path through which developmental processes could have evolved to overcome inherent biochemical stochasticity to achieve robust spatial patterning.

The robustness in the choice of biologically reasonable parameters needed to achieve such a mechanism indicates that it is an intrinsic property of the GRN topology. The core feature underlying this is a coupling assumption that arises naturally with the existence of intermediate states under mass action assumptions. Furthermore, we obtain similar results with five separate models, all with this same coupling. This suggests that MINC may be a general phenomenon related to the existence of intermediate states and probably exists in other GRNs. Computational simulations of epithelial-mesenchymal transitions (EMTs) recently has shown that increased numbers of intermediate states attenuate noise in cellular fate decisions (Ta et al., 2016), which is analogous to our predictions of how pooling in the intermediate state reduces noise. Primed lineages in hematopoietic stem cells can be represented as an intermediate state with reversible changes, which could explain data showing mean-independent noise attenuation due to a lineage commitment factor (van Galen et al., 2014). However, we note that the spatial control results were only tested in the areas of the RA gradient where the concentration is sufficiently high. There may be other factors required for robust noise control and boundary sharpening in anterior rhombomeres (r1-3) where the RA concentration is low and the gradient is shallow. In addition, our models do not take into account the effects of cell proliferation. While the cell cycle rapidly increases to 4 hr around the time of boundary sharpening (Kimmel et al., 1994), which is a much slower timescale than the noise processes, further research could better quantify the effects of cell divisions on the spatial noise.

Furthermore, our methods uncover a novel relationship between the noise source, the network topology, and the relationship between the mean and the variance using the perturbation data in experiments. Our analysis suggests intrinsic noise due to RA as the most likely dominant noise source in the zebrafish hindbrain given the known GRN topology and mean-variance relationship. With the increasing precision in experimental quantification of variance changes, this methodology could be used to identify noise sources and provide further evidence for or against GRN topologies. For example, microfluidic measurements have accurately measured noise dynamics in individual aging yeast cells (Liu et al., 2017) and a dynamic analysis similar to the one shown here could restrain the possible GRNs by requiring that not only the mean but also the noise dynamics match the data. Most importantly, this approach may be used to distinguish between models that have similar qualitative behavior with respect to the mean, thereby providing a new way to uncover details of biochemical networks from the noise in gene knockdown experiments.

Methods

All methods can be found in the accompanying Transparent Methods supplemental file.