A mechanism for hunchback promoters to readout morphogenetic positional information in less than a minute

Cell fate decisions in the fly embryo are rapid: hunchback genes decide in minutes whether nuclei follow the anterior/posterior developmental blueprint by reading out positional information in the Bicoid morphogen. This developmental system is a prototype of regulatory decision processes that combine speed and accuracy. Traditional arguments based on fixed-time sampling of Bicoid concentration indicate that an accurate readout is impossible within the experimental times. This raises the general issue of how speed-accuracy tradeoffs are achieved. Here, we compare fixed-time to on-the-fly decisions, based on comparing the likelihoods of anterior/posterior locations. We found that these more efficient schemes complete reliable cell fate decisions within the short embryological timescales. We discuss the influence of promoter architectures on decision times and error rates, present concrete examples that rapidly readout the morphogen, and predictions for new experiments. Lastly, we suggest a simple mechanism for RNA production and degradation that approximates the log-likelihood function.


Introduction
From development to chemotaxis and immune response, living organisms make precise decisions based on limited information cues and intrinsically noisy molecular processes, such as the readout of ligand concentrations by specialized genes or receptors (Houchmandzadeh et al., 2002;Perry et al., 2012;Takeda et al., 2012;Marcelletti and Katz, 1992;Bowsher and Swain, 2014). Selective pressure in biological decision-making is often strong, for reasons that range from predator evasion to growth maximization or fast immune clearance. In development, early embryogenesis of insects and amphibians unfolds outside of the mother, which arguably imposes selective pressure for speed to limit the risks of predation and infection by parasitoids (O'Farrell, 2015). In Drosophila embryos, the first 13 cycles of DNA replication and mitosis occur without cytokinesis, resulting in a multinucleated syncytium containing about 6000 nuclei (O'Farrell et al., 2004). Speed is witnessed both by the rapid and synchronous cleavage divisions observed over the cycles, and the successive fast decisions on the choice of differentiation blueprints, which are made in less than 3 min (Lucas et al., 2018).
In the early fly embryo, the map of the future body structures is set by the segmentation gene hierarchy (Nüsslein-Volhard et al., 1984;Houchmandzadeh et al., 2002;Jaeger, 2011). The definition of the positional map starts by the emergence of two (anterior and posterior) regions of distinct hunchback (hb) expression, which are driven by the readout of the maternal Bicoid (Bcd) morphogen gradient (Houchmandzadeh et al., 2002, Figure 1a). hunchback spatial profiles are sharp and the variance in hunchback expression of nuclei at similar positions along the AP axis is small (Desponds et al., 2016;Lucas et al., 2018). Taken together, these observations imply that the short-time readout is accurate and has a low error. Accuracy ensures spatial resolution and the correct positioning of future organs and body structures, while low errors ensure reproducibility and homogeneity among spatially close nuclei. The amount of positional information available at the transcriptional locus is close to the minimal amount necessary to achieve the required precision (Gregor et al., 2007b;Porcher et al., 2010;Garcia et al., 2013;Petkova et al., 2019). Furthermore, part of the morphogenetic information is not accessible for reading by downstream mechanisms (Tikhonov et al., 2015), as information is channeled and lost through subsequent cascades of gene activity. In spite of that, by the end of nuclear cycle 14 the positional information encoded in the gap gene readouts is sufficient to correctly predict the position of each nucleus within 2% of the egg length (Petkova et al., 2019). Adding to the time constraints, mitosis resets the binding of Anterior-posterior axis 2 % egg length The cartoon shows a projection on one plane of the embryo at nuclear cycles 10-13, when nuclei (red dots) have migrated to the surface of the embryo (O'Farrell, 2015). The activity of the hunchback gene decreases along the Anterior-Posterior (AP) axis. The green dots represent active transcription loci. The average concentration L(x) of maternal Bicoid decreases exponentially along the AP axis by about a factor five from the anterior (left) to the posterior (right) ends. Between the blue lines lies the boundary region. Its width d x is 2% of the egg length. hunchback activity decreases along the AP axis and undergoes a sharp drop around the boundary region. The half-maximum Bcd concentration position in WT embryos is shifted by about 5% of the egg-length toward the anterior with respect to the mid-embryo position (Struhl et al., 1992). We consider this hunchback half maximum expression as a reference point when describing the AP axis of the embryo. The hunchback readout defines the cell fate decision whether each nucleus will follow an anterior or the posterior gene expression program. (b) A typical promoter structure contains six binding sites for Bicoid molecules present at concentration L. k ¼ 6 indicates that the gene is active only when all binding sites are occupied, defining an all-or-nothing promoter architecture.
Other forms and details of the promoter structure will be discussed in Figure 2 and Figure 3. (c) The average number of nuclei making a mistake in the decision process as a function of the egg length position at cycle 11. For a fixed-time decision process completed within T ¼ 270 seconds (and k ¼ 6, that is, all-or-nothing activation scheme), a large number of nuclei make the wrong decision (full blue bars). T ¼ 270 seconds is the duration of the interphase of nuclear cycle 11 . For nuclei located in the boundary region either answer is correct so that we leave these bars unfilled. Most errors happen close to the boundary, as intuitively expected. See Appendix 1 for a detailed description of how the error is computed for the fixed-time decision strategy. (d) The time needed to reach the standard error probability of 32% (Gregor et al., 2007a) for the same process as in panel c (see also the subsection 'How many nuclei make a mistake?') as a function of egg length position. Decisions are easy away from the center but the time required for an accurate decision soars close to the boundary up to 50 min -much longer than the embryological times. Parameters for panels (c) and (d) are six binding sites that bind and unbind Bicoid without cooperaivity and a diffusion limited on rate per site max L ¼ 0:124s À1 , and an unbinding rate per site n 1 ¼ 0:0154s À1 that lead to half activation in the boundary region.
transcription factors (TF) during the phase of synchronous divisions (Lucas et al., 2018), suggesting that the decision about the nuclei's position is made by using information accessible within one nuclear cycle. Experiments additionally show that during the nuclear cycles 10-13 the positional information encoded by the Bicoid gradient is read out by hunchback promoters precisely and within 3 min (Lucas et al., 2018). Effective speed-accuracy tradeoffs are not specific to developmental processes, but are shared by a large number of sensing processes (Rinberg et al., 2006;Heitz and Schall, 2012;Chittka et al., 2009). This generality has triggered interest in quantitative limits and mechanisms for accuracy. Berg and Purcell derived the seminal tradeoff between integration time and readout accuracy for a receptor evaluating the concentration of a ligand (Berg and Purcell, 1977) based on its average binding occupancy. Later studies showed that this limit takes more complex forms when rebinding events of detached ligands (Bialek and Setayeshgar, 2005;Kaizu et al., 2014) or spatial gradients (Endres and Wingreen, 2008) are accounted for. The accuracy of the averaging method in Berg and Purcell, 1977 can be improved by computing the maximum likelihood estimate of the  Figure 2. The relation between promoter structure and on-the-fly decision-making. (a). Using six Bicoid binding sites, the promoter decides between two hypothetical Bcd concentrations L ¼ L 1 and L ¼ L 2 , given the actual (unknown) concentration L in the nucleus. The number of occupied Bicoid sites fluctuates with time (b) and we assume the gene is expressed (c) when the number of occupied Bicoid binding sites on the promoter is ! k (green dashed line, here k ¼ 2). The gene activity defines a non-Markovian telegraph process. The ratio of the likelihoods that the time trace of this telegraph process is generated by L ¼ L 1 vs L ¼ L 2 is the log-likelihood ratio used for decision-making (d). The log-likelihood ratio undergoes random excursions until it reaches one of the two decision boundaries (K, ÀK). In d. the actual concentration is L ¼ L 1 and the log-likelihood ratio hits the upper barrier and makes the right decision. When L ¼ L 2 , less Bicoid-binding sites are occupied (e) and the gene is less likely to be expressed (f), resulting in a negative drift in the log-likelihood ratio, which directs the random walk to the lower boundary ÀK and the L ¼ L 2 decision (g). We consider that all six binding sites bind Bicoid independently and are identical with binding rate per site max L ¼ 0:07s À1 and unbinding rate per site n 1 ¼ 0:08s À1 , e ¼ 0:2, time series of receptor occupancy for a given model (Endres and Wingreen, 2009;Mora and Wingreen, 2010). However, none of these approaches can result in a precise anterior-posterior (AP) decision for the hunchback promoter in the short time of the early nuclear cycles, which has led to the conclusion that there is not enough time to apply the fixed-time Berg-Purcell strategy with the desired accuracy (Gregor et al., 2007a). Additional mechanisms to increase precision (including internuclear diffusion) do yield a speed-up (Erdmann et al., 2009;Aquino et al., 2016), yet they are not sufficient to meet the 3-min challenge. The issue of the embryological speed-accuracy tradeoff remains open. The boxes indicate the states in which the gene is expressed for the 2-or-more activation rule (red box and red in panels b-d) where the gene is active when 2-or-more TF are bound and the 4or-more activation rule (blue box and blue in panels b-d) where the gene is active when 4-or-more TF are bound. (b) The dynamics of TF binding translates into bursting and inactive periods of gene activity. The OFF and ON time distributions are different for the two hypothetical concentrations (blue boxes for k ¼ 4 and red boxes for k ¼ 2). The Kullback-Leibler divergence between the distributions for the two hypothetical concentrations (D KL ) sets the decision time and is related to the difference in the area below the two distributions. For the k ¼ 4 activation rule, the OFF time distributions are similar for the two hypothetical concentrations but the ON times distributions are very different. The ON times are more informative for the k ¼ 4 activation rule than the k ¼ 2 activation rule (c) The drift V of the log-likelihood ratio characterizes the deterministic bias in the decision process. The differences in (b) translate into larger drift for k ¼ 4 for the same binding/unbinding dynamics. (d) The distribution of decision times (calculated as the first-passage time of the log-likelihood random walk) decays exponentially for long times. Higher drift leads to on average faster decisions than for the k ¼ 4 activation rule (mean decision times are shown in dashed lines). For all panels the six binding sites are independent and identical with L ¼ L 1 ¼ 5:88 m À3 , L 2 ¼ 5:32 m À3 , e ¼ 0:1, max L ¼ 0:14 s À1 and n ¼ 0:08 s À1 for all binding sites.
The approaches described above are fixed-time strategies of decision-making, that is, they assume that decisions are made at a pre-defined deterministic time T that is set long enough to achieve the desired level of error and accuracy. As a matter of fact, fixing the decision time is not optimal because the amount of available information depends on the specific statistical realization of the noisy signal that is being sensed. The time of decision should vary accordingly and therefore depend on the realization. This intuition was formalized by A. Wald by his Sequential Probability Ratio Test (SPRT) (Wald, 1945a). SPRT achieves optimality in the sense that it ensures the fastest decision strategy for a given level of expected error. The adaption of the method to biological sensing posits that the cell discriminates between two hypothetical concentrations by accumulating information through binding events, and by computing on the fly the ratio of the likelihoods (or appropriate proxies) for the two concentrations to be discriminated (Siggia and Vergassola, 2013). When the ratio 'strongly' favors one of the two hypotheses, a decision is triggered. The strength of the evidence required for decision-making depends on the desired level of error. For a given level of precision, the average decision time can be substantially shorter than for traditional averaging algorithms (Siggia and Vergassola, 2013). SPRT has also been proposed as an efficient model of decision-making in the fields of social interactions and neuroscience (Gold and Shadlen, 2007;Marshall et al., 2009) and its connections with non-equilibrium thermodynamics are discussed in Roldán et al., 2015.
Wald's approach is particularly appealing for biological concentration readouts since many of them, including the anterior-posterior decision faced by the hunchback promoter, appear to be binary decisions. Our first goal here is to specifically consider the paradigmatic example of the hunchback promoter and elucidate the degree of speed-up that can be achieved by decisions on the fly. Second, we investigate specific implementations of the decision strategy in the form of possible hunchback promoter architectures. We specifically ask how cooperative TF binding affects the sensing limits. Our results have implications beyond fly development and are generally relevant to regulatory processes. We identify promoter architectures that, by approximating Wald's strategy, do satisfy several key experimental constraints and reach the experimentally observed level of accuracy of hunchback expression within the (apparently) very stringent time limits of the early nuclear cycles.

Methodological setup
The decision process of the anterior vs posterior hunchback expression The problem faced by nuclei in their decision of anterior vs posterior developmental fate is sketched in Figure 1a. By decision we mean that nuclei commit to a cell fate through a process that is mainly irreversible leading to one of two classes of cell states that correspond to either the anterior or the posterior regions of the embryo, based on positional information acquired through gene activity. We limit our investigation of promoter architectures to the six experimentally identified Bicoid-binding sites (Figure 1b). We do not consider the known Hunchback-binding sites because before nuclear cycle 13 there is little time to produce sufficient concentrations of zygotic proteins for a significant feedback effect and the measured maternal hunchback profile has not been shown to alter anterior-posterior decision-making. Following the observation that Bcd readout is the leading factor in nuclei fate determination (Ochoa-Espinosa et al., 2009), we also neglect the role of other maternal gradients, for example Caudal, Zelda or Capicua (Jiménez et al., 2000;Sokolowski et al., 2012;Tran et al., 2018;Lucas et al., 2018), since the readout of these morphogens can only contribute additional information and decrease the decision time. We focus on the proximal promoter since no active enhancers have been identified for the hunchback locus in nuclear cycles 11-13 (Perry et al., 2011). Our results can be generalized to enhancers (Hannon et al., 2017), the addition of which only further improves the speed-accuracy efficacy, as we explicitly show for a simple model of Bicoid activated enhancers in the section 'Joint dynamics of Bicoid enhancer and promoter'. Since our goal is to show that accurate decisions can be made rapidly, we focus on the worst case decision-making scenario: the positional information (Wolpert et al., 2015) is gathered through a readout of the Bicoid concentration only, and the decision is assumed to be made independently in each nucleus. Having additional information available and/or coupling among nuclei can only strengthen our conclusion.
The profile of the average concentration of the maternal morphogen Bicoid LðxÞ is well represented by an exponential function that decreases from the anterior toward the posterior of the embryo : LðxÞ ¼ L 0 e À5ðxÀx0Þ=100 , where x is the position along the anterior-posterior axis measured in terms of percentage egg-length (EL), and x 0 is the position of half maximum hb expression corresponding to L 0 Bcd concentration. The decay length l ¼ 5 corresponds to 20% EL (Gregor et al., 2007a). Nuclei convert the graded Bicoid gradient into a sharp border of hunchback expression (Figure 1a), with high and low expressions of the hunchback promoter at the left and the right of the border respectively (Driever and Nüsslein-Volhard, 1988;Struhl et al., 1992;Crauk and Dostatni, 2005;Gregor et al., 2007b;Houchmandzadeh et al., 2002;Porcher et al., 2010;Garcia et al., 2013;Lucas et al., 2013;Tran et al., 2018;Lucas et al., 2018). We define the border region of width d x symmetrically around x 0 by the dashed lines in Figure 1a. d x is related to the positional resolution (Erdmann et al., 2009;Tran et al., 2018) of the anterior-posterior decision: it is the minimal distance measured in percentages of egg-length between two nuclei's positions, at which the nuclei can distinguish the Bcd concentrations. Although this value is not known exactly, a lower bound is estimated as d x~2 % EL (Gregor et al., 2007a), which corresponds to the width of one nucleus.
We denote the Bcd concentration at the anterior (respectively posterior) boundary of the border region by L 1 (respectively L 2 ) (see Figure 1a). At each position x, nuclei compare the probability that the local concentration LðxÞ is greater than L 1 (anterior) or smaller than L 2 (posterior). By using current best estimates of the parameters (see Appendix 1), a classic fixed-time-decision integration process and an integration time of 270 s (the duration of the interphase in nuclear cycle 11 ), we compute in Figure 1c the probability of error per nucleus for each position in the embryo (see Appendix 2 for details). As expected, errors occur overwhelmingly in the vicinity of the border region, where the decision is the hardest (Figure 1c). For nuclei located within the border region, both anterior and posterior decisions are correct since the nuclei lie close to both regions. It follows that, although the error rate can formally be computed in this region, the errors do not describe positional resolution mistakes and do not contribute to the total error (white zone in Figure 1c).
In view of Figure 1c and to simplify further analysis we shall focus on the boundaries of the border region : each nucleus discriminates between hypothesis 1 -the Bcd concentration is L ¼ L 1 , and hypothesis 2 -the Bcd concentration is L ¼ L 2 . To achieve a positional resolution of d x ¼ 2% EL, nuclei need to be able to discriminate between differences in Bcd concentrations on the order of 10%. In addition to the variation in Bcd concentration estimates that are due to biological precision, concentration estimated using many trials follows a statistical distribution. The central limit theorem suggests that this distribution is approximately Gaussian. This assumption means that the probability that the Bcd concentration estimate deviates from the actual concentration by more than the prescribed 10% positional resolution is 32% (see the subsection 'How many nuclei make a mistake?' for variations on the value and arguments on the error rate). In Figure 1d, we show that the time required under a fixed-time-decision strategy for a promoter with six binding sites to estimate the Bcd concentration within the 32% Gaussian error rate (Gregor et al., 2007a) close to the boundary is much longer than 270 s, ' 40 minutes (see Appendix 2 for details of the calculation). The activation rule for the promoter architecture in the figure is that all binding sites need to be bound for transcription initiation.

Identifying fast decision promoter architectures The promoter model
We model the hb promoter as six Bcd binding sites (Schrö der et al., 1988;Driever et al., 1989;Struhl et al., 1989;Ochoa-Espinosa et al., 2005) that determine the activity of the gene (Figure 2a). Bcd molecules bind to and unbind from each of the i ¼ 1; :::; 6 sites with rates m i and n i , which are allowed to be different for each site. For simplicity, the gene can only take two states : either it is silenced (OFF) and mRNA is not produced, or the gene is activated (ON) and mRNA is produced at full speed. While models that involve different levels of polymerase loading are biologically relevant and interesting, the simplified model allows us to gain more intuition and follows the worst-case scenario logic that we discussed in the previous subsection 'The decision process of the anterior vs posterior hunchback expression'. The same remark applies for the wide variety of promoter architectures considered in previous works (Estrada et al., 2016;Tran et al., 2018). In particular, we assume that only the number of sites that are bound matters for gene activation (and not the specific identity of the sites). Such architectures are again a subset of the range of architectures considered in Estrada et al., 2016;Tran et al., 2018. The dynamics of our model is a Markov chain with seven states with probability P i corresponding to the number of sites occupied: from all sites unoccupied (probability P 0 ) to all six sites bound by Bcd molecules (probability P 6 ). The minimum number k of bound Bicoid sites required to activate the gene divides this chain into the two disjoint and complementary subsets of active states (P ON ¼ P 6 i¼k P i , for which the gene is activated) and inactive states (P OFF ¼ P kÀ1 i¼0 P i , for which the gene is silenced) as illustrated in Figure 2b and d.
As Bicoid ligands bind and unbind the promoter (Figure 2b), the gene is successively activated and silenced (Figure 2c). This binding/unbinding dynamics results in a series of OFF and ON activation times that constitute all the information about the Bcd concentration available to downstream processes to make a decision. We note that the idea of translating the statistics of binding-unbinding times into a decision remains the same as in the Berg-Purcell approach, where only the activation times are translated into a decision (and not the deactivation times). The promoter architecture determines the relationship between Bcd concentration and the statistics of the ON-OFF activation time series, which makes it a key feature of the positional information decision process. Following (Siggia and Vergassola, 2013), we model the decision process as a Sequential Probability Ratio Test (SPRT) based on the time series of gene activation. At each point in time, SPRT computes the likelihood P of the observed time series under both hypotheses PðL 1 Þ and PðL 2 Þ and takes their ratio : RðtÞ ¼ PðL 1 Þ=PðL 2 Þ. The logarithm of RðtÞ undergoes stochastic changes until it reaches one of the two decision threshold boundaries K or ÀK (symmetric boundaries are used here for simplicity) ( Figure 2d). The decision threshold boundaries K are set by the error rate e for making the wrong decision between the hypothetical concentrations : K ¼ log ð1 À eÞ=e ð Þ (see Siggia and Vergassola, 2013 and Appendix 3). The choice of K or e depends on the level of reproducibility desired for the decision process. We set K ' 0:75, corresponding to the widely used error rate e ' 0:32 of being more than one standard deviation away from the mean of the estimate for the concentration assumed to be unbiased in a Gaussian model (see the subsection 'The decision process of the anterior vs posterior hunchback expression'). The statistics of the fluctuations in likelihood space are controlled by the values of the Bcd concentrations: when Bcd concentration is low, small numbers of Bicoid ligands bind to the promoter ( Figure 2e) and the hb gene spends little time in the active expression state (Figure 2f), which leads to a negative drift in the process and favors the lower one of the two possible concentrations (Figure 2g).

Mean decision time: connecting drift-diffusion and Wald's approaches
In this section, we develop new methods to determine the statistics of gene switches between the OFF and ON expression states. Namely, by relating Wald's approach (Wald, 1945a) with drift-diffusion, we establish the equality between the drift and diffusion coefficients in decision making space for difficult decision problems, that is, when the discrimination is hard, we elucidate the reason underlying the equality. That allows us to effectively determine long-term properties of the likelihood log-ratio and compute mean decision times for complex architectures.
A gene architecture consists of N binding sites and is represented by N + 1 Markov states corresponding to the number of bound TF, and the rates at which they bind or unbind (Figure 3a). For a given architecture, the dynamics of binding/unbinding events and the rules for activation define the two probability distributions P OFF ðt; LÞ and P ON ðs; LÞ for the duration of the OFF and ON times, respectively ( Figure 3b). The two series are denoted ft i g 1 i J þ and fs j g 1 j J À , where J þ and J À are the number of switching events in time t from OFF to ON and vice versa. For those cases where the two concentrations L 1 and L 2 are close and the discrimination problem is difficult (which is the case of the Drosophila embryo), an accurate decision requires sampling over many activation and deactivation events to achieve discrimination. The logarithm of the ratio RðtÞ can then be approximated by where V is the constant drift, that is, the bias in favor of one of the two hypotheses, D is the diffusion constant and h a standard Gaussian white noise with zero mean and delta-correlated in time (Wald, 1945a;Siggia and Vergassola, 2013). The decision time for the case of symmetric boundaries K ¼ ÀK À ¼ K þ , is given by the mean first-passage time for this biased random walk in the log-likelihood space (Redner, 2001;Siggia and Vergassola, 2013): Note that in this approximation all the details of the promoter architecture are subsumed into the specific forms of the drift V and the diffusion D.
Drift. We assume for simplicity that the time series of OFF and ON times are independent variables (when this assumption is relaxed, see Appendix 8). This assumption is in particular always true when gene activation only depends on the number of bound Bicoid molecules. Under these assumptions, we can apply Wald's equality (Wald, 1945b;Durrett, 2010) to the log-likelihood ratio, log RðtÞ. Wald considered the sum of a random number M of independent and identically distributed (i.i.d.) variables. The equality that he derived states that if M is independent of the outcome of variables with higher indices ðX i Þ i>M (i.e. M is a stopping time), then the average of the sum is the product hMi hX i i.
Wald's equality applies to our likelihood sum ( P M i log R i of the likelihoods, where M is the number of ON and OFF times before a given (large) time t). We conclude the drift of the log-likelihood ratio, log RðtÞ, is inversely proportional to ðt ON þ t OFF Þ, where t ON is the mean of the distribution of ON times P ON ðt; LÞ and t OFF is the mean of the distribution of OFF times P OFF ðs; LÞ. The term ðt ON þ t OFF Þ determines the average speed at which the system completes an activation/deactivation cycle, while the average hlog R i i describes how much deterministic bias the system acquires on average per activation/deactivation cycle. The latter can be re-expressed in terms of the Kullback-Leibler divergence D KL ðf jjgÞ ¼ R ¥ 0 dt 0 f ðt 0 Þ log f ðt 0 Þ=gðt 0 Þ ð Þ between the distributions of the OFF and ON times calculated for the actual concentration L and each one of the two hypotheses, L 1 and L 2 : Equation 2 quantifies the intuition that the drift favors the hypothetical concentration with the time distribution which is the closest to that of the real concentration L (Figure 3b).
Diffusivity : Why it is more involved to calculate and how we circumvent it. While the drift has the closed simple form in Equation 2, the diffusion term is not immediately expressed as an integral. The qualitative reason is as follows. Computing the likelihood of the two hypotheses requires computing a sum where the addends are stochastic (ratios of likelihoods) and the number of terms is also stochastic (the number of switching events). These two random variables are correlated: if the number of switching events is large, then the times are short and the likelihood is probably higher for large concentrations. While the drift is linear in the above sum (so that the average of the sum can be treated as shown above), the diffusivity depends on the square of the sum. The diffusivity involves then the correlation of times and ratios (Carballo-Pacheco et al., 2019), which is harder to obtain as it depends a priori on the details of the binding site model (see the subsection 'Equality between drift and diffusivity' and the subsection 'When are correlations between the times of events leading to decision important?' of Appendix 3 for details).
We circumvent the calculation of the diffusivity by noting that the same methods used to derive Equation (1) also yield the probability of first absorption at one of the two boundaries, say +K (see the subsection 'Equality between drift and diffusivity' of Appendix 3): By imposing P K ¼ 1 À e, we obtain VK=D ¼ log ð1 À eÞ=e ð Þ and the comparison with the expression of K ¼ log ð1 À eÞ=e ð Þ leads to the equality D ¼ V. The above equality is expected to hold for difficult decisions only. Indeed, drift-diffusion is based on the continuity of the log-likelihood process and Wald's arguments assume the absence of substantial jumps in the log-likelihood over a cycle. In other words, the two approaches overlap if the hypotheses to be discriminated are close. For very distinct hypotheses (easy discrimination problems), the two approaches may differ from the actual discrete process of decision and among themselves. We expect then that V ¼ D holds only for hypotheses that are close enough, which is verified by explicit examples (see the subsection 'Equality between drift and diffusivity' of Appendix 3). The Appendix subsection also verifies V ¼ D by expanding the general expression of V and D for close hypotheses. The origin of the equality is discussed below.
Using V ¼ D, we can reduce the general formula Equation (1) to which is formula 4.8 in Wald, 1945a and it is the expression that we shall be using (unless stated otherwise) in the remainder of the paper. The additional consequence of the equality V ¼ D is that the argument VK=2D of the hyperbolic . It follows that for any problem where the error e ( 1, the argument of the hyperbolic tangent is large and the decision time is weakly dependent on deviations to V ¼ D that occur when the two hypotheses differ substantially. A concrete illustration is provided in the subsection 'The first passage time to decision' of Appendix 3. Single binding site example. As an example of the above equations, we consider the simplest possible architecture with a single binding site (N ¼ 1), where the gene activation and de-activation processes are Markovian. In this case, the de-activation rate n is independent of TF concentration and the activation rate is exponentially distributed P OFF ðL; tÞ ¼ kLe kLt . We can explicitly calculate the drift V ¼ nkL=ðn þ kLÞðlogðL 1 =L 2 Þ þ ðL 2 À L 1 Þ=LÞ (Equation 2) and expand it for L 2 ¼ L and L 1 ¼ L þ dL, at leading order in dL. Inserting the resulting expression into Equation 4, we conclude that decreases with increasing relative TF concentration difference dL=L and gives a very good approximation of the complete formula (see Appendix 3-figure 2a-c with different values of dL=L). Equations 2 and 4 greatly reduce the complexity of evaluating the performance of architectures, especially when the number of binding sites is large. Alternatively, computing the correlation of times and log-likelihoods would be increasingly demanding as the size of the gene architecture transfer matrices increase. As an illustration, Figure 3 compares the performance of different activation strategies : the 2-or-more rule (k ¼ 2), which requires at least two Bcd-binding sites to be occupied for hb promoter activation (Figure 3a-d in blue), and the 4-or-more rule (k ¼ 4) (Figure 3a-d in red) for fixed binding and unbinding parameters. Figure 3c shows that stronger drifts lead to faster decisions. The full decision time probability distribution is computed from the explicit formula for its Laplace transform (Siggia and Vergassola, 2013, Figure 3d). With the rates chosen for Figure 3, the k ¼ 4 rule leads to an ON time distribution that varies strongly with the concentration, making it easier to discriminate between similar concentrations: it results in a stronger average drift that leads to a faster decision than k ¼ 2 ( Figure 3d).
What is the origin of the V ¼ D equality? The special feature of the SPRT random process is that it pertains to a log-likelihood. This is at the core of the V ¼ D equality that we found above. First, note that the equality is dimensionally correct because log-likelihoods have no physical dimensions so that both V and D have units of time À1 . Second, and more important, log-likelihoods are built by Bayesian updating, which constrains their possible variations. In particular, given the current likelihoods P 1 ðtÞ ¼ P j P ON ðt j ; L 1 ÞP OFF ðs j ; L 1 Þ and P 2 ðtÞ ¼ P j P ON ðt j ; L 2 ÞP OFF ðs j ; L 2 Þ at time t for the two concentrations L 1 and L 2 and the respective probabilities Q 1 ðtÞ ¼ P 1 =ðP 1 þ P 2 Þ and Q 2 ðtÞ ¼ 1 À Q 1 of the two hypotheses, it must be true that the expected values after a certain time Dt remain the same if the expectation is taken with respect to the current P i ðtÞ (see, e.g. Reddy et al., 2016). In formulae, this implies that the average variation of the probability Q 2 over a given time Dt that is should vanish (see the subsection 'Equality between drift and diffusivity' of Appendix 3 for a derivation). Here, hDQ 2 i 1 is the expected variation of Q 2 under the assumption that hypothesis 1 is true and hDQ 2 i 2 is the same quantity but under the assumption that hypothesis 2 is true. We notice now that Q 2 ðtÞ ¼ 1 1þe LðtÞ , where L ¼ log R is the log-likelihood, and that the drift-diffusion of the log-likelihood implies that hDLi 1 ¼ VDt, hDLi 2 ¼ ÀVDt and h DL À hDLi 1 ð Þ 2 i 1 ¼ h DL À hDLi 2 ð Þ 2 i 2 ¼ 2DDt. By using that dQ 2 =dL ¼ ÀQ 1 Q 2 and d 2 Q 2 =dL 2 ¼ ÀQ 1 Q 2 ðQ 2 À Q 1 Þ, we finally obtain that and imposing hDQ 2 i ¼ 0 yields the equality V ¼ D. Note that the above derivation holds only for close hypotheses, otherwise the velocity and the diffusivity under the two hypotheses do not coincide.

Additional embryological constraints on promoter architectures
In addition to the requirements imposed by their performance in the decision process (green dashed line in Figure 4a), promoter architectures are constrained by experimental observations and properties that limit the space of viable promoter candidates for the fly embryo. A discussion about their possible function and their relation to downstream decoding processes is deferred to the final section.
First, we require that the average probability for a nucleus to be active in the boundary region is equal to 0.5, as it is experimentally observed (Lucas et al., 2013;Figure 1a). This requirement mainly impacts and constrains the ratio between binding rates i and unbinding rates n i .
Second, there is no experimental evidence for active search mechanisms of Bicoid molecules for its targets. It follows that, even in the best case scenario of a Bcd ligand in the vicinity of the promoter always binding to the target, the binding rate is equal to the diffusion limited arrival rate max L ' 0:124s À1 (Appendix 1). As a result, the binding rates i are limited by diffusion arrivals and the number of available binding sites: i max ð7 À iÞ (black dashed line in Figure 4b), where L is the concentration of Bicoid. This sets the timescale for binding events. In Appendix 1, we explore the different measured values and estimates of parameters defining the diffusion limit max L and their influence on the decision time (see Appendix 1-table 1 for all the predictions).
Third, as shown in Figure 1, the hunchback response is sharp, as quantified by fitting a Hill function to the expression level vs position along the egg length. Specifically, the hunchback expression (in arbitrary units) f hb is well approximated as a function of the Bicoid concentration LðxÞ by the Hill function f hb ðxÞ ' LðxÞ H =ðLðxÞ H þ L H 0 Þ, where L 0 is the Bcd concentration at the half-maximum hb expression point and H is the Hill coefficient ( Figure 1a). Experimentally, the measured Hill coefficient for mRNA expression from the WT hb promoter is H~7 À 8 (Lucas et al., 2018;Tran et al., 2018). Recent work  suggests that these high values might not be achieved by Bicoid-binding sites only. Given current parameter estimates and an equilibrium binding model,  shows that a Hill coefficient of 7 is not achievable within the duration of an early nuclear cycle (' 5 min). That points at the contribution of other mechanisms to pattern steepness. Given these reasons (and the fact that we limit ourselves only to a model with six equilibrium Bcdbinding sites only), we shall explore the space of possible equilibrium promoter architectures limiting the steepness of our profiles to Hill coefficients H~4 À 5.

Numerical procedure for identifying fast decision-making architectures
Using Equations 2 and 4, we explore possible hb promoter architectures and activation rules to find the ones that minimize the time required for an accurate decision, given the constraints listed in the paragraph 'Additional embryological constraints on promoter architectures'. We optimize over all possible binding rates ð i Þ 1 i 6 (m 1 is the binding rate of the first Bcd ligand and 6 the binding rate of the last Bcd ligand when 5 Bcd ligands are already bound to the promoter), and the unbinding rates ðn i Þ 1 i 6 (n 1 is the unbinding rate of a single Bcd ligand bound to the promoter and n 6 is the unbinding rate of all Bcd ligands when all six Bcd-binding sites are occupied). We also explore different activation rules by varying the minimal number of Bcd ligands k required for activation in the kor-more activation rule (Estrada et al., 2016;Tran et al., 2018). We use the most recent estimates of biological constants for the hb promoter and Bcd diffusion (see Appendix 1) and set the error rate at the border to 32% (Gregor et al., 2007b;Petkova et al., 2019). Reasons for this choice were given in the subsection 'The decision process of the anterior vs posterior hunchback expression' and will be revisited in the subsection 'How many nuclei make a mistake?', where we shall introduce some embryological considerations on the number of nuclei involved in the decision process and determine the error probability accordingly. The optimization procedure that minimized the Figure 4. Performance, constraints and statistics of fastest decision-making architectures. (a) Mean decision time for discriminating between two concentrations with jL 2 À L 1 j ¼ 0:1L and e ¼ 0:32. Results shown for the fastest decision-making architectures for different activation rules and steepness constraints. For a given activation rule (k), we optimize over all values of ON rates i and OFF rates n i (see Figure 3a) within the diffusion limit (0:124 s À1 per site), constraining the steepness H and probability of nuclei to be active at the boundary (see the paragraph 'Additional embryological constraints on promoter architectures'). The green lines denotes the interphase duration of nuclear cycle 11 and even for the strongest constraints (H>5) we identify architectures that make an accurate decision within this time limit. (b) The unbinding rates (blue) and binding rates (red) of the fastest decisionmaking architectures with H>5 -all these regulatory systems require cooperativity in TF binding to the promoter-binding sites. The dashed line on the ON rates plots shows the upper bound set by the diffusion limit. (c) Histogram of the probability distribution of the time spent in different Bcd-binding site occupancy states for the fastest decision-making architecture for k ¼ 3 and no constraints on the slope (blue), H>4 (red) and H>5 (black). (d) Probability distribution of the time spent bound to th DNA by Bicoid molecules for the fastest decision-making architecture with H>5 and k ¼ 2. Our prediction is compared to the exponential distribution with parameters fit by Mir et al., 2017, for the specific binding at the boundary. While the distributions are close, our simulated distribution is not exponential, as expected for the 6-binding site architecture. The non-exponential behavior in the experimental curve is likely masked by the convolution with non-specific binding. We use the boundary region concentration L ¼ 5:6 m À3 (see panel b, k ¼ 2 for rates). average decision time for different values of k and H is implemented using a mixed strategy of multiple random starting points and steepest gradient descent (Figure 4a).

Logic and properties of the identified fast decision architectures
The main conclusion we reach using the methodology presented in the 'Methodological setup' section is that there exist promoter architectures that reach the required precision within a few minutes and satisfy all the additional embryological constraints that were discussed previously ( Figure 4a). The fastest promoters (blue crosses in Figure 4a) reach a decision within the time of nuclear cycle 11 (green line in Figure 4a) for a wide range of activation rules. Even imposing steep readouts (H>4) allows us to identify relatively fast promoters, although imposing the nuclear cycle time limit, pushes the activation rule to smaller k. Interestingly, we find that the fastest architectures identified perform well over a range of high enough concentrations (Appendix 4-figure 1). The optimal architectures differ mainly by the distribution of their unbinding rates ( Figure 4b). We shall now discuss their properties, namely the binding times of Bicoid molecules to the DNA binding sites, and the dependence of the promoter activity on various features, such as activation rules and the number of binding sites in detail. Together, these results elucidate the logic underlying the process of fast decisionmaking.

How many nuclei make a mistake?
The precision of a stochastic readout process is defined by two parameters: the resolution of the readout dx, and the probability of error, which sets the reproducibility of the readout. In Figure 4, we have used the statistical Gaussian error level (32%) to obtain our results. However, the error level sets a crucial quantity for a developing organism and it is important to connect it with the embryological process, namely how many nuclei across the embryo will fail to properly decide (whether they are positioned in the anterior or in the posterior part of the embryo). To make this connection, we compute this number for a given average decision time t and we integrate the error probability along the AP axis to obtain the error per nucleus e s . The expected number of nuclei that fail to correctly identify their position is given by hn error i ¼ e s 2 cÀ1 , where c is the nuclear cycle and we have neglected the loss due to yolk nuclei remaining in the bulk and arresting their divisions after cycle 10 (Foe and Alberts, 1983). Assuming a 270 s readout time -the total interphase time of nuclear cycle 11 ) -for the fastest architecture identified above and an error rate of 32%, we find that hn error i ' 0:3, that is an essentially fail-proof mechanism. This number can be compared with >30 nuclei in the embryo that make an error in a hTi ¼ 270s read-out in a Berg-Purcell-like fixed-time scheme (integrated blue area in Figure 1c).
Conversely, for a given architecture, reducing the error level increases drastically the mean firstpassage time to decision: the mean time for decision as a function of the error rate for the fastest architecture identified with H>5 and k ¼ 1 is shown in Appendix 2-figure 1. The decision can be made in about a minute for e ¼ 32% but requires on average 10 min for e ¼ 10% (Appendix 2-figure 1). Note that, because the mean first-passage depends simply on the inverse of the drift per cycle (Equation 4), the relative performance of two architectures is the same for any error rate so that the fastest architectures identified in Figure 4a are valid for all error levels.
Just like for the fixed-time strategy (Figure 1c and d), nuclei located in the mid-embryo region are more likely to make mistakes and take longer on average to trigger a decision (Appendix 4-figure 2).

Residence times among the various states
As shown in Tran et al., 2018, high Hill coefficients in the hunchback response are associated with frequent visits of the extreme expression states where available binding sites are either all empty (state 0), or all occupied (state 6). Figure 4c provides a concrete illustration by showing the distribution of residence times for the promoter architectures that yield the fastest decision times for k ¼ 3 and no constraints (blue bars), H>4 (red bars) and H>5 (black bars). When there are no constraints on the slope of the hunchback response, the most frequently occupied states are close to the ON-OFF transition (2 and 3 occupied binding sites in Figure 4c) to allow for fast back and forth between the active and inactive states of the gene and thereby gather information more rapidly by reducing t ON þ t OFF (see formulae 2 and 4).
We notice that for higher Hill coefficients, the system transits quickly through the central states (in particular states with 3 and 4 occupied Bcd sites, Figure 4 red and black bars). As expected for high Hill coefficients, such dynamics requires high cooperativity. Cooperativity helps the recruitment of extra transcription factors once one or two of them are already bound and thus speeds up the transitioning through the states with 2, 3 and 4 occupied binding sites. An even higher level of cooperativity is required to make TF DNA binding more stable when 5 or 6 of them are bound, reducing the OFF rates n 5 or n 6 ( Figure 4b).
The (short) binding times of Bicoid on DNA The distribution of times spent bound to DNA of individual Bicoid molecules is shown in Figure 4d obtained from Monte Carlo simulations using rates from the fastest architecture with H>5 and k ¼ 2. We find an exponential decay, an average bound time of about 7.1 s and a median around 0.5 s. Our median-time-bound prediction is of the same order of magnitude as the observed bound times seen in recent experiments by Mir et al., 2017;Mir et al., 2018, who found short (mean~0.629 s and median~0.436 s based on exponential fits), yet quite distinguishable from background, bound times to DNA. These results were considered surprising because it seemed unclear how such short binding events could be consistent with the processing of ON and OFF gene switching events. Our results show that such short binding times may actually be instrumental in achieving the tradeoff between accuracy and speed, and rationalize how longer activation events are still achieved despite the fast binding and unbinding. High cooperativity architectures lead to non-exponential bound times to DNA (Figure 4d) for which the typical bound time (median) is short but the tail of the distribution includes slower dynamics that can explain longer activation events (the mean is much larger than the median). This result suggests that cells can use the bursty nature of promoter architectures to better discriminate between TF concentrations.
In Mir et al., 2017, the raw distribution comprises both non-specific and specific binding and cannot be directly compared to simulation results. Instead, we use the largest of the two exponents fit for the boundary region (Mir et al., 2017), which should correspond to specific binding. The agreement between the distributions in Figure 4d is overall good, and we ascribe discrepancies to the fact that (Mir et al., 2017) fit two exponential distributions assuming the observed times were the convolution of exponential specific and non-specific binding times. Yet the true specific binding time distribution is likely not exponential, e.g. due to the effect of binding sites having different binding affinities. We show in Appendix 5-figure 1 that the two distributions are very similar and hard to distinguish once they are mixed with the non-specific part of the distribution.

Activation rules
In the parameter range of the early fly embryo, the fastest decision-making architectures share the one-or-more (k ¼ 1) activation rule : the promoter switches rapidly between the ON and OFF expression states and the extra binding sites are used for increasing the size of the target rather than building a more complex signal. Architectures with k ¼ 2 and k ¼ 3 activation rules can make decisions in less than 270 s and satisfy all the required biological constraints. Generally speaking, our analysis predicts that fast decisions require a small number of Bicoid-binding sites (less than three) to be occupied for the gene to be active. The advantage of the k ¼ 2 or k ¼ 3 activation rules is that the ON and OFF times are on average longer than for k ¼ 1, which makes the downstream processing of the readout easier. We do not find any architecture satisfying all the conditions for the k ¼ 4; 5; 6 activation rules, although we cannot exclude there could be some architectures outside of the subset that we managed to sample, especially for the k ¼ 4 activation rules where we did identify some promoter structures that are close to the time constraint.
Activation rules with higher k can give higher information per cycle for the ON rate, yet they do not seem to lead to faster decisions because of the much longer duration of the cycles. To gain insight on how the tradeoff between fast cycles and information affects the efficiency of activation rules, we consider architectures with only two binding sites, which lend to analytical understanding (Figure 5a and b). Both of these considered architectures are out of equilibrium and require energy consumption (as opposed to the two equilibrium architectures of Figure 5c and d). We compare promoters of different complexity: the all-or-nothing k ¼ 2 out-of-equilibrium model (a), the 1-or-more k ¼ 1 out-of-equilibrium model (b), the two binding site all-or-nothing k ¼ 2 equilibrium model (c) and the one binding site equilibrium model (d). (e) Comparison of the mean decision time between k ¼ 2 (a) and k ¼ 1 (b) activation schemes for the two binding site out-of-equilibrium models as a function of the unbinding rate n and binding rate m 1 . The binding rate m 2 is fixed to 2 ¼ 1 =2. The fastest decision-making solution associates the second binding with slower variables to maximize V. Along the line L 1 ¼ n the activation rules k ¼ 1 and k ¼ 2 perform at the same speed. e ¼ 10%. (f) Comparison of the mean decision time for equilibrium architectures with one (d) and two (c) equilibrium TFbinding sites. 1 ; 2 ; n 1 are optimized at fixed n 2 , e ¼ 5%. We set a maximum value of 5 m 3 s À1 for m 1 and m 2 , corresponding to the diffusion limited arrival at the binding site. For n 1 , the maximum value of 5 s À1 corresponds to the inverse minimum time required to read the presence of a ligand, or to differentiate it from unspecific binding of other proteins. Additional binding sites are beneficial at high ligand concentrations and for small unbinding rates. In the blue region, the fastest mean decision time for a fixed accuracy assuming equilibrium binding, comes from a two binding site architecture with a non-zero first unbinding rate. In the white region, one of the binding rates !0 (see inset), which reduces to a one binding site model. e ¼ 5%. (g) Weaker binding sites can lead to faster decision times within a range of parameters (gray stripe). We consider the k ¼ 1 activation scheme with two binding sites (b). For fixed L (x axis), n (y axis) and 1 ¼ 0:2 m 3 s À1 , we optimize over m 2 while setting 2 < 1 (context of diffusion limited first and second bindings). The gray regions corresponds to parameters for which the optimal second unbinding rate Ã 2 < 1 and the second binding is weak. In the white region When is 1-or-more faster than all-or-nothing activation?
A first model has the promoter consisting of two binding sites with the all-or-nothing rule k ¼ 2 ( Figure 5a). We consider the mathematically simpler, although biologically more demanding, situation where TFs cannot unbind independently from the intermediate states -once one TF binds, all the binding sites need to be occupied before the promoter is freed by the unbinding with rate n of the entire complex of TFs. This situation can be formulated in terms of a non-equilibrium cycle model, depicted for two binding sites in Figure 5a. The activation time is a convolution of the exponential distributions P OFF ðt; LÞ ¼ 12L 2À1 e À1Lt À e À2Lt ð Þ. In the simple case, when the two binding rates are similar ( 1 ' 2 ), the OFF times follow a Gamma distribution and the drift and diffusion can be computed analytically (see Appendix 4). When the two binding rates are not similar the drift and diffusion must be obtained by numerical integration (see Appendix 4). In the first model described above (Figure 5a), deactivation times are independent of the concentration and do not contribute to the information gained per cycle and, as a result, to V=ðt ON þ t OFF Þ. To explore the effect of deactivation time statistics on decision times, we consider a cycle model where the gene is activated by the binding of the first TF (the 1-or-more k ¼ 1 rule) and deactivation occurs by complete unbinding of the TFs complex ( Figure 5b). The resulting activation times are exponentially distributed and contribute to drift and diffusion as in the simple two state promoter model (Figure 5d). The deactivation times are a convolution of the concentrationdependent second binding and the concentration-independent unbinding of the complex and their probability distribution is P ON ðt; LÞ ¼ n2L nÀ2L e À2Lt À e Ànt ð Þ. Drift and diffusion can be obtained analytically (Appendix 4). The concentration-dependent deactivation times prove informative for reducing the mean decision time at low TF concentrations but increase the decision time at high TF concentrations compared to the simplest irreversible binding model. In the limit of unbinding times of the complex (1=n) much larger than the second binding time (1= 2 L), no information is gained from deactivation times. In the limit of 1 L=n ! ¥, the k ¼ 1 model reduces to a one binding site exponential model and the two architectures (Figure 5b and d) have the same asymptotic performance.
Within the irreversible schemes of Figure 5a and Figure 5b, the average time of one activation/ deactivation cycle is the same for the all-or-nothing k ¼ 2 and 1-or-more k ¼ 1 activation schemes. The difference in the schemes comes from the information gained in the drift term V=ðt ON þ t OFF Þ, which begs the question : is it more efficient to deconvolve the second binding event from the first one within the all-or-nothing k ¼ 2 activation scheme, or from the deactivation event in the k ¼ 1 activation scheme?
In general, the convolution of two concentration-dependent events is less informative than two equivalent independent events, and more informative than a single binding event. For small concentrations L, activation events are much longer than deactivation events. In the k ¼ 1 scheme, OFF times are dominated by the concentration-dependent step 2 L and the two activation events can be read independently. This regime of parameters favors the k ¼ 1 rule (Figure 5e). However, when the concentration L is very large the two binding events happen very fast and for 2 L>>n, in the k ¼ 1 scheme, it is hard to disentangle the binding and the unbinding events. The information gained in the second binding event goes to 0 as L ! ¥ and the one-or-more k ¼ 1 activation scheme (Figure 5b) effectively becomes equivalent to a single binding site promoter (Figure 5d), making the all-or-nothing k ¼ 2 activation (Figure 5a) scheme more informative (Figure 5e). The fastest decision time architecture systematically convolves the second binding event with the slowest of the other reactions (Figure 5e), with the transition between the two activation schemes when the other reactions have exactly equal rates ( 1 L ¼ n line in Figure 5e) (see Appendix 6 for a derivation).

How the number of binding sites affects decisions
The above results have been obtained with six binding sites. Motivated by the possibility of building synthetic promoters (Park et al., 2019) or the existence of yet undiscovered binding sites, we investigate here the role of the number of binding sites. Our results suggest that the main effect of additional binding sites in the fly embryo is to increase the size of the target (and possibly to allow for higher cooperativity and Hill coefficients). To better understand the influence of the number of binding sites on performance at the diffusion limit, we compare a model with one binding site ( Figure 5d) to a reversible model with two binding sites where the gene is activated only when both binding sites are bound (all-or-nothing k ¼ 2, Figure 5c). Just like for the six binding site architectures, we describe this two binding site reversible model by using the transition matrix of the N þ 1 Markov chain and calculate the total activation time P OFF ðt; LÞ.
For fixed values of the real concentration L, the two hypothetical concentrations L 1 and L 2 , the error e and the second off-rate n 2 , we optimize the remaining parameters m 1 , m 2 and n 1 for the shortest average decision time.
For high gene deactivation rates n 2 , the fastest decision time is achieved by a promoter with one binding site (Figure 5f): once one ligand has bound, the promoter never goes back to being completely unbound (n Ã 1 =0 in Figure 5c) but toggles between one and two bound TF (Figure 5d with n ¼ n 2 and ¼ 2 ). For lower values of gene deactivation rates n 2 , there is a sharp transition to a minimal hTi solution using both binding sites. In the all-or-nothing activation scheme that is used here, the distribution of deactivation times is ligand independent and the concentration is measured only through the distribution of activation times, which is the convolution of the distributions of times spent in the 0 and 1 states before activation in the two state. For very small deactivation n 2 rates, it is more informative to 'measure' the ligand concentration by accumulating two binding events every time the gene has to go through the slow step of deactivating ( Figure 5c). However, for large deactivation rates little time is 'lost' in the uninformative expressing state and there is no need to try and deconvolve the binding events but rather use direct independent activation/deactivation statistics from a single binding site promoter (Figure 5d, see Appendix 7 for a more detailed calculation).

The role of weak binding sites
An important observation about the strength of the binding sites that emerge from our search is that the binding rates are often below the diffusion limit max L 0 ' 0:124s À1 (see black dashed line in Figure 4b) : some of the ligands reach the receptor, they could potentially bind but the decision time decreases if they do not. In other words, binding sites are 'weak' and, since this is also a feature of many experimental promoters (Gertz et al., 2009), the purpose of this section is to investigate the rationale for this observation by using the models described in Figure 5.
Naively, it would seem that increasing the binding rate can only increase the quality of the readout. This statement is only true in certain parameter regimes, and weaker binding sites can be advantageous for a fast and precise readout. To provide concrete examples, we fix the deactivation rate n and the first binding rate m 1 in the 1-or-more irreversible binding model of Figure 5b and we look for the unbinding rate Ã 2 that leads to the fastest decision. We consider a situation where the two binding sites are not interchangeable and binding must happen in a specific order. In this case, the diffusion limit states that 2 1 if the first binding is strong and happens at the diffusion limit. We optimize the mean decision time for 0 2 1 (see Appendix 9-figure 1 for an example) and find a range of parameters where the fastest-decision value Ã 2 < 1 is not as fast as parameter range allows (Figure 5g). We note that this weaker binding site that results in fast decision times can only exist within a promoter structure that features cooperativity. In this specific case, the first binding site needs to be occupied for the second one to be available. If the two binding sites are independent, then the diffusion limit is 2 1 and the fastest hTi solution always has the fastest possible binding rates.

Predictions for Bicoid-binding sites mutants
In addition to results for wild type embryos, our approach also yields predictions that could be tested experimentally by using synthetic hb promoters with a variable numbers of Bicoid-binding sites (Figure 6a). For any of the fast decision-making architectures identified and activation rules chosen, we can compute the effects of reducing the number of binding sites. Specifically, our predictions for the k ¼ 3 activation rule and H>4 in Figure 6b can be compared to FISH or fluorescent live imaging measurements of the fraction of active loci at a given position along the anterior-posterior axis. Bcd-binding site mutants of the WT promoter have been measured by immunostaining in cycle 14 (Park et al., 2019), although mRNA experiments in earlier cell cycles of well characterized mutants are needed to provide for a more quantitative comparison.
An important consideration for the comparison to experimental data is that there is a priori no reason for the hb promoter to have an optimal architecture. We do find indeed many architectures that satisfy all the experimental constraints and are not the fastest decision-making but 'good  Egg length boundary displacement in EL% Two or more rule Three or more rule enough' hb promoters. A relevant question then is whether or not similarity in performance is associated with similarity in the microscopic architecture. This point is addressed in Figure 6c, where we compare the fraction of active loci along the AP axis using several constraint-conforming architectures for H>4 and the k ¼ 2 and k ¼ 3 activation rules. The plot shows that solutions corresponding to the same activation rule are clustered together and quite distinguishable from the rest. This result suggests that the precise values of the binding and unbinding constants are not important for satisfying the constraints, that many solutions are possible, and that FISH or MS2 imaging experiments can be used to distinguish between different activation rules. The fraction of active loci in the boundary region is an even simpler variable that can differentiate between different activation rules ( Figure 6d). Lastly, we make a prediction for the displacement of the anterior-posterior boundary in mutants, showing that a reduced numbers of Bcd sites results in a strong anterior displacement of the hb expression border compared to six binding sites, regardless of the activation rule ( Figure 6e). Error bars in Figure 6e, that correspond to different close-to-fastest architectures, confirm that these share similar properties and different activation rules are distinguishable.

Joint dynamics of Bicoid enhancer and promoter
The Bicoid transcription factor has been shown to target more than a thousand enhancer loci in the Drosophila embryo with a wide concentration range of sensitivities (Driever and Nüsslein-Volhard, 1988;Struhl et al., 1989;Hannon et al., 2017). Enhancers are of special interest because they can be located far away from the promoter (Ribeiro et al., 2010;Krivega and Dean, 2012) and perform a statistically independent sample of the concentration that is later combined with that of the promoter. Evidence suggests that promoter-based conformational changes can be stable over long times (Fukaya et al., 2016), which mimics information storage during a process of signal integration.
To explore these effects, we consider a simple model of enhancer dynamics where a Bicoid-specific enhancer switches between two states ON and OFF independently of the promoter dynamics. We assume a simple rule for the gene activity: the gene is transcribed when both the promoter and the enhancer are ON. As an example, we consider an enhancer made of one binding site so that the ON rate e ON of the enhancer is limited by the diffusion rate max L. As an example, we perform a parameter search for the promoter activation rule k ¼ 2 (see Appendix 12), while still assuming that about half the nuclei are active at the boundary and a required Hill coefficient greater than 4, looking for architectures yielding the shortest decision time for an error rate of 32% for the 10% relative concentration difference discrimination problem. We find that the enhancer improves the performance of the readout, reducing the time to decision by about ' 6%. We find that adding extra binding sites to the promoter increases the computing power of the enhancer-promoter system and can reduce the time to decision to about 60% of the performance of the best architectures without enhancers.

Estimating the log-likelihood function with RNA concentrations
To illustrate how a biochemical network can approximate the calculation of the log-likelihood, we consider the case of the fastest architecture identified in the paragraph 'Numerical procedure for identifying fast decision-making architectures' with k ¼ 1 and H>4. In Figure 7a, we show the contributions of the OFF times (blue) and the ON times (red) to the log-likelihood for this architecture. We notice that the behavior of the log-likelihood contributions at long times is simply linear in time with a positive rate for ON times and a negative rate for OFF times. Conversely, short ON times contribute negatively to the log-likelihood while short OFF times contribute positively to the log-likelihood (Figure 7a). This observation suggests a simple model of RNA production with delay to approximate the computation of the log-likelihood. We consider a model of RNA production with five parameters (Figure 7b, details of the model are given in Appendix 11). We assume that when the promoter is in the ON state, polymerase is loaded and RNA is transcribed at a constant rate r ON while when the promoter is in the OFF state, RNA is produced at a lower basal rate r b . In order to approximate the shown result in the fastest decisions for H>4 and k ¼ 2 and k ¼ 3. Error bar width is the standard deviation of the various architectures that are close to minimal hTi. For all panels, LðxÞ has an exponentially decreasing profile with decay length one fifth of total egg length with L 0 ¼ 5:6 m À3 at the boundary. Parameters are given in Appendix 10.
linear behavior of the log-likelihood function at long times we suggest the existence of an enzyme actively degrading hunchback RNA (Chanfreau, 2017) and include in our model of RNA production the delay d ON and inertia d OFF associated with promoter dynamics and polymerase loading. RNA levels fluctuate and trigger decisions when the concentration ½RNA is greater than a threshold c 1 ðtÞ (anterior decision) or lower than a threshold c 2 ðtÞ (posterior decision). Since many forms of switch have already been presented in the literature (Goldbeter and Koshland, 1981;Tyson et al., 2003;Ozbudak et al., 2004;Siggia and Vergassola, 2013;Sandefur et al., 2013), we shall concentrate on the log-likelihood calculation and refer to previous references for the implementation of a decision when reaching a threshold. We look for parameters that satisfy both a high speed and high accuracy requirement for a decision between two points located 2% egg lengths apart across the mid-embryo boundary. For the fastest architecture identified with k ¼ 1 and H>4, we identify parameters that satisfy e<32% and a mean decision time T<3 min (see Appendix 11-figure 1). We check that this model produces a profile of RNA that is consistent with the observed high Hill coefficient (Figure 7c). Interestingly, we find that for this particular set of parameters the RNA profile Hill coefficient is increased by the delayed transcription dynamics and the active degradation from ' 4 (blue line in Figure 7c) up to ' 5:2 (red line in Figure 7c, details of the calculation are given in Appendix 11). This result could shed new light on the fundamental limits to Hill coefficients in the context of cooperative TF binding (Estrada et al., 2016;Tran et al., 2018) and provide a possible mechanism to explain how mRNA profiles can reach higher steepness than the corresponding TF activities do. We also looked for parameter sets that approximate the log-likelihood for the optimal architecture identified for k ¼ 2 and H>4 and find several candidates that fall close to the requirement of speed and accuracy (Appendix 11-figure 1). Together these results show that implementing the log-likelihood using a  Mir et al., 2018), RNA is being produced at rate r ON . When OFF, after a delay d OFF , the gene switches to basal rate r b . Hunchback RNA is being degraded actively by an enzyme at rate r OFF . The RNA is in excess for this reaction. (c) The model of RNA production with delay that yields an error of less than 32% in less than 3 min produces a profile of RNA production with high Hill coefficient ' 5:2 (red lines, renormalized RNA profile) that is higher than the Hill coefficient of renormalized gene activity H ' 4 (blue line). Parameters for promoter activity are those of the fastest architecture identified with k ¼ 1 and H>4, parameters for RNA production are d ON ¼ 0:047s, d OFF ¼ 1:6s, r b ¼ 0:2s -1 , r OFF ¼ 0:5s -1 , r ON ¼ 0:805s -1 . molecular circuit with a hb promoter is possible. They do not show this is what is happening in the embryo itself.

Discussion
The issue of precision in the Bicoid readout by the hunchback promoter has a long history (Nüsslein-Volhard and Wieschaus, 1980;Tautz, 1988). Recent interest was sparked by the argument that the amount of information at the hunchback locus available during one nuclear cycle is too small for the observed 2% EL distance between neighboring nuclei that are able to make reproducible distinct decisions (Gregor et al., 2007b). By using updated estimates of the biophysical parameters Tran et al., 2018), and the Berg-Purcell error estimation, we confirm that establishing a boundary with 2% variability between neighbouring nuclei would take at least about 13.4 min -roughly the non-transient expression time in nuclear cycle 14 (Lucas et al., 2018;Tran et al., 2018) (Appendix 1). This holds for a single Bicoid-binding site. An intuitive way to achieve a speed up is to increase the number of binding sites: multiple occupancy time traces are thereby made available, which provides a priori more information on the Bicoid concentration. Possible advantages of multiple sites are not so easy to exploit, though. First, the various sites are close and their respective bindings are correlated (Kaizu et al., 2014), so that their respective occupancy time traces are not independent. That reduces the gain in the amount of information. Second, if the activation of gene expression requires the joint binding of multiple sites, the transition to the active configuration takes time. The overall process may therefore be slowed down with respect to a single binding site model, in spite of the additional information. Third, and most importantly, information is conveyed downstream via the expression level of the gene, which is again a single time trace. This channeling of the multiple sites' occupancy traces into the single time trace of gene expression makes gene activation a real information bottleneck for concentration readout. All these factors can combine and even lead to an increase in the decision time. To wit, an all-or-nothing equilibrium activation model with six identical binding sites functioning at the diffusion limit and no cooperativity takes about 38 min to achieve the same above accuracy. In sum, the binding site kinetics and the gene activation rules are essential to harness the potential advantage of multiple binding sites.
Our work addresses the question of which multisite promoters architecture minimize the effects of the activation bottleneck. Specifically, we have shown that decision schemes based on continuous updating and variable decision times significantly improve speed while maintaining the desired high readout accuracy. This should be contrasted to previously considered fixed-time integration strategies. In the case of the hunchback promoter in the fly embryo, the continuous update schemes achieve the 2% EL positional resolution in less than 1 min, always outperforming fixed-time integration strategies for the same promoter architecture (see Appendix 1-table 1). While 1 min is even beyond what is required for the fly embryo, this margin in speed allows to accommodate additional constraints, viz. steep spatial boundary and biophysical constraints on kinetic parameters. Our approach ultimately yields many promoter architectures that are consistent with experimental observables in fly embryos, and results in decision times that are compatible with a precise readout even for the fast nuclear cycle 11 (Lucas et al., 2018;Tran et al., 2018).
Several arguments have been brought forward to suggest that the duration of a nuclear cycle is the limiting time period for the readout of Bicoid concentration gradient. The first one concerns the reset of gene activation and transcription factor binding during mitosis. In that sense, any information that was stored in the form of Bicoid already bound to the gene is lost. The second argument is that the hunchback response integrated over a single nuclear cycle is already extremely precise. However, none of these imply that the hunchback decision is made at a fixed-time (corresponding to mitosis) so that strategies involving variable decision times are quite legitimate and consistent with all the known phenomenology.
We have performed our calculations in a worst-case scenario. First, we did not consider averaging of the readout between neighbouring nuclei. While both protein (Gregor et al., 2007a) and mRNA concentrations (Little et al., 2013) are definitely averaged, and it has been shown theoretically that averaging can both increase and decrease (Erdmann et al., 2009) readout variability between nuclei, we do not take advantage of this option. The fact that we achieve less than 3 min in nuclear cycle 11, demonstrates that averaging is a priori dispensable. Second, we demand that the hunchback promoter results in a readout that gives the positional resolution observed in nuclear cycle 14, in the time that the hunchback expression profile is established in nuclear cycle 11. The reason for this choice is twofold. On the one hand, we meant to show that such a task is possible, making feasible also less constrained set-ups. On the other hand, the hunchback expression border established in nuclear cycle 11 does not move significantly in later nuclear cycles in the WT embryo, suggesting that the positional resolution in nuclear cycle 11 is already sufficient to reach the precision of later nuclear cycles. The positional resolution that can be observed in nuclear cycle 11 at the gene expression level is~10% EL , but this is also due to smaller nuclear density.
Two main factors generally affect the efficiency of decisions: how information is transmitted and how available information is decoded and exploited. Decoding depends on the representation of available information. Our calculations have considered the issue of how to convey information across the bottleneck of gene activation, under the constraint of a given Hill coefficient. The latter is our empirical way of taking into account the constraints imposed by the decoding process. High Hill coefficients are a very convenient way to package and represent positional information: decoding reduces to the detection of a sharp transition, an edge in the limit of very high coefficients. The interpretation of the Hill coefficient as a decoding constraint is consistent with our results that an increase in the coefficient slows down the decision time. The resulting picture is that promoter architecture results from a balance between the constraints imposed by a quick and accurate readout and those stemming from the ease of its decoding. The very possibility of a balance is allowed by the main conclusion demonstrated here that promoter structures can go significantly below the time limit imposed by the duration of the early nuclear cycles. That leaves room for accommodating other features without jeopardising the readout timescale. While the constraint of a fixed Hill coefficient is an effective way to take into account constraints on decoding, it will be of interest to explore in future work if and how one can go beyond this empirical approach. That will require developing a joint description for transmission and decoding via an explicit modeling of the mechanisms downstream of the activation bottleneck.
Recent work has shed light on the role of out of equilibrium architectures on steepness of response (Estrada et al., 2016) and gradient establishment Park et al., 2019).
Here, we showed that equilibrium architectures perform very well and achieve short decision times, and that out of equilibrium architectures do not seem to significantly improve the performance of promoters, except for making some switches from gene states a bit faster. Non-equilibrium effect can, however, increase the Hill coefficient of the response without adding extra binding sites, which is useful for the downstream readout of positional information that we formulated above as decoding.
We also showed how short bound times of Bicoid molecules to the DNA (Mir et al., 2017;Mir et al., 2018) are translated into accurate and fast decisions. Our fast decision-making architectures also display short DNA-bound times. However, the constraint of high cooperativity means that the distribution of bound times to the DNA is non-exponential and the rare long binding times that occur during the bursty binding process are exploited during the read-out. The combination of high cooperativity and high temporal variance due to bursty dynamics is a possible recipe for an accurate readout.
At the technical level, we developed new methods for the mean decision time of complex gene architectures within the framework of variable time decision-making (SPRT). This allowed us to establish the equality V ¼ D between drift and diffusion of the log-likelihood between two close hypotheses. Its underlying reason is the martingale property that the conditional expectation of probabilities for two hypotheses, given all prior history, is equal to their present value. The methodology developed here will be useful for the broad range of decision processes where SPRT is relevant, including neuroscience (Gold and Shadlen, 2007;Bitzer et al., 2014) and synthetic biology (O'Brien and Murugan, 2019;Pittayakanchit et al., 2018).
We made predictions about how promoter architectures with different activation schemes can be compared in synthetic embryos with different numbers of Bcd binding sites. Furthermore, experiments that change the composition of the syncytial medium would influence the diffusion constant and assay the assumption of diffusion-limited activation. Our model predicts that these changes would result in modifications of hunchback activation profiles: higher or lower diffusion rates slide the hunchback profile towards the anterior or the posterior end of the embryo, respectively, similarly to an increase or a decrease of the number of Bicoid binding sites. Any of the above experiments would greatly advance our understanding of the molecular control of spatial patterning in Drosophila embryo and, more generally, of regulatory processes.
Finally, we showed how a simple model of RNA production with delay and active degradation could easily approximate the seemingly complex log-likelihood calculation. The specific implementation that we described is tentative and aimed at simplicity, yet it illustrates several points, namely that the log-likelihood can be approximated by preserving a high Hill coefficients, observed characteristics of transcription dynamics and ensuring speed-accuracy limits. Future experiments could identify candidates for the enzyme responsible for the active degradation of RNA, or image the formation and dissolution of clusters using super-resolution imaging methods. This mechanism is also very close to the multifactor clusters (Mir et al., 2018) observed recently in early Drosophila in conjunction with active transcription sites. We suggest this mechanism as an example of implementation of the log-likelihood calculation, but note that the added complexity of the computation could happen at different levels of the expression machinery, including upstream of promoter activation through enhancer dynamics and chromatin arrangement.

Additional information
Competing interests Aleksandra M Walczak: Reviewing editor, eLife. The other authors declare that no competing interests exist. The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication. . Transparent reporting form

Data availability
All data analyzed during this study were previously published in the literature and references are included in the paper.

Biological parameters in the embryo
To build a model of the promoter, we combine parameters from recent experimental work. The embryo at nuclear cycle n is modelled as having 2 nÀ1 nuclei/cells. For simplicity, we assume they are equidistributed on the periphery of the embryo and across the embryo's length and do not take into account the effect of the geometry of the embryo because the curvature of the embryo is small (the embryo is about 500 mm-long along the anterior-posterior axis and only about 100 mmlong along the dorso-ventral axis). We also neglect the few nuclei forming pole cells and remaining in the bulk (Foe and Alberts, 1983).
The Bicoid concentration in the embryo is given by LðxÞ ¼ L 0 e ÀðxÀx0Þ=l , where x is the position along the anterior-posterior (AP) embryo axis measured in % of the egg length, l is the decay length also measured in % of the egg length (l ' 100m which is roughly 20% of the egg length (Gregor et al., 2007b) and x 0 is the position of half-maximum hunchback expression (x 0 is of the order of 250 mm and varies slightly depending on cell cycle, usually close to 45% egg length for the WT hunchback promoter . L 0 ' 5:6m À3 is the concentration of free Bicoid molecules at the AP boundary  (that also corresponds to the point of largest hunchback expression slope . To compute the diffusion limited arrival rate at the locus, we use the following parameters: diffusivity D ' 7:4m 2 s À1 Abu-Arish et al., 2010), concentration of free Bicoid molecules at the AP boundary L 0 , size of the binding target a ' 3nm (Gregor et al., 2007a), which leads to an effective max L 0 ¼ DaL 0 ' 0:124s À1 at the boundary. This value is an upper bound, assuming that every encounter between a transcription factor and a binding site results in successful binding. We note in Main Text Figure 4b that most of the ON rates are close to the diffusion limit. We conclude that in this parameter regime, the most efficient strategy is to have ON events that are as fast as possible. The only reason to reduce them is to achieve the required Hill coefficient. That can be done by either adjusting the ON or the OFF rates.
The above estimate max L 0 ¼ 0:124s À1 may be inaccurate for various reasons and we ought to explore the sensitivity of results to those uncertainties. Appendix 1-table 1 recapitulates the timeperformance of different strategies for different choices of the parameters. A first source of uncertainty is the value of the diffusivity, which is estimated to vary between 4 and 7m 2 s À1 Abu-Arish et al., 2010). We consider then two possible values for the diffusivity from Porcher et al., 2010: D 1 ¼ 4:5m 2 s À1 and the aforementioned value D 2 ¼ 7:4m 2 s À1 . A second source of uncertainty is that the actual Bicoid concentration at the boundary could vary by up to a factor two depending on estimates of the concentration and its decay length (Gregor et al., 2007b;Tran et al., 2018;Abu-Arish et al., 2010). We consider then two possible value for the concentration: the aforementioned L ð1Þ 0 ¼ 5:6 molecules per m 3 , and L ð2Þ 0 ¼ 11:2 molecules per m 3 . Finally, we assumed above that the size of the target is the full Bicoid operator site, which we took about ten base pairs following Gregor et al., 2007a. However, assuming that the TF must reach a specific position on the promoter-binding site could reduce the size of the target to a single base pair, that is a by a factor 10. In terms of parameters, we consider then two possible values for a : either a 1 ' 3:10 À4 m, or a 2 ¼ 3:10 À3 m. All in all, taking into account the various sources of uncertainty max L 0 can range in the interval ½0:0076; 0:25s À1 .
We consider four possible decision strategies. The first one is a single binding site making a fixed-time decision. This computation is made using the original Berg-Purcell formula (Berg and Purcell, 1977). In the Berg-Purcell strategy, the concentration of the ligand is inferred based on the total time that the receptor or the binding site has spent occupied by ligands. Due to averaging, the relative error of the concentration readout is inversely proportional to the number of independent measurements of the concentration that can be made within the total fixed time T, that is, to the arrival rate of new Bicoid molecules at the binding site, multiplied by the probability to find the binding site empty (in our case, at the boundary, the probability is roughly one half). Since the rate of arrival of new Bicoid molecules to the binding site is max L 0 ¼ DaL 0 , the relative error of the concentration readout is given by where L 0 is the estimate of the concentration, T is the integration time and n is the probability that the binding site is full.
In the second strategy, there are two sets of six binding sites being read independently. In that case, the information from each binding site can be accessed individually and their contributions averaged to give a more precise estimate. This calculation again can be made using the original Berg-Purcell formula (Berg and Purcell, 1977) for several receptors, dividing the relative error by the square root of the total number of binding sites in Equation 8: where N is the total number of binding sites (in our case, N ¼ 12).
For the third possibility, we consider a decision made at a fixed time using the fastest architecture identified (Main Text Figure 4) without constraint on the slope (activation rule k ¼ 1). We compute the decision time using the drift-diffusion approximation with fixed time (see Appendix 2).
Finally, we consider the fastest architecture identified with a random decision time and the Sequential Probability Ratio Test (SPRT) strategy.
The result of the above calculations is that for a single receptor estimating Bcd concentration with 10% precision within a fixed-time Berg-Purcell type calculation (see Appendix 2), decisions take between 6 min for the fastest binding rates and~4 hr for the slowest estimates. Conversely, by using the on-the-fly SPRT decision-making process and the one-or-more k ¼ 1 scheme at equilibrium, the time needed to make decisions with 10% precision and an error rate of 32% at the boundary is ' 30s for the fastest rates and ' 17 min for the slowest rates. For all sets of parameters, the on-the-fly SPRT decision-making process gives a~3.5-fold faster decision time than the N ¼ 12 Berg-Purcell estimate and a >10-fold faster decision making time than the one-binding-site Berg-Purcell estimate. For the fastest rates, a decision with an error rate of less than 5% can be achieved in about 5 min within the SPRT scheme.
Appendix 1-table 1. Mean decision times for various choices of parameters and four different decision processes (see sections 'Biological parameters in the embryo' and 'Error rate and decision time for the fixed-time decision strategy'). For the optimal architectures identified (third and fourth lines of the table) we take the fastest architectures without any constraints on the slope. These architectures systematically use the activation rule k ¼ 1. Highlighted in red are the results for the range of parameters presented in the text. All calculations are made with e ¼ 32%, L ¼ L 1 ¼ 1:05 Á L 0 , L 2 ¼ 0:95 Á L 0 . The diffusion limited ON rate is max L 0 ¼ DaL 0 . For the two Berg-Purcell architectures, both the ON rate and OFF rate per site are equal to max . For the optimal architectures, we have i L 0 ¼ max L 0 Á ð7 À iÞ and n i ¼ i Á max L 0 Á 0:5 1=6 =ð1 À 0:5 1=6 Þ to keep half the genes active at the boundary.

Error rate and decision time for the fixed-time decision strategy
In this section, we describe how we compute the decision time for a fixed-time strategy (or 'Berg-Purcell type decision') and a complex promoter architecture. The classic Berg-Purcell calculation is based on the idea of averaging the time spent by the ligand bound to a receptor (or, in our case, a binding site). The original calculation assumed that the waiting times between binding and unbinding are exponential and that the bound times are not informative about the concentration. Neither of these assumptions hold in the case of the hunchback promoter. Indeed, in the context of a complex promoter architecture, the waiting times that are available to the nucleus or cell downstream are the gene's ON and OFF switching times. They are not exponentially distributed, and, depending on the activation rule, the OFF times can be just as informative about the concentration as the ON times. For these two reasons, we ought to readapt the Berg-Purcell idea to compute the mean decision time.
To that purpose, we consider a decision with a given rate of error e, fix a time of decision T and choose the concentration that has the highest likelihood between the two options L 1 and L 2 . In other words, if log RðTÞ ¼ log PðL 1 Þ=PðL 2 Þ>0 then the nucleus chooses L ¼ L 1 , while if log RðTÞ<0 then it chooses L ¼ L 2 . For instance, if the actual concentration L ¼ L 1 , then the probability of error at time T is given by Pðlog RðtÞ<0Þ.
To calculate the above error, we use the drift-diffusion approximation for log RðtÞ, compute the drift V and diffusivity D from Main Text Equations 1-4 and approximate the distribution of log RðTÞ by a normal distribution with mean VT and standard deviation ffiffiffiffiffiffiffi DT p . We compute the error rate e for the fixed-time decision process based on the Gaussian approximation. Finally, to find the mean decision time for a given error rate, we perform a quick one-dimensional search over T, which yields the value of the fixed decision time appropriate for the prescribed error e.

The mean first-passage time for the decision-making process
To investigate the role of promoter architectures in decision-making, we apply the SPRT approach (Sequential Probability Ratio Test) (Wald, 1945a;Siggia and Vergassola, 2013). In the simplest formulation of this approach, a decision is made between two hypotheses about the concentration of a surrounding TF: either the TF is at concentration L 1 or it is at concentration L 2 . The decision is based on the binding and unbinding events of the TF to a promoter. At each point in time, the error of committing to a given concentration (scenario) is estimated by computing the ratio RðtÞ of the probability of one hypothesis, PðL 1 Þ (the surrounding concentration is L 1 ) over the other (the concentration is L 2 ), PðL 2 Þ: Technically, the time dependent likelihood ratio, RðtÞ undergoes a random walk as TFs bind and unbind to the promoter, activating and deactivating the gene. The process is terminated and a decision is made when the likelihood ratio falls below the desired error level, which is expressed in terms of absorbing boundaries at K + (the concentration is L 1 ) and K À (the concentration is L 2 ): log RðtÞ K À or log RðtÞ ! K þ (see Main Text Figure 2). The boundaries K AE can be expressed in terms of the error e and for symmetric errors, K þ ¼ ÀK À ¼ ð1 À eÞ=e (Siggia and Vergassola, 2013). In our case e ¼ 32%. The mean time for decision can be computed for each discrimination task as the mean firstpassage time of a random walk (in the limit of long decision times). In Appendix 2-figure 1, we show the mean decision-time for different values of e.
Assuming the gene has two levels of activation ON and OFF, the information available to downstream mechanisms is a series of ON times s i and OFF times t j of gene activity. The probability of a concentration hypothesis is then PðL m Þ ¼ Pðft i g; fs j gjL m Þ. If successive ON and OFF times are independent then the probabilities factorize. The log ratio is then written as a function of the ON (respectively OFF) time probability distribution P ON ðt; LÞ (respectively P OFF ðt; LÞ) where J À is the number of ON to OFF switching events and J þ the number of OFF to ON switching events (Siggia and Vergassola, 2013). To understand the dynamics of the log ratio, it is necessary to compute the distributions P ON ðt; LÞ and P OFF ðt; LÞ based on the promoter dynamics, which is the focus of the following subsection.

From binding to gene activation
In this subsection, we describe how the high-dimensional complete state of the promoter is mapped onto the low-dimensional activity of the gene. We use a formalism similar to that of Desponds et al., 2016;Tran et al., 2018. The promoter features N binding sites : its complete state at time T is described by an N dimensional vectorBðtÞ, where B i ¼ 0=1 if the binding site i is empty/full. We assume the gene has two levels of activity: either it is OFF and no mRNA is produced, or it is ON and mRNA is produced at a fixed rate. Activity is then described by a simple Boolean variable aðtÞ equal to 0/1, corresponding to the gene being OFF/ON. We also assume that the activity of the gene depends only on the total number of transcription factors bound to the promoter so that there is an integer 1 k N such that aðtÞ ¼ 1 1 P i Bi!k (where 1 1 is the indicator function). From the point of view of gene activity, we are only interested in the dynamics of BðtÞ ¼ P i B i ðtÞ. We make another simplifying assumption: the probabilities of BðtÞ increasing or decreasing only depend on the value of BðtÞ and not on which binding sites specifically are bound or unbound, that is BðtÞ itself has a Markovian dynamics in f0; 1; :::Ng.
At a given time t, we describe the stochastic state of the promoter as a vector of probabilityXðtÞ where X i ðtÞ is the probability of having B ¼ i À 1 and P i X i ðtÞ ¼ 1. The Markovian dynamics of B translate into an ðN þ 1Þ Â ðN þ 1Þ transition matrix M forX :Xðt þ sÞ ¼ e MsX ðtÞ. M describes the dynamics of the promoter from the point of view of gene activity. If transcription factors do not dimerize or form complexes, then M i;j 6 ¼ 0 only if ji À jj 1 since the probability of two of them binding or unbinding decreases as the square of time for short times.
Let us now relate the statistics of the switching times for the gene activity a to the transition matrix M. Starting from a distribution of OFF statesX 0 we compute the cumulative distribution function of the waiting time t before switching to the ON state where Trans is the transpose,X ON is a vector of 1 for states corresponding to active genes and 0 for states corresponding to inactive genes, and M Ã OFF is a modified version of M where all ON states act as 'sinks' (i.e. all the transition rates from ON states have been set to 0). The expression in Equation 12 computes the cumulative probability of transitions to an ON state before time t, that is P OFF ðtÞ ¼ dPðt tÞ=dt.
For simplicity, we restrict ourselves to promoter structures where there is only one point of entry into the ON or the OFF states (i.e. any switching event from ON to OFF or vice versa will end with the promoter having a specific number of binding sites). Under this hypothesis, the probability vec-torX 0 in Equation 12 is uniquely determined and independent of the specific dynamics of the previous ON or OFF times. We relax these constraints on the promoter structure in Appendix 8. We denote by t OFF (respectively t ON ) the average of the OFF (respectively ON) times for the real concentration L.

Random walk of the log ratio
When the two hypothesized concentrations are very similar to each other as compared to the concentration scale L set by the actual concentration value, jL 1 À L 2 j<<L, the discrimination is hard and requires a large number of binding and unbinding events. In this limit, the random walk log RðtÞ can be approximated by a drift-diffusion process with drift V and diffusion D: where h is a Gaussian white noise. The decision time for the case of symmetric boundaries K þ ¼ ÀK À is a random variable T with mean given by Equation 1 in the main text (Siggia and Vergassola, 2013). hTi depends on the details of the biochemical sensing of the TF concentration and promoter architecture via the drift and diffusivity. Computing V and D in Equation 13 is enough to derive the mean first-passage time to the decision and even its full distribution by computing its Laplace transform as a solution of the drift-diffusion equation using standard techniques of first-passage for Gaussian processes (Siggia and Vergassola, 2013).

Drift
For many switching events J þ >>1 and J À >>1, the sums in Equation 11 can be replaced by continuous averages over binding time distributions. Since J À and J þ differ at most by one, we can also replace the two values J À ' J þ by a unique value J, which is the number of ON-OFF cycles. The times are distributed according to the ON and OFF probability distributions for the real concentration L. At a given large time t, the number J of terms in the sum and the log-ratios appearing in Equation 11 are a priori correlated but Wald's equality (Wald, 1945a) ensures that the average of the sum is the product of the two averages. Since the expected number of cycles grows linearly in time as hJi / t= t ON þ t OFF ð Þ, we conclude that the drift term in Equation 13 is given by: where D KL ðPjjQÞ ¼ R ¥ 0 dsPðsÞ log PðsÞ=QðsÞ ð Þ is the Kullback-Leibler divergence between the two distributions P and Q.
In sum, the drift is determined by how informative the distribution of waiting times in the ON and OFF states are about the concentration differences, with an average rate equal to the inverse of the cycle time t OFF þ t ON . The Kullback-Leibler divergence appearing in Equation 14 is intuitive : it represents the distance between the real concentration and the hypothetical concentration in the space of probabilistic models of switching times. The larger that distance, the easier it becomes to tell the difference between the two distributions through random sampling.

Expansion of the drift for small concentration differences
When the two candidate concentrations are similar, L ' L 1 ' L 2 , the quantities computed in the previous subsection can be expanded at first order in the differences in concentrations dL 1 ¼ L 1 À L and dL 2 ¼ L 2 À L. Starting from Equation 14, we expand the drift V at first and second orders: where ,! PON means that the same operations and integrations are performed for P ON . Due to conservation of probability, the integral of the first and second L-derivatives of P OFF ðt; LÞ vanish. The first-order expansion of V vanishes and we have at first non-vanishing order in jL i À Lj If for instance jL 2 À Lj>jL 1 À Lj then the drift is positive, favouring the concentration L 1 , closer to the real concentration.

Equality between drift and diffusivity
In this subsection, we present different ways of proving the equality between V and D in SPRT when the two hypotheses are close by connecting different approaches. We show that this equality is a general property of random walks in Bayesian belief space. We check that these results are true in the controlled case of one binding site in Appendix 3-figure 1.
First approach: the exit points of the decision process As proved by Wald in the original paper where he introduced sequential probability ratio tests (Wald, 1944), the nature of the test (the ratio between the likelihood of two hypotheses) requires a specific relationship between the error and the boundaries that define the regions of decision. In our case, we assume the same probability of calling L 1 when L 2 is true and calling L 2 when L 1 is true, leading to the definition of only one error level e and two symmetric boundaries K and ÀK. Specifically, Wald shows that K ¼ log 1Àe e (see also Siggia and Vergassola, 2013). When the two hypotheses are close enough that the variations of the log-likelihood can be approximated by a Gaussian process, one can compute the exit probabilities in terms of the drift V and diffusivity D. The equation for the probability of absorption PðXÞ at the upper boundary K is with the boundary condition PðKÞ ¼ 1 and PðÀKÞ ¼ 0. The corresponding solution is as shown in the supplementary material of Siggia and Vergassola, 2013. Setting x ¼ 0 in Equation 18, we find that We note that this probability of absorption is also 1 À e by definition of the error (assuming for instance that L ¼ L 1 ) leading to which in turn gives And so we find that in the limit of close hypotheses, V ¼ D.

Second approach: expansion for small concentration differences
Let us consider the SPRT process between two hypotheses L 1 and L 2 and assume for simplicity that L ¼ L 2 . In this version of the proof, we consider that the difference dL ¼ L 1 À L is small compared to L (i.e dL=L<<1) and expand the expressions for drift and diffusion in increasing orders of dL=L to find that they match. We have shown in the subsection 'Expansion of the drift for small concentration differences' of Appendix 3 that the first non-vanishing term in the expansion of the drift is of order 2 in dL=L and is given by Equation 16. An integral formula for the second moment of the log-likelihood is given in Equation A15 of Carballo-Pacheco et al., 2019 from which we get that the diffusivity of the log ratio is the sum of four terms. The first term is given by D 1 is proportional to V 2 and so is of order ðdL=LÞ 4 . The second term is given by The prefactor V is of order ðdL=LÞ 2 . Expanding the first two terms in the brackets gives the same type of terms as in Equation 15. Because of probability conservation we find that these terms are of the same order as V (i.e ðdL=LÞ 2 ). The last two terms have prefactors s and t respectively which break the argument for vanishing first order terms in Equation 15 and these terms are of order dL=L. Putting the pieces together, we find that D 2 is of order ðdL=LÞ 3 . The third term is given by Because ON and OFF times are independent, cross terms in D 3 are products of two hlog Pð:; L 1 Þ=Pð:; L 2 Þi and are of subleading order ðdL=LÞ 4 . Using the symmetry between P ON and P OFF , we expand one of the square terms in D 3 The fourth term is ÀV 2 which is of order ðdL=LÞ 4 . So we find that, at leading order, D ' D 3 which has the same expansion as V (see Equation 16). And so we recover that for small concentration differences, V ¼ D.
In the case of one binding site with ON rate L and OFF rate n we check that the formula from Carballo-Pacheco et al., 2019 gives the exact expression computed in Siggia and Vergassola, 2013 Replacing L 2 with L in Equation 26 and expanding in dL=L we have We identify the drift computed for the one binding site example in the paragraph 'Mean decision time : connecting drift-diffusion and Wald's approaches' of the main text and find that at first order drift and diffusivity are equal.
Appendix 3-figure 1. The equality V ¼ D holds for close hypotheses. We compare the exact drift and diffusion computed from the formulae derived in Siggia and Vergassola, 2013 for the case of one binding site. We find that drift and diffusion are approximately equal for close concentrations (i. e small dL=L). Parameters are n ¼ 1 s

Third approach: general properties of Bayesian random walks
To understand the origin of the equality between drift and diffusion, we go back to the definition of SPRT as the update of the Bayesian beliefs in two competing hypotheses. We no longer condition the process on one hypothesis being true but rather on the probability that each one is true given the value of the log ratio a certain time point. A very general property of posterior distributions updated with evidence is that their average does not change. This martingale property of belief update can be shown in the following way: assume that the nucleus receives extra evidence with probability pðÞ. Before that evidence comes, the probability that hypothesis 1 is true p 1 has a certain distribution Pðp 1 Þ. Then we have that the probability that hypothesis 1 is true varies as : where we used Bayes rule. The integral of pðjp 1 Þ over sums up to one and we are left with In the context of two hypothetical concentrations, we consider the normalized likelihood of the two hypotheses at time t: Q 2 ðtÞ ¼ P 2 ðtÞ=ðP 1 ðtÞ þ P 2 ðtÞÞ and Q 1 ðtÞ ¼ P 1 ðtÞ=ðP 1 ðtÞ þ P 2 ðtÞÞ, where P i ðtÞ is the likelihood that hypothesis i is true based on the evidence accumulated up to time t. The martingale property translates as where hi a are averages taken assuming hypothesis a is true and DQ is the variation of Q through time (or accumulated evidence). We note that Q 1 ¼ e LðtÞ =ð1 þ e LðtÞ Þ and Q 2 ¼ 1=ð1 þ e LðtÞÞ , where LðtÞ ¼ log P 1 ðtÞ=P 2 ðtÞ is the log-likelihood ratio at a given time. We express the variations of the probabilities in terms of the log-likelihood so that we can connect it to the drift-diffusion parameters. The drift and diffusivity that depend a priori on the real concentration are the average and variance of the variation of L, respectively. We have hDLðtÞi 1 ¼ V ð1Þ Dt, hDLðtÞi 2 ¼ V ð2Þ Dt, h DLðtÞ À hDLi 1 ð Þ 2 i 1 ¼ 2D ð1Þ Dt and h DLðtÞ À hDLi 2 ð Þ 2 i 2 ¼ 2D ð2Þ Dt (where D is the total change over time Dt). When the two concentrations are close, the statistics of the acquisition of evidence are symmetric for L 1 and L 2 and we get that D ð1Þ ¼ D ð2Þ ¼ D and V ð1Þ ¼ ÀV ð2Þ ¼ V. As explained in the paragraph 'Mean decision time : connecting drift-diffusion and Wald's approaches' of the main text, computing the derivatives of Q i with respect to L is straightforward. Plugging these relationships into Equation 30 and gathering terms in V and D we find that from which we derive that V ¼ D.

The first passage time to decision
The first exit time of a drift-diffusion process starting from 0 with two symmetric boundaries at +K and ÀK is given by hTi ¼ K tanhðVK=2DÞ=V (see for instance [Redner, 2001]). Plugging in V ¼ D, we recover Equation 4 of the main text: We check that this approximation is correct for small concentration differences in Appendix 3figure 2a,b and c. As mentioned in the main text, because the hyperbolic tangent is very flat for high values of its argument, we find that the mean decision time depends weakly on corrections to V ¼ D when the error rate e is small (which is equivalent to K large). We check that this result in true in Appendix 3-figure 2d. We note that from K ¼ logð1 À eÞ=e we have tanhðK=2Þ ¼ 1 À 2e. And so we recover the second part of Equation 4 from the main text In one of the original papers about SPRT (Wald, 1945a), Wald derives a similar formula for the total number of events before the decision: where J exit is the number of ON-OFF events when decision happens. Combining Equation 34 with Wald's equality (Wald, 1945b;Durrett, 2010) that states that hTi ¼ hJ exit iðt ON þ t OFF Þ, we recover Equation 33.
Pacheco et al., 2019, but it is sufficient to grasp the consequences of the correlations introduced by a global constraint on the duration of the process, that is, the Q function. Note that such constraint breaks the statistical independence of the various cycles, and introduces dependencies even though the P ON and P OFF factorize. We rewrite the Heaviside function using the Laplace transform form One can remark that the expression is factorized, that is, f is raised to the power n, and even f itself can be interpreted as 'expectation value' of e l pONðs;L2ÞÀpONðs;L1Þ ½ (and same with p OFF ðt; L a Þ) over the 'distribution' e pONðs;LiÞþqOFFðt;LiÞÀsðtþsÞ , that is, the constraint introduces an exponential factor that distorts the weight of the various durations. These elements are consistent with the idea of a random variable per cycle that can be averaged to compute the drift and diffusivity. However, the idea of ignoring correlations between duration of events is not valid because s is a function of l. Indeed, summing the series over n gives f ðs; lÞ=ð1 À f ðs; lÞÞ and the asymptotic behavior is then determined by the first singularity encountered as g is moved to the left for the inverse Laplace transform. For each value of l, this determines a value of s, that is s s ðlÞ (see Carballo-Pacheco et al., 2019 for the complete calculation).
What the above says more qualitatively is that the variables per cycle are correlated because of the global constraint and there is not a single-cycle effective probability distribution that can account for that because of the dependence of s on l. In other words, different moments involve different configurations and distortions of the weights for the durations of the cycles. The only exception is when ds s ðlÞ=dl vanishes, which is the case for two close hypotheses. Indeed, from f ðs s ðlÞ; lÞ ¼ 1, one obtains ds s =dl / qf =ql and it also holds qf =qlðl ¼ 0Þ ¼ h p ON ðs; L 2 Þ À p ON ðs; L 1 Þ ð Þ þ p OFF ðt; L 2 Þ À p ON ðt; L 1 Þ ð Þi. We checked explicitly in the subsection 'Equality between drift and diffusivity' of Appendix 3 that in that limit in the formula from Carballo-Pacheco et al., 2019 the first two terms are negligible and the diffusivity reduces to the variance of the log-likelihoods only.
Indeed, when the two hypothetical concentrations are close and the correlations can be ignored, a diffusivity per cycle can be defined naturally as where D 3 is the average of the squared log-likelihood variation (see Equation 24). We check in Appendix 3-figure 3 that this formula is a good approximation for small concentration differences in the context of one binding site. In that context we have where n is the OFF rate and L is the ON rate. We find that this approximation is better than the V ¼ D approximation only when the cycle times depend weakly on the concentration (i.e, for one binding site, when n is small, see Appendix 3- figure 3).

The mean first-passage time for different architectures
Activation by multiple TF -irreversible binding all-or-nothing models We first consider the all-or-nothing k ¼ 2 two binding site cycle model depicted in Main Text Figure 5a.
In this model, ON events end when the complex formed by two copies of the TF unbinds. It follows that ON times are exponentially distributed and independent of the concentration P ON ðtÞ ¼ ne Ànt . Because they are independent of the concentration, they will not contribute to the log ratios. The mean ON time is given by t ON ¼ 1=n.
The activation time (OFF time) requires two copies of the TF to bind and is the sum of the times it takes for each of them to bind. From a probabilistic point of view, this sum becomes a convolution and the OFF times are given by The average OFF time is t OFF ¼ 1= 1 L þ 1= 2 L.
We can now compute the drift (the ON times do not contribute to information) The calculations above are only valid when the two binding rates are different 1 6 ¼ 2 . Let us now assume that 1 ¼ 2 . The ON time distribution is unchanged but the OFF time distribution is now P OFF ðt; LÞ ¼ Z t 0 ds 1 Le À1Ls 2 Le À1LðtÀsÞ ¼ 2 1 L 2 te À1Lt ¼ gð2; 1= 1 LÞ; where we have identified gð2; 1= 1 LÞ, the standard Gamma distribution with exponent 2 and parameter 1= 1 L. The expression of the drift simplifies as: When the two binding rates m 1 and m 2 are similar but not exactly equal the activation time distribution can be approximated by the Gamma distribution P OFF ðt; LÞ » g 2; 12L 2þ1 . It is then convenient to use Equation 43 with 12L 2þ1 as a parameter.
Activation by multiple TF -irreversible binding with 1-or-more k ¼ 1 activation We now consider the non-equilibrium two binding site 1-or-more k ¼ 1 activation model depicted in Main Text Figure 5b. The OFF times are exponentially distributed P OFF ðt; LÞ ¼ 1 Le À1Lt . They will contribute as in the one binding-site exponential model. The ON times are now a convolution of a concentration-dependent exponential step with rate 2 L and a concentration-independent unbinding with rate n. Their distribution is P ON ðt; LÞ ¼ Z t 0 ds 2 Le À2Ls ne ÀnðtÀsÞ ¼ 2 Ln n À 2 L e À2Lt À e Ànt À Á : The expression is valid for 2 L 6 ¼ n and, as previously, reduces to a g distribution with exponent two in case of equality.
The corresponding drift is e À2Lt À e Ànt À Á log 2 L 1 n n À 2 L 1 n À 2 L 2 2 L 2 n e À2L1t À e Ànt ð Þ e À2L2t À e Ànt ð Þ Á log e À2L1t À e Ànt e À2L2t À e Ànt ! : Appendix 4-figure 1. Mean time to decision across the embryo based on local Bicoid concentration difference. For each position across the embryo, we compute the performance of the fastest identified architecture (k ¼ 1, no Hill coefficient constraint) if it was implemented to distinguish between two concentrations with 10% relative variation, and an average equal to the Bicoid concentration at the AP position denoted on the x-axes. We find that our architecture performs well in the anterior region but is slow in the posterior region, due to low concentrations.
All-or-nothing vs 1-or-more activation In this section, we give a more detailed derivation of the result given in the subsection 'Activation rules' of the main text. Specifically, we compare two activation rules for the same promoter architecture: two binding sites with binding rates m 1 and m 2 . We consider that the two binding sites bind TF independently so that 1 ¼ 2 2 (there are two free targets for the first binding). The two copies of the TF unbind together with unbinding rate n. It is a cycle model. We consider the 1-or-more k ¼ 1 activation rule (Main Text Figure 5b) and the all-or-nothing activation rule (Main Text Figure 5a). We will prove that the fastest activation scheme is the all-or-nothing scheme when 1 L>n and 1-ormore when 1 L<n. This result corresponds to the following intuitive idea : given a choice to associate the second binding with either the first binding or the unbinding for the readout by the cell, it is more informative to convolve it with the fastest of the rates to maximize the information extracted from it. This is in general true if the time per cycle is the same for the different architectures considered.
The total time it takes to go through a cycle is the same for the two architectures: t ON þ t OFF ¼ 1 1 þ 1= 2 þ 1=n, so the difference in performance will come from the amount of information extracted per cycle. We are interested in the limit of similar concentrations, that is, L ¼ L 2 and L 1 ¼ L þ dL, dL<<L. Since D ¼ V and the decision time is inversely proportional to V, it is enough to compare V ðallÞ (the drift in the all-or-nothing scheme) with V ðoneÞ (the drift in the 1-or-more scheme) to determine the fastest architecture.
To compute V ðallÞ , we start from Equation 41, plugging in L 2 ¼ L, L 1 ¼ L þ dL and 1 ¼ 2 2 . We expand for dL=L ! 0: dt2 2 L e À2Lt À e À22Lt À Á log e À22ðLþdLÞt À e À2ðLþdLÞt e À22Lt À e À2Lt ¼ dL L À dL 2 2L 2 þ Z dt2 2 L e À2Lt À e À22Lt À Á log e À2dLt À e À2tðLþ2dLÞ 1 À e À2Lt ' dL L À dL 2 2L 2 þ R dt2 2 L e À2Lt À e À22Lt ð Þ log 1 À dL 2 t 1À2e À 2 Lt 1Àe À 2 Lt þ dL 2 2 t 2 2 1À4e À 2 Lt 1Àe À 2 Lt h i ' dL L À dL 2 2L 2 þ Z dt 2 Le À2Lt À2dL 2 t 1 À 2e À2Lt À Á Â þdL 2 2 2 t 2 1 À 4e À2Lt À Á À dL 2 2 2 t 2 1 À 2e À2Lt ð Þ 2 1 À e À2Lt where z is the Riemann zeta function. Eventually we have V all t ON þ t OFF ð Þ ¼ dL 2 2L 2 4zð3Þ À 5 ½ : While the expansion of V one does not take a simple form in general, it can be computed for the specific value of 1 L ¼ n. When n ¼ L 1 ¼ 2L 2 , L ¼ L 2 and L 1 ¼ L þ dL, expanding to second order Equation 45 becomes dt2 2 L e À2Lt À e À2L2t À Á log e À2ðLþdLÞt À e À2L2t À Á e À2Lt À e À2L2t ð Þ " # ' 2dL L À dL 2 2L 2 þ Z þ¥ 0 dt2 2 L ÀdL 2 te ÀL2t À dL 2 2 2 t 2 2 e ÀL2t 1 À e 2Lt ¼ 2dL L À dL 2 2L 2 À 2dL L À 2 dL 2 L 2 þ 2zð3Þ dL 2 L 2 ¼ V ðallÞ ðt ON þ t OFF Þ: We can now use the above equality and the following arguments to reach the conclusion stated at the beginning of the section. V ðallÞ is independent of n because no information is gained about the concentrations from the step involving unbinding. Conversely, V ðoneÞ is an increasing function of n: as n decreases, the unbinding part takes over in the convolution for the ON times. Since the difference between the two convolutions can only decrease as n decreases, it becomes harder and harder to differentiate the two ON time distributions, their Kullback-Leibler divergence becomes smaller, and the drift term V ðoneÞ is reduced. In sum, for n ¼ 1 L, V ðallÞ ¼ V ðoneÞ ; V ðallÞ is independent of n whilst V ðoneÞ increases with n. We conclude that for n> 1 L, V ðoneÞ >V ðallÞ and the 1-or-more scheme is preferred. Conversely, for n< 1 L, V ðoneÞ <V ðallÞ and the all-or-nothing scheme is preferred.

Comparing one and two binding-site architectures
In this section we compare the performance of a one binding-site architecture (Main Text Figure 5d) to that of an equilibrium all-or-nothing k ¼ 2 two binding-site architecture (Main Text Figure 5c). The motivation is to provide detailed explanations for the results given in subsection 'How the number of binding sites affects decisions' of the main text. To rank their performance, we compare, as in the previous appendix 6, the drift for the two architectures.
To do the comparison, we optimize the two binding site architecture rates m 1 , m 2 and n 1 for a fixed value of the error rate e, the real concentration L, the hypothetical concentrations L 1 and L 2 and the second unbinding rate n 2 . The fixed second unbinding rate sets a time scale for the process. For concreteness, we suppose that the real concentration L ¼ L 2 and decision is hard L 1 ¼ L þ dL, dL=L<<1.
Decisions can be trivially sped up by increasing indefinitely all the rates to very high values. To avoid this, in the optimization process we set an upper bound to be 5s À1 . We choose this upper bound to be smaller than the largest values of n 2 considered (9s À1 ) and larger than the smaller values of n 2 considered (1s À1 ). From a biological point of view, this upper bound for the ON rates can correspond to the diffusion limited arrival rate. For the OFF rates, it can correspond to a minimum bound time required to trigger a downstream mechanism (for instance, the activation of the gene). We check in Appendix 7-figure 2 that curves and effects similar to those discussed below are obtained if the upper bounds are modified.
For all values of parameters, we find that the optimal ON rates are maximal ( 1 ¼ 2 ¼ 5s À1 ). We observed and discussed in subsection 'How the number of binding sites affects decision' of the main text that for certain values of the parameters (n and L), the optimal first unbinding rate n 1 reaches the upper bound 5s À1 while for smaller values the optimal first unbinding rate remains at 0 (Main Text Figure 5f inset). The transition between the two regions is sharp. Setting n 1 to 0 is equivalent to having an effective one binding site model because the promoter never goes back to the 0 state. For that reason, we want to compare the performance of a two binding site model with parameters 1 , 2 , n 1 and n 2 to that of a one binding site model with parameters 2 and n 2 . Since we are in the limit of small concentration differences, the speed of decision is set by V and t ¼ t ON þ t OFF . We denote the mean drifts of the one and the two binding-site models as V ð1Þ and V ð2Þ , respectively. Similarly, t 1 and t 2 are the mean times per cycle of the two models.
The one binding-site architecture activates with rate 2 L and deactivates with rate n 2 . Both waiting times for ON and OFF expression states are exponential. Following results in the paragraph 'Mean decision time : connecting drift-diffusion and Wald's approaches' of the main text, we have t 1 ¼ 1=n 2 þ 1=ð 2 LÞ and V ð1Þ ¼ ð1=n 2 þ 1= 2 LÞ logðL 1 =L 2 Þ þ ðL 2 À L 1 Þ=L ½ . We conclude that : In the two binding site architecture, the ON times are exponentially distributed with rate n 2 . The OFF times are more complex as they can result from several cycles from a promoter with a binding site occupied to an empty promoter before finally switching to two full binding sites and gene activation. We compute P OFF ðt; LÞ using the modified transition matrix of the Markov chain with the ON states acting as sinks as described in Appendix 3. The OFF time distribution is given by the derivative of Equation 12.
In the two binding site model, the transition matrix is a 3 by 3 matrix and can be diagonalized analytically to compute explicitly the exponential in Equation 12 . We set 1 L eq ¼ 2 L eq ¼ n 1 since it is always the case in the identified optimal architectures for the two binding site model (where L eq is the specific value of L for which upper bounds on n 1 and L are equal, L eq ¼ 1 m 3 in Main Text Figure 5f). We compute the distribution of the OFF times explicitly and find P OFF ðL; tÞ ¼ n 1 r 2 þ 8r e À 1 2 n1tð1þ2rþ ffiffiffiffiffiffiffi ffi where r ¼ L=L eq .
Computing the mean of the distribution in Equation 50 , we find that the average time for a cycle in the two binding site architecture is t ð2Þ ¼ 1 þ r r 1 n 1 r þ 1 n 2 : We can now numerically integrate V ð2Þ ¼ t ð2Þ R P OFF ðt; LÞ log P OFF ðt; L 1 Þ=P OFF ðt; LÞ ð Þ. We use this simplified formula to compare the performance of the two architectures and recover the results of Main Text Figure 5f . We show an example for a specific value of L, varying n 2 in Appendix 7-figure 1.
Appendix 7-figure 1. Comparison of the drift for the one and two binding site equilibrium architectures. We compare the drift V ð1Þ of the one binding-site architecture of Main Text Figure 5d and the drift V ð2Þ for the two binding site equilibrium architecture of Main Text Figure 5c. Both have the same parameters L ¼ L 2 ¼ 1m À3 , L 1 ¼ 1:1m À3 , 1 ¼ 2 ¼ 5m 3 s À1 , n 1 ¼ 5s À1 and we vary n 2 . The fastest architecture is the one with the highest absolute value of the drift (lowest on the graph). We recover the transition observed in Main Text Figure 5f between the optimal architectures. when L<0:95L 0 . For such architectures, the ON and OFF production levels will roughly balance each other at the boundary. To ensure a positive RNA level at the boundary, we redefine the ON rate: r Ã ON ¼ r ON À r b >0 so that the basal rate is now a factor proportional to time in both ON time and OFF time RNA contributions. In that sense, the basal rate is a systematic drift of the RNA level that is proportional to time and can be removed from the equations by defining moving boundaries c 1 ðtÞ ¼ c 1 ð0Þ þ r b t and c 2 ðtÞ ¼ c 2 ð0Þ þ r b t. With these redefined rates we now look for architectures such that the drift at the anterior boundary is positive: Specifically, we look for architectures that satisfy these conditions and for which the ratio of the drift over the variance contributions to the log-likelihood is as large as possible to avoid dynamics that are dominated by noise. Finally we look for a set of threshold concentrations that give the most favorable speed-accuracy tradeoff. We show the results of our non-exhaustive search in Appendix 11-figure 1. For one of the architectures that fall below the error level of 32% and the decision time threshold of 3 min, we check the RNA expression profile (Main text Figure 7). For RNA degradation, we use the rate r OFF ½RNA=ð1 þ ½RNAÞ that reduces to r OFF close to the boundary and yields the desired linear in time degradation level, but that is proportional to ½RNA for small values of the concentration. For this function, we compute the average number of RNA molecules in the cell after 3 min of dynamics (Main text Figure 7) and find that the RNA profile has a Hill coefficient of 5.2.