A bistable autoregulatory module in the developing embryo commits cells to binary fates

Positive autoregulation has been repeatedly proposed as a mechanism for cells to adopt binary fates during embryonic development through bistability. However, without quantitatively determining their parameters, it is unclear whether the plethora of positive autoregulatory modules found within developmental gene regulatory networks are actually bistable. Here, we combine in vivo live imaging with mathematical modeling to dissect the binary cell fate dynamics of the fruit fly pair-rule gene fushi tarazu (ftz), which is regulated by two known enhancers: the early (non-autoregulating) element and the autoregulatory element. Live imaging of transcription and protein concentration in the blastoderm revealed that binary Ftz cell states are achieved as ftz expression rapidly transitions from being dictated by the early element to the autoregulatory element. Moreover, we discovered that Ftz concentration alone is insufficient to activate the autoregulatory element, and that this element only becomes responsive to Ftz at a prescribed developmental time. Based on these observations, we developed a dynamical systems model, and quantitated its kinetic parameters directly from experimental measurements. Our model demonstrated that the ftz autoregulatory module is indeed bistable and that the early element transiently establishes the content of the binary cell fate decision to which the autoregulatory module then commits. Further analysis in silico revealed that the autoregulatory element locks the Ftz expression fate quickly, within 35 min of exposure to the transient signal of the early element. Overall, our work confirms the widely held hypothesis that autoregulation can establish developmental fates through bistability and, most importantly, provides a framework for the quantitative dissection of cellular decision-making based on systems dynamics models and real-time measurements of transcriptional and protein dynamics.


. Bistability
A simplified dynamical systems model of the ftz autoregulatory element that ignores the dynamics of mRNA production describes the rate of change in Ftz concentration over time as where is Ftz concentration, ( ) is the gene regulatory (or input-output) function that describes how input Ftz concentration controls the rate of Ftz production, and is the Ftz degradation rate. Since Ftz promotes its own production, ( ) increases with . In the long term, Ftz concentrations will tend toward stable steady states, or attractors, for which (by definition) the change in concentration over time goes to 0 (i.e., = 0). Attractors are stable, meaning that if Ftz concentration is perturbed slightly away from the attractor, it will eventually return to the attractor. However, it is also possible to have unstable steady states for which, after a small perturbation, Ftz concentration will evolve away from the steady state. For the system in Eq. 1, a steady state * will solve * = ( * ), ( meaning that all steady states can be found graphically as the intersections between a line of slope and the function . In our case, the gene regulatory function is a sigmoidal function (red curves in Figure B1), such that there are between 1 and 3 intersections of the total degradation rate with the Ftz production rate ( ). Then the autoregulatory module is either monostable, meaning it possesses one attractor, or bistable, meaning it possesses two attractors and one unstable steady state between these attractors. Since the number of intersections depends on the shape of and the slope , the exact parameter values are crucial for determining the possible behaviors that can be exhibited by the autoregulatory module. We will also consider controlled systems, i.e., dynamical systems that have an additional regulatory input. In this case, this input corresponds to the Ftz concentration produced by the early element, , such that the dynamics of Ftz production by the autoregulatory element are given by ( ) = ( ( ) + ( )) − ( ).
Since the early element shuts down over time (that is, (∞) = 0), the steady states in Eq. 3 are the same as those for in Eq. 1. However, the transient dynamics of shape the trajectory will take and determine which steady state Ftz will ultimately reach. . 151 We observed that the early element already drives a relatively constant gene expression level 152 around 20 min prior to gastrulation (Figure 3B and C). Then, at 15 min before gastrulation, its tran-153 scriptional activity decreases significantly, resulting in a 60% reduction within the next 20 min of 154 development ( Figure 3B and C). Conversely, autoregulation is initiated 20 min prior to gastrulation, 155 with its activity increasing until gastrulation starts ( Figure 3C). This transition between the early 156 and autoregulatory elements occurs while binary cell states are being established (Figure 2). Since 2.2 ftz autoactivation is triggered at a specific developmental time 160 The tight transition between the early and autoregulatory elements that occurs within 20 min (Fig-161 ure 3) could be indicative of autoactivation being initiated by the increase in Ftz concentration 162 driven by the early element. However, autoactivation could also be triggered by upstream factors 163 at a specific developmental time, regardless of the Ftz level at that time point.

164
To distinguish between these two scenarios, we measured the gene regulatory function-the  We measured the gene regulatory function of the autoregulatory element by constructing an -20 -18 -16 -14 -12 -10 -8 -6 -4 -2 0 time before gastrulation (min) [Ftz] (au) autoregulatory output (au) Ftz protein transcription C Figure 4. Two color live-imaging reveals that ftz autoregulation is initiated at a specific developmental time. (A) Two-color tagging permits in vivo simultaneous visualization of input Ftz protein concentration using a LlamaTag and output autoregulatory transcriptional dynamics using a reporter carrying the MS2 system. (B) Representative frames from live-imaging data. Green and magenta channels correspond to Ftz concentration and the transcriptional output from the ftz autoregulatory element, respectively. Dark gray and light gray outlines correspond to cells expressing high or low levels of Ftz protein, respectively. (C) Illustrative single-cell trace of Ftz protein and autoregulatory activity. Green and magenta lines correspond to the Ftz protein and transcriptional activity of the autoregulatory element, respectively. Both protein and MS2 traces are smoothened using a moving average of 1 min. (D) Experimentally measured gene regulatory function of the ftz autoregulatory element between -10 min to -5 min relative to gastrulation. Grey points correspond to simultaneous measurements of Ftz and MS2 fluorescence at individual time points from single-cell traces at the anterior boundary of stripe 4 ( = 211 nuclei, example traces shown in Figure S1). These points were grouped into quantiles, and a Hill function (red line) was fit to the quantile means. (E) The autoregulatory input-output function evolves over time, as the ftz autoregulatory element transitions from an unresponsive to a responsive state within 15 min, indicating that Ftz autoregulation is initiated through a developmental time-based mechanism. Error bars shown indicate standard errors. All data are from = 7 embryos. response clearly following the stripe boundary ( Figure 4B). 181 To calculate the regulatory function, we restricted our analysis to the cells at the anterior bound-182 ary of stripe 4 as a means to minimize the influence of other position-dependent transcription fac-183 tors that might also contribute to ftz autoactivation (Schier and Gehring, 1993). We first extracted 184 two rows (high and low) of boundary cells. Then, we separated the input Ftz concentration and 185 output transcription from the autoregulatory element in the data corresponding to each individ-186 ual cell (Figure 4C; Figure S1) into ten quantiles and fit a Hill function to the quantile averages to get 187 the gene regulatory function of ftz autoregulation within a defined temporal range (for example,

188
-10 min to -5 min for Figure 4D). Our analysis revealed a sharp regulatory relationship between Ftz 189 protein and autoregulatory response, with a Hill coefficient of 3.3±0.8 ( Figure 4D). Such Hill coeffi-  (Gregor et al., 2007). 193 We repeated the process described above at multiple developmental times to analyze how 194 the regulatory function for ftz autoactivation evolves over time. The results, shown in Figure 4E, 195 revealed that the regulatory function is clearly distinct at different time points. Specifically, initially

203
As we argued in Figure 1B, the fact that Ftz can exhibit a high state at one point in time does not 204 necessarily imply that this high state will persist in the absence of upstream regulation. Determin-205 ing whether the autoregulatory module is bistable and hence possesses developmental memory 206 requires turning the schematic shown in Figure 1A into an explicit mathematical model with empir-207 ically determined parameter values. To that end, we first developed a dynamical systems model 208 for the full Ftz regulatory system (including both the early and autoregulatory elements) to verify 209 whether we could accurately recapitulate experimental results in silico. 210 Our measurements revealed that the autoregulatory element is unresponsive to Ftz concentra-211 tion until about = −20 min ( Figure 4E). Once the element becomes responsive, we can describe 212 the ftz autoregulatory module using the dynamical systems model given by where ( ) and ( ) are the mRNA and protein concentrations produced by the autoregulatory 214 module, , are the decay rates of mRNA and protein, respectively, and is the translation rate.  Example trajectories for the mRNA production rate ( ), mRNA concentration ( ), and protein concentration ( ) from the early element, as described by Eq. 5. (B) Representative traces for three nuclei comparing empirical total Ftz trajectories (dark green) to simulated trajectories (green) comprising contributions from both the early (blue) and autoregulatory (red) elements. Gray dashed line marks the experimentally determined threshold used to classify cells as high or low Ftz expression states at gastrulation ( Figure 2D). The plotted nuclei correspond to the circled and labeled points in panel C. Note that the simulations extend well past the time that we can obtain experimental measurements. (C) Each nucleus has a set of initial conditions ( ( ) = 0 , ( ) = 0 , ( ) = 0 ) for the dynamics of the early element. The blue surface separates those nuclei (circles) that are predicted to express low levels of Ftz at gastrulation (below the surface) from those that are predicted to express high levels of Ftz (above the surface). Color intensity corresponds to the Ftz concentration at gastrulation. Downward-facing red triangles indicate false positives (predicted high at gastrulation, but experimentally determined to be low) and upward-facing red triangles indicate false negatives (predicted low, but experimentally determined to be high). Results are shown for = 118 nuclei from 3 embryos at the anterior boundary of stripe 4, with model parameters given in Table S1. the protein resulting from the autoregulatory element through the translation of the mRNA pro-225 duced by this element (first term on the right-hand side) and degradation (second term on the 226 right-hand side). Finally, as indicated in Eq. 4, we assume that the autoregulatory element is only 227 active for ≥ . We note, however, that our model produced nearly identical results whether we 228 assumed that this transition to full responsiveness occurred instantaneously at time , or whether 229 we assumed a gradual increase in responsiveness over time (Section S 3.1). 230 After time , in addition to mathematically describing the expression dynamics of the autoreg-231 ulatory element, we can also model the contribution of the early element to ftz expression as a Ftz protein and the corresponding autoregulatory response as previously shown in Figure 4D.  ond, the translation rate was found by simultaneously measuring the rate of transcription and 251 the resulting protein concentration in a transgenic construct where ftz mRNA was labeled with MS2 252 and the resulting Ftz protein tagged with a LlamaTag ( Figure S3). Third, the decay in the transcrip-253 tion rate of the early element over time was determined by fitting an exponential function to 254 the transcriptional dynamics of the early element construct ( Figure S4). Fourth, the scaling factor 255 was inferred by systematically comparing simulated to measured traces in the regime where the 256 gene regulatory function was saturated and mRNA production rate was at its maximum ( Figure S5). 257 Finally, the mRNA and protein decay rates and were drawn from existing measurements in 258 the literature. All parameter values are reported in Table S1. 259 If our model of the ftz regulatory system is accurate, then for each nucleus along the anterior  Hence, for a given set of initial conditions, we can simulate the full dynamical system starting at such developmental memory is at play, we need to further analyze our theoretical model.

296
Having established that our model accurately predicts transient Ftz expression state at gastrulation, 297 we next analyzed the behavior of the early and autoregulatory modules separately, with the goal 298 of uncovering whether the autoregulatory module is bistable and, as a result, retains a long-term 299 binary memory of these Ftz levels. 300 We first performed a test of the autoregulatory module to ascertain whether it is bistable in The intersection between Ftz degradation rate and the effective Ftz production rate (i.e., late Ftz production rate adjusted for late ftz mRNA production and decay rates) reveals that the autoregulatory module is bistable, as described in Box 1 and Section S 1.3. (B) Most cells that express high Ftz state at gastrulation reach the high fate at steady state (cell 3), but some transiently express high Ftz at gastrulation before ultimately reaching the low fate (cell 2). The threshold (dashed gray line) is only used to determine state at gastrulation; fate is decided by which of the two stable steady states is approached by the system in simulation at long times. (C) The early module dictates the transient dynamics of Ftz as well as the fate it is predicted to adopt at infinite time once the autoregulatory element has reached steady state. The blue surface is repeated from Figure 5 and separates nuclei into low (below surface) and high (above surface) Ftz state at gastrulation. The red surface (switching separatrix) separates nuclei predicted to adopt the low Ftz fate (black circles) from those predicted to adopt the high Ftz fate (green circles), where by "fate" we mean expression level at infinite time (steady state). Nuclei between the blue and red surface (blue circles) are considered transiently high in that they express high Ftz at gastrulation but are predicted to adopt the low Ftz fate. Data are plotted for = 118 nuclei from 3 embryos. Red triangles indicate false negatives (upward) and false positives (downward) as determined from classification at gastrulation. (D) Results from a representative embryo show that the experimentally measured stripe pattern (left) is recapitulated by simulation (middle). A stripe pattern is still evident at gastrulation even from the predicted early Ftz concentration alone (right). Nuclear intensities at all time points are normalized to the predicted steady-state high Ftz concentration. Red "x"s denote nuclei with a transiently high Ftz state at gastrulation; triangles denote false positives (downward) and false negatives (upward). Parameters for simulation are as given in Table S1. Figure 6B. The remaining 6.3% of nuclei (5 of 79), despite having a high Ftz expression state at gas-328 trulation, were predicted to drop to a low Ftz fate in steady state as shown by cell 2 in Figure 6B. 329 No nuclei classified as low Ftz expression state at gastrulation were predicted to express high Ftz 330 after gastrulation. Thus, the autoregulatory element ensures that the vast majority of cells adopt 331 a fate matching the transient state at gastrulation.

332
While it is clear that the autoregulatory element establishes developmental memory by fixing 333 the Ftz expression fate in steady state, we wondered whether this memory was already at play 334 at gastrulation, or whether the early module is principally responsible for setting Ftz state at gas-335 trulation. Our simulations show that a stripe pattern is already evident at gastrulation from the 336 contribution of the early protein alone, ignoring the autoregulatory contribution ( Figure 6D). Thus,337 it appears that the anterior boundary of stripe 4 at gastrulation is defined by the regulatory activ-338 ity of upstream factors binding the early element in a manner that is largely independent of the 339 activity of the autoregulatory element. Therefore, this result supports the conclusion that the au-340 toregulatory module is bistable and that it primarily serves to commit cells to fates predetermined 341 by the early element. and, if so, to adopt a trajectory destined for a high steady-state expression fate. Thus, for the 359 results reported in this section, we restricted our analysis to a subset of the best predicted nuclei, 360 the nuclei that were correctly classified at gastrulation and have relatively low cumulative error 361 between simulation and experiment over time (Section S 3.1), that were also predicted to adopt 362 the high fate according to the switching separatrix analysis ( = 21; Figure 6B). 363 We define the commitment window as − , where indicates the start of autoregulatory In this work, we utilized live imaging approaches to quantitatively characterize the dynamics of 407 the fruit fly ftz regulatory system in vivo. We elucidated tight temporal coordination between the 408 two enhancer elements that regulate ftz expression (Figure 3), and combined dynamical systems , 412 we speculate that the approach employed by the Ftz system to decide cell fate is not limited to fruit 413 flies, but might also be widely adopted during development in other organisms.

414
One of our central discoveries is that ftz autoregulation is triggered at a specific developmental  (Soluri et al., 2020), and two, Caudal and Dichaete, have been shown to 420 bind directly to the eve autoregulatory element (MacArthur et al., 2009). We speculate that timer 421 genes might also bind the ftz autoregulatory element to trigger its responsiveness to Ftz protein. behavior of the whole network can then be predicted from the behavior of the modules in isolation.

430
How to define the modules, or equivalently how to break the network down into parts, has been 431 the subject of much discussion (Hartwell et al., 1999). As we have argued in this work, topology, represses its own production.

458
Throughout developmental biology, the concept of a commitment window has been repeatedly 459 utilized to describe the amount of time cells need to be exposed to upstream signals in order to de-  The fly lines used in this study were generated by inserting transgenic reporters into the fly genome 486 or by CRISPR-Cas9 genome editing, as described below. See Table S2 for detailed information on 487 the plasmid sequences used in this study.   (Cosentino and Bates, 2011). Here, unspecified upstream factors regulate the early element module, which produces output ( ). This, in turn, acts as the input to the autoregulatory module, which has been schematized as a latch with a clock signal (clk) to show that the element only becomes responsive at a particular time. The output ( ) can be summed with ( ) to recover total Ftz concentration. Relative to the more traditional schematic in (A), this representation includes dynamic signals ( ( ) and ( )), a timing element that triggers autoregulatory responsiveness, and a clear separation of inputs and outputs from each module. 569   Table S1. Quantitative 575 We use ( ) to describe the concentration of ftz mRNA transcribed from the early element, which 576 is translated into protein ( ). We define ( ) as the ftz mRNA transcribed from the autoregu-

Supplementary Tables
where and are the decay rates of mRNA and protein respectively, is the translation rate, 581 and is a scaling factor equivalent to the ratio of maximum production rate from the endogenous 582 locus vs. the transgene (wherẽis measured). Note that̃( , ⋅), the gene regulatory function for 583 the autoregulatory element, is time dependent. Unless otherwise stated, we will assume where is the sigmoid describing the autoregulatory relationship at maximum amplitude. Using the defi-586 nition of ( ) we will then write Eq. S1 as From our empirical measurements, we observed that the production rate ( ) of early ftz mRNA is 592 well approximated by an exponential decay (Figure S4), allowing us to model ( ) and ( ) 593 through the dynamical system As before, is the translation rate and and are the decay rates of mRNA and protein.

595
Because Eq. S5 is linear, it can be equivalently written as The analytical solution is then given by where 1 , 2 , 3 are the eigenvalues of corresponding to eigenvectors ⃗ 1 , ⃗ 2 , ⃗ 3 respectively. In 599 particular, Thus, the input ( ) to the autoregulatory element is completely characterized by three param-602 eters ( ( ), ( ), ( )) corresponding to the initial conditions for Eq. S5. Since every term 603 in the solution is multiplied by an exponential that decays in time, all state variables will tend to 0 604 as goes to infinity. 606 In Box 1 in the main text we give an example of how to identify stable steady states in a model that 607 only acknowledges protein concentration and ignores mRNA dynamics. Here, we show how to use 608 the graphical intersection test on systems including both mRNA and protein.

609
To reiterate, our dynamical system is described by Eq. S4 is given by which has no input and can therefore be analyzed for steady states by standard methods. (S15) * = ( * ). (S16) Hence, the intersections of a line of slope with the right-hand side give the steady-state late 619 protein concentrations * , from which we can recover * through Eq. S14. For a plot of the 620 intersection test, see Figure 6A in the main text. We calculated the regulatory function ( ( )) of the ftz autoregulatory element (shown in Figure 4D) 624 for the anterior boundary of stripe 4. To make this possible, we identified the boundary in a man- to -15, -15 to -10, -10 to -5, and -5 to 0 min) to obtain the trend shown in Figure 4E).

S 2.2 Translation rate 635
To calculate the translation rate, we simultaneously imaged the ftz transcription rate and the re- where ( ) is the transcription rate (i.e., the mRNA production rate reported by MS2 fluorescence), 643 and =0.099 min −1 (see Table S1). An example of a resulting prediction for the amount of mRNA 644 as a function of time is shown in Figure S3C. 645 Next, we performed a parameter sweep for the translation rate and, for each value of , we 646 integrated ( ) to predict the protein dynamics ( ) using the following equation (see Figure S3D 647 for sample nuclear traces) The translation rate that results in the best fit ( Figure S3D, red line) is recorded for each nu-649 cleus. We then calculated values for each embryo by averaging best-fitted for each single cell 650 ( Figure S3E). These values are averaged across = 129 cells (embryo 1) and = 119 cells (em-651 bryo 2), respectively. Then we averaged the resulting value of between two embryos, giving us 652 = 0.082±0.004 protein AU (mRNA AU min) -1 (see Table S1), which is used in our dynamical systems 653 model in the main text. Note that values larger than 0.2 au are omitted for accuracy ( = 2 nuclei). We first calculated the average for each individual embryo and then computed the average between these two embryos, resulting in = 0.082 ± 0.004 protein AU (mRNA AU min) -1 , which is used in our dynamical systems model.  Figure S4. Inferring the decay rate of the early element transcriptional activity . (A) Imaging transcriptional dynamics of the early element using the MS2 system. The construct is the same as the one described in Figure 3A. (B) Exponential fit from two embryos. Gray dots represent the averaged early element transcription rate at individual time points from a single embryo. Red line is the exponential fit that is used to estimate the early element decay rate . 655 We calculated the decay in the mRNA production rate of early element, described by the decay 656 rate , from fluorescence measurements of the early element MS2 reporter construct ( Figure S4A). 657 We performed an exponential fit to the average trajectory of each embryo (see Figure S4B for an 658 example fit) and averaged the resulting decay rates to obtain the mean value = 0.048 ± 0.0021 659 (see Table S1, = 2 embryos) that we used in the dynamical systems model. 661 We estimated , the ratio of maximum production rate from the endogenous locus vs. the trans-662 gene, by approximating the solution to the dynamical system in Eq. S4 between two time points 663 for which we have empirical data. In particular, we began from the experimental observation that 664̃( , ) plateaus shortly before gastrulation for nuclei with total Ftz levels above about ℎ ℎ ∞ = 665 1.5 × 10 6 a.u. Therefore, for nuclei that satisfy ( ) ≥ ℎ ℎ ∞ during this time (Figure S5A), we can 666 approximate the nonlinear system for the autoregulatory element by the following linear system

S 2.4 Scaling factor
where ( ℎ ℎ ∞ ) is now a constant. This system of equations has an analytical solution given by ) . (S20) Rearranging these expressions gives Note that here = 0 is assigned to the beginning of the time window over which the simulation is 670 performed.

671
We calculated the solution to Eq. S21 for the individual boundary nuclei in the same embryos 672 as used to fit the gene regulatory function (Figure 4D). We restricted our estimations to the 3 min 673 before gastrulation based on personal observations that the prediction accuracy of the simulation 674 of this linearized system tended to fall after ∼3 min. We derived the initial conditions ( (0), (0)) 675 from estimates of the late protein obtained by subtracting simulated early protein from the total 676 protein trace (where the autoregulatory element was assumed to begin contributing at -20 min 677 before gastrulation). We estimated in a windowed approach whereby, for each nucleus, we simu- From a histogram of the errors (Figure S6A), which was roughly bimodal, we identified a threshold 690 of 7 × 10 7 to identify the 79 best predicted nuclei out of 118 nuclei total. Some sample traces from 691 these best predicted nuclei are shown in Figure S6B while traces for nuclei with high cumulative 692 error are shown in Figure S6C. 693 We also calculated the cumulative error when simulations were conducted for a gradual in-694 crease in autoregulatory responsiveness described by where ( ) is as defined in Eq. S3 and = 0.14 min -1 corresponds to a half-life of 5 min (from 696 the observation that the autoregulatory element transitions from unresponsive to fully responsive 697 over 10 min; see Figure 4E). The resulting sample traces are shown in Figure S7B. Compared to 698 the case where we assume that the autoregulatory element becomes active instantaneously at 699 (Figure S6A), there was a slight shift in the distribution toward lower error (88 best predicted with 700 the same cutoff as before in Figure S6A), but the binary classification accuracy at gastrulation was 701 the same. Therefore we opted to use ( ) rather than̂( , ) for the main analysis. Simulations were conducted assuming that the autoregulatory element becomes instantaneously responsive, as in Eq. S3 and the main text. (B) Sample traces of best predicted nuclei. In each plot, the blue curve is the measured trace and the black line is the corresponding predicted trace. The dashed gray line is the threshold for classifying a nucleus as "on" at gastrulation. (C) As for (B), but for nuclei that did not qualify to be best predicted.  Figure 6C. Pixel intensities are normalized to the steady-state high Ftz level. 703 We were curious about how accurately our model, which is based on measurements at the anterior rates. We examined these questions using stochastic differential equations (SDEs), assuming that 745 the noise in our data arises solely from stochastic dynamics within the cells rather than from mea-746 surement noise. Generally, we expect this method to overestimate the error. We can perform the above analysis on individual measured traces to produce a large number 765 of sample points of ( ) and ( ) . With these data we will aim to estimate ( ) 766 and ( ). We assume the noise characteristics are time invariant, which allows us to pool all 767 samples at all time points and bin them by the corresponding or . For example, for protein, 768 we treat each protein concentration bin as a population of samples of (̄) where769 is the mean of the samples in bin (Figure S11A). Since we assume is normally distributed 770 with variance , we divide all sample values by √ and fit a normal distribution to the resulting 771 distribution within each protein concentration bin (Figure S11B). We found that the variances ( ) 2 772 are quite well approximated by a linear relation ( ) 2 = + (Figure S11C, right). The mRNA 773 variance was estimated similarly and also found to fit a linear relation (Figure S11C, left). Noise 774 was estimated from all available trajectories, regardless of whether they were part of the stripe 4 775 anterior boundary.

776
Having estimated the noise, we investigated whether stochasticity could explain the error rate where we compare the predicted outcomes from deterministic simulations to the "ground truth" 792 of the stochastic simulations.

793
In Figure S13 we plot the cumulative distribution of error rates for our simulations (black), bro- rate up to that rate. If the system is really described by the stochastic dynamics we have inferred, 798 then the most likely error rates are those that intersect the curves where their slope is highest. 799 We found that the empirical error rate across 3 embryos was roughly twice that of the most 800 likely error rates from our stochastic simulations corresponding to the middle of calculated cumu-801 lative distribution functions (Figure S13, left), with the majority of errors being false negatives. We 802 Figure S13. Error rates resulting from deterministic predictions made when gene expression is stochastic. Solid lines, cumulative distribution functions for error rates in the prediction of results for = 100 simulated stochastic experiments (using the stochastic differential equation in Eq. S24), when using a deterministic simulation to make the predictions. For each simulation, we predicted the trajectories of boundary nuclei evenly partitioned between starting states above and below the surface for classifying high or low Ftz at gastrulation (Figure 5). Dotted vertical lines correspond to empirical error rates for three embryos with = 118 nuclei (left) across 3 embryos, or excluding one embryo that had a high number of deterministic false negatives (Section 2.3) for = 78 nuclei (right) across 2 embryos.
knew from observation that one embryo had a large number of false negative predictions (Sec- depending upon how stringently one defines a threshold for switching.

814
Both trends indicate that switching from high to low occurs at a much faster rate than low to   Figure S14. First-passage times suggest that the stochastic switching rate between Ftz fates at steady state occurs on a timescale of hours. Left, time to first passage of a stochastic trace starting from the high (black) or low (red) fate to a given protein value. Right, time to first passage of a stochastic trajectory with a given Euclidean distance from the opposite steady state. In both plots, red denotes simulations beginning in the low Ftz state and black denotes simulations beginning in the high Ftz state.