# Supplementary Information:
Deriving a genetic regulatory network from an optimization principle

Thomas R. Sokolowski, Thomas Gregor, William Bialek & Gašper Tkačik

###### Contents

  1. 1 Spatial-stochastic gene expression model
    1. 1.1 Principal ansatz
    2. 1.2 Mean expression levels
    3. 1.3 Covariances
      1. 1.3.1 Local inter-gene covariance
      2. 1.3.2 Local variance
      3. 1.3.3 Inter-volume covariance for the same gene
      4. 1.3.4 A single self-regulating gene without diffusive coupling
      5. 1.3.5 Short-correlations assumption
    4. 1.4 Noise powers
    5. 1.5 Regulatory function
    6. 1.6 Combined mRNA and protein model
      1. 1.6.1 Incorporating mRNA into the model
      2. 1.6.2 Normalized differential equation for the means of mRNA and protein
      3. 1.6.3 Approximations for computing noise covariances
      4. 1.6.4 Normalized noise covariances with mRNA
  2. 2 Numerical computation of means and covariance matrix of gap gene expression
    1. 2.1 Forward integration of the means
    2. 2.2 Numerical computation of the covariance matrix
  3. 3 Calculation of positional information
  4. 4 Optimization procedure
    1. 4.1 Parameter changes and acceptance
    2. 4.2 Cooling protocol
    3. 4.3 Final gradient descent
    4. 4.4 Exploration of Drosophila-like patterns via constrained optimization
  5. 5 Parameters, simulation geometry and optimization constraints
    1. 5.1 Fixed parameters and geometric constraints
    2. 5.2 Optimization parameters
    3. 5.3 Resource constraints
    4. 5.4 Pattern stability constraint
  6. 6 Details of data analysis
    1. 6.1 Computing positional error profiles
    2. 6.2 Determining number of slopes
    3. 6.3 Measuring resource allocation
    4. 6.4 Determining regulatory interactions
  7. 7 Additional analyses of the optimal ensemble
  8. 8 Optimizations with altered system components
    1. 8.1 Varying gap protein diffusion constant
    2. 8.2 Altered number of gap genes and maternal inputs
    3. 8.3 Bounding maximal regulatory strength

###### List of Supplementary Figures

  1. 1 Estimating extrinsic noise from data.
  2. 2 Approximate computation of positional information via optimal decoding.
  3. 3 Optimality check for the Drosophila-like solution.
  4. 4 Detailed map of regulatory contributions for the Drosophila-like solution.
  5. 5 Additional patterning phenotypes in the optimal and random ensembles.
  6. 6 _Ab initio_ optimized solutions similar to the Drosophila WT expression patterns.
  7. 7 Optimal solutions with 3 gap genes have localized positional defects.
  8. 8 Altered number of gap genes and maternal morphogen gradients.
  9. 9 Altered number of gap genes and maternal morphogen gradients, cont’d.
  10. 10 Feed-forward regulation by maternal morphogen gradients only.

##  1 Spatial-stochastic gene expression model

###  1.1 Principal ansatz

We model the stochastic gene expression dynamics of the gap gene system using
a set of coupled stochastic differential equations describing the time
evolution of gene products of all gap gene species $\alpha$ present in the
system in all nuclear volumes $i$. These dynamics are governed by a production
term that depends on the local levels of maternal inputs and gap gene proteins
in volume $i$, a degradation term, a diffusion term describing spatial
exchange with neighboring volumes, and several noise sources that we specify
in detail later. For clarity here we first derive a model that conceptually
considers protein production only and skips over the intermediate mRNA
production step. Later, in Sec. 1.6, we construct an extended model that
incorporates mRNA in an approximate fashion and conserves the principal
structure of the ansatz below, such that all derivations made in this section
also apply to the extended model.

Let us denote the normalized protein copy number of gene $\alpha$ in volume
$i$ by $g_{i}^{\alpha}$ and non-normalized protein copy numbers by
$G_{i}^{\alpha}$; we will explain the choice of the normalization factor
$N_{\rm m}$ further below. The dynamics of proteins of gene species $\alpha$
in volume $i$ then are governed by:

$\displaystyle\partial_{t}g_{i}^{\alpha}$
$\displaystyle=\frac{r_{\alpha}(\mathbf{g}_{i},\mathbf{c}_{i})}{N_{\rm
m}}-\frac{1}{\tau_{\alpha}}g_{i}^{\alpha}-h_{\alpha}\sum_{n}\left(g_{i}^{\alpha}-g_{n}^{\alpha}\right)+\sum_{k}\Gamma_{ik}^{\alpha}\eta_{k}$
(1)

where the production function $r_{\alpha}(\mathbf{g}_{i},\mathbf{c}_{i})$
depends on the set of all maternal and gap gene protein levels,
$\mathbf{c}_{i}=(c_{i}^{1},c_{i}^{2},c_{i}^{3})$ and
$\mathbf{g}_{i}=(g_{i}^{1},g_{i}^{2},g_{i}^{3},g_{i}^{4})$, respectively, in
the volume; note that the ansatz above remains valid for arbitrary numbers of
gap genes and maternal inputs. $\tau_{\alpha}$ is a (gene-specific) product
lifetime, while $h_{\alpha}\equiv D_{\alpha}/\Delta^{2}$ denotes a gene-
specific “hopping rate” for spatial exchange of gene products between
neighboring nuclei; herein $D_{\alpha}$ is the diffusion coefficient and
$\Delta$ the distance between neighboring nuclear volumes, assumed to be the
same for all nuclei (uniform spacing). The sum of the diffusion term runs over
all $2d$ nearest neighbors $n$ of $i$, whose number varies with the dimension
$d$ of the nuclear lattice; we considered a cylindrical lattice in $d=2$ and a
uniformly spaced one-dimensional lattice ($d=1$), but keep $d$ arbitrary in
the following derivations. The $k$-sum runs over all noise sources in the
system that affect the copy number of gene $\alpha$ in volume $i$. In this
sum, the noise powers $\Gamma_{ik}^{\alpha}$ are understood to be normalized
via $N_{\rm m}$ as well, and will be specified later, in Sec. 1.4. All noise
sources $\eta_{k}$ are assumed to be zero-mean Gaussian white noise, i.e.

$\displaystyle\langle\eta_{k}(t)\rangle$ $\displaystyle=0$
$\displaystyle\langle\eta_{k}(t)\eta_{k}(t+dt)\rangle$
$\displaystyle=\delta(dt)$ (2)

By default we assume that all genes can be expressed with the same maximal
mean copy number $N_{\rm m}$, i.e., each nucleus can potentially produce
$N_{\rm m}$ gene products in stationary state on average, when all stochastic
fluctuations are averaged out. Our model straightforwardly extends to the more
general case in which the maximal mean copy numbers of the zygotically
expressed species $\zeta$ can be different. We can then take the largest mean
copy number as a reference:

$\displaystyle N_{\rm m}\equiv\max_{\zeta}\left\langle G_{\zeta}\right\rangle$
(3)

All other copy numbers then will be smaller or equal to $N_{\rm m}$, and we
can account for this by introducing additional, species-specific production
rates in the production function

$\displaystyle\frac{r_{\alpha}(\mathbf{g}_{i},\mathbf{c}_{i})}{N_{\rm m}}$
$\displaystyle=\frac{1}{N_{\rm m}}\left[r_{\rm b}^{\alpha}+(r_{\rm
on}^{\alpha}-r_{\rm
b}^{\alpha})f_{\alpha}(\mathbf{g}_{i},\mathbf{c}_{i})\right]$
$\displaystyle=\frac{1}{\tau}\left[\frac{r_{\rm b}^{\alpha}}{r_{\rm
on}}+\frac{r_{\rm on}^{\alpha}-r_{\rm b}^{\alpha}}{r_{\rm
on}}f_{\alpha}(\mathbf{g}_{i},\mathbf{c}_{i})\right]$
$\displaystyle\equiv\frac{1}{\tau}\left[\hat{\rho}_{\rm
b}^{\alpha}+(\hat{\rho}_{\rm on}^{\alpha}-\hat{\rho}_{\rm
b}^{\alpha})f_{\alpha}(\mathbf{g}_{i},\mathbf{c}_{i})\right]$ (4)

where we denote the production rate and lifetime of the species with the
largest copy number $N_{\rm m}$ by $r_{\rm on}$ and $\tau$, respectively, i.e.
$N_{\rm m}\equiv r_{\rm on}\tau$. Rate $r_{\rm b}^{\alpha}$ is a small basal
expression rate, while $\hat{\rho}_{\rm b}^{\alpha}$ denotes the same rate
measured in units of $r_{\rm on}$. Note that the “normalized” rate
$\hat{\rho}_{\rm on}^{\alpha}\equiv r_{\rm on}^{\alpha}/r_{\rm on}$ then can
be smaller or larger than 1, because the lifetimes $\tau$ and $\tau_{\alpha}$
are not taken into account in the normalization.

The regulation function $f_{\alpha}(\mathbf{g}_{i},\mathbf{c}_{i})$ also can
be imposed later, such that in the following calculations we will leave the
specific form of $f_{\alpha}$ open for tractability. In our model we opted for
a regulatory function based on the MWC model, described in detail in 1.5. In
the following we will also operate with the function

$\displaystyle
F_{\alpha}(\mathbf{g}_{i},\mathbf{c}_{i})\equiv\frac{1}{\tau}\left[\hat{\rho}_{\rm
b}^{\alpha}+(\hat{\rho}_{\rm on}^{\alpha}-\hat{\rho}_{\rm
b}^{\alpha})f_{\alpha}(\mathbf{g}_{i},\mathbf{c}_{i})\right]-\frac{1}{\tau_{\alpha}}g_{i}^{\alpha}$
(5)

grouping the production and degradation terms for $g_{i}^{\alpha}$.

We will now proceed with the calculation of the means and covariances of the
spatially coupled gene expression model, needed for later computation of the
positional information contained in the expression pattern.

###  1.2 Mean expression levels

The stationary mean expression levels
$\bar{g}_{i}^{\alpha}\equiv\left\langle\bar{g}_{i}^{\alpha}\right\rangle$ in
principle could be obtained from Eq. (1) by setting the time derivatives on
the left-hand side to zero and averaging over the right-hand side.
Practically, however, this is only feasible if the production term
$\left\langle r_{\alpha}(\mathbf{g}_{i},\mathbf{c}_{i})\right\rangle$ and the
regulatory function $f_{\alpha}(\mathbf{g}_{i},\mathbf{c}_{i})$ in particular
could be approximated by (the leading terms of) its Taylor-expansion around
$\mathbf{\bar{g}}_{i}$ [see Eq. (10) below]. Then, since $\forall
k:\left\langle\eta_{k}\right\rangle=0$, this would yield the following coupled
system

$\displaystyle\partial_{t}\bar{g}_{i}^{\alpha}=\frac{r_{\alpha}(\mathbf{\bar{g}}_{i},\mathbf{c}_{i})}{N_{\rm
m}}-\frac{1}{\tau_{\alpha}}\bar{g}_{i}^{\alpha}-h_{\alpha}\sum_{n}\left(\bar{g}_{i}^{\alpha}-\bar{g}_{n}^{\alpha}\right)\equiv
0\quad\Leftrightarrow$
$\displaystyle\left(2dh_{\alpha}+\frac{1}{\tau_{\alpha}}\right)\bar{g}_{i}^{\alpha}=\frac{r_{\alpha}(\mathbf{\bar{g}}_{i},\mathbf{c}_{i})}{N_{\rm
m}}+h_{\alpha}\sum_{n}\bar{g}_{n}^{\alpha}\quad\Leftrightarrow$
$\displaystyle\bar{g}_{i}^{\alpha}=\frac{\tau_{\alpha}}{\tau}\frac{\hat{\rho}_{\rm
b}^{\alpha}+(\hat{\rho}_{\rm on}^{\alpha}-\hat{\rho}_{\rm
b}^{\alpha})f_{\alpha}(\mathbf{\bar{g}}_{i},\mathbf{c}_{i})}{1+2d\cdot
h_{\alpha}\tau_{\alpha}}+\frac{h_{\alpha}\tau_{\alpha}}{1+2d\cdot
h_{\alpha}\tau_{\alpha}}\sum_{n}\bar{g}_{n}^{\alpha}$

where $d$ denotes the dimensionality of the spatial lattice. Here the
prefactors have the same form and interpretation as in the single-
input/single-output system studied earlier by us [1], and reflect the fact
that with increasing diffusive coupling (increasing $h_{\alpha}$) the
expression level in volume $i$ is increasingly determined by the spatial
average of the neighboring expression levels (second term on the right), such
that neighboring levels increasingly equilibrate in spite of (potentially)
different local production rates (first term on the right).

In practice the coupled system defined above is hard to solve because the
regulatory function $f_{\alpha}(\mathbf{g}_{i},\mathbf{c}_{i})$ (yet to
impose) can be highly nonlinear. Moreover, it is not a priori clear around
which expression levels one should linearize, as these levels arise as a
consequence of the non-linear dynamics. We therefore opted for numerical
forward-integration of the means, which do not require linearizing
$f_{\alpha}(\mathbf{g}_{i},\mathbf{c}_{i})$ and can be carried out with
sufficient computational efficiency for our needs; the details are described
in Sec. 2.1 along with the other numerical methods that we used for solving
and optimizing the spatial-stochastic model derived here.

###  1.3 Covariances

Assuming that we have determined the full set of mean expression levels in the
system, let us now compute the steady-state covariances between the product
copy number $g_{i}^{\alpha}$ of gene $\alpha$ in volume $i$ and the product
copy number $g_{j}^{\beta}$ of gene $\beta$ in volume $j$, defined as:

$\displaystyle C_{ij}^{\alpha\beta}\equiv\left\langle\delta
g_{i}^{\alpha}\delta g_{j}^{\beta}\right\rangle$ (7)

Here the fluctuations $\delta g_{i}^{\alpha}$ are assumed to be linear and
symmetric around the mean, i.e. $g_{i}^{\alpha}=\bar{g}_{i}^{\alpha}+\delta
g_{i}^{\alpha}$. Note the symmetry
$C_{ij}^{\alpha\beta}=C_{ji}^{\beta\alpha}$.

We compute $C_{ij}^{\alpha\beta}$ using the same formula as in our earlier
work on a single-input/single-output system [1], which correlates the
fluctuations of local expresssion level $g_{i}^{\alpha}$ with the “forces”
(production, degradation and diffusion) driving expression level
$g_{j}^{\beta}$ and vice versa, and all noise sources affecting the respective
expression levels:

$\displaystyle\partial_{t}\left\langle\delta g_{i}^{\alpha}\delta
g_{j}^{\beta}\right\rangle$
$\displaystyle=\underset{\mathfrak{F}_{\beta}}{\underbrace{\left\langle\delta
g_{i}^{\alpha}F_{\beta}(\mathbf{g}_{j},\mathbf{c}_{j})\vphantom{{}^{\beta}}\right\rangle}}+\underset{\mathfrak{F}_{\alpha}}{\underbrace{\left\langle\delta
g_{j}^{\beta}F_{\alpha}(\mathbf{g}_{i},\mathbf{c}_{i})\right\rangle}}$
$\displaystyle\quad+\underset{\mathfrak{D}_{\beta}}{\underbrace{h_{\beta}\left\langle\delta
g_{i}^{\alpha}\sum_{n_{j}}\left(g_{n_{j}}^{\beta}-g_{j}^{\beta}\right)\right\rangle}}+\underset{\mathfrak{D}_{\alpha}}{\underbrace{h_{\alpha}\left\langle\delta
g_{j}^{\beta}\sum_{n_{i}}\left(g_{n_{i}}^{\alpha}-g_{i}^{\alpha}\right)\right\rangle}}$
$\displaystyle\quad+\underset{\mathfrak{G}}{\underbrace{\sum_{k}\left\langle\Gamma_{ik}^{\alpha}\Gamma_{jk}^{\beta}\right\rangle}}$
(8)

We will now compute the different terms denoted by capital blackletters
separately. Following our usual strategy, we will specify and compute the
noise term $\mathfrak{G}$ only at the end.

Using simple algebraic steps, we can write out $\mathfrak{F}_{\alpha}$ as:

$\displaystyle\mathfrak{F}_{\beta}$ $\displaystyle=\frac{\hat{\rho}_{\rm
on}^{\beta}-\hat{\rho}_{\rm b}^{\beta}}{\tau}\left\langle\delta
g_{i}^{\alpha}f_{\beta}(\mathbf{g}_{j},\mathbf{c}_{j})\right\rangle-\frac{1}{\tau_{\beta}}\left\langle\delta
g_{i}^{\alpha}\delta g_{j}^{\beta}\right\rangle$ (9)

To proceed further, we have to compute how fluctuations in all the involved
gene copy numbers $\delta\mathbf{g}_{j}$ propagate through the regulatory
function $f_{\beta}$. Here, in contrast to the computation of the mean
expression levels (Sec. 1.2), we pursue a linearization ansatz assuming the
fluctuations to be small compared to the allowed expression range. We can
expand around the mean $\mathbf{\bar{g}}_{j}$, retaining terms up to first
order, as follows:

$\displaystyle f_{\beta}(\mathbf{g}_{j},\mathbf{c}_{j})$
$\displaystyle=f_{\beta}(\mathbf{\bar{g}}_{j}+\delta\mathbf{g}_{j},\mathbf{c}_{j})\simeq
f_{\beta}(\mathbf{\bar{g}}_{j},\mathbf{c}_{j})+\sum_{\zeta}\left(\frac{\partial
f_{\beta}}{\partial
g_{j}^{\zeta}}\right)_{\mathbf{\bar{g}}_{j},\mathbf{c}_{j}}\delta
g_{j}^{\zeta}$ (10)

Herein $\zeta$ loops over all zygotic (gap gene) species in the system. From
this it follows that

$\displaystyle\left\langle\delta
g_{i}^{\alpha}f_{\beta}(\mathbf{g}_{j},\mathbf{c}_{j})\right\rangle$
$\displaystyle=\sum_{\zeta}\left(\frac{\partial f_{\beta}}{\partial
g_{j}^{\zeta}}\right)_{\mathbf{\bar{g}}_{j},\mathbf{c}_{j}}\left\langle\delta
g_{i}^{\alpha}\delta
g_{j}^{\zeta}\right\rangle\equiv\sum_{\zeta}f_{\beta}^{\prime\zeta}(\mathbf{\bar{g}}_{j},\mathbf{c}_{j})C_{ij}^{\alpha\zeta}$
(11)

where we introduce the short notation $f_{\beta}^{{}^{\prime}\zeta}$ for
indicating that the derivative of the regulatory function $f_{\beta}$ is taken
with respect to zygotic species $\zeta$. Note that here we allow for
$\zeta=\beta$, which corresponds to self-regulation of species $\beta$. With
simple intermediate steps we obtain

$\displaystyle\mathfrak{F}_{\beta}$ $\displaystyle=\frac{\hat{\rho}_{\rm
on}^{\beta}-\hat{\rho}_{\rm
b}^{\beta}}{\tau}\sum_{\zeta}f_{\beta}^{\prime\zeta}(\mathbf{\bar{g}}_{j},\mathbf{c}_{j})C_{ij}^{\alpha\zeta}-\frac{1}{\tau_{\beta}}C_{ij}^{\alpha\beta}$
$\displaystyle\mathfrak{F}_{\alpha}$ $\displaystyle=\frac{\hat{\rho}_{\rm
on}^{\alpha}-\hat{\rho}_{\rm
b}^{\alpha}}{\tau}\sum_{\zeta}f_{\alpha}^{\prime\zeta}(\mathbf{\bar{g}}_{i},\mathbf{c}_{i})C_{ji}^{\beta\zeta}-\frac{1}{\tau_{\alpha}}C_{ij}^{\alpha\beta}$
(12)

where the calculation for $\mathfrak{F}_{\alpha}$ is completely analogous, and
in the last step we have used $C_{ji}^{\beta\alpha}=C_{ij}^{\alpha\beta}$.

For the diffusive coupling terms $\mathfrak{D}_{\beta}$ and
$\mathfrak{D}_{\alpha}$, after taking the average, we get:

$\displaystyle\mathfrak{D}_{\beta}$
$\displaystyle=h_{\beta}\sum_{n_{j}}C_{in_{j}}^{\alpha\beta}-2dh_{\beta}C_{ij}^{\alpha\beta}$
$\displaystyle\mathfrak{D}_{\alpha}$
$\displaystyle=h_{\alpha}\sum_{n_{i}}C_{jn_{i}}^{\beta\alpha}-2dh_{\alpha}C_{ij}^{\alpha\beta}$
(13)

Here, $d$ again is the lattice dimension.

Putting all terms together, setting $\partial_{t}[...]=0$ in Eq. (8), and
collecting parts that only contain covariance $C_{ij}^{\alpha\beta}$ on one
side, we obtain the following:

$\displaystyle\Bigg{[}\underset{\text{from
}\mathfrak{F}_{\alpha}}{\underbrace{\frac{1}{\tau_{\alpha}}\vphantom{\frac{1}{\tau_{\beta}}}}}+\underset{\text{from
}\mathfrak{F}_{\beta}}{\underbrace{\frac{1}{\tau_{\beta}}\vphantom{\frac{1}{\tau_{\beta}}}}}+\underset{\text{from
}\mathfrak{D}_{\alpha}\text{ and
}\mathfrak{D}_{\beta}}{\underbrace{2d(h_{\alpha}+h_{\beta})\vphantom{\frac{1}{\tau_{\beta}}}}}\Bigg{]}C_{ij}^{\alpha\beta}=$
$\displaystyle\quad\frac{\hat{\rho}_{\rm on}^{\beta}-\hat{\rho}_{\rm
b}^{\beta}}{\tau}\sum_{\zeta}f_{\beta}^{\prime\zeta}(\mathbf{\bar{g}}_{j},\mathbf{c}_{j})C_{ij}^{\alpha\zeta}+\frac{\hat{\rho}_{\rm
on}^{\alpha}-\hat{\rho}_{\rm
b}^{\alpha}}{\tau}\sum_{\zeta}f_{\alpha}^{\prime\zeta}(\mathbf{\bar{g}}_{i},\mathbf{c}_{i})C_{ji}^{\beta\zeta}$
$\displaystyle\quad+h_{\beta}\sum_{n_{j}}C_{in_{j}}^{\alpha\beta}+h_{\alpha}\sum_{n_{i}}C_{jn_{i}}^{\beta\alpha}+\mathfrak{G}$
(14)

With the definitions

$\displaystyle T_{\alpha\beta}$ $\displaystyle\equiv
2\left[\frac{1}{\tau_{\alpha}}+\frac{1}{\tau_{\beta}}+2d(h_{\alpha}+h_{\beta})\right]^{-1}~{},$
(15) $\displaystyle\Lambda^{2}_{\alpha}$ $\displaystyle\equiv
h_{\alpha}T_{\alpha\beta}~{},\quad\quad\quad\Lambda^{2}_{\beta}\equiv
h_{\beta}T_{\alpha\beta}~{},$ (16)

which can be interpreted as an effective averaging time and effective
diffusive coupling lengths, respectively, we can rewrite the equation above in
the more compact final form:

$\displaystyle C_{ij}^{\alpha\beta}$
$\displaystyle=\frac{T_{\alpha\beta}}{2}\Bigg{[}\frac{\hat{\rho}_{\rm
on}^{\beta}-\hat{\rho}_{\rm
b}^{\beta}}{\tau}\sum_{\zeta}f_{\beta}^{\prime\zeta}(\mathbf{\bar{g}}_{j},\mathbf{c}_{j})C_{ij}^{\alpha\zeta}$
$\displaystyle\quad\quad\quad\quad\quad+\frac{\hat{\rho}_{\rm
on}^{\alpha}-\hat{\rho}_{\rm
b}^{\alpha}}{\tau}\sum_{\zeta}f_{\alpha}^{\prime\zeta}(\mathbf{\bar{g}}_{i},\mathbf{c}_{i})C_{ji}^{\beta\zeta}\Bigg{]}$
$\displaystyle\quad~{}+\frac{1}{2}\Bigg{[}\Lambda^{2}_{\beta}\sum_{n_{j}}C_{in_{j}}^{\alpha\beta}+\Lambda^{2}_{\alpha}\sum_{n_{i}}C_{jn_{i}}^{\beta\alpha}\Bigg{]}+\frac{T_{\alpha\beta}}{2}\mathfrak{G}$
(17)

This defines a linear system for the whole set of inter-gene/inter-position
covariances (or entries of the four-dimensional covariance matrix of the
system) $C_{ij}^{\alpha\beta}$, which can be solved after specifying the noise
term $\mathfrak{G}$ and imposing suitable boundary conditions. The equation
above has the same structure as in the single-input/single-output model
studied in our earlier work [1], with additional terms introduced by the
dependencies of the regulatory functions on the additional gap genes in the
system. Accordingly, by setting $\alpha=\beta$ and imposing that the
regulatory functions only depend on a single maternal input $c_{i}$, implying
$\forall\zeta:~{}f_{\alpha}^{\prime\zeta}(\mathbf{\bar{g}}_{j},\mathbf{c}_{j})=0=f_{\beta}^{\prime\zeta}(\mathbf{\bar{g}}_{j},\mathbf{c}_{j})$,
we recover the formula for the single gene system (see Sec. 1.3.3 below). We
will now consider this and other special cases of the general formula above,
including a simplified case which reproduces a known formula from our earlier
work as a consistency test. We will then proceed by imposing a “short-
correlation assumption”, which assumes that significant covariance
correlations only arise between nearby neighbor volumes and markedly
simplifies the (numerical) computation of the covariance matrix.

####  1.3.1 Local inter-gene covariance

By setting $i=j$ in Eq. (17) we obtain the formula for covariance between
$g_{\alpha}$ and $g_{\beta}$ in nuclear volume i:

$\displaystyle C_{ii}^{\alpha\beta}$
$\displaystyle=\frac{T_{\alpha\beta}}{2}\Bigg{[}\frac{\hat{\rho}_{\rm
on}^{\beta}-\hat{\rho}_{\rm
b}^{\beta}}{\tau}\sum_{\zeta}f_{\beta}^{\prime\zeta}(\mathbf{\bar{g}}_{i},\mathbf{c}_{i})C_{ii}^{\alpha\zeta}$
$\displaystyle\quad\quad\quad\quad\quad+\frac{\hat{\rho}_{\rm
on}^{\alpha}-\hat{\rho}_{\rm
b}^{\alpha}}{\tau}\sum_{\zeta}f_{\alpha}^{\prime\zeta}(\mathbf{\bar{g}}_{i},\mathbf{c}_{i})C_{ii}^{\beta\zeta}\Bigg{]}$
$\displaystyle\quad~{}+\frac{1}{2}\sum_{n_{i}}\Bigg{[}\Lambda^{2}_{\beta}C_{in_{i}}^{\alpha\beta}+\Lambda^{2}_{\alpha}C_{in_{i}}^{\beta\alpha}\Bigg{]}+\frac{T_{\alpha\beta}}{2}\sum_{k}\left\langle\Gamma_{ik}^{\alpha}\Gamma_{ik}^{\beta}\right\rangle$
(18)

Assuming equal parameters for all genes (which is not the same as setting
$\alpha=\beta$), this simplifies further to:

$\displaystyle C_{ii}^{\alpha\beta}$
$\displaystyle=\frac{T}{2}\frac{1-\hat{\rho}_{\rm
b}}{\tau}\sum_{\zeta}f^{\prime\zeta}(\mathbf{\bar{g}}_{i},\mathbf{c}_{i})\Big{(}C_{ii}^{\alpha\zeta}+C_{ii}^{\beta\zeta}\Big{)}$
$\displaystyle\quad~{}+\frac{\Lambda^{2}}{2}\sum_{n_{i}}\Big{(}C_{in_{i}}^{\alpha\beta}+C_{in_{i}}^{\beta\alpha}\Big{)}+\frac{T}{2}\sum_{k}\left\langle\Gamma_{ik}^{\alpha}\Gamma_{ik}^{\beta}\right\rangle$
(19)

Obviously, in this case $\hat{\rho}_{\rm on}^{\alpha}=\hat{\rho}_{\rm
on}^{\beta}\equiv 1$ because $r_{\rm on}^{\alpha}=r_{\rm on}^{\beta}\equiv
r_{\rm on}$ is the same for all genes.

####  1.3.2 Local variance

If we locally consider even the same gene, i.e. set $\alpha=\beta$ in Eq.
(18), we get the formula for the local variance $\sigma^{2}_{i,\alpha}$ of
expression level $g_{\alpha}$:

$\displaystyle\sigma^{2}_{i,\alpha}\equiv C_{ii}^{\alpha\alpha}$
$\displaystyle=T_{\alpha\alpha}\frac{\hat{\rho}_{\rm
on}^{\alpha}-\hat{\rho}_{\rm
b}^{\alpha}}{\tau_{\alpha}}\sum_{\zeta}f_{\alpha}^{\prime\zeta}(\mathbf{\bar{g}}_{i},\mathbf{c}_{i})C_{ii}^{\alpha\zeta}$
$\displaystyle\quad~{}+\Lambda^{2}_{\alpha}\sum_{n_{i}}C_{in_{i}}^{\alpha\alpha}+\frac{T_{\alpha\alpha}}{2}\sum_{k}\left\langle\big{[}\Gamma_{ik}^{\alpha}\big{]}^{2}\right\rangle$
(20)

If the parameters for all genes again are equal, we can also directly write
this as:

$\displaystyle\sigma^{2}_{i,\alpha}$ $\displaystyle=T\frac{1-\hat{\rho}_{\rm
b}}{\tau}\sum_{\zeta}f^{\prime\zeta}(\mathbf{\bar{g}}_{i},\mathbf{c}_{i})C_{ii}^{\alpha\zeta}$
$\displaystyle\quad~{}+\Lambda^{2}\sum_{n_{i}}C_{in_{i}}^{\alpha\alpha}+\frac{T}{2}\sum_{k}\left\langle\big{[}\Gamma_{ik}^{\alpha}\big{]}^{2}\right\rangle$
(21)

####  1.3.3 Inter-volume covariance for the same gene

This special case is obtained by setting $\alpha=\beta$ in Eq. (17) directly:

$\displaystyle C_{ij}^{\alpha\alpha}=$
$\displaystyle\quad\frac{T_{\alpha\alpha}}{2}\frac{\hat{\rho}_{\rm
on}^{\alpha}-\hat{\rho}_{\rm
b}^{\alpha}}{\tau}\sum_{\zeta}\Big{(}f_{\alpha}^{\prime\zeta}(\mathbf{\bar{g}}_{j},\mathbf{c}_{j})C_{ij}^{\alpha\zeta}+f_{\alpha}^{\prime\zeta}(\mathbf{\bar{g}}_{i},\mathbf{c}_{i})C_{ji}^{\alpha\zeta}\Big{)}$
$\displaystyle\quad\quad~{}+\frac{\Lambda^{2}_{\alpha}}{2}\Bigg{[}\sum_{n_{j}}C_{in_{j}}^{\alpha\alpha}+\sum_{n_{i}}C_{jn_{i}}^{\alpha\alpha}\Bigg{]}+\frac{T_{\alpha\alpha}}{2}\sum_{k}\left\langle\Gamma_{ik}^{\alpha}\Gamma_{jk}^{\alpha}\right\rangle$
(22)

This formula again reduces to the single-input/single-output formula from [1]
when setting $\forall\zeta:~{}f_{\alpha}^{\prime\zeta}=0$ (limit of no mutual
and no self-regulation).

####  1.3.4 A single self-regulating gene without diffusive coupling

As a consistency test, we can check whether the extended spatial framework
derived here correctly reproduces the formula for the variance of the self-
interacting gene without spatial coupling, derived earlier in [2].

If there is only a single species $g_{\alpha}\equiv g$ which is self-
regulating, Eq. (21) for the local variance reduces to:

$\displaystyle\sigma^{2}_{i,\alpha}$ $\displaystyle=T\frac{1-\hat{\rho}_{\rm
b}}{\tau}\left(\frac{\partial f_{\alpha}}{\partial
g_{\alpha}}\right)_{\bar{g}^{\alpha}_{i},\mathbf{c}_{i}}\underset{\sigma^{2}_{i,\alpha}}{\underbrace{C_{ii}^{\alpha\alpha}}}$
$\displaystyle\quad~{}+\Lambda^{2}\sum_{n_{i}}C_{in_{i}}^{\alpha\alpha}+\frac{T}{2}\sum_{k}\left\langle\big{[}\Gamma_{ik}^{\alpha}\big{]}^{2}\right\rangle$
(23)

For the last part of this section, we will drop the index $\alpha$. In the
limiting case without diffusion we have

$\displaystyle\Lambda^{2}$ $\displaystyle=0~{},$ $\displaystyle T$
$\displaystyle=\tau~{}.$ (24)

This leads to

$\displaystyle\sigma^{2}_{i}$ $\displaystyle=(1-\hat{\rho}_{\rm
b})\left(\frac{\partial f}{\partial
g}\right)_{\bar{g}_{i},\mathbf{c}_{i}}\sigma^{2}_{i}+\frac{\tau}{2}\sum_{k}\left\langle\Gamma_{ik}^{2}\right\rangle$
$\displaystyle\Leftrightarrow\quad\quad\sigma^{2}_{i}$
$\displaystyle=\frac{\tau}{1-(1-\hat{\rho}_{\rm b})\left(\frac{\partial
f}{\partial
g}\right)_{\bar{g}_{i},\mathbf{c}_{i}}}\frac{1}{2}\sum_{k}\left\langle\Gamma_{ik}^{2}\right\rangle$
$\displaystyle=\frac{1}{\frac{1}{\tau}-\frac{r_{\rm on}-r_{\rm b}}{N_{\rm
m}}\left(\frac{\partial f}{\partial
g}\right)_{\bar{g}_{i},\mathbf{c}_{i}}}\frac{1}{2}\sum_{k}\left\langle\Gamma_{ik}^{2}\right\rangle$
$\displaystyle\overset{(r_{\rm b}=0)}{=}\frac{1}{\frac{1}{\tau}-\frac{r_{\rm
on}}{\Omega}\left(\frac{\partial
f}{\partial\left(G/\Omega\right)}\right)_{\bar{G}_{i},\mathbf{c}_{i}}}\frac{1}{2}\sum_{k}\left\langle\Gamma_{ik}^{2}\right\rangle$
(25)

where $\Omega$ is the volume. This recovers the formula obtained from the
fluctuation-dissipation approach in [2], in which the variance is attenuated
by a “Langevin damping force”.

####  1.3.5 Short-correlations assumption

In the short-correlations assumption we assume that only nearest-neighbor
volumes have significant correlations. We thus set

$\displaystyle\forall\alpha\forall\beta:\quad C_{ij}^{\alpha\beta}$
$\displaystyle=0\quad\text{if }j\notin N(i)\cup\\{i\\}$ (26)

where $N(i)$ is the set of indices of the nearest-neighbor volumes of $i$.

In formula (17), this simplifies the diffusive coupling terms originating from
parts $\mathfrak{D}_{\alpha}$ and $\mathfrak{D}_{\beta}$. In particular, the
terms $C^{\alpha\beta}_{in_{j}}$, which denote correlations between volume $i$
and the neighbors of $j$ (!) now only have nonzero values when $i$ is a
nearest neighbor volume of $j$, i.e. for $n_{j}=i$, while all other terms
under the sum vanish. This correspondingly holds for the terms
$C^{\beta\alpha}_{jn_{i}}$, such that only $C^{\beta\alpha}_{jj}\geq 0$ when
$n_{i}=j$. We thus obtain the markedly simpler formula

$\displaystyle C_{ij}^{\alpha\beta}$
$\displaystyle=\frac{T_{\alpha\beta}}{2}\Bigg{[}\frac{\hat{\rho}_{\rm
on}^{\beta}-\hat{\rho}_{\rm
b}^{\beta}}{\tau}\sum_{\zeta}f_{\beta}^{\prime\zeta}(\mathbf{\bar{g}}_{j},c_{j})C_{ij}^{\alpha\zeta}$
$\displaystyle\quad\quad\quad\quad\quad+\frac{\hat{\rho}_{\rm
on}^{\alpha}-\hat{\rho}_{\rm
b}^{\alpha}}{\tau}\sum_{\zeta}f_{\alpha}^{\prime\zeta}(\mathbf{\bar{g}}_{i},\mathbf{c}_{i})C_{ji}^{\beta\zeta}\Bigg{]}$
$\displaystyle\quad~{}+\frac{1}{2}\Bigg{[}\Lambda^{2}_{\beta}C_{ii}^{\alpha\beta}+\Lambda^{2}_{\alpha}C_{jj}^{\beta\alpha}\Bigg{]}+\frac{T_{\alpha\beta}}{2}\sum_{k}\left\langle\Gamma_{ik}^{\alpha}\Gamma_{jk}^{\beta}\right\rangle$
(27)

which now is restricted to $j\in N(i)$, i.e. volume $j$ being a nearest
neighbor volume of volume $i$ (while $C_{ij}^{\alpha\beta}=0$ by covention if
this is not the case).

Note that in the case $i=j$ the terms in the coupling sums do _n_ ot vanish,
because then the terms $C^{\beta\alpha}_{jn_{i}}=C^{\beta\alpha}_{in_{i}}$ and
$C^{\beta\alpha}_{jn_{i}}=C^{\beta\alpha}_{in_{i}}$ all are nonzero again.
However, in this case we obtain the same formula as in Eq. (18). Eqs. (27) and
(18) together define a complete equation system for all nearest-neighbor
covariances $C_{ij}^{\alpha\beta}$ and local (inter-gene) covariances
$C_{ii}^{\alpha\beta}$ that contains significantly less unknowns than the
system for general $C_{ij}^{\alpha\beta}$ (where $j$ can also denote a non-
neighboring volume) defined by Eq. (17).

###  1.4 Noise powers

Building on our earlier work [1], our model comprises three noise sources
affecting the expression of the gap genes in the nuclei: (1.) “input noise”
from diffusive arrival of transcription factors (TFs); (2.) “output noise”
(birth-death noise) coming from the production and degradation of gap gene
products following gene activation; (3.) “diffusion noise” from the stochastic
hopping of gene products between neighboring volumes.

The diffusion noise not only contributes to the summed noise sources affecting
the copy number of gene $\alpha$ in volume $i$,
$\Big{\langle}\left[\mathcal{N}^{\alpha}_{ii}\right]^{2}\Big{\rangle}$; it
also induces anticorrelated fluctuations in neighboring volumes that we have
to include in the nearest-neighbor noise correlations
$\Big{\langle}\left[\mathcal{N}^{\alpha}_{i\nu}\right]^{2}\Big{\rangle}\equiv\Big{\langle}\Gamma_{ik}^{\alpha}\Gamma_{\nu
k}^{\alpha}\Big{\rangle}=-\Big{\langle}\left[\mathcal{D}^{\alpha}_{[i\rightarrow\nu]}\right]^{2}\Big{\rangle}-\Big{\langle}\left[\mathcal{D}^{\alpha}_{[\nu\rightarrow
i]}\right]^{2}\Big{\rangle}$, where $\nu$ indicates one of the nearest
neighbor volumes of $i$. Here and throughout we will assume that diffusive
processes of different protein species are uncorrelated, i.e. the
corresponding term $\Big{\langle}\Gamma_{ik}^{\alpha}\Gamma_{\nu
k}^{\beta}\Big{\rangle}\equiv 0$. We calculate the detailed form of the
diffusion noise contributions further below.

We will now first compute the normalized “local” noise term
$\Big{\langle}\left[\mathcal{N}^{\alpha}_{ii}\right]^{2}\Big{\rangle}$, which
in more detail we can write as follows:

$\displaystyle\left\langle\left[\mathcal{N}^{\alpha}_{ii}\right]^{2}\right\rangle$
$\displaystyle\equiv\underset{\text{input
noise}}{\underbrace{\left\langle\left[\mathcal{J}^{\alpha}_{i}\right]^{2}\right\rangle}}+\underset{\text{output
noise}}{\underbrace{\left\langle\left[\mathcal{O}^{\alpha}_{i}\right]^{2}\right\rangle}}+\underset{\text{diffusion/hopping
noise}}{\underbrace{\sum_{\nu}\left\langle\left[\mathcal{D}^{\alpha}_{[i\rightarrow\nu]}\right]^{2}+\left[\mathcal{D}^{\alpha}_{[\nu\rightarrow
i]}\right]^{2}\right\rangle}}$ (28)

Let us first consider the input noise sources. Here we assume that the
combined diffusive search and binding processes responsible for this noise
will be uncorrelated among the different genes. This “crosstalk-free”
situation is expected to hold when the copy numbers of all transcription
factors are comparably low, and the accessibility of a TF-specific binding
site does not depend on the presence of a different TF type. Then we can
simply add up the input noise contributions contributed by all the different
TFs that regulate gene $\alpha$:

$\displaystyle\left\langle\left[\mathcal{J}^{\alpha}_{i}\right]^{2}\right\rangle=\underset{\equiv\sum_{\kappa}\mathcal{J}_{c}^{\kappa}}{\underbrace{\frac{1}{N_{\rm
m}^{2}}(r_{\rm on}^{\alpha}-r_{\rm
b}^{\alpha})^{2}\sum_{\kappa}\left[\left(\frac{\partial f_{\alpha}}{\partial
c^{\kappa}}\right)^{2}\frac{2c^{\kappa}~{}\Phi_{\kappa}^{\alpha}}{D_{\kappa}l_{\alpha}}\right]_{c_{i}^{\kappa}}}}$
$\displaystyle+\underset{\equiv\sum_{\zeta}\mathcal{J}_{g}^{\zeta}}{\underbrace{\frac{1}{N_{\rm
m}^{2}}(r_{\rm on}^{\alpha}-r_{\rm
b}^{\alpha})^{2}\sum_{\zeta}\left[\left(\frac{\partial
f_{\alpha}}{\partial(G^{\zeta}/\Omega)}\right)^{2}\frac{2(G^{\zeta}/\Omega)\Phi_{\zeta}^{\alpha}}{D_{\zeta}l_{\alpha}}\right]_{\left\langle
G_{i}^{\zeta}\right\rangle}}}$ (29)

where $G_{j}^{\zeta}$ are absolute copy numbers of the zygotic species,
$\Omega$ the reaction (nuclear) volume, and $c_{i}^{\kappa}$ the concentration
of maternal input species $\kappa$, and $l_{\alpha}$ the spatial size of the
TF binding sites, which we assume to be the same for all TF species and simply
denote by $l$. $D_{\kappa}$ and $D_{\zeta}$ denote the internal diffusion
constants of the maternal and zygotic TF species, which we assume to be equal
as well for all TF species, denoting it by $D_{\text{int}}$ from now on.
$\Phi_{\kappa}^{\alpha}$ and $\Phi_{\zeta}^{\alpha}$ are functions of the
internal structure of the regulatory region of gene $\alpha$, which may depend
on the average occupancies. While throughout this work we set all
$\Phi_{\kappa(\zeta)}^{\alpha}\equiv 1$, we included these factors in the
derivations for potential future applications of the model.

The first terms for the input noise from maternal inputs can (for a particular
$\kappa$) be simplified further as

$\displaystyle\mathcal{J}_{c}^{\kappa}$ $\displaystyle=\frac{1}{N_{\rm
m}\tau}\frac{(r_{\rm on}^{\alpha}-r_{\rm b}^{\alpha})^{2}}{r_{\rm
on}^{2}}\left[\left(\frac{\partial f_{\alpha}}{\partial
c^{\kappa}}\right)^{2}\frac{2c^{\kappa}N_{\rm
m}}{D_{\text{int}}l\tau}\Phi_{\kappa}^{\alpha}\right]_{c_{i}^{\kappa}}$
$\displaystyle=\frac{1}{N_{\rm m}\tau}(\hat{\rho}_{\rm
on}^{\alpha}-\hat{\rho}_{\rm b}^{\alpha})^{2}\left[\left(\frac{\partial
f_{\alpha}}{\partial
c^{\kappa}}\right)^{2}2c^{\kappa}c_{0}\Phi_{\kappa}^{\alpha}\right]_{c_{i}^{\kappa}}$

where $c_{0}\equiv N_{\rm m}/D_{\text{int}}l\tau$ defines a natural
concentration scale, as in our previous work [3, 4, 2, 1].

For the other ($G_{i}^{\zeta}$-dependent) terms, we can similarly write, for
each $\zeta$ separately:

$\displaystyle\mathcal{J}_{g}^{\zeta}$ $\displaystyle=\frac{1}{N_{\rm
m}^{2}}(r_{\rm on}^{\alpha}-r_{\rm b}^{\alpha})^{2}\left[\left(\frac{\partial
f_{\alpha}}{\partial\left(\frac{G^{\zeta}}{N_{\rm
m}}\right)}\right)^{2}\frac{\Omega}{N_{\rm
m}}\frac{2\left(\frac{G^{\zeta}}{N_{\rm
m}}\right)\Phi_{\zeta}^{\alpha}}{D_{\text{int}}l}\right]_{\left\langle
G_{i}^{\zeta}\right\rangle}$ $\displaystyle=\frac{1}{N_{\rm
m}\tau}\frac{(r_{\rm on}^{\alpha}-r_{\rm b}^{\alpha})^{2}}{r_{\rm
on}^{2}}\left[\left(\frac{\partial f_{\alpha}}{\partial
g^{\zeta}}\right)^{2}\frac{\Omega}{N_{\rm m}}\frac{2g^{\zeta}N_{\rm
m}}{D_{\text{int}}l\tau}\Phi_{\zeta}^{\alpha}\right]_{\bar{g}_{i}^{\zeta}}$
$\displaystyle=\frac{1}{N_{\rm m}\tau}\frac{1}{\gamma_{\rm m}}(\hat{\rho}_{\rm
on}^{\alpha}-\hat{\rho}_{\rm b}^{\alpha})^{2}\left[\left(\frac{\partial
f_{\alpha}}{\partial
g^{\zeta}}\right)^{2}2g^{\zeta}c_{0}\Phi_{\zeta}^{\alpha}\right]_{\bar{g}_{i}^{\zeta}}$

where $\gamma_{\rm m}\equiv N_{\rm m}/\Omega$ is the maximal mean
concentration.

The output noise is simple birth-death shot noise, and given by the sum of the
average birth and death rates:

$\displaystyle\left\langle\left[\mathcal{O}^{\alpha}_{i}\right]^{2}\right\rangle$
$\displaystyle=\frac{1}{N_{\rm m}^{2}}\left[\left(r_{\rm b}^{\alpha}+(r_{\rm
on}^{\alpha}-r_{\rm
b}^{\alpha})f_{\alpha}(\mathbf{\bar{g}}_{i},\mathbf{c}_{i})\right)+\frac{1}{\tau_{\alpha}}\left\langle
G_{i}^{\alpha}\right\rangle\right]$ $\displaystyle=\frac{1}{N_{\rm
m}}\left[\frac{1}{\tau}\left(\hat{r}_{\rm b}^{\alpha}+(\hat{\rho}_{\rm
on}^{\alpha}-\hat{r}_{\rm
b}^{\alpha})f_{\alpha}(\mathbf{\bar{g}}_{i},\mathbf{c}_{i})\right)+\frac{1}{\tau_{\alpha}}\bar{g}_{i}^{\alpha}\right]$

Finally, for the diffusion noise sources we have:

$\displaystyle\left\langle\left[\mathcal{D}^{\alpha}_{[i\rightarrow\nu]}\right]^{2}\right\rangle$
$\displaystyle=\frac{1}{N_{\rm m}^{2}}h_{\alpha}\left\langle
G_{i}^{\alpha}\right\rangle=\frac{1}{N_{\rm m}}h_{\alpha}\bar{g}_{i}^{\alpha}$
$\displaystyle\left\langle\left[\mathcal{D}^{\alpha}_{[\nu\rightarrow
i]}\right]^{2}\right\rangle$ $\displaystyle=\frac{1}{N_{\rm
m}^{2}}h_{\alpha}\left\langle G_{\nu}^{\alpha}\right\rangle=\frac{1}{N_{\rm
m}}h_{\alpha}\bar{g}_{\nu}^{\alpha}$ (33)

These terms not only appear in
$\Big{\langle}\left[\mathcal{N}^{\alpha}_{ii}\right]^{2}\Big{\rangle}$; in
their negated form, they also constitute the nearest-neighbor noise covariance
$\Big{\langle}\left[\mathcal{N}^{\alpha}_{i\nu}\right]^{2}\Big{\rangle}$,
listed in the final summary of all noise contributions below.

Taken together, assuming equal internal TF diffusion constants and binding
size lengths, the noise terms read

$\displaystyle\left\langle\left[\mathcal{N}^{\alpha}_{ii}\right]^{2}\right\rangle=\frac{1}{N_{\rm
m}}\Bigg{\\{}$ $\displaystyle\quad\quad\frac{1}{\tau}\left[\hat{r}_{\rm
b}^{\alpha}+(\hat{\rho}_{\rm on}^{\alpha}-\hat{r}_{\rm
b}^{\alpha})f_{\alpha}(\mathbf{\bar{g}}_{i},\mathbf{c}_{i})\right]+\frac{1}{\tau_{\alpha}}\bar{g}_{i}^{\alpha}$
“output noise” $\displaystyle\quad+\frac{(\hat{\rho}_{\rm
on}^{\alpha}-\hat{\rho}_{\rm
b}^{\alpha})^{2}}{\tau}\sum_{\kappa}\left[\left(\frac{\partial
f_{\alpha}}{\partial
c^{\kappa}}\right)^{2}2c^{\kappa}c_{0}\Phi_{\kappa}^{\alpha}\right]_{c_{i}^{\kappa}}$
“maternal input noise” $\displaystyle\quad+\frac{1}{\gamma_{\rm
m}}\frac{(\hat{\rho}_{\rm on}^{\alpha}-\hat{\rho}_{\rm
b}^{\alpha})^{2}}{\tau}\sum_{\zeta}\left[\left(\frac{\partial
f_{\alpha}}{\partial
g^{\zeta}}\right)^{2}2g^{\zeta}c_{0}\Phi_{\zeta}^{\alpha}\right]_{\bar{g}_{i}^{\zeta}}$
“zygotic input noise” $\displaystyle\quad+2d\cdot
h_{\alpha}\bar{g}_{i}^{\alpha}+\sum_{\nu}h_{\alpha}\bar{g}_{\nu}^{\alpha}~{}\Bigg{\\}}$
“diffusion noise” (34)
$\displaystyle\left\langle\left[\mathcal{N}^{\alpha}_{i\nu}\right]^{2}\right\rangle=-\frac{h_{\alpha}}{N_{\rm
m}}\left(\bar{g}_{i}^{\alpha}+\bar{g}_{\nu}^{\alpha}\right)$ “diffusion noise”
(35)

where $\kappa$ and $\zeta$, respectively, run over all maternal input and gap
gene output species in the system, and $\nu\in N(i)$ over all nearest
neighbors of volume $i$.

###  1.5 Regulatory function

Until now we could carry out all calculations without specifying the nonlinear
regulatory function $f_{\alpha}(\mathbf{g}_{i},\mathbf{c}_{i})$ that converts
the local maternal input and gap protein levels into an average activity level
of the regulated gap gene. While multiple choices are possible here, such as
Hill-type regulatory functions, in this work we opted for the Monod-Wyman-
Changeux (MWC) regulation model, which allows exploration of a richer set of
solutions. Our parametrization of the model is based on the approximation used
in [4], which we briefly recapitulate here once more for completeness.

The MWC model was originally introduced for allosteric enzymes that can bind
ligands with different affinities in two different conformational states named
$T$ (“tense”) and $R$ (“relaxed”). In the absence of any ligand, the occupancy
of the respective states is determined by a single equilibrium constant
$L=[T]/[R]$. The ligands can bind to the enzyme in both states, but with
different affinities, characterized by different dissociation constants
$K_{T}$ and $K_{R}$. Since the affinity of the ligand to binding sites in the
$R$ state is higher, the equilibrium is shifted towards the $R$ state upon
ligand binding. One notable limitation of the MWC model is that it requires
all binding subunits to be in the same state; as such, the whole enzyme
switches state at once, partial switching of subunits is not accounted for.
With these assumptions, formulas for two quantities can be derived: the (mean)
fractional occupancy of (all) ligand binding sites, $\bar{Y}$, and the
fraction of enzymes in the $R$ state, $\bar{R}$:

$\displaystyle\bar{Y}$
$\displaystyle=\frac{\xi(1+\xi)^{n-1}+LC\xi(1+C\xi)^{n-1}}{(1+\xi)^{n}+L(1+C\xi)^{n}}$
(36) $\displaystyle\bar{R}$
$\displaystyle=\frac{(1+\xi)^{n}}{(1+\xi)^{n}+L(1+C\xi)^{n}}$ (37)

Here $C=K_{R}/K_{T}$ is the state affinity ratio, $n$ is the number of
(identical) subunits that can bind the ligand, and $\xi=[X]/K_{R}$ is the
ligand concentration $[X]$ normalized to the affinity of the $R$ state.

This model can be employed as a phenomenological model of gene regulation, as
in previous published works. This model is phenomenological in that the gene
regulation is likely not a process at thermal equilibrium; nevertheless, the
model provides a convenient and versatile parametrization of the regulatory
functions, as we will show in what follows. Applied to gene regulation, a gene
can exist in “off” or “on” expression states (mimicking the balance of the $T$
and $R$ states in the original model application to haemoglobin), and the
binding of transcription factors (TFs) (which take the role of ligands in the
original model) shifts the effective balance between the two gene expression
states. For a specific gene $\alpha$, let us denote the constant describing
the ratio between the “on” and “off” states by $L_{\alpha}$, the TF
concentration by $c$, and the dissociation constants for TF binding in the
“on” and “off” states by $K_{\rm on}^{\alpha 0}$ and $K_{\rm off}^{\alpha 0}$,
respectively; the indices $\alpha 0$ are used to indicate that there is only
one regulator protein species, which we extend to the case of multiple
regulators further below. Then we can rewrite Eq. (37), which now will
describe the activity $f_{\alpha}(c)$ of the gene as a function of the TF
concentration $c$, as follows:

$\displaystyle f_{\alpha}(c)$ $\displaystyle=\frac{\left(1+\frac{c}{K_{\rm
on}^{\alpha 0}}\right)^{H^{\alpha 0}_{\rm
reg}}}{L_{\alpha}\left(1+\frac{c}{K_{\rm off}^{\alpha 0}}\right)^{H^{\alpha
0}_{\rm reg}}+\left(1+\frac{c}{K_{\rm on}^{\alpha 0}}\right)^{H^{\alpha
0}_{\rm reg}}}$ (38)

Here the parameter $n$ (which could be interpreted as the number of binding
sites for the TF) determines the strength of regulatory action of input $c$ on
species $\alpha$, denoted by $H^{\alpha 0}_{\rm reg}$. Notice that the
constant $C$ in Eq. (37) now is absorbed into the different denominators
dividing the TF concentration $c$. For further clarity, we can also rewrite
the formula as follows

$\displaystyle f_{\alpha}(c)$ $\displaystyle=\frac{e^{-\beta\varepsilon_{\rm
on}}\left(1+\frac{c}{K_{\rm on}^{\alpha 0}}\right)^{H^{\alpha 0}_{\rm
reg}}}{e^{-\beta\varepsilon_{\rm off}}\left(1+\frac{c}{K_{\rm off}^{\alpha
0}}\right)^{H^{\alpha 0}_{\rm reg}}+e^{-\beta\varepsilon_{\rm
on}}\left(1+\frac{c}{K_{\rm on}^{\alpha 0}}\right)^{H^{\alpha 0}_{\rm reg}}}$
(39)

where $L_{\alpha}=e^{\beta\Delta\varepsilon}\equiv e^{\beta(\varepsilon_{\rm
on}-\varepsilon_{\rm off})}$ now corresponds to the ratio of the Boltzmann-
like weights of the (unregulated/unbound) “on” and “off” states. As it will
become clear later, the difference in the exponential corresponds to an
effective “base energy” appearing in the denominator of $f_{\alpha}(c)$,
controlling to the probability of spontaneous (“basal”) activity.

In reference [4] the constant $L_{\alpha}$ is written as
$L_{\alpha}=\tilde{L}_{\alpha}^{H^{\alpha 0}_{\rm reg}}$, which makes it
easier to combine it with the other terms. Let us now consider the limit in
which spontaneous activity is improbable (compared to regulated activity),
i.e. $\varepsilon_{\rm on}\gg\varepsilon_{\rm off}$ or $L_{\alpha}\gg 1$, and
in which binding to the off-state is improbable, too, $K_{\rm off}^{\alpha
0}\gg K_{\rm on}^{\alpha 0}$. Then we can ignore the term describing the
probability of the bound-but-inactive state and write:

$\displaystyle f_{\alpha}(c)$
$\displaystyle\simeq\frac{\left(1+\frac{c}{K_{\rm on}^{\alpha
0}}\right)^{H^{\alpha 0}_{\rm reg}}}{\tilde{L}_{\alpha}^{H^{\alpha 0}_{\rm
reg}}+\left(1+\frac{c}{K_{\rm on}^{\alpha 0}}\right)^{H^{\alpha 0}_{\rm
reg}}}=\frac{1}{1+\tilde{L}_{\alpha}^{H^{\alpha 0}_{\rm
reg}}\left/\left(1+\frac{c}{K_{\rm on}^{\alpha 0}}\right)^{H^{\alpha 0}_{\rm
reg}}\right.}$ $\displaystyle=\frac{1}{1+\exp\left[-H^{\alpha 0}_{\rm
reg}\ln\left(\frac{1+c/K_{\rm on}^{\alpha 0}}{1+\tilde{c}/K_{\rm on}^{\alpha
0}}\right)\right]}=\frac{1}{1+\exp\left[-H^{\alpha 0}_{\rm
reg}\ln\left(1+\frac{c}{K_{\rm on}^{\alpha
0}}\right)+F_{\alpha}^{(0)}\right]}$ (40)

Here, $\tilde{c}\equiv(\tilde{L}_{\alpha}-1)K_{\rm on}^{\alpha 0}$ has the
interpretation of an effective half-activation (half-regulation) threshold
concentration, while we see that $F_{\alpha}^{(0)}=\ln L_{\alpha}$ indeed
corresponds to a basal free-energy-like term.

We can straightforwardly generalize this to the multi-gene case. When gene
$\alpha$ is regulated by multiple TF concentrations $g^{\zeta}$ we can
initially write

$\displaystyle f_{\alpha}(\\{g^{\zeta}\\})$
$\displaystyle=\frac{\prod_{\zeta}\left(1+\frac{g^{\zeta}}{K_{\rm
on}^{\alpha\zeta}}\right)^{H^{\alpha\zeta}_{\rm
reg}}}{L~{}\prod_{\zeta}\left(1+\frac{g^{\zeta}}{K_{\rm
off}^{\alpha\zeta}}\right)^{H^{\alpha\zeta}_{\rm
reg}}+\prod_{\zeta}\left(1+\frac{g^{\zeta}}{K_{\rm
on}^{\alpha\zeta}}\right)^{H^{\alpha\zeta}_{\rm reg}}}$ (41)

where $\zeta$ now loops over all regulatory input species. Performing the same
steps as before for each regulating input separately and regrouping yields

$\displaystyle
f_{\alpha}(\\{g^{\zeta}\\})=\frac{1}{1+e^{-F_{\alpha}(\\{g^{\zeta}\\})}}~{},$
(42)

where the free-energy-like term $F_{\alpha}$ now reads

$\displaystyle F_{\alpha}(\\{g^{\zeta}\\})=\sum_{\zeta}H^{\alpha\zeta}_{\rm
reg}\ln\left(1+g^{\zeta}/K_{\rm on}^{\alpha\zeta}\right)-F^{(0)}_{\alpha}$
(43)

with $F^{(0)}_{\alpha}=\ln L$ again.

While the $F$ term above groups all regulatory interactions (including self-
regulation) under one sum, in our model we want to distinguish maternal
regulation from the zygotic interactions, and therefore use the following more
specific notation

$\displaystyle f_{\alpha}(\mathbf{g},\mathbf{c})$
$\displaystyle=\frac{1}{1+e^{-F_{\alpha}(\mathbf{g},\mathbf{c})}}$

with

$\displaystyle F_{\alpha}(\mathbf{g},\mathbf{c})$
$\displaystyle=\sum_{\zeta}H_{G}^{\alpha\zeta}\ln\left(1+g^{\zeta}/K_{G}^{\alpha\zeta}\right)+\sum_{\kappa}H_{M}^{\alpha\kappa}\ln\left(1+c^{\kappa}/K_{M}^{\alpha\kappa}\right)-F^{(0)}_{\alpha}$
(45)

where $\kappa$ loops over all maternal and $\zeta$ over all zygotic inputs,
respectively. In our standard model with 3 maternal inputs and 4 gap genes
expressed downstream of them, $\zeta\in\\{1,2,3,4\\}$ and
$\kappa\in\\{A,P,T\\}$, where the letters stand for the anterior (A),
posterior (P) and terminal (T) maternal systems.

###  1.6 Combined mRNA and protein model

In our original generic stochastic differential equation (SDE) ansatz in Sec.
1.1 we did not take into account the intermediate step of mRNA production.
However, it is known that the translation process can induce additional noise
in the output, in particular when protein production occurs in bursts. In this
section we will extend our original framework to effectively incorporate
bursty protein production by adapting the original set of SDEs and the noise
covariance terms.

####  1.6.1 Incorporating mRNA into the model

Let us start by extending the set of differential equations in Sec. 1.1 in
order to incorporate the mRNA populations, which we will call
$M_{i}^{\alpha}$. In the following we will distinguish mRNA rates and
lifetimes from protein rates and lifetimes by the tilde, i.e. whenever
$\tau_{\alpha}$ will refer to the protein lifetime, $\tilde{\tau}_{\alpha}$
will denote the lifetime of the corresponding mRNA. The dynamics for a mRNA-
protein population pair for gene species $\alpha$ then can be written as
(abbreviating noise terms as $\mathcal{N}_{M}$ and $\mathcal{N}_{G}$,
respectively):

$\displaystyle\partial_{t}M_{i}^{\alpha}$
$\displaystyle=r_{\alpha}(\mathbf{G}_{i},\mathbf{c}_{i})-\frac{1}{\tilde{\tau}_{\alpha}}M_{i}^{\alpha}+\mathcal{N}_{M}$
(46) $\displaystyle\partial_{t}G_{i}^{\alpha}$
$\displaystyle=k_{\alpha}M_{i}^{\alpha}-\frac{1}{\tau_{\alpha}}G_{i}^{\alpha}-h_{\alpha}\sum_{\nu}\left(G_{i}^{\alpha}-G_{\nu}^{\alpha}\right)+\mathcal{N}_{G}$
(47)

Herein, $h_{\alpha}=D_{\alpha}/\Delta^{2}$ is the protein hopping rate between
neighboring volumes at distance $\Delta$ and $\nu$ runs over all nearest-
neighbor volumes, as before. These equations are still unnormalized, i.e.
given in absolute numbers. Note that the equation for the protein levels
$G_{i}^{\alpha}$ closely follows the original ansatz.

Let us first look at the means of these equations, denoted by bars over the
stochastic quantities:

$\displaystyle\bar{M}_{i}^{\alpha}$
$\displaystyle=\tilde{\tau}_{\alpha}~{}r_{\alpha}(\bar{\mathbf{G}}_{i},\mathbf{c}_{i})$
(48) $\displaystyle\bar{G}_{i}^{\alpha}$
$\displaystyle=\tau_{\alpha}k_{\alpha}\bar{M}_{i}^{\alpha}-\tau_{\alpha}h_{\alpha}\sum_{\nu}\left(\bar{G}_{i}^{\alpha}-\bar{G}_{\nu}^{\alpha}\right)$
(49)

The second equation illustrates that we cannot simply define the burst size as
$\beta_{\alpha}$ as a ratio $\bar{G}_{i}^{\alpha}/\bar{M}_{i}^{\alpha}$,
because $\bar{G}_{i}^{\alpha}$ is also affected by diffusion. We therefore opt
for defining the burst size in our model as:

$\displaystyle\beta_{\alpha}\equiv k_{\alpha}\tau_{\alpha}$ (50)

####  1.6.2 Normalized differential equation for the means of mRNA and protein

In order to normalize the SDEs introduced in the previous subsection let us
define

$\displaystyle M_{\rm max}$ $\displaystyle\equiv r_{\rm on}\tilde{\tau}$
$\displaystyle G_{\rm max}$ $\displaystyle\equiv k_{\rm max}\tau M_{\rm
max}=k_{\rm max}r_{\rm on}\tau\tilde{\tau}$ (51)

where $r_{\rm on}$ is the fastest mRNA transcription rate (among all species
$\alpha$), $k_{\rm max}$ is the (fastest) protein translation rate from mRNA,
$\tilde{\tau}$ is the longest mRNA lifetime, and $\tau$ the longest protein
lifetime of all the species. This means that $M_{\rm max}$ and $G_{\rm max}$
are the highest possible (per-species) mRNA and protein numbers in the system.
These rates also imply a maximal burst size $\beta_{\rm m}\equiv k_{\rm
max}\tau$; while the species-specific burst sizes $\beta_{\alpha}$ in
principle may be different, they will all be smaller than or equal to
$\beta_{\rm m}$, and the corresponding ratios $\beta_{\alpha}/\beta_{\rm m}$
will appear in the normalized noise covariances further below.

Using the maximal protein copy number $G_{\max}$ for normalization, $N_{\rm
m}\equiv G_{\max}=k_{\rm max}r_{\rm on}\tau\tilde{\tau}$, and averaging over
Eqs. (46) and (47) , we directly obtain the set of coupled ordinary
differential equations describing the time-evolution of the normalized mean
expression levels of mRNA and protein populations in our extended model:

$\displaystyle\partial_{t}\bar{m}_{i}^{\alpha}$
$\displaystyle=\frac{r_{\alpha}(\bar{\mathbf{g}}_{i},\mathbf{c}_{i})}{N_{\rm
m}}-\frac{1}{\tilde{\tau}_{\alpha}}\bar{m}_{i}^{\alpha}$ (52)
$\displaystyle\partial_{t}\bar{g}_{i}^{\alpha}$
$\displaystyle=\frac{\beta_{\alpha}}{\tau_{\alpha}}\bar{m}_{i}^{\alpha}-\frac{1}{\tau_{\alpha}}\bar{g}_{i}^{\alpha}-h_{\alpha}\sum_{\nu}\left(\bar{g}_{i}^{\alpha}-\bar{g}_{\nu}^{\alpha}\right)$
(53)

where $\bar{m}_{i}^{\alpha}\equiv\langle M_{i}^{\alpha}\rangle/N_{\rm m}$ and
$\bar{g}_{i}^{\alpha}\equiv\langle G_{i}^{\alpha}\rangle/N_{\rm m}$.

By further carrying out the normalization of the production function
$r_{\alpha}(\mathbf{g}_{i},\mathbf{c}_{i})$ we recover almost the same
structure as in the original model:

$\displaystyle\frac{r_{\alpha}(\mathbf{g}_{i},\mathbf{c}_{i})}{N_{\rm m}}$
$\displaystyle=\frac{1}{\beta_{\rm m}r_{\rm on}\tilde{\tau}}\left[r_{\rm
b}^{\alpha}+(r_{\rm on}^{\alpha}-r_{\rm
b}^{\alpha})f^{\alpha}(\mathbf{g}_{i},\mathbf{c}_{i})\right]$
$\displaystyle=\frac{1}{\beta_{\rm m}\tilde{\tau}}\left[\hat{\rho}_{\rm
b}^{\alpha}+(\hat{\rho}_{\rm on}^{\alpha}-\hat{\rho}_{\rm
b}^{\alpha})f^{\alpha}(\mathbf{g}_{i},\mathbf{c}_{i})\right]\equiv\hat{r}(\mathbf{g}_{i},\mathbf{c}_{i})$
(54)

Note, however, that now the normalized rates $\hat{\rho}_{\rm b}^{\alpha}$ and
$\hat{\rho}_{\rm on}^{\alpha}$ have a slightly different meaning, as they
refer to the mRNA basal production rate and the mRNA production rate at full
induction, respectively; the additional prefactor $1/\beta_{\rm m}$ accounts
for our choice of normalizing by the maximal protein level.

Eqs. (46) and (47) allow us to forward-integrate all coupled mRNA and protein
levels in the system in order to obtain predictions of the mean expression
levels at any desired time point. Obtaining a similar prediction for the
covariance matrix of the coupled system is technically more challenging, but
simplifies considerably when the cases of significantly differing mRNA and
protein time scales are considered, as described in the next subsection.

####  1.6.3 Approximations for computing noise covariances

In order to avoid the necessity of repeating the derivations for the
covariance matrix with explicit inclusion of the additionally arising
covariances of the kind $\left\langle\delta G_{i}^{\alpha}\delta
M_{j}^{\beta}\right\rangle$, we consider the extended model in the regime of
separated mRNA and protein time scales, which allows us to reuse all previous
derivations.

If the mRNA dynamics were much faster than the protein dynamics, mRNA
fluctuations would be largely averaged out at the protein level and the mRNA
levels would be mainly determined by the protein levels propagated through the
regulatory function. We could then approximate

$\displaystyle
M_{i}^{\alpha}\simeq\tilde{\tau}_{\alpha}r_{\alpha}(\mathbf{G}_{i},\mathbf{c}_{i})$
(55)

such that

$\displaystyle\partial_{t}G_{i}^{\alpha}$ $\displaystyle\simeq
k_{\alpha}\tilde{\tau}_{\alpha}~{}r_{\alpha}(\mathbf{G}_{i},\mathbf{c}_{i})-\frac{1}{\tau_{\alpha}}G_{i}^{\alpha}-h_{\alpha}\sum_{\nu}\left(G_{i}^{\alpha}-G_{\nu}^{\alpha}\right)+\mathcal{N}_{G}$
(56)

which recovers the original set of equations without explicit treatment of
mRNA in Sec. 1.1.

In the opposite regime, i.e. $\tilde{\tau}_{\alpha}\gg\tau_{\alpha}$, which we
consider in our model, we can make another approximation assuming that the
mRNA is the relevant species and that the protein populations equilibrate on a
timescale faster than the mRNA populations are changing. In this limit, the
protein populations effectively “track” the mRNA populations and we can make
the ansatz:

$\displaystyle\tilde{G}_{i}^{\alpha}\equiv\beta_{\alpha}M_{i}^{\alpha}\equiv
k_{\alpha}\tau_{\alpha}M_{i}^{\alpha}$ (57)

This allows us to write down an effective SDE for the production and
degradation of $\tilde{G}_{i}^{\alpha}$ by simply multiplying the
corresponding SDE for the mRNA population by the protein burst size
$\beta_{\alpha}$:

$\displaystyle\partial_{t}\tilde{G}_{i}^{\alpha}$
$\displaystyle=\beta_{\alpha}r_{\alpha}(\mathbf{\tilde{G}}_{i},\mathbf{c}_{i})-\frac{1}{\tilde{\tau}_{\alpha}}\tilde{G}_{i}^{\alpha}+\beta_{\alpha}\mathcal{N}_{\tilde{G}}$
(58)

The equation above does not incorporate protein diffusion yet. But since
$\tilde{G}$ is a very good proxy for the protein population $G$ we can assume
that the $\tilde{G}$ population is altered by diffusion in the same way as the
$G$ population, and that the noise contributions originating from diffusive
exchange, which we denote by $\mathcal{D}_{\tilde{G}}$, can be simply added on
top of the other noise sources that affect $\tilde{G}$ (and actually propagate
from the mRNA to the protein levels). We can therefore incorporate diffusion
into our ansatz as follows:

$\displaystyle\partial_{t}\tilde{G}_{i}^{\alpha}$
$\displaystyle=\beta_{\alpha}r_{\alpha}(\mathbf{\tilde{G}}_{i},\mathbf{c}_{i})-\frac{1}{\tilde{\tau}_{\alpha}}\tilde{G}_{i}^{\alpha}-h_{\alpha}\sum_{\nu}\left(\tilde{G}_{i}^{\alpha}-\tilde{G}_{\nu}^{\alpha}\right)+\beta_{\alpha}\mathcal{N}_{\tilde{G}}+\mathcal{D}_{\tilde{G}}$
(59)

This equation again recovers the structure of our ansatz in Sec. 1.1. The
relevant differences are the burst-size prefactors in the production and noise
term; while the former will disappear after normalization, the second
prefactor will lead to an additional noise contribution, as shown further
below. Also note that the relevant integration time scale in the equation
above is the mRNA life time $\tilde{\tau}_{\alpha}$.

We can now again normalize the approximated SDE above by the normalization
constant $N_{\rm m}=G_{\rm max}$ and defining
$\tilde{g}_{i}\equiv\tilde{G}_{i}/N_{\rm m}$:

$\displaystyle\partial_{t}\tilde{g}_{i}^{\alpha}$
$\displaystyle\simeq\beta_{\alpha}~{}\hat{r}_{\alpha}(\mathbf{\tilde{g}}_{i},\mathbf{c}_{i})-\frac{1}{\tilde{\tau}_{\alpha}}\tilde{g}_{i}^{\alpha}-h_{\alpha}\sum_{\nu}\left(\tilde{g}_{i}^{\alpha}-\tilde{g}_{\nu}^{\alpha}\right)+\frac{\beta_{\alpha}}{N_{\rm
m}}\mathcal{N}_{\tilde{G}}+\frac{1}{N_{\rm m}}\mathcal{D}_{\tilde{G}}$ (60)

The equation above implies that in the considered limit we can approximate the
noise powers affecting the protein levels by the noise powers affecting the
(absolute) mRNA levels multiplied by prefactor $\beta_{\alpha}/N_{\rm m}$ and
supplemented by the protein-level (normalized) diffusion noise
$\mathcal{D}_{\tilde{G}}/N_{\rm m}$, which has the same structure as in our
original model since we normalize by the maximal protein copy number as
before. We calculate the explicit normalized forms of all the modified noise
covariances comprised in the extended mRNA-protein model in the subsequent
section.

####  1.6.4 Normalized noise covariances with mRNA

As in the simpler model without separate treatment of mRNA and protein, for
which the noise covariances are treated in Sec. 1.4, the noise term
$\frac{\beta_{\alpha}}{N_{\rm m}}\mathcal{N}_{\tilde{G}}$ in Eq. (60) consists
of three contributions that propagate into or arise at the mRNA level: (1) the
input noise from maternal regulation ($\mathcal{J}_{c}$), (2) the input noise
from regulation by gap proteins acting as transcription factors to its own or
other gap genes’ regulatory region ($\mathcal{J}_{g}$), and (3) the output
noise originating from production and degradation, denoted by ($\mathcal{O}$).
The covariance terms describing these noise contributions are the same ones
already used in Sec. 1.4, only that we have to work with the rates and
lifetimes relevant to the mRNA populations now, and account for the additional
prefactor $\frac{\beta_{\alpha}}{N_{\rm m}}$ by multiplying with its square.
For the maternal input noise contributions (cf. Eq.
(LABEL:eqInputNoiseMaternal) in the model without mRNA) this yields (using
$N_{\rm m}=\beta_{\rm m}r_{\rm on}\tilde{\tau}$):

$\displaystyle\mathcal{J}_{c}^{\kappa}$
$\displaystyle=\frac{\beta_{\alpha}^{2}}{N_{\rm m}^{2}}(r_{\rm
on}^{\alpha}-r_{\rm b}^{\alpha})^{2}\left[\left(\frac{\partial
f_{\alpha}}{\partial
c^{\kappa}}\right)^{2}\frac{2c^{\kappa}}{D_{\text{int}}l}\Phi_{\kappa}^{\alpha}\right]_{c_{i}^{\kappa}}$
$\displaystyle=\frac{\beta_{\alpha}^{2}}{N_{\rm m}\tilde{\tau}}\frac{(r_{\rm
on}^{\alpha}-r_{\rm b}^{\alpha})^{2}}{\beta_{\rm m}^{2}r_{\rm
on}^{2}}\left[\left(\frac{\partial f_{\alpha}}{\partial
c^{\kappa}}\right)^{2}\frac{2c^{\kappa}N_{\rm
m}}{D_{\text{int}}l\tilde{\tau}}\Phi_{\kappa}^{\alpha}\right]_{c_{i}^{\kappa}}$
$\displaystyle=\frac{1}{N_{\rm m}\tilde{\tau}}(\hat{\rho}_{\rm
on}^{\alpha}-\hat{\rho}_{\rm
b}^{\alpha})^{2}\left(\frac{\beta_{\alpha}}{\beta_{\rm
m}}\right)^{2}\left[\left(\frac{\partial f_{\alpha}}{\partial
c^{\kappa}}\right)^{2}2c^{\kappa}c_{0}\Phi_{\kappa}^{\alpha}\right]_{c_{i}^{\kappa}}$
(61)

where $c_{0}\equiv N_{\rm m}/(D_{\text{int}}l\tilde{\tau})$ now is defined
based on the longest mRNA lifetime $\tilde{\tau}$. This recovers almost
exactly the same formula as before, Eq. (LABEL:eqInputNoiseMaternal), apart
from the prefactor $(\beta_{\alpha}/\beta_{\rm m})^{2}$, which vanishes when
all burst sizes $\beta_{\alpha}$ are equal.

A completely analogous calculation yields the corresponding term for the other
(zygotic) input noise contributions (Eq. (LABEL:eqInputNoiseZygotic) in the
model without mRNA):

$\displaystyle\mathcal{J}_{g}^{\zeta}$ $\displaystyle=\frac{1}{N_{\rm
m}\tilde{\tau}}\frac{1}{\gamma_{\rm m}}(\hat{\rho}_{\rm
on}^{\alpha}-\hat{\rho}_{\rm
b}^{\alpha})^{2}\left(\frac{\beta_{\alpha}}{\beta_{\rm
m}}\right)^{2}\left[\left(\frac{\partial f_{\alpha}}{\partial
g^{\zeta}}\right)^{2}2g^{\zeta}c_{0}\Phi_{\zeta}\right]_{\bar{g}_{i}^{\zeta}}$
(62)

The output noise follows the same ansatz as in Sec. 1.4 and is given by the
sum of the average mRNA production and degradation rate again multiplied with
prefactor $(\beta_{\alpha}/N_{\rm m})^{2}$, but care has to put into its
correct normalization. Repeating the original calculation with the required
modifications yields

$\displaystyle\left\langle\left[\mathcal{O}^{\alpha}_{i}\right]^{2}\right\rangle$
$\displaystyle=\frac{\beta_{\alpha}^{2}}{N_{\rm m}^{2}}\left[\left(r_{\rm
b}^{\alpha}+(r_{\rm on}^{\alpha}-r_{\rm
b}^{\alpha})f_{\alpha}(\mathbf{\bar{g}}_{i},\mathbf{c}_{i})\right)+\frac{1}{\tilde{\tau}_{\alpha}}\left\langle
M_{i}^{\alpha}\right\rangle\right]$
$\displaystyle=\frac{\beta_{\alpha}}{N_{\rm
m}}\left[\frac{\beta_{\alpha}}{\beta_{\rm m}r_{\rm
on}\tilde{\tau}}\left(r_{\rm b}^{\alpha}+(r_{\rm on}^{\alpha}-r_{\rm
b}^{\alpha})f_{\alpha}(\mathbf{\bar{g}}_{i},\mathbf{c}_{i})\right)+\frac{1}{\tilde{\tau}_{\alpha}}\frac{\beta_{\alpha}\left\langle
M_{i}^{\alpha}\right\rangle}{N_{\rm m}}\right]$
$\displaystyle=\beta_{\alpha}\times\frac{1}{N_{\rm
m}}\left[\frac{1}{\tilde{\tau}}\left(\frac{\beta_{\alpha}}{\beta_{\rm
m}}\right)\left(\hat{r}_{\rm b}^{\alpha}+(\hat{\rho}_{\rm
on}^{\alpha}-\hat{r}_{\rm
b}^{\alpha})f_{\alpha}(\mathbf{\bar{g}}_{i},\mathbf{c}_{i})\right)+\frac{1}{\tilde{\tau}_{\alpha}}\bar{g}_{i}^{\alpha}\right]$
(63)

where in the last line we make use of $\beta_{\alpha}\left\langle
M_{i}^{\alpha}\right\rangle/N_{\rm m}=\left\langle
G_{i}^{\alpha}\right\rangle/G_{\rm max}=\bar{g}_{i}^{\alpha}$. While for equal
burst sizes $\beta_{\alpha}$ the prefactor $\beta_{\alpha}/\beta_{\rm m}$
preceding the production term again vanishes, the overall prefactor
$\beta_{\alpha}$ that multiplies the whole output noise always remains. As
known from simpler expression models, this precisely reflects the
superpoissonian increase of the output noise due to bursty protein production.

Taken together the normalized noise covariances in the model with mRNA read:

$\displaystyle\left\langle\left[\mathcal{N}^{\alpha}_{ii}\right]^{2}\right\rangle=\frac{1}{N_{\rm
m}}\Bigg{\\{}$
$\displaystyle\quad\quad\frac{\beta_{\alpha}}{\tilde{\tau}}\left(\frac{\beta_{\alpha}}{\beta_{\rm
m}}\right)\left[\hat{r}_{\rm b}^{\alpha}+(\hat{\rho}_{\rm
on}^{\alpha}-\hat{r}_{\rm
b}^{\alpha})f_{\alpha}(\mathbf{\bar{g}}_{i},\mathbf{c}_{i})\right]+\frac{1}{\tilde{\tau}_{\alpha}}\bar{g}_{i}^{\alpha}$
“output noise” $\displaystyle\quad+\frac{(\hat{\rho}_{\rm
on}^{\alpha}-\hat{\rho}_{\rm
b}^{\alpha})^{2}}{\tilde{\tau}}\left(\frac{\beta_{\alpha}}{\beta_{\rm
m}}\right)^{2}\sum_{\kappa}\left[\left(\frac{\partial f_{\alpha}}{\partial
c^{\kappa}}\right)^{2}2c^{\kappa}c_{0}\Phi_{\kappa}^{\alpha}\right]_{c_{i}^{\kappa}}$
“maternal input noise” $\displaystyle\quad+\frac{1}{\gamma_{\rm
m}}\frac{(\hat{\rho}_{\rm on}^{\alpha}-\hat{\rho}_{\rm
b}^{\alpha})^{2}}{\tilde{\tau}}\left(\frac{\beta_{\alpha}}{\beta_{\rm
m}}\right)^{2}\sum_{\zeta}\left[\left(\frac{\partial f_{\alpha}}{\partial
g^{\zeta}}\right)^{2}2g^{\zeta}c_{0}\Phi_{\zeta}^{\alpha}\right]_{\bar{g}_{i}^{\zeta}}$
“zygotic input noise” $\displaystyle\quad+2d\cdot
h_{\alpha}\bar{g}_{i}^{\alpha}+\sum_{\nu}h_{\alpha}\bar{g}_{\nu}^{\alpha}~{}\Bigg{\\}}$
“diffusion noise” (64)
$\displaystyle\left\langle\left[\mathcal{N}^{\alpha}_{i\nu}\right]^{2}\right\rangle=-\frac{h_{\alpha}}{N_{\rm
m}}\left(\bar{g}_{i}^{\alpha}+\bar{g}_{\nu}^{\alpha}\right)$ “diffusion noise”
(65)

In our optimizations we assume all burst sizes $\beta_{\alpha}$ to be equal,
such that the increase of the output noise is the only (but significant)
additional noise contribution arising from incorporation of mRNA and bursty
protein production. Note that since all noise covariances derived in this
section and Eq. (60) recover the structure of the corresponding formulae
derived in the model without mRNA, all derivations presented in Sec. 1.3 also
apply to the extended model, such that the computation of the covariances can
be carried out in exactly the same way, with the slight change that protein
lifetimes have to be replaced by mRNA lifetimes.

##  2 Numerical computation of means and covariance matrix of gap gene
expression

The spatial-stochastic model of gap gene expression derived in the previous
sections establishes two sets of equations for the means and covariances of
the gap gene expression levels, respectively, which we solve numerically. Note
that we need to solve for the means first because they appear in the
definition of the coupled linear system relating the entries of the covariance
matrix.

###  2.1 Forward integration of the means

As already outlined in Sec. 1.2, Eqs. (46) and (47) for the mean expression
levels must be integrated numerically. We start with zero expression levels at
$t=0$ and use the built-in MATLAB ODE solver ode23t. We account for nuclear
divisions that result in a doubling of production capacities upon every
division by ramping up the production rate at the start of every new nuclear
cycle $c$ via a factor $\phi(c)=2^{c}/2^{14}$, meaning that the full
production rate is only reached at the beginning of nuclear cycle 14
($\phi(14)=1$). The division times (end times of the respective nuclear cycle)
used in our model are listed in Table 1 and follow the timeline described in
[5].

Nuclear cycle | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13
---|---|---|---|---|---|---|---|---|---|---|---|---|---
Division time [min] | 12 | 20 | 28 | 36 | 44 | 52 | 60 | 68 | 76 | 84 | 94 | 110 | 126
Table 1: Nuclear division times / end times of nuclear cycles employed in our
model, following the timeline described in [5].

###  2.2 Numerical computation of the covariance matrix

Equation 17 defines a coupled linear system for computing the entries of the
generalized stationary covariance matrix that comprises all covariances
$C_{ij}^{\alpha\beta}$ across both all the gene species and all positions.
This linear system can itself be represented as a matrix equation

$\displaystyle\mathcal{M}\mathbf{c}=\mathbf{G}$ (66)

where $\mathbf{c}$ is a vector containing the entries $C_{ij}^{\alpha\beta}$
in a defined order, $\mathcal{M}$ a matrix collecting the coupling terms
relating different entries $C_{ij}^{\alpha\beta}$ and $C_{kl}^{\gamma\delta}$
with each other, and $\mathbf{G}$ a vector containing the respective right-
hand side values that do not depend on $C_{ij}^{\alpha\beta}$, in our case the
noise covariance terms. Thus, by constructing matrix $\mathcal{M}$ and vector
$\mathbf{c}$ with a defined bookkeeping (order of $C_{ij}$ entries in
$\mathbf{c}$) and imposing the corresponding right-hand side noise covariances
vector $\mathbf{G}$ we can obtain $\mathbf{c}$ by numerial matrix inversion,

$\displaystyle\mathbf{c}=\mathcal{M}^{-1}\mathbf{G}$ (67)

and use (map back) its entries for constructing the covariance matrix
$\mathcal{C}=(C_{ij}^{\alpha\beta})$.

Note that the matrix $\mathcal{M}$ typically is large and grows quadratically
in the system size variables. Already for a one-dimensional system of typical
size, i.e. for $N_{\rm x}=70$ and $N_{\rm y}=1$ nuclei in the respective
dimensions, and 4 gap genes modeled, $\mathcal{M}$ is a $78400\times 78400$
quadratic matrix with more than 6 billion entries ($70\times 70\times 4\times
4=78400$). However, when the short-correlations assumption (Sec. 1.3.5) is
imposed, the vast majority of entries of $\mathcal{M}$ becomes zero. This
sparse matrix can be inverted considerably faster and for the system size
parameters above can be carried out in less than one second of computation
time.

##  3 Calculation of positional information

Our aim is to compute the positional information $I\equiv I(\mathbf{g};x)$,
which quantifies how much information the local noisy expression levels of the
gap genes ($\mathbf{g}$) carry about the position at some fixed developmental
time when their readout occurs. Based on and validated by our earlier work [6,
7], we assume that the local distributions of expression levels,
$P(\mathbf{g}|x_{i})$, can be written as multivariate Gaussian distributions
and therefore are completely characterized by the set of mean expression
levels $\bar{g}(x_{i})$ and the local covariance matrices
$C^{\alpha\beta}(x_{i})=C^{\alpha\beta}_{ii}$ at position $x_{i}$ predicted by
our spatial-stochastic model. However, even with this approximation it remains
difficult to compute the positional information by direct integration, because
this technically requires integration over the four-dimensional distribution
of all possible combinations of local expression levels. This is numerically
unfeasible even with Monte-Carlo integration if no further approximations are
made [6]. Therefore, here we opted for an approach which, instead of directly
computing the positional information $I(\mathbf{g};x)$, computes a lower bound
on this quantity from an averaged decoding map (confusion matrix) at the
relevant time point ($T_{\text{meas}})$; the corresponding decoding maps, when
averaged, then can be integrated, bounding $I(\mathbf{g};x)$ from below [8].
The tightest lower bound results from employing Bayes’ optimal decoder [9],
which we also opt for here. Additionally, for the case of Drosophila gap
genes, we have previously shown that the decoding-based bound that we employ
here is tight to within experimental and numerical errors for the wild-type
[10], and moreover predictive about the expression profiles in mutants with
knocked-out maternal inputs or gap genes [9].

In order to construct the optimal decoder, we first apply Bayes’ rule as
follows:

$\displaystyle P(x^{*}_{j}|\mathbf{g})P(\mathbf{g})$
$\displaystyle=P(\mathbf{g}|x^{*}_{j})P(x^{*}_{j})$
$\displaystyle\Leftrightarrow\quad\quad\quad P(x^{*}_{j}|\mathbf{g})$
$\displaystyle=\frac{P(\mathbf{g}|x^{*}_{j})P(x^{*}_{j})}{P(\mathbf{g})}=\frac{P(\mathbf{g}|x^{*}_{j})}{\sum_{j}P(\mathbf{g}|x^{*}_{j})}$
(68)

Here the distribution on the left-hand side of the second line is the
posterior distribution, i.e. the distribution over decoded positions
$x^{*}_{j}$ implied by reading the set of expression levels $\mathbf{g}$.
$P(x^{*}_{j})=\frac{1}{N_{\rm x}}$ is the prior distribution over positions
which in our case is uniform (reflecting initial complete ignorance of the
nuclei about their position and the nuclei being uniformly spaced along the AP
axis) and therefore cancels in the last formula. Equation (68) tells us that,
for any fixed set of expression levels $\mathbf{g}$, we can obtain the
posterior distribution by tabulating the values $P(\mathbf{g}|x^{*}_{j})$ for
all $x^{*}_{j}$ and normalizing them by their sum over all positions.

For any specific set of local expression levels $\mathbf{g}_{n}(x_{i})$, e.g.
reflecting the local expression levels at the readout time $T_{\text{meas}}$
from one individual embryo $n$, we can now construct a decoding map
$P_{n}(x^{*}_{j}|x_{i})$ as follows:

$\displaystyle
P_{n}(x^{*}_{j}|x_{i})=P(x^{*}_{j}|\mathbf{g}_{n}(x_{i}))=\frac{P(\mathbf{g}_{n}(x_{i})|x^{*}_{j})}{\sum_{j}P(\mathbf{g}_{n}(x_{i})|x^{*}_{j})}.$
(69)

$P_{n}(x^{*}_{j}|x_{i})$ establishes the mapping between the actual position
$x_{i}$ and the locally implied positions $x^{*}_{j}$ for a single decoding
process at all positions along the embryo axis.

Asking for the general performance of the decoding process, i.e. the amount of
positional information that the gap genes can encode on average, we can repeat
all the steps leading up to the construction of an individual decoding map
$P_{n}(x^{*}_{j}|x_{i})$ for a set of $n=1,..,N_{\rm e}$ embryo samples with
corresponding expression levels $\mathbf{g}_{n}(x_{i})$ sampled from the
Gaussian output distributions $P(\mathbf{g}|x_{i})$, and then obtain the
average decoding map:

$\displaystyle\langle P\rangle(x^{*}_{j}|x_{i})\equiv\langle
P_{n}(x^{*}_{j}|x_{i})\rangle_{N_{\rm e}}=\frac{1}{N_{\rm
e}}\sum_{n}P_{n}(x^{*}_{j}|x_{i})$ (70)

The decoding capacity can be computed from any (normalized) decoding map
$P(x^{*}|x)$ via the known formula for mutual information:

$\displaystyle I(x^{*};x)$
$\displaystyle=\int_{0}^{L}\int_{0}^{L}P(x^{*}|x)\log_{2}\frac{P(x^{*}|x)}{P_{x}(x^{*})P_{x^{*}}(x)}dx^{*}dx$
(71)

where by convention the integrand equals zero when $P(x,x^{*})=0$ and $P_{x}$
and $P_{x^{*}}$ denote the marginal distributions:

$\displaystyle P_{x}(x^{*})$ $\displaystyle=\int_{0}^{L}P(x^{*}|x)dx$ (72)
$\displaystyle P_{x^{*}}(x)$ $\displaystyle=\int_{0}^{L}P(x^{*}|x)dx^{*}$ (73)

These integrals run over the whole domain of support of $x$ and $x^{*}$, i.e.
the whole embryo length. In our discrete system the integrals are sums over
all nuclear positions, while the summation is carried out over the average
decoding map $\langle P\rangle(x^{*}_{j}|x_{i})$:

$\displaystyle I(x^{*};x)$ $\displaystyle=\sum_{i}\sum_{j}\langle
P\rangle(x^{*}_{j}|x_{i})\log_{2}\frac{\langle
P\rangle(x^{*}_{j}|x_{i})}{P_{x}(x^{*}_{j})P_{x^{*}}(x_{i})}$ $\displaystyle
P_{x}(x^{*}_{j})$ $\displaystyle=\sum_{i}\langle P\rangle(x^{*}_{j}|x_{i})$
$\displaystyle P_{x^{*}}(x_{i})$ $\displaystyle=\sum_{j}\langle
P\rangle(x^{*}_{j}|x_{i})$ (74)

We found empirically that averaging the individual decoding maps over $N_{\rm
e}=50$ Monte Carlo embryo “samples” (i.e., using $N_{\rm x}\times N_{\rm e}$
sampled expression vectors $\mathbf{g}_{n}(x_{i}))$) is sufficient for
reproducible computation of $I(x^{*};x)$, which serves as an input to the
stochastic optimization algorithm described below.

Note that the portrayed way of estimating the positional information closely
follows the biological decoding process, where the developing nuclei have to
interpret the experienced gap expression levels in a short time frame, thus
converting a _single_ “snapshot” of the noisy pattern into a positional
estimate. Thus, while as described above, the computation of $I(x^{*};x)$ (as
a tractable lower bound on the positional information $I(\mathbf{g};x)$) is
technically motivated, it is intriguing to think that $I(x^{*};x)$ itself is
the quantity of primary biological relevance.

##  4 Optimization procedure

We use an individually customized simulated annealing algorithm for random
optimization of the positional information $I$ over the set of regulatory
parameters and (in some cases) the gap protein diffusion constant, followed by
a final gradient descent optimization run after random optimization has
settled into a minimum. The optimization procedure starts with a totally
random set of parameters uniformely sampled on the prescribed parameter bounds
(see Table 3 and Sec. 5.2). The parameters then are iteratively changed in a
random fashion (described in detail below). Upon each iteration, we solve our
spatial-stochastic embryo model for the altered set of parameters, and
(re)compute the positional information $I$ of the pattern as described in Sec.
3. The parameter change is always accepted if it leads to an increase of $I$;
if the parameter change decreases $I$, the change is still accepted with a
finite probability according to the Metropolis-Hastings algorithm, with an
acceptance probability (effective temperature) that is lowered with every
subsequent optimization step. The detailed cooling protocol is specified
further below in Sec. 4.2.

###  4.1 Parameter changes and acceptance

For each parameter $n$ we initially define a minimal and maximal bound,
$P^{\text{min}}_{n}$ and $P^{\text{max}}_{n}$, which in most optimizations are
very generous and cover the full range of biophysically relevant values for
the given parameter type (see Sec. 5.2 for details). The implied parameter
interval $\Delta P_{n}=P^{\text{max}}_{n}-P^{\text{min}}_{n}$ then is
subdivided into $N_{\delta}$ equidistant subintervals $\delta P_{n}$. While in
classical approaches the subinterval boundaries together with
$P^{\text{min}}_{n}$ and $P^{\text{max}}_{n}$ could be used as a lattice of
($N_{\delta}+1$) discrete parameter values on which the optimization can
explore the range $\Delta P_{n}$, we found that such discretization strongly
restricts the set of possible expression patterns, and therefore opted for the
following approach in which all parameter combinations in the bounded
hyperparameter space are accessible: At the start of the optimization the
parameter value $P_{n}^{0}$ for each optimization parameter $n$ is randomly
chosen on the interval $\Delta P_{n}$ with uniform probability
($p(P_{n}^{0})=1/\Delta P_{n}$); this defines the initial parameter vector
$\mathbf{\Theta}^{0}=(P^{0}_{1},...,P^{0}_{N_{p}})$, where $N_{p}$ is the
total number of optimized parameters. At each iteration $i\geq 1$ of the
random optimization, first we randomly choose which parameter to change by
uniformly sampling a random integer number $n_{i}\in[1,N_{p}]$. Second, a new
proposed value for parameter $n_{i}$ is chosen by uniformly sampling it within
the $2\delta P_{n_{i}}$ interval around the current value $P_{n_{i}}^{i-1}$,
i.e.

$\displaystyle r_{i}$ $\displaystyle=\mathcal{R}_{i}([-1,1])$ $\displaystyle
P_{n_{i}}^{i}$ $\displaystyle=P_{n_{i}}^{i-1}+r_{i}\delta P_{n_{i}}$ (75)

where $\mathcal{R}_{i}([-1,1])$ denotes a continuous random number uniformly
sampled on the interval $[-1,1]$. If the new value $P_{n_{i}}^{i}$ reaches
beyond the bounds $P^{\text{min}}_{n}$ or $P^{\text{max}}_{n}$, it is set
(truncated) to the respective boundary value. The proposed parameter vector
$\mathbf{\Theta}_{i}^{\text{prop}}$ for the current optimization iteration
then is constructed by replacing the $n_{i}$-th entry of the current vector
$\mathbf{\Theta}_{i-1}$ by $P_{n_{i}}^{i}$.

We then compute the mean gap gene expression levels and their full covariance
matrix for the proposed parameter vector $\mathbf{\Theta}_{i}^{\text{prop}}$
and the corresponding positional information $I_{i}^{\text{prop}}$ (see Sec.
3). The new value of the objective function $I_{i}^{\text{prop}}$ then is used
to carry out the “accept or reject” test according to the Metropolis-Hastings
algorithm as follows:

  1. 1.

Compute objective function difference:

$\displaystyle\Delta I=I_{i}^{\text{prop}}-I_{i-1}$

  2. 2.

Compute the acceptance probability, given current temperature $T=T_{i}$:

$\displaystyle p_{\text{acc}}=\min(1,\exp(\Delta I/T))$

  3. 3.

Metropolis-Hastings test:

     * •

Sample a uniformly distributed random number $\mathcal{R}\in[0,1]$.

     * •

If $\mathcal{R}\leq p_{\text{acc}}$, accept the change:

       * –

Update parameter vector,
$\mathbf{\Theta_{i}}\leftarrow\mathbf{\Theta}_{i}^{\text{prop}}$.

       * –

Update objective function, $I_{i}\leftarrow I_{i}^{\text{prop}}$.

     * •

Else reject:

       * –

Parameter vector remains unchanged,
$\mathbf{\Theta_{i}}\leftarrow\mathbf{\Theta}_{i-1}$.

       * –

Objective function value remains unchanged, $I_{i}\leftarrow I_{i-1}$.

  4. 4.

Decrease temperature according to cooling protocol (see below), $T\leftarrow
T_{i+1}$

  5. 5.

Iterate, $i\leftarrow i+1$

The random optimization procedure ends when $N_{\rm opt}$ optimization steps
have been attempted. By default we set $N_{\rm opt}=250*N_{\rm x}*N_{\rm g}$,
thus scaling the number of optimization steps with the system size (number of
positions $N_{\rm x}$ and number of simulated gap genes $N_{\rm g}$). We found
that with our model and optimization procedure, this number of optimization
steps is sufficient for finding optimized patterns whose positional
information is comparable to or even slightly larger than the $\sim
4.2~{}\text{bits}$ of the measured gap gene profiles [7]. Finally, the local
maximum found by random optimization is taken as the starting point for a
final gradient descent optimization run as described in Sec. 4.3.

###  4.2 Cooling protocol

The temperature is decreased according to an initially defined exponentially
decaying temperature ramp as follows:

$\displaystyle T_{i}=T_{\text{max}}\exp(-t_{i})$ (76)

where the $t_{i}$ are $N_{\rm opt}$ numbers uniformly spaced on the interval
$[0,\log(T_{\text{max}}/T_{\text{min}})]$, such that $T_{0}=T_{\text{max}}$
and $T_{N_{\rm opt}}=T_{\text{min}}$.

After analysing the typical changes in positional information $\Delta I$ upon
changing the optimization parameters we empirically explored different maximal
and minimal temperature bounds $T_{\text{max}}$ and $T_{\text{min}}$, and
found that the choice of $T_{\text{max}}=0.2$ and $T_{\text{min}}=0.01$ leads
to robust optimization behavior, i.e. access of local information maxima in
parameter space.

###  4.3 Final gradient descent

To ensure that the optimization fully settles into the local maximum
$I^{*}_{\text{rand}}$ of the objective function (positional information) found
by random optimization, we start an additional gradient-descent optimization
run from $I^{*}_{\text{rand}}$. For this we used the built-in MATLAB fmincon
optimizer with negative positional information as an objective function and
$N_{\rm opt}/10$ optimization steps, and the same parameter bounds as in
random optimization. In most cases, the final gradient descent leads to only
small changes in the parameter values, reflecting that random optimization
already settled to a point close to the actual local maximum.

###  4.4 Exploration of Drosophila-like patterns via constrained optimization

We observed that optimized solutions obtained by starting from random
parameter values can be degenerate: we found different gap gene expression
patterns that encoded similarly high amounts of positional information. To
systematically test whether any of these solutions lies in the vicinity of the
observed fruit fly gap gene expression pattern, we employed a recently
developed statistical framework which we reported in a separate publication
[11]; we summarize the workflow below.

We first tested whether our model is expressive enough to be able to generate
expression patterns akin to the measured Drosophila gap gene expression
pattern. This was done by minimizing the pattern distance (defined further
below) between the measured mean expression profiles and the mean expression
profiles generated by our model, without maximizing the encoded positional
information (i.e., by traditional “fitting”). These tests confirmed that such
Drosophila-like patterns can be generated by our model (with the exception of
the posterior Hunchback bump, which we attribute to the likely presence of
additional Hb enhancers in the WT system, not included in our model). We then
carried out constrained optimization runs, in which we maximized positional
information but supplemented by the pattern distance as a Lagrange multiplier
term, which generates a small “pulling force” that biases the stochastic
search in parameter space towards expression patterns that are similar to the
observed expression pattern. Note that there is no a priori guarantee that an
information-optimizing solution close to the Drosophila WT exists. As a last
step, we chose the most similar of these optimal, Drosophila-like patterns and
confirmed that they remain a maximum of positional information when the
pulling force is taken away again. The key idea behind this workflow is that
while the small “pulling force” helps break the degeneracy among information-
maximizing solutions to guide the stochastic search towards patterns similar
to the Drosophila wild-type, the force needs to be small compared to the
information term, so that optimization remains dominated by parameter changes
that increase the information. This “guided optimization” is very different in
the landscape and outcomes compared to pure fitting, which can be recovered in
the limit $\lambda\rightarrow\infty$ in Eq.(77) in the function defined below;
see also main paper Fig. 2F.

The extended objective function implementing optimization of positional
information guided by similarity with the measured expression pattern is
defined as

$\displaystyle\mathcal{L}\equiv I(\mathbf{g};x)-\lambda
D(\mathbf{g},\hat{\mathbf{g}})$ (77)

where $\mathbf{g}=(g_{i}^{\alpha})$ represents represents the pattern being
optimized and $\hat{\mathbf{g}}=(\hat{g}_{i}^{\alpha})$ the measured pattern,
with expression levels normalized to the range $[0,1]$ in both patterns. The
pattern distance function $D(\mathbf{g},\hat{\mathbf{g}})$ is defined as a
finite p-norm,

$\displaystyle
D(\mathbf{g},\hat{\mathbf{g}})\equiv\underset{\pi}{\min}\left(\frac{1}{N}\sum_{\alpha}\sum_{i}|g_{i}^{\alpha}-\hat{g}_{i}^{\mathcal{P}_{\pi}(\alpha)}|^{p}\right)^{1/p}$
(78)

where $\alpha=1..N_{\rm g}$ runs over all gap gene species, $i=1..N$ over all
positions/volumes in the model, and index $\pi$ runs over all permutation
functions $\mathcal{P}_{\pi}(\alpha)$ that map a gene species in $\mathbf{g}$
to a gene species in $\hat{\mathbf{g}}$. This accounts for the fact that in
the optimized pattern the gene identities are not a priori specified, such
that even for identical patterns the distance measure would be nonzero if the
distances are not computed on the “matching” genes that minimize the measure.
Interestingly, we found that the pattern distance measure performs better for
$p\geq 3$ than for $p=2$, such that we set $p=3$ in all our distance
computations. Moreover, we found that the choice of $\lambda=6$ for the
Lagrange mulitplier can successfully filter out high-information solutions
that are indeed similar to the measured Drosophila pattern.

We verified that the extended objective function successfully maximized
$I(\mathbf{g};x)$ in two ways. Firstly, we carried out parameter perturbations
away from the optimal parameters which lead to reduction of $I(\mathbf{g};x)$
or negligible changes, as described in SI Fig. 3. Secondly, we computed
$I(\mathbf{g};x)$ on “purely fitted” solutions in which only $\lambda
D(\mathbf{g},\hat{\mathbf{g}})$ was minimized (see Fig. 2E of the main text).
While the resulting mean expression patterns are similar to the observed
pattern, the positional information for these solutions is considerably lower,
highlighting that optimization of $I(\mathbf{g};x)$ is crucial for obtaining
the parameter set that reproduces the high information content of the measured
gap gene pattern.

##  5 Parameters, simulation geometry and optimization constraints

We discriminate between parameters that are optimized in our optimization runs
as described in Sec. 4 and those that are set to fixed values, which include
the parameters defining the simulated embryo geometry. The division of
parameters into fixed and optimized is not arbitrary: fixed parameters
represent biophysical constraints on the system that set the intrinsic levels
of noise and thus limit positional information. We will describe the fixed
parameters first and then proceed to the optimized parameters.

###  5.1 Fixed parameters and geometric constraints

Table 2 gives an overview of the model parameters that we keep fixed in our
optimizations, although some of them are varied between different ensembles of
optimization runs (such as the number of output / gap genes $N_{\rm g}$ or the
number of input gradients $N_{\rm i}$, see Sec. 8 for details). The geometric
parameters are kept fixed in all optimization runs and are chosen to reproduce
the geometry of the Drosophila embryo in nuclear cycle 14. While our model in
principle allows to study arbitrary spatial arrangements of the embryo nuclei,
we restricted our simulatons to two cases: (1.) A 2d model in which the embryo
is abstracted as a cylinder with $N_{\rm x}$ equally spaced nuclei along the
axial direction and $N_{\rm y}$ equally spaced nuclei along the circumference
with periodic boundary conditions along the $y$-dimension, i.e. the nuclei
with indices $y=1$ are coupled to the nuclei at $y=N_{\rm y}$. (2.) A reduced
1d model in which the coupling along the $y$-dimension is skipped; this
arrangement is representative of a single “string” of nuclei along the half-
perimeter of the ellipsoidal section of the embryo. To our surprise we found
in early simulations that there is little difference in optimal values of
positional information between the 1d and the 2d model, whereas the
computational effort for solving the model is considerably higher in the 2d
setting. Moreover we found that the 1d model already is capable of reproducing
the measured features of the gap gene system very well. For the optimization
runs presented in this work, we therefore opted for the 1d model in order to
boost computational efficiency. Note that a reduction of the simulated lattice
of positions has a profound effect on the computational cost, because the
numerical computation of the cross-covariances in the system requires the
inversion of a matrix that grows quadratically with the number of considered
positions (see Sec. 2.2). For the same reason, we reduced the spatial size of
the system further by excluding the very posterior positions of the embryo
where the wild-type embryo does not show any expression of gap genes. This
“posterior cut-off” amounted to the $N_{\rm c}=5$ posteriormost nuclei.

Our model includes three morphogen input gradients: an anterior gradient (A),
representative of Bicoid, a posterior gradient (P), representative of Nanos,
and a terminal system (T), representative of Torso-like, that emerges
symetrically from the opposite poles of the embryo. We assume an exponential
shape for the gradients, and that the posterior gradient is a perfect mirror
image of the anterior gradient with respect to mid-embryo. The length scale
$\lambda_{\rm AP}$ is chosen based on the observed properties of the Bicoid
gradient, $\lambda_{\rm AP}=20\%~{}{\rm EL}$; the length scale of the terminal
system gradients is chosen four times shorter, $\lambda_{\rm T}=5\%~{}{\rm
EL}$. Restrictions imposed on the possible regulatory action of the gradients
are discussed further below.

While our model is derived such that it allows for different protein and mRNA
production rates and lifetimes, we make the simplifying assumption that these
parameters are equal for all gap gene species, as indicated in Table 2. The
same holds for the gap protein burst size and the internal diffusion constant
of the transcription factor proteins (both zygotic and maternal). We further
assume that the binding sites of the regulated genes all have the same target
size (binding region length) $l$.

We set the value of the inter-nuclear diffusion constant $D$, assumed to be
equal for all gap genes, to $D=0.5~{}\mu m^{2}/s$ by default, and later varied
it accross several optimization ensembles (cf. Fig.3E in the main text).

Name | Symbol | Value
---|---|---
No. of output genes | $N_{\rm g}$ | 1–5, default = 4
No. of maternal input morphogens | $N_{\rm i}$ | 1–3, default = 3
No. of nuclei along embryo axis | $N_{\rm x}$ | 70
Posterior cut-off | $N_{\rm c}$ | 5
Corr. length cut-off (SCA) | $N_{\rm SCA}$ | variable, default = 1
Nuclear radius | $R$ | 3.25 $\mu m$
Nuclear volume | $\Omega$ | 144 $\mu m^{3}$
Inter-nuclear distance | $\Delta$ | 8.5 $\mu m$
Max. gap gene mRNA copy no. (single species) | $M_{\rm m}$ | 540
Max. protein burst size | $\beta_{\rm m}$ | 12
Max. no. of gap proteins (single species) | $N_{\rm m}$ | 6480
Max. input concentration | $c_{\rm m}$ | 150 nM = 90/$\mu m^{3}$
Input gradient length (anterior and posterior) | $\lambda_{\rm AP}$ | 0.2 EL
Input gradient length (terminal system) | $\lambda_{\rm T}$ | 0.05 EL
mRNA life times | $\tilde{\tau}_{\alpha}\equiv\tilde{\tau}$ | 20 min
Protein life times | $\tau_{\alpha}\equiv\tau$ | 10 min
mRNA prod. rates at full induction | $r_{\rm on}^{\alpha}\equiv r_{\rm on}$ | 27/min
Basal mRNA prod. rates | $r_{\rm b}^{\alpha}\equiv r_{\rm b}$ | 0.027/min
Protein burst size | $\beta_{\alpha}\equiv\beta_{\rm m}$ | 12
Internal TF diff. const. | $D_{\text{int}}$ | 10 $\mu m^{2}/s$
TF binding site length | $l$ | 10 nm

Table 2: Overview of fixed biophysical and geometric parameters.

###  5.2 Optimization parameters

An overview of the model parameters that are optimized for maximizing
positional information is shown in Table 3. The core set consists of the
regulatory parameters that determine the regulation of the gap genes by the
maternal inputs, $H_{M}^{\alpha\kappa}$ and $K_{M}^{\alpha\kappa}$, and their
mutual and self-regulation, $H_{G}^{\alpha\zeta}$ and $K_{G}^{\alpha\zeta}$;
note that while we separate maternal from zygotic regulation notationally, $K$
and $H$ describe the same quantities (regulation threshold and regulation
strength, respectively) in both cases (see Sec. 1.5 for details). These
parameters were optimized on the bounds indicated in the respective column of
Table 3, which were chosen such that all genes have the possibility to be
fully expressed or completely inactive throughout the whole system, i.e.
generous enough to guarantee that the expression patterns emerging are not
constrained by the parameter bounds but by the biophysical laws imposed by our
model (specifically, by the structure of the regulatory function and the
diffusive coupling).

While keeping the regulatory interactions as generic as possible, we imposed
the following additional constraints aiding the optimization performance.
Firstly, optimization looks for the best (if any) activating regulatory
parameters for the anterior and posterior morphogen gradients, and for the
best (if any) repressive regulatory parameters for the terminal morphogen
gradient. Secondly, optimization looks for the best (if any) repressive
interactions between gap genes, and the best (if any) self-activating auto-
regulation. These two constraints are informed by our prior work [4, 2], which
derived optimal topologies for feed-forward and interacting networks, as well
as for self-interaction. In this prior work we had analytic control over the
optimization scenarios, providing a solid theoretical foundation for the
nature of optimal solutions in the more complicated case we analyze
numerically here. For performance reasons, we allow only two of the gap genes
to couple to the posterior gradient in the production runs, since that
decreases the number of optimization parameters; we verified that this
constraint has little influence on the optimal positional information values
reached in our standard ensemble described here. Taken together, these
constraints do not change the nature of the optimal solutions that we find,
but speed up the search time such that multiple optimization scenarios can be
tractably explored. Relaxing these constraints could produce additional
solutions that locally maximize positional information, but it will not remove
the solutions we do find.

For two optimal ensembles (solid circles Fig.3E in the main text), we also
optimized the inter-nuclear diffusion constant $D$ within the bounds indicated
in Table 3. Interestingly, the average optimal value found, $D_{\rm
opt}=0.4~{}\mu m^{2}/s$, is very close to our initially chosen default value
$D=0.5~{}\mu m^{2}/s$ (see also Sec. 8.1).

Name | Symbol | No. contained | Unit | Default | Default
---|---|---|---|---|---
|  | in model |  | bounds | value
Maternal regulation thresholds | $K_{M}^{\alpha\kappa}$ | $N_{\rm g}\times N_{\rm i}$ | $1/\mu m^{3}$ | $[0.4,2000]$ |
Maternal regulation strength | $H_{M}^{\alpha\kappa}$ | $N_{\rm g}\times N_{\rm i}$ |  | $[0.01,50]$ |
Zygotic regulation thresholds | $K_{G}^{\alpha\zeta}$ | $N_{\rm g}\times N_{\rm g}$ | $N_{\rm m}$ | $[5\cdot 10^{-7},2000]$ |
Zygotic regulation strength | $H_{G}^{\alpha\zeta}$ | $N_{\rm g}\times N_{\rm g}$ |  | $[0.01,50]$ |
Protein inter-nuclear diffusion const. | $D$ | 1 | $\mu m^{2}/s$ | $[0.01,10]$ | $0.5$
Total no. |  | $2N_{\rm g}^{2}+2N_{\rm g}N_{\rm i}+1$ |  |  |
\- for $N_{\rm g}=2$ |  | 21 |  |  |
\- for $N_{\rm g}=3$ |  | 37 |  |  |
\- for $N_{\rm g}=4$ |  | 57 |  |  |
\- for $N_{\rm g}=5$ |  | 81 |  |  |

Table 3:  Overview of optimized biophysical parameters. The regulatory
parameters are always optimized starting from random values uniformly sampled
within the indicated parameter bounds. For the analysis summarized in Fig. 3E
and F of the main article, we also optimized the inter-nuclear diffusion
constant $D$; otherwise the default value was used. The bounds indicated for
the regulatory strength terms $H_{M}^{\alpha\kappa}$ and $H_{G}^{\alpha\zeta}$
refer to the absolute value of these quantities, i.e. whenever
$H_{M}^{\alpha\kappa}$ or $H_{G}^{\alpha\zeta}$ were restricted to be negative
(representing repressive interactions) we used the indicated bounds with a
negative sign. By default, the upper bound for the regulatory strength was
$H_{\rm max}=50$, but reduced to smaller values in some optimization runs (see
Sec. 8.3). In the lower part we show how the total number of optimized
parameters (including $D$) scales with the number of gap genes simulated,
$N_{\rm g}$; these numbers reduce by $2(N_{\rm g}-2)$ when the constraint that
only two genes can couple to the posterior gradient is imposed, as explained
in Sec. 5.2.

###  5.3 Resource constraints

One of the most important constraints impacting on the positional information
encoded by a noisy gap gene expression pattern is the total number of
molecules that establish the pattern. If this number could be made arbitrarily
high, we could, in principle, remove all intrinsic noise in the system and
attain the maximal positional information possible. Note that in our discrete
model this maximal value is given by the logarithm of the number of nuclei
along the embryo axis, $I_{\rm max}=\log_{2}N_{x}\simeq 6~{}{\rm bits}$ – this
bound corresponds to a perfect (error-free) identifiability of each nucleus
based on a single readout of local gap gene expression patterns. In our model,
the constraint on molecular copy numbers is imposed on two levels.

Firstly, we limit the maximal mRNA and protein copy numbers, $M_{\rm m}$ and
$N_{\rm m}$, that can be reached at full induction of gene expression at each
position by imposing the maximal production rates and burst size indicated in
Table 2. These values were chosen such that (1.) maximal mRNA copy numbers
approximately correspond to experimentally determined values in nuclear cycle
14 [12], and (2.) maximal protein copy numbers correspond to current order-of-
magnitude estimates ($10^{3}-10^{4}$ proteins at maximal induction per
nucleus). With this choice, the gene expression variance in optimized systems
turns out to be close close to the experimentally measured variance, which we
computed from existing experimental data [9].

Secondly, we impose a constraint on the “resource utilization” (RU), which we
define to be the fraction of the total maximal protein number throughout the
whole system (i.e., summing over all positions) that is available for creating
the pattern (see Sec. 6.3 for the corresponding formula). ${\rm RU}<1$ means
that the pattern in which all genes are fully induced at all positions cannot
be created, and the lower RU, the more positions in the embryo must have at
least some gap genes below maximal induction. By default, we set the resource
utilization constraint to be equal to the corresponding value in the measured
normalized mean expression pattern of Drosophila, RU = 0.2, which allows us to
compare theoretically optimal solutions under the same molecular cost
constraint that the fly embryo experiences.

Note that the experimental value of RU is well below RU = 0.5. RU = 0.5 would
lead to shifted “counter patterns” in which the genes are expressed at 50% of
the available positions. Such patterns are indeed optimal for a binary code
[13], but in wild-type Drosophila the RU is significantly lower, since the
embryo can employ a richer code that is not restricted to binary “ON”/“OFF”
expression levels and which can productively use intermediate expression
levels (also see Sec. 8.2). This is shown in Fig. 3C in the main text, where
we explored higher and lower values of RU in dedicated optimization runs, and
how this influences the positional information values reached by optimal
solutions. Note that even when the constraint on resource utilization is
removed, the positional information values of the optimized patterns still
remain $\lesssim 5~{}{\rm bits}$. The only way to increase positional
information beyond this bound would be to allow higher protein copy numbers
(i.e., higher mRNA and protein production rates) – this truly is a fundamental
biophysical constraint that can only be removed at a higher time- or metabolic
cost.

The resource utilization (RU) constraint is implemented in the optimization
runs by simply rejecting simulated patterns that exceed the imposed RU value.

###  5.4 Pattern stability constraint

When exploring the role of the inter-nuclear diffusion constant $D$ (see Sec.
8.1), we observed that at lower $D$ optimal solutions tend to compensate for
the loss of spatial averaging by producing patterns that have more wiggled
spatial profiles–while at the same time exhibiting a larger amount of non-
stationarity of patterns at read-out time $T_{\rm meas}$. As explained in the
main text, achieving pattern stability around the read-out time might be
essential to the patterning system architecture, motivating us to explore
optimizations with an upper bound on the allowed pattern rate-of-change (RoC).
This constraint was imposed in the same way as the constraint on resource
utilization, i.e. by rejecting patterns that exceed the imposed maximal RoC
value.

We defined the pattern rate-of-change as the temporal derivative of the
normalized gap protein expression levels evaluated at the read-out time
$T_{\rm meas}$ and averaged over all positions and gap gene species, i.e.,:

$\displaystyle{\rm RoC}(\mathbf{g},T_{\rm
meas})=\left\langle\left\langle\frac{dg_{i}^{\alpha}(T_{\rm
meas})}{dt}\right\rangle_{\alpha}\right\rangle_{i}$ (79)

where $\alpha$ runs over $1..N_{\rm g}$ and $i$ over $1..N_{\rm x}$.

##  6 Details of data analysis

###  6.1 Computing positional error profiles

The local value of the positional error, $\sigma_{x}(x_{i})$, was computed
from the first and second moments of the posterior decoding-map distributions.
Recall that the average decoding map is $\langle P\rangle(x_{j}^{*}|x_{i})$,
which is defined in the section on positional information, Sec. 3. The
formulae for the moments and the positional error therefore read:

$\displaystyle M_{1}(x_{i})$ $\displaystyle=\frac{1}{N_{\rm
x}}\sum_{j}x_{j}^{*}\langle P\rangle(x_{j}^{*}|x_{i})$ (80) $\displaystyle
M_{2}(x_{i})$ $\displaystyle=\frac{1}{N_{\rm x}}\sum_{j}(x_{j}^{*})^{2}\langle
P\rangle(x_{j}^{*}|x_{i})$ (81) $\displaystyle\sigma_{x}(x_{i})$
$\displaystyle=\sqrt{M_{2}(x_{i})-M_{1}^{2}(x_{i})}.$ (82)

###  6.2 Determining number of slopes

In order to numerically determine the number of slopes we first computed (per
gene) the numerical derivatives of the normalized expression profiles with
respect to space, and then determined local extrema by comparing the sign of
successive derivative values
$(\bar{g}^{\alpha}_{i})^{\prime}\equiv\frac{d\bar{g}^{\alpha}}{dx}|_{x=x_{i}}=\frac{\bar{g}^{\alpha}_{i+1}-\bar{g}^{\alpha}_{i}}{\Delta}$.
In order to avoid falsely calling small local fluctuations full profile
slopes, we restricted the test for derivative slope change to positions with
significantly elevated gene expression. Moreover, only absolute derivative
values surpassing a required minimum were considered to constitute “slopes.”
This excluded, e.g., patterns in which the expression level reduced only
slightly along the embryo axis and never stretched from very high to very low
levels. Since the derivative is a local quantity defined between two
successive lattice points, we averaged both the expression levels and
derivative values before testing whether they are above the required minimal
values. Mathematically, these additional conditions read:

$\displaystyle(\bar{g}_{i-1}^{\alpha}+\bar{g}_{i}^{\alpha})/2$
$\displaystyle\geq\bar{g}_{\rm min}$
$\displaystyle(|(\bar{g}_{i-1}^{\alpha})^{\prime}|+|(\bar{g}_{i}^{\alpha})^{\prime}|)/2$
$\displaystyle\geq\bar{g}^{\prime}_{\rm min}$ (83)

In practice, $\bar{g}_{\rm min}=0.1$ and $\bar{g}^{\prime}_{\rm min}=1/EL$
(corresponding to a change from full to zero expression over the whole embryo
length) proved suitable choices. Overall this simple scheme has proven to be
very robust in discarding occasional small random bumps in the profiles while
calling only larger changes in expression levels significant slopes.

###  6.3 Measuring resource allocation

For any given expression pattern, we define its “resource utilization” as the
average expression level across all gene species and positions,

$\displaystyle RU\equiv\frac{1}{N_{\rm g}N_{\rm
x}}\sum_{\alpha}\sum_{i}\bar{g}_{i}^{\alpha}\quad,$ (84)

where the mean expression levels $\bar{g}_{i}$ already are normalized between
0 and 1, as explained in Sec. 1.1; the prefactor $\frac{1}{N_{\rm g}N_{\rm
x}}$ ensures that the resource utilization is also normalized between 0 and 1,
with 1 corresponding to the situation in which all genes are expressed with
the maximal expression rate at all positions.

The partial “resource utilization per gene” is defined as

$\displaystyle RU^{\alpha}\equiv\frac{1}{N_{\rm g}N_{\rm
x}}\sum_{i}\bar{g}_{i}^{\alpha}\quad,$ (85)

where we retain the normalization factor $\frac{1}{N_{\rm g}N_{\rm x}}$ for
comparability with $RU$; notice that the coefficient of variation of this
quantity reported for different patterns in the main article (Fig. 3A) is
independent of the specific choice of the normalization prefactor.

Similarly, we define the “resource utilization per position” as:

$\displaystyle RU_{i}\equiv\frac{1}{N_{\rm g}N_{\rm
x}}\sum_{\alpha}\bar{g}_{i}^{\alpha}$ (86)

Here we locally sum the expression levels of all gene species.

Both $RU^{\alpha}$ and $RU_{i}$ define discrete distributions for which we
compute the coefficient of variation (CV), by normalizing the standard
deviation of the distributions by their respective means (see Fig. 3A and SI
Fig. 5C).

The resource bias shown in SI Fig. 5D was computed according to the formula

$\displaystyle B=\frac{R_{A}-R_{P}}{R_{A}+R_{P}}$ (87)

where $R_{A}$ is the total number of proteins expressed in the anterior half
of the embryo, and $R_{P}$ the remaining total number of proteins expressed in
the posterior half, respectively. For an uneven number of nuclei along the
embryo axis, the resources of the central position were assigned to the
anterior ($R_{A}$) and posterior ($R_{P}$) with weight 1/2, respectively.

###  6.4 Determining regulatory interactions

We visualized the regulatory interactions of the optimized patterns in two
ways. Firstly, we computed the local “regulatory contributions” by inserting
the local concentrations of maternal inputs into the free-energy-like terms
that enter the exponent of the MWC regulatory function (see Sec. 1.5), i.e.,
$H_{M}^{\alpha\kappa}\log(1+c_{i}/K_{M}^{\alpha\kappa})$ for maternal and
$H_{G}^{\alpha\zeta}\log(1+\bar{g}_{i}/K_{G}^{\alpha\zeta})$ for zygotic
interactions. This resulted in detailed spatial profiles of the regulatory
inputs into each of the simulated genes along the embryo axis. Stacked example
profiles for our best Drosophila-like solution (Fig. 2 of the main article)
are shown in SI Fig. 4.

Secondly, in order to obtain an even simpler and stereotyped representation in
terms of “regulatory arrows,” which ignores spatial dependency of the
regulatory contributions, we reconstructed regulatory network “cartoons” from
the zygotic regulation parameters ($H_{G}^{\alpha\zeta}$ and
$K_{G}^{\alpha\zeta}$), by evaluating the logarithmic terms (see above) at the
highest possible zygotic expression levels (which in our model all are equal
to 1). Terms resulting in values larger than 10 were deemed to be strong
interactions depicted by fat arrow symbols in the cartoons; terms resulting in
values smaller than 1 were deemed insignificant and therefore not shown in the
cartoons. Networks visualized in this way are shown in SI Fig. 6 and Fig. 2D
of the main article.

##  7 Additional analyses of the optimal ensemble

In addition to quantities discussed in Fig. 3A and B of the main text, we
extended our comparison between the optimal ensemble and the random ensemble
to further quantities introduced below; comparison results are shown in SI
Fig. 5.

SI Fig. 5A shows the pattern rate-of-change (RoC) distributions for the
solutions of the optimal (graded colors) and random (grey) ensembles; for the
optimal ensemble, the data shown is the same as used Fig. 3B of the main text.
As in the main figure, SI Fig. 5A shows that higher positional information
(yellow bullets) tends to lead to more stable patterns, i.e., a slower pattern
RoC. Interestingly, the RoC is distributed over a similar range in the random
ensemble. This is due to the fact that a considerable fraction of the random
ensemble consists of solutions with little or no cross- and self-regulation,
which tend to have slow dynamics (but also carry little information).

Since we observed that for some solutions the RoC itself may either increase
or decrease as the patterns build up, we asked whether such “RoC slowdown” or
“RoC speedup” is correlated with the ability to encode more positional
information in the pattern. We defined “RoC slowdown” as the ratio between the
RoC at 30 min before the readout time $T_{\rm meas}$, i.e. at
$t=135~{}\text{min}$, and the RoC at $T_{\rm meas}=165~{}\text{min}$; in that
respect, a RoC slowdown smaller than one is a RoC speedup. Our RoC estimate is
thus a discrete approximation to the second-derivative of the gap gene
expression dynamics.

We observe that high positional information tends to correlate with a dynamics
slowdown when the patterns progress towards $T_{\rm meas}$; as shown in SI
Fig. 5B, the optimal patterns tend to reduce their RoC by a factor of 1.6 on
average during the last 30 min of patterning. Comparison to random patterns is
not trivial. The random patterns also tend to slow down, as a consequence of a
generic anticorrelation between the RoC slowdown and the RoC itself; in other
words, patterns with a small RoC at $T_{\rm meas}$ tend to have a larger
relative reduction of the RoC before. Since random patterns tend to have a
smaller RoC in general, for the reasons mentioned in the last paragraph, they
also tend to display a more pronounced slowdown. However, while for the
optimal solutions the RoC slowdown is positively correlated with the
positional information of the pattern (${\rm PC}\simeq 0.4$, $p<10^{-13}$),
such correlation is very weak (and has an inverted sign) for random patterns
(${\rm PC}\simeq-0.12$, $p=0.026$).

Panel C of SI Fig. 5 shows, for the standard optimal ensemble and the random
ensemble, how uniformly the patterns tend to distribute the available
resources (gap gene proteins) along the embryo axis, quantified by the average
(across-gene) coefficient of variation (CV) of protein distributions, computed
across the length of the embryo. Here the median CV is significantly higher
(0.46) for the optimized patterns as compared to patterns from the random
ensemble (0.13). This reflects the fact that upon optimization the variety of
combinations of locally expressed genes increases, meaning that locally more
resources have to be devoted in order to obtain a larger number of significant
signals (expression levels). In contrast, randomly generated patterns without
significant cross-regulation tend to locally express one, or rarely at most
two genes simultaneously, which reduces the variance in the necessary
resources.

In SI Fig. 5D we analyse how the available protein resources are distributed
between the anterior and posterior half of the embryo. The “resource bias” is
defined such that a value of 0 means equal distribution of resources between
the anterior and posterior halves of the embryo, while a value of $\pm 100$
implies a complete attribution of resources to the anterior (+) or posterior
(-) half, respectively (see Sec. 6.3 for the corresponding formula). We
observe that optimized solutions, on average, tend to attribute somewhat more
resources (bias of $14\%$) to the anterior half, in part because posterior-
most 5 nuclei of the system, as described in Sec. 5.1, do not contain gap gene
patterns. This effectively leads to a stronger anterior maternal signal and,
consequently, optimal solutions tend to allocate more resources there.

##  8 Optimizations with altered system components

In addition to the standard optimized ensemble of Fig. 3 in the main text, we
also carried out _ab-initio_ optimizations for hypothetical alternative
settings in which important components and parameters held fixed in the
standard ensemble were changed. In particular, we varied: (1.) the inter-
nuclear diffusion constant $D$; (2.) the number of available gap genes
expressed downstream of maternal inputs (1–5 genes); (3.) the number of
maternal input gradients regulating the gap genes, i.e. systems with all three
maternal gradients present vs. systems in which one or two of these are
lacking; (4.) the maximal regulatory strength for maternal and zygotic
regulatory interactions.

In the subsections below we describe these four altered optimiziation settings
and their results.

###  8.1 Varying gap protein diffusion constant

In our standard optimization ensemble we used a fixed value of the inter-
nuclear diffusion constant of $D=0.5~{}\mu{\rm m}^{2}/{\rm s}$ for gap gene
proteins, which is in the range reported for Bcd [14]. This choice generates
mean expression profiles that agree very well with the experimental data. We
later systematically relaxed this choice and carried out additional
optimizations in which (1.) other fixed values of $D=10^{-2}~{}\mu{\rm
m}^{2}/{\rm s}~{}\cdots~{}10~{}\mu{\rm m}^{2}/{\rm s}$ were used, and (2.) $D$
itself was considered an additional optimization parameter and optimized for
maximizing positional information jointly with the regulatory parameters. In
both cases, the additional optimizations were carried out both without and
with a bound on the maximal pattern rate-of-change (RoC).

The corresponding results are summarized in Fig. 3E and Fig. 3F of the main
text, where in panel E we show how the mean positional information $I$ and the
mean pattern RoC vary with the imposed diffusion constant $D$, while panel F
shows example expression patterns for two largely different values of $D$. For
the cases in which $D$ itself is optimized, the $I$ and RoC values are plotted
at the single value of $D$ which is computed as the average optimal $D$ over
optimization runs.

In accordance with previous findings [1], the positional information of
optimal expression patterns decreases markedly for high diffusion constants,
$D\gtrsim 1~{}\mu{\rm m}^{2}/{\rm s}$. In this regime, though spatial
averaging is very efficient in removing non-Poissonian noise components, the
gap gene expression patterns become increasingly spread-out and spatially
flatter, such that a smaller number of slopes (or combinations thereof) can be
accommodated along the embryo axis. This markedly reduces the information
content of the “positional code”; in short, the ability to create sufficiently
(but not infinitely!) steep boundaries is essential for good patterning, and
high $D$ values prevent the such boundaries from being set up. An example of a
suboptimal pattern at large diffusion constant values is shown in the lower
part of Fig. 3F of the main text.

Interestingly, in the opposite regime (and without any bound on the pattern
RoC), optima with comparably high positional information ($I\gtrsim 4~{}{\rm
bits}$) are found even for very low values of $D\sim 10^{-2}~{}\mu{\rm
m}^{2}/{\rm s}$. This is due to the fact that in this regime, patterns with
very wiggled (high spatial derivatives), irregular shapes can be created,
which nominally increases positional information. However, this comes at the
cost of decreasing pattern stability, resulting in a high rate-of-change (see
red curve in lower part of Fig. 3E of the main text). Accordingly, when the
bound on the pattern RoC is also imposed, the optimal solutions at low $D$ do
not reach the plateau of high information any more, and an optimal regime
emerges at intermediate values of $D\simeq 0.1~{}\mu{\rm m}^{2}/{\rm
s}~{}\cdots~{}1~{}\mu{\rm m}^{2}/{\rm s}$. As expected, this optimal regime is
found around the average value of optimized $D_{\rm opt}\simeq 0.4~{}\mu{\rm
m}^{2}/{\rm s}$, which happens to be very close to our originally chosen fixed
value of $D_{\rm fix}=0.5~{}\mu{\rm m}^{2}/{\rm s}$.

###  8.2 Altered number of gap genes and maternal inputs

An intriguing question to ask is why the Drosophila system features exactly
four gap genes expressed as the first layer of zygotic downstream genes
activated by maternal inputs. Could evolution have realized a similar
patterning precision with a smaller number, e.g. three, gap genes?
Alternatively, could there be any significant benefit of using more than four
gap genes? Similarly one can ask whether all three of the maternal input
systems (anterior gradient, posterior gradient, and terminal system) are
necessary for generating patterns with high positional information (and if so,
why or why not).

To address these questions, we carried out hundreds of optimization runs with
altered numbers of downstream (gap) genes and different combinations of
maternal input gradients present in the system (“APT” = all present, “AP” =
anterior and posterior gradients only, “AT” = anterior gradient and terminal
system only, “A” = anterior gradient only).

In the runs with altered numbers of gap genes, we disallowed the expression of
more than the prescribed number of genes, $N_{g}<4$, and in the $N_{g}=5$ case
we introduced one additional hypothetical gap gene species for simulating the
expression of five-gene patterns, allowing the new gene to freely couple to
the available morphogen gradients and other gap genes. In order to compare the
optimized ensembles at equal resource utilization, the respective bounds on
resource utilization were scaled with the number of gap genes, such that the
total amount of expressed molecules was the same for all gene numbers, i.e.:
for 3 genes, the resource utilization bound was scaled up by factor 4/3, for 5
genes it was scaled down by factor 4/5, and so forth for the other gap gene
numbers. In short, we imagine that there is a fixed total (across gap genes)
budget of gap gene expression, which can be split across different numbers of
gap gene species: with fewer different gap genes, each can be expressed at
more positions in space; with more different gap genes, each must be expressed
in a more localized setting. In all these runs, the hard biophysical
constraints (i.e., the gap gene product degradation rates and maximal
transcription and translation rates) were held fixed.

In the optimization runs with reduced number of maternal inputs the
concentrations of the respective gradients “knocked out” were simply set to
zero, and the optimization parameters pertaining to the coupling of downstream
genes to these gradients were taken out of the set of optimization parameters.

Early in the analysis of alternative optimal ensembles we realized that the
variety of optimal solutions markedly decreases with decreasing number of
available downstream genes, and also with decreasing number of maternal
gradients. This, in itself, is an interesting finding. We therefore carried
out a pattern similarity analysis by computing the mutual distances between
the optimal patterns found in the respective ensemble (using the distance
function defined in Eq. (78) of Sec. 4.4), and then resorting the distance
matrix by hierarchical clustering (via standard MATLAB routines).

We summarize the results of the similarity analyses in SI Fig. 8 and SI Fig.
9. In these figures, we show the sorted distance matrices for the combinations
of gene numbers and maternal input gradients (as indicated in the titles above
the matrices). The number of expressed downstream genes increases along the
rows. The left column panels all feature the full set of maternal inputs
(“APT”), whereas the right column panels show optimizations with the anterior
input gradient only (“A”). Selected segments of the distance similarity
matrices in which mutual distances are low – which define “clusters” of
similar patterns – are highlighted by red boxes in the matrices; the boxes
next to the matrices linked to these clusters show several randomly chosen
example patterns from the respective clusters.

We observe the lowest variety of optimal patterns in the case with one
downstream gene only (SI Fig. 8A and B), although that variety is not limited
to one single optimal pattern when all three maternal inputs are present (“1
GENE / APT”, SI Fig. 8A). We find two optimal shapes, symmetric along the
anterior-posterior axis, but, surprisingly, also one solution which prefers to
create one bump shallowly increasing towards the posterior and then sharply
decreasing again. This can be attributed to the fact that in our system the
posterior gradient reaches slightly lower maximal levels due to the “cut-off”
described in Sec. 5.1. In contrast, the situation with one gap gene and the
anterior input gradient only (“1 GENE / A”, SI Fig. 8B) is simpler: here we
identify one single cluster of stereotypical solutions, in accordance with our
previous theoretical findings [3, 1]. Notice that while the optimal patterns
with one downstream gene nominally are permitted to be fully expressed at 80%
of the positions (to match the resource utilization of the WT), these patterns
in fact prefer to use less of the available resources in order to create more
varying pattern shapes.

In sum, stochastic optimization runs for one gap gene converge on a single
solution (for one morphogen input) or a pair of qualitatively different
solutions (for three morphogens inputs), showing that optimizing positional
information in this setup strongly constraints the space of optimal expression
patterns.

Solutions optimized with two gap genes show the increase in the variety of
optimal solutions, but we can still clearly identify structural similarities
(SI Fig. 8C and D). Optimal solutions for two gap genes responding to a single
maternal morphogen form a small number of discrete clusters (SI Fig. 8D, “2
GENES / A”), which are furthermore qualitatively consistent with our previous
work [3, 4]; the stochastic optimization repeatedly finds multiple solutions
belonging to these stereotyped clusters. Amongst all these solutions, the
available resources are predominantly assigned to the anterior half of the
embryo.

For optimal solutions with two gap genes responding to three maternal
morphogens, in addition to almost symmetric solutions in which the two genes
are expressed in opposite halves of the embryo, we also found asymmetric
solutions in which the resources are assigned predominantly either to the
anterior or posterior of the embryo (see example patterns shown in SI Fig.
8C). Notice that in these solutions one of the two gap genes follows a profile
that constitutes an optimal solution for the optimizations with one gene (SI
Fig. 8C), while the other gap gene profile is “added on top.”

The observed trend of increasing solution diversity continues as the number of
available gap genes is increased beyond two. With three gap genes and all
three maternal morphogens, there are more clusters, each containing more
degenerate solutions, and the optimal patterns within each cluster are also
less similar to each other (SI Fig. 8E). Patterns remain more stereotypic and
similar to each other, though, within the clusters in the optimal ensemble
with 3 gap genes and anterior maternal morphogen only (SI Fig. 8F).

With 4 and 5 gap genes and all three maternal gradients we still can sometimes
find clusters of structurally similar solutions, but as a rule the optimal
solutions now are dissimilar (SI Fig. 9A and C). Once again the clusters of
similar solutions are larger with one maternal gradient only, and the
solutions within these clusters also are more similar to each other (SI Fig.
9B and D). Interestingly, some of these solutions do not make use of all
available genes but prefer to assign all available resources to a smaller set
of genes (i.e., the pattern consists only of nonzero expression for 3 out of
4, or 4 out of 5, available genes).

Overall, the similarity between the patterns decreases less when going from 3
to 4 genes as compared to going from 2 to 3 genes. This is best seen in the
comparison of the pattern distance histograms across the ensembles with
different gene numbers, SI Fig. 9E and F (upper panels). Here the ensembles
with 1 or 2 gap genes extend out towards very low pattern distances, while
starting from 3 genes a characteristic distance distribution forms at larger
pattern distances. The data for 5 genes was omitted from these plots because
it overlaps almost perfectly with the data for 4 genes. Below the distance
histograms we also show the corresponding histograms of positional information
encoded by the optimized patterns. Here the largest jumps occur between 1 and
2, and between 2 and 3 gap genes. 3-gap-gene patterns can reach surprisingly
high values of positional information (maximal value $I_{3}^{\rm
max}=4.56~{}{\rm bits}$), but 4-gap-gene patterns still reach about 0.2 bits
more ($I_{4}^{\rm max}=4.77~{}{\rm bits}$); the increase between 4- and 5-gene
patterns is smaller ($I_{5}^{\rm max}=4.84~{}{\rm bits}$). Note that while the
nominal increases of these maximal $I$-values may seem minuscule, due to the
logarithmic nature of the information-theoretic quantities, these information
values correspond to an increase from 23 to 27 (and 28, respectively)
perfectly distinguishable states (positions).

###  8.3 Bounding maximal regulatory strength

We also asked about how the maximal allowed strength of regulatory
interactions affects the nature of optimal solutions. By default, we set the
bounds on the the maximal regulatory interaction strength to be very generous,
$H_{\rm max}=50$. As we found in our earlier work for a single output gene
under the control of a single maternal input gradient (but using the same type
of regulatory function, see Sec. 1.5), this maximal value is high enough to
produce a sharp, step-like transition of the activation function from full to
zero expression within two successive nuclei along the anterior-posterior axis
[1].

In order to assess the influence of this generous choice, we also generated
optimization ensembles with lower values of $H_{\rm max}\in[3,4,5,10]$.
Moreover, we extended this analysis to two additional hypothetical
restrictions on the possible regulatory interactions: (1.) regulation with
different $H_{\rm max}$ and anterior gradient only (“A only”); (2.) “feed-
forward regulation only” (“FF only”) at different $H_{\rm max}$, meaning that
zygotic regulatory interactions, i.e. mutual regulation and self-regulation,
are not possible at all, while maternal regulation is constrained by lowering
$H_{\rm max}$. For the highest value $H_{\rm max}=50$, we also generated the
“FF only” optimal ensembles with two (anterior + posterior = “AP”, and
anterior + terminal = “AT”) gradients.

Fig. 4B and C of the main text summarize how these extra constraints impact
the positional information encoded by optimized patterns. While in all three
hypothetical settings the positional information decreases with decreasing
maximal strength of regulatory interactions $H_{\rm max}$, mutual and self-
regulation can partly compensate for this decrease compared to the “FF only”
case. This is particularly pertinent at very low $H_{\rm max}=3$.
Qualitatively, the same conclusion holds with one morphogen input gradient