Research Article

Physics of Living Systems

Predictive modeling reveals that higher-order cooperativity drives transcriptional repression in a synthetic developmental enhancer

Chan Zuckerberg Biohub, United States
Department of Chemical Biology, University of California, Berkeley, United States
Department of Physics, University of California, Berkeley, United States
Department of Biology, Ecole Polytechnique, France
Biophysics Graduate Group, University of California, Berkeley, United States
Department of Biology, New York University, United States
Department of Molecular and Cell Biology, University of California, Berkeley, United States
Institute for Quantitative Biosciences–QB3, University of California at Berkeley, United States

Dec 12, 2022

https://doi.org/10.7554/eLife.73395

Version of Record

Accepted for publication after peer review and revision.

Download
Cite
Share
CommentOpen annotations (there are currently 0 annotations on this page).

Version of Record published: January 12, 2023 (This version)
Accepted Manuscript published: December 12, 2022 (Go to version)
Accepted: December 9, 2022
Received: August 27, 2021
Preprint posted: July 29, 2021 (Go to version)

1. Of interest
The chemorepellent, SLIT2, bolsters innate immunity against Staphylococcus aureus

Vikrant K Bhosle, Chunxiang Sun ... Lisa A Robinson

Research Article Sep 29, 2023
Further reading

Abstract
Editor's evaluation
Introduction
Results
Discussion
Materials and methods
Appendix 1
Data availability
References
Decision letter
Author response
Article and author information
Metrics

Abstract

A challenge in quantitative biology is to predict output patterns of gene expression from knowledge of input transcription factor patterns and from the arrangement of binding sites for these transcription factors on regulatory DNA. We tested whether widespread thermodynamic models could be used to infer parameters describing simple regulatory architectures that inform parameter-free predictions of more complex enhancers in the context of transcriptional repression by Runt in the early fruit fly embryo. By modulating the number and placement of Runt binding sites within an enhancer, and quantifying the resulting transcriptional activity using live imaging, we discovered that thermodynamic models call for higher-order cooperativity between multiple molecular players. This higher-order cooperativity captures the combinatorial complexity underlying eukaryotic transcriptional regulation and cannot be determined from simpler regulatory architectures, highlighting the challenges in reaching a predictive understanding of transcriptional regulation in eukaryotes and calling for approaches that quantitatively dissect their molecular nature.

Editor's evaluation

The work by Kim et al., used synthetic constructs in Drosophila to examine the relationship between regulators and transcription initiation. By measuring regulator concentrations and the corresponding RNA polymerase initiation rates in different synthetic constructs and using a thermodynamic model, the authors concluded that higher-order cooperativities between the repressor on adjacent binding sites, and that between the repressor and RNA polymerase are needed to explain the observed response curves in RNA polymerase loading rate. This work targets a challenging question in eukaryotic transcription regulation, where higher-order cooperativity between different molecular components, in addition to simple transcription factor binding and unbinding, is often necessary to account for observed promoter behaviors when multiple elements (repressors, mediators, activators) exist.

https://doi.org/10.7554/eLife.73395.sa0

Introduction

During embryonic development, transcription factors bind stretches of regulatory DNA termed enhancers to dictate the spatiotemporal dynamics of gene expression patterns that will lay out the future body plan of multicellular organisms (Spitz and Furlong, 2012; Small and Arnosti, 2020). One of the greatest challenges in quantitative developmental biology is to predict these patterns from knowledge of the number, placement, and affinity of transcription factor binding sites within enhancers. The early embryo of the fruit fly Drosophila melanogaster has become one of the main workhorses in this attempt to achieve a predictive understanding of cellular decision-making in development due to its well-characterized gene regulatory network and transcription factor binding motifs, and the ease with which its development can be quantified using live imaging (Garcia et al., 2020; Small and Arnosti, 2020; Rivera et al., 2019).

Predictive understanding calls for the derivation of theoretical models that generate quantitative and experimentally testable predictions. Thermodynamic models based on equilibrium statistical mechanics have emerged as a widespread theoretical framework to achieve this goal (Ackers et al., 1982; Vilar and Leibler, 2003; Bolouri and Davidson, 2003; Bintu et al., 2005b; Bintu et al., 2005a; Segal et al., 2008; Fakhouri et al., 2010; Sayal et al., 2016; Phillips et al., 2019; Eck et al., 2020). For instance, over the last decade, a dialogue between these thermodynamic models and experiments demonstrated the capacity to quantitatively predict bacterial transcriptional regulation from knowledge of the DNA regulatory architecture (He et al., 2010; Garcia and Phillips, 2011; Brewster et al., 2014; Garcia et al., 2012; Sepúlveda et al., 2016).

The predictive power of these models is evident when inferring model parameters from simple regulatory architectures and using those parameters to make parameter-free predictions of more complex architectures (Boedicker et al., 2013a; Boedicker et al., 2013b, Razo-Mejia et al., 2018; Phillips et al., 2019). Consider, for example, that RNA polymerase II (RNAP)—which we take as a proxy for the whole basal transcriptional machinery—binds to a promoter with a dissociation constant $K_{p}$ . When RNAP is bound, transcription is initiated at a rate $R$ (Figure 1A). In the absence of any regulation, a thermodynamic model will only have $K_{p}$ and $R$ as its free parameters which can be experimentally determined by, for example, measuring mRNA distributions (Razo-Mejia et al., 2020). Now, we assume that the parameters $K_{p}$ and $R$ inferred in this step do not just enable a fit to the data, but that their values represent physical quantities that remain unaltered as more complex regulatory architectures are iteratively considered. As a result, when we consider the case where a single repressor molecule can bind, our model calls for only two new free parameters: a dissociation constant for repressor to its binding motif $K_{r}$ , and a negative cooperativity between repressor and RNAP, $ω_{r p}$ , that makes the recruitment of RNAP to the DNA less favorable when the repressor is bound to its binding site (Figure 1B). Once again, after determining $K_{r}$ and $ω_{r p}$ experimentally (Phillips et al., 2019), we consider the case where two repressors can bind simultaneously (Figure 1C). If the repressors interact with RNAP independently of each other, then our model has no remaining free parameters such that we will have reached complete predictive power. However, protein-protein interactions between repressors could exist or even higher-order interactions giving rise to a repressor-repressor-RNAP ternary complex might be present. This extra complexity would require yet another round of experimentation to quantify these interactions represented by $ω_{r r}$ and $ω_{r r p}$ in Figure 1C, respectively. Even after quantifying these parameters, predictive power might not be reached if, after adding yet another repressor binding site, a complex between all three repressors and RNAP can be formed (Figure 1D).

Figure 1

Download asset Open asset

Building up predictive models of transcriptional repression.

(A) In the absence of repressor binding, gene expression can be characterized by a dissociation constant between RNAP and the promoter $K_{p}$ and the rate of transcription initiation when the promoter is bound by RNAP $R$ . (B) In the presence of a single repressor binding site, models need to account for two additional parameters describing the repressor dissociation constant $K_{r}$ and a repressor-RNAP interaction term $ω_{r p}$ . (C) For two-repressor architectures, parameters accounting for repressor-repressor interactions $ω_{r r}$ and for interactions giving rise to a repressor-repressor-RNAP complex could also have to be incorporated. (D) For the case of three repressor binding sites, additional parameters $ω_{r r r}$ and $ω_{r r r p}$ capturing the higher-order cooperativity between three repressor molecules and between three Runt molecules and RNAP, respectively, could be necessary. Note the nomenclature shown below each construct, which indicates which Runt binding sites are present in each construct.

While protein-protein cooperativity captured by $ω_{r r}$ has been studied both in bacteria (Ackers et al., 1982; Ptashne and Gann, 2002) and eukaryotes (Giniger and Ptashne, 1988; Ma et al., 1996; Lebrecht et al., 2005; Parker et al., 2011; Fakhouri et al., 2010; Sayal et al., 2016), the necessity of accounting for higher-order interactions such as those described in our example by the $ω_{r r p}$ and $ω_{r r r p}$ terms had only been demonstrated in archeae (Peeters et al., 2013) and bacteria (Dodd et al., 2004). The need to invoke this higher-order cooperativity in eukaryotes only became apparent in the last few years (Estrada et al., 2016b; Park et al., 2019; Biddle et al., 2020). These higher-order cooperativities might be necessary in order to account for the complex interactions mediated by, for example, the recruitment of co-repressors (Courey and Jia, 2001; Walrad et al., 2011), mediator complex (Park et al., 2019), or any other element of the transcriptional machinery. As a result, while posing a challenge to reaching a parameter-free predictive understanding of transcriptional regulation, higher-order cooperativity provides an avenue for quantifying the complexity of the molecular processes underlying eukaryotic cellular decision-making.

In this paper, we sought to test whether an iterative and predictive approach, such as that outlined in Figure 1, was possible for transcriptional repression in the early embryo of the fruit fly Drosophila melanogaster or whether it is necessary to invoke higher-order cooperativities that challenge the reach of our predictive models as we add more complexity to the system. To make this possible, we engineered binding sites for the Runt repressor into the Bicoid-activated hunchback P2 minimal enhancer. We systematically varied the number and placement of Runt binding sites within this enhancer (Chen et al., 2012) in order to determine whether model fits to real-time transcriptional measurements from the enhancer constructs containing only one-Runt binding site could accurately predict repression in two- and three-Runt binding site constructs (Figure 1). We found that a thermodynamic model can recapitulate all our data. However, we also discovered that, while the model could describe repression by a single Runt repressor, protein-protein and higher-order cooperativities had to be invoked in order to quantitatively account for regulation by two or more repressor molecules. While these higher-order cooperativities limit the iterative bottom-up discourse between theory and experiment that has been successful in bacteria (Phillips et al., 2009), they also provide a concrete theoretical framework for quantifying the complexities behind eukaryotic transcriptional control, and call for the development of new theories and experiments specifically conceived to uncover the the molecular underpinnings of this complexity.

Results

Predicting transcription rate using a thermodynamic model of Bicoid activation and Runt repression

Inspired by the theory-experiment dialogue leading to predictive understanding of the lac operon in E. coli over the last four decades (Phillips et al., 2019; Razo-Mejia et al., 2018; Garcia and Phillips, 2011; Garcia et al., 2012; Ackers et al., 1982; Buchler et al., 2003), we built a predictive model of Runt repression on the Bicoid-activated hunchback P2 enhancer using the thermodynamic model framework (Phillips et al., 2019; Bintu et al., 2005b; Bintu et al., 2005a) with the goal of predicting the rate of transcription initiation as a function of input transcription factor concentration, and the number and placement of Runt repressor binding sites. Our model rests on the ‘occupancy hypothesis’ that states that the rate of mRNA production, $d [m R N A] / d t$ , is proportional to the probability of the promoter being bound by RNA polymerase II (RNAP), $p_{b o u n d}$ , such that

\frac{d [m R N A]}{d t} = R p_{b o u n d},

where $R$ is the rate of mRNA production when the promoter is occupied by RNAP. Note that, throughout this study, we treat the rate of transcription initiation and the rate of RNAP loading interchangeably.

To generate intuition, we start by modeling the case of hunchback P2 with one Runt binding site. Figure 2A illustrates the possible states the system can be found in. Each state has an associated statistical weight which can be calculated as prescribed by equilibrium statistical mechanics (Bintu et al., 2005b; Bintu et al., 2005a). Here, we assume that there are six Bicoid binding sites with the same dissociation constant given by $K_{b}$ , one Runt binding site with a dissociation constant specified by $K_{r}$ , and a promoter with a dissociation constant for RNAP prescribed by $K_{p}$ . In the absence of Runt, we consider four states as shown in the top two rows of Figure 2A. Here, we assume that Bicoid-Bicoid cooperativity is so strong that the enhancer can either be unoccupied or completely bound by Bicoid molecules (Gregor et al., 2007; Park et al., 2019). Further, we consider an interaction between Bicoid and RNAP given by $ω_{b p}$ . For simplicity, we use the dimensionless parameters $b = [B i c o i d] / K_{b}$ , $r = [R u n t] / K_{r}$ and $p = [R N A P] / K_{p}$ . These assumptions lead to a functional form reminiscent of a Hill function that explains the sharp step-like expression pattern along the embryo’s anterior-posterior axis of the hunchback gene (Gregor et al., 2007; Park et al., 2019; Driever and Nüsslein-Volhard, 1988; Driever et al., 1989). A full thermodynamic model in which we do not make this assumption of high Bicoid-Bicoid cooperativity is discussed in detail in Section ‘Derivation of the general thermodynamic model for the hunchback P2 enhancer’ and Section ‘Derivation of the general and simpler thermodynamic model for the hunchback P2 enhancer with one Runt binding site’.

Figure 2 with 2 supplements see all

Download asset Open asset

Thermodynamic model of transcriptional regulation by Bicoid activator and Runt repressor.

(A) States and statistical weights for the regulation of *hunchback* P2 with one Runt binding site in the limit of strong Bicoid-Bicoid cooperativity. Here, we use the dimensionless parameters $b = [B i c o i d] / K_{b}$ , $r = [R u n t] / K_{r}$ , and $p = [R N A P] / K_{p}$ , where $K_{b}$ , $K_{r}$ , and $K_{p}$ are the dissociation constants of Bicoid, Runt, and RNAP, respectively. $ω_{b p}$ represents the cooperativity between Bicoid and RNAP, $ω_{r p}$ captures the cooperativity between Runt and RNAP, and $R$ represents the rate of transcription when the promoter is occupied by RNAP. The top two rows correspond to states where only Bicoid and RNAP act, while the bottom two rows represent repression by Runt. (B) Representative prediction of RNAP loading rate as a function of Bicoid and Runt concentrations for $ω_{b p} = 3, ω_{r p} = 0.001, p = 0.001, R = 1 (A U / m i n)$ .

The molecular mechanism by which Runt downregulates transcription of its target genes remains unclear (Chen et al., 2012; Hang and Gergen, 2017; Koromila and Stathopoulos, 2017; Koromila and Stathopoulos, 2019). Here, we assume the so-called ‘direct repression’ model (Gray et al., 1994) that posits that Runt operates by inhibiting RNAP binding to the promoter through a direct Runt-RNAP interaction term given by $ω_{r p} < 1$ independently of Bicoid. As a result, in the presence of Runt, we consider four additional states as shown in the bottom two rows of Figure 2A. Other potential mechanisms of Runt repression are further discussed in Supplementary Section ‘Comparison of different modes of repression’, where we also show that the choice of specific mechanism does not change our conclusions.

Given these assumptions, we arrive at the microstates and corresponding statistical weights shown in Figure 2A. The probability of finding RNAP bound to the promoter, $p_{b o u n d}$ , is calculated by dividing the sum of all statistical weights featuring RNAP by the sum of the weights of all possible microstates. The calculation of $p_{b o u n d}$ combined with Equation 1 leads to the expression

R a t e = R p_{b o u n d} = R \frac{p + b^{6} p ω_{b p} + r p ω_{r p} + b^{6} r p ω_{b p} ω_{r p}}{1 + b^{6} + r + b^{6} r + p + b^{6} p ω_{b p} + r p ω_{r p} + b^{6} r p ω_{b p} ω_{r p}},

which makes it possible to predict the output rate of mRNA production as a function of the input concentrations of Bicoid and Runt (Figure 2B). With this theoretical framework in hand, we experimentally tested the predictions of this model.

Measuring transcriptional input-output to test model predictions

The transcriptional input-output function in Figure 2B indicates that, in order to predict the rate of RNAP loading and to test our theoretical model, we need to first measure the concentration of the input Bicoid and Runt transcription factors. In order to quantify the concentration profile of Bicoid, we used an established eGFP-Bicoid line (Gregor et al., 2007) and measured mean Bicoid nuclear concentration dynamics along the anterior-posterior axis of the embryo over nuclear cycles 13 and 14 (nc13 and nc14, respectively) as shown in Movie Figure 3—video 1 (Eck et al., 2020). An example snapshot and time trace of Bicoid nuclear concentration dynamics at 40% of the embryo length appear in Figure 3A and B.

Figure 3 with 11 supplements see all

Download asset Open asset

Measurement of input transcription factor concentrations and output rate of transcription to test model predictions.

(A) Snapshot of an embryo expressing eGFP-Bicoid spanning 20–60% of the embryo length. (For a full time-lapse movie, see Movie Figure 3—video 1) (B) Bicoid nuclear fluorescence dynamics taken at 40% of the embryo. (C) Snapshot of an embryo expressing eGFP:LlamaTag-Runt spanning 20–60% of the embryo length. (For a full time-lapse movie, see Movie Figure 3—video 2) (D) Runt nuclear concentration dynamics in males and females. (E) Measured transcription factor concentration profiles along the anterior-posterior axis of the embryo. The concentration profiles are averaged over the gray shaded regions shown in (B) and (D) which corresponds to a time window between 5 and 10min into nc14. (F) Predicted RNAP loading rate for *hunchback* P2 with one Runt binding site over the anterior-posterior axis generated for a reasonable set of model parameters $K_{b} = 30$ AU, $K_{r} = 100$ AU, $ω_{b p} = 100$ , $p = 0.001$ , and $R = 1$ AU/min for varying values of the Runt-RNAP interaction term $ω_{r p} = [10^{- 2}, 1]$ . (G) Schematic of the MS2 system where 24 repeats of the MS2 loop sequence are inserted downstream of the promoter followed by the *lacZ* gene. The MS2 coat protein (MCP) fused to GFP binds the MS2 loops. (H) Example snapshot of an embryo expressing MCP-GFP and Histone-RFP. Green spots correspond to active transcriptional loci and red circles correspond to nuclei. Spot intensities are proportional to the number of actively transcribing RNAP molecules. (I) Representative MS2 fluorescence averaged over a narrow window (2.5% of the embryo length) along the anterior-posterior axis of the embryo. The initial rate of RNAP loading was obtained by fitting a line (brown) to the initial rise of the data and the x-intercept is defined as the onset of transcription ( $T_{O N}$ ). (J) Measured initial rate of RNAP loading (over a spatial bin of 2.5% of the embryo length) across the anterior-posterior axis of the embryo, from the *hunchback* P2 enhancer. (B, D, E, and J, error bars represent standard error of the mean over $\geq 3$ embryos; I, error bars represent standard error of the mean over the spatial averaging corresponding to roughly ten nuclei; A, C, and H, white scale bars represent 20 μm.).

Quantification of the Runt concentration using standard fluorescent protein fusions is not possible due to the slow maturation times of these proteins (Bothma et al., 2018). We therefore measured Runt concentration dynamics using our recently developed LlamaTags, which are devoid of such maturation dynamics artifacts (Bothma et al., 2018). Specifically, we generated a new fly line harboring a fusion of a LlamaTag against eGFP to the endogenous runt gene using CRISPR/Cas9-mediated homology-directed repair (Materials and Methods; Harrison et al., 2010, Gratz et al., 2015).

Using this LlamaTag fusion, we measured the mean Runt nuclear fluorescence along the anterior-posterior axis of the embryo over nc13 and nc14 (Materials and Methods; Figure 3B; Movie Figure 3—video 2). As expected due to the location of the runt gene on the X chromosome (Lott et al., 2011), there is a sex dependence in the nuclear concentration levels in nc13, with males displaying lower Runt levels than females; this difference is compensated by early nc14 (Figure 3C and D). As a result, for ease of analysis, we focused subsequent quantitative dissection on nc14.

We used the measured input protein concentration profiles to predict the output transcription rate. To make this possible, we invoked previous observations stating that the concentration dynamics of input transcription factors does not significantly affect the initial rate of RNAP loading (Garcia et al., 2013; Eck et al., 2020). As a result, we decided to use the time-averaged concentration dynamics of Bicoid and Runt over a time window spanning 5 min after the 13th anaphase to 10 min after this anaphase (gray shaded region in Figure 3B and D) as inputs to our model, resulting in the static spatial concentration profiles shown in Figure 3E. We then used these time-averaged concentration profiles of input transcription factors to calculate the time-averaged rate of transcription initiation over the same time window. In the Supplementary Information Section ‘Comparing using static versus dynamic transcription factor concentrations as model inputs’ we compare this methodology with one that acknowledges input transcription factor concentration dynamics and show that the prediction stemming from both approaches leads to equivalent theoretical predictions. Specifically, the time-averaged rate of transcription predicted by the dynamic inputs was similar to the rate of transcription predicted by the static inputs.

Along the anterior-posterior axis of the embryo, the measured Bicoid and Runt concentration profiles define a trajectory through the input-output function (Figure 2B). Given a set of parameters, this trajectory predicts the initial rate of RNAP loading. This quantitative prediction can be directly compared with experimentally measured transcription initiation rates. For example, given the concentration profiles shown in Figure 3E, we calculate the RNAP loading rate as a function of the position along the embryo for different values of the Runt-RNAP interaction, captured by $ω_{r p}$ . Figure 3F illustrates how $ω_{r p}$ shapes the predicted profiles for the RNAP loading rate. As expected, the prediction shows that the rate of transcription decreases as the strength of the Runt-RNAP interaction decreases.

Next, we sought to experimentally test these predictions by measuring the rate of RNAP loading using the MS2 system (Bertrand et al., 1998; Lucas et al., 2013; Garcia et al., 2013). Here, we inserted 24 repeats of the MS2 loop sequence following the hunchback P2 enhancer and even-skipped promoter in our reporter construct, which leads to the fluorescent labeling of sites of active transcription in living embryos (Figure 3G and H; Movie Figure 3—video 3). The fluorescence intensity of each MS2 spot is proportional to the number of actively transcribing RNAP molecules (Garcia et al., 2013). In order to quantify the transcriptional activity reported by MS2, we measured the mean MS2 spot fluorescence over nuclei in a narrow spatial window (Figure 3I; Garcia et al., 2013; Eck et al., 2020). To measure the initial rate of RNAP loading, we obtained the slope of the initial rise in the number of actively transcribing RNAP molecules over the same time window used to average input transcription factor concentration (Figure 3I, brown line). The resulting RNAP loading rate plotted over the anterior-posterior axis is in qualitative agreement with the classic pattern driven by the hunchback P2 minimal enhancer (Figure 3J; Garcia et al., 2013, Chen et al., 2012, Park et al., 2019).

While we chose the initial rate of transcription as the experimental measurable to confront against our model predictions, the MS2 technique can also report on other dynamical features of transcription such as the time window over which transcription occurs and the fraction of loci that engage in transcription at any point over the nuclear cycle. Although these two quantities have been shown to be relevant in shaping gene expression patterns in other regulatory contexts (Garcia et al., 2013; Lammers et al., 2020; Eck et al., 2020; Dufourt et al., 2018; Reimer et al., 2021), we found that the transcription time window was not significantly regulated in the presence of Runt (Figure 3—figure supplement 3). As described in Section ‘Quantitative interpretation of MS2 signals’, we did find some modulation of the fraction of transcriptionally engaged loci for a subset of our synthetic enhancer constructs but, as we could not detect a clear trend in how this fraction of active loci was modulated, we did not pursue a theoretical dissection of the control of this quantity by Runt.

Enhancer sequence dictates unrepressed transcription rates by determining RNAP-promoter interactions

With these theoretical models and our experimental platform in hand, we designed a set of synthetic enhancer constructs with differing number and placement of Runt binding sites as shown in Figure 4A (top) , and Figure 4—figure supplement 1. Our enhancer sequences are identical to those created and validated by Chen et al., 2012, which kept the length of the enhancer sequence consistent and inserted experimentally validated Runt binding sites (Melnikova et al., 1993; Lewis et al., 1999; Chen et al., 2012; Koromila and Stathopoulos, 2017) by mutating the base pairs within the enhancer that are not mapped to binding sites for any known transcripiton factor in the early fruit fly embryo (Hertz et al., 1990; Hertz and Stormo, 1999).

Figure 4 with 2 supplements see all

Download asset Open asset

Enchancer-to-enhancer variability in the unrepressed transcription level stems from unique RNAP-dependent parameters.

(A) Measured initial rates of RNAP loading across the anterior-posterior axis of the embryo for all synthetic enhancer constructs in the *absence* of Runt protein. (The [111] synthetic enhancer construct with the position of Bicoid (red) and Runt (green) binding sites is shown in genomic length scale on top as a reference.) (B) Representative best MCMC fit and (C) associated corner plot for the [001] construct in the *runt* null background. (D) Inferred model parameters for all synthetic enhancers in the absence of Runt repressor. Note the large spread in $ω_{b p}$ , consistent with the corner plot shown in (C), which indicates that our model does not constrain this parameter well compared to the other parameters. (E) Coefficient of variation of inferred parameters. (**A, B**), shaded regions represent the standard error of the mean over>3 embryos; (B) error bars from MCMC fit represent 95% confidence interval; (D) error bars represent standard deviations calculated from the MCMC posterior chains; (E) error bars are calculated by propagating the standard deviation of individual parameters from their MCMC chains.

A major assumption of our theoretical approach is that the model parameters obtained from simple regulatory architectures can be used as inputs for more complex constructs. For instance, we assume that the Runt-independent model parameters for Bicoid and RNAP action— $K_{b}$ , $ω_{b p}$ , $p$ and $R$ (Figure 2A)—are conserved for all constructs containing Runt binding sites regardless of their number and placement in the enhancer. If model parameters can be shared across constructs, then our model should predict the same profile for the rate of transcription across all synthetic enhancer constructs.

To test this assumption, we measured the initial rate of RNAP loading in all of our reporter constructs, in runt null embryos (Materials and Methods). Notably, unrepressed transcription rates varied significantly across synthetic enhancers (Figure 4A). For example, despite no Runt being present, the [001] construct had almost twice the unrepressed rate of [000].

This large construct-to-construct variability in unrepressed transcription rates likely originates from the Runt binding site sequences interfering with some combination of Bicoid and RNAP function. To uncover the mechanistic effect of these Runt binding sites sequences on unrepressed activity, we sought to determine which parameters in our thermodynamic model varied across constructs. In the absence of Runt repressor, only four states remain corresponding to the two top rows of Figure 2A. In this limit, the predicted rate of transcription is given by

R a t e = R \frac{p + {(\frac{[B i c o i d]}{K_{b}})}^{6} p ω_{b p}}{1 + p + {(\frac{[B i c o i d]}{K_{b}})}^{6} + {(\frac{[B i c o i d]}{K_{b}})}^{6} p ω_{b p}},

where we have invoked the same parameters as in Figure 2 and Equation 2. For clarity, the free parameters in this equation are marked using the red color.

To obtain the model parameters for each construct measured in Figure 4A, we invoked the Bayesian inference technique of Markov Chain Monte Carlo (MCMC) sampling that has been widely used for inferring the biophysical parameters from theoretical models (Liu et al., 2021, Razo-Mejia et al., 2018, Geyer and Thompson, 1992; Supplementary Section ‘Markov Chain Monte Carlo inference protocol’). A representative comparison of the MCMC fit to the experimental data reveals a good agreement between theory and experiment (Figure 4B). MCMC sampling also gives the distribution of the posterior probability for each parameter as well as their cross-correlation (Figure 4C). These corner plots reveal relatively unimodal posterior distributions, suggesting that a unique set of parameters can explain the data.

Note that, while the Bicoid dissociation constant $K_{b}$ and the Bicoid-RNAP interaction term $ω_{b p}$ remain largely unchanged regardless of enhancer sequence, there is considerable variability in the inferred mean RNAP-dependent parameters $p$ and $R$ (Figure 4D). This variability can be further quantified by examining the coefficient of variation,

C V = \frac{σ}{μ},

where $σ$ and μ are the standard deviation and the mean of each parameter, respectively, calculated over all constructs. The coefficients of variation for the RNAP and promoter-dependent parameters are much higher than those for Bicoid-dependent parameters (≈ 40% versus < 10%; Figure 4E). This suggests that the variability in unrepressed transcription rates due to the presence of Runt binding sites stems from differences in the behavior of RNAP at the promoter rather than differences in Bicoid binding or activation. As a result, as we consider increasingly more complex regulatory architectures, we associated each construct with its own specific Bicoid- and RNAP-dependent parameters as inferred in Figure 4D. In contrast, as we will show below, we will conserve Runt-dependent parameters as we consider increasingly more complex constructs featuring more Runt binding sites.

The thermodynamic model recapitulates repression by one Runt binding site

Next, we asked whether our model recapitulates gene expression for the hunchback P2 enhancer with a one-Runt binding site in the presence of Runt repressor as predicted by Equation 2. We posited that, since the binding site sequence remains unaltered throughout our constructs (Figure 4—figure supplement 1), the value of the Runt dissociation constant $K_{r}$ would also remain unchanged across these enhancers regardless of Runt binding site position; however, we assumed that, as the distance between Runt and the promoter varied, so could the Runt-RNAP interaction term $ω_{r p}$ .

We measured the initial rate of transcription along the embryo for all our constructs containing one Runt binding site in the presence of Runt protein. In this case of a single Runt binding site, Equation 2 predicts that the initial rate of RNAP loading will be given by

R a t e = R p_{b o u n d} = R \frac{p + b^{6} p ω_{b p} + \frac{[R u n t]}{K_{r}} p ω_{r p} + b^{6} \frac{[R u n t]}{K_{r}} p ω_{b p} ω_{r p}}{1 + b^{6} + \frac{[R u n t]}{K_{r}} + b^{6} r + p + b^{6} p ω_{b p} + \frac{[R u n t]}{K_{r}} p ω_{r p} + b^{6} \frac{[R u n t]}{K_{r}} p ω_{b p} ω_{r p}} .

Here, we have have rewritten Equation 2 to clarify which parameters are fixed and which parameters are inferred using color coding. Specifically, we took Runt-independent parameters ( $K_{b}$ , $ω_{b p}$ , $p$ and $R$ ), shown in black, as given by the inference from our previous experiments in the absence of Runt (Figure 4). Further, Runt-dependent parameters ( $K_{r}$ and $ω_{r p}$ ) which we will infer, are shown in red. We then used MCMC sampling to infer these Runt-dependent parameters for each of our constructs while retaining the mean values of Runt-independent parameters.

The resulting MCMC fits show significant agreement with the experimental data (Figure 5A), confirming that, within our model, the same dissociation constant $K_{r}$ can be used for all Runt binding sites regardless of their position within the enhancer. Further, the corner plot yielded a unimodal distribution of posterior probability of the inferred parameters (Figure 5B), indicating the existence of a unique set of most-likely model parameters. We challenged our assumption of constant $K_{r}$ across our constructs in Section Figure 5—figure supplement 3, where we show that, even if we posit that each construct has a different Runt dissociation constant, the obtained $K_{r}$ values are comparable.

Figure 5 with 3 supplements see all

Download asset Open asset

Testing the direct repression model in the presence of one Runt binding site.

(A) Initial transcription rate as a function of position along the embryo for the three constructs containing one Runt binding site in the presence and absence of Runt repressor, together with their best MCMC fits. (B) Corner plots from MCMC inference for all constructs with one Runt binding site. (C) Inferred $ω_{r p}$ value as a function of distance between the promoter and the Runt binding site. (A, data points represent mean and standard error of the mean over the embryos and shaded error bars represent 95% confidence intervals for the best MCMC fits for Runt WT datasets; C, data and error bars represent the mean and standard deviation of the posterior chains, respectively.).

The observed trend in the Runt-RNAP interaction captured by $ω_{r p}$ qualitatively agrees with the “direct repression” model. Specifically, because the model assumes that Runt interacts directly with RNAP, it predicts that, the farther apart Runt and the promoter are, the lower this interaction should be (Gray et al., 1994). In agreement with this prediction, the mean value of $ω_{r p}$ obtained from our fits changes from high repression ( $ω_{r p} \approx 0.1$ ) in the [001] construct to almost no repression ( $ω_{r p} \approx 1$ ) in the [100] construct as the Runt site is moved away from the promoter (Figure 5C). Thus, the direct repression model recapitulates repression by a single Runt molecule using the the same dissociation constant regardless of Runt binding site position, and displays the expected dependence of the Runt-RNAP interaction term on the distance between these two molecules.

Predicting repression by two-Runt binding sites requires both Runt-Runt and Runt-Runt-RNAP higher-order cooperativity

Could the parameters inferred in the preceding section be used to accurately predict repression in the presence of two Runt binding sites? An extra Runt binding site enables new protein-protein interactions between Runt molecules and RNAP (Figure 6A). First, we considered individual Runt-RNAP interaction terms, $ω_{r p 1}$ and $ω_{r p 2}$ , whose values were already inferred from the one-Runt binding site constructs as $ω_{r p_{[001]}}, ω_{r p_{[010]}}, a n d ω_{r p_{[100]}}$ (Figure 5D). Second, we considered protein-protein interactions (positive or negative) between two Runt molecules, $ω_{r r}$ . Third, following recent studies of Bicoid activation of the hunchback P2 minimal enhancer (Estrada et al., 2016a; Park et al., 2019), we also posited the existence of simultaneous Runt-Runt-RNAP higher-order cooperativity $ω_{r r p}$ . Given these different cooperativities, and as shown in detail in Figure 6—figure supplement 6B, the predicted rate of transcription is

\begin{array}{ll} R a t e = & R (p + b^{6} p ω_{b p} + r p (ω_{r p 1} + ω_{r p 2}) + r^{2} p ω_{r p 1} ω_{r p 2} ω_{r r} ω_{r r p} + b^{6} r p ω_{b p} (ω_{r p 1} + ω_{r p 2}) + \\ b^{6} r^{2} p ω_{b p} ω_{r p 1} ω_{r p 2} ω_{r r} ω_{r r p}) (1 + b^{6} (1 + 2 r + p ω_{b p}) + 2 r + p + r p (ω_{r p 1} + ω_{r p 2}) + r^{2} (ω_{r r} \\ {+ p ω_{r p 1} ω_{r p 2} ω_{r r} ω_{r r p}) + b^{6} r p ω_{b p} (ω_{r p 1} + ω_{r p 2}) + b^{6} r^{2} ω_{r r} + b^{6} r^{2} p ω_{b p} ω_{r p 1} ω_{r p 2} ω_{r r} ω_{r r p})}^{- 1} . \end{array}

Figure 6 with 7 supplements see all

Download asset Open asset

Prediction for the transcription initiation rate of *hunchback* P2 with two-Runt binding sites under different models of cooperativity.

(A) Direct repression model for *hunchback* P2 with two Runt binding sites featuring Runt-RNAP interaction terms given by $ω_{r p 1}$ and, $ω_{r p 2}$ Runt-Runt cooperativity captured by $ω_{r r}$ , and Runt-Runt-RNAP higher-order cooperativity accounted for by $ω_{r r p}$ . (B) Parameter-free model prediction for two Runt binding sites when the two Runt molecules bind the DNA and interact with RNAP independently of each other. (**C,D,E**) Best MCMC fits for the data for two-Runt binding site constructs for models with various combinations of cooperativity parameters. (C) Model incorporating Runt-Runt cooperativity. (D) Model incorporating Runt-Runt-RNAP higher-order cooperativity. (E) Model accounting for both Runt-Runt cooperativity and Runt-Runt-RNAP higher-order cooperativity. (F) Fixed or inferred parameters $ω_{r r}$ and $ω_{r r p}$ for all two-Runt binding site constructs. Note that $ω_{r r}$ is fixed to 1 for [011] and [101] constructs due to the fact that no Runt-Runt cooperativity is necessary to quantitatively describe the expression driven by these constructs; only the [110] construct is used to infer both $ω_{r r}$ and $ω_{r r p}$ . The horizontal line of $ω = 1$ denotes the case of no cooperativity. (G) Akaike Information Criterion (AIC) for all four scenarios of different free parameters shown throughout (**B–E**). (B-E, data points represent mean and standard error of the mean over the embryos. C-E, shaded error bars represent 95% confidence intervals for the best MCMC fits for the Runt WT datasets; F, data and error bars represent the mean and standard deviation of the posterior chain, while the standard deviation for the fixed $ω_{r r}$ is set to 0.).

Here, once again, we have color-coded parameters to be inferred in red to differentiate them from fixed parameters that were already inferred in previous sections. Despite the complexity of this equation, note that its only free parameters are the cooperativity parameters $ω_{r r}$ and $ω_{r r p}$ . As a result, we sought to determine whether the Runt-RNAP cooperativity terms, $ω_{r p 1}$ and $ω_{r p 2}$ , are sufficient to predict repression by two Runt molecules, or whether the cooperativities given by $ω_{r r}$ and $ω_{r r p}$ also need to be invoked.

Consider the simplest case where two Runt molecules bind and interact with RNAP independently from each other. Here, $ω_{r r} = 1$ , and $ω_{r r p} = 1$ . This model has no free parameters; all parameters have already been determined by the inferences performed on Runt null datasets and one-Runt binding site constructs (Figure 4 and Figure 5, respectively). While there was some agreement between the model and the data for the [101] construct (Figure 6B, center), significant deviations from the prediction occurred for the other two constructs. These deviations ranged from less repression than predicted for [011] (Figure 6B, left) to more repression than predicted for [110] (Figure 6B, right). Thus, this simple model of Runt independent repression is not supported by the experimental data, suggesting additional regulatory interactions between the Runt molecules and RNAP.

A first alternative to the independent repression model is the consideration of Runt-Runt cooperative interactions such as those that characterize many transcription factors (Park et al., 2019; Estrada et al., 2016b; He et al., 2010; Segal et al., 2008; Ptashne, 2004). However, adding a Runt-Runt cooperativity term, $ω_{r r}$ , was insufficient to account for the observed regulatory behavior (Figure 6C; Figure 6—figure supplement 4 more thoroughly analyzes this discrepancy). A second alternative consists in incorporating a Runt-Runt-RNAP higher-order cooperativity term, $ω_{r r p}$ . While the best MCMC fits revealed significant improvements in predictive power, important deviations still existed for the [110] construct (Figure 6D, right; Figure 6—figure supplement 5 more thoroughly analyzes the MCMC inference results).

Not surprisingly, given the agreement of the higher-order cooperativity model with the data for the [011] and [101] constructs (Figure 6D, left and center), this agreement persisted when both Runt-Runt cooperativity and Runt-Runt-RNAP higher-order cooperativity were considered (Figure 6E, left and center). However, including these two cooperativities also significantly improved the ability of the model at explaining the [110] experimental data (Figure 6E, right). Thus, while higher-order cooperativity is the main interaction necessary to quantitatively describe repression by two Runt repressors, pairwise cooperativity also needs to be invoked. This conclusion is supported by our MCMC sampling: posterior distributions for the Runt-Runt cooperativity term are not well constrained for the [011] or [101] constructs, whereas Runt-Runt-RNAP higher-order cooperativity is constrained very well across all constructs (Figure 6—figure supplement 6D; Figure 6—figure supplement 6 more thoroughly analyzes the MCMC inference results). As a result, accounting for both pairwise and higher-order cooperativity is necessary for the model to explain the observed rate of RNAP loading of all three constructs.

The higher-order cooperativity revealed by our analysis can lead to more or less repression than predicted by the independent repression model, motivating us to determine the magnitude of this cooperativity across constructs. To make this possible, we inferred the magnitude of the Runt-Runt cooperativity $ω_{r r}$ and the Runt-Runt-RNAP higher-order cooperativity $ω_{r r p}$ . As shown in Figure 6F, depending on the spatial arrangement of Runt binding sites, the Runt-Runt-RNAP higher-order cooperativity term $ω_{r r p}$ can be below or above 1. Note that, in doing these fits, we first set the Runt-Runt cooperativity, $ω_{r r}$ , values for [011] and [101] to 1 because, as we had demonstrated in Figure 6D, only the higher-order Runt-Runt-RNAP cooperativity was necessary. Thus, different placements of Runt molecules on the enhancer lead to distinct higher-order interactions with RNAP which, in turn, can result in less or more repression than predicted by a model where Runt molecules act independently of each other.

Repression by three-Runt binding sites also requires higher-order cooperativity

Building on our success in deploying thermodynamic models to explain repression by one- and two-Runt binding sites, we investigated repression by three-Runt binding sites. First, we accounted for pairwise interactions between Runt and RNAP, which were inferred from measurements of the one-Runt binding site constructs (Figure 1B), yielding $ω_{r p_{[001]}}, ω_{r p_{[010]}}$ , and $ω_{r p_{[100]}}$ from [001], [010], and [100]. Second, we considered pairwise protein-protein interactions between Runt molecules (Figure 1C), which were inferred from the two-Runt binding sites constructs through the parameters $ω_{r r_{[011]}}, ω_{r r_{[101]}}$ , and $ω_{r r_{[110]}}$ . Finally, we incorporated Runt-Runt-RNAP higher-order cooperativity acquired from the two-Runt binding sites constructs (Figure 1C) captured by $ω_{r r p_{[011]}}, ω_{r r p_{[101]}}$ , and $ω_{r r p_{[110]}}$ . we tested our model predictions using a similar scheme to that described in the previous section: we generated a parameter-free prediction for the initial rate of transcription by using the inferred parameters from the one- and two-Runt binding sites constructs, including the pairwise and higher-order interactions described above.

Figure 7A shows the resulting parameter-free prediction. As seen in the figure, our model could not qualitatively recapitulate the experimental data as it predicted too much repression. Such disagreement suggests that additional regulatory interactions are at play. Building on the need for higher-order cooperativity in the two-Runt binding site case, we propose the existence of higher-order cooperativities necessary to describe regulation by three Runt molecules—Runt-Runt-Runt higher-order cooperativity, $ω_{r r r}$ and Runt-Runt-Runt-RNAP higher-order cooperativity, $ω_{r r r p}$ (Figure 1D). The resulting expression for the predicted rate of transcription in the presence of all these sources of cooperativity is shown in Equation S10 in Section ‘Derivation of the general and simpler thermodynamic model for the hunchback P2 enhancer with one Runt binding site’. For simplicity, we assumed that the Runt-Runt-Runt cooperativity is one, and only determined the Runt-Runt-Runt-RNAP higher-order cooperativity. By including only a Runt-Runt-Runt-RNAP higher-order cooperativity parameter, our model recapitulated the experimental data (Figure 7B). Thus, our results further support the view in which the addition of Runt repressor binding motifs in an enhancer calls for the incorporation of cooperativities of increasingly higher-order.

Figure 7

Download asset Open asset

Prediction for *hunchback* P2 with three-Runt binding sites and multiple sources of cooperativity.

(A) Prediction using previously inferred Runt-RNAP, Runt-Runt, and Runt-Runt-RNAP cooperativity parameters. (B) Best MCMC fit obtained by incorporating an additional Runt-Runt-Runt-RNAP higher-order cooperativity parameter of $ω_{r r r p} = 857$ , corresponding to roughly $7 k_{B} T$ of free energy. (A,B, data points represent mean and standard error of the mean over>3 embryos; B, shaded regions represent 95% confidence intervals for the best MCMC fit.).

Discussion

One of the challenges in generating predictions to probe thermodynamic models is that, often, these models are contrasted against experimental data from endogenous regulatory regions (Segal et al., 2008; Sayal et al., 2016; Park et al., 2019). Here, the presence of multiple binding sites for several transcription factors—known and unknown (Vincent et al., 2016)—leads to models with a combinatorial explosion of free parameters. Like the proverbial elephant that can be fit with four parameters (Mayer et al., 2010), experiments with endogenous enhancers typically contain enough parameters to render it possible to explain away apparent disagreement between theory and experiment (Garcia et al., 2020).

To close this gap, synthetic minimal enhancers have emerged as an attractive alternative to endogenous enhancers (Fakhouri et al., 2010; Sayal et al., 2016; Park et al., 2019; Crocker et al., 2016). Here, the presence of only a handful of transcription factor binding sites and the ability to systematically control their placement and affinity dramatically reduce the number of free parameters in the model (Garcia et al., 2020). Inferences performed on these synthetic constructs could then inform model parameters that would make it possible to quantitatively predict transcriptional output of de novo enhancers (Sayal et al., 2016).

Building on these works, we sought to predict how the Runt repressor, which counteracts activation by Bicoid along the anterior-posterior axis of the early fly embryo (Chen et al., 2012), dictates output levels of transcription. To dissect repression, a strong and detectable level of expression in the absence of the repressor was needed, prompting us to choose a simple system of synthetic enhancers based on the strong hunchback P2 minimal enhancer (Garcia et al., 2013; Chen et al., 2012). This enhancer has been carefully studied in terms of its activator Bicoid and the pioneer-like transcription factor Zelda in the early embryo (Driever and Nüsslein-Volhard, 1988; Garcia et al., 2013; Park et al., 2019; Eck et al., 2020), making it easier to identify neutral sequences within the enhancer for introducing Runt binding sites (Chen et al., 2012). Further, when inserted into hunchback P2, Runt binding site number determines the level of transcription incrementally (Chen et al., 2012). Thus, hunchback P2 provided an ideal scaffold for quantitatively and systematically dissecting repression by Runt.

Previous studies using synthetic enhancers relied on measurements of input transcription factor patterns using fluorescence immunostaining, and of cytoplasmic mRNA patterns using fluorescence in situ hybridization (FISH) or single-molecule FISH. These fixed-tissue techniques have key differences from the live-imaging approach adopted here. First, given the dynamical nature of development, it is necessary to know when data were acquired. Doing so with high temporal resolution using FISH is challenging, although it can be accomplished to some degree by synchronizing embryo deposition before fixation (Park et al., 2019). Second, while most transcription factors directly dictate the rate of RNAP loading, and hence the rate of mRNA production (Spitz and Furlong, 2012; Garcia et al., 2013; Eck et al., 2020), typical FISH measurements report on the accumulated mRNA in the cytoplasm, which is a convolution of all processes of the transcription cycle—initiation, elongation, and termination (Liu et al., 2021; Alberts et al., 2015)—as well as mRNA nuclear export dynamics, diffusion, and degradation. These processes could be modulated in space and time, potentially confounding measurements. Here, we overcame these challenges by using the MS2 technique to precisely time our embryos and acquire the rate of transcription initiation. Of course, despite the ease of measuring the rate of transcription initiation using MS2, the accumulated mRNA is presumably a more relevant quantity for predicting downstream cellular decision making. Previous studies have shown that the MS2-MCP technique can also be used to quantify such patterns of accumulated mRNA, and that this quantification leads to results comparable to those obtained by smFISH (Garcia et al., 2013; Lammers et al., 2020). Following the same quantification method, we assessed the relationship between the initial rate of RNAP loading and the accumulated mRNA (Figure 3—figure supplement 5, Figure 3—figure supplement 6) by plotting them against each other. Reassuringly, as shown in Figure 3—figure supplement 8, our analysis revealed a strong correlation (with Pearson’s correlation coefficient of 0.90), supporting our claim that higher-order cooperativity is essential for explaining the action of multiple transcription factors during the development.

Interestingly, our initial dissection of constructs containing various combinations of Runt binding sites, but in the absence of Runt protein, revealed that unrepressed gene expression levels depend strongly on the number and placement of the binding sites within the enhancer (Figure 4A). These results challenge previous assumptions that unregulated gene expression levels stay unchanged as enhancer architecture is modulated (Sayal et al., 2016; Fakhouri et al., 2010; Barr et al., 2017), but they are in accordance with observations in bacterial systems (Garcia et al., 2012). As a result, our measurements call for accounting for unregulated levels in future quantitative dissections of eukaryotic enhancers, or to study relative magnitudes such as the fold-change in gene expression that has driven the dissection of bacterial transcriptional regulation (Phillips et al., 2019).

Using the thermodynamic model shown in Equation 3, we determined that the Bicoid-dependent parameters remain constant while RNAP-dependent parameters vary across these synthetic enhancer constructs. We speculate that the overall enhancer sequence, which changed as a result of the placement of different combinations of Runt binding sites within it, might affect the binding of the transcriptional machinery. Specifically, since the enhancer is proximal to the promoter, the transcriptional machinery might see slightly different DNA sequences in the vicinity of the promoter as suggested by published structures of the transcriptional machinery assembled on DNA (Louder et al., 2016).

Once we accounted for this difference in unrepressed gene expression levels, we determined that the repression profiles obtained for constructs bearing one-Runt binding site could be described by a simple thermodynamic model (Figure 2). Specifically, we showed that the same dissociation constant described Runt binding regardless of the position of its binding site along the enhancer (Figure 5A). Further, the Runt-RNAP interaction terms describing repressor action decreased as the binding site was placed farther from the promoter (Figure 5C), qualitatively consistent with a ‘direct repression’ model in which Runt needs to physically contact RNAP in order to realize its function (Jaynes and O’Farrell, 1991; Gray et al., 1994; Hewitt et al., 1999).

Although our model recapitulated repression by a one-Runt binding site, the inferred parameters were insufficient to quantitatively predict repression by two-Runt binding sites (Figure 6B). These results suggest that multiple repressors do not act independently of each other. Instead, new parameters describing both Runt-Runt cooperativity and Runt-Runt-RNAP higher-order cooperativity had to be incorporated into our models to quantitatively describe Runt action in these constructs (Figure 6—figure supplement 1C–E). An examination of the various cooperativity values inferred in the language of interaction energies (Table 1) revealed that these energies were of a magnitude comparable to protein-protein interaction energies previously measured in bacterial systems (Dodd et al., 2004; Bintu et al., 2005a; Amit et al., 2011). Interestingly, these interaction energies were both positive and negative, suggesting that both cooperativity or anti-cooperativity are at play depending on enhancer architecture (Amit et al., 2011). Additionally, the [101] construct showed a closer agreement with the parameter-free prediction, without invoking higher-order cooperativity, than the other two constructs ([110] or [011]). This further supports a picture where higher-order cooperativity is sensitive to the placement and orientation of transcription factor binding sites within regulatory regions.

Table 1

Interaction energies for the Runt-related cooperativity parameters from one-, two-, and three-Runt sites constructs.

Note that we used the Boltzmann relation of $ω = e x p (- E / (k_{B} T))$ , where the $E$ is the interaction energy, $k_{B}$ is the Boltzmann constant, and $T$ is the temperature.

Interaction energies for the Runt-related cooperativity parameters
model parameter	construct	interaction energy (K_BT)
Runt-RNAP interaction, $ω_{r p}$	[001]	2.34 ± 0.63
	[010]	1.36 ± 0.36
	[100]	0.18 ± 0.24
Runt-Runt interaction, $ω_{r r}$	[011]	0 (manually set)
	[110]	-0.95 ± 0.12
	[101]	0 (manually set)
Runt-Runt-RNAP interaction, $ω_{r r p}$	[011]	-2.09 ± 0.27
	[110]	4.15 ± 1.14
	[101]	1.12 ± 0.51
Runt-Runt-Runt-RNAP interaction, $ω_{r r r p}$	[111]	-2.12 ± 0.14

While we have long known about protein-protein cooperative interactions (Ackers et al., 1982), in the last few years it has become clear that higher-order cooperativity can also be at play in eukaryotic systems (Estrada et al., 2016a; Park et al., 2019; Biddle et al., 2020) as well as in bacteria (Dodd et al., 2004) and archaea (Peeters et al., 2013). The existence of this higher-order cooperativity suggests that, to predict gene expression from DNA sequence, it might be necessary to build an understanding of the many simultaneous interactions that precede transcriptional initiation. Our discovery of higher-order cooperativity in the action of multiple Runt molecules opens up new avenues to uncover the molecular nature of this phenomenon. For example, following an approach developed in Park et al., 2019, it could be possible to determine whether and how these cooperativity parameters are modulated upon perturbation of molecular players such as the Groucho or CtBP co-repressors, Big-brother, a co-factor facilitating the Runt binding to DNA, and components of the mediator complex (Park et al., 2019; Courey and Jia, 2001; Walrad et al., 2011). Indeed, Park et al., 2019 recently showed that co-activators and mediator units are involved in dictating the magnitude of similar higher-order cooperativity terms in activation by Bicoid. Thus, our thermodynamic models provide a lens through which to dissect the molecular underpinnings of Runt interactions with itself and with the transcriptional machinery.

Notably, the need to invoke cooperative interactions as more Runt binding sites are being added opposes our goal of predicting complex regulatory architectures from experiments with simpler architectures without the need to invoke new parameters. However, it will be interesting to determine whether more parameters need to be invoked as the number of Runt binding sites increases beyond three, or whether the parameters already inferred are sufficient to endow our models with parameter-free predictive power.

Importantly, while our model adopted a ‘direct repression’ view of the mechanism of Runt action, other mechanisms of repression such as ‘quenching’ could also describe the data. While all such models call for higher-order cooperativity to describe the data (Supplementary Section ‘Comparison of different modes of repression’), our data cannot differentiate among those models. Thus, we did not attempt to distinguish different molecular mechanisms of Runt transcriptional repression.

Finally, even though the work presented here has relied exclusively on thermodynamic models, it is important to note that a much more general approach based on kinetic models that are not in thermodynamic equilibrium could also be appropriate for describing our data. Indeed, an increasing body of work over the last few years has provided evidence for the necessity of invoking these more complex models in the context of transcriptional regulation in eukaryotes (Estrada et al., 2016a; Li et al., 2018; Park et al., 2019; Eck et al., 2020). In future work, it will be interesting to determine whether, when our data is viewed through the lens of these non-equilibrium models, invoking higher-order cooperativity is still necessary or whether, instead, simple pairwise protein-protein interactions suffice to reach an agreement between theory and experiment.

Overall, the work presented here establishes a framework for systematically and quantitatively studying repression in the early fly embryo. As showcased here, synthetic enhancers based on the hunchback P2 minimal enhancer constitute an ideal scaffold for the study of other repressors in early fly embryos. For example, we envision that this approach could be used to dissect repression by other transcription factors such as Capicua or Krüppel (Löhr et al., 2009; Sauer and Jäckle, 1991; Papagianni et al., 2018; Chen et al., 2012), and to probe observations of multiple repressors working together to oppose activation by Bicoid in establishing gene expression patterns along the anterior-posterior axis (Chen et al., 2012; Briscoe and Small, 2015). We anticipate that a similar approach could be used to dissect repression along the dorso-ventral axis of the embryo, by for example, adding repressor binding sites to well-established reporter constructs that are only regulated by the Dorsal activator (Jiang and Levine, 1993). Critically, we need to understand not only how one species of repressor works in concert with an activator, but also how multiple species of repressors work together as a system. The approach presented here provides a way forward for predictively understanding the complex gene regulatory network that shapes gene expression patterns in the early fly embryo.

Materials and methods

Generation of synthetic enhancers with MS2 reporter

Request a detailed protocol

The synthetic enhancer constructs used in this study are based off of Chen et al., 2012. In summary, the hunchback P2 enhancer was used as a scaffold to introduce Runt binding sites at different positions that are thought to be neutral (i.e. these Runt binding sites do not interfere with any other obvious binding sites for other transcription factors in the early Drosophila embryos as shown in Figure 4—figure supplement 1). For the three positions chosen to introduce Runt binding sites in Chen et al., 2012, the Gene Synthesis service from Genscript was used to generate synthetic enhancers with all possible configurations of zero-, one-, two-, and three-Runt binding sites in hunchback P2 as shown in Figure 1A. The enhancer sequences were placed into the original plasmid pIB backbone (Chen et al., 2012) using the Gene Fragment Synthesis service in Genscript, followed by the even-skipped promoter, and 24 repeats of the MS2v5 loop (Wu et al., 2015), the lacZ coding sequence, and the $α$ -Tubulin 3’UTR sequence (Chen et al., 2012) as shown in Table 2. These plasmids were injected into the 38F1 landing site using the RMCE method (Bateman et al., 2006) by BestGene Inc Flies were screened by selecting for white eye color and made homozygous. The orientation of the insertion was determined by genomic PCR to ensure a consistent orientation across all of our constructs. Specifically, we used two sets of primers that each amplified one of these two possible orientations: ‘Upward’, where the forward primer binds to a genomic location outside of 38F1 (TTCTAGTTCCAGTGAAATCCAAGCA) and the reverse primer binds to a location in our reporter transgene (ACGCCAGGGTTTTCCCAG), and ‘Downward’, where the forward primer remains the same as the ‘Upward’ set and the reverse primer binds to a location in our reporter transgene (CTCTGTTCTCGCTATTATTCCAACC) when the insertion is the opposite orientation to the ‘Upward’ orientation. As a result, only amplicons from either one of the orientations of insertion in the 38F1 landing site can be obtained. We chose the ‘Downward’ orientation for all our constructs.

Table 2

List of plasmids used to create the transgenic fly lines used in this study.

Plasmids
Name (hyperlinked to Benchling)	Function
pIB-hbP2-evePr-MS2v5-LacZ-Tub3UTR	[000]-MS2v5 reporter construct
pIB-hbP2+r 1 -far-evePr-MS2v5-LacZ-Tub3UTR	[100]-MS2v5 reporter construct
pIB-hbP2+r 1 -mid-evePr-MS2v5-LacZ-Tub3UTR	[010]-MS2v5 reporter construct
pIB-hbP2+r 1 -close-evePr-MS2v5-LacZ-Tub3UTR	[001]-MS2v5 reporter construct
pIB-hbP2+r 2-2 +3 -evePr-MS2v5-LacZ-Tub3UTR	[011]-MS2v5 reporter construct
pIB-hbP2+r 2-1 +3 -evePr-MS2v5-LacZ-Tub3UTR	[101]-MS2v5 reporter construct
pIB-hbP2+r 2-1 +2 -evePr-MS2v5-LacZ-Tub3UTR	[110]-MS2v5 reporter construct
pIB-hbP2+r 3 -evePr-MS2v5-LacZ-Tub3UTR	[111]-MS2v5 reporter construct
pHD-scarless-LlamaTag-Runt	Donor plasmid for LlamaTag-Runt CRISPR knock-in fusion for the N-terminal
pU6:3-gRNA(Runt-N-2)	gRNA plasmid for LlamaTag-Runt CRISPR knock-in fusion for the N-terminal
pCasper-vasa-eGFP	vasa maternal driver for ubiquitous eGFP expression in the early embryo

CRISPR-Cas9 knock-in of the green LlamaTag in the endogenous runt locus

Request a detailed protocol

We used CRISPR-Cas9 mediated Homology Directed Repair (HDR) to insert the LlamaTag against eGFP into the N-terminal of the runt endogenous locus (Bothma et al., 2018; Gratz et al., 2015). The donor plasmid was constructed by stitching individual fragments—PCR amplified left/right homology arms from the endogenous runt locus roughly 1 kb in length each, LlamaTag, and pHD-scarless vector—using Gibson assembly (Gratz et al., 2015). The PAM sites in the donor plasmid were mutated such that the Cas9 only cleaved the endogenous locus, not the donor plasmid, without changing the amino acid sequence of the Runt protein. The final donor plasmid contained the 3xP3-dsRed marker such that dsRed is expressed in the fly eye and ocelli for screening. Positive transformant flies were screened using a fluorescence dissection scope and set up for single fly crosses to establish individual lines that were then verified with PCR amplification and Sanger sequencing (UC Berkeley Sequencing Facility). Importantly, this llamaTag-runt allele rescues development to adulthood as a homozygous. Thus we concluded that the LlamaTag-Runt allele can be used to monitor the behavior of endogenous Runt protein.

Fly strains

Request a detailed protocol

Transcription from the synthetic enhancer reporter constructs was measured by using embryos from crossing yw;his2av-mRFP1;MCP-eGFP(2) females and yw;synthetic enhancer-MS2v5-lacZ;+ males as described in Garcia et al., 2013; Eck et al., 2020; Lammers et al., 2020.

eGFP-Bicoid measurements were performed using the fly line from Gregor et al., 2007. The LlamaTag-Runt measurements were done using the fly line LlamaTag-Runt; +; vasa-eGFP, His2Av-iRFP illustrated in Table 3. Briefly, eGFP was supplied by a vasa maternal driver. Females carrying both the LlamaTag-Runt and the vasa-driven eGFP were crossed with males carrying the LlamaTag-Runt, the progeny from this cross were imaged and then recovered to determine the embryo’s sex using PCR. PCR was run with three sets of primers: Y chr1 (Forward: CGATCCAGCCCAATCTCTCATATCACTA, Reverse: ATCGTCGGTAATGTGTCCTCCGTAATTT), Y chr2 (Forward: AACGTAACCTAGTCGGATTGCAAATGGT, Reverse: GAGGCGTACAATTTCCTTTCTCATGTCA), and Auto1 (Forward: GATTCGATGCACACTCACATTCTTCTCC, Reverse: GCTCAGCGCGAAACTAACATGAAAAACT). Two of primers in the set (Y chr1 and Y chr2) bind to the Y chromosome while the other one (Auto1) binds to one of the autosomes and constitutes a positive control (Lott et al., 2011). The Histone-iRFP fly line was from Pan et al., 2022, and was used for nuclei segmentation to extract nuclear flourescence from the eGFP channel.

Table 3

List of fly lines used in this study and their experimental usage.

Fly lines
Genotype	Use
LlamaTag-Runt; +; vasa-eGFP, His2Av-iRFP	Visualize LlamaTagged Runt protein and label nuclei
LlamaTag-Runt; +; MCP-eGFP(4F), His2Av-iRFP	Visualize LlamaTagged Runt protein, nascent transcripts and label nuclei
run3/FM6; +; +	Visualize LlamaTagged Runt protein, nascent transcripts and label nuclei
yw; His2Av-mRFP; MCP-eGFP	Females to label nascent RNA and nuclei
yw; [000]-MS2v5; +	Males carrying the MS2 reporter transgene
yw; [100]-MS2v5; +	Males carrying the MS2 reporter transgene
yw; [010]-MS2v5; +	Males carrying the MS2 reporter transgene
yw; [001]-MS2v5; +	Males carrying the MS2 reporter transgene
yw; [011]-MS2v5; +	Males carrying the MS2 reporter transgene
yw; [101]-MS2v5; +	Males carrying the MS2 reporter transgene
yw; [110]-MS2v5; +	Males carrying the MS2 reporter transgene
yw; [111]-MS2v5; +	Males carrying the MS2 reporter transgene

To generate the embryos that are zygotic null for the runt allele, we used a fly cross scheme consisting of two crosses. In the first generation, we crossed LlamaTag-Runt;+;+ males with run3/FM6;+;MCP-eGFP(4 F),his2av-mRFP1 females. run3 is the null allele for runt, missing around 5 kb including the coding sequence of the runt locus (Gergen and Butler, 1988; Chen et al., 2012). The MCP-eGFP(4 F) transgene expresses approximately twice the amount of MCP protein than the MCP-eGFP(2) (Garcia et al., 2013; Eck et al., 2020) and thus results in similar levels of MCP to those of MCP-eGFP(2) in the trans-heterozygotes. The female progeny from this cross, LlamaTag-Runt/run3;+;MCP-eGFP(4 F),his2av-mRFP1/+ was then crossed with males whose genotype was LlamaTag-Runt/Y;synthetic enhancer-MS2v5-lacZ;+ to produce the embryos that we used for live imaging. The resulting embryos carried maternally supplied MCP-eGFP and His-RFP for visualization of nascent transcripts and nuclei. The X chromosome contained a LlamaTag-Runt allele or run3 null allele. We could differentiate between these two genotypes because, when the embryo had the Runt allele, a stripe pattern would appear in late nc14. We imaged all embryos until late nc14 to make sure that we were capturing the nulls.

Sample preparation and data collection

Request a detailed protocol

Sample preparation was done following the protocols described in Garcia et al., 2013. Briefly, embryos were collected, dechorionated with bleach for 1–2 min, and then mounted between a semipermeable membrane (Lumox film, Starstedt, Germany) and a coverslip while embedded in Halocarbon 27 oil (Sigma-Aldrich). Live imaging was performed using a Leica SP8 scanning confocal microscope, a White Light Laser and HyD dectectors (Leica Microsystems, Biberach, Germany). Imaging settings for the MS2 experiments with the presence of MCP-eGFP and Histone-RFP were the same as in Eck et al., 2020 except that we used a 1024x245 pixel format to image a wider field of view along the anterior-posterior axis. The settings for the eGFP-Bicoid measurements were the same as described in Eck et al., 2020.

The settings for the eGFP:LlamaTag-Runt measurements were similar to that of eGFP-Bicoid except for the following. To increase our imaging throughput, we utilized the ‘Mark and Position’ functionality in the LASX software (Leica SP8) to image 5–6 embryos simultaneously. To account for the decreased time resolution, we lowered the z-stack size from 10 μm to 2.5 μm, keeping the 0.5 μm z-step. By doing this, we could maintain 1-min frame rate for each imaged embryo. Additionally, these flies expressed Histone-iRFP, instead of Histone-RFP as in Eck et al., 2020, so that we used a 670 nm laser at 40 μW (measured at a 10x objective) for excitation of the histone channel, and the HyD detector was set to a 680 nm-800 nm spectral window (Figure 3—figure supplement 7).

Image analysis

Request a detailed protocol

Images were analyzed using custom-written software (MATLAB, mRNA Dynamics Github repository; Garcia Lab @ UC Berkeley, 2022) following the protocol in Garcia et al., 2013 and Eck et al., 2020. Briefly, this procedure involved segmentation and tracking of nuclei and transcription spots. First, segmentation and tracking of individual nuclei were done using the histone channel as a nuclear mask. Second, segmentation of each transcription spot was done based on its fluorescence intensity and existence over multiple z-stacks. The intensity of each MCP-GFP transcriptional spot was calculated by integrating pixel intensity values in a small window around the spot and subtracting the background fluorescence measured outside of the active transcriptional locus. When there was no detectable transcriptional activity, we assigned NaN values for the intensity. The tracking of transcriptional spots was done by using the nuclear tracking and proximity of transcriptional spots between consecutive time points. The nuclear protein fluorescence intensities from the eGFP-Bicoid and LlamaTag-Runt fly lines, which we use as a proxy for the protein nuclear concentration, were calculated as follows. Using the nuclear mask generated from the histone channel, we performed the same nuclear segmentation and tracking as described above for the MS2 spots. Then, for every z-section, we extracted the integrated fluorescence over a $2 μ m$ diameter circle on the xy-plane centered on each nucleus. For each nucleus, the recorded fluorescence corresponded to the z-position where the fluorescence was maximal. This resulted in an average nuclear concentration as a function of time for each single nucleus. These concentrations from individual nuclei were then averaged over a narrow spatial window (2.5% of the embryo length) to generate the spatially averaged protein concentration reported in the main text. For the eGFP:LlamaTag-Runt datasets, we had to subtract the background eGFP fluorescence due to the presence of an unbound eGFP population (Bothma et al., 2018). We used the same protocol described in Bothma et al., 2018 and in the Supplementary Section ‘Quantifying the nuclear concentration of LlamaTag-Runt’ to extract this background.

Bayesian inference procedure: Markov Chain Monte Carlo sampling

Request a detailed protocol

Parameter inference was done using the Markov Chain Monte Carlo (MCMC) method. We used a well-established package MCMCstat that uses an adaptive MCMC algorithm (Haario et al., 2006; Haario et al., 2001). A detailed description on how we performed the MCMC parameter inference, for example setting the priors and bounds for parameters, can be found in Supplementary Section ‘Markov Chain Monte Carlo inference protocol’.

Appendix 1

S1 Derivation of the general thermodynamic model for the hunchback P2 enhancer

In this section, we rederive the thermodynamic model presented in the main text, now without the assumption of strong Bicoid-Bicoid cooperativity. The equilibrium thermodynamic modeling framework that we used in this paper is described in more detail in Bintu et al., 2005a; Bintu et al., 2005b.

We start by modeling the case of hunchback P2 without any Runt binding sites, which is believed to have at least six Bicoid binding sites (Park et al., 2019; Driever et al., 1989). As shown by the states and weights presented in Figure 2—figure supplement 1A, in our thermodynamic model, we assume that the six Bicoid binding sites have the same dissociation constant given by $K_{b}$ , and we posit that RNAP-promoter binding is governed by a dissociation constant given by $K_{p}$ . We also assume pairwise cooperativity between Bicoid molecules given by $ω_{b}$ , and cooperativity between each Bicoid molecule and RNAP given by $ω_{b p}$ . For simplicity, we will use the dimensionless parameters $b = [B i c o i d] / K_{b}$ and $p = [R N A P] / K_{p}$ , where $[B i c o i d]$ , and $[R N A P]$ are the concentrations of Bicoid and RNAP, respectively, and $K_{b}$ and $K_{p}$ are their corresponding dissociation constants.

We factor the total partition function into two categories: $Z_{b}$ corresponding to states that only have Bicoid bound, and $Z_{b p}$ describing states with both Bicoid and RNAP bound. Then, we calculate each component separately. The sum of microstates for $Z_{b}$ is

Z_{b} = 1 + 6 b + 15 b^{2} ω_{b} + \dots + b^{6} ω_{b}^{5} = 1 + \sum_{i = 1}^{6} (\binom{6}{i}) b^{i} ω_{b}^{i - 1} .

Using the binomial theorem, we can simplify Equation S1 leading to

Z_{b} = 1 + \sum_{i = 1}^{6} (\binom{6}{i}) b^{i} ω_{b}^{i - 1} = 1 + \frac{1}{ω_{b}} [(1 + b ω_{b})^{6} - 1] .

Using the same logic, we obtain $Z_{b p}$ such that

Z_{b p} = (p + p \sum_{i = 1}^{6} (\binom{6}{i}) b^{i} ω_{b}^{i - 1} ω_{b p}^{i}) = p + \frac{p}{ω_{b}} [(1 + b ω_{b} ω_{b p})^{6} - 1] .

Using these two partition functions, we then calculate the probability of the promoter being bound by RNAP, $p_{b o u n d}$ as

P_{b o u n d} = \frac{Z_{b p}}{Z_{b} + Z_{b p}} = \frac{p + \frac{p}{ω_{b}} [{(1 + b ω_{b} ω_{b p})}^{6} - 1]}{1 + \frac{1}{ω_{b}} [{(1 + b ω_{b})}^{6} - 1] + p + \frac{p}{ω_{b}} [{(1 + b ω_{b} ω_{b p})}^{6} - 1]} .

Following recent work (Gregor et al., 2007; Park et al., 2019), we now assume that the Bicoid-Bicoid pairwise cooperativity is very strong ( $ω_{b} ≫ 1$ ). We can then simplify Equation S4 to obtain

P_{b o u n d} = \frac{p + p b^{6} ω_{b}^{5} ω_{b p}^{6}}{1 + p + b^{6} ω_{b}^{5} + p b^{6} ω_{b}^{5} ω_{b p}^{6}} .

If we now define a new binding constant for Bicoid, $K_{b}^{'} = K_{b} * {(\frac{1}{ω_{b}})}^{\frac{5}{6}}$ , such that $b^{'} = b ω_{b}^{\frac{5}{6}}$ , and a new cooperativity term between Bicoid and RNAP given by $ω_{b p}^{'} = ω_{b p}^{6}$ , we can then rewrite Equation S5 as

P_{b o u n d} = \frac{p + b^{', 6} p ω_{b p}^{'}}{1 + p + b^{', 6} + b^{', 6} p ω_{b p}^{'}},

which is the expression we use throughout the main text. Thus, strong pairwise cooperativity between Bicoid molecules leads to a functional form where only the state with all Bicoid molecules bound remain (six in this case). This strong cooperativity can explain the sharp step-like expression pattern along the embryo’s anterior-posterior axis of the hunchback gene (Figure 3J; Gregor et al., 2007, Park et al., 2019, Driever and Nüsslein-Volhard, 1988; Driever and Nüsslein-Volhard, 1989).

S2 Derivation of the general and simpler thermodynamic model for the hunchback P2 enhancer with one Runt binding site

Having derived the equation for the strong cooperative binding of Bicoid to the wild-type hunchback P2 enhancer, we will now extend that model to the case of hunchback P2 with one Runt binding site. The corresponding states and weights of our full model are shown in Figure 2—figure supplement 2A.

Using a similar logic for calculating the partition functions as described in the previous section, we can compute the probability of the promoter being bound by RNAP as

\begin{array}{ll} p_{b o u n d} = \frac{\overset{Bicoid and RNAP}{\overset{⏞}{(p + p \sum_{i = 1}^{6} (\binom{6}{i}) b^{i} ω_{b}^{i - 1} ω_{b p}^{i})}} + \overset{Bicoid, Runt, and RNAP}{\overset{⏞}{(r p ω_{r p} + r p ω_{r p} \sum_{i = 1}^{6} (\binom{6}{i}) b^{i} ω_{b}^{i - 1} ω_{b p}^{i})}}}{\underset{Bicoid only}{\underset{⏟}{(1 + \sum_{i = 1}^{6} (\binom{6}{i}) b^{i} ω_{b}^{i - 1})}} + \underset{Bicoid and RNAP}{\underset{⏟}{(p + p \sum_{i = 1}^{6} (\binom{6}{i}) b^{i} ω_{b}^{i - 1} ω_{b p}^{i})}} + \underset{Bicoid and Runt}{\underset{⏟}{(r + r \sum_{i = 1}^{6} (\binom{6}{i}) b^{i} ω_{b}^{i - 1})}} + \underset{Bicoid, Runt, and RNAP}{\underset{⏟}{(r p ω_{r p} + r p ω_{r p} \sum_{i = 1}^{6} (\binom{6}{i}) b^{i} ω_{b}^{i - 1} ω_{b p}^{i})}}}, \end{array}

where, in addition to the parameters defined in the above section for the wild-type hunchback P2 case in the absence of Runt, we have added two parameters: the dissociation constant for Runt given by $K_{r}$ , and a Runt-RNAP interaction term (an anti-cooperativity), $ω_{r p}$ . Using the binomial theorem as in Equation S2, we can simplify Equation S7 to obtain

p_{b o u n d} = \frac{p + \frac{p}{ω_{b}} [{(1 + b ω_{b} ω b p)}^{6} - 1] + r p ω_{r p} + \frac{r p ω_{r p}}{ω_{b}} [{(1 + b ω_{b} ω_{b p})}^{6} - 1]}{1 + \frac{1}{ω_{b}} [{(1 + b ω_{b})}^{6} - 1] + p + \frac{p}{ω_{b}} [{(1 + b ω_{b} ω b p)}^{6} - 1] + r + \frac{r}{ω_{b}} [{(1 + b ω_{b})}^{6} - 1] + r p ω_{r p} + \frac{r p ω_{r p}}{ω_{b}} [{(1 + b ω_{b} ω_{b p})}^{6} - 1]} .

We now again assume that Bicoid-Bicoid cooperativity is very strong such that $ω_{b} ≫ 1$ . Then, we can combine Equation S8 with Equation 1 to obtain

R a t e = R p_{b o u n d} = R \frac{p + b^{', 6} p ω_{b p} + r p ω_{r p} + b^{', 6} r p ω_{b p}^{'} ω_{r p}}{1 + b^{', 6} + r + b^{', 6} r + p + b^{', 6} p ω_{b p}^{'} + r p ω_{r p} + b^{', 6} r p ω_{b p}^{'} ω_{r p}},

where the new parameters, $b^{'}$ and $ω_{b p}^{'}$ are defined in the same way as in Equation S6. The effective states and weights remaining after taking this limit are shown in Figure 2—figure supplement 2B. Similarly, we can derive expressions for $p_{b o u n d}$ in the presence of two and three Runt binding sites, and in the strong Bicoid-Bicoid cooperativity limit in order to obtain the predictions used throughout this text. We show this expression for two Runt binding sites in Equation 6. Further, for the case of repression by three Runt binding sites, the rate of transcription is given by

\begin{array}{ll} R a t e = & R (p + b^{6} p ω_{b p} + r p (ω_{r p 1} + ω_{r p 2} + ω_{r p 3}) + b^{6} r p ω_{b p} (ω_{r p 1} + ω_{r p 2} + ω_{r p 3}) + \\ r^{2} p (ω_{r p 1} ω_{r p 2} ω_{r r 1} ω_{r r p 1} + ω_{r p 2} ω_{r p 3} ω_{r r 2} ω_{r r p 2} + ω_{r p 3} ω_{r p 1} ω_{r r 3} ω_{r r p 3}) + \\ r^{3} p ω_{r p 1} ω_{r p 2} ω_{r p 3} ω_{r r 1} ω_{r r 2} ω_{r r 3} ω_{r r p 1} ω_{r r p 2} ω_{r r p 3} ω_{r r r} ω_{r r r p} + \\ b^{6} r^{2} p ω_{b p} (ω_{r p 1} ω_{r p 2} ω_{r r 1} ω_{r r p 1} + ω_{r p 2} ω_{r p 3} ω_{r r 2} ω_{r r p 2} + ω_{r p 3} ω_{r p 1} ω_{r r 3} ω_{r r p 3}) + \\ b^{6} r^{3} p ω_{b p} ω_{r p 1} ω_{r p 2} ω_{r p 3} ω_{r r 1} ω_{r r 2} ω_{r r 3} ω_{r r p 1} ω_{r r p 2} ω_{r r p 3} ω_{r r r} ω_{r r r p}) \\ (1 + b^{6} (1 + 3 r + p ω_{b p}) + 3 r + p + r p (ω_{r p 1} + ω_{r p 2} + ω_{r p 3}) + \\ r^{2} p (ω_{r p 1} ω_{r p 2} ω_{r r 1} ω_{r r p 1} + ω_{r p 2} ω_{r p 3} ω_{r r 2} ω_{r r p 2} ω_{r p 3} ω_{r p 1} ω_{r r 3} ω_{r r p 3}) + b^{6} r^{2} (ω_{r r 1} + ω_{r r 2} + ω_{r r 3}) + \\ b^{6} r^{3} ω_{r r 1} ω_{r r 2} ω_{r r 3} ω_{r r r} + r^{2} (ω_{r r 1} + ω_{r r 2} + ω_{r r 3}) + r^{3} ω_{r r 1} ω_{r r 2} ω_{r r 3} + \\ r^{3} p ω_{r p 1} ω_{r p 2} ω_{r p 3} ω_{r r 1} ω_{r r 2} ω_{r r 3} ω_{r r p 1} ω_{r r p 2} ω_{r r p 3} ω_{r r r} ω_{r r r p} + \\ b^{6} r p ω_{b p} (ω_{r p 1} + ω_{r p 2} + ω_{r p 3}) + b^{6} r^{2} p ω_{b p} (ω_{r p 1} ω_{r p 2} ω_{r r 1} ω_{r r p 1} + ω_{r p 2} ω_{r p 3} ω_{r r 2} ω_{r r p 2} + ω_{r p 3} ω_{r p 1} ω_{r r 3} ω_{r r p 3}) + \\ {b^{6} r^{3} p ω_{b p} ω_{r p 1} ω_{r p 2} ω_{r p 3} ω_{r r 1} ω_{r r 2} ω_{r r 3} ω_{r r p 1} ω_{r r p 2} ω_{r r p 3} ω_{r r r} ω_{r r r p})}^{- 1}, \end{array}

where the parameters are defined as in Figure 1 and Section ‘Repression by three-Runt binding sites also requires higher-order cooperativity’.

S3 Comparing using static versus dynamic transcription factor concentrations as model inputs

In this section, we tested whether using static, time-averaged transcription factor concentration profiles yielded comparable theoretical predictions than when instead acknowledging the fact that input transcription factor concentration changes over time. Briefly, we compared the predicted rate of transcription calculated in two ways: (1) time-averaging the instantaneous rate from the dynamic transcription factor concentration profiles over a specified time window (from 5 to 10 minutes from the 13th anaphase) and (2) using static input transcription factors already time-averaged over the same time window.

As a concrete example, we focused on the hunchback P2 enhancer with one Runt binding site. We calculated the predicted rate of transcription using the thermodynamic model given by Equation 2. First, we performed this calculation using the dynamic concentration profiles of Bicoid and Runt shown in Figure 3B and D, respectively. Briefly, the terns $b$ and $r$ in Equation 2 now become functions of time such that

R a t e (t) = R \frac{p + b^{6} (t) p ω_{b p} + r (t) p ω_{r p} + b^{6} (t) r (t) p ω_{b p} ω_{r p}}{1 + b^{6} (t) + r (t) + b^{6} (t) r (t) + p + b^{6} (t) p ω_{b p} + r (t) p ω_{r p} + b^{6} (t) r (t) p ω_{b p} ω_{r p}},

where $b (t) = [B i c o i d] (t) / K_{b}$ and $r (t) = [R u n t] (t) / K_{r}$ . We choose a set of reasonable values for the model parameters to illustrate the calculation of $R a t e (t)$ at 30% of the embryo length. The resulting dynamic rate of transcription profile is shown in Figure 3—figure supplement 1A (blue curve). We then use this profile to calculate the time-averaged rate of transcription over the time window of 5–10 minutes from the 13th anaphase, resulting in the green area shown in Figure 3—figure supplement 1A.

The predicted average rate of RNAP loading given dynamic input transcription factors can be compared to the predicted rate of RNAP loading given the average input concentrations that we used throughout the main text (Figure 3E). Specifically, we plug the static concentration profiles of Bicoid and Runt shown in Figure 3E into Equation 2 to obtain the red area shown in Figure 3—figure supplement 1A. As shown in the figure, the predicted rate of transcription obtained by these two analysis methodologies are equivalent within error.

Finally, we performed this comparison between different approaches to calcualate the rate of transcription as a function of position along the embryo (from 20% to 70% of the embryo length). As shown in Figure 3—figure supplement 1B, the resulting spatial profiles are comparable within error. Thus, we have shown that our approach of using time-averaged, static transcription factor concentrations as inputs to our model yield quantitatively equivalent result as accounting for the dynamic concentration profiles of these transcription factors.

S4 Quantitative interpretation of MS2 signals

The MS2 signal reports on three features of transcriptional dynamics: (1) the initial RNAP loading rate, (2) the duration of transcription, and (3) the fraction of loci that engage in transcription at any time point in the nuclear cycle. In this section, we will explain in further detail how we extract these features from the MS2 signal over nuclear cycle 14.

S4.1 Extracting the initial RNAP loading rate

The initial rate of RNAP loading corresponds the average transcription rate observed after transcriptional onset and until the MS2 signal reaches its peak value during nuclear cycle 14. In order to measure this rate, we followed the protocol described in Garcia et al., 2013. Briefly, as shown in Figure 3—figure supplement 2A, we fitted a line to the MS2 time trace (averaged over nuclei within a spatial window of 2.5% of the embryo length) within the time window of 5–10 minutes after the 13th anaphase. The slope of this line reported on the initial rate of RNAP loading (Figure 3G). The spatial profiles of this initial rate of RNAP loading across all our synthetic enhancer constructs and genotypes are shown in Figure 3—figure supplement 2B.

S4.2 Extracting the duration of transcription

In the main text, we focused on the theoretical prediction of the initial rate of transcription. However, the length of the time window over which transcription occurs (Lammers et al., 2020) is another regulatory knob that, in principle, Runt could modulate to dictate gene expression patterns. We sought to determine the duration of time over which transcription occurs to assess whether Runt affects not only the initial rate of transcription, but also the time window over which transcription could initiate. To quantify the effective duration of transcription initiation, we resorted to the analysis methodology developed in Garcia et al., 2013. Briefly, we parametrized the MS2 signal decay regime—after transcription reaches its peak and becomes slower than the unloading rate (Garcia et al., 2013)—as an exponential decay (Figure 3—figure supplement 3A). Thus, we can describe the MS2 spot fluorescence trace in the decay regime as

F l u o (t) = F l u o_{m a x} e^{- (t - T_{p e a k}) / τ},

where $T_{p e a k}$ represents the time point where the MS2 spot fluorescence reaches its peaks, $F l u o_{m a x}$ is the maximum level of fluorescence, and $τ$ is the decay time.

Given the sometimes noisy MS2 traces (data not shown), we fitted an exponential curve to the more robust integral of the MS2 spot fluorescence over time from $T_{p e a k}$ to the end of nuclear cycle 14 as shown in Figure 3—figure supplement 3B. This quantity is proportional to the amount of mRNA produced between the integration bounds (Garcia et al., 2013). The resulting accumulated mRNA time trace is then fitted to the integrated form of Equation S12, which is given by

m R N A (t) = m R N A_{m a x} (1 - e^{(t - T_{p e a k}) / τ}),

where $m R N A_{m a x}$ is the accumulated mRNA at the end of nuclear cycle 14.

The resulting profiles of the duration of transcription along the embryo for our all synthetic enhancer constructs are illustrated in Figure 3—figure supplement 3C in the presence and absence of Runt protein. As shown in the figure, this duration time is not significantly modulated by Runt repressor.

S4.3 Calculation of the fraction of competent nuclei

Another quantity that could be modulated by Runt repressor is the fraction of loci that ever engage in transcription during a given nuclear cycle, which we termed as the “fraction of competent loci”. As demonstrated by Garcia et al., 2013, Dufourt et al., 2018, Lammers et al., 2020 and Eck et al., 2020, this fraction of transcriptionally competent loci is modulated along the anterior-posterior axis, presumably due to the action of transcription factor gradients.

To show a concrete example of how this quantity is calculated, we take data from one construct ([000]) showing the MS2 spot fluorescence time traces from individual loci of transcription as shown in Figure 3—figure supplement 4A. Here, columns represent time points during nuclear cycle 14, and rows represent individual transcriptional loci. As shown in the figure, roughly 80% of the loci, labeled as “competent loci”, show active transcription during nuclear cycle 14. However, the remaining 20% of the loci never engage in transcription, which we termed as “incompetent loci”. Because these two populations exhibit wildly different behaviors, we define the fraction of competent loci as

fraction of competent loci = \frac{number of competent loci}{number of total loci} .

Thus, in this example in Figure 3—figure supplement 4A, the fraction of competent loci is approximately 0.8.

Figure 3—figure supplement 4B shows the measured fraction of active loci for all synthetic enhancer constructs in the presence and absence of Runt repressor. As seen in the figure, although this quantity can be modulated by the presence of Runt repressor, this is not always the case (e.g., [011] and [111]). Moreover, we could not find a trend for how the fraction of competent loci is modulated by different combinations of Runt binding sites. For example, the [100] construct alone did show a change in the fraction of active loci in the presence of Runt, whereas the [010] construct did not. When these two binding sites were combined as the [110], there was no significant modulation of the fraction of competent loci when adding Runt repressor. In another example, the [001] construct showed a mild modulation of the fraction of competent loci. However, when this Runt binding site was combined with the [010], which did not show any modulation, the [011] construct showed a much bigger modulation of the fraction of competent loci than the [001]. Thus, the [010] Runt binding site could drive more or less modulation of the fraction of competent loci when combined with different Runt binding sites in a context-dependent manner. As a result of our failure to uncover an apparent trend in terms of which regulatory architectures lead to a stronger modulation of the fraction of active loci, we did not attempt to theoretically explain the regulation of this fraction of active loci in this study.

S5 Design of synthetic enhancer constructs based on the hunchback P2 enhancer

The Runt binding sites were introduced into the hunchback P2 minimal enhancers at the positions determined by Chen et al., 2012. To make this possible, the authors chose positions containing presumed neutral DNA sequences, meaning that these DNA locations did not contain obvious motifs for Bicoid or Zelda, the major input transcription factors that regulate this enhancer. Then, these DNA sequences were mutated to turn them into Runt binding sites.

To ensure that this process did not perturb the binding sites for Bicoid and Zelda we resorted to the Advanced PATSER entry form (Hertz et al., 1990; Hertz and Stormo, 1999) which identifies the location of transcription factor binding sites from a sequence of DNA based on position weight matrices. We used position weight matrices for Bicoid and Zelda from Park et al., 2019. PATSER was run with the settings described in Eck et al., 2020 for both the hunchback P2 enhancer and the hunchback P2 enchancer with three Runt binding sites (from Chen et al., 2012) for Bicoid and Zelda, respectively. The result of this analysis for these two constructs is shown for each transcription factor in Figure 4—figure supplement 1A. Here, we took a PATSER score cutoff—for considering a given sequence to be a binding site—of 3 as in Eck et al., 2020. We observed that the recognized binding motifs for both Bicoid and Zelda were identical between the two constructs, meaning that we did not add additional Bicoid or Zelda binding sites by introducing the Runt motifs. The resulting synthetic enhancer with three Runt binding sites with mapped binding sites for Bicoid, Zelda (Figure 4—figure supplement 1A), and Runt (Chen et al., 2012) is shown in Figure 4—figure supplement 1B as a reference. The position of the Runt binding sites are noted from their distance from the promoter (which is marked as 0).

S6 Markov Chain Monte Carlo inference protocol

Markov Chain Monte Carlo (MCMC) sampling is a widely used technique for robust parameter estimation using Bayesian statistics (Geyer and Thompson, 1992; Sivia and Skilling, 2006). We used the MATLAB package MCMCstat, an adaptive MCMC technique, which we could directly implement downstream of our data analysis pipeline (Haario et al., 2006; Haario et al., 2001). Detailed instructions on how to implement the MCMCstat package can be found in https://mjlaine.github.io/mcmcstat/, (copy archived at swh:1:rev:ca0cbd288f03c7a29050b3d6698a96b45ccfa4b2; Laine, 2021).

MCMC allows for an estimation of the set of parameter values of a model that best explain the experimental data along with their associated errors. In this work, we used MCMC to infer the set of best fit values of the parameters in our thermodynamic models given the observed profile of the rate of transcription initiation along the anterior-posterior axis of the embryo.

MCMC calculates a Bayesian posterior probability distribution of each free parameter given the data by stochastically sampling different parameter values. For a given set of observations $D$ and a model with parameters $θ$ , the posterior probability distribution of a particular set of values is given by Bayes’ theorem

\underset{posterior}{\underset{⏟}{p (θ | D)}} \propto \underset{likelihood}{\underset{⏟}{p (D | θ)}} \underset{prior}{\underset{⏟}{p (θ)}} .

The prior function represents the a priori assumption about the probability distribution of parameter values $θ$ . Here, we assumed a uniform prior distribution for all parameters to reflect our ignorance about the model parameters within the following intervals:

$K_{b}$ : $[0, 100]$ AU
$ω_{b p}$ : $[0, 200]$
$p$ : $[0, 1]$
$R$ : $[0, 400]$ AU/min
$K_{r}$ : $[0, 100]$ AU
$ω_{r p}$ : $[0, 1.2]$
$ω_{r r}$ : $[0, 100]$
$ω_{r r p}$ : $[0, 100]$

These intervals were justified using the following arguments.

First, because we observed a gradual modulation of the rate of transcription by both Bicoid and Runt in the middle region of the embryo we reasoned that the binding sites for these transcription factors were not saturated. As a result, we posited that the real dissociation constant should be between the minimum and maximum measured values of Bicoid and Runt (Figure 3—figure supplement 2). Our measurements of Bicoid and Runt concentration yield fluorescence values over the 0–100 AU range for the embryo region that we used for contrasting our model and experimental data (20–50% of the embryo length), such that the dissociation constants ( $K_{b}$ and $K_{r}$ ) should not exceed the maximum value of the Bicoid or Runt concentration.

Second, $ω_{b p}$ represents the cooperativity between Bicoid complex and RNAP. In the statistical mechanics framework, this cooperativity can be expressed using the interaction energy between Bicoid and RNAP, $Δ ϵ_{b p}$ , such that $ω_{b p} = e x p (- β Δ ϵ_{b p})$ , where $β = \frac{1}{k_{B} T}$ , $k_{B}$ is the Boltzmann constant and $T$ is the temperature. There is not much known about in vivo interaction energies between Bicoid and RNAP complex, thus we tried several different bounds until we found a narrow enough parameter bound with unimodal distribution of the posterior chain. As we could see from the corner plots in Figure 4C, there is a positive correlation between $K_{b}$ and $ω_{b p}$ . Thus, we constrained the $ω_{b p}$ intervals by finding an interval that gives both well-constrained $K_{b}$ and $ω_{b p}$ (Figure 4C).

Third, $R$ represents the rate of RNAP loading when the promoter is occupied, thus it is constrained by the maximum observed rate of RNAP loading (Figure 3—figure supplement 2).

Fourth, $p = ([R N A P] / K_{p})$ represents the concentration of RNAP divided by its dissociation constant. Recall that the predicted rate of transcription from hunchback P2 in the limit where the Bicoid concentration reaches zero is given by

R a t e ([B i c o i d] \to 0) = R \frac{p}{1 + p} .

This rate of transcription at the posterior region, where Bicoid reaches zero, is much lower than that at the anterior region where Bicoid saturates given by $R$ (Figure 3—figure supplement 2). As a result, we can write the inequality

R \frac{p}{1 + p} ≪ R .

such that

\frac{p}{1 + p} ≪ 1,

which holds if $p ≪ 1$ .

Finally, we did not have good estimates for the intervals of either Runt-Runt cooperativity, $ω_{r r}$ , or higher-order cooperativity, $ω_{r r p}$ . Thus, we initially started with an interval of $[0, 100]$ , of the same order as the interval we used $ω_{b p}$ . We then explored whether this parameter bound was sufficient to give us constrained values of $ω_{r r}$ and $ω_{r r p}$ . As we showed in Figure 6—figure supplement 6D, this interval gives reasonably constrained values of $ω_{r r}$ and $ω_{r r p}$ . As shown in Figure 6 and Figure 6—figure supplement 6, we posit that the $ω_{r r}$ parameter is not well-constrained not because of its width of the interval, but because it is not as essential for the model fit to the data as it is to include $ω_{r r p}$ into the model. Overall, our MCMC inference results as well as the corner plots shown demonstrate that our parameter intervals chosen were reasonable.

S6.1 Calculation of the Akaike Information Criterion

The Akaike Information Criterion (AIC) is a quantitative metric to assess how well a model performs compared to other models (Akaike, 1974). AIC is defined by two terms that account for: (1) how well the model explains the data (log-likelihood), and (2) how many parameters are used (to penalize the model). The AIC is expressed by

AIC = - 2 \log (likelihood) + 2 p,

where $p$ is the number of free parameters in the model. As a result, given a set of models, the model with minimal AIC value is preferred (Akaike, 1974).

First, we assume that our measurement error follows a Gaussian probability distribution. As a result, we can write down the probability of measuring a k-th datum, $x_{k}$ from a Gaussian distribution whose mean is $μ_{k}$ and standard deviation is $σ_{k}$ , termed $p r o b (x_{k} | μ_{k}, σ_{k})$ (Sivia and Skilling, 2006)

p r o b (x_{k} | μ_{k}, σ_{k}) = \frac{1}{σ_{k} \sqrt{2 π}} e x p (- \frac{{(x_{k} - μ_{k})}^{2}}{2 σ_{k}^{2}}) .

We can write the likelihood $L$ of a given set of $N$ measurements by multiplying the probabilities of each independent measurement such that

L = p r o b ({x_{1}, x_{2}, \dots, x_{N}} | {μ_{1}, μ_{2}, \dots, μ_{N}}, {σ_{1}, σ_{2}, \dots, σ_{N}}) = Π_{k = 1}^{N} \frac{1}{σ_{k} \sqrt{2 π}} e x p (- \frac{{(x_{k} - μ_{k})}^{2}}{2 σ_{k}^{2}}),

where x_k is the k-th datum with standard deviation of $σ_{k}$ from an expected value of $μ_{k}$ from a model (Sivia and Skilling, 2006).

By taking the logarithm of the likelihood (L) we get

\log (L) = - \sum_{k = 1}^{N} \log (σ_{k} \sqrt{2 π}) - \sum_{k = 1}^{N} \frac{(x_{k} - μ_{k})^{2}}{2 σ_{k}^{2}} .

We can now plug in this expression for the log-likelihood into the definition of the AIC from Equation S19, resulting in

AIC = 2 k + 2 \sum_{k = 1}^{N} (σ_{k} \sqrt{2 π}) + 2 \sum_{k = 1}^{N} \frac{(x_{k} - μ_{k})^{2}}{2 σ_{k}^{2}} .

All of the Akaike Information Criteria values shown in Figure 6—figure supplement 1 and Figure 5—figure supplement 3 are calculated using Equation S23.

S7 Comparison of different modes of repression

Transcriptional repressors have been classified into two broad categories: short-range and long-range, depending on the genomic length scale that they act on Courey and Jia, 2001; Li and Gilmour, 2011. Long-range repression is realized by the recruitment of chromatin modifiers. In contrast, short-range repressors act within 100–150 bp by interacting with nearby transcription factors or with the promoter (Li and Gilmour, 2011). Traditionally, the molecular mechanism of short-range repressors, such as Runt, have been further classified into three categories: “direct repression”, “competition”, and “quenching” (Gray et al., 1994; Jaynes and O’Farrell, 1991; Arnosti et al., 1996; Kulkarni and Arnosti, 2005). In “direct repression”, the repressor inhibits the binding of RNAP to the promoter (Figure 5—figure supplement 1A). “Competition” denotes a repressor that competes with an activator for the same DNA binding location (Figure 5—figure supplement 1B). This molecular mechanisms has been proposed for the action of Giant and Krüppel repressors on the even-skipped stripe 2 enhancer, where some activator and repressor binding sites partially overlap (Small et al., 1992). Lastly, “quenching” corresponds to the case where the repressor and activator do not interact with each other directly. Instead, the repressor inhibits the activators’ action of recruiting the RNAP (Figure 5—figure supplement 1C).

Despite several classic studies of the molecular mechanism of repressors in the early fly embryo (Gray et al., 1994; Ip et al., 1992; Bothma et al., 2011; Jaynes and O’Farrell, 1991), the mechanisms of many repressors remain unknown. Note that, even for the same repressor, the mode of repression might not be the same depending on, for example, its sequence context (Koromila and Stathopoulos, 2019; Hang and Gergen, 2017). For example, it has been proposed that Runt repressor acts with different mechanisms in different regulatory elements of the sloppy-paired gene (Hang and Gergen, 2017). In this section, we derive a thermodynamic model from each mode of repression and compare their explanatory power in the context of our data stemming from the hunchback P2 enhancer containing one Runt binding site. Note that, in the main text, we already developed a thermodynamic model for the “direct repression” scenario (Section ‘Derivation of the general and simpler thermodynamic model for the hunchback P2 enhancer with one Runt binding site’). As a result, in this section, we focus on deriving the thermodynamic models for the “competition” and “quenching” scenarios, but repeat the result of the derivation for the “direct repression” here for ease comparison between different models.

S7.1 Derivation of models for each scenario of repression for hunchback P2 with one Runt binding site

S7.1.1 Modeling repression for hunchback P2 with one Runt binding site: direct repression

For completeness, we repeat the expression for the direct repression scenario as shown in Section ‘Derivation of the general and simpler thermodynamic model for the hunchback P2 enhancer with one Runt binding site’ and Figure 5—figure supplement 1A. The probability of finding RNAP bound to the promoter, $p_{b o u n d}$ , is calculated by dividing the sum of all statistical weights featuring RNAP by the sum of the weights of all possible microstates. The calculation of $p_{b o u n d}$ , combined with Equation 1, leads to the expression

R a t e = R p_{b o u n d} = R \frac{p + b^{6} p ω_{b p} + r p ω_{r p} + b^{6} r p ω_{b p} ω_{r p}}{1 + b^{6} + r + b^{6} r + p + b^{6} p ω_{b p} + r p ω_{r p} + b^{6} r p ω_{b p} ω_{r p}},

where the parameters are as defined in Figure 2.

S7.1.2 Modeling repression for hunchback P2 with one Runt binding site: competition

In the competition scenario, Runt binding makes Bicoid binding less likely. This mechanism can be captured by an interaction term between Bicoid and Runt given by $ω_{b r}$ . Building on our assumption of strong Bicoid-Bicoid cooperativity, we posit that Runt disfavors the state with six bound Bicoid molecules. We can enumerate the states and weights from Figure 5—figure supplement 1B to calculate the $R a t e$ ( $\propto p_{b o u n d}$ ), which leads to

R a t e = R \frac{p + b^{6} p ω_{b p} + r p + b^{6} r p ω_{b p} ω_{b r}}{1 + p + b^{6} + r + b^{6} r ω_{b r} + b^{6} p ω_{b p} + r p + b^{6} r p ω_{b p} ω_{b r}} .

S7.1.3 Modeling repression for hunchback P2 with one Runt binding site: quenching

In the quenching scenario, Runt reduces the magnitude of the cooperativity between the Bicoid complex and RNAP by a factor $ω_{b r p}$ . We can enumerate the states and weights from Figure 5—figure supplement 1C, leading to a rate of transcription given by

R a t e = R \frac{p + b^{6} p ω_{b p} + r p + b^{6} r p ω_{b p} ω_{b r p}}{1 + p + b^{6} + r + b^{6} r + b^{6} p ω_{b p} + r p + b^{6} r p ω_{b p} ω_{b r p}} .

With these expressions for each repression mechanism in hand, we can now compare how each model fares against our experimental data.

S7.2 Comparing the three models of repression with the one-Runt binding site data

We used the MCMC sampling to fit each model to our experimentally measured initial rate of transcription over the anterior-posterior axis of the embryo. As shown in Figure 5—figure supplement 2A, B and C, we see that all three models can explain the [100] and [010] construct data relatively well. However, the competition model resulted in a qualitatively poor fit to the [001] construct as shown by the lack of saturation in the most anterior region of the embryo (Figure 5—figure supplement 2C), (ii). The direct repression and quenching models showed equally good fits to the data stemming from this construct.

S7.3 Predicting two-Runt binding sites data for each mode of repression

We further tested these different models of repression by using the parameters inferred from the one-Runt binding site constructs to predict the rate of initiation for the two-Runt binding sites constructs. As reasoned in the main text, we began by assuming that the two Runt molecules act independently of each other such that there are no interactions between Runt molecules. Figure 6—figure supplement 1 shows this parameter-free prediction for our two-Runt binding sites constructs for all three modes of repression. As shown in the figure, none of the models can explain the data, suggesting the need to invoke additional interactions between the molecular players of our model.

Next, we considered whether Runt-Runt pairwise or higher-order cooperativities had to be invoked in order to explain the two-Runt binding sites data for both the competition and quenching mechanisms. For the competition model, we considered Runt-Runt cooperativity, $ω_{r r}$ , and Runt-Runt-Bicoid higher-order cooperativity, $ω_{b r r}$ in addition to the Runt-Bicoid interaction term $ω_{b r}$ . In the quenching scenario, we accounted for Runt-Runt cooperativity, $ω_{r r}$ , and Runt-Runt-Bicoid-RNAP higher-order cooperativity, $ω_{b r r p}$ . For both the competition (Figure 6—figure supplement 2) and quenching (Figure 6—figure supplement 3) mechanisms, we observed a qualitatively similar trend to that observed for direct repression (Figure 6). Specifically, as shown in Figure 6—figure supplements 2C and 3C, considering pairwise cooperativity did not significantly improve the MCMC fits to the data for either model considered. Further, considering only the higher-order cooperativity also did not improve the fits for both competition and quenching mechanisms as shown in Figure 6—figure supplement 2D and Figure 6—figure supplement 3D. Invoking both Runt-Runt cooperativity and higher-order cooperativity improved the fits qualitatively for both competition and quenching mechanisms as shown in Figure 6—figure supplement 2E and Figure 6—figure supplement 3E.

While the quenching model showed almost equally good MCMC fits to the data as the direct repression model, the competition model showed qualitatively poor fits in any combination of cooperativities. In particular, there was a significant mismatch in the most anterior region of the embryo, where Bicoid is thought to saturate hunchback expression. While we do not view these fits as conclusive evidence to support one mechanism over the other, an exercise that would require a new round of experimentation, we conclude that higher-order cooperativity is required to explain the data from the two-Runt binding sites constructs regardless of the choice of mechanism of Runt.

S8 Quantifying the nuclear concentration of LlamaTag-Runt

The major caveat in the eGFP:LlamaTag-Runt fluorescence measurements is that the raw nuclear fluorescence that we measured consists of two populations: eGFP bound to the LlamaTag-Runt, and free, unbound eGFP. Thus, in order to measure nuclear Runt concentration, we need to factor out the contribution from free eGFP to the overall fluorescence.

We followed the procedure described in Bothma et al., 2018 which consists of using cytoplasmic fluorescence to calculate the free nuclear eGFP under two assumptions. First, we posit that most of the transcription factors reside in the nucleus such that the cytoplasmic fluorescence mostly reports on free cytoplasmic eGFP. Second, we assume that the nucleus-to-cytoplasm ratio of free eGFP is kept constant at a measured chemical equilibrium of $K_{G} = G F P_{C} / G F P_{N} = 0.8$ , where $G F P_{C}$ and $G F P_{N}$ are the eGFP fluorescence in nuclei and cytoplasm in the absence of LlamaTag (Bothma et al., 2018).

As shown in Bothma et al., 2018, the nuclear concentration of the GFP-tagged transcription factor, $G F P - T F_{N}$ , is given by

G F P - T F_{N} = F l u o_{N} - \frac{F l u o_{C}}{K_{G}},

where $F l u o_{N}$ and $F l u o_{C}$ are the eGFP fluorescence in nuclei and cytoplasm, respectively, that we measured in the embryos with both eGFP and LlamaTagged Runt. The resulting nuclear concentration of LlamaTag-Runt is shown in Figure 3B.

Data availability

All data (both input transcription factor concentration and output transcription from all synthetic enhancers, both pre- and post-processed data) have been deposited in Dryad under the DOI https://doi.org/10.5061/dryad.7sqv9s4sv.

The following data sets were generated

1. Kim Y
2. Rhee K
3. Liu J
4. Jeammet S
5. Turner M
6. Small S
7. Garcia HG
(2021) Dryad Digital Repository
Predictive modeling reveals that higher-order cooperativity drives transcriptional repression in a synthetic developmental enhancer.
https://doi.org/10.5061/dryad.7sqv9s4sv

References

(1982) Quantitative model for gene regulation by lambda phage repressor
PNAS 79:1129–1133.
https://doi.org/10.1073/pnas.79.4.1129
- PubMed
- Google Scholar
1. Akaike H
(1974) A new look at the statistical model identification
IEEE Transactions on Automatic Control 19:716–723.
https://doi.org/10.1109/TAC.1974.1100705
- Google Scholar
Book
1. Alberts B
2. Johnson A
3. Lewis J
4. Raff M
5. Roberts K
6. Walter P
(2015)
Molecular Biology of the Cell

New York, NY: Taylor and Francis Group.
- Google Scholar
1. Amit R
2. Garcia HG
3. Phillips R
4. Fraser SE
(2011) Building enhancers from the ground up: a synthetic biology approach
Cell 146:105–118.
https://doi.org/10.1016/j.cell.2011.06.024
- PubMed
- Google Scholar
1. Arnosti DN
2. Barolo S
3. Levine M
4. Small S
(1996) The eve stripe 2 enhancer employs multiple modes of transcriptional synergy
Development 122:205–214.
https://doi.org/10.1242/dev.122.1.205
- PubMed
- Google Scholar
1. Barr KA
2. Martinez C
3. Moran JR
4. Kim AR
5. Ramos AF
6. Reinitz J
(2017) Synthetic enhancer design by in silico compensatory evolution reveals flexibility and constraint in cis-regulation
BMC Systems Biology 11:116.
https://doi.org/10.1186/s12918-017-0485-2
- PubMed
- Google Scholar
1. Bateman JR
2. Lee AM
3. Wu C
(2006) Site-Specific transformation of Drosophila via phiC31 integrase-mediated cassette exchange
Genetics 173:769–777.
https://doi.org/10.1534/genetics.106.056945
- PubMed
- Google Scholar
1. Bertrand E
2. Chartrand P
3. Schaefer M
4. Shenoy SM
5. Singer RH
6. Long RM
(1998) Localization of ash1 mRNA particles in living yeast
Molecular Cell 2:437–445.
https://doi.org/10.1016/s1097-2765(00)80143-4
- PubMed
- Google Scholar
Preprint
(2020) Allosteric Conformational Ensembles Have Unlimited Capacity for Integrating Information
bioRxiv.
https://doi.org/10.1101/2020.12.10.420117
- Google Scholar
1. Bintu L
2. Buchler NE
3. Garcia HG
4. Gerland U
5. Hwa T
6. Kondev J
7. Kuhlman T
8. Phillips R
(2005a) Transcriptional regulation by the numbers: applications
Current Opinion in Genetics & Development 15:125–135.
https://doi.org/10.1016/j.gde.2005.02.006
- PubMed
- Google Scholar
1. Bintu L
2. Buchler NE
3. Garcia HG
4. Gerland U
5. Hwa T
6. Kondev J
7. Phillips R
(2005b) Transcriptional regulation by the numbers: models
Current Opinion in Genetics & Development 15:116–124.
https://doi.org/10.1016/j.gde.2005.02.007
- Google Scholar
(2013a) DNA sequence-dependent mechanics and protein-assisted bending in repressor-mediated loop formation
Physical Biology 10:066005.
https://doi.org/10.1088/1478-3975/10/6/066005
- PubMed
- Google Scholar
(2013b) Theoretical and experimental dissection of DNA loop-mediated repression
Physical Review Letters 110:018101.
https://doi.org/10.1103/PhysRevLett.110.018101
- PubMed
- Google Scholar
1. Bolouri H
2. Davidson EH
(2003) Transcriptional regulatory cascades in development: initial rates, not steady state, determine network kinetics
PNAS 100:9371–9376.
https://doi.org/10.1073/pnas.1533293100
- PubMed
- Google Scholar
(2011) The snail repressor inhibits release, not elongation, of paused Pol II in the Drosophila embryo
Current Biology 21:1571–1577.
https://doi.org/10.1016/j.cub.2011.08.019
- PubMed
- Google Scholar
(2018) LlamaTags: a versatile tool to image transcription factor dynamics in live embryos
Cell 173:1810–1822.
https://doi.org/10.1016/j.cell.2018.03.069
- PubMed
- Google Scholar
(2014) The transcription factor titration effect dictates level of gene expression
Cell 156:1312–1323.
https://doi.org/10.1016/j.cell.2014.02.022
- PubMed
- Google Scholar
1. Briscoe J
2. Small S
(2015) Morphogen rules: design principles of gradient-mediated embryo patterning
Development 142:3996–4009.
https://doi.org/10.1242/dev.129452
- PubMed
- Google Scholar
(2003) On schemes of combinatorial transcription logic
PNAS 100:5136–5141.
https://doi.org/10.1073/pnas.0930314100
- PubMed
- Google Scholar
1. Chen H
2. Xu Z
3. Mei C
4. Yu D
5. Small S
(2012) A system of repressor gradients spatially organizes the boundaries of bicoid-dependent target genes
Cell 149:618–629.
https://doi.org/10.1016/j.cell.2012.03.018
- PubMed
- Google Scholar
1. Courey AJ
2. Jia S
(2001) Transcriptional repression: the long and the short of it
Genes & Development 15:2786–2796.
https://doi.org/10.1101/gad.939601
- PubMed
- Google Scholar
(2016) Quantitatively predictable control of Drosophila transcriptional enhancers in vivo with engineered transcription factors
Nature Genetics 48:292–298.
https://doi.org/10.1038/ng.3509
- PubMed
- Google Scholar
1. Dodd IB
2. Shearwin KE
3. Perkins AJ
4. Burr T
5. Hochschild A
6. Egan JB
(2004) Cooperativity in long-range gene regulation by the lambda cI repressor
Genes & Development 18:344–354.
https://doi.org/10.1101/gad.1167904
- PubMed
- Google Scholar
1. Driever W
2. Nüsslein-Volhard C
(1988) The bicoid protein determines position in the Drosophila embryo in a concentration-dependent manner
Cell 54:95–104.
https://doi.org/10.1016/0092-8674(88)90183-3
- PubMed
- Google Scholar
1. Driever W
2. Nüsslein-Volhard C
(1989) The bicoid protein is a positive regulator of hunchback transcription in the early Drosophila embryo
Nature 337:138–143.
https://doi.org/10.1038/337138a0
- PubMed
- Google Scholar
(1989) Determination of spatial domains of zygotic gene expression in the Drosophila embryo by the affinity of binding sites for the bicoid morphogen
Nature 340:363–367.
https://doi.org/10.1038/340363a0
- PubMed
- Google Scholar
1. Dufourt J
2. Trullo A
3. Hunter J
4. Fernandez C
5. Lazaro J
6. Dejean M
7. Morales L
8. Nait-Amer S
9. Schulz KN
10. Harrison MM
11. Favard C
12. Radulescu O
13. Lagha M
(2018) Temporal control of gene expression by the pioneer factor Zelda through transient interactions in hubs
Nature Communications 9:5194.
https://doi.org/10.1038/s41467-018-07613-z
- PubMed
- Google Scholar
(2020) Quantitative dissection of transcription in development yields evidence for transcription-factor-driven chromatin accessibility
eLife 9:e56429.
https://doi.org/10.7554/eLife.56429
- PubMed
- Google Scholar
(2016a) SiteOut: an online tool to design binding site-free DNA sequences
PLOS ONE 11:e0151740.
https://doi.org/10.1371/journal.pone.0151740
- PubMed
- Google Scholar
(2016b) Information integration and energy expenditure in gene regulation
Cell 166:234–244.
https://doi.org/10.1016/j.cell.2016.06.012
- PubMed
- Google Scholar
1. Fakhouri WD
2. Ay A
3. Sayal R
4. Dresch J
5. Dayringer E
6. Arnosti DN
(2010) Deciphering a transcriptional regulatory code: modeling short-range repression in the Drosophila embryo
Molecular Systems Biology 6:341.
https://doi.org/10.1038/msb.2009.97
- PubMed
- Google Scholar
1. Garcia HG
2. Phillips R
(2011) Quantitative dissection of the simple repression input-output function
PNAS 108:12173–12178.
https://doi.org/10.1073/pnas.1015616108
- PubMed
- Google Scholar
1. Garcia HG
2. Sanchez A
3. Boedicker JQ
4. Osborne M
5. Gelles J
6. Kondev J
7. Phillips R
(2012) Operator sequence alters gene expression independently of transcription factor occupancy in bacteria
Cell Reports 2:150–161.
https://doi.org/10.1016/j.celrep.2012.06.004
- PubMed
- Google Scholar
1. Garcia HG
2. Tikhonov M
3. Lin A
4. Gregor T
(2013) Quantitative imaging of transcription in living Drosophila embryos links polymerase activity to patterning
Current Biology 23:2140–2145.
https://doi.org/10.1016/j.cub.2013.08.054
- PubMed
- Google Scholar
1. Garcia HG
2. Berrocal A
3. Kim YJ
4. Martini G
5. Zhao J
(2020) Lighting up the central dogma for predictive developmental biology
Current Topics in Developmental Biology 137:1–35.
https://doi.org/10.1016/bs.ctdb.2019.10.010
- PubMed
- Google Scholar
Software
1. Garcia Lab @ UC Berkeley
(2022) MRNADynamics, version 1ca47dd
Github.

https://github.com/GarciaLab/mRNADynamics
1. Gergen JP
2. Butler BA
(1988) Isolation of the Drosophila segmentation gene runt and analysis of its expression during embryogenesis
Genes & Development 2:1179–1193.
https://doi.org/10.1101/gad.2.9.1179
- PubMed
- Google Scholar
1. Geyer CJ
2. Thompson EA
(1992) Constrained monte carlo maximum likelihood for dependent data
Journal of the Royal Statistical Society 54:657–683.
https://doi.org/10.1111/j.2517-6161.1992.tb01443.x
- Google Scholar
1. Giniger E
2. Ptashne M
(1988) Cooperative DNA binding of the yeast transcriptional activator GAL4
PNAS 85:382–386.
https://doi.org/10.1073/pnas.85.2.382
- PubMed
- Google Scholar
(2015) Crispr-cas9 genome editing in Drosophila
Current Protocols in Molecular Biology 111:31.
https://doi.org/10.1002/0471142727.mb3102s111
- PubMed
- Google Scholar
(1994) Short-Range repression permits multiple enhancers to function autonomously within a complex promoter
Genes & Development 8:1829–1838.
https://doi.org/10.1101/gad.8.15.1829
- PubMed
- Google Scholar
(2007) Probing the limits to positional information
Cell 130:153–164.
https://doi.org/10.1016/j.cell.2007.05.025
- PubMed
- Google Scholar
(2001) An adaptive metropolis algorithm
Bernoulli 7:223.
https://doi.org/10.2307/3318737
- Google Scholar
1. Haario H
2. Laine M
3. Mira A
4. Saksman E
(2006) Dram: efficient adaptive MCMC
Statistics and Computing 16:339–354.
https://doi.org/10.1007/s11222-006-9438-0
- Google Scholar
1. Hang S
2. Gergen JP
(2017) Different modes of enhancer-specific regulation by runt and even-skipped during Drosophila segmentation
Molecular Biology of the Cell 28:681–691.
https://doi.org/10.1091/mbc.E16-09-0630
- PubMed
- Google Scholar
(2010) Grainyhead and Zelda compete for binding to the promoters of the earliest-expressed Drosophila genes
Developmental Biology 345:248–255.
https://doi.org/10.1016/j.ydbio.2010.06.026
- PubMed
- Google Scholar
1. He HH
2. Meyer CA
3. Shin H
4. Bailey ST
5. Wei G
6. Wang Q
7. Zhang Y
8. Xu K
9. Ni M
10. Lupien M
11. Mieczkowski P
12. Lieb JD
13. Zhao K
14. Brown M
15. Liu XS
(2010) Nucleosome dynamics define transcriptional enhancers
Nature Genetics 42:343–347.
https://doi.org/10.1038/ng.545
- PubMed
- Google Scholar
(1990) Identification of consensus patterns in unaligned DNA sequences known to be functionally related
Bioinformatics 6:81–92.
https://doi.org/10.1093/bioinformatics/6.2.81
- Google Scholar
1. Hertz GZ
2. Stormo GD
(1999) Identifying DNA and protein patterns with statistically significant alignments of multiple sequences
Bioinformatics 15:563–577.
https://doi.org/10.1093/bioinformatics/15.7.563
- PubMed
- Google Scholar
1. Hewitt GF
2. Strunk BS
3. Margulies C
4. Priputin T
5. Wang XD
6. Amey R
7. Pabst BA
8. Kosman D
9. Reinitz J
10. Arnosti DN
(1999) Transcriptional repression by the Drosophila giant protein: cis element positioning provides an alternative means of interpreting an effector gradient
Development 126:1201–1210.
https://doi.org/10.1242/dev.126.6.1201
- PubMed
- Google Scholar
1. Ip YT
2. Park RE
3. Kosman D
4. Yazdanbakhsh K
5. Levine M
(1992) Dorsal-twist interactions establish snail expression in the presumptive mesoderm of the Drosophila embryo
Genes & Development 6:1518–1530.
https://doi.org/10.1101/gad.6.8.1518
- PubMed
- Google Scholar
1. Jaynes JB
2. O’Farrell PH
(1991) Active repression of transcription by the engrailed homeodomain protein
The EMBO Journal 10:1427–1433.
https://doi.org/10.1002/j.1460-2075.1991.tb07663.x
- PubMed
- Google Scholar
1. Jiang J
2. Levine M
(1993) Binding affinities and cooperative interactions with bHLH activators delimit threshold responses to the dorsal gradient morphogen
Cell 72:741–752.
https://doi.org/10.1016/0092-8674(93)90402-c
- PubMed
- Google Scholar
1. Koromila T
2. Stathopoulos A
(2017) Broadly expressed repressors integrate patterning across orthogonal axes in embryos
PNAS 114:8295–8300.
https://doi.org/10.1073/pnas.1703001114
- PubMed
- Google Scholar
1. Koromila T
2. Stathopoulos A
(2019) Distinct roles of broadly expressed repressors support dynamic enhancer action and change in time
Cell Reports 28:855–863.
https://doi.org/10.1016/j.celrep.2019.06.063
- PubMed
- Google Scholar
1. Kulkarni MM
2. Arnosti DN
(2005) Cis-Regulatory logic of short-range transcriptional repression in Drosophila melanogaster
Molecular and Cellular Biology 25:3411–3420.
https://doi.org/10.1128/MCB.25.9.3411-3420.2005
- PubMed
- Google Scholar
Software
1. Laine M
(2021) Mcmcstat, version swh:1:rev:ca0cbd288f03c7a29050b3d6698a96b45ccfa4b2
Software Heritage.

https://archive.softwareheritage.org/swh:1:dir:6dfd72aa7998045fb1e4725081b8ffdefba354c7;origin=https://github.com/mjlaine/mcmcstat;visit=swh:1:snp:d91631a580497fff421f83ef3e706159a0d67602;anchor=swh:1:rev:ca0cbd288f03c7a29050b3d6698a96b45ccfa4b2
1. Lammers NC
2. Galstyan V
3. Reimer A
4. Medin SA
5. Wiggins CH
6. Garcia HG
(2020) Multimodal transcriptional control of pattern formation in embryonic development
PNAS 117:836–847.
https://doi.org/10.1073/pnas.1912500117
- PubMed
- Google Scholar
1. Lebrecht D
2. Foehr M
3. Smith E
4. Lopes FJP
5. Vanario-Alonso CE
6. Reinitz J
7. Burz DS
8. Hanes SD
(2005) Bicoid cooperative DNA binding is critical for embryonic patterning in Drosophila
PNAS 102:13176–13181.
https://doi.org/10.1073/pnas.0506462102
- PubMed
- Google Scholar
(1999) Core-Binding factor influences the disease specificity of Moloney murine leukemia virus
Journal of Virology 73:5535–5547.
https://doi.org/10.1128/JVI.73.7.5535-5547.1999
- PubMed
- Google Scholar
1. Li J
2. Gilmour DS
(2011) Promoter proximal pausing and the control of gene expression
Current Opinion in Genetics & Development 21:231–235.
https://doi.org/10.1016/j.gde.2011.01.010
- PubMed
- Google Scholar
1. Li C
2. Cesbron F
3. Oehler M
4. Brunner M
5. Höfer T
(2018) Frequency modulation of transcriptional bursting enables sensitive and rapid gene regulation
Cell Systems 6:409–423.
https://doi.org/10.1016/j.cels.2018.01.012
- PubMed
- Google Scholar
Preprint
1. Liu J
2. Hansen D
3. Eck E
4. Kim YJ
5. Turner M
6. Alamos S
7. Garcia HG
(2021) Real-Time Single-Cell Characterization of the Eukaryotic Transcription Cycle Reveals Correlations between RNA Initiation, Elongation, and Cleavage
bioRxiv.
https://doi.org/10.1101/2020.08.29.273474
- Google Scholar
1. Löhr U
2. Chung HR
3. Beller M
4. Jäckle H
(2009) Antagonistic action of bicoid and the repressor capicua determines the spatial limits of Drosophila head gene expression domains
PNAS 106:21695–21700.
https://doi.org/10.1073/pnas.0910225106
- PubMed
- Google Scholar
1. Lott SE
2. Villalta JE
3. Schroth GP
4. Luo S
5. Tonkin LA
6. Eisen MB
(2011) Noncanonical compensation of zygotic X transcription in early Drosophila melanogaster development revealed through single-embryo RNA-seq
PLOS Biology 9:e1000590.
https://doi.org/10.1371/journal.pbio.1000590
- PubMed
- Google Scholar
1. Louder RK
2. He Y
3. López-Blanco JR
4. Fang J
5. Chacón P
6. Nogales E
(2016) Structure of promoter-bound TFIID and model of human pre-initiation complex assembly
Nature 531:604–609.
https://doi.org/10.1038/nature17394
- PubMed
- Google Scholar
(2013) Live imaging of bicoid-dependent transcription in Drosophila embryos
Current Biology 23:2135–2139.
https://doi.org/10.1016/j.cub.2013.08.053
- PubMed
- Google Scholar
1. Ma D
2. Alberti M
3. Lynch C
4. Nikaido H
5. Hearst JE
(1996) The local repressor AcrR plays a modulating role in the regulation of acrab genes of Escherichia coli by global stress signals
Molecular Microbiology 19:101–112.
https://doi.org/10.1046/j.1365-2958.1996.357881.x
- PubMed
- Google Scholar
(2010) Drawing an elephant with four complex parameters
American Journal of Physics 78:648–649.
https://doi.org/10.1119/1.3254017
- Google Scholar
1. Melnikova IN
2. Crute BE
3. Wang S
4. Speck NA
(1993) Sequence specificity of the core-binding factor
Journal of Virology 67:2408–2411.
https://doi.org/10.1128/JVI.67.4.2408-2411.1993
- PubMed
- Google Scholar
(2022) Near-infrared nuclear markers for Drosophila imaging
MicroPubl Biol 2022:000639.
https://doi.org/10.17912/micropub.biology.000639
- Google Scholar
1. Papagianni A
2. Forés M
3. Shao W
4. He S
5. Koenecke N
6. Andreu MJ
7. Samper N
8. Paroush Z
9. González-Crespo S
10. Zeitlinger J
11. Jiménez G
(2018) Capicua controls toll/IL-1 signaling targets independently of RTK regulation
PNAS 115:1807–1812.
https://doi.org/10.1073/pnas.1713930115
- PubMed
- Google Scholar
1. Park J
2. Estrada J
3. Johnson G
4. Vincent BJ
5. Ricci-Tam C
6. Bragdon MD
7. Shulgina Y
8. Cha A
9. Wunderlich Z
10. Gunawardena J
11. DePace AH
(2019) Dissecting the sharp response of a canonical developmental enhancer reveals multiple sources of cooperativity
eLife 8:e41266.
https://doi.org/10.7554/eLife.41266
- PubMed
- Google Scholar
1. Parker DS
2. White MA
3. Ramos AI
4. Cohen BA
5. Barolo S
(2011) The cis-regulatory logic of hedgehog gradient responses: key roles for gli binding affinity, competition, and cooperativity
Science Signaling 4:176.
https://doi.org/10.1126/scisignal.2002077
- PubMed
- Google Scholar
(2013) A thermodynamic model of the cooperative interaction between the archaeal transcription factor ss-lrpb and its tripartite operator DNA
Gene 524:330–340.
https://doi.org/10.1016/j.gene.2013.03.118
- PubMed
- Google Scholar
Book
(2009)
Physical Biology of the Cell

New York, United States: Garland Science.
- Google Scholar
(2019) Figure 1 theory meets figure 2 experiments in the study of gene expression
Annual Review of Biophysics 48:121–163.
https://doi.org/10.1146/annurev-biophys-052118-115525
- PubMed
- Google Scholar
Book
1. Ptashne M
2. Gann A
(2002)
Genes and Signals

New York: Cold Spring Harbor Laboratory Press.
- Google Scholar
Book
1. Ptashne M
(2004)
A Genetic Switch: Phage Lambda Revisited

Cold Spring Harbor, N.Y: Cold Spring Harbor Laboratory Press.
- Google Scholar
1. Razo-Mejia M
2. Barnes SL
3. Belliveau NM
4. Chure G
5. Einav T
6. Lewis M
7. Phillips R
(2018) Tuning transcriptional regulation through signaling: a predictive theory of allosteric induction
Cell Systems 6:456–469.
https://doi.org/10.1016/j.cels.2018.02.004
- PubMed
- Google Scholar
1. Razo-Mejia M
2. Marzen S
3. Chure G
4. Taubman R
5. Morrison M
6. Phillips R
(2020) First-principles prediction of the information processing capacity of a simple genetic circuit
Physical Review. E 102:022404.
https://doi.org/10.1103/PhysRevE.102.022404
- PubMed
- Google Scholar
Preprint
1. Reimer A
2. Alamos S
3. Westrum C
4. Turner MA
5. Talledo P
6. Zhao J
7. Garcia HG
(2021) Minimal Synthetic Enhancers Reveal Control of the Probability of Transcriptional Engagement and Its Timing by a Morphogen Gradient
bioRxiv.
https://doi.org/10.1101/2021.07.10.451524
- Google Scholar
(2019) REDfly: the transcriptional regulatory element database for Drosophila
Nucleic Acids Research 47:D828–D834.
https://doi.org/10.1093/nar/gky957
- PubMed
- Google Scholar
1. Sauer F
2. Jäckle H
(1991) Concentration-Dependent transcriptional activation or repression by Krüppel from a single binding site
Nature 353:563–566.
https://doi.org/10.1038/353563a0
- PubMed
- Google Scholar
1. Sayal R
2. Dresch JM
3. Pushel I
4. Taylor BR
5. Arnosti DN
(2016) Quantitative perturbation-based analysis of gene expression predicts enhancer activity in early Drosophila embryo
eLife 5:e08445.
https://doi.org/10.7554/eLife.08445
- PubMed
- Google Scholar
(2008) Predicting expression patterns from regulatory sequence in Drosophila segmentation
Nature 451:535–540.
https://doi.org/10.1038/nature06496
- PubMed
- Google Scholar
1. Sepúlveda LA
2. Xu H
3. Zhang J
4. Wang M
5. Golding I
(2016) Measurement of gene regulation in individual cells reveals rapid switching between promoter states
Science 351:1218–1222.
https://doi.org/10.1126/science.aad0635
- PubMed
- Google Scholar
Book
1. Sivia D
2. Skilling J
(2006)
Data Analysis: A Bayesian Tutorial

Oxford, United Kingdom: Oxford University Press.
- Google Scholar
(1992) Regulation of even-skipped stripe 2 in the Drosophila embryo
The EMBO Journal 11:4047–4057.
https://doi.org/10.1002/j.1460-2075.1992.tb05498.x
- PubMed
- Google Scholar
1. Small S
2. Arnosti DN
(2020) Transcriptional enhancers in Drosophila
Genetics 216:1370.
https://doi.org/10.1534/genetics.120.301370
- Google Scholar
1. Spitz F
2. Furlong EEM
(2012) Transcription factors: from enhancer binding to developmental control
Nature Reviews. Genetics 13:613–626.
https://doi.org/10.1038/nrg3207
- PubMed
- Google Scholar
1. Vilar JMG
2. Leibler S
(2003) Dna looping and physical constraints on transcription regulation
Journal of Molecular Biology 331:981–989.
https://doi.org/10.1016/s0022-2836(03)00764-2
- PubMed
- Google Scholar
(2016) The appeasement of doug: a synthetic approach to enhancer biology
Integrative Biology 8:475–484.
https://doi.org/10.1039/c5ib00321k
- PubMed
- Google Scholar
(2011) Hairless is a cofactor for runt-dependent transcriptional regulation
Molecular Biology of the Cell 22:1364–1374.
https://doi.org/10.1091/mbc.E10-06-0483
- PubMed
- Google Scholar
1. Wu B
2. Miskolci V
3. Sato H
4. Tutucci E
5. Kenworthy CA
6. Donnelly SK
7. Yoon YJ
8. Cox D
9. Singer RH
10. Hodgson L
(2015) Synonymous modification results in high-fidelity gene expression of repetitive protein and nucleotide sequences
Genes & Development 29:876–886.
https://doi.org/10.1101/gad.259358.115
- PubMed
- Google Scholar

Decision letter

Jie Xiao

Reviewing Editor; Johns Hopkins University, United States
Naama Barkai

Senior Editor; Weizmann Institute of Science, Israel

Our editorial process produces two outputs: (i) public reviews designed to be posted alongside the preprint for the benefit of readers; (ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.

Decision letter after peer review:

Thank you for submitting your article "Predictive modeling reveals that higher-order cooperativity drives transcriptional repression in a synthetic developmental enhancer" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Naama Barkai as the Senior Editor. The reviewers have opted to remain anonymous.

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission. While the work is potentially important, the uniqueness of the model and experimental verification of model predictions need to be further justified.

The essential revisions are:

1. Modeling: Please provide a better articulation/construction of the model regarding the fitting parameters (which to keep constant and which to vary, the number of parameters to use, error bars in the experimental data), whether an alternative model should be considered, whether the provided model can indeed account quantitatively observed data (not just to demonstrate that previous models do not work), and whether statistically the model fits experimental data adequately (not by overfitting, examine using Akaike Information Criterion AIC or Bayesian Information Criterion BIC). See reviewers #1 pt 1, #2 pt 1-5, and #3's overall review, and pt1.

2. Experiments: please provide experimental measurement for a few conditions where the outcome of transcription regulation by Runt is measured by smFISH or protein concentrations to demonstrate that the measured promoter activity using RNAP loading rate is equivalent or better (without complications in mRNA/protein degradation).

Reviewer #1 (Recommendations for the authors):

1. Modeling: I would like to see a better articulation of what the alternative is to a 'thermodynamic model.' I understand the authors meaning to be equilibrium models arrived at through partition functions. But kinetic models based on detailed balance are also thermodynamic. I don't want to get twisted around in the language, but the over-emphasis on the ideological purity of thermodynamics was distracting and seems to refer to an ongoing dialog / conversation which is likely lost on the general reader.

a. One example is when they assume that the Runt dissociation constant remains unchanged with different positions (K_r) because the sequence is the same. The problem stems from their previous section where they found that there was an unsuspected relationship between the position of the Runt Binding sites and the parameter "p=[RNAP]/K_p" --- the point being that the sequence is the same in each construct, and therefore, based on their model K_p is changing. So why is it preferrable to assume that K_r is constant with position, but K_p changes?

b. Expanding on the previous point, they then proceed to state that an unchanging K_r is supported by the fact that by only varying the w_rp they were able to fit their data. Considering the similar role of K_r with w_rp within the equation, I wonder if one could reach the same result varying the K_r and holding w_rp constant? Furthermore, who's to say that both don't vary between the constructs? Could the authors address this concern by varying K_r and further discussing this possibility and what it means for the interpretations of the model? In summary, the authors have a model which explains their data. Great. I'm not sure this agreement can be used as an argument for the propriety of the 'thermodynamic approach.'

2. A simple plot showing mRNA/cell as a function of repressor concentration for the different configurations in Figure 6 would be much appreciated. Perhaps I missed the inference, but the ultimate goal of an input/output understanding is to map concentration of a regulator to expression level of the target gene.

Reviewer #2 (Recommendations for the authors):

The manuscript has several issues as noted below.

1) For MCMC fitting of experimental data, the absolute error or relative error of all the results should be given. One would like to know the confidence degree of the fitting results.

2) Formula (5) in the main text is used to fit the data, and the results are shown in figure6. However, the number of the parameters in the fitting formula is almost the same as the amount of the experimental data, and the formula is too complicated. The corresponding parameter space is quite large. Therefore, it is necessary to further analyze the accuracy of the fitting and the necessity of adding the cooperativity.

3) What is source of error bar of the initial RNAP loading rate? Need to give detailed explanation. If the source of error bar is from the data of multiple traces, why not fit multiple trajectories separately?

4) The fitting equations (1), (2) and (5), are all the fitting function for the experimental data. Are the parameter b, p and r the fitting parameters? If they are the fitting parameters, the fitting result needs to be given. These parameters are related to the weights of different binding states in your model, which need to be listed in detail, and compared the differences of fitting results under variant conditions.

5) In this paper, a variety of data is shown, which are Runt binding to different sites such as [001][010],[011] and so on. Please elaborate on how Runt is experimentally controlled or determined to bind to the specified sites.

6) The fluorescence of nucleus and cytoplasm is quantitatively analyzed, how to distinguish nucleus and cytoplasm? Please give the process and corresponding results.

7) It is shown in Fiugre6b that if the cooperativity is not included, a good fitting result cannot be obtained. Has the fitting covered a large enough parameter space search? The fitting results in figure 6 are quite different from the experiment results regardless of the cooperativity. If so, does it mean that the cooperativity plays a major role in the fitting, or even that the individual interaction between Runt and binding sites can be ignored? Could you separate out the part of cooperativity in the fitting results so as to confirm that the cooperativity plays a key role in the experiment?

Reviewer #3 (Recommendations for the authors):

1) Results of Figure 5 should be accompanied by reports of fitting alternative models and suitable model comparisons performed. Also, a baseline model with constant K_R and \omega_RP (same for all three constructs) should be evaluated and compared to. Such comparison should ideally involve model comparison techniques such as likelihood ratio and AIC/BIC.

2) Results of Figure 6 should be accompanied by statistical assessment of improvements due to higher-order cooperativity parameters.

3) The fact that fitting of more complex models typically relied on setting different values for the additional parameters for each construct (rather than the same values for all tested constructs of that category) should be examined more closely. How well can the data be modeled without this flexibility?

4) The authors should address the concern that alternative model structures different from that used in the thermodynamics-based models here might have explained the data without requiring higher-order cooperativity.

Other suggestions:

One of the first experimental tests reported is that of predicted expression profile in the hunchback promoter with one Runt binding site. Figure 3J is supposed to be a qualitative match to the predictions in Figure 3F. Could the authors elaborate on what reasonable regimes of the model parameters would provide predictions that do not qualitatively match the observation? It seems like with an activator (Bicoid) progressively decreasing and a repressor (Runt) progressively increasing from anterior to posterior, the enhancer expression will obviously decrease in some form. Perhaps the authors might wish to make it clear that their goal here was to establish some basic parameter estimates to be used in the remaining analyses, and perhaps also to assess the uncertainties associated with those model parameters.

The results presented in Figure 4 are puzzling. Figure 4D suggests that the differences among the constructs arises purely from differences in the p and R parameters, and not from the bicoid or bicoid-RNAP related parameters. Could the authors add a possible mechanistic speculation here regarding how making changes solely in the enhancer sequence might result in biochemical changes solely in the behavior of RNAP at the promoter (without any change in the protein-DNA interactions at the enhancer itself)?

[Editors' note: further revisions were suggested prior to acceptance, as described below.]

Thank you for resubmitting your work entitled "Predictive modeling reveals that higher-order cooperativity drives transcriptional repression in a synthetic developmental enhancer" for further consideration by eLife. Your revised article has been evaluated by Naama Barkai (Senior Editor) and a Reviewing Editor.

The manuscript has been significantly improved but please further address one issue on the validation of transcription rate measurement, the reviewers asked to use an orthogonal method such as smFISH to validate the measurement using MS2, which was not done in the revision. Please provide at least a thorough comparison and contrast of the current method with other orthogonal methods in the Discussion to equip readers with a full context of the quantifications that are central to the work.

https://doi.org/10.7554/eLife.73395.sa1

Author response

The essential revisions are:

1. Modeling: Please provide a better articulation/construction of the model regarding the fitting parameters (which to keep constant and which to vary, the number of parameters to use, error bars in the experimental data), whether an alternative model should be considered, whether the provided model can indeed account quantitatively observed data (not just to demonstrate that previous models do not work), and whether statistically the model fits experimental data adequately (not by overfitting, examine using Akaike Information Criterion AIC or Bayesian Information Criterion BIC). See reviewers #1 pt 1, #2 pt 1-5, and #3's overall review, and pt1.

We apologize for the lack of clarity in our original manuscript regarding the model construction as well as the statistical interpretation of the MCMC inference. In our revised manuscript, we colored the free parameters in red to emphasize which parameters are being inferred from a given model equation. We also included error bars representing 95% of confidence intervals for all of our MCMC inference results in the main Figures (Figure 4, 5, 6, and 7).

For the question of whether an alternative model should be considered, specifically for one-Runt binding site constructs, we tested alternative models under four different scenarios: (1) K_r and w_rp being constant across constructs, (2) K_r remains constant and w_rp varies, (3) K_r varies and w_rp remains constant, and (4) both K_r and w_rp vary across constructs. The results are given in Figure S16, where we can see that our original assumption of keeping K_r constant and varying w_rp indeed explains the data better than the other models. This result is further supported by the quantification of model performance using the Akaike Information Criterion (as suggested by the reviewers) shown in Figure S16F. Finally, the inferred parameters from the case where we let both K_r and w_rp vary across constructs are shown in Figure S16G. These results suggest that the K_r values do not vary significantly across the constructs whereas w_rp values vary significantly between the constructs. From the aforementioned evidence, we argue that our assumption of constant K_r across constructs is well supported.

To assuage the reviewers’ concern on whether the improved model with higher-order cooperativities considered in the main text indeed explains the data better than the others, we computed Akaike Information Criterion for all models presented in Figure 6, resulting in a new panel, Figure 6G. The result indeed supports that both Runt-Runt cooperativity and higher-order cooperativity (emphasis on the higher-order cooperativity) are required to explain the experimental data without overfitting.

Finally, in addition to the “direct repression” that was assumed throughout the main text, we tested alternative models based on different modes of repression that have been proposed in the literature (Gray et al., Genes and Development 8:1829, 1994,). Specifically, we tested two additional modes of repression, “competition” and “quenching” in Figure S7 and S8, respectively, and summarized the results in Section S5. In short, our goal was to see if the pairwise cooperativities are sufficient to explain the two-Runt binding sites cases in these alternate models of repression, or whether we had to invoke higher-order cooperativities as in the “direct repression” mechanism. We concluded that both modes of repression required higher-order cooperativities to explain the two-Runt binding sites cases (even though the competition mechanism generally gave worse fits than direct repression or quenching). Thus, we believe that our conclusion of the necessity of higher-order cooperativity to account for our experimental data holds regardless of the specific model of repression that is assumed.

2. Experiments: please provide experimental measurement for a few conditions where the outcome of transcription regulation by Runt is measured by smFISH or protein concentrations to demonstrate that the measured promoter activity using RNAP loading rate is equivalent or better (without complications in mRNA/protein degradation).

We have revised our manuscript to include the results for accumulated mRNA patterns. Here is the summary of our major revisions.

First, we quantified the accumulated mRNA level (pattern) by integrating the MS2 time traces during nuclear cycle 14. This method was previously validated to give quantitatively comparable results to single-molecule fluorescence in situ hybridization (smFISH) (Garcia et al., Current Biology, 23:2140, 2013). The patterns of accumulated mRNA corresponding to each construct with and without the Runt protein are given in Figure S17. The figure clearly shows that the effect of repression is present at the level of accumulated mRNA.

Second, to demonstrate that the measured promoter activity using RNAP loading rate is equivalent to measures of accumulated mRNA, we compared the initial rate of RNAP loading to the accumulated mRNA at each position along the anterior-posterior axis for all of our constructs. The result is now included in a new figure, Figure S19. Interestingly, there is substantial correlation between these two metrics, with an average Pearson correlation coefficient of 0.9. Thus, we argue that the initial rate of RNAP loading, the quantity that we focused on in this paper, is comparable to the accumulated mRNA, suggesting that higher-order cooperativity is necessary to predict patterns of accumulated mRNA.

Regardless, we would like to emphasize that, even though there was a good correlation between the rate of transcription and the accumulated mRNA as shown in Figure S19, the objective of our work was not to predict accumulated mRNA patterns. Instead, our goal was to uncover the molecular mechanisms by which multiple transcription factors work together and to quantitatively predict the consequences of these mechanisms. We argue that the action of transcription factors is more accurately captured by monitoring the rate of transcription, rather than the accumulated mRNA, which could be potentially confounded by the regulation of different aspects of the transcription cycle (initiation, elongation and termination), as well as degradation and diffusion. We now emphasized this point in the Discussion (Lines 417-422).

Note that the reviewer requested to make this comparison using orthogonal techniques such as smFISH. However, we hope that the reviewers will agree that the already established relationship between the integral of the MS2 signal and the accumulated mRNA makes invoking such orthogonal techniques (which lack the temporal resolution of the MS2 system) unnecessary.

Reviewer #1 (Recommendations for the authors):

1. Modeling: I would like to see a better articulation of what the alternative is to a 'thermodynamic model.' I understand the authors meaning to be equilibrium models arrived at through partition functions. But kinetic models based on detailed balance are also thermodynamic. I don't want to get twisted around in the language, but the over-emphasis on the ideological purity of thermodynamics was distracting and seems to refer to an ongoing dialog / conversation which is likely lost on the general reader.

We agree with the reviewer’s comment. Kinetic models that account for the transition rates between the different binding states the enhancer become thermodynamic models in the extreme case where the rates that characterize these transitions are much faster than the rate of mRNA production. In this limit, we can invoke quasi-equilibrium, making it possible to forego rates altogether and focus on the probabilities of each state. As the reviewer suggests, we could have entertained such a more general kinetic model. However, our goal was to ask whether the simple thermodynamic model, which has been very successful in bacteria, could also predict transcriptional regulation in the animal context. We have now included a discussion of alternative models in the main text (Lines 107-108) and Discussion (Lines 492-493) to clarify our points.

a. One example is when they assume that the Runt dissociation constant remains unchanged with different positions (K_r) because the sequence is the same. The problem stems from their previous section where they found that there was an unsuspected relationship between the position of the Runt Binding sites and the parameter "p=[RNAP]/K_p" --- the point being that the sequence is the same in each construct, and therefore, based on their model K_p is changing. So why is it preferrable to assume that K_r is constant with position, but K_p changes?

We agree that the justification for the assumption that the Runt dissociation constant K_r is equivalent throughout our constructs should be clarified further.

In Section 2.3 and Figure 4, we attempted to determine whether the simple thermodynamic model derived in Equation 3 could explain our data, and if so, which parameters were consistent across our synthetic enhancer constructs, and which parameters varied. As the reviewer pointed out, the parameter “p=[RNAP]/K_p” changes between constructs, thus suggesting that K_p changes even though the sequence of the promoter is exactly the same across constructs. We do not know why this is the case, but we could speculate on a potential scenario.

We speculate that, since the enhancer is proximal to the promoter, the overall enhancer sequence might affect the binding of the transcriptional machinery itself. Therefore, the transcriptional machinery might see slightly different DNA sequences, a combination of the promoter and part of synthetic enhancer for each construct. Indeed, Louder et al. (Nature 531:604, 2018), showed using cryo-EM that promoter-bound general transcription factor IID (TFIID) occupies roughly 60 bp of DNA (from -20 to +40, where 0 represents the start of the transcription site in the promoter) when it forms pre-initiation complex. Thus, it is conceivable that the general transcriptional machinery would interact with the general DNA sequence of an enhancer located proximally upstream from the promoter. In contrast, the binding site sequence for Runt is small enough (8 bp) compared to the enhancer sequence (~300 bp) that Runt protein might see a similar local DNA sequence for each binding site. Thus, our results could be explained if the general transcriptional machinery interacts with much longer DNA sequences than Runt repressor. This discussion is now included in the Discussion (Lines 433-440)

Regardless of the reason behind the variation in K_p between reporter constructs, our assumption that the Runt dissociation constant remains equivalent across all binding sites is supported by relaxing this assumption and performing MCMC inference on an unconstrained model. In a new figure, Figure S16G, we show the inferred parameters corresponding to the case where we let both K_r and w_rp vary. In Author response image 1, it is clear that the inferred K_r values do overlap between constructs within error bars, whereas the p=[RNAP]/K_p varies across synthetic enhancer constructs. Thus, we believe that our assumption of K_r remaining constant across all of our constructs is valid, and chose to continue making this assumption throughout the main text. We have incorporated this evidence that the inferred K_r values do overlap between constructs within error bars in Figure S16.

Author response image 1

Download asset Open asset

Inferred values of K_r and [RNAP]/K_p across all synthetic enhancer constructs.

K_r is inferred from one-Runt binding site constructs, and [RNAP]/K_p is inferred from all constructs in the absence of Runt protein.

b. Expanding on the previous point, they then proceed to state that an unchanging K_r is supported by the fact that by only varying the w_rp they were able to fit their data. Considering the similar role of K_r with w_rp within the equation, I wonder if one could reach the same result varying the K_r and holding w_rp constant? Furthermore, who's to say that both don't vary between the constructs? Could the authors address this concern by varying K_r and further discussing this possibility and what it means for the interpretations of the model? In summary, the authors have a model which explains their data. Great. I'm not sure this agreement can be used as an argument for the propriety of the 'thermodynamic approach.'

This is a great suggestion. We tested alternative models under four different scenarios: (1) K_r and w_rp being constant across constructs, (2) K_r remains constant and w_rp varies, (3) K_r varies and w_rp remains constant, and (4) both K_r and w_rp vary across constructs. The results are given in Figure S16, where we can see that our original assumption of keeping K_r constant and varying w_rp indeed explains the data better than the other models. Briefly, the best MCMC fits to our data look qualitatively similar for the [001] and [010] constructs regardless of our assumptions, but differ for the [100] construct. For the [100] construct, the best MCMC fits predict that there should be a repression in the presence of Runt protein for the (1) and (3) scenarios. However, no repression is present in the data, suggesting that the assumption of keeping K_r constant is essential in explaining all three constructs, especially the [100] construct. This result is further supported by the quantification of model performance given by the Akaike Information Criterion shown in Figure S16F. The Akaike Information Criterion is a metric that reports the performance of the model. Briefly, Akaike Information Criterion quantifies the deviation between the model and the data, then penalizes the model by the number of free parameters. As a result, the smaller the Akaike Information Criterion is, the better the model performs. Our analysis suggests that the assumption of K_r remaining constant and w_rp varying is favorable. Finally, the inferred parameters from the case where we let both K_r and w_rp vary across constructs are shown in Figure S16G, showing that the inferred K_r values are comparable between constructs (within error bars), whereas w_rp values vary.

2. A simple plot showing mRNA/cell as a function of repressor concentration for the different configurations in Figure 6 would be much appreciated. Perhaps I missed the inference, but the ultimate goal of an input/output understanding is to map concentration of a regulator to expression level of the target gene.

In the new Figure S18, we now present the accumulated mRNA as a function of repressor concentration for all of our synthetic enhancers. We would like to emphasize that this visualization could be misleading since Bicoid and Runt both vary along the anterior posterior axis of the embryo. Thus, to avoid confusion regarding the simultaneous variation of Bicoid and Runt across the embryo, in the main text we favor our current visualization where we plot our transcriptional output along the embryo’s body axis (anterior-posterior axis).

Reviewer #2 (Recommendations for the authors):

The manuscript has several issues as noted below.

1) For MCMC fitting of experimental data, the absolute error or relative error of all the results should be given. One would like to know the confidence degree of the fitting results.

Thank you for pointing this out. We agree that we missed this important information in our first submission. We have now included the 95% confidence interval for all our MCMC fits and prediction plots in our revised submission. As a result, Figures 4, 5, 6, and 7 now contain proper confidence intervals from the MCMC inference.

2) Formula (5) in the main text is used to fit the data, and the results are shown in figure6. However, the number of the parameters in the fitting formula is almost the same as the amount of the experimental data, and the formula is too complicated. The corresponding parameter space is quite large. Therefore, it is necessary to further analyze the accuracy of the fitting and the necessity of adding the cooperativity.

We agree that what is now Equation 6 is quite complicated and that we did not clearly present which parameters are free and which parameters are fixed because they were determined using previous measurements of reporter constructs with simpler regulatory architectures. We revised the way we display the formula to emphasize that only w_rr and w_rrp are the free parameters by coloring these two parameters in red. Thus, the number of free parameters is indeed smaller than the number of data points.

3) What is source of error bar of the initial RNAP loading rate? Need to give detailed explanation. If the source of error bar is from the data of multiple traces, why not fit multiple trajectories separately?

We thank the reviewer for pointing out this lack of clarification in our manuscript. Figure 3I shows how the initial rate of RNAP loading is measured by extracting the slope resulting from a linear fit to the averaged MS2 time traces at the beginning of nuclear cycle 14. Here, the MS2 traces were averaged over a set of nuclei in an individual embryo that are in the same spatial bin along the anterior-posterior axis of the embryo (2.5% of the embryo length) as previously described (Garcia et al., Current Biology 23:2140, 2013; Eck and Liu et al., eLife 9:e56429, 2020). Then, for all figures throughout this paper, we averaged these initial RNAP loading rates over multiple embryos (>3 embryos) at each spatial position, with an error bar corresponding to the standard error of the mean. Thus, the source of the error bar in the initial RNAP loading rate stems from embryo-to-embryo variability. This information has now been included throughout the paper, especially in Figure 3 caption.

To address the reviewer’s concern of whether our parameter inference would change if we inferred the model parameters from individual embryos, instead of the mean initial RNAP loading rates averaged over multiple embryos, we now included a new figure, Figure S20. In Figure S20B, we show that the initial RNAP loading rate obtained from individual embryos and then averaged over these embryos does not differ significantly from the initial RNAP loading rate calculated from an averaged embryo. Additionally, for a construct with the largest embryo-to-embryo variability, [010], we compared the results of performing MCMC inference using these two different approaches to calculating the initial RNAP loading rate. In Figure S20C, we now show that the values of parameters inferred using both approaches are comparable.

4) The fitting equations (1), (2) and (5), are all the fitting function for the experimental data. Are the parameter b, p and r the fitting parameters? If they are the fitting parameters, the fitting result needs to be given. These parameters are related to the weights of different binding states in your model, which need to be listed in detail, and compared the differences of fitting results under variant conditions.

We apologize for the lack of clarity in the presentation of these equations in our first draft. The parameters b, p, and r are b=[Bicoid]/K_b, p=[RNAP]/K_p, and r=[Runt]/K_r, respectively. Thus, the fitting parameters are K_b, p, and K_r. However, K_b and p were already fitted in Figure 4 in the absence of Runt protein. These inferred parameters were taken as given in order to determine K_r in Figure 5. Thus, parameters are not inferred simultaneously, but sequentially. To clarify further, we now mark all free parameters in our equations in red in our manuscript to more easily distinguish free and fixed parameters.

5) In this paper, a variety of data is shown, which are Runt binding to different sites such as [001][010],[011] and so on. Please elaborate on how Runt is experimentally controlled or determined to bind to the specified sites.

We now elaborate on our synthetic enhancer construct design in Section 2.3 (Lines 219-226). Briefly, we used the same sequence motif that was identified in a study by Chen et al. (Cell 149:618, 2012). In that work, using synthetic enhancers, the binding motif used in our manuscript was validated to recruit Runt to repress the transcription of a lacZ reporter. Specifically, the resulting expression patterns were quantified using in situ hybridization in the presence and absence of Runt protein. We used the exact same sequence and positions of Runt binding sites across all our synthetic enhancer constructs, and we also observed transcriptional repression for all constructs by comparing the presence and absence of Runt protein. Additionally, this binding motif contains the canonical binding motifs for Run domain proteins, which have conserved amino acid sequences for the DNA binding domains between Drosophila Run and its human homologue Runx1 (Chen et al., Cell 149:618, 2012; Koromila et al., Cell Reports 28:855, 2019). Several studies using ChIP-chip, or electric mobility shift assay (EMSA) have identified this consensus sequence to be directly bound by Runx1 protein’s DNA binding domain (Melnikova et al., Journal of Virology 67:2408, 1993; Lewis et al., Journal of Virology 73:5535, 1999). We have now included these references supporting Runt’s direct binding to its binding sites where we introduced our synthetic enhancer constructs in the main text.

6) The fluorescence of nucleus and cytoplasm is quantitatively analyzed, how to distinguish nucleus and cytoplasm? Please give the process and corresponding results.

We segmented nuclei by using the signal stemming from a histone fusion to infrared fluorescent protein (His-iRFP). We have included a new figure (Figure S22) to clearly show this. We have now clarified this point in the Materials and methods.

7) It is shown in Figure 6b that if the cooperativity is not included, a good fitting result cannot be obtained. Has the fitting covered a large enough parameter space search? The fitting results in figure 6 are quite different from the experiment results regardless of the cooperativity. If so, does it mean that the cooperativity plays a major role in the fitting, or even that the individual interaction between Runt and binding sites can be ignored? Could you separate out the part of cooperativity in the fitting results so as to confirm that the cooperativity plays a key role in the experiment?

We agree with the reviewer that this point has to be clarified further. In brief, the curves in Figure 6B did not have any free parameters, as we used the values for the parameters K_r and w_rp that were inferred from Figure 5 for each Runt binding site (w_rr = 1, and w_rrp=1 in Equation 6 in this particular case).

To address the reviewer’s point about whether the higher-order cooperativity is sufficient to get a good fit without the interaction between Runt and its binding site, we performed a sensitivity analysis for the initial RNAP loading rate as a function of the Runt dissociation constant for its binding site, K_r. Briefly, we fixed the K_r value to the inferred value from one-Runt binding site constructs shown in Figure 5, varied the K_r values by 100-fold (0.1x or 10x K_r) and then repeated the MCMC inference for different scenarios of cooperativities to see if the higher-order cooperativity alone could give equivalently good fits to the data, regardless of the magnitude of K_r. The result is shown in our new figure, Figure S21. Figure S21B-E, center column shows the case of a 10 times larger value of K_r, where we found that the best MCMC fit cannot recapitulate our data even in the presence of both Runt-Runt cooperativity and higher-order cooperativity. On the other hand, Figure S21B-E, right column shows the case of K_r/10, the best MCMC fit could recapitulate the data quite well if we invoked higher-order cooperativity, but ignored Runt-Runt pairwise cooperativity. In conclusion, we found that the Runt dissociation constant, K_r, plays a pivotal role in fitting as well as the combination of different cooperativities, but that, even as we vary this dissociation constant, higher-order cooperativity needs to be invoked to obtain a quantitative agreement between theory and experiment.

Reviewer #3 (Recommendations for the authors):

1) Results of Figure 5 should be accompanied by reports of fitting alternative models and suitable model comparisons performed. Also, a baseline model with constant K_R and \omega_RP (same for all three constructs) should be evaluated and compared to. Such comparison should ideally involve model comparison techniques such as likelihood ratio and AIC/BIC.

We have included a new figure, Figure S16 to address this point about alternative model assumptions in the context of the one-Runt binding site constructs. We tested alternative models with four different scenarios: (1) K_r and w_rp being constant across constructs, (2) K_r remains constant and w_rp varies, (3) K_r varies and w_rp remains constant, and (4) both K_r and w_rp vary across constructs. The results are given in Figure S16, where we show that our original assumption of keeping K_r constant and varying w_rp indeed explains the data better than the other models with the minimal set of free parameters. These results are further supported by a new calculation of the Akaike Information Criterion shown in Figure S16F.

2) Results of Figure 6 should be accompanied by statistical assessment of improvements due to higher order cooperativity parameters.

We again thank the reviewer for bringing up this point. We now include an Akaike Information Criterion (AIC) assessment for all models presented in Figure 6 (Figure 6G) to show how the model improves as we change our assumptions.

3) The fact that fitting of more complex models typically relied on setting different values for the additional parameters for each construct (rather than the same values for all tested constructs of that category) should be examined more closely. How well can the data be modeled without this flexibility?

Indeed, the simplest, baseline model would have been to consider equal values of additional parameters across the synthetic enhancer constructs such as when we introduce K_r and w_rp in Figure 5 to explain the repression by Runt protein for the one-Runt binding site cases. We considered the reviewer’s suggestion and performed MCMC inference keeping K_r and w_rp constant. However, the result was quite off from the experimental data as shown in Figure S16B.

4) The authors should address the concern that alternative model structures different from that used in the thermodynamics-based models here might have explained the data without requiring higher order cooperativity.

As the reviewer pointed out, our objective was to test whether thermodynamic models could predict Runt repressor action for a wide array of regulatory architectures. Within the thermodynamic modeling framework, we have tested different modes of repression, such as “competition” and “quenching'' as shown in Figure S7 and S8, and shown that they all require higher-order cooperativities in order to quantitatively recapitulate the data.

Of course, it is interesting to consider what would happen if we expanded our models to, for example, the non-equilibrium regime. While we are very interested in this approach, we believe that testing the non-equilibrium models is beyond the scope of this manuscript, and now address this point explicitly in the Discussion (Lines 491-499).

Other suggestions:

One of the first experimental tests reported is that of predicted expression profile in the hunchback promoter with one Runt binding site. Figure 3J is supposed to be a qualitative match to the predictions in Figure 3F. Could the authors elaborate on what reasonable regimes of the model parameters would provide predictions that do not qualitatively match the observation? It seems like with an activator (Bicoid) progressively decreasing and a repressor (Runt) progressively increasing from anterior to posterior, the enhancer expression will obviously decrease in some form. Perhaps the authors might wish to make it clear that their goal here was to establish some basic parameter estimates to be used in the remaining analyses, and perhaps also to assess the uncertainties associated with those model parameters.

Figure 3J was given as an illustrative example of one of our measured profiles of the initial RNAP loading. Our objective with this figure was not to go into the detail of the reasonable parameter regimes that would be consistent with these data. Indeed, as the reviewer points out, the role of Figure 3J is to confirm the qualitative prediction that, given the Bicoid and Runt profiles, the initial rate of RNAP loading should decrease towards the posterior of the embryo. Figure 5 then attempts to quantitatively confront this type of data with our models to explore the reasonable regime of the model parameters as well as uncertainties. We clarified this in the main text, Section 2.2 (Lines 189-190) and Figure 3 caption.

The results presented in Figure 4 are puzzling. Figure 4D suggests that the differences among the constructs arises purely from differences in the p and R parameters, and not from the bicoid or bicoid-RNAP related parameters. Could the authors add a possible mechanistic speculation here regarding how making changes solely in the enhancer sequence might result in biochemical changes solely in the behavior of RNAP at the promoter (without any change in the protein-DNA interactions at the enhancer itself)?

Indeed, this result was puzzling to us. We think that our result reveals an important potential caveat in the study of synthetic enhancers: the basal level of expression should be measured experimentally in the absence of the transcription factor of interest. This can, of course, be experimentally challenging in cases where there are more than two transcription factors of interest.

As pointed out by the reviewer, the MCMC inference results shown in Figure 4 lead us to conclude that the different basal activities of our synthetic enhancer constructs in the absence of Runt protein is due to the variability in p=[RNAP]/K_p and the maximum rate of transcription R. Currently, we do not have a satisfactory explanation for this observation. We would like to leave this for future investigations using experimental methods that can directly address these possibilities, but want to offer some speculations which we repeat from our response to Reviewer #1.

We speculate that, since the enhancer is proximal to the promoter, the overall enhancer sequence might affect the binding of the transcriptional machinery itself. Therefore, the transcriptional machinery might see slightly different DNA sequences, a combination of the promoter and part of synthetic enhancer for each construct. Indeed, Louder et al. (Nature 531:604, 2018), used cryo-EM to show that promoter-bound general transcription factor IID (TFIID) occupies roughly 60 bp of DNA (from -20 to +40, where 0 represents the start of the transcription site in the promoter) when it forms pre-initiation complex. Thus, it is conceivable that the general transcriptional machinery would interact with the general DNA sequence of an enhancer located proximally upstream from the promoter. In contrast, the binding site sequence for Runt is small enough (8 bp) compared to the enhancer sequence (~300 bp) that Runt protein might see a similar local DNA sequence for each binding site. Thus, our results could be explained if the general transcriptional machinery interacts with much longer DNA sequences than Runt repressor. This discussion is now included in the Discussion (Lines 433-440)

[Editors' note: further revisions were suggested prior to acceptance, as described below.]

The manuscript has been significantly improved but please further address one issue on the validation of transcription rate measurement-- the reviewers asked to use an orthogonal method such as smFISH to validate the measurement using MS2, which was not done in the revision. Please provide at least a thorough comparison and contrast of the current method with other orthogonal methods in the Discussion to equip readers with a full context of the quantifications that are central to the work.

Once again, we thank the editors and reviewers for their insightful comments and suggestions on our manuscript. We further address the issue raised by the editors on the validation of transcription rate measurement via comparison to orthogonal methods, such as smFISH below.

First, we quantified the accumulated mRNA level (pattern) by integrating the MS2 time traces during nuclear cycle 14. This method was previously validated to give quantitatively comparable results to single-molecule fluorescence in situ hybridization (smFISH) (Garcia et al., Current Biology 23:2140, 2013; Lammers et al., PNAS 17:836, 2020). Briefly, Garcia et al. compared the patterns of accumulated mRNA along the anterior-posterior axis of the embryo at late nuclear cycle 14 as quantified via MS2 and smFISH. These two methods yielded quantitatively comparable results as shown in Figure S3H from Garcia et al., Current Biology, 23:2140, 2013. We have elaborated further on this validation study in the Discussion (Lines 419-429).

Furthermore, we have compared two publicly available datasets from the same reporter construct driven by hunchback P2 enhancer by fluorescence in situ hybridization (FISH; Park et al., eLife 8: e41266, 2019) and MS2 ( Eck and Liu et al., eLife 9:e56429, 2020). Again, we integrated the MS2 time traces during nuclear cycle 13 for the datasets from Eck and Liu et al., and then compared that to the patterns of accumulated mRNA quantified via fluorescence in situ hybridization at early nuclear cycle 14 from Park et al. (shown in Author response image 2). By normalizing the accumulated mRNA profiles, we could see a good match along the anterior-posterior axis. From these comparisons to orthogonal methods widely used in the field (both smFISH and FISH), we claim that the MS2 technique is an accurate method for quantifying the patterns of accumulated mRNA as much as smFISH or FISH.

Note that the reviewer requested we repeat these experiments in the context of our work. However, we hope that the reviewers will agree that the already established relationship between the integral of the MS2 signal and the accumulated mRNA makes repeating these measurements unnecessary.

Author response image 2

Download asset Open asset

Comparison of accumulated mRNA profiles measured by FISH and MS2 for a reporter construct driven by the *hunchback* P2 enhancer.

Normalized profiles of accumulated mRNA averaged over embryos from FISH (blue, Park et al., *eLife* 8:e41266, 2019) and MS2 (red, Eck and Liu et al., *eLife* 9:e56429). The blue and red curves show the mean and standard error of the mean over 16 and 9 embryos, respectively. The accumulated mRNA profiles were normalized to the average level within 25%-35% of the embryo length.

https://doi.org/10.7554/eLife.73395.sa2

Article and author information

Author details

Yang Joon Kim

Chan Zuckerberg Biohub, San Francisco, United States

Contribution
Conceptualization, Resources, Data curation, Software, Formal analysis, Investigation, Visualization, Methodology, Writing – original draft, Writing – review and editing

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0003-1742-5657
Kaitlin Rhee

Department of Chemical Biology, University of California, Berkeley, Berkeley, United States

Contribution
Data curation, Formal analysis, Investigation

Competing interests
No competing interests declared
Jonathan Liu

Department of Physics, University of California, Berkeley, Berkeley, United States

Contribution
Software, Investigation, Writing – review and editing

Competing interests
No competing interests declared
Selene Jeammet

Department of Biology, Ecole Polytechnique, Paris, France

Contribution
Investigation

Competing interests
No competing interests declared
Meghan A Turner

Biophysics Graduate Group, University of California, Berkeley, Berkeley, United States

Contribution
Resources

Competing interests
No competing interests declared
Stephen J Small

Department of Biology, New York University, New York, United States

Contribution
Conceptualization, Resources, Supervision

Competing interests
No competing interests declared
Hernan G Garcia
1. Chan Zuckerberg Biohub, San Francisco, United States
2. Department of Physics, University of California, Berkeley, Berkeley, United States
3. Biophysics Graduate Group, University of California, Berkeley, Berkeley, United States
4. Department of Molecular and Cell Biology, University of California, Berkeley, Berkeley, United States
5. Institute for Quantitative Biosciences–QB3, University of California at Berkeley, Berkeley, United States
Contribution
Conceptualization, Supervision, Funding acquisition, Investigation, Visualization, Writing – original draft, Project administration, Writing – review and editing

For correspondence
hggarcia@berkeley.edu

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-5212-3649

Funding

Burroughs Wellcome Fund (Career Award)

Hernan G Garcia

Sloan Research Foundation

Hernan G Garcia

Human Frontier Science Program

Hernan G Garcia

Searle Scholars Program

Hernan G Garcia

Shurl and Kay Curci Foundation

Hernan G Garcia

Hellman Foundation

Hernan G Garcia

National Institute of Health (DP2 OD024541-01)

Hernan G Garcia

National Science Foundation (1652236)

Hernan G Garcia

Korea Foundation for Advanced Studies (Graduate Student Fellowship)

Yang Joon Kim

Chan Zuckerberg Biohub (Chan Zuckerberg Biohub Investigator Award)

Hernan G Garcia

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We are grateful to Armando Reimer, Brandon Schlomann, Elizabeth Eck, Gabriella Martini, Jacques Bothma, Jeehae Park, Jia Ling, Matthew Norstad, Michael Eisen, Nicholas Lammers, Nipam Patel, Rob Phillips, Simon Alamos, Xavier Darzacq, and Yasemin Kirişçioğlu for their guidance and comments on our manuscript. This work was supported by the Burroughs Wellcome Fund Career Award at the Scientific Interface, the Sloan Research Foundation, the Human Frontiers Science Program, the Searle Scholars Program, the Shurl and Kay Curci Foundation, the Hellman Foundation, the NIH Director’s New Innovator Award (DP2 OD024541-01), NSF CAREER Award (1652236), and an NIH R01 Award (R01GM139913) to HGG, and a KFAS scholarship to YJK. HGG is also a Chan Zuckerberg Biohub investigator.

Senior Editor

Naama Barkai, Weizmann Institute of Science, Israel

Reviewing Editor

Jie Xiao, Johns Hopkins University, United States

Version history

Preprint posted: July 29, 2021 (view preprint)
Received: August 27, 2021
Accepted: December 9, 2022
Accepted Manuscript published: December 12, 2022 (version 1)
Version of Record published: January 12, 2023 (version 2)

Copyright

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.