hunchback promoters can readout morphogenetic positional information in less than a minute

We model the hb promoter as six Bcd binding sites [41, 42, 43, 44] that determine the activity of the gene (Fig. 2 a). Bcd molecules bind to and unbind from each of the $i=1,...,6$ sites with rates $\mu_{i}$ and $\nu_{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 section 1.1. The same remark applies for the wide variety of promoter architectures considered in previous works [45, 34]. 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 [45, 34].

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_{\rm ON}=\sum_{i=k}^{6}P_{i}$, for which the gene is activated) and inactive states ($P_{\rm OFF}=\sum_{i=0}^{k-1}P_{i}$, for which the gene is silenced) as illustrated in Fig. 2 b and d.

As Bicoid ligands bind and unbind the promoter (Fig. 2 b), the gene is successively activated and silenced (Fig. 2 c). 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 Ref. [27], 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) (Fig. 2 d). The decision threshold boundaries $K$ are set by the error rate $e$ for making the wrong decision between the hypothetical concentrations : $K=\log\left((1-e)/e\right)$ (see [27] and Appendix C.5). The choice of $K$ or $e$ depends on the level of reproducibility desired for the decision process. We set $K\simeq 0.75$, corresponding to the widely used error rate $e\simeq 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 Section 1.1 above). 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 (Fig. 2 e) and the hb gene spends little time in the active expression state (Fig. 2 f), which leads to a negative drift in the process and favors the lower one of the two possible concentrations (Fig. 2 g).

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 [26] with drift-diffusion, we establish the equality between the drift and diffusion coefficients in decision making space for difficult decision problems, i.e. 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 (Fig. 3 a). For a given architecture, the dynamics of binding/unbinding events and the rules for activation define the two probability distributions $P_{\rm OFF}(t,L)$ and $P_{\rm ON}(s,L)$ for the duration of the OFF and ON times, respectively (Fig. 3b). The two series are denoted $\{t_{i}\}_{1\leq i\leq J^{+}}$ and $\{s_{j}\}_{1\leq j\leq 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 a drift–diffusion equation: $d\log R(t)/dt=V+\sqrt{2D}\eta$, where $V$ is the constant drift, i.e. the bias in favor of one of the two hypotheses, $D$ is the diffusion constant and $\eta$ a standard Gaussian white noise with zero mean and delta-correlated in time [26, 27]. 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 [46, 27]:

$$\langle T\rangle=\frac{K\tanh(VK/2D)}{V}\,.$$

(1)

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 G). 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 [47, 48] 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 $\langle M\rangle\,\langle X_{i}\rangle$.

Wald’s equality applies to our likelihood sum ($\sum_{i}^{M}\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 $(\tau_{\rm ON}+\tau_{\rm OFF})$, where $\tau_{\rm ON}$ is the mean of the distribution of ON times $P_{\rm ON}(t,L)$ and $\tau_{\rm OFF}$ is the mean of the distribution of OFF times $P_{\rm OFF}(s,L)$. The term $(\tau_{\rm ON}+\tau_{\rm OFF})$ determines the average speed at which the system completes an activation/deactivation cycle, while the average $\langle\log R_{i}\rangle$ 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||g)=\int^{\infty}_{0}dt^{\prime}f(t^{\prime})\log\left(f(t^{\prime})/g(t^{\prime})\right)$ 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}$ :

	$\displaystyle V$	$\displaystyle=\frac{1}{(\tau_{\rm ON}+\tau_{\rm OFF})}\big{[}D_{KL}(P_{\rm OFF}(.,L)||P_{\rm OFF}(.,L_{2}))-D_{KL}(P_{\rm OFF}(.,L)||P_{\rm OFF}(.,L_{1}))$		(2)
		$\displaystyle+D_{KL}(P_{\rm ON}(.,L)||P_{\rm ON}(.,L_{2}))-D_{KL}(P_{\rm ON}(.,L)||P_{\rm ON}(.,L_{1}))\big{]}\,.$

Eq. 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$ (Fig. 3 b).

Diffusivity : Why it is more involved to calculate and how we circumvent it. While the drift has the closed simple form in Eq. 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 [49], which is harder to obtain as it depends a priori on the details of the binding site model (see Appendices C.5 and C.7 for details).

We circumvent the calculation of the diffusivity by noting that the same methods used to derive Eq. (1) also yield the probability of first absorption at one of the two boundaries, say $+K$ (see Appendix C.5) :

$$\Pi_{K}=\frac{e^{VK/D}}{1+e^{\frac{KV}{D}}}\,.$$

(3)

By imposing $\Pi_{K}=1-e$, we obtain $VK/D=\log\left((1-e)/e\right)$ and the comparison with the expression of $K=\log\left((1-e)/e\right)$ 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 Appendix C.5). The Appendix 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 Eq. (1) to

$$\langle T\rangle=\frac{K\tanh(K/2)}{V}=\frac{K}{V}\left(1-2e\right)\,,$$

(4)

which is formula $4.8$ in [26] 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 tangent in Eq. (1) is $\simeq K/2=\ln\left(\sqrt{(1-e)/e}\right)$. It follows that for any problem where the error $e\ll 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 Appendix C.6.

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 $\nu$ is independent of TF concentration and the activation rate is exponentially distributed $P_{\rm OFF}(L,t)=kLe^{\-kLt}$. We can explicitly calculate the drift $V=\nu kL/(\nu+kL)(\log(L_{1}/L_{2})+(L_{2}-L_{1})/L)$ (Eq. 2) and expand it for $L_{2}=L$ and $L_{1}=L+\delta L$, at leading order in $\delta L$. Inserting the resulting expression into Eq. 4, we conclude that

$$\langle T\rangle=\frac{\nu+kL}{\nu kL}\frac{2KL^{2}}{\delta L^{2}}\tanh\left(\frac{K}{2}\right),$$

(5)

decreases with increasing relative TF concentration difference $\delta L/L$ and gives a very good approximation of the complete formula (see SI Fig 9 a-c with different values of $\delta L/L$).

Eqs. 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, Fig. 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 (Fig. 3a-d in blue), and the $4$-or-more rule ($k=4$) (Fig. 3a-d in red) for fixed binding and unbinding parameters. Fig. 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 [27] (Fig. 3d). With the rates chosen for Fig. 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$ (Fig. 3 d).

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)=\Pi_{j}P_{\rm ON}(t_{j},L_{1})P_{\rm OFF}(s_{j},L_{1})$ and $P_{2}(t)=\Pi_{j}P_{\rm ON}(t_{j},L_{2})P_{\rm 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 $\Delta t$ remain the same if the expectation is taken with respect to the current $P_{i}(t)$ (see, e.g. [50]). In formulae, this implies that the average variation of the probability $Q_{2}$ over a given time $\Delta t$ that is

$$\langle\Delta Q_{2}\rangle=Q_{1}\langle\Delta Q_{2}\rangle_{1}+Q_{2}\langle\Delta Q_{2}\rangle_{2}\,,$$

(6)

should vanish (see Appendix C.5 for a derivation). Here, $\langle\Delta Q_{2}\rangle_{1}$ is the expected variation of $Q_{2}$ under the assumption that hypothesis $1$ is true and $\langle\Delta Q_{2}\rangle_{2}$ is the same quantity but under the assumption that hypothesis $2$ is true. We notice now that $Q_{2}(t)=\frac{1}{1+e^{{\cal L}(t)}}$, where ${\cal L}=\log R$ is the log-likelihood, and that the drift-diffusion of the log-likelihood implies that $\langle\Delta{\cal L}\rangle_{1}=V\Delta t$, $\langle\Delta{\cal L}\rangle_{2}=-V\Delta t$ and $\langle\left(\Delta{\cal L}-\langle\Delta\mathcal{L}\rangle_{1}\right)^{2}\rangle_{1}=\langle\left(\Delta{\cal L}-\langle\Delta\mathcal{L}\rangle_{2}\right)^{2}\rangle_{2}=2D\Delta t$. By using that $dQ_{2}/d{\cal L}=-Q_{1}Q_{2}$ and $d^{2}Q_{2}/d{\cal L}^{2}=-Q_{1}Q_{2}(Q_{2}-Q_{1})$, we finally obtain that

$$\langle\Delta Q_{2}\rangle=Q_{1}Q_{2}\left(Q_{2}-Q_{1}\right)\Delta t\left[V-D\right]\,,$$

(7)

and imposing $\langle\Delta Q_{2}\rangle=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.

In addition to the requirements imposed by their performance in the decision process (green dashed line in Fig. 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 [40] (Fig. 1a). This requirement mainly impacts and constrains the ratio between binding rates $\mu_{i}$ and unbinding rates $\nu_{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 $\mu_{\rm max}L\simeq 0.124s^{-1}$ (Appendix A). As a result, the binding rates $\mu_{i}$ are limited by diffusion arrivals and the number of available binding sites: $\mu_{i}\leq\mu_{\rm max}(7-i)$ (black dashed line in Fig. 4b), where $L$ is the concentration of Bicoid. This sets the timescale for binding events. In Appendix A, we explore the different measured values and estimates of parameters defining the diffusion limit $\mu_{\rm max}L$ and their influence on the decision time (see Table 1 for all the predictions).

Third, as shown in Fig. 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_{\rm hb}$ is well approximated as a function of the Bicoid concentration $L(x)$ by the Hill function $f_{\rm hb}(x)\simeq L(x)^{H}/(L(x)^{H}+L_{0}^{H})$, where $L_{0}$ is the Bcd concentration at the half-maximum $hb$ expression point and $H$ is the Hill coefficient (Fig. 1a). Experimentally, the measured Hill coefficient for mRNA expression from the WT hb promoter is $H\sim 7-8$ [8, 34]. Recent work [34] suggests that these high values might not be achieved by Bicoid binding sites only. Given current parameter estimates and an equilibrium binding model, [34] shows that a Hill coefficient of $7$ is not achievable within the duration of an early nuclear cycle ($\simeq 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 $6$ equilibrium Bcd binding sites only), we shall explore the space of possible equilibrium promoter architectures limiting the steepness of our profiles to Hill coefficients $H\sim 4-5$.

Using Eqs. 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 Section 1.2.3. We optimize over all possible binding rates $(\mu_{i})_{1\leq i\leq 6}$ ($\mu_{1}$ is the binding rate of the first Bcd ligand and $\mu_{6}$ the binding rate of the last Bcd ligand when $5$ Bcd ligands are already bound to the promoter), and the unbinding rates $(\nu_{i})_{1\leq i\leq 6}$ ($\nu_{1}$ is the unbinding rate of a single Bcd ligand bound to the promoter and $\nu_{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 $k$-or-more activation rule [45, 34]. We use the most recent estimates of biological constants for the hb promoter and Bcd diffusion (see Appendix A) and set the error rate at the border to $32\%$[12, 15]. Reasons for this choice were given in subsection 1.1 and will be revisited in subsection 2.1.1, 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 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 (Fig. 4a).

The precision of a stochastic readout process is defined by two parameters: the resolution of the readout $\delta x$, and the probability of error, which sets the reproducibility of the readout. In Fig. 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 $\langle n_{\rm error}\rangle=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$ [51]. Assuming a $270$s readout time – the total interphase time of nuclear cycle $11$ ([34]) – for the fastest architecture identified above and an error rate of $32\%$, we find that $\langle n_{\rm error}\rangle\simeq 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 $\langle{T}\rangle{}=270$s read-out in a Berg-Purcell-like fixed time scheme (integrated blue area in Fig. 1c).

Conversely, for a given architecture, reducing the error level increases drastically the mean first-passage 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 Fig. 7. The decision can be made in about a minute for $e=32\%$ but requires on average $10$ minutes for $e=10\%$ (Fig. 7). Note that, because the mean first-passage depends simply on the inverse of the drift per cycle (Eq. 4), the relative performance of two architectures is the same for any error rate so that the fastest architectures identified in Fig. 4a are valid for all error levels.

As shown in Ref. [34], 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$). Fig. 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 Fig. 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 $\tau_{\rm ON}+\tau_{\rm 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, Fig. 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 $\nu_{5}$ or $\nu_{6}$ (Fig. 4b).

The distribution of times spent bound to DNA of individual Bicoid molecules is shown in Fig. 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 [52], who found short (mean $\sim 0.629$s and median $\sim 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 (Fig. 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 Ref. [52], 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 [52], which should correspond to specific binding. The agreement between the distributions in Fig. 4d is overall good, and we ascribe discrepancies to the fact that Mir et al. [52] 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 Supplementary Fig. 13 that the two distributions are very similar and hard to distinguish once they are mixed with the non-specific part of the distribution.

In the parameter range of the early fly embryo, the fastest decision-making architectures share the $1$-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$ seconds 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 (Fig. 5 a and b). Both of these considered architectures are out of equilibrium and require energy consumption (as opposed to the two equilibrium architectures of Fig. 5 c and d).

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$ (Fig. 5 a). 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 $\nu$ 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 Fig. 5 a. The activation time is a convolution of the exponential distributions $P_{\rm OFF}(t,L)=\frac{\mu_{1}\mu_{2}L}{\mu_{2}-\mu_{1}}\left(e^{-\mu_{1}Lt}-e^{-\mu_{2}Lt}\right)$. In the simple case when the two binding rates are similar ($\mu_{1}\simeq\mu_{2}$), the OFF times follow a Gamma distribution and the drift and diffusion can be computed analytically (see Appendix D). When the two binding rates are not similar the drift and diffusion must be obtained by numerical integration (see Appendix D).

In the first model described above (Fig. 5 a), deactivation times are independent of the concentration and do not contribute to the information gained per cycle and, as a result, to $V/(\tau_{\rm ON}+\tau_{\rm 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 (Fig. 5 b). The resulting activation times are exponentially distributed and contribute to drift and diffusion as in the simple two state promoter model (Fig. 5 d). The deactivation times are a convolution of the concentration-dependent second binding and the concentration-independent unbinding of the complex and their probability distribution is $P_{\rm ON}(t,L)=\frac{\nu\mu_{2}L}{\nu-\mu_{2}L}\left(e^{-\mu_{2}Lt}-e^{-\nu t}\right)$. Drift and diffusion can be obtained analytically (Appendix D). 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/\nu$) much larger than the second binding time ($1/\mu_{2}L$), no information is gained from deactivation times. In the limit of $\mu_{1}L/\nu\to\infty$, the $k=1$ model reduces to a one binding site exponential model and the two architectures (Fig. 5 b and d) have the same asymptotic performance.

Within the irreversible schemes of Fig. 5 a and Fig. 5 b, 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/(\tau_{\rm ON}+\tau_{\rm 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 $\mu_{2}L$ and the two activation events can be read independently. This regime of parameters favors the $k=1$ rule (Fig.5 e). However, when the concentration $L$ is very large the two binding events happen very fast and for $\mu_{2}L>>\nu$, 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\rightarrow\infty$ and the $1$-or-more $k=1$ activation scheme (Fig.5 b) effectively becomes equivalent to a single binding site promoter (Fig.5 d), making the all-or-nothing $k=2$ activation (Fig.5 a) scheme more informative (Fig.5 e). The fastest decision time architecture systematically convolves the second binding event with the slowest of the other reactions (Fig.5 e), with the transition between the two activation schemes when the other reactions have exactly equal rates ($\mu_{1}L=\nu$ line in Fig. 5 e) (see Appendix E for a derivation).

The above results have been obtained with six binding sites. Motivated by the possibility of building synthetic promoters [53] 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 (Fig. 5 d) 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$, Fig. 5 c). Just like for the $6$ 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_{\rm 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 $\nu_{2}$, we optimize the remaining parameters $\mu_{1}$, $\mu_{2}$ and $\nu_{1}$ for the shortest average decision time.

For high gene deactivation rates $\nu_{2}$, the fastest decision time is achieved by a promoter with one binding site (Fig. 5 f): once one ligand has bound, the promoter never goes back to being completely unbound ($\nu_{1}^{*}$=0 in Fig. 5 c) but toggles between one and two bound TF (Fig. 5 d with $\nu=\nu_{2}$ and $\mu=\mu_{2}$). For lower values of gene deactivation rates $\nu_{2}$, there is a sharp transition to a minimal $\langle{T}\rangle{}$ 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 $2$ state. For very small deactivation $\nu_{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 (Fig. 5 c). 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 (Fig. 5 d, see Appendix F for a more detailed calculation).

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 $\mu_{\rm max}L_{0}\simeq 0.124s^{-1}$ (see black dashed line in Fig. 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 [54], the purpose of this section is to investigate the rationale for this observation by using the models described in Fig. 5.

Naïvely, 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 $\nu$ and the first binding rate $\mu_{1}$ in the $1$-or-more irreversible binding model of Fig. 5 b and we look for the unbinding rate $\mu_{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 $\mu_{2}\leq\mu_{1}$ if the first binding is strong and happens at the diffusion limit. We optimize the mean decision time for $0\leq\mu_{2}\leq\mu_{1}$ (see Fig. 14 for an example) and find a range of parameters where the fastest-decision value $\mu_{2}^{*}<\mu_{1}$ is not as fast as parameter range allows (Fig. 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 $\mu_{2}\leq\mu_{1}$ and the fastest $\langle{T}\rangle{}$ solution always has the fastest possible binding rates.