1.1 Fold-changes are indistinguishable across models

As a first stop on our search for the “right” model of simple repression, let us consider what we can learn from theory and experimental measurements on the average level of gene expression in a population of cells. One experimental strategy that has been particularly useful (if incomplete since it misses out on gene expression dynamics) is to measure the fold-change in mean expression [25]. The fold-change FC is defined as

where angle brackets 〈⋅〉 denote the average over a population of cells and mean mRNA 〈m〉 and mean protein 〈p〉 are viewed as a function of repressor copy number R. What this means is that the fold-change in gene expression is a relative measurement of the effect of the transcriptional repressor (R > 0) on the gene expression level compared to an unregulated promoter (R = 0). The third equality in Eq 2 follows from assuming that the translation efficiency, i.e., the number of proteins translated per mRNA, is the same in both conditions. In other words, we assume that mean protein level is proportional to mean mRNA level, and that the proportionality constant is the same in both conditions and therefore cancels out in the ratio. This is reasonable since the cells in the two conditions are identical except for the presence of the transcription factor, and the model assumes that the transcription factor has no direct effect on translation.

Fold-change has proven a very convenient observable in past work [41–44]. Part of its utility in dissecting transcriptional regulation is its ratiometric nature, which removes many secondary effects that are present when making an absolute gene expression measurement. Also, by measuring otherwise identical cells with and without a transcription factor present, any biological noise common to both conditions can be made to cancel out. Fig 1B and 1C depicts a smorgasbord of mathematicized cartoons for simple repression using both thermodynamic and kinetic models, respectively, that have appeared in previous literature. For each cartoon, we calculate the fold-change in mean gene expression as predicted by that model, deferring most algebraic details to the S1 Supporting Information. What we will find is that for all cartoons the fold-change can be written as a Fermi function of the form

where the effective free energy contains two terms: the parameters ΔF_R, an effective free energy parametrizing the repressor-DNA interaction, and ρ, a term derived from the level of coarse-graining used to model all repressor-free states. In other words, the effective free energy of the Fermi function can be written as the additive effect of the regulation given by the repressor via ΔF_R, and the kinetic scheme used to describe the steps that lead to a transcriptional event via log(ρ) (See Fig 1D, left panel). This implies all models collapse to a single master curve as shown in Fig 1D. We will offer some intuition for why this master curve exists and discuss why at the level of the mean expression, we are unable to discriminate “right” from “wrong” cartoons given only measurements of fold-changes in expression.

1.1.1 Two- and three-state thermodynamic models

We begin our analysis with models 1 and 2 in Fig 1B. In each of these models the promoter is idealized as existing in a set of discrete states; the difference being whether or not the RNAP bound state is included or not. Gene expression is then assumed to be proportional to the probability of the promoter being in either the empty state (model 1) or the RNAP-bound state (model (2)). We direct the reader to the S1 Supporting Information for details on the derivation of the fold-change. For our purposes here, it suffices to state that the functional form of the fold-change for model 1 is

where R is the number of repressors per cell, N_NS is the number of non-specific binding sites where the repressor can bind, Δε_R is the repressor-operator binding energy, and β ≡ (k_BT)⁻¹. This equation matches the form of the master curve in Fig 1D with ρ = 1 and ΔF_R = βΔε_R − log(R/N_NS). For model 2 we have a similar situation. The fold-change takes the form where P is the number of RNAP per cell, and Δε_P is the RNAP-promoter binding energy. For this model we have ΔF_R = βΔε_R − log(R/N_NS) and $$ \rho =1+\frac{P}{N_{NS}}{\mathrm{e}}^{-\beta \Delta {\epsilon }_P}$$. Thus far, we see that the two thermodynamic models, despite making different coarse-graining commitments, result in the same functional form for the fold-change in mean gene expression. We now explore how kinetic models fare when faced with computing the same observable.

1.1.2 Kinetic models

One of the main difference between models shown in Fig 1C, cast in the language of chemical master equations, compared with the thermodynamic models discussed in the previous section is the probability space over which they are built. Rather than keeping track only of the microstate of the promoter, and assuming that gene expression is proportional to the probability of the promoter being in a certain microstate, chemical master equation models are built on the entire probability state of both the promoter microstate, and the current mRNA count. Therefore, in order to compute the fold-change, we must compute the mean mRNA count on each of the promoter microstates, and add them all together [32].

Again, we consign all details of the derivation to the S1 Supporting Information. Here we just highlight the general findings for all five kinetic models. As already shown in Fig 1C and 1D, all the kinetic models explored can be collapsed onto the master curve. Given that the repressor-bound state only connects to the rest of the promoter dynamics via its binding and unbinding rates, $$ k_R^{+}$$ and $$ k_R^{-}$$ respectively, all models can effectively be separated into two categories: a single repressor-bound state, and all other promoter states with different levels of coarse graining. This structure then guarantees that, at steady-state, detailed balance between these two groups is satisfied. What this implies is that the steady-state distribution of each of the non-repressor states has the same functional form with or without the repressor, allowing us to write the fold-change as a product of the ratio of the binding and unbinding rates of the promoter, and the promoter details. This results in a fold-change of the form

where $$ \Delta F_R\equiv -\text{log}(k_R^{+}/k_R^{-})$$, and the functional forms of ρ for each model change as shown in Fig 1C. Another intuitive way to think about these two terms is as follows: in all kinetic models shown in Fig 1C the repressor-bound state can only be reached from a single repressor-free state. The ratio of these two states --repressor-bound and adjacent repressor-free state-- must remain the same for all models, regardless of the details included in other promoter states if ΔF_R represents an effective free energy of the repressor binding the DNA operator. The presence of other states then draws probability density from the promoter being in either of these two states, making the ratio between the repressor-bound state and all repressor-free states different. The log difference in this ratio is given by log(ρ). Since model 1 and model 5 of Fig 1C consist of a single repressor-free state, ρ is then necessarily 1 (See the S1 Supporting Information for further details).

The key outcome of our analysis of the models in Fig 1 is the existence of a master curve shown in Fig 1D to which the fold-change predictions of all the models collapse. This master curve is parametrized by only two effective parameters: ΔF_R, which characterizes the number of repressors and their binding strength to the DNA, and ρ, which characterizes all other features of the promoter architecture. The key assumption underpinning this result is that no transcription occurs when a repressor is bound to its operator. Given this outcome, i.e., the degeneracy of the different models at the level of fold-change, a mean-based metric such as the fold-change that can be readily measured experimentally is insufficient to discern between these different levels of coarse-graining. The natural extension that the field has followed for the most part is to explore higher moments of the gene expression distribution in order to establish if those contain the key insights into the mechanistic nature of the gene transcription process [24, 35]. Following a similar trend, in the next section we extend the analysis of the models to higher moments of the mRNA distribution as we continue to examine the discriminatory power of these different models.

2.1 Kinetic models for unregulated promoter noise

Before we can tackle simple repression, we need an adequate phenomenological model of constitutive expression. The literature abounds with options from which we can choose, and we show several potential kinetic models for constitutive promoters in Fig 2A. Let us consider the suitability of each model for our purposes in turn.

2.1.1 Poisson noise promoter

The simplest model of constitutive expression that we can imagine is shown as model 1 in Fig 2A and assumes that transcripts are produced as a Poisson process from a single promoter state. This is the picture from Jones et. al. [33] that was used to interpret a systematic study of gene expression noise over a series of promoters designed to have different strengths, but no regulation. This model insists that the “true” steady-state mRNA distribution is Poisson, implying the Fano factor ν must be 1. In [33], the authors carefully attribute measured deviations from Fano = 1 to intensity variability in fluorescence measurements, gene copy number variation, and copy number fluctuations of the transcription machinery, e.g., RNAP itself. In this picture, all the corrections to Poisson behavior are derived as additive corrections to the Fano factor. This picture is appealing in its simplicity, with only two parameters, the initiation rate r and degradation rate γ. In other words, the model is not excessively complex for the data at hand. But for many interesting questions, for instance in the recent work [47], attributing all deviations from the model to extrinsic noise sources, limits the kinds of predictions that can be done. To make progress then we need a (slightly) more complex model than model 1 that would allow us to incorporate the non-Poissonian features of constitutive promoters directly into a master equation formulation.

2.1.2 Sub-Poissoninan noise promoters with RNAP escape

A natural extension of the one-state promoter studied in the previous section is to explicitly include an empty promoter state. This state allows for single RNAP to bind and unbind from the promoter with rates $$ k_P^{+}$$ and $$ k_P^{-}$$, respectively, before engaging in a transcriptional event. Once the RNAP irreversibly escapes from the promoter to transcribe the mRNA, the promoter goes back to this empty state. Models 2 and 3 from Fig 2A show two different versions of this process. In model 2 the polymerase can directly engage in a transcription event with rate r after binding to the promoter, while model 3 includes an irreversible step with rate k_O that coarse-grains the known multi-step process involved in transcriptional initiation [49].

Although physically intuitive given our understanding of how a single transcript is produced, we are forced to discard both of these models on the basis that the Fano factor for both of them is strictly < 1 as shown in Fig 2A (See S1 Supporting Information for detailed derivations). This is in disagreement with the experimental data from [33], reproduced in Fig 2B. All of these unregulated promoters –which only differ in their promoter sequence– show strictly super-Poissonian noise (Fano factor > 1). In fact, we suspect any model in which transcription proceeds through a multistep cycle must necessarily have ν < 1. The intuitive argument compares the waiting time distribution to traverse the cycle with the waiting time for a Poisson promoter (model 1) with the same mean time. The latter is simply an exponential distribution. The former is a convolution of multiple exponentials, and intuitively the waiting time distribution for a multistep process should be more peaked with a smaller fractional width than a single exponential with the same mean. A less disperse waiting time distribution means transcription initiations are more uniformly distributed in time relative to a Poisson process. Hence the distribution of mRNA over a population of cells should be less variable compared to Poisson, giving ν < 1. (In the S1 Supporting Information we present a more precise version of the intuitive arguments in this paragraph).

2.1.3 Noise in a two-state promoter with “active” and “inactive” states

Inspired by [47], we next revisit an analog of model 2 in Fig 2A, but as with the analogous models considered in Section 1, the interpretation of the two states is changed. Rather than explicitly viewing them as RNAP bound and unbound, we view them as “active” and “inactive,” which are able and unable to initiate transcripts, respectively. We are noncommittal as to the microscopic details of these states.

One interpretation [37, 38, 50] for the active and inactive states is that they represent the promoter’s supercoiling state: transitions to the inactive state are caused by accumulation of positive supercoiling, which inhibits transcription, and transitions back to “active” are caused by gyrase or other topoisomerases relieving the supercoiling. This is an interesting possibility because it would mean the timescale for promoter state transitions is driven by topoisomerase kinetics, not by RNAP kinetics. From in vitro measurements, the former are suggested to be of order minutes [37]. Contrast this with model 2, where the state transitions are assumed to be governed by RNAP, which, assuming a copy number per cell of order 10³, has a diffusion-limited association rate k_on ∼ 10² s⁻¹ to a target promoter. Combined with known K_d’s of order μM, this gives an RNAP dissociation rate k_off of order 10² s⁻¹. As we will show below, however, there are some lingering puzzles with interpreting this supercoiling hypothesis, so we leave it as a speculative hypothesis and refrain from assigning definite physical meaning to the two states in this model.

Intuitively one might expect that, since transcripts are produced as a Poisson process only when the promoter is in one of the two states in this model, transcription initiations should now be “bunched” in time, in contrast to the “anti-bunching” of models 2 and 3 above. One might further guess that this bunching would lead to super-Poissonian noise in the mRNA distribution over a population of cells. Indeed, as shown in the S1 Supporting Information, a calculation of the Fano factor produces

which is strictly greater than 1, verifying the above intuition. Note we have dropped the P label on the promoter switching rates to emphasize that these very likely do not represent kinetics of RNAP itself. This calculation can also be sidestepped by noting that the model is mathematically equivalent to the simple repression model from [33], with states and rates relabeled and reinterpreted.

How does this model compare to model 1 above? In model 1, all non-Poisson features of the mRNA distribution were handled as extrinsic corrections. By contrast, here the 3 parameter model is used to fit the full mRNA distribution as measured in mRNA FISH experiments. In essence, all variability in the mRNA distribution is regarded as “intrinsic,” arising either from stochastic initiation or from switching between the two coarse-grained promoter states. The advantage of this approach is that it fits neatly into the master equation picture, and the parameters thus inferred can be used as input for more complicated models with regulation by transcription factors.

While this seems promising, there is a major drawback for our purposes which was already uncovered by the authors of [47]: the statistical inference problem is nonidentifiable, in the sense described in Section 4.3 of [51]. What this means is that it is impossible to infer the parameters r and k⁻ from the single-cell mRNA counts data of [33] (as shown in S2 Fig of [47]). Rather, only the ratio r/k⁻ could be inferred. In that work, the problem was worked around with an informative prior on the ratio k⁻/k⁺. That approach is unlikely to work here, as, recall, our entire goal in modeling constitutive expression is to use it as the basis for a yet more complicated model, when we add on repression. But adding more complexity to a model that is already poorly identified is a fool’s errand, so we will explore one more potential model.

2.1.4 Noise model for one-state promoter with explicit bursts

The final model we consider is inspired by the failure mode of model 4. The key observation above was that, as found in [47], only two parameters, k⁺ and the ratio r/k⁻, could be directly inferred from the single-cell mRNA data from [33]. So let us take this seriously and imagine a model where these are the only two model parameters. What would this model look like?

To develop some intuition, consider model 4 in the limit k⁺ ≪ k⁻ ≲ r, which is roughly satisfied by the parameters inferred in [47]. In this limit, the system spends the majority of its time in the inactive state, occasionally becoming active and making a burst of transcripts. This should call to mind the well-known phenomenon of transcriptional bursting, as observed in, e.g., [37, 38, 52–54]. Let us make this correspondence more precise. The mean dwell time in the active state is 1/k⁻. While in this state, transcripts are produced at a rate r per unit time. So on average, r/k⁻ transcripts are produced before the system switches to the inactive state. Once in the inactive state, the system dwells there for an average time 1/k⁺ before returning to the active state and repeating the process. r/k⁻ resembles an average burst size, and 1/k⁺ resembles the time interval between burst events. More precisely, 1/k⁺ is the mean time between the end of one burst and the start of the next, whereas 1/k⁺ + 1/k⁻ would be the mean interval between the start of two successive burst events, but in the limit k⁺ ≪ k⁻, 1/k⁺ + 1/k⁻ ≈ 1/k⁺. Note that this limit ensures that the waiting time between bursts is approximately exponentially distributed, with mean set by the only timescale left in the problem, 1/k⁺. If instead it were the case that k⁺ ∼ k⁻, then the waiting time 1/k⁺ + 1/k⁻ would have a peaked distribution, but this does not appear to be the case for any of datasets from [33].

Let us now verify this intuition with a precise derivation to check that r/k⁻ is in fact the mean burst size and to obtain the full burst size distribution. Consider first a constant, known dwell time T in the active state. Transcripts are produced at a rate r per unit time, so the number of transcripts n produced during T is Poisson distributed with mean rT, i.e.,

Since the dwell time T is unobservable from single molecule mRNA counts, we actually want P(n), the dwell time distribution with no conditioning on T. Basic rules of probability theory tell us we can write P(n) in terms of P(n ∣ T) as But we know the dwell time distribution P(T), which is exponentially distributed according to so P(n) can be written as A standard integral table shows $$ {\int }_0^{\infty }{x}^n{e}^{-ax}dx=n!/a^{n+1}$$, so which is exactly the geometric distribution with standard parameter θ ≡ k⁻/(k⁻ + r) and domain n ∈ {0, 1, 2, …}. The mean of the geometric distribution, with this convention, is exactly as we guessed intuitively above.

So in taking the limit r, k⁻ → ∞, r/k⁻ ≡ b, we obtain a model which effectively has only a single promoter state, which initiates bursts at rate k⁺ (transitions to the active state, in the model 4 picture). The master equation for mRNA copy number m as derived in the S1 Supporting Information takes the form

where we use k_i to denote the burst initiation rate, Geom(n;b) is the geometric distribution with mean b, i.e., $$ \text{Geom}(n;b)=\frac{1}{1+b}{\left(\frac{b}{1+b}\right)}^n$$ (with domain over nonnegative integers as above). The first two terms are the usual mRNA degradation terms. The third term enumerates all ways the system can produce a burst of transcripts and arrive at copy number m, given that it had copy number m′ before the burst. The fourth term allows the system to start with copy number m, then produce a burst and end with copy number m′. In fact this last sum has trivial m′ dependence and simply enforces normalization of the geometric distribution. Carrying it out we have We direct readers again to the S1 Supporting Information for further details. This improves on model 4 in that now the parameters are easily inferred, as we will see later, and have clean interpretations. The non-Poissonian features are attributed to the empirically well-established phenomenological picture of bursty transcription.

The big approximation in going from model 4 to 5 is that a burst is produced instantaneously rather than over a finite time. If the true burst duration is not short compared to transcription factor kinetic timescales, this could be a problem in that mean burst size in the presence and absence of repressors could change, rendering parameter inferences from the constitutive case inappropriate. Let us make some simple estimates of this.

Consider the time delay between the first and final RNAPs in a burst initiating transcription (not the time to complete transcripts, which potentially could be much longer.) If this timescale is short compared to the typical search timescale of transcription factors, then all is well. The estimates from deHaseth et. al. [49] put RNAP’s diffusion-limited on rate around ∼ few × 10⁻² nM⁻¹ s⁻¹ and polymerase loading as high as 1 s⁻¹. Then for reasonable burst sizes of < 10, it is reasonable to guess that bursts might finish initiating on a timescale of tens of seconds or less (with another 30-60 sec to finish elongation, but that does not matter here). A transcription factor with typical copy number of order 10 (or less) would have a diffusion-limited association rate of order (10 sec)⁻¹ [55]. Higher copy number TFs tend to have many binding sites over the genome, which should serve to pull them out of circulation and keep their effective association rates from rising too large [56]. Therefore, there is perhaps a timescale separation possible between transcription factor association rates and burst durations, but this assumption could very well break down, so we will have to keep it in mind when we infer repressor rates from the Jones et. al. single-cell mRNA counts data later [33].

In reflecting on these 5 models, the reader may feel that exploring a multitude of potential models just to return to a very minimal phenomenological model of bursty transcription may seem highly pedantic. But the purpose of the exercise was to examine a host of models from the literature and understand why they are insufficient, one way or another, for our purposes. Along the way we have learned that the detailed kinetics of RNAP binding and initiating transcription are probably irrelevant for setting the population distribution of mRNA. The timescales are simply too fast, and as we will see later in Figs 3 and 4, the noise seems to be governed by slower timescales. Perhaps in hindsight this is not surprising: intuitively, the degradation rate γ sets the fundamental timescale for mRNA dynamics, and any other processes that substantially modulate the mRNA distribution should not differ from γ by orders of magnitude.

3.1 Parameter inference for constitutive promoters

From consideration of Fano factors in the previous section, we suspect that model 5 in Fig 2A, a one-state bursty model of constitutive promoters, achieves the right balance of complexity and simplicity, by allowing both Fano factor ν > 1, but also by remedying, by design, the problems of parameter degeneracy that model 4 in Fig 2A suffered [47]. Does this stand up to closer scrutiny, namely, comparison to full mRNA distributions rather than simply their moments? We will test this thoroughly on single-cell mRNA counts for different unregulated promoters from Jones et. al. [33].

It will be instructive, however, to first consider the Poisson promoter, model 1 in Fig 2. As we alluded to earlier, since the Poisson distribution has a Fano factor ν strictly equal to 1, and all of the observed data in Fig 2B has Fano factor ν > 1, we might already suspect that this model is incapable of fitting the data. We will verify that this is in fact the case. Using the same argument we can immediately rule out models 2 and 3 from Fig 2A. These models have Fano factors ν ≤ 1 meaning they are underdispersed relative to the Poisson distribution. We will also not explicitly consider model 4 from Fig 2A since it was already thoroughly analyzed in [47], and since model 5 can be viewed as a special case of it.

Our objective for this section will then be to assess whether or not model 5 is quantitatively able to reproduce experimental data. In other words, if our claim is that the level of coarse graining in this model is capable of capturing the relevant features of the data, then we should be able to find values for the model parameters that can match theoretical predictions with single-molecule mRNA count distributions. A natural language for this parameter inference problem is that of Bayesian probability. We will then build a Bayesian inference pipeline to fit the model parameters to data. We will use the full dataset of single-cell mRNA counts from [33] used in Fig 2B. To gain intuition on how this analysis is done we analyzed the “wrong” model 1 in Fig 2A. For this analysis, we direct the reader to the S1 Supporting Information as an intuitive introduction to the Bayesian framework with a simpler model.

3.1.1 Model 5—bursty promoter

Let us now consider the problem of parameter inference for model five. As derived in the S1 Supporting Information, the steady-state mRNA distribution in this model is a negative binomial distribution [36], given by

where b is the mean burst size and k_i is the burst rate in units of the mRNA degradation rate γ. As sketched earlier, to think of the negative binomial distribution in terms of an intuitive “story,” in the precise meaning of [57], we imagine the arrival of bursts as a Poisson process with rate k_i, where each burst has a geometrically-distributed size with mean size b.

As for the Poisson promoter model, this expression for the steady-state mRNA distribution is exactly the likelihood we want to use when stating Bayes theorem. Denoting the single-cell mRNA count data as D = {m₁, m₂, …, m_N}, Bayes’ theorem takes the form

We already have our likelihood --the product of N negative binomials as Eq 18—so we only need to choose priors on k_i and b. For the datasets from [33] that we are analyzing, as for the Poisson promoter model above, we are still data-rich so the prior’s influence remains weak, but not nearly as weak because the dimensionality of our model has increased from one parameter to two. Details on the arguments behind our prior distribution selection are left for the S1 Supporting Information. We state here that the natural scale to explore these parameters is logarithmic. This is commonly the case for parameters for which our previous knowledge based on our domain expertise spans several orders of magnitude. For this we chose log-normal distributions for both k_i and b. Details on the mean and variance of these distributions can be found in the S1 Supporting Information.

We carried out Markov-Chain Monte Carlo (MCMC) sampling on the posterior of this model, starting with the constitutive lacUV5 dataset from [33]. The resulting MCMC samples are shown in Fig 3A. In contrast to the active/inactive constitutive model considered in [47] (kinetic model 4 in Fig 2A), this model is well-identified with both parameters determined to a fractional uncertainty of 5-10%. The strong correlation reflects the fact that their product sets the mean of the mRNA distribution, which is tightly constrained by the data, but there is weak “sloppiness” [58] along a set of values with a similar product.

Having found the model’s posterior to be well-identified as with the Poisson promoter (S1 Supporting Information), the next step is to compare both models with experimental data. To do this for the case of the bursty promoter, for each of the parameter samples shown in Fig 3A we generated negative bionomial-distributed mRNA counts. As MCMC samples parameter space proportionally to the posterior distribution, this set of random samples span the range of possible values that we would expect given the correspondence between our theoretical model and the experimental data. A similar procedure can be applied to the Poisson promoter. To compare so many samples with the actual observed data, we can use empirical cumulative distribution functions (ECDF) of the distribution quantiles. This representation is shown in Fig 3B. In this example, the median for each possible mRNA count for the Poisson distribution is shown as a dark green line, while the lighter green contains 95% of the randomly generated samples. This way of representing the fit of the model to the data gives us a sense of the range of data we might consider plausible, under the assumption that the model is true. For this case, as we expected given our premise of the Poisson promoter being wrong, it is quite obvious that the observed data, plotted in black is not consistent with the Poisson promoter model. An equivalent plot for the bursty promoter model is shown in blue. Again the darker tone shows the median, while the lighter color encompasses 95% of the randomly generated samples. Unlike the Poisson promoter model, the experimental ECDF closely tracks the posterior predictive ECDF, indicating this model is actually able to generate the observed data and increasing our confidence that this model is sufficient to parametrize the physical reality of the system.

The commonly used promoter sequence lacUV5 is our primary target here, since it forms the core of all the simple repression constructs of [33] that we consider in Section 3.2. Nevertheless, we thought it wise to apply our bursty promoter model to the other 17 unregulated promoters available in the single-cell mRNA count dataset from [33] as a test that the model is capturing the essential phenomenology. If the model fit well to all the different promoters, this would increase our confidence that it would serve well as a foundation for inferring repressor kinetics later in Section 3.2. Conversely, were the model to fail on more than a couple of the other promoters, it would give us pause.

Fig 3C shows the results, plotting the posterior distribution from individually MCMC sampling all 18 constitutive promoter datasets from [33]. To aid visualization, rather than plotting samples for each promoter’s posterior as in Fig 3A, for each posterior we find and plot the curve that surrounds the 95% highest probability density region. What this means is that each contour encloses approximately 95% of the samples, and thus 95% of the probability mass, of its posterior distribution. Theory-experiment comparisons, shown in the S1 Supporting Information, display a similar level of agreement between data and predictive samples as for the bursty model with lacUV5 in Fig 3B.

One interesting feature from Fig 3C is that burst rate varies far more widely, over a range of ∼ 10², than burst size, confined to a range of ≲ 10¹ (and with the exception of promoter 6, just a span of 3 to 5-fold). This suggests that k_i, not b, is the key dynamic variable that promoter sequence tunes.

3.1.2 Connecting inferred parameters to prior work

It is interesting to connect these inferences on k_i and b to the work of [48], where these same 18 promoters were considered through the lens of the three-state thermodynamic model (model 2 in Fig 1B) and binding energies Δε_P were predicted from an energy matrix model derived from [59]. As previously discussed the thermodynamic models of gene regulation can only make statements about the mean gene expression. This implies that we can draw the connection between both frameworks by equating the mean mRNA 〈m〉. This results in

By taking the weak promoter approximation for the thermodynamic model (P/N_NS exp(−βΔε_P) ≪ 1) results in [48] valid for all the binding energies considered here.

Given this result, how are the two coarse-grainings related? A quick estimate can shed some light. Consider for instance the lacUV5 promoter, which we see from Fig 3A has k_i/γ ∼ b ∼ few, from Fig 3B has 〈m〉 ∼ 20, and from [48] has βΔε_P ∼ −6.5. Further we generally assume P/N_NS ∼ 10⁻³ since N_NS ≈ 4.6 × 10⁶ and P ∼ 10³. After some guess-and-check with these values, one finds the only association that makes dimensional sense and produces the correct order-of-magnitude for the known parameters is to take

and Fig 3D shows that this linear scaling between ln k_i and −βΔε_P is approximately true for all 18 constitutive promoters considered. The plotted line is simply Eq 22 and assumes P ≈ 5000.

While the associations represented by Eqs 22 and 23 appear to be borne out by the data in Fig 3, we do not find the association of parameters they imply to be intuitive. We are also cautious to ascribe too much physical reality to the parameters. Indeed, part of our point in comparing the various constitutive promoter models is to demonstrate that these models each provide an internally self-consistent framework that adequately describes the data, but attempting to translate between models reveals the dubious physical interpretation of their parameters.

Our remaining task in this work is a determination of the physical reality of thermodynamic binding energies in Fig 1, as codified by the thermodynamic-kinetic equivalence of Eq 3. For our phenomenological needs here model 5 in Fig 2A is more than adequate: the posterior distributions in Fig 3C are cleanly identifiable and the predictive checks in the S1 Supporting Information indicate no discrepancies between the model and the mRNA single-molecule count data of [33]. Of the models we have considered it is unique in satisfying both these requirements. As a result we will use it as a foundation to build upon in the next section when we add regulation.

3.2 Transcription factor kinetics can be inferred from single-cell mRNA distribution measurements

3.2.1 Building the model and performing parameter inference

Now that we have a satisfactory model in hand for constitutive promoters, we would like to return to the main thread: can we reconcile the thermodynamic and kinetic models by putting to the test Eq 3, the correspondence between indirectly inferred thermodynamic binding energies and kinetic rates? To make this comparison, is it possible to infer repressor binding and unbinding rates from mRNA distributions over a population of cells as measured by single-molecule Fluorescence in situ Hybridization in [33]? If so, how do these inferred rates compare to direct single-molecule measurements such as from [55] and to binding energies such as from [41] and [43], which were inferred under the assumptions of the thermodynamic models in Fig 1B? And can this comparison shed light on the unreasonable effectiveness of the thermodynamic models, for instance, in their application in [44, 60]?

As we found in Section 2, for our purposes the “right” model of a constitutive promoter is the bursty picture, model five in Fig 2A. Therefore our starting point here is the analogous model with repressor added, model 5 in Fig 1C. For a given repressor binding site and copy number, this model has four rate parameters to be inferred: the repressor binding and unbinding rates $$ k_R^{+}$$, and $$ k_R^{-}$$, the initiation rate of bursts, k_i, and the mean burst size b (we nondimensionalize all of these by the mRNA degradation rate γ).

Before showing the mathematical formulation of our statistical inference model, we would like to sketch the intuitive structure. The dataset from [33] we consider consists of single-cell mRNA counts data of nine different conditions, spanning several combinations of three unique repressor binding sites (the so-called Oid, O1, and O2 operators) and four unique repressor copy numbers. We assume that the values of k_i and b are known, since we have already cleanly inferred them from constitutive promoter data, and further we assume that these values are the same across datasets with different repressor binding sites and copy numbers. In other words, we assume that the regulation of the transcription factor does not affect the mean burst size nor the burst initiation rate. The regulation occurs as the promoter is taken away from the transcriptionally active state when the promoter is bound by repressor. We assume that there is one unbinding rate parameter for each repressor binding site, and likewise one binding rate for each unique repressor copy number. This makes our model seven dimensional, or nine if one counts k_i and b as well. Note that we use only a subset of the datasets from Jones et. al. [33], as discussed more in the S1 Supporting Information.

Formally now, denote the set of seven repressor rates to be inferred as

where subscripts for dissociation rates k⁻ indicate the different repressor binding sites, and subscripts for association rates k⁺ indicate the concentration of the small-molecule that controlled the expression of the LacI repressor (see the S1 Supporting Information). This is because for this particular dataset the repressor copy numbers were not measured directly, but it is safe to assume that a given concentration of the inducer resulted in a specific mean repressor copy number [60]. Also note that the authors of [33] report estimates of LacI copy number per cell rather than direct measurements. However, these estimates were made assuming the validity of the thermodynamic models in Fig 1, and since testing these models is our present goal, it would be circular logic if we were to make the same assumption. Therefore we will make no assumptions about the LacI copy number for a given inducer concentration.

Having stated the problem, Bayes’ theorem reads

where D is again the set of all N observed single-cell mRNA counts across the various conditions. We assume that individual single-cell measurements are independent so that the likelihood factorizes as where $$ k_j^{\pm }$$ represent the appropriate binding and unbinding rates out of $$ \overrightarrow{k}$$ for the j-th measured cell. The probability $$ p(m\mid k_j^{+},k_j^{-},k_i,b)$$ appearing in the last expression is exactly Eq. S205, the steady-state distribution for our bursty model with repression derived in the S1 Supporting Information, which for completeness we reproduce here as where ₂F₁ is the confluent hypergeometric function of the second kind and α and β, defined for notational convenience, are

This likelihood is rather inscrutable. We did not find any of the known analytical approximations for ₂F₁ useful in gaining intuition, so we instead resorted to numerics. One insight we found was that for very strong or very weak repression, the distribution in Eq 27 is well approximated by a negative binomial with burst size b and burst rate k_i equal to their constitutive lacUV5 values, except with k_i multiplied by the fold-change $$ {(1+k_R^{+}/k_R^{-})}^{-1}$$. In other words, once again only the ratio $$ k_R^{+}/k_R^{-}$$ was detectable. But for intermediate repression, the distribution was visibly broadened with Fano factor greater than 1 + b, the value for the corresponding constitutive case. This indicates that the repressor rates had left an imprint on the distribution, and perhaps intuitively, this intermediate regime occurs for values of $$ k_R^{\pm }$$ comparable to the burst rate k_i. Put another way, if the repressor rates are much faster or much slower than k_i, then there is a timescale separation and effectively only one timescale remains, $$ k_i{(1+k_R^{+}/k_R^{-})}^{-1}$$. Only when all three rates in the problem are comparable does the mRNA distribution retain detectable information about them.

Next we specify priors. As for the constitutive model, weakly informative log-normal priors are a natural choice for all our rates. We found that if the priors were too weak, our MCMC sampler would often become stuck in regions of parameter space with very low probability density, unable to move. We struck a balance in choosing our prior widths between helping the sampler run while simultaneously verifying that the marginal posteriors for each parameter were not artificially constrained or distorted by the presence of the prior. All details for our prior distributions are listed in the S1 Supporting Information.

We ran MCMC sampling on the full nine dimensional posterior specified by this model. To attempt to visualize this object, in Fig 4A we plot several two-dimensional slices as contour plots, analogous to Fig 3C. Each of these nine slices corresponds to the $$ (k_R^{+},k_R^{-})$$ pair of rates for one of the conditions from the dataset used to fit the model and gives a sense of the uncertainty and correlations in the posterior. We note that the 95% uncertainties of all the rates span about ∼ 0.3 log units, or about a factor of two, with the exception of $$ k_{0.5}^{+}$$, the association rate for the lowest repressor copy number which is somewhat larger.

3.2.2 Comparison with prior measurements of repressor binding energies

Our primary goal in this work is to reconcile the kinetic and thermodynamic pictures of simple repression. Towards this end we would like to compare the repressor kinetic rates we have inferred with the repressor binding energies inferred through multiple methods in [41] and [43]. If the agreement is close, then it suggests that the thermodynamic models are not wrong and the repressor binding energies they contain correspond to physically real free energies, not mere fit parameters.

Fig 4B shows both comparisons, with the top panel comparing to thermodynamic binding energies and the bottom panel comparing to single-molecule measurements. First consider the top panel and its comparison between repressor kinetic rates and binding energies. As described in section 1, if the equilibrium binding energies from [41] and [43] indeed are the physically real binding energies we believe them to be, then they should be related to the repressor kinetic rates via Eq 3, which we restate here,

Assuming mass action kinetics implies that $$ k_R^{+}$$ is proportional to repressor copy number R, or more precisely, it can be thought of as repressor copy number times some intrinsic per molecule association rate. But since R is not directly known for our data from [33], we cannot use this equation directly. Instead we can consider two different repressor binding sites and compute the difference in binding energy between them, since this difference depends only on the unbinding rates and not on the binding rates. This can be seen by evaluating Eq 29 for two different repressor binding sites, labeled (1) and (2), but with the same repressor copy number R, and taking the difference to find or simply The left and right hand sides of this equation are exactly the horizontal and vertical axes of the top panel of Fig 4. Since we inferred rates for three repressor binding sites (O1, O2, and Oid), there are only two independent differences that can be constructed, and we arbitrarily chose to plot O2-O1 and O1-Oid in Fig 4B. Numerically, we compute values of $$ k_{O1}^{-}/k_{Oid}^{-}$$ and $$ k_{O2}^{-}/k_{O1}^{-}$$ directly from our full posterior samples, which conveniently provides uncertainties as well, as detailed in the S1 Supporting Information. We then compare these log ratios of rates to the binding energy differences Δε_O1 − Δε_Oid and from Δε_O2 − Δε_O1 as computed from the values from both [41] and [43]. Three of the four values are within ∼ 0.5 k_BT of the diagonal representing perfect agreement, which is comparable to the ∼ few × 0.1 k_BT variability between the independent determinations of the same quantities between [41] and [43]. The only outlier involves Oid measurements from [43], and as the authors of [43] note, this is a difficult measurement of low fluorescence signal against high background since Oid represses so strongly. We are therefore inclined to regard the failure of this point to fall near the diagonal as a testament to the difficulty of the measurement and not as a failure of our theory.

On the whole then, we find these to be results very suggestive. Although not conclusively, our results point at the possible interpretation of the phenomenological free energy of transcription factor binding as the log of a ratio of transition rates. This did not need to be the case; the free energy was inferred from bulk measurements of mean gene expression relative to an unregulated strain –the fold-change– while the rates were extracted from fitting single-molecule counts to the full mRNA steady state distribution. Further single molecule experiments along the lines of [55] could shed further light into this idea.

3.2.3 Comparison with prior measurements of repressor kinetics

In the previous section we established the equivalence between the equilibrium models’ binding energies and the repressor kinetics we infer from mRNA population distributions. But one might worry that the repressor rates we infer from mRNA distributions are themselves merely fitting parameters and that they do not actually correspond to the binding and unbinding rates of the repressor in vivo. To verify that this is not the case, we next compare our kinetic rates with a different measurement of the same rates using a radically different method: single molecule measurements as performed in Hammar et. al. [55]. This is plotted in the lower panel of Fig 4B.

Since we do not have access to repressor copy number for either the single-cell mRNA data from [33] or the single-molecule data from [55], we cannot make an apples-to-apples comparison of association rates $$ k_R^{+}$$. Further, while Hammar et. al. directly measure the dissociation rates $$ k_R^{-}$$, our inference procedure returns $$ k_R^{-}/\gamma $$, i.e., the repressor dissociation rate nondimensionalized by the mRNA degradation rate γ. So to make the comparison, we must make an assumption for the value of γ since it was not directly measured. For most mRNAs in E. coli, quoted values for the typical mRNA lifetime γ⁻¹ range between about 2.5 min [61] to 8 min. We chose γ⁻¹ = 3 min and γ⁻¹ = 5 min as representative values and plot a comparison of $$ k_{O1}^{-}$$ and $$ k_{Oid}^{-}$$ from our inference with corresponding values reported in [55] for both these choices of γ.

The degree of quantitative agreement in the lower panel of Fig 4B clearly depends on the precise choice of γ. Nevertheless we find this comparison very satisfying, when two wildly different approaches to a measurement of the same quantity yield broadly compatible results. We emphasize the agreement between our rates and the rates reported in [55] for any reasonable γ: values differ by at most a factor of 2 and possibly agree to within our uncertainties of 10-20%. From this we feel confident asserting that the parameters we have inferred from Jones et. al.’s single-cell mRNA counts data do in fact correspond to repressor binding and unbinding rates, and therefore our conclusions on the agreement of these rates with binding energies from [41] and [43] are valid.

3.2.4 Model checking

In Fig 3B we saw that the simple Poisson model of a constitutive promoter, despite having a well behaved posterior, was clearly insufficient to describe the data. It behooves us to carry out a similar check for our model of simple repression, codified by Eq 27 for the steady-state mRNA copy number distribution. As derived in Sections 1 and 2, we have compelling theoretical reasons to believe it is a good model, but if it nevertheless turned out to be badly contradicted by the data we should like to know.

The details are deferred to the S1 Supporting Information, and here we only attempt to summarize the intuitive ideas, as detailed at greater length by Jaynes [62] as well as Gelman and coauthors [51, 63]. From our samples of the posterior distribution, plotted in Fig 4A, we generate many replicate data using a random number generator. In Fig 4C, we plot empirical cumulative distribution functions of the middle 95% quantiles of these replicate data with the actual experimental data from Jones et. al. [33] overlaid, covering all ten experimental conditions spanning repressor binding sites and copy numbers (as well as the constitutive baseline UV5).

The purpose of Fig 4C is simply a graphical, qualitative assessment of the model: do the experimental data systematically disagree with the simulated data, which would suggest that our model is missing important features? A further question is not just whether there is a detectable difference between simulated and experimental data, but whether this difference is likely to materially affect the conclusions we draw from the posterior in Fig 4A. More rigorous and quantitative statistical tests are possible [51], but their quantitativeness does not necessarily make them more useful. As stated in [63], we often find this graphical comparison more enlightening because it better engages our intuition for the model, not merely telling if the model is wrong but suggesting how the model may be incomplete.

Our broad brush takeaway from Fig 4C is overall of good agreement. There some oddities, in particular the long tails in the data for Oid, 1 ng/mL, and O2, 0.5 ng/mL. The latter is especially odd since it extends beyond the tail of the unregulated UV5 distribution. This is a relatively small number of cells, however, so whether this is a peculiarity of the experimental data, a statistical fluke of small numbers, or a real biological effect is unclear. It is conceivable that there is some very slow timescale switching dynamics that could cause this bimodality, although it is unclear why it would only appear for specific repressor copy numbers. There is also a small offset between experiment and simulation for O2 at the higher repressor copy numbers, especially at 2 and 10 ng/mL. From the estimate of repressor copy numbers from [33], it is possible that the repressor copy numbers here are becoming large enough to partially invalidate our assumption of a separation of timescales between burst duration and repressor association rate. Another possibility is that the very large number of zero mRNA counts for Oid, 2 ng/mL is skewing its partner datasets through the shared association rate. None of these fairly minor potential caveats cause us to seriously doubt the overall correctness of our model, which further validates its use to compare the thermodynamic models’ binding energies to the kinetic models’ repressor binding and unbinding rates, as we originally set out to do.