Quantitative dissection of transcription in development yields evidence for transcription factor-driven chromatin accessibility

Thermodynamic models of gene regulation can predict transcriptional regulation in bacteria, but in eukaryotes chromatin accessibility and energy expenditure may call for a different framework. Here we systematically tested the predictive power of models of DNA accessibility based on the Monod-Wyman-Changeux (MWC) model of allostery, which posits that chromatin fluctuates between accessible and inaccessible states. We dissected the regulatory dynamics of hunchback by the activator Bicoid and the pioneer-like transcription factor Zelda in living Drosophila embryos and showed that no thermodynamic or non-equilibrium MWC model can recapitulate hunchback transcription. Therefore, we explored a model where DNA accessibility is not the result of thermal fluctuations but is catalyzed by Bicoid and Zelda, possibly through histone acetylation, and found that this model can predict hunchback dynamics. Thus, our theory-experiment dialogue uncovered potential molecular mechanisms of transcriptional regulatory dynamics, a key step toward reaching a predictive understanding of developmental decision-making.


Introduction
Over the last decade, hopeful analogies between genetic and electronic circuits have posed the challenge of predicting the output gene expression of a DNA regulatory sequence in much the same way that the output current of an electronic circuit can be predicted from its wiring diagram (Endy, 2005). 5 This challenge has been met with a plethora of theoretical works, including thermodynamic models, which use equilibrium statistical mechanics to calculate the probability of finding transcription factors bound to DNA and to relate this probability to the output rate of mRNA production (Ackers et al., 1982;Buchler et al., 2003;Vilar and Leibler, 2003;Bolouri and Davidson, 2003;Bintu 10 et al., 2005b,a;Sherman and Cohen, 2012). Thermodynamic models of bacterial transcription launched a dialogue between theory and experiments that has largely confirmed their predictive power for several operons (Ackers et al., 1982;Bakk et al., 2004;Zeng et al., 2010;He et al., 2010;Garcia and Phillips, 2011;Brewster et al., 2012;Cui et al., 2013;Brewster et al., 2014;Sepulveda et al., 15 2016; Razo-Mejia et al., 2018) with a few potential exceptions (Garcia et al., 2012;Hammar et al., 2014).
Following these successes, thermodynamic models have been widely applied to eukaryotes to describe transcriptional regulation in yeast (Segal et al., 2006;Gertz et al., 2009;Sharon et al., 2012;Zeigler and Cohen, 2014), human cells 20 (Giorgetti et al., 2010), and the fruit fly Drosophila melanogaster (Jaeger et al., 2004a;Zinzen et al., 2006;Segal et al., 2008;Fakhouri et al., 2010;Parker et al., 2011;White et al., 2012;Samee et al., 2015;Sayal et al., 2016). However, two key differences between bacteria and eukaryotes cast doubt on the applicability of thermodynamic models to predict transcriptional regulation in the latter. 25 First, in eukaryotes, DNA is tightly packed in nucleosomes and must become accessible in order for transcription factor binding and transcription to occur (Polach and Widom, 1995;Levine, 2010;Schulze and Wallrath, 2007;Lam et al., 2008;Li et al., 2011;Fussner et al., 2011;Bai et al., 2011;Li et al., 2014a;Hansen and O'Shea, 2015). Second, recent reports have speculated that, unlike 30 in bacteria, the equilibrium framework may be insufficient to account for the energy-expending steps involved in eukaryotic transcriptional regulation, such as histone modifications and nucleosome remodeling, calling for non-equilibrium models of transcriptional regulation (Kim and O'Shea, 2008;Estrada et al., 2016;Li et al., 2018;Park et al., 2019). This thermodynamic MWC model assumes that chromatin rapidly transitions between accessible and inaccessible states via thermal fluctuations, and that the binding of transcription factors to accessible DNA shifts this equilibrium toward the accessible state. Like all thermodynamic models, this model relies on the "occupancy hypothesis" (Hammar et al., 2014;Garcia et al., 2012;Phillips et al., 2019): the probability p bound of finding RNA polymerase (RNAP) bound to the promoter, a quantity that can be easily computed, is linearly related to the rate of mRNA production dmRNA dt , a quantity that can be experimentally measured, Here, R is the rate of mRNA production when the system is in an RNAPbound state (see Section S1.1 for a more detailed overview). Additionally, in all thermodynamic models, the transitions between states are assumed to be much faster than both the rate of transcriptional initiation and changes in transcription factor concentrations. This separation of time scales, combined with a lack 40 of energy dissipation in the process of regulation, makes it possible to consider the states to be in equilibrium such that the probability of each state can be computed using its Boltzmann weight (Garcia et al., 2007).
Despite the predictive power of thermodynamic models, eukaryotic transcription may not adhere to the requirements imposed by the thermodynamic Here, we performed a systematic dissection of the predictive power of these MWC models of DNA allostery in the embryonic development of the fruit fly Drosophila melanogaster in the context of the step-like activation of the hunchback gene by the Bicoid activator and the pioneer-like transcription factor Zelda 60 (Driever et al., 1989;Nien et al., 2011;Xu et al., 2014). Specifically, we compared the predictions from these MWC models against dynamical measurements of input Bicoid and Zelda concentrations and output hunchback transcriptional activity. Using this approach, we discovered that no thermodynamic or non- ular steps that make DNA accessible, to systematically test our model of transcription factor-driven chromatin accessibility, and to make progress toward a predictive understanding of transcriptional regulation in development.

A thermodynamic MWC model of activation and chromatin accessibility by Bicoid and Zelda
During the first two hours of embryonic development, the hunchback P2 minimal enhancer (Margolis et al., 1995;Driever et al., 1989;Perry et al., 2012;80 Park et al., 2019) is believed to be devoid of significant input signals other than activation by Bicoid and regulation of chromatin accessibility by both Bicoid and Zelda (Perry et al., 2012;Xu et al., 2014;Hannon et al., 2017). As a result, the early regulation of hunchback provides an ideal scaffold for a stringent test of simple theoretical models of eukaryotic transcriptional regulation. 85 Our implementation of the thermodynamic MWC model (Fig. 1A) in the context of hunchback states that in the inaccessible state, neither Bicoid nor Zelda can bind DNA. In the accessible state, DNA is unwrapped and the binding sites become accessible to these transcription factors. Due to the energetic cost of opening the chromatin (∆ε chrom ), the accessible state is less likely to occur 90 than the inaccessible one. However, the binding of Bicoid or Zelda can shift the equilibrium toward the accessible state (Adams and Workman, 1995;Miller and Widom, 2003;Mirny, 2010;Narula and Igoshin, 2010;Marzen et al., 2013).
In our model, we assume that all binding sites for a given molecular species have the same binding affinity. Relaxing this assumption does not affect any 95 of our conclusions (as we will see below in Sections 2.3 and 2.4). Bicoid upregulates transcription by recruiting RNAP through a protein-protein interaction characterized by the parameter ω bp . We allow cooperative protein-protein interactions between Bicoid molecules, described by ω b . However, since to our knowledge there is no evidence of direct interaction between Zelda and any 100 other proteins, we assume no interaction between Zelda and Bicoid, or between Zelda and RNAP.
In Fig. 2A, we illustrate the simplified case of two Bicoid binding sites and one Zelda binding site, plus the corresponding statistical weights of each state given by their Boltzmann factors . Note that the actual model utilized through-105 out this work accounts for at least six Bicoid binding sites and ten Zelda binding sites that have been identified within the hunchback P2 enhancer (Section 4.1;Driever and Nusslein-Volhard, 1988b;Driever et al., 1989;Park et al., 2019). This general model is described in detail in Section S1.2.
The probability of finding RNAP bound to the promoter is calculated by dividing the sum of all statistical weights featuring RNAP by the sum of the weights corresponding to all possible system states. This leads to Bicoid and RNAP binding , [Bicoid], [Zelda], and [RN AP ] being the concentrations of Bicoid, Zelda, and RNAP, respectively, and K b , K z , and K p their dissociation constants (see Sections S1.1 and S1.2 for a detailed derivation). Given a set of model parameters, plugging p bound into Equation 1 predicts the rate of RNAP loading as a function of Bicoid and Zelda concentrations as shown in Fig. 2B. Note that in this work, 115 we treat the rate of transcriptional initiation and the rate of RNAP loading interchangeably.

Dynamical prediction and measurement of input-output functions in development
In order to experimentally test the theoretical model in Fig. 2, it is nec-120 essary to measure both the inputs -the concentrations of Bicoid and Zeldaas well as the output rate of RNAP loading. Typically, when testing models of transcriptional regulation in bacteria and eukaryotes, input transcriptionfactor concentrations are assumed to not be modulated in time: regulation is in steady state (Ackers et al., 1982;Bakk et al., 2004;Segal et al., 2008;Garcia 125 and Phillips, 2011;Sherman and Cohen, 2012;Cui et al., 2013;Little et al., 2013;Raveh-Sadka et al., 2009;Sharon et al., 2012;Zeigler and Cohen, 2014;Xu et al., 2015;Sepulveda et al., 2016;Estrada et al., 2016;Razo-Mejia et al., 2018;Zoller et al., 2018;Park et al., 2019). However, embryonic development is a highly dynamic process in which the concentrations of transcription factors are 130 constantly changing due to their nuclear import and export dynamics, and due to protein production, diffusion, and degradation (Edgar and Schubiger, 1986;Edgar et al., 1987;Jaeger et al., 2004b;Gregor et al., 2007b Given the high reproducibility of the concentration dynamics of Bicoid and 155 Zelda (Fig. S3), we combined measurements from multiple embryos by synchronizing their anaphase in order to create an "averaged embryo" (Section S2.1), an approach that has been repeatedly used to describe protein and transcriptional dynamics in the early fly embryo (Garcia et al., 2013;Bothma et al., 2014Bothma et al., , 2015Berrocal et al., 2018;Lammers et al., 2019).  In order to test the predictions of the thermodynamic MWC model ( Fig. 3E and F), we used the MS2 system (Bertrand et al., 1998;Garcia et al., 2013;Lucas et al., 2013). Here, 24 repeats of the MS2 loop are inserted in the To compare theory and experiment, we next obtained the initial RNAP loading rates (Fig. 4C, points) and t on (Fig. 4D, points) from the experimental 200 data (Section S2.3; Fig. S5B). The step-like shape of the RNAP loading rate ( Fig. 4C) agrees with previous measurements performed on this same reporter construct (Garcia et al., 2013). The plateaus at the extreme anterior and posterior positions were used to constrain the maximum and minimum theoretically allowed values in the model (Section S1.3). With these constraints in place, we 205 attempted to simultaneously fit the thermodynamic MWC model to both the initial rate of RNAP loading and t on . For a given set of model parameters, the measurements of Bicoid and Zelda concentration dynamics predicted a corresponding initial rate of RNAP loading and t on ( Fig. 3E and F). The model parameters were then iterated using standard curve-fitting techniques (Section 4.6) 210 until the best fit to the experimental data was achieved ( Fig. 4C and D, lines).
Although the model accounted for the initial rate of RNAP loading (Fig. 4C, line), it produced transcriptional onset times that were much lower than those that we experimentally observed (Fig. S6B, purple line). We hypothesized that this disagreement was due to our model not accounting for mitotic repression, when the transcriptional machinery appears to be silent immediately after cell division (Shermoen and O'Farrell, 1991;Gottesfeld and Forbes, 1997;Parsons and Tg, 1997;Garcia et al., 2013). Thus, we modified the thermodynamic MWC model to include a mitotic repression window term, implemented as a time window at the start of the nuclear cycle during which no transcription could occur; the rate of mRNA production is thus given by where R and p bound are as defined in Eqns. 1 and 2, respectively, and t M itRep is the mitotic repression time window over which no transcription can take place after anaphase (Sections S1.2 and S3.1). After incorporating mitotic repression, the thermodynamic MWC model successfully fit both the rates of RNAP loading 215 and t on ( Fig. 4C and D, lines, Fig. S6A and B).
Given this success, we next challenged the model to perform the simpler task of explaining Bicoid-mediated regulation in the absence of Zelda. This scenario corresponds to setting the concentration of Zelda to zero in the models in Section S1.2 and Fig. 2. In order to test this seemingly simpler model, we re-220 peated our measurements in embryos devoid of Zelda protein (Video S4). These zelda − embryos were created by inducing clones of non-functional zelda mutant (zelda 294 ) germ cells in female adults (Sections 4.2, 4.3;Liang et al., 2008). All embryos from these mothers lack maternally deposited Zelda; female embryos still have a functional copy of zelda from their father, but this copy is not tran-225 scribed until after the maternal-to-zygotic transition, during nuclear cycle 14 (Liang et al., 2008). We confirmed that the absence of Zelda did not have a substantial effect on the spatiotemporal pattern of Bicoid (Section S4.1;Xu et al., 2014).
While close to 100% of nuclei in wild-type embryos exhibited transcription  positions in the mutant embryos that did exhibit transcription in at least 50% of observed nuclei, we extracted the initial rate of RNAP loading and t on as a function of position. Interestingly, these RNAP loading rates were comparable to the corresponding rates in wild-type embryos (Fig. 4F, points). However, unlike in the wild-type case (Fig. 4D, points), t on was not constant in the zelda −  to the next mitosis have already started and repressed transcriptional activity.
Next, we attempted to simultaneously fit the model to the initial rates of RNAP loading and t on in the zelda − mutant background. Although the model recapitulated the observed initial RNAP loading rates, we noticed a discrepancy between the observed and fitted transcriptional onset times of up to ∼5 min 255 ( Fig. 4F and G). While the mutant data exhibited a substantial delay in more posterior nuclei, the model did not produce significant delays (Fig. 4G, line). Further, our model could not account for the lack of transcriptional activity posterior to 40% of the embryo length in the zelda − mutant (Fig. 4H).
These discrepancies suggest that the thermodynamic MWC model cannot 260 fully describe the transcriptional regulation of the hunchback promoter by Bicoid and Zelda. However, the attempted fits in Fig. 4F and G correspond to a particular set of model parameters and therefore do not completely rule out the possibility that there exists some parameter set of the thermodynamic MWC model capable of recapitulating the zelda − data.

265
In order to determine whether this model is at all capable of accounting for the zelda − transcriptional behavior, we systematically explored how its parameters dictate its predictions. To characterize and visualize the limits of our model, we examined a mathematically convenient metric: the average transcriptional onset delay along the anterior-posterior axis (Fig. 5A). This quantity is defined as the area under the curve of t on versus embryo position, from 20% to 37.5% along the embryo (the positions where the zelda − data display transcription), divided by the corresponding distance along the embryo where x is the position along the embryo and the value of t on at 20% along the embryo was chosen as the offset with respect to which to define the zero of this integral (Section S5.1). The average t on delay corresponding to the wildtype data is close to zero, and is substantially different from the larger value obtained from measurements in the zelda − background within experimental and ω bp , energy to make the DNA accessible ∆ε chrom , and length of the mi-275 totic repression window t M itRep ), and to calculate the corresponding average t on delay for each parameter set. Fig. 5B features three specific realizations of this parameter search; for each combination of parameters considered, the predicted t on is calculated and the corresponding average t on delay computed. Although the wild-type data were contained within the thermodynamic MWC model re-280 gion, the range of the average t on delay predicted by the model (Fig. 5C, green rectangle) did not overlap with the zelda − data. We concluded that our thermodynamic MWC model is not sufficient to explain the regulation of hunchback by Bicoid and Zelda.
where p inacc and P r,i are arbitrary weights describing the states in our gener- Qualitatively, the reason for the failure of thermodynamic models to predict 305 hunchback transcriptional is revealed by comparing Bicoid and Zelda concentration dynamics to those of the MS2 output signal (Fig. S10). The thermodynamic models investigated in this work have assumed that the system responds instantaneously to any changes in input transcription factor concentration. As a result, since Bicoid and Zelda are imported into the nucleus by around 3 min into the 310 nuclear cycle ( Fig. 3A and B), these models always predict that transcription will ensue at approximately that time. Thus, thermodynamic models cannot accommodate delays in the t on such as those revealed by the zelda − data (see Section S6.3 for a more detailed explanation). Rather than further complicating our thermodynamic models with additional molecular players to attempt 315 to describe the data, we instead decided to examine the broader purview of non-equilibrium models to attempt to reach an agreement between theory and experiment.
2.5. A non-equilibrium MWC model also fails to describe the zelda − data Thermodynamic models based on equilibrium statistical mechanics can be Since this model can operate out of steady state, we calculate the probabilities of each state as a function of time by solving the system of coupled ordinary differential equations (ODEs) associated with the system shown in

385
We therefore proposed a model of transcription factor-driven chromatin accessibility in which, in order for the DNA to become accessible and transcription to ensue, the system slowly and irreversibly transitions through m transcriptionally silent states (Section S8.1; Fig. 7A). We assume that the transitions between these states are all governed by the same rate constant π. Finally, in a stark deviation from the MWC framework, we posit that these transitions can be catalyzed by the presence of Bicoid and Zelda such that Here, π describes the rate (in units of inverse time) of each irreversible step, expressed as a sum of rates that depend separately on the concentrations of Bicoid and Zelda, and c b and c z are rate constants that scale the relative contribution of each transcription factor to the overall rate.
In this model of transcription factor-driven chromatin accessibility, once the 390 DNA becomes irreversibly accessible after transitioning through the m nonproductive states, we assume that, for the rest of the nuclear cycle, the system equilibrates rapidly such that the probability of it occupying any of its possible (B-D, error bars represent standard error of the mean over multiple embryos).

Discussion
For four decades, thermodynamic models rooted in equilibrium statistical mechanics have constituted the null theoretical model for calculating how the We conclude that the MWC thermodynamic model cannot accurately predict hunchback transcriptional dynamics.
To explore non-equilibrium models, we retained the MWC mechanism of chromatin accessibility, but did not demand that the accessible and inaccessible These models could be tested with future experiments that implement these 570 modulations in reporter constructs.
In sum, here we engaged in a theory-experiment dialogue to respond to the theoretical challenges of proposing a passive MWC mechanism for chromatin accessibility in eukaryotes (Mirny, 2010;Narula and Igoshin, 2010;Marzen et al., 2013); we also questioned the suitability of thermodynamic models in the con-  new insights afforded by this dialogue will undoubtedly refine theoretical descriptions of transcriptional regulation as a further step toward a predictive understanding of cellular decision-making in development.

Acknowledgments
We are grateful to Jack Bateman, Jacques Bothma, Mike Eisen, Jeremy  , 1990;Hertz and Stormo, 1999). PATSER was run with setting "Seq. Alphabet and Normalization" as "a:t 3 g:c 2" to provide the approximate background frequencies as annotated in the Berkeley Drosophila Genome Project (BDGP)/Celera Release 1. Reverse complementary sequences were also scored. In order to image transcription in embryos lacking maternally de-620 posited Zelda protein, we crossed mother flies whose germline was w, his2av-mrfp1,zelda (294) bicoid transgene used for imaging nuclear Bicoid in a wildtype background.

Zelda germline clones 635
In order to generate mother flies containing a germline homozygous null for zelda, we first crossed virgin females of w,his2av-mrfp1,zelda (294)

+ [P ]
Kp Here, [P ] and [A] are the concentrations of RNAP and activator, respectively.
K p and K a are their corresponding dissociation constants, and ω indicates an 1075 interaction between activator and RNAP: ω > 1 corresponds to cooperativity, whereas 0 < ω < 1 corresponds to anti-cooperativity.
Using p bound , we write the subsequent rate of mRNA production by assuming the occupancy hypothesis, which states that Figure S1: Equilibrium thermodynamic model of simple activation. A promoter region with one binding site for an activator molecule has four possible microstates, each with its corresponding statistical weight and rate of RNAP loading.
where R is an underlying rate of transcriptional initiation (usually interpreted as the rate of loading RNAP from the promoter-bound state). In the case of simple activation illustrated in Fig. S1  Hill equation is invoked, equilibrium thermodynamics is implicitly used, bringing with it all of the underlying assumptions described in Section S6.4. This highlights the importance of rigorously grounding the assumptions made in any model of transcription, to better discriminate between the effects of equilibrium and non-equilibrium processes.

S1.2. Thermodynamic MWC model
In the thermodynamic MWC model, we consider a system with six Bicoid binding sites and ten Zelda binding sites. In addition, we allow for RNAP binding to the promoter.
In our model, the DNA can be in either an accessible or an inaccessible state.
The difference in free energy between the two states is given by −∆ε chrom , where ∆ε chrom is defined as Here, ε accessible and ε accessible are the energies of the accessible and inaccessi-1090 ble states, respectively. A positive ∆ε chrom signifies that the inaccessible state is at a lower energy level, and therefore more probable, than the accessible state.
We assume that all binding sites for a given molecular species have the same binding affinity, and that all accessible states exist at the same energy level compared to the inaccessible state. Thus, the total number of states is determined 1095 by the combinations of occupancy states of the three types of binding sites as well as the presence of the inaccessible, unbound state. We choose to not allow any transcription factor or RNAP binding when the DNA is inaccessible.
In this equilibrium model, the statistical weight of each accessible microstate is given by the thermodynamic dissociation constants K b , K z , and K p of Bicoid, 1100 Zelda, and RNAP respectively. The statistical weight for the inaccessible state is e ∆ε chrom . We allow for a protein-protein interaction term ω b between nearestneighbor Bicoid molecules, as well as a pairwise cooperativity ω bp between Bicoid and RNAP. However, we posit that Zelda does not interact directly with either Bicoid or RNAP. For notational convenience, we express the statistical weights  Figure S2: States, weights, and rate of RNAP loading diagram for the thermodynamic MWC model, containing six Bicoid binding sites, ten Zelda binding sites, and a promoter.
shows the states and statistical weights for this thermodynamic MWC model, with all the associated parameters.
Incorporating all the microstates, we can calculate a statistical mechanical partition function, the sum of all possible weights, which is given by (1 + z) 10 Zelda binding Bicoid and RNAP binding .
Using the binomial theorem Eq. S7 can be expressed more compactly as From this partition function, we can calculate p bound , the probability of being in an RNAP-bound state. This term is given by the sum of the statistical weights of the RNAP-bound states divided by the partition function In this model, we once again assume that the transcription associated with each microstate is zero unless RNAP is bound, in which case the associated rate is R. Then, the overall transcriptional initiation rate is given by the product of (S10) Note that since the MS2 technology only measures nascent transcripts, we can 1110 ignore the effects of mRNA degradation and focus on transcriptional initiation.

S1.3. Constraining model parameters
The transcription rate R of the RNAP-bound states can be experimentally constrained by making use of the fact that the hunchback minimal reporter used in this work produces a step-like pattern of transcription across the length of the fly embryo (Fig. 4C). Since in the anterior end of the embryo, the observed transcription appears to level out to a maximum value, we assume that Bicoid binding is saturated in this anterior end of the embryo such that In this limit, Eq. S10 can be written as where R max is the maximum possible transcription rate. Importantly, R max is an experimentally observed quantity rather than a free parameter. As a result, the model parameter R is determined by experimentally measurable quantity 1115 R max .
The value of p can also be constrained by measuring the transcription rate in the embryo's posterior, where we assume Bicoid concentration to be negligible. Here, the observed transcription bottoms out to a minimum level R min (Fig. 4C), which we can connect with the model's theoretical minimum rate.
Specifically, in this limit, b approaches zero in Eq. S10 such that all Bicoiddependent terms drop out, resulting in where we have replaced R with R max as described above. Next, we can express p in terms of the other parameters such that Thus, p is no longer a free parameter, but is instead constrained by the experimentally observed maximum and minimum rates of transcription R max and R min , as well as our choices of K z and ∆ε chrom . In our analysis, R max and R min are calculated by taking the mean RNAP loading rate across all embryos 1120 from the anterior and posterior of the embryo respectively, extrapolated using the trapezoidal fitting scheme described in Section S2.3.
Finally, we expand this thermodynamic MWC model to also account for suppression of transcription in the beginning of the nuclear cycle via mechanisms such as mitotic repression (Section S3.1). To make this possible, we include a 1125 trigger time term t M itRep , before which we posit that no readout of Bicoid or Zelda by hunchback is possible and the rate of RNAP loading is fixed at 0. For times t > t M itRep , the system behaves according to Eq. S10. Thus, given the constraints stemming from direct measurements of R max and R min , the model has six free parameters: ∆ε chrom , ω b , ω bp , K b , K z , and t M itRep . The final 1130 calculated transcription rate is then integrated in time to produce a predicted MS2 fluorescence as a function of time (Section S2.2).
For subsequent parameter exploration of this model (Section S5.1), constraints were placed on the parameters to ensure sensible results. Each parameter was constrained to be strictly positive such that:  , 2008). Intra-embryo variability in Zelda-sfGFP nuclear fluorescence was very low (less than 10%), whereas inter-embryo variability was relatively higher, up to 20% (Fig. S3F, red and black, respectively). Three separate Zelda-sfGFP embryos were measured.
Due to the consistency of Zelda-sfGFP nuclear fluorescence, we assumed 1180 the Zelda profile to be spatially uniform in our analysis, and thus created a mean Zelda-sfGFP measurement for each individual embryo by averaging all mean nuclear fluorescence traces in space across the anterior-posterior axis of the embryo (Fig. S3D, inset). This mean measurement was used as an input in the theoretical models. However, we still retained inter-embryo variability in 1185 Zelda, as described below.
To combine multiple embryo datasets as inputs to the models explored throughout this work, the fluorescence traces corresponding to each dataset were aligned at the start of nuclear cycle 13, defined as the start of anaphase.
Because each embryo may have possessed slightly different nuclear cycle lengths (S15) Finally, we assume that upon reaching the end of the reporter gene, the RNAP 1215 molecules terminate and disappear instantly such that they no longer contribute to spot fluorescence.
The MS2 fluorescence signal reports on the number of RNAP molecules actively occupying the gene at any given time and, under the assumptions outlined above, is given by the integral where F (t) is the predicted fluorescence value, R(t) is the RNAP loading rate predicted by each specific model, R(t − t elon ) is the time-shifted loading rate that accounts for RNAP molecules finishing transcription at the end of the gene, 1220 and α is an arbitrary scaling factor to convert from absolute numbers of RNAP molecules to arbitrary fluorescence units. The predicted value F (t) was scaled by α to match the experimental data.
The final predicted MS2 signal was modified in a few additional ways. First, any RNAP molecule that had not yet reached the position of the MS2 stem loops had its fluorescence value set to zero (Fig. S4, i), since only RNAP molecules downstream of the MS2 stem loop sequence exhibit a fluorescent signal. Second, RNAP molecules that were only partially done elongating the MS2 stem loops contributed a partial fluorescence intensity, given by the ratio of the distance traversed through the stem loops to the total length of the stem loops Figure S4: MS2 fluorescence calculation protocol. RNAP molecules load onto the reporter gene at a time-dependent rate R(t), after which they elongate at a constant velocity v. Upon reaching the end of the gene after a length L has been transcribed, they are assumed to terminate and disappear instantly, given by the time-shifted rate R(t − L v ). The time-dependent MS2 fluorescence is calculated by summing the contributions of RNAP molecules that are located before, within, or after the MS2 stem loop sequence (i, ii, and iii, respectively).
where F partial is the partial fluorescence contributed by an RNAP molecule within the stem loop sequence region, L partial is the distance within the stem 1225 loop sequence traversed, and L loops is the length of the stem loop sequence (Fig. S4, ii). For this reporter construct, the length of the stem loops was approximately L loops = 1.28 kb. RNAP molecules that had finished transcribing the MS2 stem loops contributed the full amount of fluorescence (Fig. S4, iii).
Finally, to make this simulation compatible with the trapezoidal fitting scheme 1230 in Section S2.3, we included a falling signal at the end of the nuclear cycle, achieved by setting R(t) = 0 after 17 min into the nuclear cycle and thus preventing new transcription initiation events.
Given the predicted MS2 fluorescence trace, the rate of RNAP loading and t on were extracted with the fitting procedure used on the experimental data 1235 (Section S2.3).

S2.3. Extracting initial RNAP loading rate and transcriptional onset time
To extract the initial rate of RNAP loading and the transcriptional onset time t on used in the data analysis, we fit both the experimental and calculated The trapezoidal model provides a heuristic fit of the main features of the MS2 signal by assuming that the RNAP loading rate is either zero or some constant value r (Fig. S5A). At time t on , the loading rate switches from zero to this constant value r, producing a linear rise in the MS2 signal. After the  With this fit, we can extract the initial rate of RNAP loading (given by the 1255 initial slope) as well as t on (given by the intercept of the fit onto the x-axis).
As a consistency check, the t on values extrapolated from the trapezoidal fit of the data were compared with the experimental time points at which the first MS2 spots were observed for both the wild-type and zelda − mutant experiments (Fig. S5C). Due to the detection limit of the microscope, this latter 1260 method reports on the time at which a few RNAP molecules have already begun transcribing the reporter gene, rather than a "true" transcriptional onset time. Using the first frame of spot detection yields similar trends to the trapezoidal fits, except that the measured first frame times are systematically larger by 3-5 min. Additionally, utilizing the first frame of detection to measure t on 1265 appears to be a noisier method, likely because the actual MS2 spots cannot be observed below a finite signal-detection limit, whereas the extrapolated t on from the trapezoidal fit corresponds to a "true" onset time below the signal-detection limit.

1270
S3.1. Mitotic repression is necessary to recapitulate Bicoid-and Zelda-mediated regulation of hunchback using the thermodynamic MWC model As described in Section 2.3 of the main text, a mitotic repression window was incorporated into the thermodynamic MWC model (Section S1.2) in order to explain the observed transcriptional onset times of hunchback. Here, we justify 1275 and explain this theoretical modification in greater detail. Fig. S6A and B depicts the experimentally observed initial rates of RNAP loading and t on across the length of the embryo (blue points) for the wildtype background. After constraining the maximum and minimum theoretically allowed rates of RNAP loading (Section S1.3), we attempted to simultaneously 1280 fit the thermodynamic MWC model to both the rate of RNAP loading and t on .
The fit results demonstrate that while the thermodynamic MWC model can recapitulate the measured step-like rate of RNAP loading at hunchback Bicoid and Zelda concentration measurements at 45% along the embryo, shown for a single embryo in Fig. S6C, are used in conjunction with the previously mentioned best-fit model parameters to predict the output MS2 signal at the same position. This prediction can then be directly compared with experimental data (Fig. S6D, purple line vs. black points, respectively). Whereas the model 1295 predicts that transcription will commence around 1 min after anaphase due to the concurrent increase in the Bicoid and Zelda concentrations, the observed MS2 signal begins to increase around 3 min after anaphase (Fig. S6D). As a result, the predicted transcriptional dynamics in Fig. S6D are systematically shifted in time with respect to the observed data.

1300
The observed disagreement in t on suggests that in this model, transcription is prevented from starting at the time dictated solely by the increase of Bicoid and Zelda concentrations. While we speculate that this effect could stem from processes such as RNAP escape from the promoter, DNA replication at the start of the cell cycle, and post-mitotic nucleosome clearance from the promoter, we 1305 choose not to commit to a detailed molecular picture and instead ascribe this transcriptional refractory period at the beginning of the nuclear cycle to mitotic repression, the observation that the transcriptional machinery cannot operate during mitosis (Shermoen and O'Farrell, 1991;Gottesfeld and Forbes, 1997;Parsons and Tg, 1997;Garcia et al., 2013). To account for this phenomenon, 15.46% ± 1.64%. This value is within the range of the inter-embryo variability of eGFP Bicoid in wild-type background embryos (Fig. S3C). Measuring the decay length of the eGFP-Bicoid profile in the zelda − background also yielded a slightly different result: 20.8% +/-1.4% of the total embryo length, as opposed to 23.5% +/-0.6% in the wild-type background (dashed curves).

1345
Nevertheless, these differences would have a negligible effect on our overall conclusions. In the context of our models, an overall rescaling in the magnitude of the Bicoid gradient between the wild-type and zelda − backgrounds can be compensated by a corresponding rescaling in the dissociation constant of Bicoid, K b . Because our systematic exploration of theoretical models considers many 1350 possible parameter values (Section S5.1), this rescaling has no effect on our conclusion that the equilibrium models are insufficient to explain the zelda − data. As a result and given that our statistics for the wild-type eGFP-Bicoid data consisted of more embryos than the data for the zelda − background, we used this wild-type data in our analyses as an input to both the wild-type and 1355 zelda − model calculations. wild-type (+zelda) mutant (-zelda) rescaled mutant exponential fit (wild-type) exponential fit (mutant) λ = 23.5% +/-0.6% λ = 20.8% +/-1.4% Figure S7: eGFP-Bicoid measurements in wild-type (blue) and zelda − mutant embryos (red), along with rescaled mutant profiles (black). Fits to an exponentially decaying function yield decay lengths in each background (blue and black dashed curves). A total of n=3 embryos were measured. All error bars are standard error of mean across embryos.

S5.1. General methodology
To help visualize the limits of our models, we collapsed our observations onto a two-dimensional state space, following a method similar to that described 1360 in Estrada et al. (2016). In this space, the x-axis is the average t on delay. This magnitude was computed by integrating the t on across 20% to 37.5% of the embryo length, corresponding to the range in which both wild-type and zelda − experiments exhibited transcription in at least 50% of observed nuclei (Figs. S8A and 5A; Eq. 4). We first subtracted all values of t on by the value at 1365 20% of the embryo length, removing the "baseline" from the calculation of the integral. After computing the integral, the resulting value was then normalized by dividing by the distance along the embryo integrated over (see Section 2.3 for more detail).
The y-axis in our state space exploration plane is given by the position along 1370 the embryo at which the rate of RNAP loading reaches its midpoint (Fig. S8B) when fitted to a Hill function of Hill-coefficient six, the number of known Bi-coid binding sites on the hunchback P2 enhancer (Driever and Nusslein-Volhard, 1988a;Park et al., 2019). The maximum and minimum values of the RNAP loading rate are fixed to R max and R min (Section S1.3), motivated by the obser-1375 vation that the wild-type RNAP loading rates exhibit maximum and minimum plateaus at either end of the embryo (Fig. S6A). Note that due to the large space of possible model behaviors, many predicted rate profiles did not fit well to this Hill function (e.g. Fig. S8B). Instead, this process should be interpreted as a way to extract, for a particular model realization, the best-fit midpoint 1380 position of the observed wild-type RNAP loading rate profile. If a particular model is capable of recapitulating the experimental data, then the best-fit midpoint position under this protocol should agree with that of the data. Hereafter we refer to this quantity as the rate midpoint position.
Combined, the average t on delay and the rate midpoint position provide 1385 a simplified description of our data as well as of our theoretical predictions.
Each theoretical model inhabits a finite region in this two-dimensional state space, which we can calculate by systematically varying model parameters. Due to the large number of parameters of each model explored, the corresponding state-space boundaries were generated by efficient sampling of the underlying high-dimensional parameter space (Fig. S8C). The methodology is similar to the one described in Estrada et al. (2016). Briefly, a set of 100 ini-1395 tial points was generated, each with a randomized set of initial parameters, the specifics of which depended on the model being tested (Fig. S8C, i). The state space was sectioned into 10 horizontal and vertical slices (Fig. S8C, ii).
The most extremal points in each slice were found, resulting in two extremal anterior-posterior positions for the mid-point of the rate of RNAP loading and 1400 two extremal average t on delay points (Fig. S8C, iii). For each of these points, a new set of five points was generated using random parameters within a small neighborhood of the seed points determined by the extremal points of the pre-vious iteration (Fig. S8C, iv). These new points were plotted; some of these points may be more extreme than the previous set of points. Steps ii-iv were 1405 iterated, resulting in a growing boundary over time (Fig. S8C, v). Constraints imposed by the data were used to filter unrealistic results. First, if the simulated t on at the most anterior position (20% of the embryo length) was greater than 5.5 min or less than 1.5 min, the point was filtered out. This removal was justified experimentally, since none of the observed transcriptional onset values in 1410 this position in either experiment lay outside of these bounds ( Fig. 4D and G).
Second, if the simulated rate midpoint position was smaller than 16% or larger than 100% of the embryo length, then this point was also excluded. This exclusion was also justified experimentally, since the rate midpoints values of both experiments were within this interval ( Fig. 4C and F). Points that fulfilled these 1415 constraints were retained for the next iteration of the algorithm. This process was repeated until the resulting space of points no longer grew appreciably, resulting in an estimate of the size and shape of the state space for each of the models presented in Sections S1.2, S6.1, S7.1, and S8.1.
To determine whether the algorithm had indeed converged, the total area of 1420 each model's region in parameter space was tracked with each iteration number.
If the algorithm worked well, then this area would approach some saturation value. Fig. S8D shows the area of the state space corresponding to each model (normalized by the final area at the final iteration number) as a function of the iteration number. Each model converged to a finite value, indicating that the 1425 parameter space occupied by the models had been thoroughly explored.   The partition function in this generalized thermodynamic model is given by the polynomial

S5.2. State space exploration with the thermodynamic MWC model
where p inacc is the weight of the inaccessible state and P r,n is the weight of the accessible state with r RNAP molecules bound and n Bicoid molecules bound.
The overall transcriptional initiation rate is now where P 1,n is the statistical weight of each RNAP-bound state and R is the corresponding rate of transcriptional initiation. Note that, as described above, R is still equal to R max , the constraint described in Section S1.3, but we no 1480 longer use the R min constraint.
The resulting rate of transcriptional initiation is integrated over time to produce a simulated MS2 fluorescence trace using the same procedure as for the models presented in Sections S1.2, S7.1, and S8.1 (see Section S2.2 for details).
As with the thermodynamic MWC model, we allow for a mitotic repression time 1485 window to account for the lack of transcription early in the nuclear cycle.

S6.2. Generalized thermodynamic model state space exploration
Due to the high-dimensional parameter space of the generalized thermodynamic model, constraints were necessary to efficiently explore this parameter space (Section S5.1). These constraints were placed on the values of the in-1490 dividual microstate weights P r,n , based on dimensional analysis and heuristic arguments. Specifically, each weight P r,n is derived from a product of bind- which has bounds 0 < P 1,1 < (1000) 2 (1000) = 10 9 (S20) and thus provide a bound for the possible values that the weight P 1,1 can take.
In general, this process can be applied to enforce bounds on any microstate weight P r,n through constraining the possible values of p, b, ω bp , and ω b . As a result, the weight of a microstate with more Bicoid bound (i.e. higher values of n) will have a more generous dynamic range, due to the larger powers of b and 1505 ω b . In this way, exploration of parameter space can be made more constrained by restricting the possible values of the microstate weights P r,n .
The generalized thermodynamic model of transcription encompasses a larger area of the explored state space than the thermodynamic MWC model (Fig. S9A). However, as evident by the projection onto the x-axis, this model 1510 fails to capture the delays observed in zelda − data (Fig. S9B, yellow rectangle).

S6.3. Investigation of the failure of thermodynamic models
Here, we provide an intuitive explanation for why equilibrium thermodynamic models fail to recapitulate the delay in t on for zelda − embryos. The combination of the occupancy hypothesis and the assumption of separation of 1515 times scales described in Section S6.4 imply that the rate of transcriptional initiation at any moment in time is an instantaneous readout of the Bicoid concentration at that time point. Thus, any thermodynamic model is memoryless.
Intuitively, this means that a thermodynamic model requires transcription to begin as soon as the Bicoid concentration crosses a certain "threshold" since 1520 time delays between input and output require some sense of memory. Examination of the dynamic measurements of MS2 output in zelda − embryos reveals that no matter what "threshold" concentration of Bicoid is assigned for the start of transcription, the model cannot simultaneously describe two values of t on corresponding to different positions along the anterior-posterior axis 1525 ( Fig. S10A and B). Another self-consistency check of a thermodynamic model is to examine the concentration of Bicoid at t on for various positions along the embryo. Due to the memoryless nature of thermal equilibrium, a valid thermodynamic model predicts that, at different positions along the embryo, t on will occur when Bicoid 1530 reaches the same threshold value. For the zelda − data, however, the level of Bicoid at each anterior-posterior position's t on value actually decreases with increasing t on , suggesting the failure of the thermodynamic model (Fig. S10C).
Thus, the strong position-dependent delay in t on for the zelda − data cannot be explained by an instantaneous Bicoid readout mechanism.

S6.4. Re-examining thermodynamic models of transcriptional regulation
Thermodynamic models based on equilibrium statistical mechanics can be seen as limiting cases of more general kinetic models. For example, consider simple activation, where an activator whose concentration is modulated in time regulates transcription by binding to a single site (Fig. S11). In this generic 1540 model, the presence of activator can modulate the rates of activator and RNAP binding and unbinding through the parameters α, β, γ, and δ.
In order to reduce kinetic models to thermodynamic models where the probabilities of each state are dictated by Boltzmann weights such as those in Fig. 2, four conditions must be fulfilled. First, the rate of mRNA production must be linearly related to the probability of finding RNAP bound to the promoter ( Fig. S11i). This occupancy hypothesis is necessary for Eq. S2 to hold. Second, the time scales of binding and unbinding of RNAP and transcription factors must be much faster than the time scales of the concentration dynamics of these proteins (Fig. S11ii). Third, these time scales must also be much faster than the rate of transcriptional initiation and mRNA production (Fig. S11iii).
Under these conditions of separation of time scales, the binding and unbinding of proteins quickly reaches steady state while the overall concentrations of these molecular players are modulated (Segel and Slemrod, 1989). Fourth, there must be no energy input into the system (Fig. S11iv). This condition demands "detailed balance" (Vilar and Leibler, 2003;Ahsendorf et al., 2014;Hill, 1985): the product of state transition rates in the clockwise direction over a closed loop is equal to the product going in the counterclockwise direction, a constraint known as the cycle condition (Estrada et al., 2016). In the case of Fig. S11, this requirement implies that  Figure S11: A simple kinetic model of transcriptional activation in which activator molecules influence RNAP binding kinetics. The assumptions that make it possible to turn this kinetic model into a thermodynamic one are (i) the occupancy hypothesis, (ii, iii) a separation of time scales between binding and unbinding rates, and activator and mRNA production dynamics, respectively, and (iv) no energy expenditure (detailed balance).

S7.1. Non-equilibrium MWC model
The non-equilibrium MWC model is an extension of the thermodynamic MWC model presented in Section S1.2, where we now relax the assumption of separation of time scales (Fig. S11 ii and iii) and make it possible to assume, 1550 for example, that the system responds instantaneously to changes in activator concentration. Here, we explicitly simulate the full system of ordinary differential equations (ODEs) that describe the dynamics of the system out of steady state. Additionally, we allow for energy to be expended and thus do not enforce detailed balance through the cycle condition (Fig. S11iv). We still employ a 1555 mitotic repression window term, before which no transcription is allowed.
We consider a generic model with n Bicoid binding sites, and again ignore Zelda since we are only interested in recapitulating the zelda − mutant data. As a result, this new model has n + 1 total binding sites which, together with the closed chromatin state, results in a total of 2 n+1 + 1 = N microstates. In the 1560 case of six Bicoid binding sites, this results in N = 129 total microstates. We assign each microstate x i a label i and describe the transition rate from state j to state i using k ij , where i, j range from 0 to N − 1, inclusive.
In matrix notation, we write the system of ODEs as where X is a vector containing the fractional occupancy of each microstate x i and K is a matrix containing all the transition rates k ij . Normalizing such that 1565 the sum of all the components in the vector X is unity, we now have a vector representing the instantaneous probability of being in each microstate.
Here, bin() indicates taking the base 2 value of the binary label in the parentheses. The closed chromatin state is added manually and assigned to the last position in our binary label, x N −1 . This convention allows us to intuitively de-1575 fine each unique label for the system's microstates and provides a way to map the physical contents of a microstate with its associated label i.
In general, the overall transition matrix K can be very complex. However, we benefit from the fact that the only non-zero transitions k ij are the ones that correspond to physical processes: modifying the open/closed chromatin 1580 state, and binding and unbinding of Bicoid or RNAP molecules. In this binary notation, these constraints imply that the only nonzero transitions are the ones that represent individual flips between 0 and 1, as well as between the open and closed states 0 and N − 1. The transition matrix K is then easier to write, since it is clear from the binary representation which transitions must be nonzero.

1585
Finally, diagonal elements k ii are entirely constrained because they represent probability loss from a particular state i, and must be equal to the negative of the rest of the column i, such that the sum over each column in K is zero.
Given that the Bicoid concentration changes as a function of time and that we assume first-order binding kinetics, whichever rates k ij correspond to Bicoid To model transcription specifically, we assumed that at the beginning of the 1600 nuclear cycle, the system is in the closed chromatin state: x i (t = 0) = 0 except for the closed chromatin state x N −1 (t = 0) = 1. We simulated the full trajectory of all the microstates x i over time by solving the system of ODEs given in Eq. S22. Finally, we calculated p bound by summing the x i 's that correspond to RNAP-bound states, and then computed the subsequent transcriptional ini-1605 tiation rate by multiplying p bound with the transcription rate R. Here, R is the same R max as in Sections S1.2 and S6.1 but again we do not constrain the model using R min , just as in Section S6.1. Fig. S12A shows an example of this model for a system with only one Bicoid binding site and no closed chromatin state, for simplicity, resulting in a fourstate network. The binary indexing labels (shown beneath each state in light pink) can be converted into the base-10 labels (light teal) ranging from 0 to 3.
The connection matrix for this system is and the corresponding transition rate matrix K is where, in this example, k 02 represents the transition rate from state j to state i.
The diagonal elements k ii are equal to the negative of the sum of the elements 1610 in the rest of the column in order to preserve conservation of probability. For example, k 00 = −(k 10 + k 20 + k 30 ).
With all this information in hand, we solve for the occupancy of each of the four states using the matrix ODE k 00 k 01 k 02 0 k 10 k 11 0 k 13 k 20 0 k 22 k 23 0 k 31 k 32 k 33 In this case, the occupancy hypothesis relates p bound to the overall transcription rate, resulting in This model can produce time-dependent behavior not found in the thermodynamic models. Fig. S12B contains an example of a hypothetical input Bicoid activator concentration that switches instantaneously from zero to a finite value.

1615
In the thermodynamic models, the predicted transcriptional initiation rate also responds instantaneously (Fig. S12B, top). In contrast, for a suitable set of parameters, the non-equilibrium MWC model predicts a slow response over time (Fig. S12B, bottom).
To produce a simulated MS2 fluorescence trace, the resulting rate of mRNA 1620 production is integrated over time using the same procedure (Section S2.2) as the models presented in Sections S1.2, S6.1, and S8.1. As with the thermodynamic MWC model, we allow for a time window of mitotic repression to account for the lack of transcription early in the nuclear cycle. into a base-10 label for each state (light teal). The transition rates k ij are defined as the transition rate from state i into state j using this labeling system. (B) For an example input activator concentration temporal profile that is a step function, the time-dependent response is compared for the cases of separation of time scales and lack thereof. In the former, the transcriptional initiation rate responds instantaneously to the increase in activator input, while the response is slower in the latter.

1625
In the parameter exploration of this model (Section S5.1), the transition rates k ij were constrained with minimum and maximum values of k min = 1 and k max = 10 5 respectively, in units of inverse minutes. These bounds were conservatively chosen using the following estimates. First, we estimate the values of the possible unbinding rates k of f . We assume that RNAP and Bicoid obey the same unbinding kinetics. Estimates of in vivo single-molecule binding kinetics inferred from Mir et al. (2018) indicate that the lifetime of Bicoid on DNA is on the order of 3 s −1 . Second, we estimate the values of the possible on-rates k on using the classic Berg-Purcell equation for the case of a diffusion-limited binding to a perfectly absorbing spherical receptor (Berg and Purcell, 1977). In this case, the on-rate of molecule binding is given by where D is the diffusion coefficient of the molecule, a is the estimated size of the spherical receptor, and c 0 is the background concentration of the molecular species. Since here we are talking about transcription factor binding to a Bicoid binding site, we assume a to be on the order of 5 nm. We assume that RNAP and Thus, our maximum and minimum transition rate bounds of k min = 1 min −1 and k max = 10 5 min −1 lie outside these estimated binding and unbinding rates.
One caveat of the state-space exploration approach is that the high dimensionality of the non-equilibrium MWC model prevented us from calculating the full state-space boundary using six Bicoid binding sites. Due to computational 1630 costs, we were only able to accurately produce a state-space boundary for this model (Section S7.1) using five Bicoid binding sites. Running the exploration for a model with six Bicoid binding sites took over two weeks on our own server, and the algorithm had not noticeably converged in the end.
The results of the state space exploration for the non-equilibrium MWC 1635 model using five Bicoid binding sites resulted in larger average t on delays than the thermodynamic models (Sections S1.2 and S6.1). However, this model, like those, failed to reproduce the delays observed in zelda − data (Fig. S9B, gray bar).
S8. Transcription factor-driven model of chromatin accessibility 1640 S8.1. Transcription factor-driven model of chromatin accessibility The transcription factor-driven model of chromatin accessibility is a slight modification of the thermodynamic MWC model (Section S1.2) that replaces the MWC mechanism of chromatin transitions with a direct driving action due to Bicoid and Zelda. Here, we retain the idea of inaccessible vs. accessible 1645 states, but no longer demand that these states be in thermodynamic equilibrium. Instead, the system begins in the inaccessible state and undergoes a series of m identical, slow, and effectively irreversible transitions to the accessible state. Once these transitions into the accessible state occur, the system can rapidly and reversibly transition into all of its accessible microstates such 1650 that the probability of the system being in any of these microstates is described by thermodynamic equilibrium. The accessible states are governed by the same rules and parameters as the thermodynamic MWC model (Section S1.2), albeit without the ∆ε chrom parameter since now the transition from the inaccessible to accessible state is unidirectional.

1655
Additionally, we consider two possible contributions for these irreversible transitions: a Bicoid-dependent pathway and a Zelda-dependent pathway. We assume the transition rates to be first-order in Bicoid and Zelda, respectively, such that and Here, π b is the Bicoid-dependent contribution to the transition rates and π z is the corresponding Zelda-dependent contribution. There are two input parameters c b and c z that give the relative speed of each transition rate contribution.
The overall rate π of each irreversible transition is given by the sum Calculating the overall RNAP loading rate then simply corresponds to rescaling p bound with the overall probability p m (t) of being in the accessible state: where R is the same maximum rate used in Section S1.2. Note that p m (t) is a time-dependent quantity that changes over time. To calculate p m (t), we solve the corresponding system of ODEs that describes the time evolution of P d P dt = Π P , where Π is the transition rate matrix describing the time evolution of the system. Π, by definition, is a square matrix with dimension m + 1. Given the initial condition that the system begins in the inaccessible state the system of ODEs can be solved to find the probability of being in the ac- where π is given by Eq. S31.
For simplicity, the time evolution of P was solved using MATLAB's ode15s solver.
With the probability p m (t) of the system being in the accessible state calculated, we now calculate the probability p bound of RNAP bound to the promoter in the accessible states, which lie in thermodynamic equilibrium with each other.
Because we now only have accessible states, the partition function is where z, p, and b correspond to the non-dimensionalized concentrations of Zelda, RNAP, and Bicoid, respectively, and ω b and ω bp are the cooperativities between Bicoid molecules and between Bicoid and RNAP, respectively. Thus, the overall transcriptional initiation rate is given by Due to the lack of the inaccessible state in the partition function and because we assume that Zelda does not directly interact with Bicoid or RNAP, now the 1660 presence of Zelda mathematically separates out so that only Bicoid influences transcription. The calculation above is a standard equilibrium statistical mechanical calculation, except that we have weighted the final result with p m (t), the probability of being in the accessible states. The resulting rate is integrated to produce a simulated MS2 fluorescence trace using the same procedure 1665 (Section S2.2) as the models presented in Sections S1.2, S6.1, and S7.1.
Interestingly, we found that a mitotic repression term was not necessary to recapitulate the data, since the presence of intermediary states produced the necessary delay to explain the experimentally observed t on values in the data ( Fig. 4D and G, points).

1670
In order to sufficiently explain the data, we found that a minimum of m = 5 irreversible steps was necessary. Fig. S13A

S8.2. Transcription factor-driven model of chromatin accessibility state space exploration
In the parameter exploration of this model (Section S5.1), the parameters were constrained as The parameters shared with the thermodynamic MWC model retained the constraints described in Section S1.3.