Abstract

Recent measurements of morphogen distribution have allowed us to subject this hypothesis to rigorous physical testing. In the early Drosophila embryo, measurements of the morphogen Dorsal, which is a transcription factor responsible for initiating the earliest zygotic patterns along the dorsal-ventral axis, have revealed a gradient that is too narrow to pattern the entire axis. In this study, we use a mathematical model of Dorsal dynamics, fit to experimental data, to determine the ability of the Dorsal gradient to regulate gene expression across the entire dor-sal-ventral axis. We found that two assumptions are required for the model to match experimental data in both Dorsal distribution and gene expression patterns. First, we assume that Cactus, an inhibitor that binds to Dorsal and prevents it from entering the nuclei, must itself be present in the nuclei. And second, we assume that fluorescence measurements of Dorsal reflect both free Dorsal and Cactus-bound Dorsal. Our model explains the dynamic behavior of the Dorsal gradient at lateral and dorsal positions of the embryo, the ability of Dorsal to regulate gene expression across the entire dorsal-ventral axis, and the robustness of gene expression to stochastic effects. Our results have a general implication for interpreting fluorescence-based measurements of signaling molecules.

Author Summary

Across the fruit fly embryo, proteins known as morphogens cause chain reactions that eventually lead to the growth of many different body parts. One specific morphogen, known as Dorsal, is distributed in a gradient across the dorsal-ventral (DV) aXis of the embryo. The Dorsal gradient has recently been measured in live embryos using a fluorescently-tagged version of Dorsal, but the gradient measured by fluorescence microscopy seems to be too narrow to pattern all of the DV axis. Using a mathematical model of the Drosophila embryo, we have proposed a solution to this outstanding problem: namely that Cactus, the inhibitor to Dorsal, is present with Dorsal in nuclei across the embryo, which creates a disparity between the gradient measured by fluorescence and the gradient measured by gene expression. By using distinct model equations for active and inactive pools of Dorsal, we were able to recreate the dynamics of the Dorsal gradient and the eXpression patterns of its target genes with a high level of accuracy, showing that mathematical models may be critical for properly interpreting fluorescence data.

Introduction

In the early (1—3 hr old) Drosophila embryo, the transcription factor Dorsal (dl; Fly-Base ID: FBgn0260632) acts as a morphogen to pattern the embryo’s dorsal-ventral (DV) axis ([1, 2]; reviewed in [3, 4]). dl, a homologue of mammalian NF-KB, is expressed ubiquitously throughout the syncitial blastoderm, and is sequestered to the cytoplasm by its association with the inhibitor Cactus (Cact; FlyBase: FBgn0000250), a homologue of mammalian IKB [5]. Signaling through Toll receptors on the ventral side of the embryo causes the dissociation of the dl/Cact complex, and free dl accumulates in the ventral nuclei [5—7] to create a spatial gradient that causes differential gene expression based on multiple gene expression thresholds.

Types I, II and 111+ (such as sna, vnd and 50g; FlyBase: FBgn0003448, FBgn0261930, FBgn0003463) are activated by high, moderate, and low levels of nuclear dl, respectively, and thus form boundaries at roughly 20%, 33%, and 50% DV position (with 0% DV being the ventral midline). Type 111- genes (such as dpp and zen; FlyBase: FBgn0000490, FBgn0004053) are repressed by dl, and thus are expressed in nuclei on the dorsal half of the embryo where there is little or no d1 [9, 10].

Quantitative measurements in fixed embryos revealed a dl gradient that becomes flat by ~ 40% DV position, calling into question how such a short-ran-ged gradient could pattern the Type III genes [12]. However, these results were in contrast to classical enhancer trap studies that suggested that dl directly sets the border of the Type III-gene zen [13]. One proposed explanation to rectify this difference is the dl gradient width may be dynamic, as in the case of Tribolium [14].

This “saw tooth” observation was later confirmed by both detailed measurements in fixed embryos [12], as well as modeling work [15]. In contrast, dl levels in the dorsal-most nuclei were found to slowly decrease during interphase, and recover back at the start of the next interphase [10]. However, while the nuclear concentration of dl was found to be highly dynamic all along the DV axis, the seemingly narrow width of the dl gradient was measured to be constant [10, 12].

To determine the answers to these questions, we analyze a model of dl/Cact dynamics. Our model is based on previously-published work [15] , and makes two distinctive assumptions: first, that Cact can be present in the nuclei to regulate dl; and second, that fluorescence measurements (from either dl immunofluorescence experiments or from fluorescent protein-tagged dl) comprise both free and Cact-bound dl. Under these two assumptions, our model, when fit to the spatiotemporal measurements in live embryos [10], can eXplain fine details of the dynamics of the dl gradient. Using this result as a starting point, we simulated dl-dependent gene eXpression, taking into account the fact that only the free (and not Cact-bound) nuclear dl can regulate gene eXpression. We found that nuclei can accurately interpret their position along the DV aXis even in regions of the embryo where stochasticity causes high levels of gradient readout error. Thus, our model eXplains both the dynamics of the dl gradient as well as its spatial range. Our results have implications on how fluorescence measurements should be interpreted, in particular when there is a binding partner that modulates the activity of a signaling molecule.

Methods

dI/Cact model formulation

Here we sketch the essentials for understanding the model. (For full details, see 81 Text.) To simulate the dynamics of dl and Cact during NC10-NC14, a cross section of the embryo was modeled as a linear array of rectangular prism-shaped compartments that each contain a single nucleus (Fig. 1a,b). Each compartment and each nucleus are well mixed, with slow exchange between neighboring compartments [11, 17]. Because the embryo is approximately symmetric about the DV axis, only one half of the axis is simulated and no-exchange boundary conditions are assumed at both the ventral and dorsal midlines. The number of compartments, as well as their dimensions, depends on the number of nuclei, which increases (to the nearest integer) by a factor of

1b). The length and width of each compartment is calculated as the length of the simulated region, L, divided by the number of nuclei in interphase i, 11,-; the height, H, remains constant. The simulation begins at the onset of NC10 interphase with the initial conditions for each molecular species uniform in space. The normalized concentration of nuclear/cytoplasmic Cact is at steady state; the concentration of dl/Cact complex is unity; and the concentration of free dl is zero.

(For a description of the initial model used—depicted in Fig. 2—see 81 Text.) During interphase, the nuclear and cytoplasmic concentrations are governed by the parameters associated with nuclear import/ export, intercompartmental exchange, production/ degradation of Cact, and binding/unbinding of dl/Cact complex. At the start of mitosis the nuclei break down and the contents of each nucleus are mixed with that of the surrounding cytoplasmic compartments (Fig. 1c, center panel). The equations for the mitosis phase are subsequently solved for the appropriate duration, and the next interphase begins with nuclei and cytoplasmic compartments having the same concentration (Fig. 1c, right panel). This allows for the presence of all three molecular species (dl, Cact and dl/Cact complex) to exist in the nuclei at the start of interphase. The surface area and volume of each nucleus, while dynamic over the course of the entire simulation, are considered static for the duration of each NC for simplicity. Nuclear dimensions are calculated using published measurements ([18]; for full description see 81 Text). Interphase and mitosis dynamics are represented by six and three non-dimensionalized differential equations, respectively, consisting of up to 15 free parameters (see Box 1; see 81 Table and 82 Table for dimensional analysis). During interphase, each of the three molecular species is represented by two equations, one for the nuclear concentration and one for the cytoplasmic concentration. During mitosis the nuclei are undefined, so only the equations representing cytoplasmic concentrations are solved. When appropriate, the import/export rates (C 1-, 51-), intercompartmental exchange rates (1,), dl/Cact binding (7/) and dissociation (fig) rates, and degradation of Cact (a) are all represented by free parameters. The import and export rates of each species are approximated as a first-order process, and exchange is formulated as mass transport flux between neighboring compartments. In addition, the Toll-mediated dissociation rate of dl/Cact complex is modeled as a Gaussian curve centered on the ventral midline (fie_%($) ), Which represents the distribution of active Toll receptors. The intensity and range of Toll-mediated dissociation is specified by the amplitude and Width parameters (fl, gb), respectively.

Gene expression model

For simplicity, f is formulated as a hard-threshold on/off sWitch for each species (Hill function With Hill exponent of n H = 100), Where production is unity When the concentration of dl is above (for sna/sog/vnd) or below (for zen) a threshold BdlszNAi. In addition, production of vnd and 50g is repressed by sna above a threshold esmszNAi. To simulate the error in gradient interpretation in each nucleus due to the stochasticity of the arrival of dl at the enhancer site (Which we Will hereafter refer to as signal noise), the effective concentration of dl in nucleus 11 is calculated as Us? 2 Uh —|— nN(0, 1) x / U51“, Where MO, 1) is a random number selected from the standard normal distribution and 17 is a tunable constant. (Note: we assume that random fluctuations in the actual concentration of nuclear dl are small compared to the fluctuations in effective concentration due to the random arrival of dl molecules at enhancer sites.) The value of 17 is chosen to be between 0.02 and 0.5. Values of Us? that dip below zero are set to zero.

Optimization

The following is a brief description (for full detail, see 81 Text). We start each optimization run With A = 500 “in-diViduals” that possess random assignments for each of the 15 parameters, chosen from a discrete set of numbers ({1 2 3 4 . . . 9} x 10{'2 '1 0}). (Based on our preliminary simulations, the starting intercompartmental exchange parameters and the dl/Cact dissociation rate (fig) are divided by 100 to avoid starting the search With unrealistically high values of these parameters.) The simulation results for nucleus 11 at time indeX k, Yfic(PN), corresponding to a set of parameters PN 2 [p1 p2 p3. . .p15]N, are evaluated and compared to the data using a residual sum of squares calculation: where th = dl-Venus data, clth 2 measurement uncertainty—both corresponding to simulated nucleus h and timepoint k—and SN is the ordinary least squares estimate of the scale factor that minimizes the difference between X and Y across all time and space coordinates and for parameter set N. Um“; and Wnuc are the dimensionless versions of nuclear dl and nuclear dl/ Cact complex, respectively (see Equations (1) and (3)). The dl-Venus data are taken directly from Reeves et al. [10], in which the concentration of nuclear Venus-tagged dl was imaged in live embryos and quantified in space and time and averaged.

The 500 “children” are then run through the simulation and ranked, and the top 100 individuals from the combined pool of parents + children are used to generate the next generation, and the process repeats. We have observed this process to stagnate at about 25 generations. After the 25th generation, we keep the top 100 parameter sets as the end product of each evolutionary optimization run.

For each optimization run, we use a set of parameters that are randomly chosen from the top 10% of those collected to simulate the dl/Cact dynamics. We use the dl/Cact dynamics associated with this set of parameters as an input to the gene eXpression model equations, and allow only the gene eXpression parameters to evolve for 9 different values of the noise parameter, 1], between 0.02 and 0.5. Since the concentration range of the simulated dl gradient is known to be between approximately 0 and 3 prior to optimization, the initial threshold (91-) parameters are restricted to random numbers between 0 and 3. For simplicity, the lifetime parameters (T1) are also restricted to this range, in minutes, but often evolve to values outside of that range. Again, we have observed this process to stagnate at about 25 generations. After the 25th generation, we keep the top 50 parameter sets as the end product of each evolutionary optimization run.

Results

Simulations of dI/Cact dynamics

To obtain an accurate picture of dl/Cact dynamics, we formulated a model of the dl gradient based on previous work ([15]; see 81 Text for details) and fit the model to a set of detailed spa-tiotemporal data for the dl gradient, obtained by live imaging of dl tagged with YFP variant Venus [10].

Specifically, we were able to simulate the overall increase in nuclear dl seen along the ventral midline (Fig. 2a), and the overall spatiotemporal shape of the simulated data set was in qualitative agreement with the dl-Venus data set (Fig. 2b; compare with Fig. 3a). However, our model could not simulate the decrease in nuclear dl seen in the dorsal-most nuclei (see Fig. 2a).

In the model, however, the initial conditions for each interphase include newly-formed nuclei that are devoid of any dl, Cact, or dl/Cact complex. This simplifying assumption precludes the possibility that dl/Cact complex can ever be found in the nuclei. Therefore, we formulated an “extended model” in which we removed this simplifying assumption but left the kinetic interactions the same. Specifically, our new model assumed that as nuclear envelopes reform in the cytoplasm, the contents of the nucleus would simply reflect that of the cytoplasm at the start of interphase (a straightforward assumption, given the porosity of nuclear envelopes during mitosis [23, 24]). This entailed modeling not only nuclear dl, but also nuclear Cact and nuclear dl/Cact complex, as each nascent nuclear compartment envelops all cytoplasmic species at the beginning of interphase. In other words, even if the nuclear import rates of Cact and dl/Cact complex were zero, the initial conditions for each interphase would include nonzero levels of these species in the nuclei.

If dl/Cact complex contributes to fluorescence and not gene expression, this result also provides a straightforward solution to the previously-noted difference between the dl gradient’s spatial range as measured by fluorescence and its spatial range as assayed by its ability to specify target gene locations (for example, see [10]). Therefore, under this extended model, we take nuclear dl fluorescence intensity to result from the sum of contributions from free dl and dl/Cact complex (see Equation 11). We thus fit the sum of free nuclear d1 (1) and nuclear dl/Cact complex (3) in our extended model to the dl-Venus data set in both space and time, revealing parameter sets that show an excellent fit to the data (see Fig. 3d). In particular, the extended model captures the decreasing basal levels, which affirms our assumption that nuclei begin interphase with the same concentration profile as the cytoplasm.

3a-c). Our model predicts that in the lateral and dorsal regions of the embryo, the observed levels of dl-Venus fluorescence, including the nonzero basal levels, are predominantly composed of dl/Cact complexes, and that the free nuclear dl gradient decays to near zero levels.

Parameter analysis of extended dI/Cact model

We represent each evolutionary optimization run by an average parameter set, in which we calculate the mean and standard deviation (weighted by RSS error; see Methods and 81 Text) for each parameter, across all 100 sets. We collected 254 such runs, resulting in 254 average parameter sets that are nearly indistinguishable in terms of their average RSS error (see Methods and 81 Text), yet vary in the values of their parameters, as has been observed previously for biological models [25, 26]. Despite the variation in parameters, several trends among the parameter sets have surfaced, most notably among the nuclear import/ export rates of the three molecular species.

4a). Deviations from this cluster are highly correlated between these two parameters, with a sparse, high-RSS-error tail. Thus, the nuclear import equilibrium constant (the ratio of import to export) for dl is largely contained within one order of magnitude and, perhaps unsurprisingly, favors import (Fig. 4b).

4c). However, it is clear from this plot that these two parameters are highly correlated. Indeed, their ratio, representing the nuclear import equilibrium constant for dl/Cact complex, is largely constrained to one order of magnitude, which slightly favors export (Fig. 4d).

4e,f). This is relatively unsurprising as, in our analysis, free Cact is not explicitly fit to any experimental measurements.

Gene expression simulations

While it was understood there were likely other factors involved in the expression of dl target genes, the modeling work revealed the extent to which the measured dl fluorescence alone could specify gene expression patterns. In particular, these simulations resulted in acceptable fits for Type I genes, adequate fits for Type II genes, and mediocre fits for Type III genes (which have dl-dependent boundaries past 40% DV axis; [10, 13]). The mediocre fits for the Type III genes are due to the measured spatial range of the dl nuclear fluorescence gradient being too narrow to account for the Type III genes. In light of our modeling work here, the other factors involved in DV gene regulation may include modulation of the dl nuclear gradient by nuclear Cact. Indeed, the difference in using fluorescence as a direct measurement of active dl (vs. total d1) becomes greater further away from the ventral midline, so that Type III genes such as 50g and zen are the most sensitive to this interpretation of the data. Therefore, in order to investigate how our proposed interpretation of the dl fluorescence data could impact our understanding of dl-dependent patterning, we built a model of gene expression.

To account for the error in gradient readout inherent to the stochastic process of dl arriving at the enhancer site, we added noise to the concentration of dl to calculate an effective dl concentration seen by each gene’s enhancer. For each nucleus 11, we cal-standard normal distribution and 17 is a tunable “noise” constant. A separate amount of noise was added for the calculation of each gene.

5a). In particular, our results indicate that the dl nuclear gradient can indeed pattern Type III genes, perhaps due to the fact that the nuclear concentration of free (active) dl drops by another order of magnitude past 40% DV aXis (Fig. 5b), in constrast to total dl, which drops by only 20% in the same range (Fig. 5d). As shown in Fig. 5b, the gene eXpression thresholds cross the NC14 dl gradient at about the same location as the half-maX of the corresponding gene’s boundary.

In those simulations, we observe that Type III genes cannot be properly specified by the total dl gradient (Fig. 5c). This is the same trend as seen previously [10]. Several results arise from this formulation that show the model is highly susceptible to noise. First, the locations where the predicted gene eXpression thresholds cross the dl gradient in this case do not align with the putative locations of the gene eXpression boundaries (Fig. 5d). In other words, much of gene expression does not result from a straightforward interpretation of the gradient, but instead may rely on local fluctuations in the effective concentration. Second, the dynamic range (defined as C(x1)/c(x2), for any choice of x1 < x2) of d1 gradient interpretation is very narrow, with genes of Types 1, II, and 111+ having thresholds all within a threefold range (Fig. 5d). Such a low dynamic range could prevent the different cell fates from being properly established. Third, the best-fit threshold for Type 111- genes is so low that the basal levels never drop below this threshold, and Type 111- genes are only eXpressed as a result of dl gradient read error (noise). This phenomenon is also observed in the Type III+ genes, as these simulations show they are eXpressed all the way to the dorsal midline.

Sensitivity analysis of gene expression model

To determine the sensitivity of our results using free dl with respect to changes in model parameters, we took the best fit parameters for both free dl and total dl and varied them by i 10%. In both of these scenarios, we used 17 = 0.2.

mRNA, TmRNA, n) for each gene tested (Fig. 6). In contrast, the Type III genes in the model using total dl are highly sensitive to the value of 9611:,” RNA as well as 1]. This result emphasizes that our understanding of how the dl gradient specifies genes far away from the ventral midline is contingent on the proper interpretation of dl fluorescence measurements. In our gene expression model, we idealize the relationship between dl concentration and gene expression rate as a hard-threshold (Hill coefficient 11 H = 100). Even with such a strict phenomenology, our extended model can reproduce the non-sharp boundaries of gene expression, as exhibited by the Type 111 genes and, to a lesser extent, the Type 11 genes [10]. Specifically, our results indicate that the graded nature of the boundaries of each gene may be a direct result of the amount of gradient read error (noise), which confirms preVious results based on using fluorescence directly as a measure of active dl [10]. However, if we model gene

Background subtraction mechanism allows the oil gradient to overcome noise

5), subtracting the inactive component of the dl gradient from that measured by fluorescence reveals an active dl gradient that is able to convey spatial information over a greater proportion of the DV aXis due to its eXpanded dynamic range. One may ask how this “background subtraction” mechanism achieves this. As illustrated in Fig. 7a, background subtraction via Cact increases the relative difference in active dl concentration between the ventral and dorsal midlines. This effect is achieved even if the dl/Cact profile is flat, in which case the shape and slope of the dl activity gradient is the same as that of the total dl gradient, but the dynamic range is greatly improved. According to our analysis, this becomes most important in the issue of overcoming noise. To quantify this effect, we calculate the potential for nuclei to misinterpret their position in the x-direction due to noise (Fig. 7b).

Suppose these nuclei are in a concentration gradient C(x), and a background subtracted gradient cb(x) = C(x) — B, where B is a constant less than min{c(x)}. The relative difference in concentration between the nuclei is E m l E | Ax = i E | an. In the same general region of the embryo, the relative amount of c 6 dx dx noise is % = In order for these two nuclei to perceive different enough dl concentrations to _ IdCID , which is much easier to overcome. a n In other words, for nuclei separated as k* 2 lg]; dx :4 the smaller k* need be, as long as the slope remains constant. This effect is plotted in Fig. 7b, in Which the solid lines are a representative dl gradient With more and more of the basal levels subtracted uniformly (red to blue), and the dotted lines represent the corresponding k*, the number of nuclei misinterpreting their position along the axis. Past 40% DV position, the value of k* is dramatically smaller With a larger background subtraction. Therefore, nuclei patterned nearest neighbors, the lower the value of c, by a “background subtracted” dl gradient, one in Which the nuclear dl levels are closer to zero (but the slopes are the same), Will be better able to distinguish themselves from their closely-separated neighbors.

Discussion

However, questions remain as to how the spatial information carried by a morphogen gradient results in gene expression boundaries. In this study, we attempt to address this question by analyzing a mechanistic model of dl/Cact dynamics in the early embryo in light of recent live imaging measurements of nuclear dl-Venus fluorescence [10]. In comparing our modeling results against the experimental data, we found that only when our model includes nuclear Cact and nuclear dl/Cact complex can it account for experimental observations such as the declining basal levels of dl-Venus fluorescence (Fig. 8a-c). In such a model, the most straightforward assumption is that fluorescence measurements are a combination of both free nuclear dl and dl/Cact complex (Fig. 8d). This assumption implies that the dl nuclear concentration gradient, as measured by fluorescence, represents total nuclear dl, and not necessarily the gradient of dl activity. This is especially true on the lateral and dorsal side of the embryo, where dl/Cact complex becomes the dominant dl-containing species.

When we distinguish between free nuclear dl and nuclear dl/Cact complex, the resulting active gradient possesses a larger dynamic range than when measured by fluorescence alone (see Fig. 7a), and thus is capable of generating expression patterns that accurately reflect what is measured in vivo. Indeed, this result is borne out in our simulations of gene expression patterns.

6; see also Sl—S4 Figs.). Using raw fluorescence becomes especially problematic in lateral and dorsal regions of the embryo, as the difference in dl-fluorescence measurements between neighboring nuclei becomes vanishingly small, resulting in lateral and dorsal nuclei becoming almost indistinguishable in their predicted gene expression (see Fig. 7b). For example, Fig. 6b shows that when 50g expression best fits the FISH data near 40% DV axis, the nuclei beyond 40% DV almost uniformly express 50g, creating a “tail” of 50g expression that does not match the data. However, raising the threshold value (Bdlzsog) in order to eliminate the “tail” results in a very poor fit in which the 50g domain collapses to more closely resemble the vnd domain, as described previously [10]. In addition, our optimization results found that, in the raw fluorescence case, the threshold for zen repression is actually below the concentration of dl, meaning that expression of zen is only permitted because signal noise in the dorsal half of the embryo causes the effective concentration of dl to occasionally dip below the zen threshold.

As speculated recently, this noise, plus a noise-filtering mechanism, may be the mechanism by which we get graded gene expression boundaries, even if the input-output response phenomenology is infinitely sharp [10]. This result is consistent with the observation that gene eXpression boundaries become more graded with increasing distance from the ventral midline, which is precisely what would be eXpected were noise to play a part in the boundary shape.

Including such factors would be beneficial in quantitatively understanding gene eXpression patterns, including timing and sharpness [30—32]. However, we do not anticipate these factors influencing our conclusion regarding the regulation of Type 111 genes, such as 50g and zen, as Twist is a short-ranged transcriptional partner of d1 [33], and Zelda eXpression is spatially uniform [34]. On the other hand, including Zelda in a gene eXpression model may help eXplain how a changing dl gradient can result in early gene eXpression [31].

However, our modeling work demonstrates that the presence of Cact in the nucleus is a straightforward conclusion from the observational data. This implies Cact may have a previously unappreciated function in the nucleus to regulate the effective concentration of (tran-scriptionally-active) dl. In mammalian systems, IKB does indeed enter the nucleus to regulate NF-KB molecules [38—41]. In those cases, it is newly-synthesized IKB, expressed as a result of NF-KB activation, that appears to regulate nuclear NF-KB. Circumstantial evidence suggests the cact locus may be zygotically regulated by dl activity, as binding sites for Twist, a transcrip-tional partner of dl, are functional within the cact locus [42, 43].

First, experiments could be performed to directly detect Cact within nuclei, such as by Western blot of nuclear extracts or fluorescence correlation spectroscopy. Second, a careful measurement of dl nuclear fluorescence in mutants with uniform levels of dl (and thus, uniform gene expression) would allow us to directly correlate dl nuclear fluorescence levels to gene expression. These measurements could be compared to those in wild type embryos to infer whether dl/Cact complex is a major contributor to dl fluorescence.

Such issues can be seen in other systems, notably in Dpp signaling in the Drosophila larval wing disc [44—46]. As a recent example, work in the Bicoid system has revealed a cofactor for Bicoid, called Dampened, that potentiates its activity [47]. Given the results of the present study, it is possible that Bicoid fluorescence readings do not accurately reflect its signaling gradient, similar to the case for dl presented here. Further experimental data is needed to verify the predictions made by our model; however, the modeling work presented here suggests that, in some cases, an accurate mathematical framework may be necessary to properly interpret fluorescence-based data.

Supporting Information

Detailed description of model formulation. (PDF)

(Left to right) Increasing the nosie parameter, 1] from 0 to 1 shows that the slopes of the gene expression boundaries approach infinity at 17 = O, and become very noise above 17 = 0.2. (Note: each run is an average of 10 runs for each parameter adjustment to reduce randomness in the plot due to noise. This comports With the experimental data, Which are the average of 10+ embryos. The same is true for 82—84 Figs.)

Using free d1 as the input to the gene expression model, a sensitivity analysis shows little sensitivity to changes in the dl threshold parameters (BdlmRNA), lifetime parameters (Ti), and noise parameter (1]) for our genes of interest. (Hill coefficient 11H 2 100.)

Using both free dl and dl/Cact complex as the input to the gene expression model, a sensitivity analysis shows high sensitivity to changes in the dl threshold parameters (BdlszNA) for both Type 111 genes (50g and zen; green and yellow, respectively), and little sensitivity to changes in lifetime and noise parameters. (Hill coefficient 11 H = 100). not change the conclusions of our sensitivity analysis. (TIF) Acknowledgments The authors are grateful to Sophia Carrell (NC State) for critical reading of the manuscript.

Author Contributions