When $k=0$, the term inside the sum on the right-hand side cancels with the
other $m^{x}$, so we can simplify to

$\mathbf{R}_{m}\mathbf{\left\langle
p^{y}\left[(m+1)^{x}-m^{x}\right]\right\rangle}=\mathbf{R}_{m}\mathbf{\left\langle
p^{y}\left[\sum_{k=1}^{x}{x\choose k}m^{x-k}\right]\right\rangle}.$ (S71)

Once the sum is expanded we can see that the highest moment in this sum is
given by $\mathbf{\left\langle m^{(x-1)}p^{y}\right\rangle}$ which satisfies
both of the conditions on Eq. S65.

For the mRNA degradation term we similarly have

$\mathbf{\Gamma}_{m}\mathbf{\left\langle
mp^{y}\left[(m-1)^{x}-m^{x}\right]\right\rangle}=\mathbf{\Gamma}_{m}\mathbf{\left\langle
mp^{y}\left[\sum_{k=0}^{x}{x\choose
k}(-1)^{k}m^{x-k}-m^{x}\right]\right\rangle}.$ (S72)

Simplifying terms we obtain

$\mathbf{\Gamma}_{m}\mathbf{\left\langle mp^{y}\left[\sum_{k=0}^{x}{x\choose
k}(-1)^{k}m^{x-k}-m^{x}\right]\right\rangle}=\mathbf{\Gamma}_{m}\mathbf{\left\langle
p^{y}\left[\sum_{k=1}^{x}{x\choose k}(-1)^{k}m^{x+1-k}\right]\right\rangle}.$
(S73)

The largest moment in this case is $\mathbf{\left\langle
m^{x}p^{y}\right\rangle}$, which again satisfies the conditions on Eq. S65.

The protein production term gives

$\mathbf{R}_{p}\mathbf{\left\langle
m^{(x+1)}\left[(p+1)^{y}-p^{y}\right]\right\rangle}=\mathbf{R}_{p}\mathbf{\left\langle
m^{(x+1)}\left[\sum_{k=0}^{y}{y\choose
k}(-1)^{k}p^{y-k}-p^{y}\right]\right\rangle}.$ (S74)

Upon simplification we obtain

$\mathbf{R}_{p}\mathbf{\left\langle m^{(x+1)}\left[\sum_{k=0}^{y}{y\choose
k}(-1)^{k}p^{y-k}-p^{y}\right]\right\rangle}=\mathbf{R}_{p}\mathbf{\left\langle
m^{(x+1)}\left[\sum_{k=1}^{y}{y\choose
k}(-1)^{k}p^{y-k}\right]\right\rangle}.$ (S75)

Here the largest moment is given by $\mathbf{\left\langle
m^{x+1}p^{y-1}\right\rangle}$, that again satisfies both of our conditions.
For the last term, for protein degradation we have

$Rp\mathbf{\left\langle
m^{(x+1)}\left[(p+1)^{y}-p^{y}\right]\right\rangle}=Rp\mathbf{\left\langle
m^{(x+1)}\left[\sum_{k=1}^{y}{y\choose
k}(-1^{k})p^{y-k}\right]\right\rangle}.$ (S76)

The largest moment involved in this term is therefore $\mathbf{\left\langle
m^{x}p^{y-1}\right\rangle}$. With this, we show that the four terms involved
in our general moment equation depend only on lower moments that satisfy Eq.
S65.

As a reminder, what we showed in this section is that the kinetic model
introduced in Fig. 2(A) has no moment-closure problem. In other words, moments
of the joint mRNA and protein distribution can be computed just from knowledge
of lower moments. This allows us to cleanly integrate the moments of the
distribution dynamics as cells progress through the cell cycle.

###  S3.3 Computing single promoter steady-state moments

(Note: The Python code used for the calculations presented in this section can
be found in the following link as an annotated Jupyter notebook)

As discussed in Section 1.4, one of the main factors contributing to cell-to-
cell variability in gene expression is the change in gene copy number during
the cell cycle as cells replicate their genome before cell division. Our
minimal model accounts for this variability by considering the time trajectory
of the distribution moments as given by Eq. S66. These predictions will be
contrasted with the predictions from a kinetic model that doesn’t account for
changes in gene copy number during the cell cycle in Appendix S5.

If we do not account for change in gene copy number during the cell cycle or
for the partition of proteins during division, the dynamics of the moments of
the distribution described in this section will reach a steady state. In order
to compute the steady-state moments of the kinetic model with a single gene
across the cell cycle, we use the moment closure property of our master
equation. By equating Eq. S66 to zero for a given $\mathbf{x}$ and
$\mathbf{y}$, we can solve the resulting linear system and obtain a solution
for $\mathbf{\left\langle m^{x}p^{y}\right\rangle}$ at steady state as a
function of moments $\mathbf{\left\langle
m^{x^{\prime}}p^{y^{\prime}}\right\rangle}$ that satisfy Eq. S65. Then, by
solving for the zero${}^{\text{th}}$ moment $\mathbf{\left\langle
m^{0}p^{0}\right\rangle}$ subject to the constraint that the probability of
the promoter being in any state should add up to one, we can substitute back
all of the solutions in terms of moments $\mathbf{\left\langle
m^{x^{\prime}}p^{y^{\prime}}\right\rangle}$ with solutions in terms of the
rates shown in Fig. 2. In other words, through an iterative process, we can
get at the value of any moment of the distribution. We start by solving for
the zero${}^{\text{th}}$ moment. Since all higher moments, depend on lower
moments we can use the solution of the zero${}^{\text{th}}$ moment to compute
the first mRNA moment. This solution is then used for higher moments in a
hierarchical iterative process.

##  S4 Accounting for the variability in gene copy number during the cell
cycle

(Note: The Python code used for the calculations presented in this section can
be found in the following link as an annotated Jupyter notebook)

When growing in rich media, bacteria can double every $\approx$ 20 minutes.
With two replication forks each traveling at $\approx$ 1000 bp per second, and
a genome of $\approx$ 5 Mbp for E. coli [61], a cell would need $\approx$ 40
minutes to replicate its genome. The apparent paradox of growth rates faster
than one division per 40 minutes is solved by the fact that cells have
multiple replisomes, i.e. molecular machines that replicate the genome running
in parallel. Cells can have up to 8 copies of the genome being replicated
simultaneously depending on the growth rate [37].

This observation implies that during the cell cycle gene copy number varies.
This variation depends on the growth rate and the relative position of the
gene with respect to the replication origin, having genes close to the
replication origin spending more time with multiple copies compare to genes
closer to the replication termination site. This change in gene dosage has a
direct effect on the cell-to-cell variability in gene expression [25, 38].

###  S4.1 Numerical integration of moment equations

(Note: The Python code used for the calculations presented in this section can
be found in the following link as an annotated Jupyter notebook)

For our specific locus (galK) and a doubling time of $\approx$ 60 min for our
experimental conditions, cells have on average 1.66 copies of the reporter
gene during the cell cycle [25]. What this means is that cells spend 60% of
the time having one copy of the gene and 40% of the time with two copies. To
account for this variability in gene copy number across the cell cycle we
numerically integrate the moment equations derived in Appendix S3 for a time
$t=[0,t_{s}]$ with an mRNA production rate $r_{m}$, where $t_{s}$ is the time
point at which the replication fork reaches our specific locus. For the
remaining time before the cell division $t=[t_{s},t_{d}]$ that the cell spends
with two promoters, we assume that the only parameter that changes is the mRNA
production rate from $r_{m}$ to $2r_{m}$. This simplifying assumption ignores
potential changes in protein translation rate $r_{p}$ or changes in the
repressor copy number that would be reflected in changes on the repressor on
rate $k^{(r)}_{\text{on}}$.

####  S4.1.1 Computing distribution moments after cell division

(Note: The Python code used for the calculations presented in this section can
be found in the following link as an annotated Jupyter notebook)

We have already solved a general form for the dynamics of the moments of the
distribution, i.e. we wrote differential equations for the moments
${d\left\langle m^{x}p^{y}\right\rangle\over dt}$. Given that we know all
parameters for our model we can simply integrate these equations numerically
to compute how the moments of the distribution evolve as cells progress
through their cell cycle. Once the cell reaches a time $t_{d}$ when is going
to divide the mRNA and proteins that we are interested in undergo a binomial
partitioning between the two daughter cells. In other words, each molecule
flips a coin and decides whether to go to either daughter. The question then
becomes given that we have a value for the moment $\left\langle
m^{x}p^{y}\right\rangle_{t_{d}}$ at a time before the cell division, what
would the value of this moment be after the cell division takes place
$\left\langle m^{x}p^{y}\right\rangle_{t_{o}}$?

The probability distribution of mRNA and protein after the cell division
$P_{t_{o}}(m,p)$ must satisfy

$P_{t_{o}}(m,p)=\sum_{m^{\prime}=m}^{\infty}\sum_{p^{\prime}=p}^{\infty}P(m,p\mid
m^{\prime},p^{\prime})P_{t_{d}}(m^{\prime},p^{\prime}),$ (S77)

where we are summing over all the possibilities of having $m^{\prime}$ mRNA
and $p^{\prime}$ proteins before cell division. Note that the sums start at
$m$ and $p$; this is because for a cell to have these copy numbers before cell
division it is a requirement that the mother cell had at least such copy
number since we are not assuming that there is any production at the
instantaneous cell division time. Since we assume that the partition of mRNA
is independent from the partition of protein, the conditional probability
$P(m,p\mid m^{\prime},p^{\prime})$ is simply given by a product of two
binomial distributions, one for the mRNA and one for the protein, i.e.

$P(m,p\mid m^{\prime},p^{\prime})={m^{\prime}\choose m}\left({1\over
2}\right)^{m^{\prime}}\cdot{p^{\prime}\choose p}\left({1\over
2}\right)^{p^{\prime}}.$ (S78)

Because of these product of binomial probabilities are allowed to extend the
sum from Eq. S118 to start at $m^{\prime}=0$ and $p^{\prime}=0$ as

$P_{t_{o}}(m,p)=\sum_{m^{\prime}=0}^{\infty}\sum_{p^{\prime}=0}^{\infty}P(m,p\mid
m^{\prime},p^{\prime})P_{t_{d}}(m^{\prime},p^{\prime}),$ (S79)

since the product of the binomial distributions in Eq. S119 is zero for all
$m^{\prime}<m$ and/or $p^{\prime}<0$. So from now on in this section we will
assume that a sum of the form $\sum_{x}\equiv\sum_{x=0}^{\infty}$ to simplify
notation.

We can then compute the distribution moments after the cell division
$\left\langle m^{x}p^{y}\right\rangle_{t_{o}}$ as

$\left\langle
m^{x}p^{y}\right\rangle_{t_{o}}=\sum_{m}\sum_{p}m^{x}p^{y}P_{t_{o}}(m,p),$
(S80)

for all $x,y\in\mathbb{N}$. Substituting Eq. S118 results in

$\left\langle
m^{x}p^{y}\right\rangle_{t_{o}}=\sum_{m}\sum_{p}m^{x}p^{y}\sum_{m^{\prime}}\sum_{p^{\prime}}P(m,p\mid
m^{\prime},p^{\prime})P_{t_{d}}(m^{\prime},p^{\prime}).$ (S81)

We can rearrange the sums to be

$\left\langle
m^{x}p^{y}\right\rangle_{t_{o}}=\sum_{m^{\prime}}\sum_{p^{\prime}}P_{t_{d}}(m^{\prime},p^{\prime})\sum_{m}\sum_{p}m^{x}p^{y}P(m,p\mid
m^{\prime},p^{\prime}).$ (S82)

The fact that Eq. S119 is the product of two independent events allows us to
rewrite the joint probability $P(m,p\mid m^{\prime},p^{\prime})$ as

$P(m,p\mid m^{\prime},p^{\prime})=P(m\mid m^{\prime})\cdot P(p\mid
p^{\prime}).$ (S83)

With this we can then write the moment $\left\langle
m^{x}p^{y}\right\rangle_{t_{o}}$ as

$\left\langle
m^{x}p^{y}\right\rangle_{t_{o}}=\sum_{m^{\prime}}\sum_{p^{\prime}}P_{t_{d}}(m^{\prime},p^{\prime})\sum_{m}m^{x}P(m\mid
m^{\prime})\sum_{p}p^{y}P(p\mid p^{\prime}).$ (S84)

Notice that both terms summing over $m$ and over $p$ are the conditional
expected values, i.e.

$\sum_{z}z^{x}P(z\mid z^{\prime})\equiv\left\langle z^{x}\mid
z^{\prime}\right\rangle,\;{\text{ for }z\in\\{m,p\\}}.$ (S85)

These conditional expected values are the expected values of a binomial random
variable $z\sim\text{Bin}(z^{\prime},1/2)$, which can be easily computed as we
will show later in this section. We then rewrite the expected values after the
cell division in terms of these moments of a binomial distribution

$\left\langle
m^{x}p^{y}\right\rangle_{t_{o}}=\sum_{m^{\prime}}\sum_{p^{\prime}}\left\langle
m^{x}\mid m^{\prime}\right\rangle\left\langle p^{y}\mid
p^{\prime}\right\rangle P_{t_{d}}(m^{\prime},p^{\prime}).$ (S86)

To see how this general formula for the moments after the cell division works
let’s compute the mean protein per cell after the cell division $\left\langle
p\right\rangle_{t_{o}}$. That is setting $x=0$, and $y=1$. This results in

$\left\langle
p\right\rangle_{t_{o}}=\sum_{m^{\prime}}\sum_{p^{\prime}}\left\langle
m^{0}\mid m^{\prime}\right\rangle\left\langle p\mid p^{\prime}\right\rangle
P_{t_{d}}(m^{\prime},p^{\prime}).$ (S87)

The zeroth moment $\left\langle m^{0}\mid m^{\prime}\right\rangle$ by
definition must be one since we have

$\left\langle m^{0}\mid m^{\prime}\right\rangle=\sum_{m}m^{0}P(m\mid
m^{\prime})=\sum_{m}P(m\mid m^{\prime})=1,$ (S88)

since the probability distribution must be normalized. This leaves us then
with

$\left\langle
p\right\rangle_{t_{o}}=\sum_{m^{\prime}}\sum_{p^{\prime}}P_{t_{d}}(m^{\prime},p^{\prime})\left\langle
p\mid p^{\prime}\right\rangle.$ (S89)

If we take the sum over $m^{\prime}$ we simply compute the marginal
probability distribution
$\sum_{m^{\prime}}P_{t_{d}}(m^{\prime},p^{\prime})=P_{t_{d}}(p^{\prime})$,
then we have

$\left\langle p\right\rangle_{t_{o}}=\sum_{p^{\prime}}\left\langle p\mid
p^{\prime}\right\rangle P_{t_{d}}(p^{\prime}).$ (S90)

For the particular case of the first moment of the binomial distribution with
parameters $p^{\prime}$ and $1/2$ we know that

$\left\langle p\mid p^{\prime}\right\rangle={p^{\prime}\over 2}.$ (S91)

Therefore the moment after division is equal to

$\left\langle p\right\rangle_{t_{o}}=\sum_{p^{\prime}}{p^{\prime}\over
2}P_{t_{d}}(p^{\prime})={1\over
2}\sum_{p^{\prime}}p^{\prime}P_{t_{d}}(p^{\prime}).$ (S92)

Notice that this is just 1/2 of the expected value of $p^{\prime}$ averaging
over the distribution prior to cell division, i.e.

$\left\langle p\right\rangle_{t_{o}}={\left\langle
p^{\prime}\right\rangle_{t_{d}}\over 2},$ (S93)

where $\left\langle\cdot\right\rangle_{t_{d}}$ highlights that is the moment
of the distribution prior to the cell division. This result makes perfect
sense. What this is saying is that the mean protein copy number right after
the cell divides is half of the mean protein copy number just before the cell
division. That is exactly we would expect. So in principle to know the first
moment of either the mRNA distribution $\left\langle m\right\rangle_{t_{o}}$
or the protein distribution $\left\langle m\right\rangle_{t_{o}}$ right after
cell division it suffices to multiply the moments before the cell division
$\left\langle m\right\rangle_{t_{d}}$ or $\left\langle p\right\rangle_{t_{d}}$
by 1/2. Let’s now explore how this generalizes to any other moment
$\left\langle m^{x}p^{y}\right\rangle_{t_{o}}$.

####  S4.1.2 Computing the moments of a binomial distribution

The result from last section was dependent on us knowing the functional form
of the first moment of the binomial distribution. For higher moments we need
some systematic way to compute such moments. Luckily for us we can do so by
using the so-called moment generating function (MGF). The MGF of a random
variable $X$ is defined as

$M_{X}(t)=\left\langle e^{tX}\right\rangle,$ (S94)

where $t$ is a dummy variable. Once we know the MGF we can obtain any moment
of the distribution by simply computing

$\left\langle X^{n}\right\rangle=\left.{d^{n}\over
dt^{n}}M_{X}(t)\right|_{t=0},$ (S95)

i.e. taking the $n$-th derivative of the MGF returns the $n$-th moment of the
distribution. For the particular case of the binomial distribution
$X\sim\text{Bin}(N,q)$ it can be shown that the MGF is of the form

$M_{X}(t)=\left[(1-q)+qe^{t}\right]^{N}.$ (S96)

As an example let’s compute the first moment of this binomially distributed
variable. For this, the first derivative of the MGF results in

${dM_{X}(t)\over dt}=N[(1-q)+qe^{t}]^{N-1}qe^{t}.$ (S97)

We just need to follow Eq. S136 and set $t=0$ to obtain the first moment

$\left.{dM_{X}(t)\over dt}\right|_{t=0}=Nq,$ (S98)

which is exactly the expected value of a binomially distributed random
variable.

So according to Eq. S127 to compute any moment $\left\langle
m^{x}p^{y}\right\rangle$ after cell division we can just take the $x$-th
derivative and the $y$-th derivative of the binomial MGF to obtain
$\left\langle m^{x}\mid m^{\prime}\right\rangle$ and $\left\langle p^{y}\mid
p^{\prime}\right\rangle$, respectively, and take the expected value of the
result. Let’s follow on detail the specific case for the moment $\left\langle
mp\right\rangle$. When computing the moment after cell division $\left\langle
mp\right\rangle_{t_{o}}$ which is of the form

$\left\langle
mp\right\rangle_{t_{o}}=\sum_{m^{\prime}}\sum{p^{\prime}}\left\langle m\mid
m^{\prime}\right\rangle\left\langle p\mid p^{\prime}\right\rangle
P_{t_{d}}(m^{\prime},p^{\prime}),$ (S99)

the product $\left\langle m\mid m^{\prime}\right\rangle\left\langle p\mid
p^{\prime}\right\rangle$ is then

$\left\langle m\mid m^{\prime}\right\rangle\left\langle p\mid
p^{\prime}\right\rangle={m^{\prime}\over 2}\cdot{p^{\prime}\over 2},$ (S100)

where we used the result in Eq. S139, substituting $m$ and $p$ for $X$,
respectively, and $q$ for 1/2. Substituting this result into the moment gives

$\left\langle
mp\right\rangle_{t_{o}}=\sum_{m^{\prime}}\sum_{p^{\prime}}{m^{\prime}p^{\prime}\over
4}P_{t_{d}}(m^{\prime},p^{\prime})={\left\langle
m^{\prime}p^{\prime}\right\rangle_{t_{d}}\over 4}.$ (S101)

Therefore to compute the moment after cell division $\left\langle
mp\right\rangle_{t_{o}}$ we simply have to divide by 4 the corresponding
equivalent moment before the cell division.

Not all moments after cell division depend only on the equivalent moment
before cell division. For example if we compute the third moment of the
protein distribution $\left\langle p^{3}\right\rangle_{t_{o}}$, we find

$\left\langle p^{3}\right\rangle_{t_{o}}={\left\langle
p^{3}\right\rangle_{t_{d}}\over 8}+{3\left\langle
p^{2}\right\rangle_{t_{d}}\over 8}.$ (S102)

So for this particular case the third moment of the protein distribution
depends on the third moment and the second moment before the cell division. In
general all moments after cell division $\left\langle
m^{x}p^{y}\right\rangle_{t_{o}}$ linearly depend on moments before cell
division. Furthermore, there is “moment closure” for this specific case in the
sense that all moments after cell division depend on lower moments before cell
division. To generalize these results to all the moments computed in this work
let us then define a vector to collect all moments before the cell division up
the $\left\langle m^{x}p^{y}\right\rangle_{t_{d}}$ moment, i.e.

$\mathbf{\left\langle m^{x}p^{y}\right\rangle}_{t_{d}}=\left(\left\langle
m^{0}p^{0}\right\rangle_{t_{d}},\left\langle
m^{1}\right\rangle_{t_{d}},\ldots,\left\langle
m^{x}p^{y}\right\rangle_{t_{d}}\right).$ (S103)

Then any moment after cell division $\left\langle
m^{x^{\prime}}p^{y^{\prime}}\right\rangle_{t_{o}}$ for $x^{\prime}\leq x$ and
$y^{\prime}\leq y$ can be computed as

$\left\langle
m^{x^{\prime}}p^{y^{\prime}}\right\rangle_{t_{o}}=\mathbf{z}_{x^{\prime}y^{\prime}}\cdot\mathbf{\left\langle
m^{x}p^{y}\right\rangle}_{t_{d}},$

where we define the vector $\mathbf{z}_{x^{\prime}y^{\prime}}$ as the vector
containing all the coefficients that we obtain with the product of the two
binomial distributions. For example for the case of the third protein moment
$\left\langle p^{3}\right\rangle_{t_{o}}$ the vector
$\mathbf{z}_{x^{\prime}y^{\prime}}$ would have zeros for all entries except
for the corresponding entry for $\left\langle p^{2}\right\rangle_{t_{d}}$ and
for $\left\langle p^{3}\right\rangle_{t_{d}}$, where it would have $3/8$ and
$1/8$ accordingly.

If we want then to compute all the moments after the cell division up to
$\left\langle m^{x}p^{y}\right\rangle_{t_{o}}$ let us define an equivalent
vector

$\mathbf{\left\langle m^{x}p^{y}\right\rangle}_{t_{o}}=\left(\left\langle
m^{0}p^{0}\right\rangle_{t_{o}},\left\langle
m^{1}\right\rangle_{t_{o}},\ldots,\left\langle
m^{x}p^{y}\right\rangle_{t_{o}}\right).$ (S104)

Then we need to build a square matrix $\mathbf{Z}$ such that each row of the
matrix contains the corresponding vector $\mathbf{z}_{x^{\prime}y^{\prime}}$
for each of the moments. Having this matrix we would simply compute the
moments after the cell division as

$\mathbf{\left\langle
m^{x}p^{x}\right\rangle}_{t_{o}}=\mathbf{Z}\cdot\mathbf{\left\langle
m^{x}p^{x}\right\rangle}_{t_{d}}.$ (S105)

In other words, matrix $\mathbf{Z}$ will contain all the coefficients that we
need to multiply by the moments before the cell division in order to obtain
the moments after cell division. Matrix $\mathbf{Z}$ was then generated
automatically using Python’s analytical math library sympy [62].

Fig. S15 (adapted from Fig. 3(B)) shows how the first moment of both mRNA and
protein changes over several cell cycles. The mRNA quickly relaxes to the
steady state corresponding to the parameters for both a single and two
promoter copies. This is expected since the parameters for the mRNA production
were determined in the first place under this assumption (See Appendix S1). We
note that there is no apparent delay before reaching steady state of the mean
mRNA count after the cell divides. This is because the mean mRNA count for the
two promoters copies state is exactly twice the expected mRNA count for the
single promoter state (See Appendix S1). Therefore once the mean mRNA count is
halved after the cell division, it is already at the steady state value for
the single promoter case. On the other hand, given that the relaxation time to
steady state is determined by the degradation rate, the mean protein count
does not reach its corresponding steady state value for either promoter copy
number state. Interestingly once a couple of cell cycles have passed the first
moment has a repetitive trajectory over cell cycles. We have observed this
experimentally by tracking cells as they grow under the microscope. Comparing
cells at the beginning of the cell cycle with the daughter cells that appear
after cell division shown that on average all cells have the same amount of
protein at the beginning of the cell cycle (See Fig. 18 of [28]), suggesting
that these dynamical steady state takes place in vivo.

Figure S8: First and second moment dynamics over cell the cell cycle. Mean
$\pm$ standard deviation mRNA (upper panel) and mean $\pm$ standard deviation
protein copy number (lower panel) as the cell cycle progresses. The dark
shaded region delimits the fraction of the cell cycle that cells spend with a
single copy of the promoter. The light shaded region delimits the fraction of
the cell cycle that cells spend with two copies of the promoter. For a 100 min
doubling time at the galK locus cells spend 60% of the time with one copy of
the promoter and the rest with two copies.

In principle when measuring gene expression levels experimentally from an
asynchronous culture, cells are sampled from any time point across their
individual cell cycles. This means that the moments determined experimentally
correspond to an average over the cell cycle. In the following section we
discuss how to account for the fact that cells are not uniformly distributed
across the cell cycle in order to compute these averages.

###  S4.2 Exponentially distributed ages

As mentioned in Appendix S2, cells in exponential growth have exponentially
distributed ages across the cell cycle, having more young cells compared to
old ones. Specifically the probability of a cell being at any time point in
the cell cycle is given by [40]

$P(a)=(\ln 2)\cdot 2^{1-a},$ (S106)

where $a\in[0,1]$ is the stage of the cell cycle, with $a=0$ being the start
of the cycle and $a=1$ being the cell division. In Appendix S10 we reproduce
this derivation. It is a surprising result, but can be intuitively thought as
follows: If the culture is growing exponentially, that means that all the time
there is an increasing number of cells. That means for example that if in a
time interval $\Delta t$ $N$ “old” cells divided, these produced $2N$ “young”
cells. So at any point there is always more younger than older cells.

Our numerical integration of the moment equations gave us a time evolution of
the moments as cells progress through the cell cycle. Since experimentally we
sample asynchronous cells that follow Eq. S147, each time point along the
moment dynamic must be weighted by the probability of having sampled a cell at
such specific time point of the cell cycle. Without loss of generality let’s
focus on the first mRNA moment $\left\langle m(t)\right\rangle$ (the same can
be applied to all other moments). As mentioned before, in order to calculate
the first moment across the entire cell cycle we must weigh each time point by
the corresponding probability that a cell is found in such point of its cell
cycle. This translates to computing the integral

$\left\langle m\right\rangle_{c}=\int_{\text{beginning cell cycle}}^{\text{end
cell cycle}}\left\langle m(t)\right\rangle P(t)dt,$ (S107)

where $\left\langle m\right\rangle_{c}$ is the mean mRNA copy number averaged
over the entire cell cycle trajectory, and $P(t)$ is the probability of a cell
being at a time $t$ of its cell cycle.

If we set the time in units of the cell cycle length we can use Eq. S147 and
compute instead

$\left\langle m\right\rangle=\int_{0}^{1}\left\langle m(a)\right\rangle
P(a)da,$ (S108)

where $P(a)$ is given by Eq. S147.

What Eq. S149 implies is that in order to compute the first moment (or any
moment of the distribution) we must weigh each point in the moment dynamics by
the corresponding probability of a cell being at that point along its cell
cycle. That is why when computing a moment we take the time trajectory of a
single cell cycle as the ones shown in Fig. S15 and compute the average using
Eq. S147 to weigh each time point. We perform this integral numerically for
all moments using Simpson’s rule.

###  S4.3 Reproducing the equilibrium picture

Given the large variability of the first moments depicted in Fig. S15 it is
worth considering why a simplistic equilibrium picture has shown to be very
successful in predicting the mean expression level under diverse conditions
[20, 26, 24, 27]. In this section we compare the simple repression
thermodynamic model with this dynamical picture of the cell cycle. But before
diving into this comparison, it is worth recapping the assumptions that go
into the equilibrium model.

####  S4.3.1 Steady state under the thermodynamic model

Given the construction of the thermodynamic model of gene regulation for which
the probability of the promoter microstates rather than the probability of
mRNA or protein counts is accounted for, we are only allowed to describe the
dynamics of the first moment using this theoretical framework [36]. Again
let’s only focus on the mRNA first moment $\left\langle m\right\rangle$. The
same principles apply if we consider the protein first moment. We can write a
dynamical system of the form

${d{\left\langle m\right\rangle}\over dt}=r_{m}\cdot
p_{{}_{\text{bound}}}-\gamma_{m}\left\langle m\right\rangle,$ (S109)

where as before $r_{m}$ and $\gamma_{m}$ are the mRNA production and
degradation rates respectively, and $p_{{}_{\text{bound}}}$ is the probability
of finding the RNAP bound to the promoter [54]. This dynamical system is
predicted to have a single stable fixed point that we can find by computing
the steady state. When we solve for the mean mRNA copy number at steady state
$\left\langle m\right\rangle_{ss}$ we find

$\left\langle
m\right\rangle_{ss}={r_{m}\over\gamma_{m}}p_{{}_{\text{bound}}}.$ (S110)

Since we assume that the only effect that the repressor has over the
regulation of the promoter is exclusion of the RNAP from binding to the
promoter, we assume that only $p_{{}_{\text{bound}}}$ depends on the repressor
copy number $R$. Therefore when computing the fold-change in gene expression
we are left with

$\text{fold-change}={\left\langle m(R\neq
0)\right\rangle_{ss}\over\left\langle
m(R=0)\right\rangle_{ss}}={p_{{}_{\text{bound}}}(R\neq 0)\over
p_{{}_{\text{bound}}}(R=0)}.$ (S111)

As derived in [20] this can be written in the language of equilibrium
statistical mechanics as

$\text{fold-change}=\left(1+{R\over
N_{\text{NS}}}e^{-\beta\Delta\varepsilon_{r}}\right)^{-1},$ (S112)

where $\beta\equiv(k_{B}T)^{-1}$, $\Delta\varepsilon_{r}$ is the repressor-DNA
binding energy, and $N_{\text{NS}}$ is the number of non-specific binding
sites where the repressor can bind.

To arrive at Eq. S153 we ignore the physiological changes that occur during
the cell cycle; one of the most important being the variability in gene copy
number that we are exploring in this section. It is therefore worth thinking
about whether or not the dynamical picture exemplified in Fig. S15 can be
reconciled with the predictions made by Eq. S153 both at the mRNA and protein
level.

Fig. S16 compares the predictions of both theoretical frameworks for varying
repressor copy numbers and repressor-DNA affinities. The solid lines are
directly computed from Eq. S153. The hollow triangles and the solid circles,
represent the fold-change in mRNA and protein respectively as computed from
the moment dynamics. To compute the fold-change from the kinetic picture we
first numerically integrate the moment dynamics for both the two- and the
three-state promoter (See Fig. S15 for the unregulated case) and then average
the time series accounting for the probability of cells being sampled at each
stage of the cell cycle as defined in Eq. S149. The small systematic
deviations between both models come partly from the simplifying assumption
that the repressor copy number, and therefore the repressor on rate
$k^{(r)}_{\text{on}}$ remains constant during the cell cycle. In principle the
gene producing the repressor protein itself is also subjected to the same
duplication during the cell cycle, changing therefore the mean repressor copy
number for both stages.

Figure S9: Comparison of the equilibrium and kinetic reressor titration
predictions. The equilibrium model (solid lines) and the kinetic model with
variation over the cell cycle (solid circles and white triangles) predictions
are compared for varying repressor copy numbers and operator binding energy.
The equilibrium model is directly computed from Eq. S153 while the kinetic
model is computed by numerically integrating the moment equations over several
cell cycles, and then averaging over the extent of the cell cycle as defined
in Eq. S149.

For completeness Fig. S17 compares the kinetic and equilibrium models for the
extended model of [27] in which the inducer concentration enters into the
equation. The solid line is directly computed from Eq. 5 of [27]. The hollow
triangles and solid points follow the same procedure as for Fig. S16, where
the only effect that the inducer is assume to have in the kinetics is an
effective change in the number of active repressors, affecting therefore
$k^{(r)}_{\text{on}}$.

Figure S10: Comparison of the equilibrium and kinetic inducer titration
predictions. The equilibrium model (solid lines) and the kinetic model with
variation over the cell cycle (solid circles and white triangles) predictions
are compared for varying repressor copy numbers and inducer concentrations.
The equilibrium model is directly computed as Eq. 5 of reference [27] with
repressor-DNA binding energy $\Delta\varepsilon_{r}=-13.5\;k_{B}T$ while the
kinetic model is computed by numerically integrating the moment dynamics over
several cell cycles, and then averaging over the extent of a single cell cycle
as defined in Eq. S149.

###  S4.4 Comparison between single- and multi-promoter kinetic model

After these calculations it is worth questioning whether the inclusion of this
change in gene dosage is drastically different with respect to the simpler
picture of a kinetic model that ignores the gene copy number variability
during the cell cycle. To this end we systematically computed the average
moments for varying repressor copy number and repressor-DNA affinities. We
then compare these results with the moments obtained from a single-promoter
model and their corresponding parameters. The derivation of the steady-state
moments of the distribution for the single-promoter model are detailed in
Appendix S3.

Fig. S16 and Fig. S17 both suggest that since the dynamic multi-promoter model
can reproduce the results of the equilibrium model at the first moment level
it must then also be able to reproduce the results of the single-promoter
model at this level (See Appendix S2). The interesting comparison comes with
higher moments. A useful metric to consider for gene expression variability is
the noise in gene expression [53]. This quantity, defined as the standard
deviation divided by the mean, is a dimensionless metric of how much
variability there is with respect to the mean of a distribution. As we will
show below this quantity differs from the also commonly used metric known as
the Fano factor (variance / mean) in the sense that for experimentally
determined expression levels in fluorescent arbitrary units, the noise is a
dimensionless quantity while the Fano factor is not.

Fig. S18 shows the comparison of the predicted protein noise between the
single- (dashed lines) and the multi-promoter model (solid lines) for
different operators and repressor copy numbers. A striking difference between
both is that the single-promoter model predicts that as the inducer
concentration increases, the standard deviation grows much slower than the
mean, giving a very small noise. In comparison the multi-promoter model has a
much higher floor for the lowest value of the noise, reflecting the expected
result that the variability in gene copy number across the cell cycle should
increase the cell-to-cell variability in gene expression [38, 25]

Figure S11: Comparison of the predicted protein noise between a single- and a
multi-promoter kinetic model. Comparison of the noise (standard
deviation/mean) between a kinetic model that considers a single promoter at
all times (dashed line) and the multi-promoter model developed in this section
(solid line) for different repressor operators. (A) Operator O1,
$\Delta\varepsilon_{r}=-15.3\;k_{B}T$, (B) O2,
$\Delta\varepsilon_{r}=-13.9\;k_{B}T$, (C) O3,
$\Delta\varepsilon_{r}=-9.7\;k_{B}T$

###  S4.5 Comparison with experimental data

Having shown that the kinetic model presented in this section can not only
reproduce the results from the equilibrium picture at the mean level (See Fig.
S16 and Fig. S17), but make predictions for the cell-to-cell variability as
quantified by the noise (See Fig. S18), we can assess whether or not this
model is able to predict experimental measurements of the noise. For this we
take the single cell intensity measurements (See Methods) to compute the noise
at the protein level.

As mentioned before this metric differs from the Fano factor since for
fluorescent arbitrary units the noise is a dimensionless quantity. To see why
consider that the noise is defined as

$\text{noise}\equiv\frac{\sqrt{\left\langle p^{2}\right\rangle-\left\langle
p\right\rangle^{2}}}{\left\langle p\right\rangle}.$ (S113)

We assume that the intensity level of a cell $I$ is linearly proportional to
the absolute protein count, i.e.

$I=\alpha p,$ (S114)

where $\alpha$ is the proportionality constant between arbitrary units and
protein absolute number $p$. Substituting this definition on Eq. S154 gives

$\text{noise}={\sqrt{\left\langle(\alpha
I)^{2}\right\rangle-\left\langle\alpha
I\right\rangle^{2}}\over\left\langle\alpha I\right\rangle}.$ (S115)

Since $\alpha$ is a constant it can be taken out of the average operator
$\left\langle\cdot\right\rangle$, obtaining

$\text{noise}={\sqrt{\alpha^{2}\left(\left\langle
I^{2}\right\rangle-\left\langle
I\right\rangle^{2}\right)}\over\alpha\left\langle
I\right\rangle}={\sqrt{\left(\left\langle I^{2}\right\rangle-\left\langle
I\right\rangle^{2}\right)}\over\left\langle I\right\rangle}.$ (S116)

Notice that in Eq. S155 the linear proportionality between intensity and
protein count has no intercept. This ignores the autofluorescence that cells
without reporter would generate. To account for this, in practice we compute

$\text{noise}={\sqrt{\left(\left\langle(I-\left\langle
I_{\text{auto}}\right\rangle)^{2}\right\rangle-\left\langle I-\left\langle
I_{\text{auto}}\right\rangle\right\rangle^{2}\right)}\over\left\langle
I-\left\langle I_{\text{auto}}\right\rangle\right\rangle}.$ (S117)

where $I$ is the intensity of the strain of interest and $\left\langle
I_{\text{auto}}\right\rangle$ is the mean autofluorescence intensity, obtained
from a strain that does not carry the fluorescent reporter gene.

Fig. S19 shows the comparison between theoretical predictions and experimental
measurements for the unregulated promoter. The reason we split the data by
operator despite the fact that since these are unregulated promoters, they
should in principle have identical expression profiles is to precisely make
sure that this is the case. We have found in the past that sequences
downstream of the RNAP binding site can affect the expression level of
constitutively expressed genes. We can see that both models, the single-
promoter (gray dotted line) and the multi-promoter (black dashed line)
underestimate the experimental noise to different degrees. The single-promoter
model does a worse job at predicting the experimental data since it doesn’t
account for the differences in gene dosage during the cell cycle. But still we
can see that accounting for this variability takes us to within a factor of
two of the experimentally determined noise for these unregulated strains.

Figure S12: Protein noise of the unregulated promoter. Comparison of the
experimental noise for different operators with the theoretical predictions
for the single-promoter (gray dotted line) and the multi-promoter model (black
dashed line). Each datum represents a single date measurement of the
corresponding $\Delta lacI$ strain with $\geq 300$ cells. The points
correspond to the median, and the error bars correspond to the 95% confidence
interval as determined by 10,000 bootstrap samples.

To further test the model predictive power we compare the predictions for the
three-state regulated promoter. Fig. S20 shows the theoretical predictions for
the single- and multi-promoter model for varying repressor copy numbers and
repressor-DNA binding affinities as a function of the inducer concentration.
We can see again that our zero-parameter fits systematically underestimates
the noise for all strains and all inducer concentrations. We highlight that
the $y$-axis is shown in a log-scale to emphasize more this deviation; but, as
we will show in the next section, our predictions still fall within a factor
of two from the experimental data.

Figure S13: Protein noise of the regulated promoter. Comparison of the
experimental noise for different operators ((A) O1,
$\Delta\varepsilon_{r}=-15.3\;k_{B}T$, (B) O2,
$\Delta\varepsilon_{r}=-13.9\;k_{B}T$, (C) O3,
$\Delta\varepsilon_{r}=-9.7\;k_{B}T$) with the theoretical predictions for the
single-promoter (dashed lines) and the multi-promoter model (solid lines).
Points represent the experimental noise as computed from single-cell
fluorescence measurements of different E. coli strains under 12 different
inducer concentrations. Dotted line indicates plot in linear rather than
logarithmic scale. Each datum represents a single date measurement of the
corresponding strain and IPTG concentration with $\geq 300$ cells. The points
correspond to the median, and the error bars correspond to the 95% confidence
interval as determined by 10,000 bootstrap samples. White-filled dots are plot
at a different scale for better visualization.

####  S4.5.1 Systematic deviation of the noise in gene expression

Fig. S19 and Fig. S20 highlight that our model underestimates the cell-to-cell
variability as measured by the noise. To further explore this systematic
deviation Fig. S21 shows the theoretical vs. experimental noise both in linear
and log scale. As we can see the data is systematically above the identity
line. The data is colored by their corresponding experimental fold-change
values. The data that has the largest deviations from the identity line also
corresponds to the data with the largest error bars and the smallest fold-
change. This is because measurements with very small fold-changes correspond
to intensities very close to the autofluorescence background. Therefore
minimal changes when computing the noise are amplified given the ratio of
std/mean. In Appendix S9 we will explore empirical ways to improve the
agreement between our minimal model and the experimental data to guide future
efforts to improve the minimal.

Figure S14: Systematic comparison of theoretical vs experimental noise in gene
expression. Theoretical vs. experimental noise both in linear (left) and log
(right) scale. The dashed line shows the identity line of slope 1 and
intercept zero. All data are colored by the corresponding value of the
experimental fold-change in gene expression as indicated by the color bar.
Each datum represents a single date measurement of the corresponding strain
and IPTG concentration with $\geq 300$ cells. The points correspond to the
median, and the error bars correspond to the 95% confidence interval as
determined by 10,000 bootstrap samples.

##  S5 Accounting for the variability in gene copy number during the cell
cycle

(Note: The Python code used for the calculations presented in this section can
be found in the following link as an annotated Jupyter notebook)

When growing in rich media, bacteria can double every $\approx$ 20 minutes.
With two replication forks each traveling at $\approx$ 1000 bp per second, and
a genome of $\approx$ 5 Mbp for E. coli [61], a cell would need $\approx$ 40
minutes to replicate its genome. The apparent paradox of growth rates faster
than one division per 40 minutes is solved by the fact that cells have
multiple replisomes, i.e. molecular machines that replicate the genome running
in parallel. Cells can have up to 8 copies of the genome being replicated
simultaneously depending on the growth rate [37].

This observation implies that during the cell cycle gene copy number varies.
This variation depends on the growth rate and the relative position of the
gene with respect to the replication origin, having genes close to the
replication origin spending more time with multiple copies compare to genes
closer to the replication termination site. This change in gene dosage has a
direct effect on the cell-to-cell variability in gene expression [25, 38].

###  S5.1 Numerical integration of moment equations

(Note: The Python code used for the calculations presented in this section can
be found in the following link as an annotated Jupyter notebook)

For our specific locus (galK) and a doubling time of $\approx$ 60 min for our
experimental conditions, cells have on average 1.66 copies of the reporter
gene during the cell cycle [25]. What this means is that cells spend 60% of
the time having one copy of the gene and 40% of the time with two copies. To
account for this variability in gene copy number across the cell cycle we
numerically integrate the moment equations derived in Appendix S3 for a time
$t=[0,t_{s}]$ with an mRNA production rate $r_{m}$, where $t_{s}$ is the time
point at which the replication fork reaches our specific locus. For the
remaining time before the cell division $t=[t_{s},t_{d}]$ that the cell spends
with two promoters, we assume that the only parameter that changes is the mRNA
production rate from $r_{m}$ to $2r_{m}$. This simplifying assumption ignores
potential changes in protein translation rate $r_{p}$ or changes in the
repressor copy number that would be reflected in changes on the repressor on
rate $k^{(r)}_{\text{on}}$.

####  S5.1.1 Computing distribution moments after cell division

(Note: The Python code used for the calculations presented in this section can
be found in the following link as an annotated Jupyter notebook)

We have already solved a general form for the dynamics of the moments of the
distribution, i.e. we wrote differential equations for the moments
${d\left\langle m^{x}p^{y}\right\rangle\over dt}$. Given that we know all
parameters for our model we can simply integrate these equations numerically
to compute how the moments of the distribution evolve as cells progress
through their cell cycle. Once the cell reaches a time $t_{d}$ when is going
to divide the mRNA and proteins that we are interested in undergo a binomial
partitioning between the two daughter cells. In other words, each molecule
flips a coin and decides whether to go to either daughter. The question then
becomes given that we have a value for the moment $\left\langle
m^{x}p^{y}\right\rangle_{t_{d}}$ at a time before the cell division, what
would the value of this moment be after the cell division takes place
$\left\langle m^{x}p^{y}\right\rangle_{t_{o}}$?

The probability distribution of mRNA and protein after the cell division
$P_{t_{o}}(m,p)$ must satisfy

$P_{t_{o}}(m,p)=\sum_{m^{\prime}=m}^{\infty}\sum_{p^{\prime}=p}^{\infty}P(m,p\mid
m^{\prime},p^{\prime})P_{t_{d}}(m^{\prime},p^{\prime}),$ (S118)

where we are summing over all the possibilities of having $m^{\prime}$ mRNA
and $p^{\prime}$ proteins before cell division. Note that the sums start at
$m$ and $p$; this is because for a cell to have these copy numbers before cell
division it is a requirement that the mother cell had at least such copy
number since we are not assuming that there is any production at the
instantaneous cell division time. Since we assume that the partition of mRNA
is independent from the partition of protein, the conditional probability
$P(m,p\mid m^{\prime},p^{\prime})$ is simply given by a product of two
binomial distributions, one for the mRNA and one for the protein, i.e.

$P(m,p\mid m^{\prime},p^{\prime})={m^{\prime}\choose m}\left({1\over
2}\right)^{m^{\prime}}\cdot{p^{\prime}\choose p}\left({1\over
2}\right)^{p^{\prime}}.$ (S119)

Because of these product of binomial probabilities are allowed to extend the
sum from Eq. S118 to start at $m^{\prime}=0$ and $p^{\prime}=0$ as

$P_{t_{o}}(m,p)=\sum_{m^{\prime}=0}^{\infty}\sum_{p^{\prime}=0}^{\infty}P(m,p\mid
m^{\prime},p^{\prime})P_{t_{d}}(m^{\prime},p^{\prime}),$ (S120)

since the product of the binomial distributions in Eq. S119 is zero for all
$m^{\prime}<m$ and/or $p^{\prime}<0$. So from now on in this section we will
assume that a sum of the form $\sum_{x}\equiv\sum_{x=0}^{\infty}$ to simplify
notation.

We can then compute the distribution moments after the cell division
$\left\langle m^{x}p^{y}\right\rangle_{t_{o}}$ as

$\left\langle
m^{x}p^{y}\right\rangle_{t_{o}}=\sum_{m}\sum_{p}m^{x}p^{y}P_{t_{o}}(m,p),$
(S121)

for all $x,y\in\mathbb{N}$. Substituting Eq. S118 results in

$\left\langle
m^{x}p^{y}\right\rangle_{t_{o}}=\sum_{m}\sum_{p}m^{x}p^{y}\sum_{m^{\prime}}\sum_{p^{\prime}}P(m,p\mid
m^{\prime},p^{\prime})P_{t_{d}}(m^{\prime},p^{\prime}).$ (S122)

We can rearrange the sums to be

$\left\langle
m^{x}p^{y}\right\rangle_{t_{o}}=\sum_{m^{\prime}}\sum_{p^{\prime}}P_{t_{d}}(m^{\prime},p^{\prime})\sum_{m}\sum_{p}m^{x}p^{y}P(m,p\mid
m^{\prime},p^{\prime}).$ (S123)

The fact that Eq. S119 is the product of two independent events allows us to
rewrite the joint probability $P(m,p\mid m^{\prime},p^{\prime})$ as

$P(m,p\mid m^{\prime},p^{\prime})=P(m\mid m^{\prime})\cdot P(p\mid
p^{\prime}).$ (S124)

With this we can then write the moment $\left\langle
m^{x}p^{y}\right\rangle_{t_{o}}$ as

$\left\langle
m^{x}p^{y}\right\rangle_{t_{o}}=\sum_{m^{\prime}}\sum_{p^{\prime}}P_{t_{d}}(m^{\prime},p^{\prime})\sum_{m}m^{x}P(m\mid
m^{\prime})\sum_{p}p^{y}P(p\mid p^{\prime}).$ (S125)

Notice that both terms summing over $m$ and over $p$ are the conditional
expected values, i.e.

$\sum_{z}z^{x}P(z\mid z^{\prime})\equiv\left\langle z^{x}\mid
z^{\prime}\right\rangle,\;{\text{ for }z\in\\{m,p\\}}.$ (S126)

These conditional expected values are the expected values of a binomial random
variable $z\sim\text{Bin}(z^{\prime},1/2)$, which can be easily computed as we
will show later in this section. We then rewrite the expected values after the
cell division in terms of these moments of a binomial distribution

$\left\langle
m^{x}p^{y}\right\rangle_{t_{o}}=\sum_{m^{\prime}}\sum_{p^{\prime}}\left\langle
m^{x}\mid m^{\prime}\right\rangle\left\langle p^{y}\mid
p^{\prime}\right\rangle P_{t_{d}}(m^{\prime},p^{\prime}).$ (S127)

To see how this general formula for the moments after the cell division works
let’s compute the mean protein per cell after the cell division $\left\langle
p\right\rangle_{t_{o}}$. That is setting $x=0$, and $y=1$. This results in

$\left\langle
p\right\rangle_{t_{o}}=\sum_{m^{\prime}}\sum_{p^{\prime}}\left\langle
m^{0}\mid m^{\prime}\right\rangle\left\langle p\mid p^{\prime}\right\rangle
P_{t_{d}}(m^{\prime},p^{\prime}).$ (S128)

The zeroth moment $\left\langle m^{0}\mid m^{\prime}\right\rangle$ by
definition must be one since we have

$\left\langle m^{0}\mid m^{\prime}\right\rangle=\sum_{m}m^{0}P(m\mid
m^{\prime})=\sum_{m}P(m\mid m^{\prime})=1,$ (S129)

since the probability distribution must be normalized. This leaves us then
with

$\left\langle
p\right\rangle_{t_{o}}=\sum_{m^{\prime}}\sum_{p^{\prime}}P_{t_{d}}(m^{\prime},p^{\prime})\left\langle
p\mid p^{\prime}\right\rangle.$ (S130)

If we take the sum over $m^{\prime}$ we simply compute the marginal
probability distribution
$\sum_{m^{\prime}}P_{t_{d}}(m^{\prime},p^{\prime})=P_{t_{d}}(p^{\prime})$,
then we have

$\left\langle p\right\rangle_{t_{o}}=\sum_{p^{\prime}}\left\langle p\mid
p^{\prime}\right\rangle P_{t_{d}}(p^{\prime}).$ (S131)

For the particular case of the first moment of the binomial distribution with
parameters $p^{\prime}$ and $1/2$ we know that

$\left\langle p\mid p^{\prime}\right\rangle={p^{\prime}\over 2}.$ (S132)

Therefore the moment after division is equal to

$\left\langle p\right\rangle_{t_{o}}=\sum_{p^{\prime}}{p^{\prime}\over
2}P_{t_{d}}(p^{\prime})={1\over
2}\sum_{p^{\prime}}p^{\prime}P_{t_{d}}(p^{\prime}).$ (S133)

Notice that this is just 1/2 of the expected value of $p^{\prime}$ averaging
over the distribution prior to cell division, i.e.

$\left\langle p\right\rangle_{t_{o}}={\left\langle
p^{\prime}\right\rangle_{t_{d}}\over 2},$ (S134)

where $\left\langle\cdot\right\rangle_{t_{d}}$ highlights that is the moment
of the distribution prior to the cell division. This result makes perfect
sense. What this is saying is that the mean protein copy number right after
the cell divides is half of the mean protein copy number just before the cell
division. That is exactly we would expect. So in principle to know the first
moment of either the mRNA distribution $\left\langle m\right\rangle_{t_{o}}$
or the protein distribution $\left\langle m\right\rangle_{t_{o}}$ right after
cell division it suffices to multiply the moments before the cell division
$\left\langle m\right\rangle_{t_{d}}$ or $\left\langle p\right\rangle_{t_{d}}$
by 1/2. Let’s now explore how this generalizes to any other moment
$\left\langle m^{x}p^{y}\right\rangle_{t_{o}}$.

####  S5.1.2 Computing the moments of a binomial distribution

The result from last section was dependent on us knowing the functional form
of the first moment of the binomial distribution. For higher moments we need
some systematic way to compute such moments. Luckily for us we can do so by
using the so-called moment generating function (MGF). The MGF of a random
variable $X$ is defined as

$M_{X}(t)=\left\langle e^{tX}\right\rangle,$ (S135)

where $t$ is a dummy variable. Once we know the MGF we can obtain any moment
of the distribution by simply computing

$\left\langle X^{n}\right\rangle=\left.{d^{n}\over
dt^{n}}M_{X}(t)\right|_{t=0},$ (S136)

i.e. taking the $n$-th derivative of the MGF returns the $n$-th moment of the
distribution. For the particular case of the binomial distribution
$X\sim\text{Bin}(N,q)$ it can be shown that the MGF is of the form

$M_{X}(t)=\left[(1-q)+qe^{t}\right]^{N}.$ (S137)

As an example let’s compute the first moment of this binomially distributed
variable. For this, the first derivative of the MGF results in

${dM_{X}(t)\over dt}=N[(1-q)+qe^{t}]^{N-1}qe^{t}.$ (S138)

We just need to follow Eq. S136 and set $t=0$ to obtain the first moment

$\left.{dM_{X}(t)\over dt}\right|_{t=0}=Nq,$ (S139)

which is exactly the expected value of a binomially distributed random
variable.

So according to Eq. S127 to compute any moment $\left\langle
m^{x}p^{y}\right\rangle$ after cell division we can just take the $x$-th
derivative and the $y$-th derivative of the binomial MGF to obtain
$\left\langle m^{x}\mid m^{\prime}\right\rangle$ and $\left\langle p^{y}\mid
p^{\prime}\right\rangle$, respectively, and take the expected value of the
result. Let’s follow on detail the specific case for the moment $\left\langle
mp\right\rangle$. When computing the moment after cell division $\left\langle
mp\right\rangle_{t_{o}}$ which is of the form

$\left\langle
mp\right\rangle_{t_{o}}=\sum_{m^{\prime}}\sum{p^{\prime}}\left\langle m\mid
m^{\prime}\right\rangle\left\langle p\mid p^{\prime}\right\rangle
P_{t_{d}}(m^{\prime},p^{\prime}),$ (S140)

the product $\left\langle m\mid m^{\prime}\right\rangle\left\langle p\mid
p^{\prime}\right\rangle$ is then

$\left\langle m\mid m^{\prime}\right\rangle\left\langle p\mid
p^{\prime}\right\rangle={m^{\prime}\over 2}\cdot{p^{\prime}\over 2},$ (S141)

where we used the result in Eq. S139, substituting $m$ and $p$ for $X$,
respectively, and $q$ for 1/2. Substituting this result into the moment gives

$\left\langle
mp\right\rangle_{t_{o}}=\sum_{m^{\prime}}\sum_{p^{\prime}}{m^{\prime}p^{\prime}\over
4}P_{t_{d}}(m^{\prime},p^{\prime})={\left\langle
m^{\prime}p^{\prime}\right\rangle_{t_{d}}\over 4}.$ (S142)

Therefore to compute the moment after cell division $\left\langle
mp\right\rangle_{t_{o}}$ we simply have to divide by 4 the corresponding
equivalent moment before the cell division.

Not all moments after cell division depend only on the equivalent moment
before cell division. For example if we compute the third moment of the
protein distribution $\left\langle p^{3}\right\rangle_{t_{o}}$, we find

$\left\langle p^{3}\right\rangle_{t_{o}}={\left\langle
p^{3}\right\rangle_{t_{d}}\over 8}+{3\left\langle
p^{2}\right\rangle_{t_{d}}\over 8}.$ (S143)

So for this particular case the third moment of the protein distribution
depends on the third moment and the second moment before the cell division. In
general all moments after cell division $\left\langle
m^{x}p^{y}\right\rangle_{t_{o}}$ linearly depend on moments before cell
division. Furthermore, there is “moment closure” for this specific case in the
sense that all moments after cell division depend on lower moments before cell
division. To generalize these results to all the moments computed in this work
let us then define a vector to collect all moments before the cell division up
the $\left\langle m^{x}p^{y}\right\rangle_{t_{d}}$ moment, i.e.

$\mathbf{\left\langle m^{x}p^{y}\right\rangle}_{t_{d}}=\left(\left\langle
m^{0}p^{0}\right\rangle_{t_{d}},\left\langle
m^{1}\right\rangle_{t_{d}},\ldots,\left\langle
m^{x}p^{y}\right\rangle_{t_{d}}\right).$ (S144)

Then any moment after cell division $\left\langle
m^{x^{\prime}}p^{y^{\prime}}\right\rangle_{t_{o}}$ for $x^{\prime}\leq x$ and
$y^{\prime}\leq y$ can be computed as

$\left\langle
m^{x^{\prime}}p^{y^{\prime}}\right\rangle_{t_{o}}=\mathbf{z}_{x^{\prime}y^{\prime}}\cdot\mathbf{\left\langle
m^{x}p^{y}\right\rangle}_{t_{d}},$

where we define the vector $\mathbf{z}_{x^{\prime}y^{\prime}}$ as the vector
containing all the coefficients that we obtain with the product of the two
binomial distributions. For example for the case of the third protein moment
$\left\langle p^{3}\right\rangle_{t_{o}}$ the vector
$\mathbf{z}_{x^{\prime}y^{\prime}}$ would have zeros for all entries except
for the corresponding entry for $\left\langle p^{2}\right\rangle_{t_{d}}$ and
for $\left\langle p^{3}\right\rangle_{t_{d}}$, where it would have $3/8$ and
$1/8$ accordingly.

If we want then to compute all the moments after the cell division up to
$\left\langle m^{x}p^{y}\right\rangle_{t_{o}}$ let us define an equivalent
vector

$\mathbf{\left\langle m^{x}p^{y}\right\rangle}_{t_{o}}=\left(\left\langle
m^{0}p^{0}\right\rangle_{t_{o}},\left\langle
m^{1}\right\rangle_{t_{o}},\ldots,\left\langle
m^{x}p^{y}\right\rangle_{t_{o}}\right).$ (S145)

Then we need to build a square matrix $\mathbf{Z}$ such that each row of the
matrix contains the corresponding vector $\mathbf{z}_{x^{\prime}y^{\prime}}$
for each of the moments. Having this matrix we would simply compute the
moments after the cell division as

$\mathbf{\left\langle
m^{x}p^{x}\right\rangle}_{t_{o}}=\mathbf{Z}\cdot\mathbf{\left\langle
m^{x}p^{x}\right\rangle}_{t_{d}}.$ (S146)

In other words, matrix $\mathbf{Z}$ will contain all the coefficients that we
need to multiply by the moments before the cell division in order to obtain
the moments after cell division. Matrix $\mathbf{Z}$ was then generated
automatically using Python’s analytical math library sympy [62].

Fig. S15 (adapted from Fig. 3(B)) shows how the first moment of both mRNA and
protein changes over several cell cycles. The mRNA quickly relaxes to the
steady state corresponding to the parameters for both a single and two
promoter copies. This is expected since the parameters for the mRNA production
were determined in the first place under this assumption (See Appendix S1). We
note that there is no apparent delay before reaching steady state of the mean
mRNA count after the cell divides. This is because the mean mRNA count for the
two promoters copies state is exactly twice the expected mRNA count for the
single promoter state (See Appendix S1). Therefore once the mean mRNA count is
halved after the cell division, it is already at the steady state value for
the single promoter case. On the other hand, given that the relaxation time to
steady state is determined by the degradation rate, the mean protein count
does not reach its corresponding steady state value for either promoter copy
number state. Interestingly once a couple of cell cycles have passed the first
moment has a repetitive trajectory over cell cycles. We have observed this
experimentally by tracking cells as they grow under the microscope. Comparing
cells at the beginning of the cell cycle with the daughter cells that appear
after cell division shown that on average all cells have the same amount of
protein at the beginning of the cell cycle (See Fig. 18 of [28]), suggesting
that these dynamical steady state takes place in vivo.

Figure S15: First and second moment dynamics over cell the cell cycle. Mean
$\pm$ standard deviation mRNA (upper panel) and mean $\pm$ standard deviation
protein copy number (lower panel) as the cell cycle progresses. The dark
shaded region delimits the fraction of the cell cycle that cells spend with a
single copy of the promoter. The light shaded region delimits the fraction of
the cell cycle that cells spend with two copies of the promoter. For a 100 min
doubling time at the galK locus cells spend 60% of the time with one copy of
the promoter and the rest with two copies.

In principle when measuring gene expression levels experimentally from an
asynchronous culture, cells are sampled from any time point across their
individual cell cycles. This means that the moments determined experimentally
correspond to an average over the cell cycle. In the following section we
discuss how to account for the fact that cells are not uniformly distributed
across the cell cycle in order to compute these averages.

###  S5.2 Exponentially distributed ages

As mentioned in Appendix S2, cells in exponential growth have exponentially
distributed ages across the cell cycle, having more young cells compared to
old ones. Specifically the probability of a cell being at any time point in
the cell cycle is given by [40]

$P(a)=(\ln 2)\cdot 2^{1-a},$ (S147)

where $a\in[0,1]$ is the stage of the cell cycle, with $a=0$ being the start
of the cycle and $a=1$ being the cell division. In Appendix S10 we reproduce
this derivation. It is a surprising result, but can be intuitively thought as
follows: If the culture is growing exponentially, that means that all the time
there is an increasing number of cells. That means for example that if in a
time interval $\Delta t$ $N$ “old” cells divided, these produced $2N$ “young”
cells. So at any point there is always more younger than older cells.

Our numerical integration of the moment equations gave us a time evolution of
the moments as cells progress through the cell cycle. Since experimentally we
sample asynchronous cells that follow Eq. S147, each time point along the
moment dynamic must be weighted by the probability of having sampled a cell at
such specific time point of the cell cycle. Without loss of generality let’s
focus on the first mRNA moment $\left\langle m(t)\right\rangle$ (the same can
be applied to all other moments). As mentioned before, in order to calculate
the first moment across the entire cell cycle we must weigh each time point by
the corresponding probability that a cell is found in such point of its cell
cycle. This translates to computing the integral

$\left\langle m\right\rangle_{c}=\int_{\text{beginning cell cycle}}^{\text{end
cell cycle}}\left\langle m(t)\right\rangle P(t)dt,$ (S148)

where $\left\langle m\right\rangle_{c}$ is the mean mRNA copy number averaged
over the entire cell cycle trajectory, and $P(t)$ is the probability of a cell
being at a time $t$ of its cell cycle.

If we set the time in units of the cell cycle length we can use Eq. S147 and
compute instead

$\left\langle m\right\rangle=\int_{0}^{1}\left\langle m(a)\right\rangle
P(a)da,$ (S149)

where $P(a)$ is given by Eq. S147.

What Eq. S149 implies is that in order to compute the first moment (or any
moment of the distribution) we must weigh each point in the moment dynamics by
the corresponding probability of a cell being at that point along its cell
cycle. That is why when computing a moment we take the time trajectory of a
single cell cycle as the ones shown in Fig. S15 and compute the average using
Eq. S147 to weigh each time point. We perform this integral numerically for
all moments using Simpson’s rule.

###  S5.3 Reproducing the equilibrium picture

Given the large variability of the first moments depicted in Fig. S15 it is
worth considering why a simplistic equilibrium picture has shown to be very
successful in predicting the mean expression level under diverse conditions
[20, 26, 24, 27]. In this section we compare the simple repression
thermodynamic model with this dynamical picture of the cell cycle. But before
diving into this comparison, it is worth recapping the assumptions that go
into the equilibrium model.

####  S5.3.1 Steady state under the thermodynamic model

Given the construction of the thermodynamic model of gene regulation for which
the probability of the promoter microstates rather than the probability of
mRNA or protein counts is accounted for, we are only allowed to describe the
dynamics of the first moment using this theoretical framework [36]. Again
let’s only focus on the mRNA first moment $\left\langle m\right\rangle$. The
same principles apply if we consider the protein first moment. We can write a
dynamical system of the form

${d{\left\langle m\right\rangle}\over dt}=r_{m}\cdot
p_{{}_{\text{bound}}}-\gamma_{m}\left\langle m\right\rangle,$ (S150)

where as before $r_{m}$ and $\gamma_{m}$ are the mRNA production and
degradation rates respectively, and $p_{{}_{\text{bound}}}$ is the probability
of finding the RNAP bound to the promoter [54]. This dynamical system is
predicted to have a single stable fixed point that we can find by computing
the steady state. When we solve for the mean mRNA copy number at steady state
$\left\langle m\right\rangle_{ss}$ we find

$\left\langle
m\right\rangle_{ss}={r_{m}\over\gamma_{m}}p_{{}_{\text{bound}}}.$ (S151)

Since we assume that the only effect that the repressor has over the
regulation of the promoter is exclusion of the RNAP from binding to the
promoter, we assume that only $p_{{}_{\text{bound}}}$ depends on the repressor
copy number $R$. Therefore when computing the fold-change in gene expression
we are left with

$\text{fold-change}={\left\langle m(R\neq
0)\right\rangle_{ss}\over\left\langle
m(R=0)\right\rangle_{ss}}={p_{{}_{\text{bound}}}(R\neq 0)\over
p_{{}_{\text{bound}}}(R=0)}.$ (S152)

As derived in [20] this can be written in the language of equilibrium
statistical mechanics as

$\text{fold-change}=\left(1+{R\over
N_{\text{NS}}}e^{-\beta\Delta\varepsilon_{r}}\right)^{-1},$ (S153)

where $\beta\equiv(k_{B}T)^{-1}$, $\Delta\varepsilon_{r}$ is the repressor-DNA
binding energy, and $N_{\text{NS}}$ is the number of non-specific binding
sites where the repressor can bind.

To arrive at Eq. S153 we ignore the physiological changes that occur during
the cell cycle; one of the most important being the variability in gene copy
number that we are exploring in this section. It is therefore worth thinking
about whether or not the dynamical picture exemplified in Fig. S15 can be
reconciled with the predictions made by Eq. S153 both at the mRNA and protein
level.

Fig. S16 compares the predictions of both theoretical frameworks for varying
repressor copy numbers and repressor-DNA affinities. The solid lines are
directly computed from Eq. S153. The hollow triangles and the solid circles,
represent the fold-change in mRNA and protein respectively as computed from
the moment dynamics. To compute the fold-change from the kinetic picture we
first numerically integrate the moment dynamics for both the two- and the
three-state promoter (See Fig. S15 for the unregulated case) and then average
the time series accounting for the probability of cells being sampled at each
stage of the cell cycle as defined in Eq. S149. The small systematic
deviations between both models come partly from the simplifying assumption
that the repressor copy number, and therefore the repressor on rate
$k^{(r)}_{\text{on}}$ remains constant during the cell cycle. In principle the
gene producing the repressor protein itself is also subjected to the same
duplication during the cell cycle, changing therefore the mean repressor copy
number for both stages.

Figure S16: Comparison of the equilibrium and kinetic reressor titration
predictions. The equilibrium model (solid lines) and the kinetic model with
variation over the cell cycle (solid circles and white triangles) predictions
are compared for varying repressor copy numbers and operator binding energy.
The equilibrium model is directly computed from Eq. S153 while the kinetic
model is computed by numerically integrating the moment equations over several
cell cycles, and then averaging over the extent of the cell cycle as defined
in Eq. S149.

For completeness Fig. S17 compares the kinetic and equilibrium models for the
extended model of [27] in which the inducer concentration enters into the
equation. The solid line is directly computed from Eq. 5 of [27]. The hollow
triangles and solid points follow the same procedure as for Fig. S16, where
the only effect that the inducer is assume to have in the kinetics is an
effective change in the number of active repressors, affecting therefore
$k^{(r)}_{\text{on}}$.

Figure S17: Comparison of the equilibrium and kinetic inducer titration
predictions. The equilibrium model (solid lines) and the kinetic model with
variation over the cell cycle (solid circles and white triangles) predictions
are compared for varying repressor copy numbers and inducer concentrations.
The equilibrium model is directly computed as Eq. 5 of reference [27] with
repressor-DNA binding energy $\Delta\varepsilon_{r}=-13.5\;k_{B}T$ while the
kinetic model is computed by numerically integrating the moment dynamics over
several cell cycles, and then averaging over the extent of a single cell cycle
as defined in Eq. S149.

###  S5.4 Comparison between single- and multi-promoter kinetic model

After these calculations it is worth questioning whether the inclusion of this
change in gene dosage is drastically different with respect to the simpler
picture of a kinetic model that ignores the gene copy number variability
during the cell cycle. To this end we systematically computed the average
moments for varying repressor copy number and repressor-DNA affinities. We
then compare these results with the moments obtained from a single-promoter
model and their corresponding parameters. The derivation of the steady-state
moments of the distribution for the single-promoter model are detailed in
Appendix S3.

Fig. S16 and Fig. S17 both suggest that since the dynamic multi-promoter model
can reproduce the results of the equilibrium model at the first moment level
it must then also be able to reproduce the results of the single-promoter
model at this level (See Appendix S2). The interesting comparison comes with
higher moments. A useful metric to consider for gene expression variability is
the noise in gene expression [53]. This quantity, defined as the standard
deviation divided by the mean, is a dimensionless metric of how much
variability there is with respect to the mean of a distribution. As we will
show below this quantity differs from the also commonly used metric known as
the Fano factor (variance / mean) in the sense that for experimentally
determined expression levels in fluorescent arbitrary units, the noise is a
dimensionless quantity while the Fano factor is not.

Fig. S18 shows the comparison of the predicted protein noise between the
single- (dashed lines) and the multi-promoter model (solid lines) for
different operators and repressor copy numbers. A striking difference between
both is that the single-promoter model predicts that as the inducer
concentration increases, the standard deviation grows much slower than the
mean, giving a very small noise. In comparison the multi-promoter model has a
much higher floor for the lowest value of the noise, reflecting the expected
result that the variability in gene copy number across the cell cycle should
increase the cell-to-cell variability in gene expression [38, 25]

Figure S18: Comparison of the predicted protein noise between a single- and a
multi-promoter kinetic model. Comparison of the noise (standard
deviation/mean) between a kinetic model that considers a single promoter at
all times (dashed line) and the multi-promoter model developed in this section
(solid line) for different repressor operators. (A) Operator O1,
$\Delta\varepsilon_{r}=-15.3\;k_{B}T$, (B) O2,
$\Delta\varepsilon_{r}=-13.9\;k_{B}T$, (C) O3,
$\Delta\varepsilon_{r}=-9.7\;k_{B}T$

###  S5.5 Comparison with experimental data

Having shown that the kinetic model presented in this section can not only
reproduce the results from the equilibrium picture at the mean level (See Fig.
S16 and Fig. S17), but make predictions for the cell-to-cell variability as
quantified by the noise (See Fig. S18), we can assess whether or not this
model is able to predict experimental measurements of the noise. For this we
take the single cell intensity measurements (See Methods) to compute the noise
at the protein level.

As mentioned before this metric differs from the Fano factor since for
fluorescent arbitrary units the noise is a dimensionless quantity. To see why
consider that the noise is defined as

$\text{noise}\equiv\frac{\sqrt{\left\langle p^{2}\right\rangle-\left\langle
p\right\rangle^{2}}}{\left\langle p\right\rangle}.$ (S154)

We assume that the intensity level of a cell $I$ is linearly proportional to
the absolute protein count, i.e.

$I=\alpha p,$ (S155)

where $\alpha$ is the proportionality constant between arbitrary units and
protein absolute number $p$. Substituting this definition on Eq. S154 gives

$\text{noise}={\sqrt{\left\langle(\alpha
I)^{2}\right\rangle-\left\langle\alpha
I\right\rangle^{2}}\over\left\langle\alpha I\right\rangle}.$ (S156)

Since $\alpha$ is a constant it can be taken out of the average operator
$\left\langle\cdot\right\rangle$, obtaining

$\text{noise}={\sqrt{\alpha^{2}\left(\left\langle
I^{2}\right\rangle-\left\langle
I\right\rangle^{2}\right)}\over\alpha\left\langle
I\right\rangle}={\sqrt{\left(\left\langle I^{2}\right\rangle-\left\langle
I\right\rangle^{2}\right)}\over\left\langle I\right\rangle}.$ (S157)

Notice that in Eq. S155 the linear proportionality between intensity and
protein count has no intercept. This ignores the autofluorescence that cells
without reporter would generate. To account for this, in practice we compute

$\text{noise}={\sqrt{\left(\left\langle(I-\left\langle
I_{\text{auto}}\right\rangle)^{2}\right\rangle-\left\langle I-\left\langle
I_{\text{auto}}\right\rangle\right\rangle^{2}\right)}\over\left\langle
I-\left\langle I_{\text{auto}}\right\rangle\right\rangle}.$ (S158)

where $I$ is the intensity of the strain of interest and $\left\langle
I_{\text{auto}}\right\rangle$ is the mean autofluorescence intensity, obtained
from a strain that does not carry the fluorescent reporter gene.

Fig. S19 shows the comparison between theoretical predictions and experimental
measurements for the unregulated promoter. The reason we split the data by
operator despite the fact that since these are unregulated promoters, they
should in principle have identical expression profiles is to precisely make
sure that this is the case. We have found in the past that sequences
downstream of the RNAP binding site can affect the expression level of
constitutively expressed genes. We can see that both models, the single-
promoter (gray dotted line) and the multi-promoter (black dashed line)
underestimate the experimental noise to different degrees. The single-promoter
model does a worse job at predicting the experimental data since it doesn’t
account for the differences in gene dosage during the cell cycle. But still we
can see that accounting for this variability takes us to within a factor of
two of the experimentally determined noise for these unregulated strains.

Figure S19: Protein noise of the unregulated promoter. Comparison of the
experimental noise for different operators with the theoretical predictions
for the single-promoter (gray dotted line) and the multi-promoter model (black
dashed line). Each datum represents a single date measurement of the
corresponding $\Delta lacI$ strain with $\geq 300$ cells. The points
correspond to the median, and the error bars correspond to the 95% confidence
interval as determined by 10,000 bootstrap samples.

To further test the model predictive power we compare the predictions for the
three-state regulated promoter. Fig. S20 shows the theoretical predictions for
the single- and multi-promoter model for varying repressor copy numbers and
repressor-DNA binding affinities as a function of the inducer concentration.
We can see again that our zero-parameter fits systematically underestimates
the noise for all strains and all inducer concentrations. We highlight that
the $y$-axis is shown in a log-scale to emphasize more this deviation; but, as
we will show in the next section, our predictions still fall within a factor
of two from the experimental data.

Figure S20: Protein noise of the regulated promoter. Comparison of the
experimental noise for different operators ((A) O1,
$\Delta\varepsilon_{r}=-15.3\;k_{B}T$, (B) O2,
$\Delta\varepsilon_{r}=-13.9\;k_{B}T$, (C) O3,
$\Delta\varepsilon_{r}=-9.7\;k_{B}T$) with the theoretical predictions for the
single-promoter (dashed lines) and the multi-promoter model (solid lines).
Points represent the experimental noise as computed from single-cell
fluorescence measurements of different E. coli strains under 12 different
inducer concentrations. Dotted line indicates plot in linear rather than
logarithmic scale. Each datum represents a single date measurement of the
corresponding strain and IPTG concentration with $\geq 300$ cells. The points
correspond to the median, and the error bars correspond to the 95% confidence
interval as determined by 10,000 bootstrap samples. White-filled dots are plot
at a different scale for better visualization.

####  S5.5.1 Systematic deviation of the noise in gene expression

Fig. S19 and Fig. S20 highlight that our model underestimates the cell-to-cell
variability as measured by the noise. To further explore this systematic
deviation Fig. S21 shows the theoretical vs. experimental noise both in linear
and log scale. As we can see the data is systematically above the identity
line. The data is colored by their corresponding experimental fold-change
values. The data that has the largest deviations from the identity line also
corresponds to the data with the largest error bars and the smallest fold-
change. This is because measurements with very small fold-changes correspond
to intensities very close to the autofluorescence background. Therefore
minimal changes when computing the noise are amplified given the ratio of
std/mean. In Appendix S9 we will explore empirical ways to improve the
agreement between our minimal model and the experimental data to guide future
efforts to improve the minimal.

Figure S21: Systematic comparison of theoretical vs experimental noise in gene
expression. Theoretical vs. experimental noise both in linear (left) and log
(right) scale. The dashed line shows the identity line of slope 1 and
intercept zero. All data are colored by the corresponding value of the
experimental fold-change in gene expression as indicated by the color bar.
Each datum represents a single date measurement of the corresponding strain
and IPTG concentration with $\geq 300$ cells. The points correspond to the
median, and the error bars correspond to the 95% confidence interval as
determined by 10,000 bootstrap samples.

##  S6 Maximum entropy approximation of distributions

(Note: The Python code used for the calculations presented in this section can
be found in the following link as an annotated Jupyter notebook)

On the one hand the solution of chemical master equations like the one in
Section 1.1 represent a hard mathematical challenge. As presented in Appendix
S2 Peccoud and Ycart derived a closed-form solution for the two-state promoter
at the mRNA level [34]. In an impressive display of mathematical skills,
Shahrezaei and Swain were able to derive an approximate solution for the one-
(not considered in this work) and two-state promoter master equation at the
protein level [53]. Nevertheless both of these solutions do not give
instantaneous insights about the distributions as they involve complicated
terms such as confluent hypergeometric functions.

On the other hand there has been a great deal of work to generate methods that
can approximate the solution of these discrete state Markovian models [63, 64,
65, 66, 67]. In particular for master equations like the one that concerns us
here whose moments can be easily computed, the moment expansion method
provides a simple method to approximate the full joint distribution of mRNA
and protein [67]. In this section we will explain the principles behind this
method and show the implementation for our particular case study.

###  S6.1 The MaxEnt principle

The principle of maximum entropy (MaxEnt) first proposed by E. T. Jaynes in
1957 tackles the question of given limited information what is the least
biased inference one can make about a particular probability distribution
[42]. In particular Jaynes used this principle to show the correspondence
between statistical mechanics and information theory, demonstrating, for
example, that the Boltzmann distribution is the probability distribution that
maximizes Shannon’s entropy subject to a constraint that the average energy of
the system is fixed.

To illustrate the principle let us focus on a univariate distribution
$P_{X}(x)$. The $n^{\text{th}}$ moment of the distribution for a discrete set
of possible values of $x$ is given by

$\left\langle x^{n}\right\rangle\equiv\sum_{x}x^{n}P_{X}(x).$ (S159)

Now assume that we have knowledge of the first $m$ moments
$\mathbf{\left\langle x\right\rangle}_{m}=(\left\langle
x\right\rangle,\left\langle x^{2}\right\rangle,\ldots,\left\langle
x^{m}\right\rangle)$. The question is then how can we use this information to
build an estimator $P_{H}(x\mid\mathbf{\left\langle x\right\rangle}_{m})$ of
the distribution such that

$\lim_{m\rightarrow\infty}P_{H}(x\mid\mathbf{\left\langle
x\right\rangle}_{m})\rightarrow P_{X}(x),$ (S160)

i.e. that the more moments we add to our approximation, the more the estimator
distribution converges to the real distribution.

The MaxEnt principle tells us that our best guess for this estimator is to
build it on the base of maximizing the Shannon entropy, constrained by the
information we have about these $m$ moments. The maximization of Shannon’s
entropy guarantees that we are the least committed possible to information
that we do not posses. The Shannon entropy for an univariate discrete
distribution is given by [3]

$H(x)\equiv-\sum_{x}P_{X}(x)\log P_{X}(x).$ (S161)

For an optimization problem subject to constraints we make use of the method
of the Lagrange multipliers. For this we define the constraint equation
$\mathcal{L}(x)$ as

$\mathcal{L}(x)\equiv H(x)-\sum_{i=0}^{m}\left[\lambda_{i}\left(\left\langle
x^{i}\right\rangle-\sum_{x}x^{i}P_{X}(x)\right)\right],$ (S162)

where $\lambda_{i}$ is the Lagrange multiplier associated with the
$i^{\text{th}}$ moment. The inclusion of the zeroth moment is an additional
constraint to guarantee the normalization of the resulting distribution. Since
$P_{X}(x)$ has a finite set of discrete values, when taking the derivative of
the constraint equation with respect to $P_{X}(x)$, we chose a particular
value of $X=x$. Therefore from the sum over all possible $x$ values only a
single term survives. With this in mind we take the derivative of the
constraint equation obtaining

${d\mathcal{L}\over dP_{X}(x)}=-\log
P_{X}(x)-1-\sum_{i=0}^{m}\lambda_{i}x^{i}.$ (S163)

Equating this derivative to zero and solving for the distribution (that we now
start calling $P_{H}(x)$, our MaxEnt estimator) gives

$P_{H}(x)=\exp\left(-1-\sum_{i=0}^{m}\lambda_{i}x^{i}\right)={1\over\mathcal{Z}}\exp\left(-\sum_{i=1}^{m}\lambda_{i}x^{i}\right),$
(S164)

where $\mathcal{Z}$ is the normalization constant that can be obtained by
substituting this solution into the normalization constraint. This results in

$\mathcal{Z}\equiv\exp\left(1+\lambda_{0}\right)=\sum_{x}\exp\left(-\sum_{i=1}^{m}\lambda_{i}x^{i}\right).$
(S165)

Eq. S164 is the general form of the MaxEnt distribution for a univariate
distribution. The computational challenge then consists in finding numerical
values for the Lagrange multipliers $\\{\lambda_{i}\\}$ such that $P_{H}(x)$
satisfies our constraints. In other words, the Lagrange multipliers weight the
contribution of each term in the exponent such that when computing any of the
moments we recover the value of our constraint. Mathematically what this means
is that $P_{H}(x)$ must satisfy

$\sum_{x}x^{n}P_{H}(x)=\sum_{x}{x^{n}\over\mathcal{Z}}\exp\left(-\sum_{i=1}^{m}\lambda_{i}x^{i}\right)=\left\langle
x^{n}\right\rangle.$ (S166)

As an example of how to apply the MaxEnt principle let us use the classic
problem of a six-face die. If we are only told that after a large number of
die rolls the mean value of the face is $\left\langle x\right\rangle=4.5$
(note that a fair die has a mean of $3.5$), what would the least biased guess
for the distribution look like? The MaxEnt principle tells us that our best
guess would be of the form

$P_{H}(x)={1\over\mathcal{Z}}\exp\left(\lambda x\right).$ (S167)

Using any numerical minimization package we can easily find the value of the
Lagrange multiplier $\lambda$ that satisfies our constraint. Fig. S22 shows
two two examples of distributions that satisfy the constraint. Panel (A) shows
a distribution consistent with the 4.5 average where both 4 and 5 are equally
likely. Nevertheless in the information we got about the nature of the die it
was never stated that some of the faces were forbidden. In that sense the
distribution is committing to information about the process that we do not
posses. Panel (B) by contrast shows the MaxEnt distribution that satisfies
this constraint. Since this distribution maximizes Shannon’s entropy it is
guaranteed to be the least biased distribution given the available
information.

Figure S22: Maximum entropy distribution of six-face die. (A)biased
distribution consistent with the constraint $\left\langle x\right\rangle=4.5$.
(B) MaxEnt distribution also consistent with the constraint.

####  S6.1.1 The mRNA and protein joint distribution

The MaxEnt principle can easily be extended to multivariate distributions. For
our particular case we are interested in the mRNA and protein joint
distribution $P(m,p)$. The definition of a moment $\left\langle
m^{x}p^{y}\right\rangle$ is a natural extension of Eq. S159 of the form

$\left\langle m^{x}p^{y}\right\rangle=\sum_{m}\sum_{p}m^{x}p^{y}P(m,p).$
(S168)

As a consequence the MaxEnt joint distribution $P_{H}(m,p)$ is of the form

$P_{H}(m,p)={1\over\mathcal{Z}}\exp\left(-\sum_{(x,y)}\lambda_{(x,y)}m^{x}p^{y}\right),$
(S169)

where $\lambda_{(x,y)}$ is the Lagrange multiplier associated with the moment
$\left\langle m^{x}p^{y}\right\rangle$, and again $\mathcal{Z}$ is the
normalization constant given by

$\mathcal{Z}=\sum_{m}\sum_{p}\exp\left(-\sum_{(x,y)}\lambda_{(x,y)}m^{x}p^{y}\right).$
(S170)

Note that the sum in the exponent is taken over all available $(x,y)$ pairs
that define the moment constraints for the distribution.

###  S6.2 The Bretthorst rescaling algorithm

The determination of the Lagrange multipliers suffer from a numerical under
and overflow problem due to the difference in magnitude between the
constraints. This becomes a problem when higher moments are taken into
account. The resulting numerical values for the Lagrange multipliers end up
being separated by several orders of magnitude. For routines such as Newton-
Raphson or other minimization algorithms that can be used to find these
Lagrange multipliers these different scales become problematic.

To get around this problem we implemented a variation to the algorithm due to
G. Larry Bretthorst, E.T. Jaynes’ last student. With a very simple argument we
can show that linearly rescaling the constraints, the Lagrange multipliers and
the “rules” for how to compute each of the moments, i.e. each of the
individual products that go into the moment calculation, should converge to
the same MaxEnt distribution. In order to see this let’s consider again a
univariate distribution $P_{X}(x)$ that we are trying to reconstruct given the
first two moments $\left\langle x\right\rangle$, and $\left\langle
x^{2}\right\rangle$. The MaxEnt distribution can be written as

$P_{H}(x)={1\over\mathcal{Z}}\exp\left(-\lambda_{1}x-\lambda_{2}x^{2}\right)={1\over\mathcal{Z}}\exp\left(-\lambda_{1}x\right)\exp\left(-\lambda_{2}x^{2}\right).$
(S171)

We can always rescale the terms in any way and obtain the same result. Let’s
say that for some reason we want to rescale the quadratic terms by a factor
$a$. We can define a new Lagrange multiplier
$\lambda_{2}^{\prime}\equiv{\lambda_{2}\over a}$ that compensates for the
rescaling of the terms, obtaining

$P_{H}(x)={1\over\mathcal{Z}}\exp\left(-\lambda_{1}x\right)\exp\left(-\lambda_{2}^{\prime}ax^{2}\right).$
(S172)

Computationally it might be more efficient to find the numerical value of
$\lambda_{2}^{\prime}$ rather than $\lambda_{2}$ maybe because it is of the
same order of magnitude as $\lambda_{1}$. Then we can always multiply
$\lambda_{2}^{\prime}$ by $a$ to obtain back the constraint for our quadratic
term. What this means is that that we can always rescale the MaxEnt problem to
make it numerically more stable, then we can rescale back to obtain the value
of the Lagrange multipliers. The key to the Bretthorst algorithm lies in the
selection of what rescaling factor to choose in order to make the numerical
inference more efficient.

Bretthorst’s algorithm goes even further by further transforming the
constraints and the variables to make the constraints orthogonal, making the
computation much more effective. We now explain the implementation of the
algorithm for our joint distribution of interest $P(m,p)$.

####  S6.2.1 Algorithm implementation

Let the $M\times N$ matrix $\mathbf{A}$ contain all the factors used to
compute the moments that serve as constraints, where each entry is of the form

$A_{ij}=m_{i}^{x_{j}}\cdot p_{i}^{y_{j}}.$ (S173)

In other words, recall that to obtain any moment $\left\langle
m^{x}p^{y}\right\rangle$ we compute

$\left\langle m^{x}p^{y}\right\rangle=\sum_{m}\sum_{p}m^{x}p^{y}P(m,x).$
(S174)

If we have $M$ possible $(m,p)$ pairs in our truncated sample space (because
we can’t include the sample space up to infinity)
$\\{(m,p)_{1},(m,p)_{2},\ldots(m,p)_{N}\\}$, and we have $N$ exponent pairs
$(x,y)$ corresponding to the $N$ moments used to constraint the maximum
entropy distribution $\\{(x,y)_{1},(x,y)_{2},\ldots,(x,y)_{N}\\}$, then matrix
$\mathbf{A}$ contains all the possible $M$ by $N$ terms of the form described
in Eq. S173. Let also $\mathbf{v}$ be a vector of length $N$ containing all
the constraints with each entry of the form

$v_{j}=\left\langle m^{x_{j}}p^{y_{j}}\right\rangle,$ (S175)

i.e. the information that we have about the distribution. That means that the
constraint equation $\mathcal{L}$ to be used for this problem takes the form

$\mathcal{L}=-\sum_{i}P_{i}\ln
P_{i}+\lambda_{0}\left(1-\sum_{i}P_{i}\right)+\sum_{j>0}\lambda_{j}\left(v_{j}-\sum_{i}A_{ij}P_{i}\right),$
(S176)

where $\lambda_{0}$ is the Lagrange multiplier associated with the
normalization constraint, and $\lambda_{j}$ is the Lagrange multiplier
associated with the $j^{\text{th}}$ constraint. This constraint equation is
equivalent to Eq. S162, but now all the details of how to compute the moments
are specified in matrix $\mathbf{A}$.

With this notation in hand we now proceed to rescale the problem. The first
step consists of rescaling the terms to compute the entries of matrix
$\mathbf{A}$. As mentioned before, this is the key feature of the Bretthorst
algorithm; the particular choice of rescaling factor used in the algorithm
empirically promotes that the rescaled Lagrange multipliers are of the same
order of magnitude. The rescaling takes the form

$A_{ij}^{\prime}={A_{ij}\over G_{j}},$ (S177)

where $G_{j}$ serves to rescale the moments, providing numerical stability to
the inference problem. Bretthorst proposes an empirical rescaling that
satisfies

$G_{j}^{2}=\sum_{i}A_{ij}^{2},$ (S178)

or in terms of our particular problem

$G_{j}^{2}=\sum_{m}\sum_{p}\left(m^{x_{j}}p^{y_{j}}\right)^{2}.$ (S179)

What this indicates is that each pair $m_{i}^{x_{j}}p_{i}^{y_{j}}$ is
normalized by the square root of the sum of the all pairs of the same form
squared.

Since we rescale the factors involved in computing the constraints, the
constraints must also be rescaled simply as

$v_{j}^{\prime}=\left\langle
m^{x_{j}}p^{y_{j}}\right\rangle^{\prime}={\left\langle
m^{x_{j}}p^{y_{j}}\right\rangle\over G_{j}}.$ (S180)

The Lagrange multipliers must compensate this rescaling since at the end of
the day the probability must add up to the same value. Therefore we rescale
the $\lambda_{j}$ terms as

$\lambda_{j}^{\prime}=\lambda_{j}G_{j},$ (S181)

such that any $\lambda_{j}A_{ij}=\lambda_{j}^{\prime}A_{ij}^{\prime}$. If this
empirical value for the rescaling factor makes the rescaled Lagrange
multipliers $\lambda_{j}^{\prime}$ be of the same order of magnitude, this by
itself would already improve the algorithm convergence. Bretthorst proposes
another linear transformation to make the optimization routine even more
efficient. For this we generate orthogonal constraints that make Newton-
Raphson and similar algorithms converge faster. The transformation is as
follows

$A_{ik}^{\prime\prime}=\sum_{j}{e}_{jk}A_{ij}^{\prime},$ (S182)

for the entires of matrix $\mathbf{A}$, and

$v_{k}^{\prime\prime}=\sum_{j}{e}_{jk}u_{j}^{\prime},$ (S183)

for entires of the constraint vector $\mathbf{v}$, finally

$\lambda_{k}^{\prime\prime}=\sum_{j}{e}_{jk}\beta_{j},$ (S184)

for the Lagrange multipliers. Here ${e}_{jk}$ is the $j^{\text{th}}$ component
of the $k^{\text{th}}$ eigenvector of the matrix $\mathbf{E}$ with entries

${E}_{kj}=\sum_{i}{A}_{ik}^{\prime}{A}_{ij}^{\prime}.$ (S185)

This transformation guarantees that the matrix $\mathbf{A}^{\prime\prime}$ has
the property

$\sum_{i}A_{ij}^{\prime\prime}A_{jk}^{\prime\prime}=\beta_{j}\delta_{jk},$
(S186)

where $\beta_{j}$ is the $j^{\text{th}}$ eigenvalue of the matrix $\mathbf{E}$
and $\delta_{jk}$ is the Kronecker delta function. What this means is that, as
desired, the constraints are orthogonal to each other, improving the algorithm
convergence speed.

###  S6.3 Predicting distributions for simple repression constructs

Having explained the theoretical background along with the practical
difficulties and a workaround strategy proposed by Bretthorst, we implemented
the inference using the moments obtained from averaging over the variability
along the cell cycle (See Appendix S5). Fig. S23 and Fig. S24 present these
inferences for both mRNA and protein levels respectively for different values
of the repressor-DNA binding energy and repressor copy numbers per cell. From
these plots we can easily appreciate that despite the fact that the mean of
each distribution changes as the induction level changes, there is a lot of
overlap between distributions. This as a consequence means that at the single-
cell level cells cannot perfectly resolve between different inputs.

Figure S23: Maximum entropy mRNA distributions for simple repression
constructs. mRNA distributions for different biophysical parameters. From left
to right the repressor-DNA affinity decreases as defined by the three lacI
operators O1 ($-15.3\;k_{B}T$), O2 ($-13.9\;k_{B}T$), and O3 ($-9.7\;k_{B}T$).
From top to bottom the mean repressor copy number per cell increases. The
curves on each plot represent different IPTG concentrations. Each distribution
was fitted using the first three moments of the mRNA distribution. Figure S24:
Maximum entropy protein distributions for simple repression constructs.
Protein distributions for different biophysical parameters. From left to right
the repressor-DNA affinity decreases as defined by the three lacI operators O1
($-15.3\;k_{B}T$), O2 ($-13.9\;k_{B}T$), and O3 ($-9.7\;k_{B}T$). From top to
bottom the mean repressor copy number per cell increases. The curves on each
plot represent different IPTG concentrations. Each distribution was fitted
using the first six moments of the protein distribution.

###  S6.4 Comparison with experimental data

Now that we have reconstructed an approximation of the probability
distribution $P(m,p)$ we can compare this with our experimental measurements.
But just as detailed in Appendix S5.5 the single-cell microscopy measurements
are given in arbitrary units of fluorescence. Therefore we cannot compare
directly our predicted protein distributions with these values. To get around
this issue we use the fact that the fold-change in gene expression that we
defined as the ratio of the gene expression level in the presence of the
repressor and the expression level of a knockout strain is a non-dimensional
quantity. Therefore we normalize all of our single-cell measurements by the
mean fluorescence value of the $\Delta lacI$ strain with the proper background
fluorescence subtracted as explained in Appendix S5.5 for the noise
measurements. In the case of the theoretical predictions of the protein
distribution we also normalize each protein value by the predicted mean
protein level $\left\langle p\right\rangle$, having now non-dimensional scales
that can be directly compared. Fig. S25 shows the experimental (color curves)
and theoretical (dark dashed line) cumulative distribution functions for the
three $\Delta lacI$ strains. As in Fig. S19, we do not expect differences
between the operators, but we explicitly plot them separately to make sure
that this is the case. We can see right away that as we would expect given the
limitations of the model to accurately predict the noise and skewness of the
distribution, the model doesn’t accurately predict the data. Our model
predicts a narrower distribution compared to what we measured with single-cell
microscopy.

Figure S25: Experiment vs. theory comparison for $\Delta lacI$ strain. Example
fold-change empirical cumulative distribution functions (ECDF) for strains
with no repressors and different operators. The color curves represent single-
cell microscopy measurements while the dashed black lines represent the
theoretical distributions as reconstructed by the maximum entropy principle.
The theoretical distributions were fitted using the first six moments of the
protein distribution.

The same narrower prediction applies to the regulated promoters. Fig. S26,
shows the theory-experiment comparison of the cumulative distribution
functions for different repressor binding sites (different figures), repressor
copy numbers (rows), and inducer concentrations (columns). In general the
predictions are systematically narrower compared to the actual experimental
data.

Figure S26: Experiment vs. theory comparison for regulated promoters. Example
fold-change empirical cumulative distribution functions (ECDF) for regulated
strains with the three operators (different colors) as a function of repressor
copy numbers (rows) and inducer concentrations (columns). The color curves
represent single-cell microscopy measurements while the dashed black lines
represent the theoretical distributions as reconstructed by the maximum
entropy principle. The theoretical distributions were fitted using the first
six moments of the protein distribution

##  S7 Gillespie simulation of master equation

(Note: The Python code used for the calculations presented in this section can
be found in the following link as an annotated Jupyter notebook)

So far we have generated a way to compute an approximated form of the joint
distribution of protein and mRNA $P(m,p)$ as a function of the moments of the
distribution $\left\langle m^{x}p^{y}\right\rangle$. This is a non-
conventional form to work with the resulting distribution of the master
equation. A more conventional approach to work with master equations whose
closed-form solutions are not known or not computable is to use stochastic
simulations commonly known as Gillespie simulations. To benchmark the
performance of our approach based on distribution moments and maximum entropy
we implemented the Gillespie algorithm. Our implementation as detailed in the
corresponding Jupyter notebook makes use of just-in-time compilation as
implemented with the Python package numba.

###  S7.1 mRNA distribution with Gillespie simulations

To confirm that the implementation of the Gillespie simulation was correct we
perform the simulation at the mRNA level for which the closed-form solution of
the steady-state distribution is known as detailed in Appendix S2. Fig. S27
shows example trajectories of mRNA counts. Each of these trajectories were
computed over several cell cyles, where the cell division was implemented
generating a binomially distributed random variable that depended on the last
mRNA count before the division event.

Figure S27: Stochastic trajectories of mRNA counts. 100 stochastic
trajectories generated with the Gillespie algorithm for mRNA counts over time
for a two-state unregulated promoter. Cells spend a fraction of the cell cycle
with a single copy of the promoter (light brown) and the rest of the cell
cycle with two copies (light yellow). When trajectories reach a new cell
cycle, the mRNA counts undergo a binomial partitioning to simulate the cell
division.

To check the implementation of our stochastic algorithm we generated several
of these stochastic trajectories in order to reconstruct the mRNA steady-state
distribution. These reconstructed distributions for a single- and double-copy
of the promoter can be compared with Eq. S10 \- the steady-state distribution
for the two-state promoter. Fig. S28 shows the great agreement between the
stochastic simulation and the analytical result, confirming that our
implementation of the Gillespie simulation is correct.

Figure S28: Comparison of analytical and simulated mRNA distribution. Solid
lines show the steady-state mRNA distributions for one copy (light blue) and
two copies of the promoter (dark blue) as defined by Eq. S10. Shaded regions
represent the corresponding distribution obtained using 2500 stochastic mRNA
trajectories and taking the last cell-cyle to approximate the distribution.

###  S7.2 Protein distribution with Gillespie simulations

Having confirmed that our implementation of the Gillespie algorithm that
includes the binomial partitioning of molecules reproduces analytical results
we extended the implementation to include protein counts. Fig. S29 shows
representative trajectories for both mRNA and protein counts over several cell
cycles. Specially for the protein we can see that it takes several cell cycles
for counts to converge to the dynamical steady-state observed with the
deterministic moment equations. Once this steady-state is reached, the
ensemble of trajectories between cell cycles look very similar.

Figure S29: Stochastic trajectories of mRNA and protein counts. 2500 protein
counts over time for a two-state unregulated promoter. Cells spend a fraction
of the cell cycle with a single copy of the promoter (light brown) and the
rest of the cell cycle with two copies (light yellow). When trajectories reach
a new cell cycle, the molecule counts undergo a binomial partitioning to
simulate the cell division.

From these trajectories we can compute the protein steady-state distribution,
taking into account the cell-age distribution as detailed in Appendix S6. Fig.
S30 shows the comparison between this distribution and the one generated using
the maximum entropy algorithm. Despite the notorious differences between the
distributions, the Gillespie simulation and the maximum entropy results are
indistinguishable in terms of the mean, variance, and skewness of the
distribution. We remind the reader that the maximum entropy is an
approximation of the distribution that gets better the more moments we add. We
therefore claim that the approximation works sufficiently well for our
purpose. The enormous advantage of the maximum entropy approach comes from the
computation time. for the number of distributions that were needed for our
calculations the Gillespie algorithm proved to be a very inefficient method
given the large sample space. Our maximum entropy approach reduces the
computation time by several orders of magnitude, allowing us to extensively
explore different parameters of the regulatory model.

Figure S30: Comparison of protein distributions. Comparison of the protein
distribution generated with Gillespie stochastic simulations (blue curve) and
the maximum entropy approach presented in Appendix S6 (orange curve). The
upper panel shows the probability mass function. The lower panel compares the
cumulative distribution functions.

##  S8 Computational determination of the channel capacity

(Note: The Python code used for the calculations presented in this section can
be found in the following link as an annotated Jupyter notebook)

In this section we detail the computation of the channel capacity of the
simple genetic circuit shown in Fig. 5. As detailed in Section 1.6 the channel
capacity is defined as the mutual information between input $c$ and output $p$
maximized over all possible input distributions $P(c)$ [3]. In principle there
are an infinite number of input distributions, so the task of finding
$\hat{P}(c)$, the input distribution at channel capacity, requires an
algorithmic approach that guarantees the convergence to this distribution.
Tkačik, Callan and Bialek developed a clever analytical approximation to find
the $\hat{P}(c)$ distribution [17]. The validity of their so-called small
noise approximation requires the standard deviation of the output distribution
$P(p\mid c)$ to be much smaller than the domain of the distribution. For our
particular case such condition is not satisfied given the spread of the
inferred protein distributions shown in Fig. 4.

Fortunately there exists a numerical algorithm to approximate $\hat{P}(c)$ for
discrete distributions. In 1972 Richard Blahut and Suguru Arimoto
independently came up with an algorithm mathematically shown to converge to
$\hat{P}(c)$ [44]. To compute both the theoretical and the experimental
channel capacity shown in Fig. 5, we implemented Blahut’s algorithm. In the
following section we detail the definitions needed for the algorithm. Then we
detail how to compute the experimental channel capacity when the bins of the
distribution are not clear given the intrinsic arbitrary nature of microscopy
fluorescence measurements.

###  S8.1 Blahut’s algorithm

Following [44] we implemented the algorithm to compute the channel capacity.
We define $\mathbf{p_{c}}$ to be an array containing the probability of each
of the input inducer concentrations (twelve concentrations, See Methods). Each
entry $j$ of the array is then of the form

$p_{c}^{(j)}=P(c=c_{j}),$ (S187)

with $j\in\\{1,2,\ldots,12\\}$. The objective of the algorithm is to find the
entries $p_{c}^{(j)}$ that maximize the mutual information between inputs and
outputs. We also define $\mathbf{Q}$ to be a $|\mathbf{p_{c}}|$ by
$|\mathbf{p_{p\mid c}}|$ matrix, where $|\cdot|$ specifies the length of the
array, and $\mathbf{p_{p\mid c}}$ is an array containing the probability
distribution of an output given a specific value of the input. In other words,
the matrix $\mathbf{Q}$ recollects all of the individual output distribution
arrays $\mathbf{p_{p\mid c}}$ into a single object. Then each entry of the
matrix $\mathbf{Q}$ is of the form

$Q^{(i,j)}=P(p=p_{i}\mid c=c_{j}).$ (S188)

For the case of the theoretical predictions of the channel capacity (Solid
lines in Fig. 5) the entries of matrix $\mathbf{Q}$ are given by the inferred
maximum entropy distributions as shown in Fig. 4. In the next section we will
discuss how to define this matrix for the case of the single-cell fluorescence
measurements. Having defined these matrices we proceed to implement the
algorithm shown in Figure 1 of [44].

###  S8.2 Channel capacity from arbitrary units of fluorescence

A difficulty when computing the channel capacity between inputs and outputs
from experimental data is that ideally we would like to compute

$C(g;c)\equiv\sup_{P(c)}I(g;c),$ (S189)

where $g$ is the gene expression level, and $c$ is the inducer concentration.
But in reality we are computing

$C(f(g);c)\equiv\sup_{P(c)}I(f(g);c),$ (S190)

where $f(g)$ is a function of gene expression that has to do with our mapping
from the YFP copy number to some arbitrary fluorescent value as computed from
the images taken with the microscope. The data processing inequality, as
derived by Shannon himself, tells us that for a Markov chain of the form
$c\rightarrow g\rightarrow f(g)$ it must be true that [3]

$I(g;c)\geq I(f(g);c),$ (S191)

meaning that information can only be lost when mapping from the real
relationship between gene expression and inducer concentration to a
fluorescence value.

On top of that, given the limited number of samples that we have access to
when computing the channel capacity, there is a bias in our estimate given
this undersampling. The definition of accurate unbiased descriptors of the
mutual information is still an area of active research. For our purposes we
will use the method described in [68]. The basic idea of the method is to
write the mutual information as a series expansion in terms of inverse powers
of the sample size, i.e.

$I_{\text{biased}}=I_{\infty}+\frac{a_{1}}{N}+\frac{a_{2}}{N^{2}}+\cdots,$
(S192)

where $I_{\text{biased}}$ is the biased estimate of the mutual information as
computed from experimental data, $I_{\infty}$ is the quantity we would like to
estimate, being the unbiased mutual information when having access to infinity
number of experimental samples, and the coefficients $a_{i}$ depend on the
underlying distribution of the signal and the response. This is an empirical
choice to be tested. Intuitively this choice satisfies the limit that as the
number of samples from the distribution grows, the empirical estimate of the
mutual information $I_{\text{biased}}$ should get closer to the actual value
$I_{\infty}$.

In principle for a good number of data points the terms of higher order become
negligible. So we can write the mutual information as

$I_{\text{biased}}\approx I_{\infty}+\frac{a_{1}}{N}+\mathcal{O}(N^{-2}).$
(S193)

This means that if this particular arbitrary choice of functional form is a
good approximation, when computing the mutual information for varying number
of samples - by taking subsamples of the experimental data - we expect to find
a linear relationship as a function of the inverse of these number of data
points. From this linear relationship the intercept is a bias-corrected
estimate of the mutual information. We can therefore bootstrap the data by
taking different sample sizes and then use the Blahut-Arimoto algorithm we
implemented earlier to estimate the biased channel capacity. We can then fit a
line and extrapolate for when $1/N=0$ which corresponds to our unbiased
estimate of the channel capacity.

Let’s go through each of the steps to illustrate the method. Fig. S31 show a
typical data set for a strain with an O2 binding site
($\Delta\varepsilon_{r}=-13.9\;k_{B}T$) and $R=260$ repressors per cell. Each
of the distributions in arbitrary units is binned into a specified number of
bins to build matrix $\mathbf{Q}$.

Figure S31: Single cell fluorescence distributions for different inducer
concentrations. Fluorescence distribution histogram (A) and cumulative
distribution function (B) for a strain with 260 repressors per cell and a
binding site with binding energy $\Delta\varepsilon_{r}=-13.9\;k_{B}T$. The
different curves show the single cell fluorescence distributions under the 12
different IPTG concentrations used throughout this work. The triangles in (A)
show the mean of each of the distributions.

Given a specific number of bins used to construct $\mathbf{Q}$, we subsample a
fraction of the data and compute the channel capacity for such matrix using
the Blahut-Arimoto algorithm. Fig. S32 shows an example where 50% of the data
on each distribution from Fig. S31 was sampled and binned into 100 equal bins.
The counts on each of these bins are then normalized and used to build matrix
$\mathbf{Q}$ that is then fed to the Blahut-Arimoto algorithm. We can see that
for these 200 bootstrap samples the channel capacity varies by $\approx$ 0.1
bits. Not a significant variability, nevertheless we consider that it is
important to bootstrap the data multiple times to get a better estimate of the
channel capacity.

Figure S32: Channel capacity bootstrap for experimental data. Cumulative
distribution function of the resulting channel capacity estimates obtained by
subsampling 200 times 50% of each distribution shown in Fig. S31, binning it
into 100 bins, and feeding the resulting $\mathbf{Q}$ matrix to the Blahut-
Arimoto algorithm.

Eq. S193 tells us that if we subsample each of the distributions from Fig. S31
at different fractions, and plot them as a function of the inverse sample size
we will find a linear relationship if the expansion of the mutual information
is valid. To test this idea we repeated the bootstrap estimate of Fig. S32
sampling 10%, 20%, and so on until taking 100% of the data. We repeated this
for different number of bins since a priori for arbitrary units of
fluorescence we do not have a way to select the optimal number of bins. Fig.
S33 shows the result of these estimates. We can see that the linear
relationship proposed in Eq. S193 holds true for all number of bins selected.
We also note that the value of the intercept of the linear regression varies
depending on the number of bins.

Figure S33: Inverse sample size vs channel capacity. As indicated in Eq. S193
if the channel capacity obtained for different subsample sizes of the data is
plotted against the inverse sample size there must exist a linear relationship
between these variables. Here we perform 15 bootstrap samples of the data from
Fig. S31, bin these samples using different number of bins, and perform a
linear regression (solid lines) between the bootstrap channel capacity
estimates, and the inverse sample size.

To address the variability in the estimates of the unbiased channel capacity
$I_{\infty}$ we again follow the methodology suggested in [68]. We perform the
data subsampling and computation of the channel capacity for a varying number
of bins. As a control we perform the same procedure with shuffled data, where
the structure that connects the fluorescence distribution to the inducer
concentration input is lost. The expectation is that this control should give
a channel capacity of zero if the data is not “over-binned.” Once the number
of bins is too high, we would expect some structure to emerge in the data that
would cause the Blahut-Arimoto algorithm to return non-zero channel capacity
estimates.

Fig. S34 shows the result of the unbiased channel capacity estimates obtained
for the data shown in Fig. S31. For the blue curve we can distinguish three
phases:

  1. 1.

A rapid increment from 0 bits to about 1.5 bits as the number of bins
increases.

  2. 2.

A flat region between $\approx$ 50 and 1000 bins.

  3. 3.

A second rapid increment for large number of bins.

We can see that the randomized data presents two phases only:

  1. 1.

A flat region where there is, as expected no information being processed since
the structure of the data was lost when the data was shuffled.

  2. 2.

A region with fast growth of the channel capacity as the over-binning
generates separated peaks on the distribution, making it look like there is
structure in the data.

We take the flat region of the experimental data ($\approx$ 100 bins) to be
our best unbiased estimate of the channel capacity from this experimental
dataset.

Figure S34: Channel capacity as a function of the number of bins. Unbiased
channel capacity estimates obtained from linear regressions as in Fig. S33.
The blue curve show the estimates obtained from the data shown in Fig. S31.
The orange curve is generated from estimates where the same data is shuffled,
loosing the relationship between fluorescence distributions and inducer
concentration.

###  S8.3 Assumptions involved in the computation of the channel capacity

An interesting suggestion by Professor Gasper Tkacik was to dissect the
different physical assumptions that went into the construction of the input-
output function $P(p\mid c)$, and their relevance when comparing the
theoretical channel capacities with the experimental inferences. In what
follows we describe the relevance of four important aspects that all affect
the predictions of the information processing capacity.

##### (i) Cell cycle variability.

We think that the inclusion of the gene copy number variability during the
cell cycle and the non-Poissoninan protein degradation is a key component to
our estimation of the input-output functions and as a consequence of the
channel capacity. This variability in gene copy number is an additional source
of noise that systematically decreases the ability of the system to resolve
different inputs. The absence of the effects that the gene copy number
variability and the protein partition has on the information processing
capacity leads to an overestimate of the channel capacity as shown in Fig.
S35. Only when these noise sources are included in our inferences is that we
get to capture the experimental channel capacities with no further fit
parameters.

Figure S35: Comparison of channel capacity predictions for single- and multi-
promoter models. Channel capacity for the multi-promoter model (solid lines)
vs. the single-promoter steady state model (dot-dashed lines) as a function of
repressor copy numbers for different repressor-DNA binding energies. The
single-promoter model assumes Poissonian protein degradation ($\gamma_{p}>0$)
and steady state, while the multi-promoter model accounts for gene copy number
variability and during the cell cycle and has protein degradation as an effect
due to dilution as cells grow and divide.

##### (ii) Non-Gaussian noise distributions.

For the construction of the probability distributions used in the main text
(Fig. 4) we utilized the first 6 moments of the protein distribution. The
maximum entropy formalism tells us that the more constraints we include in the
inference, the closer the maximum entropy distribution will be to the real
distribution. But a priori there is no way of knowing how many moments should
be included in order to capture the essence of the distribution. In principle
two moments could suffice to describe the entire distribution as happens with
the Gaussian distribution. To compare the effect that including more or less
constraints on the maximum entropy inference we constructed maximum entropy
distributions using an increasing number of moments from 2 to 6. We then
computed the Kullback-Leibler divergence $D_{KL}$ of the form

$D_{KL}(P_{6}(p\mid c)||P_{i}(p\mid c))=\sum_{p}P_{6}(p\mid
c)\log_{2}{P_{6}(p\mid c)\over P_{i}(p\mid c)},$ (S194)

where $P_{i}(p\mid c)$ is the maximum entropy distribution constructed with
the first $i$ moments, $i\in\\{2,3,4,5,6\\}$. Since the Kullback-Leibler
divergence $D_{KL}(P||Q)$ can be interpreted as the amount of information lost
by assuming the incorrect distribution $Q$ when the correct distribution is
$P$, we used this metric as a way of how much information we would have lost
by using less constraints compared to the six moments used in the main text.

Fig. S36 shows this comparison for different operators and repressor copy
numbers. We can see from here that using less moments as constraints gives
basically the same result. This is because most of the values of the Kullback-
Leibler divergence are significantly smaller than 0.1 bits, and the entropy of
these distributions is in general $>10$ bits, so we would lose less than 1% of
the information contained in these distributions by utilizing only two moments
as constraints. Therefore the use of non-Gaussian noise is not an important
feature for our inferences.

Figure S36: Measuring the loss of information by using different number of
constraints. The Kullback-Leibler divergence was computed between the maximum
entropy distribution constructed using the first 6 moments of the distribution
and a variable number of moments.

##### (iii) Multi-state promoter.

This particular point is something that we are still exploring from a
theoretical perspective. We have shown that in order to capture the single-
molecule mRNA FISH data a single-state promoter wouldn’t suffice. This model
predicts a Poisson distribution as the steady-state and the data definitely
shows super Poissonian noise. Given the bursty nature of gene expression we
opt to use a two-state promoter where the states reflect effective
transcriptionally “active” and “inactive” states. We are currently exploring
alternative formulations of this model to turn it into a single state with a
geometrically distributed burst-size.

##### (iv) Optimal vs Log-flat Distributions.

The relevance of having use the Blahut-Arimoto algorithm to predict the
maximum mutual information between input and outputs was just to understand
the best case scenario. We show the comparison between theoretical and
experimental input-output functions $P(p\mid c)$ in Fig. S26. Given the good
agreement between these distributions we could compute the mutual information
$I(c;p)$ for any arbitrary input distribution $P(c)$ and obtain a good
agreement with the corresponding experimental mutual information.

The reason we opted to specifically report the mutual information at channel
capacity was to put the results in a context. By reporting the upper bound in
performance of these genetic circuits we can start to dissect how different
molecular parameters such as repressor-DNA binding affinity or repressor copy
number affect the ability of this genetic circuit to extract information from
the environmental state.

##  S9 Empirical fits to noise predictions

(Note: The Python code used for the calculations presented in this section can
be found in the following link as an annotated Jupyter notebook)

In Fig. 3(C) in the main text we show that our minimal model has a systematic
deviation on the gene expression noise predictions compared to the
experimental data. This systematics will need to be addressed on an improved
version of the minimal model presented in this work. To guide the insights
into the origins of this systematic deviation in this appendix we will explore
empirical modifications of the model to improve the agreement between theory
and experiment.

###  S9.1 Multiplicative factor for the noise

The first option we will explore is to modify our noise predictions by a
constant multiplicative factor. This means that we assume the relationship
between our minimal model predictions and the data for noise in gene
expression are of the from

$\text{noise}_{\text{exp}}=\alpha\cdot\text{noise}_{\text{theory}},$ (S195)

where $\alpha$ is a dimensionless constant to be fit from the data. The data,
especially in Fig. S19 suggests that our predictions are within a factor of
$\approx$ two from the experimental data. To further check that intuition we
performed a weighted linear regression between the experimental and
theoretical noise measurements. The weight for each datum was taken to be
proportional to the bootstrap errors in the noise estimate, this to have
poorly determined noises weigh less during the regression. The result of this
regression with no intercept shows exactly that a factor of two systematically
improves the theoretical vs. experimental predictions. Fig. S37 shows the
improved agreement when the theoretical predictions for the noise are
multiplied by $\approx 1.5$.

Figure S37: Multiplicative factor to improve theoretical vs. experimental
comparison of noise in gene expression. Theoretical vs. experimental noise
both in linear (left) and log (right) scale. The dashed line shows the
identity line of slope 1 and intercept zero. All data are colored by the
corresponding value of the experimental fold-change in gene expression as
indicated by the color bar. The $x$-axis was multiplied by a factor of
$\approx 1.5$ as determined by a linear regression from the data in Fig. S18.
Each datum represents a single date measurement of the corresponding strain
and IPTG concentration with $\geq 300$ cells. The points correspond to the
median, and the error bars correspond to the 95% confidence interval as
determined by 10,000 bootstrap samples.

For completeness Fig. S38 shows the noise in gene expression as a function of
the inducer concentration including this factor of $\approx 1.5$. It is clear
that overall a simple multiplicative factor improves the predictive power of
the model.

Figure S38: Protein noise of the regulated promoter with multiplicative
factor. Comparison of the experimental noise for different operators ((A) O1,
$\Delta\varepsilon_{r}=-15.3\;k_{B}T$, (B) O2,
$\Delta\varepsilon_{r}=-13.9\;k_{B}T$, (C) O3,
$\Delta\varepsilon_{r}=-9.7\;k_{B}T$) with the theoretical predictions for the
the multi-promoter model. A linear regression revealed that multiplying the
theoretical noise prediction by a factor of $\approx 1.5$ would improve
agreement between theory and data. Points represent the experimental noise as
computed from single-cell fluorescence measurements of different E. coli
strains under 12 different inducer concentrations. Dotted line indicates plot
in linear rather than logarithmic scale. Each datum represents a single date
measurement of the corresponding strain and IPTG concentration with $\geq 300$
cells. The points correspond to the median, and the error bars correspond to
the 95% confidence interval as determined by 10,000 bootstrap samples. White-
filled dots are plot at a different scale for better visualization.

###  S9.2 Additive factor for the noise

As an alternative way to empirically improve the predictions of our model we
will now test the idea of an additive constant. What this means is that our
minimal model underestimates the noise in gene expression as

$\text{noise}_{\text{exp}}=\beta+\text{noise}_{\text{theory}},$ (S196)

where $\beta$ is an additive constant to be determined from the data. As with
the multiplicative constant we performed a regression to determine this
empirical additive constant comparing experimental and theoretical gene
expression noise values. We use the error in the 95% bootstrap confidence
interval as a weight for the linear regression. Fig. S39 shows the resulting
theoretical vs. experimental noise where $\beta\approx 0.2$. We can see a
great improvement in the agreement between theory and experiment with this
additive constant

Figure S39: Additive factor to improve theoretical vs. experimental comparison
of noise in gene expression. Theoretical vs. experimental noise both in linear
(left) and log (right) scale. The dashed line shows the identity line of slope
1 and intercept zero. All data are colored by the corresponding value of the
experimental fold-change in gene expression as indicated by the color bar. A
value of $\approx 0.2$ was added to all values in the $x$-axis as determined
by a linear regression from the data in Fig. S18. Each datum represents a
single date measurement of the corresponding strain and IPTG concentration
with $\geq 300$ cells. The points correspond to the median, and the error bars
correspond to the 95% confidence interval as determined by 10,000 bootstrap
samples.

For completeness Fig. S40 shows the noise in gene expression as a function of
the inducer concentration including this additive factor of $\beta\approx
0.2$. If anything, the additive factor seems to improve the agreement between
theory and data even more than the multiplicative factor.

Figure S40: Protein noise of the regulated promoter with additive factor.
Comparison of the experimental noise for different operators ((A) O1,
$\Delta\varepsilon_{r}=-15.3\;k_{B}T$, (B) O2,
$\Delta\varepsilon_{r}=-13.9\;k_{B}T$, (C) O3,
$\Delta\varepsilon_{r}=-9.7\;k_{B}T$) with the theoretical predictions for the
the multi-promoter model. A linear regression revealed that an additive factor
of $\approx 0.2$ to the the theoretical noise prediction would improve
agreement between theory and data. Points represent the experimental noise as
computed from single-cell fluorescence measurements of different E. coli
strains under 12 different inducer concentrations. Dotted line indicates plot
in linear rather than logarithmic scale. Each datum represents a single date
measurement of the corresponding strain and IPTG concentration with $\geq 300$
cells. The points correspond to the median, and the error bars correspond to
the 95% confidence interval as determined by 10,000 bootstrap samples. White-
filled dots are plot at a different scale for better visualization.

###  S9.3 Correction factor for channel capacity with multiplicative factor