# Scale invariance in early embryonic development

Miloš Nikolić,a,b Victoria Antonetti,a,c Feng Liu,b,c Gentian Muhaxheri,a,d
Mariela D. Petkova,e Martin Scheeler,a Eric M. Smith,a William Bialek,a,b,f
and Thomas Gregora,b,g aJoseph Henry Laboratories of Physics and bLewis–Sigler
Institute for Integrative Genomics, Princeton University, Princeton NJ 08544
USA
cCenter for Quantitative Biology and School of Physics, Peking University,
Beijing 100871 China
dDepartment of Physics, Lehman College, City University of New York, Bronx, NY
10468 USA
eProgram in Biophysics, Harvard University, Cambridge MA 02138 USA
fInitiative for the Theoretical Sciences, The Graduate Center, City University
of New York, 365 Fifth Ave., New York, NY 10016 USA
gDepartment of Developmental and Stem Cell Biology UMR3738, Institut Pasteur,
75015 Paris, France

###### Abstract

The body plan of the fruit fly is determined by the expression of just a
handful of genes. We show that the spatial patterns of expression for several
of these genes scale precisely with the size of the embryo. Concretely,
discrete positional markers such as the peaks in striped patterns have
absolute positions along the anterior–posterior axis that are proportional to
embryo length, with better than $1\%$ accuracy. Further, the information (in
bits) that graded patterns of expression provide about position can be
decomposed into information about fractional or scaled position and
information about absolute position or embryo length; all of the available
information is about scaled position, again with $\sim 1\%$ accuracy. These
observations suggest that the underlying genetic network exhibits scale
invariance in a deeper mathematical sense. Taking this mathematical statement
seriously requires that the network dynamics have a zero mode, which connects
to many other observations on this system.

##  I Introduction

Closely related organisms can vary widely in size, but variations in their
proportions are much smaller [1, 2, 3]. There is a considerable gap between
this qualitative observation and some precise mathematical statement of
scaling, e.g., that the linear dimensions of all elements in the body plan are
in direct proportion to the linear dimensions of the organism. If correct this
scale invariance would be visible not only in the fully developed organism but
already at some earlier stages in its development.

There are many examples of “allometric scaling,” power-law relationships among
different quantities across a well-defined class of organisms [4, 5, 6]. In
some cases, these relations connect the linear dimensions of different body
parts. Nonetheless, truly precise spatial scaling in embryonic development
would be quite surprising.

We understand the mechanisms of pattern formation in a wide range of
non–biological systems, from fluid flows to crystal growth (snowflakes) and
more [7, 8, 9, 10, 11, 12], but none of these examples exhibit scale
invariance. Instead, the elements of the pattern have linear dimensions set by
microscopic parameters, and larger systems exhibit more repetitions of the
same pattern rather than expansion or contraction of pattern elements to match
the size of the system as a whole [13]. Going back to the pioneering work of
Turing [14], the mathematical structure of the equations governing these
systems is not so different from the structure of models for genetic or
biochemical networks. If we take these analogies literally, we would predict
that taller people should have more vertebrae, which is obviously wrong. Is
there a real problem here, or are we just oversimplifying?

Here we try to make the notion of scale invariance in development more
precise. We use the first hours of development in the fruit fly as an example,
following spatial patterns of morphogen concentration as they flow through
three layers of a genetic network, from maternal inputs to the gap genes to
the pair rule genes [15, 16, 17]. In the spirit of earlier work [18, 19, 20,
21, 22] we analyze discrete positional markers, such as the stripes in pair
rule-gene expression, and find that positions of these markers vary in
proportion to the length of the embryo with better than $1\%$ accuracy [23].
We then go beyond discrete markers, decomposing the information carried by
graded patterns of gap gene expression into information about fractional or
scaled position vs. information about the absolute position; we find that all
the available information is about fractional position along the
anterior–posterior axis. Information that would signal a deviation from scale
invariance is less than $1\%$ of the total.

These results provide strong evidence for scaling in a precise mathematical
sense, for both the gap genes and the pair rule genes. But at least one of the
maternal inputs, Bicoid [24, 25], does not show any sign of scale invariance:
as in well-understood non-biological pattern-forming systems, there is a
length scale that presumably is set by underlying molecular parameters and
does not adjust in response to the linear dimensions of the whole embryo. This
suggests that scale invariance is an emergent property of the gap gene
network.

We argue that true scale invariance places very specific requirements on the
dynamics of this network, independent of molecular details: it must have a
“zero mode.” This has connections to other observations on gap gene dynamics
[26, 27] and to more detailed models [28, 29].

##  II Testing for scaling

To make the notion of scaling more precise we take seriously the idea that
cell fates are determined by the concentration of particular molecules called
morphogens [30]. Since the cell fates are tied to their positions, the
concentrations of morphogens must also carry information about position along
the body axes. These ideas are especially crisp in the early fly embryo, where
we know the identities of all the relevant morphogens and rich spatial
patterns in the concentrations of these molecules are established before cells
make large-scale movements [31].

We focus on pattern formation along a single axis, which will be the
anterior–posterior axis in the analysis of fly embryos below. Then we can
measure position along this axis by a single variable $0<x<L$, where $x=0$ is
the anterior end of the embryo, $x=L$ is the posterior end, and hence $L$ is
the length of the embryo. There are multiple morphogen species, indexed by
$\rm i$, and if we neglect the discreteness of cells then their concentration
profiles are described by continuous functions $g_{\rm i}(x;L)$. The notation
emphasizes that concentration profiles may be different in embryos of
different size $L$.

True scale invariance is the statement that the concentration of morphogens
depends only on position relative to the length, that is

$g_{\rm i}(x;L)=\Phi_{\rm i}(x/L).$ (1)

If there is a fixed map from morphogen concentrations to cell fates, this
scaling behavior would guarantee that cells adopt a fate that depends on their
relative position $x/L$, and not separately on $x$ and $L$.

How do we test for scale invariance? If the concentration of the morphogen has
a single peak as a function of $x$, we can write

$g_{\rm i}(x;L)=g(x-x_{\rm p};L),$ (2)

then scale invariance as in Eq. (1) requires that all the $L$ dependence is
contained in the position of the peak.

$x_{\rm p}=\langle f_{\rm p}\rangle\cdot L+{\rm noise};$ (3)

where $f_{p}$ is the fractional or scaled peak position,
$\langle\cdots\rangle$ is the average over many embryos with different
lengths, and $\rm noise$ allows that positions jitter from embryo to embryo.
We emphasize that Eq. (3) is not just the statement that positional markers
adjust (in absolute distance) to the length of the embryo; scale invariance as
we have defined it in Eq. (1) requires that this adjustment is exactly linear
with zero intercept. There is a natural generalization to concentration
profiles that have multiple peaks, as with the pair rule genes (Fig. 1A, B).

It has been known for some time that the morphogens in the early fly embryo
carry enough information to specify scaled positions with $\sim 1\%$ precision
all along the anterior–posterior axis [32, 33]. At the same time, embryos from
the same mother, in an inbred laboratory stock, fluctuate in length with a
standard deviation of $\sigma_{L}/\langle L\rangle\sim 4\%$ (Appendix A and
[34, 35]). It would seem that to make these numbers consistent with one
another, positional signals must scale with embryo length, but this is a bit
subtle.

Imagine a hypothetical embryo in which, e.g., the peak of the morphogen
profile is perfectly anchored in absolute position relative to the anterior
pole of the embryo, with no scaling and no noise, such that $x_{\rm p}=\langle
x_{\rm p}\rangle$. Then the relative or fractional positions $f_{\rm p}=x_{\rm
p}/L$ fluctuate only because the lengths of the embryos vary,

$\displaystyle\sigma_{f_{\rm p}}^{2}(A)$ $\displaystyle\equiv$
$\displaystyle\langle(\delta f_{\rm p})^{2}\rangle=\langle x_{\rm
p}\rangle^{2}\left[\left\langle\left(\frac{1}{L}\right)^{2}\right\rangle-\left\langle\frac{1}{L}\right\rangle^{2}\right]$
(4) $\displaystyle\sim$ $\displaystyle\left({{\langle x_{\rm
p}\rangle}\over{\langle L\rangle}}\cdot{{\sigma_{L}}\over{\langle
L\rangle}}\right)^{2}.$ (5)

Thus for a marker that on average is a quarter of the way from the anterior to
posterior, $\langle x_{\rm p}\rangle=0.25\langle L\rangle$, fluctuations will
be $\sigma_{f_{\rm p}}(A)\sim 0.01$ even without scaling. Similarly, if we
have a marker anchored at some fixed absolute position relative to the
posterior then the variance in fractional position will be

$\sigma_{f_{\rm p}}^{2}(P)=\left(1-{{\langle x_{\rm p}\rangle}\over{\langle
L\rangle}}\right)^{2}\cdot\left({{\sigma_{L}}\over{\langle
L\rangle}}\right)^{2}.$ (6)

We can imagine cells combining anterior and posterior signals to reduce the
error,

${1\over{\sigma_{f_{\rm p}}^{2}(A,P)}}={1\over{\sigma_{f_{\rm
p}}^{2}(A)}}+{1\over{\sigma_{f_{\rm p}}^{2}(P)}}.$ (7)

With $\sigma_{L}/\langle L\rangle\sim 0.04$, fluctuations in fractional
position thus could be less than $\sim 1.4\%$ everywhere along the
anterior–posterior axis, even in the absence of any scaling mechanism.
Convincing ourselves that pattern formation is truly scale invariant requires
a very precise measurement and depends on the system itself being very
precise.

It is intuitive to think about scaling as the proportionality of positions to
embryo length, as in Eq. (3), but it should be possible to test the scaling of
the entire morphogen profile, as in Eq. (1), more directly. There are two
related observations. First, to compare morphogen profiles across embryos of
different lengths, we need a metric. Second, since morphogen profiles are
noisy, it is unrealistic to expect the exact equality of two functions across
all values of $x$. Fortunately, the noise level itself provides a metric for
comparison, which is made precise in the language of information theory.

The statement that morphogen profiles depend on $x$ and $L$ means that the
concentrations of these molecules provide information about the underlying
positional variables. This is quantified, uniquely, by the Shannon information
[36, 37, 38],

$I({\mathbf{g}}\rightarrow\\{x,L\\})=\int d{\mathbf{g}}\int dx\int
dL\,P\left({\mathbf{g}}|\\{x;L\\}\right)P(x,L)\log_{2}\left[{{P\left({\mathbf{g}}|\\{x,L\\}\right)}\over{P\left({\mathbf{g}}\right)}}\right]\,{\rm
bits},$ (8)

where for more compact notation we write ${\mathbf{g}}=\\{g_{\rm i}\\}$ and
$d{\mathbf{g}}=\prod_{\rm i}dg_{\rm i}$. Here
$P\left({\mathbf{g}}|\\{x,L\\}\right)$ is the probability of finding the set
of morphogen concentrations $\\{g_{\rm i}\\}$ at position $x$ in an embryo of
length $L$; $P\left({\mathbf{g}}\right)$ is the probability of finding these
concentrations averaged over all values of $x$ and $L$; and $P(x,L)$ is the
distribution of positions and lengths. It is useful to recall that this
information is mutual: the concentrations of morphogens provide cells with
information about position, and specifying position allows us to predict the
concentrations, so we write $I({\mathbf{g}};\\{x,L\\})$. Information depends
on both the mean spatial profiles of the morphogens and their noise levels.

True scale invariance would mean that all of the information conveyed by the
morphogens is about the fractional position $x/L$:

$I({\mathbf{g}};\\{x,L\\})=I({\mathbf{g}};x/L)\,\,\,\,{\rm(perfect\
scaling)}.$ (9)

Equivalently, if we want to predict the morphogen concentration, it is enough
to specify the fractional position, and no extra information is gained by
knowing $x$ and $L$ separately. We can think of the total information as
having a component about the relative position and an extra increment that
describes the deviation from scaling,

$I({\mathbf{g}};\\{x,L\\})=I({\mathbf{g}};x/L)+\Delta I,$ (10)

and we will see that with samples from a sufficiently large number of embryos,
we can make a reliable estimate of $\Delta I$. The smaller the fraction
$\Delta I/I({\mathbf{g}};x/L)$ the closer the system is to a mathematical
ideal of scaling. More explicit expressions for $\Delta I$ are developed in
Appendix B and applied to experiments in §IV.

We emphasize that true scale invariance, corresponding to $\Delta I=0$, is a
very strong condition. Different levels of evidence for scaling in embryonic
development have inspired models in which competing mechanisms can provide
some cancellation of the intrinsic length scales determined by diffusion
constants and reaction rates [39, 40, 41, 42]. These models typically allow
for scaling in the position of a single discrete positional marker (e.g., the
middle of the embryo), or for approximate scaling across a larger segment of
the relevant axes. True scale invariance would require new dynamical
mechanisms.

##  III Stripes and boundaries

In the early fly embryo, information about position along the
anterior–posterior axis flows from maternal inputs through the network of gap
genes to the pair-rule genes [17]. The pair-rule genes are expressed in
interdigitating striped patterns that provide a preview of the segmented body
plan in the fully developed organism; these stripes are visible within three
hours after the egg is laid (Fig. 1A–C). The positions of pair-rule stripes
are a clear example of the positional markers discussed above.

Here we analyze the spatial profiles of gene expression for three of the pair-
rule genes—eve, prd, and run—measured using fluorescent antibody staining of
the corresponding proteins in more than one hundred embryos that were fixed
during nuclear cycle 14 (nc14), i.e. between 2 and $3\,$h of development [33].
Our results recapitulate earlier work [23] on a larger ensemble of embryos.

As soon as the stripes are visible it is straightforward to measure their
positions $x_{\rm i}$ [33]. The time during nc14 can be measured with $\sim
1\,{\rm min}$ precision by following the progression of the invaginating
cellularization membrane [32]. The stripe positions vary systematically in
time [43, 44, 45, 19, 46] and are well described by

$\frac{x_{\rm i}(t)}{L}=\frac{x_{\rm i}(t_{0})}{L}+s_{\rm i}(t-t_{0}),$ (11)

as shown for the Eve stripes in Fig 1D. Combining data from all time points,
we shift each embryo to the reference time $t_{0}=45\,{\rm min}$,

$\frac{x_{\rm i}(t)}{L}\rightarrow\frac{x_{\rm i}(t)}{L}-s_{\rm i}(t-t_{0}).$
(12)

We use this same procedure for the Prd and Run stripes, although these become
clear only at slightly later times.

Figure 1: Precise scaling of pair-rule stripes in the Drosophila embryo. (A)
Bright-field image overlaid with fluorescent antibody staining for Eve protein
(fuschia), focusing on the mid-sagittal plane with the dorsal side up;
scalebar is $100\,\mu{\rm m}$. (B) Expression of Eve in the second half of
nuclear cycle fourteen (nc14). Solid line is the mean, and shaded region is
the standard deviation across $N_{\rm em}=108$ embryos in a time window
between 30 and $60\,$min from the start of nc14. Inset shows a single nucleus
with a white square (width $0.01L$) used to average intensities. (C) Eve
expression profiles as a function of relative position along the body axis for
12 time bins during nc14, as indicated by color. (D) Linear dynamics of Eve
peak positions during nc14, fit to Eq. (11). (E) Absolute positions of Eve
peaks measured from the anterior pole referred to $t_{0}=45\,{\rm min}$, as in
Eq. (12), plotted vs. embryo length. (F) Standard deviation of scaled stripe
positions as a function of mean position for three pair-rule genes, and for
the cephalic furrow (CF, see Appendix C). Error bars are standard deviations
from bootstrapping. Black curves with red shading (bootstrapped errors) are
estimates of precision based on anchoring in Eqs. (5–7), and $d$ is the
spacing between neighboring cells.

Figure 1E shows that the stripe positions $x_{\rm i}$ measured from the
anterior pole are proportional to the length of the embryo $L$. More
precisely, if we fit these linear relations then intercepts are zero and
slopes are equal to the mean fractional positions, as in Eq. (3), both results
with error bars of less than $1\%$ (Appendix C). This provides _prima facie_
evidence for scaling of the pair-rule stripes, reinforcing the conclusions of
earlier work [18, 19, 20, 21].

We can go beyond the mean behaviors to look at fluctuations around these
means. For each stripe $\rm i$ in each embryo $\alpha$, we can write

${{x_{\rm i}^{\alpha}}\over{L^{\alpha}}}=\langle f_{\rm i}\rangle+\delta
f_{\rm i}^{\alpha},$ (13)

where $\langle\cdots\rangle$ now is an average over all the embryos in our
sample. The variance of the relative position is $\sigma_{f_{\rm
i}}^{2}=\langle(\delta f_{\rm i})^{2}\rangle$, and Fig. 1F shows that
$\sigma_{f_{\rm i}}\leq 0.01$ for all 21 pair rule stripes that we measure.
This is consistent with previous measurements, and with the information
content of the gap gene expression patterns that feed into the generation of
pair-rule stripes [47, 33], but earlier work did not address scaling
explicitly.

As a caution, we note that the observation of scaling in fixed embryos would
be trivial if variations in embryo length were dominated by shrinkage during
fixation. Across $N_{\rm em}=609$ fixed embryos used for the analysis of gap
genes (below) we find a mean length $\langle L\rangle_{\rm fix}=455\,\mu{\rm
m}$, while across $N_{\rm em}=610$ live embryos (§V) we find $\langle
L\rangle_{\rm live}=490\,\mu{\rm m}$. Hence, shrinkage with fixation is a bit
less than $10\%$ across many different experiments. But the variations in
length are almost the same, $(\sigma_{L}/\langle L\rangle)_{\rm fix}=0.038$,
while $(\sigma_{L}/\langle L\rangle)_{\rm live}=0.037$. The small extra
variance in the length of fixed embryos cannot explain the scaling behavior
that we observe.

Figure 1F also shows that the fluctuations in fractional position are smaller
than the bound on mechanisms that have no explicit scaling, from Eq. (7). This
bound is very tight, because of the small variance in emrbyo lengths, and thus
requires extreme precision in the measurement and biological reproducibility
of the fractional positions to demonstrate scaling. To emphasize the
importance of precision, we note that the position of the cephalic furrow is
directly regulated by pair rule gene expression [48], but it has a slightly
higher relative positional variance, due to the experimental difficulty of
defining morphological features to less than the width of a single cell [49].
Here we show explicitly that the furrow position scales with embryo length
(Appendix C). Even though the precision of the CF position is almost $\sim
1\%$ in the scaled coordinates [49], this alone is not sufficient to reject
the hypothesis that positions are defined in absolute rather than relative
coordinates, as can be seen from Fig. 1F.

The pair rule stripes are shaped by input from the gap genes [50], and it is
natural to ask whether the scaling behavior that we observe is inherited from
these inputs. The gap genes were long discussed in terms of “expression
domains,” as if they were on/off switches [51, 52, 53, 54]. We now know that
this misses a substantial fraction of the positional information encoded by
these genes [47, 55, 33], but defining the boundaries of the expression
domains as positional markers (Fig. 2A–D) allows us to give a preliminary
analysis of scaling by following the same ideas as for the positions of the
pair-rule stripes.

Previous experiments have measured the expression profiles of the gap genes
[33], staining $N_{\rm em}=609$ fixed embryos in nc14 with fluorescent
antibodies directed at the proteins encoded by the gap genes (Fig. 2A–D). We
define expression boundaries as the positions where the concentrations are
half their maximum mean value, and we correct their relative positions to
$t_{0}=45$ min as above. Figure 2E shows that all thirteen of the gap gene
boundaries defined in this way have absolute positions that scale precisely
with embryo length, as with the positions of the pair rule stripes. The
accuracy of this scaling again is better than $\sim 1\%$, and this precision
is better than the limiting performance of mechanisms that do not have some
explicit sensitivity to embryo length (Fig. 2F). For the gap genes, this
procedure allows us to span almost the entire range of the anterior–posterior
axis.

In summary, stripes and boundaries of gene expression in the early fly embryo
provide discrete positional markers, and the absolute positions of these
markers are in all cases proportional to the length of the embryo. This is
consistent with previous observations [18, 19, 20, 21], but the precision of
the scaling that we observe here is surprising. This suggests that the
underlying genetic network exhibits true scale invariance, which we now test
using the information decomposition [Eq. (10)].

Figure 2: Precise scaling of gap gene expression boundaries. Expression
profiles of (A) Hunchback (Hb), (B) Giant (Gt), (C) Knirps (Kni), and (D)
Krüppel (Kr), based on immunofluorescent staining (Appendix D). Means (solid
lines) and standard deviations (shading) across embryos aligned by scaled
position $x_{s}$. Vertical lines indicate the mean positions of expression
boundaries as well as a small peak in Kni. (E) Absolute position of all gap
gene boundaries as a function of the embryo length. Dashed black line
indicates the position of the posterior of the embryo. Boundary positions are
time-corrected to $t_{0}=45\,{\rm min}$, as with the stripe positions in Fig.
1D. (F) Standard deviation of scaled boundary positions as a function of mean
position for all 13 markers. Error bars are standard deviations from
bootstrapping. Black curves with red shading (bootstrapped errors) are
estimates of precision based on anchoring in Eqs. (5–7), and $d$ is the
spacing between neighboring cells. Horizontal dashed lines denote the distance
$d$ and half-distance $d/2$, between neighboring nuclei. Dotted gray line
indicates 1% precision.

##  IV Absolute vs. scaled positional information

The concentrations of morphogens provide cells with information about their
position in the embryo. This “positional information” [30] can be measured in
bits if we have access to data on the mean and variability of spatial profiles
for the concentration of the relevant molecules [47, 55]. The local expression
levels of individual gaps genes convey roughly two bits of information about
position, twice what is possible in a model of on/off expression domains.
Taken together all four gap genes provide $\sim 4.2\,{\rm bits}$, sufficient
to specify positions with $\sim 1\%$ accuracy along the anterior–posterior
axis, as seen above. However, these earlier analyses assumed, implicitly, that
information is about the fractional or scaled position. Is this correct?

Figure 3: Expression of Hb in scaled coordinates. (A) Mean concentration of
Hb, $\langle g_{\rm Hb}(x_{s})\rangle$, vs scaled position (solid line, as in
Fig. 2A) and the conditional distribution $P(g_{\rm Hb}|x_{s})$ around this
mean (shading). Intensity bin size is 0.05 maximum $\langle g_{\rm
Hb}\rangle$. (B) A slice through the conditional distribution at $x_{s}=\
0.47$ (dashed black lines) compared with distributions estimated from embryos
in narrow bins of length, $P_{L}(g_{\rm Hb}|x_{s})$. Lengths were binned in 5
bins with an equal number of embryos in each, such that each bin contains
about 60 embryos with variations in $L$ of less than $1\%$. Mean lengths in
each bin are indicated at the upper right of each panel. Probability
distributions of $g_{\rm Hb}$ are estimated using a kernel density estimator
with a Gaussian kernel that has width $\delta g=0.07\times\max_{x_{s}}\langle
g_{\rm Hb}(x_{s})\rangle$.

Figure 4: Near zero deviation from perfect scaling, in bits. (A) The extra
information $\Delta I$ that Hb expression levels carry about absolute rather
than scaled position, defined by Eq. (10) and evaluated from Eq. (14).
Estimates are based on random choices of $N$ embryos out of the full
experimental ensemble (points; circles show means with standard deviations),
and the extrapolation $N_{\rm em}\rightarrow\infty$ follows the methods of
Appendix E (red line). The result is $\Delta I=0.00\pm 0.008\,{\rm bits}$ (red
circle with error bar). (B) The extra information $\Delta I$ conveyed by all
four gap genes together, defined as in (A) by Eq. (10) but now evaluated using
Eq. (15). Symbols as in (A); the result is $\Delta I=0.038\pm 0.039\,{\rm
bits}$. Error bars are larger because we are analyzing a multidimensional
code, but there still is no significant difference from $\Delta I=0$.

The key to separating information about scaled vs. absolute position is to
compare the variance in morphogen concentrations at a scaled position $x_{s}$
depending on whether we constrain the length of the embryo (Appendix B).
Qualitatively, if there is perfect scaling then knowing the length would not
add any information with which to predict the morphogen concentration. Since
information is mutual this would mean that all the available information is
about the scaled position. To test this quantitatively in the context of the
gap genes, we have assembled data on $N_{\rm em}=301$ embryos, in each of
which we have reliable simultaneous measurements on the spatial profiles of
expression in all four gap genes, as described in Appendix D.

Figure 3A shows the spatial profile of Hb as a function of scaled position
along the anterior–posterior axis. At each scaled position $x_{s}=x/L$ we can
visualize the distribution of expression levels, which is well approximated as
a Gaussian (Appendix D and [47]). We can then ask if this distribution changes
when we look only at embryos in a narrow range of lengths $L$, and the answer
is no (qualitatively; Fig. 3B). Quantitatively we want to estimate the
difference in entropy between these two distributions, averaged over $x_{s}$
and $L$, which will give us the deviation from scaling $\Delta I$ in Eq. (10),
as explained in Appendix B. The calculation of the entropy simplifies in the
Gaussian approximation, depending just on the variances as in Eq. (48),

$\Delta I={1\over
2}\langle\log_{2}[\sigma_{g}^{2}(x_{s})]\rangle_{x_{s}}-{1\over
2}\langle\log_{2}[\sigma_{g}^{2}(x_{s}|L)]\rangle_{x_{s},L},$ (14)

where $\sigma_{g}^{2}(x_{s}|L)$ is the variance in concentration at scaled
position $x_{s}$ across embryos of length $L$ and $\sigma_{g}^{2}(x_{s})$ is
the same variance computed across all embryos.

Applying Eq. (14) requires estimating the relevant variances and also making
bins along the $x_{s}$ and $L$ axes. For the scaled position we choose bins of
size $\Delta x_{s}=0.01$, consistent with the precision that we see in Figs. 1
and 2. To sample the range of embryo lengths we use $N_{\rm
bins}=5,\,10,\,15,$ or $20$ adaptive bins, and find the same results in all
cases (Appendix E). As is well known, estimates of entropy or information are
subject to systematic errors [56, 38]. In the present case, if we substitute
estimates of the variances into Eq. (14), we find a nonzero result for $\Delta
I$. But suppose we include different numbers of embryos in our analysis. In
that case, we see that our estimate of $\Delta I$ depends on $1/N_{\rm em}$ as
expected theoretically [56, 38], and having seen this predicted dependence we
can extrapolate $N_{\rm em}\rightarrow\infty$. In particular, if we shuffle
the data so that the true $\Delta I=0$, then our estimation procedure returns
a random number with zero mean and standard deviation equal to our quoted
error bar, demonstrating that we have control over the systematic errors.
These now standard analysis methods are explained more fully in Appendix E.

Results of this analysis for Hb are shown in Fig. 4A. Using all $N_{\rm
em}=301$ embryos in our data set produces a very small estimate of $\Delta I$,
but even this is exaggerated by systematic errors as we see by changing
$N_{\rm em}$. Our best estimate extrapolates to zero as $N_{\rm
em}\rightarrow\infty$, with an error bar smaller than $0.01\,{\rm bits}$. When
we repeat the same analyses for each of the other gap genes (i.e., Gt, Kni,
and Kr), we get the same result (Appendix E).

Figure 5: The maternal input Bicoid does not scale. (A) Measurements of Bcd
concentration in $N_{\rm em}=582$ live embryos are grouped into eight classes
by embryo length $L$ and averaged. There is only one global normalization, so
this shows that absolute concentrations have the same dependence on absolute
position $x$ across all classes. (B) The same data plotted vs. scaled position
$x_{s}=x/L$. Profiles separate, providing evidence against scaling. (C) Extra
information $\Delta I$ that Bcd concentration provides about absolute vs.
scaled position, defined by Eq. (10) and evaluated from Eq. (14). Symbols as
in Fig. 3, but the extrapolation now leads to a significantly nonzero value of
$\Delta I=0.1\pm 0.02\,{\rm bits}$. Data from [49].

We can generalize this analysis to consider all four gap genes simultaneously.
Now the role of the variance in Eq. (14) is played by the covariance matrix of
the fluctuations, as in Eq. (52):

$\Delta I={1\over
2}\langle\log_{2}\left[||\Sigma(x_{s})||\right]\rangle_{x_{s}}-{1\over
2}\langle\log_{2}\left[||\Sigma(x_{s}|L)||\right]\rangle_{x_{s},L}.$ (15)

Here $||\Sigma(x_{s}|L)||$ is the determinant of the covariance matrix
describing fluctuations in the expression levels of all four genes at scaled
position $x_{s}$ across embryos of length $L$, and $\Sigma(x_{s}|L)$ is the
covariance computed across all embryos. Because we are looking at higher
dimensional variations the impact of the finiteness of our data set is larger,
but again we see the predicted dependence on $1/N_{\rm em}$ and can
extrapolate to give $\Delta I=0.038\pm 0.039\,{\rm bits}$ (Fig. 4B). Once
again this is consistent with $\Delta I=0$: there is no significant evidence
for encoding of information about absolute, as opposed to scaled position.

Although the number of bits has meaning, it is useful to express the deviation
from perfect scaling as a fraction of the information available about scaled
position [47, 55],

${{I({\mathbf{g}}\rightarrow\\{x,L\\})-I({\mathbf{g}}\rightarrow
x/L)}\over{I({\mathbf{g}}\rightarrow x/L)}}=0.009\pm 0.009.$ (16)

Thus the patterns of gap gene expression scale with $1\%$ accuracy, not just
at discrete positional markers but across the entire range of graded spatial
variations.

##  V Maternal inputs do not scale

Having observed scaling in the pair rule stripe positions and followed this
back to the gap genes, it is natural to ask if we can go one step further and
trace the scaling behavior of the patterning system to the maternal inputs. Of
the three maternal inputs that drive patterning along the anterior–posterior
axis of the fly embryo, much attention has been given to Bicoid (Bcd) [24,
25]. The protein is present at high concentrations in the anterior, and there
is a nearly exponential decay of concentration with distance toward the
posterior; one can monitor the dynamics of Bicoid protein concentrations
quantitatively in live embryos using GFP-fusions [34].

Comparison across closely related species of flies shows that the length scale
of this exponential decay varies in proportion to the mean length of the
embryo [57]. Insertion of bicoid genes from other species into Drosophila
melanogaster produces protein concentration profiles with length scales
appropriate to the host, but these are not sufficient to rescue the embryo
from deletion of the native Bcd [58]. These results emphasize the subtlety of
comparison across species and the impact of genetic variations, leading us to
re-examine the behavior of Bcd profiles across a large number of live embryos
drawn from the same inbred laboratory strain used in the analysis of gap and
pair rule genes.

Figure 5 analyzes Bcd profiles from $N_{\rm em}=582$ live embryos [49].
Measurements are taken during a small temporal window in nuclear cycle
fourteen [34], and the only normalization (as with the gap genes) is to
subtract a common background level from all the embryos and set the highest
mean concentration to one. When we group the embryos into eight classes based
on their length $L$, we see that the average concentration profiles in all
groups are the same when plotted vs. absolute position, except for small
effects at the posterior pole (Fig. 5A). If we plot vs. scaled position the
different groups of embryos separate significantly (Fig. 5B), providing direct
evidence against scaling. We make this precise using the same information
theoretic approach as above and now find a significant nonzero value of
$\Delta I=0.1\pm 0.02\,{\rm bits}$ (Fig. 5C). This may seem like a small
number, but this is related to the $\sim 4\%$ scale of variations in embryo
length. We conclude that the maternal inputs do not scale, in agreement with
earlier suggestions [18].

We emphasize that the absence of scaling in the maternal inputs should not be
interpreted as a form of noise. Indeed, absolute concentrations of Bcd protein
are highly reproducible across embryos and this can be traced to highly
reproducible numbers of mRNA molecules [59, 49, 60]. Instead, we should think
of the maternal inputs as a nearly deterministic response to the boundary
conditions in the embryo, which also have a direct impact on the gap genes;
see Eqs. (19, 20) below.

##  VI Scaling and zero modes

The results here strongly support the view that patterns of gap gene
expression are genuinely scale invariant and that this is an emergent property
of the gap gene network. Here we take the precise mathematical notion of scale
invariance literally and explore its implications. While we do not pretend to
have a detailed model, it is useful to have in mind a class of models for how
patterns might form. As a caution we recall Turing’s introductory remarks
[14]: “This model will be a simplification and an idealization, and
consequently a falsification.”

If we focus on variations just along the anterior–posterior axis $x$, and
ignore the discreteness of nuclei, then the concentration $g_{\rm i}$ of
protein encoded by gene $\rm i$ plausibly obeys an equation of the form

${{\partial g_{\rm i}}\over{\partial t}}=D_{\rm i}{{\partial^{2}g_{\rm
i}}\over{\partial x^{2}}}+R_{\rm i}F_{\rm i}({\mathbf{g}})-{1\over{\tau_{\rm
i}}}g_{\rm i}.$ (17)

Here $D_{\rm i}$ is the diffusion constant of species $\rm i$, $R_{\rm i}$ is
the maximum rate at which these proteins can be synthesized, $\tau_{\rm i}$ is
their lifetime before being degraded, and $F_{\rm i}({\mathbf{g}})$ describes
all the potentially complex interactions by which all the proteins regulate
the expression of gene $\rm i$. We assume that the mRNA and protein dynamics
have separate time scales so that one is effectively slaved to the other and
we can write only one variable for each gene. Further, we neglect time scales
that might arise in the process of regulation itself, such as switching
between different regulatory states, so that $F_{\rm i}({\mathbf{g}})$ is an
instantaneous function of the relevant concentrations. These assumptions are
quite conventional, and other than this what we have written is very general.
For example, the function $F_{\rm i}({\mathbf{g}})$ could describe both
activating and repressive interactions, and these interactions could be
combinatorial. These equations include as special cases Turing’s original
models [14] and their intellectual descendants [61, 62].

The maximum steady state concentration of each protein is $R_{\rm i}\tau_{\rm
i}$, and we can choose units in which this is equal to one, as with the
normalized profiles of gap gene expression in Fig. 2A–D. For simplicity we
will assume that all the decay times are the same, $\tau_{\rm i}=\tau$,
although this is not essential for what follows; finally, we choose units of
time such that $\tau=1$. Then we have

${{\partial g_{\rm i}}\over{\partial t}}=\lambda_{\rm
i}^{2}{{\partial^{2}g_{\rm i}}\over{\partial x^{2}}}+F_{\rm
i}({\mathbf{g}})-g_{\rm i},$ (18)

where the length scale $\lambda_{\rm i}=\sqrt{D_{\rm i}\tau}$. This describes
an autonomous network, which is not quite realistic for the gap genes—which
are driven by maternal inputs—but should be sufficient to draw qualitative
conclusions about the implications of scale invariance.

The length of the embryo appears not in the dynamical equations but in the
boundary conditions. For most proteins, there can be no diffusive flux of
molecules into or out of the ends of the embryo, so that

$-D_{\rm i}{{\partial g_{\rm i}}\over{\partial x}}{\bigg{|}}_{x=0}=D_{\rm
i}{{\partial g_{\rm i}}\over{\partial x}}{\bigg{|}}_{x=L}=0.$ (19)

The situation for maternal inputs is different; as an example, in making the
egg the mother places mRNA for the protein Bicoid (Bcd) at the anterior end
($x=0$), and this is translated continuously, so that

$-D_{\rm Bcd}{{\partial g_{\rm Bcd}}\over{\partial x}}{\bigg{|}}_{x=0}=T_{\rm
Bcd},$ (20)

where $T_{\rm Bcd}$ is the rate of translation in appropriate units.

Let us imagine that the final pattern we observe is in steady state, so that

$0=\lambda_{\rm i}^{2}{{\partial^{2}g_{\rm i}(x;L)}\over{\partial
x^{2}}}+F_{\rm i}({\mathbf{g}})-g_{\rm i}(x;L),$ (21)

where the notation reminds us that length dependence can arise once we impose
the boundary conditions. If we have true scale invariance as in Eq. (1) then
if we make a small change in the length of the embryo, so that $L\rightarrow
L+\delta L$, the expression levels should change as

$\displaystyle g_{\rm i}(x;L)$ $\displaystyle\rightarrow$ $\displaystyle
g_{\rm i}(x;L)+{{\delta L}\over{L}}\psi_{\rm i}(x/L)$ (22)
$\displaystyle\psi_{\rm i}(x_{s})$ $\displaystyle=$ $\displaystyle-
x_{s}\Phi_{\rm i}^{\prime}(x_{s}),$ (23)

but Eq. (21) still must be true. This requires that

$\sum_{\rm j}\left[\left(\lambda_{\rm i}^{2}{{\partial^{2}}\over{\partial
x^{2}}}-1\right)\delta_{\rm ij}+{{\partial F_{\rm i}}\over{\partial g_{\rm
j}}}{\bigg{|}}_{{\mathbf{g}}={\mathbf{\Phi}}}\right]\psi_{\rm j}(x/L)=0.$ (24)

This seemingly abstract condition has a direct implication for the dynamics of
the network.

Suppose that the system is close to its steady state so that we can write

$g_{\rm i}(x;L;t)=\Phi_{\rm i}(x/L)+\delta g_{\rm i}(x;t)$ (25)

and $\delta{\mathbf{g}}$ is small. Then we can linearize the dynamics in Eq.
(18) to yield

${{\partial(\delta g_{\rm i})}\over{\partial t}}=\sum_{\rm
j}\left[\left(\lambda_{\rm i}^{2}{{\partial^{2}}\over{\partial
x^{2}}}-1\right)\delta_{\rm ij}+{{\partial F_{\rm i}}\over{\partial g_{\rm
j}}}{\bigg{|}}_{{\mathbf{g}}={\mathbf{\Phi}}}\right]\delta g_{\rm j}.$ (26)

We recognize the term in brackets as the same one that appears in Eq. (24). To
understand this connection it is useful to think of all possible spatial
patterns of gene expression as points in a high dimensional space.

Concretely we can write

$\delta g_{\rm i}(x;t)=\sum_{\mu}a_{\mu}(t)\phi_{\rm i}^{\mu}(x)$ (27)

where the functions $\\{\phi_{\rm i}^{\mu}(x)\\}$ are the spatial “modes” of
expression and the set $\\{a_{\mu}\\}$ provides the coordinates of one
expression profile in this multidimensional space. The number of modes is the
number of genes times the number of independent points along the $x$ axis,
e.g. the number of rows of cells; for the gap genes the result is that the
space has a dimensionality $d>300$. We can choose these modes as
eigenfunctions of the operator that appears in both Eqs. (24) and (26),

$\sum_{\rm j}\left[\left(\lambda_{\rm i}^{2}{{\partial^{2}}\over{\partial
x^{2}}}-1\right)\delta_{\rm ij}+{{\partial F_{\rm i}}\over{\partial g_{\rm
j}}}{\bigg{|}}_{{\mathbf{g}}={\mathbf{\Phi}}}\right]\phi_{\rm
j}^{\mu}(x)=-\lambda_{\mu}\phi_{\rm i}^{\mu}(x),$ (28)

where $\lambda_{\mu}\geq 0$ means that the steady state is stable. Then so
long as the deviations from the steady state are small, the dynamics of the
network are simple in this coordinate system,

${{da_{\mu}(t)}\over{dt}}=-\lambda_{\mu}a_{\mu}(t).$ (29)

Through Eq. (24) we see that perfect scale invariance implies a “zero mode,” a
particular mode of gene expression associated with eigenvalue
$\lambda_{\mu}=0$. Importantly this is not the steady state pattern itself,
but an additional mode.

The existence of a zero mode has several implications:

  * •

Most literally, one component in the spatial pattern of gene expression will
relax very slowly to its steady state, much more slowly than other components.
Formally the relaxation should be as a power of time rather than exponential.

  * •

The dynamics describe a “restoring force” that pulls the patterns of gene
expression toward their steady state values; the eigenvalues are the spring
constants associated with these restoring forces. Along the zero mode, there
is no (linear) restoring force, and in the presence of any finite noise, the
fluctuations along this mode will be very large compared with other modes.

  * •

Along directions with nonzero $\lambda_{\mu}$ the fluctuations in $\mathbf{g}$
will be approximately Gaussian so long as they remain small, as we see for the
gap genes. But along the zero mode, there should be some deviation from
Gaussian behavior.

There is evidence that the spatial patterns of gap gene expression can be
compressed into a lower dimensional space, consistent with the idea that a
zero mode dominates the dynamics [63]. The ($4\times 4$) covariance matrix of
fluctuations in gap gene expression is dominated by a single mode at almost
all locations along the anterior–posterior axis, this large variance mode
relaxes $\sim 10\times$ more slowly than the lower variance modes, and one can
even see hints of non–Gaussian behavior [26].

The existence of a zero mode is a statement about the linearized dynamics. If
the absence of a linear restoring force continues for finite deviations from
the steady state then there is a line of attracting spatial patterns rather
than a single stable pattern. Different points along this line are the
patterns appropriate to embryos of different lengths, and the final pattern is
selected by boundary conditions. Line attractors have long been discussed for
neural networks [64]. It has been noted that models of the gap gene network
might support such line attractors [28], and there are also suggestions that
internal dynamics of the network can generate approximate scaling [29]. The
observation of nearly perfect scale invariance in the real network leads us to
a much sharper version of these ideas.

##  VII Discussion

Scale invariance is an appealing concept. It quantifies the intuition that
organisms are built from parts that are in proportion to one another,
independent of an individual organism’s overall size. There is a long history
of searching for such scaling not just in adult organisms but at early stages
of development, and the fruit fly Drosophila melanogaster has been a
particular target for these studies [40, 20, 21, 19, 29]. If we compare
related species of flies we can see spatial patterns of gene expression that
scale, on average, across $10\times$ changes in embryo length [57, 58], and
similar results are obtained within a single species but with artificial
selection for length variation [22]. It has always been less clear whether
scaling occurs without such large genetic variations, across the natural
length variations in a single species.

We have explored scaling across many embryos from a quasi-inbred laboratory
stock, minimizing genetic variation. Across this ensemble, we see length
fluctuations with a standard deviation of $\pm 4\%$ but embryos in the tails
of the distribution have lengths $\pm 10\%$ from the mean (Fig. A1). Following
previous work, we measured the positions of discrete markers—such as the CF
position, the peaks of pair-rule stripes, and the boundaries of gap gene
domains—and found precise scaling of the absolute positions with embryo
length. This is consistent with previous results, but what is new is the
precision that we observe: markers are at positions that are scaled relative
to the embryo length with an accuracy of $\sim 1\%$ across the full extent of
the anterior–posterior axis. This observed precision excludes a broad class of
models that combine information from both ends of the embryo without explicit
scaling [39, 40, 41, 42].

There remains a gap between the positioning of discrete markers and the fuller
notion of scale invariance. The gap gets smaller as we track more markers
across a wider range of positions, but it would be attractive to address scale
invariance directly. We have introduced an information theoretic approach that
analyzes the full, graded spatial profiles of gene expression and measures
similarity in the natural units provided by the intrinsic noise levels of
these profiles. Concretely, we introduce a decomposition of the information
that morphogen concentrations provide about position into a component about
scaled position and a deviation from scaling. Applied to the gap genes in the
early fly embryo, the result is clear: the deviation from scaling is less than
one percent of the total positional information. It is perhaps surprising that
we can make such a precise statement about the functional output of a complex
network.

In contrast to the results for the gap genes and the pair-rule genes, at least
one of the maternal inputs, Bicoid, does not exhibit scaling. We can see this
“by eye,” simply plotting profiles vs. absolute or scaled position, and these
impressions are quantified by the same information theoretic approaches used
to demonstrate scaling in the gap genes. Error bars again are in the range of
$\sim 0.01\,{\rm bits}$, but the deviation from scaling now is $\sim 10\times$
as large. The conclusion is that near-perfect scale invariance is an emergent
property of the gap gene network.

If we take scale invariance as a precise mathematical statement then the
dynamics of the underlying genetic network must have a zero mode. This is
equivalent to saying that the dynamics do not have a single attractor, but
rather a line of attractors as in models for short-term memory in neural
networks [64]. Then position along this line is chosen by the boundary
conditions and hence the length of the embryo. A zero mode would provide
connections among several otherwise disparate observations on the gap genes.

Finally, recent experiments on mammalian pseudo-embryos suggest that scale
invariance may be a more universal feature of genetic networks underlying
developmental pattern formation [65]. In these self-organizing cell aggregates
derived from stem cells, scale invariance emerges without fixed boundary
conditions, but rather with boundaries that move as the aggregate grows. The
existence of a zero mode in the regulatory network becomes even more
attractive as a general mechanism for scaling.

###### Acknowledgements.

We are grateful to E. F. Wieschaus for his advice and for many inspiring
discussions. We thank M. Biggin and N. Patel for sharing the antibodies used
in these measurements. This work was supported in part by US National Science
Foundation Grant PHY–1734030 (Center for the Physics of Biological Function);
by National Institutes of Health Grants R01GM077599 and R01GM097275; by the
Simons Foundation; by the John Simon Guggenheim Memorial Foundation.

##  Appendix A Natural length variations of embryos in a laboratory strain of
flies

As described in the main text, much previous work on scaling has exploited the
natural variation in embryo lengths across the evolutionary tree or the
variations that can be selected artificially over reasonable times. Here we
use variations in length that occur within a single laboratory strain, OreR,
minimizing genetic variations. Measurements on large numbers of live embryos
are taken from Refs. [49, 35] and on fixed embryos from Ref. [33].

As an example, Fig. A1 shows the probability distribution of embryo lengths
$L$ estimated from $N_{\rm em}=610$ living dechorionated embryos (Bcd-GFP 2XA
strain in [49]). The mean length of the embryos is $\langle L\rangle=490\pm
0.76\,\mu{\rm m}$, and the standard deviation is $\sigma_{L}=18\pm
1.06\,\mu{\rm m}$. This corresponds to a fractional variation
$\sigma_{L}/\langle L\rangle=0.037$, and as noted in the main text our sample
is sufficiently large that it includes embryos $\pm 10\%$ from the mean. This
is true also in the case of fixed OreR embryos where we find
$\sigma_{L}/\langle L\rangle=0.038$ and $\sigma_{L}/\langle L\rangle=0.039$ in
the experimental ensembles used for the analysis fo the gap and pair rule
genes, respectively.

Figure A1: Distribution of live embryo lengths. Data from $N_{\rm em}=610$
embryos [49] analyzed in bins of size $\Delta L/\langle L\rangle=0.02$. The
mean embryo length $\langle L\rangle=490\pm 0.76\,\mu{\rm m}$. Overlaid is the
error bar indicating the standard deviation of $\sigma_{L}=0.0371\langle
L\rangle$, and each dot indicates the length of one of the embryos in our
sample.

##  Appendix B Decomposing information

We want to give explicit expressions that allow us to decompose positional
information, as in Eq. (10), based on estimates from real data. We give more
detail than usual in hopes of making the analysis accessible to a broader
audience.

Concentrations can depend both on the absolute position $x$ and the length of
the embryo $L$. We can rewrite this as a dependence on the scaled position
$x_{s}=x/L$ and the length $L$. Thus we have

$P\left({\mathbf{g}}|\\{x,L\\}\right)=P\left({\mathbf{g}}|\\{x/L,L\\}\right)=P_{L}(\\{g_{\rm
i}\\}|x_{s}),$ (30)

where $P_{L}$ is a distribution formed only across embryos of length $L$.
Similarly, we expect that cells are uniformly distributed along the $x$ axis
up to the length $L$, so that

$P(x,L)={1\over L}\Theta(1-x_{s})P(L),$ (31)

where $\Theta$ is the unit step function. Then we can substitute into Eq. (8):

$\displaystyle I(\\{g_{\rm i}\\}\rightarrow\\{x,L\\})$ $\displaystyle=$
$\displaystyle\int d{\mathbf{g}}\int dx\int
dL\,P\left({\mathbf{g}}|\\{x;L\\}\right)P(x,L)\log_{2}\left[{{P\left({\mathbf{g}}|\\{x,L\\}\right)}\over{P\left({\mathbf{g}}\right)}}\right]$
(32) $\displaystyle=$ $\displaystyle\int
dL\,P(L)\int_{0}^{1}dx_{s}\,P_{L}({\mathbf{g}}|x_{s})\log_{2}\left[{{P_{L}({\mathbf{g}}|x_{s})}\over{P\left({\mathbf{g}}\right)}}\right].$
(33)

Now we insert a factor of unity:

$\displaystyle\log_{2}\left[{{P_{L}({\mathbf{g}}|x_{s})}\over{P\left({\mathbf{g}}\right)}}\right]$
$\displaystyle=$
$\displaystyle\log_{2}\left[{{P_{L}({\mathbf{g}}|x_{s})}\over{P\left({\mathbf{g}}\right)}}{{P({\mathbf{g}}|x_{s})}\over{P({\mathbf{g}}|x_{s})}}\right]$
(34) $\displaystyle=$
$\displaystyle\log_{2}\left[{{P({\mathbf{g}}|x_{s})}\over{P\left({\mathbf{g}}\right)}}\right]-\log_{2}\left[P({\mathbf{g}}|x_{s})\right]+\log_{2}\left[P_{L}({\mathbf{g}}|x_{s})\right].$
(35)

Substituting, we can write

$I({\mathbf{g}}\rightarrow\\{x,L\\})=I_{1}+I_{2}+I_{3},$ (36)

where the three components are

$\displaystyle I_{1}$ $\displaystyle=$ $\displaystyle\int_{0}^{1}dx_{s}\int
d{\mathbf{g}}\,P({\mathbf{g}}|x_{s})\log_{2}\left[{{P({\mathbf{g}}|x_{s})}\over{P\left({\mathbf{g}}\right)}}\right]$
(37) $\displaystyle I_{2}$ $\displaystyle=$
$\displaystyle-\int_{0}^{1}dx_{s}\int
d{\mathbf{g}}\,P({\mathbf{g}}|x_{s})\log_{2}[P({\mathbf{g}}|x_{s})]$ (38)
$\displaystyle I_{3}$ $\displaystyle=$ $\displaystyle\int
dL\,P(L)\int_{0}^{1}dx_{s}\int
d{\mathbf{g}}\,P_{L}({\mathbf{g}}|x_{s})\log_{2}[P_{L}({\mathbf{g}}|x_{s})],$
(39)

We can identify the three terms: First, $I_{1}$ is the information that the
concentrations of morphogens provide about the scaled position,

$I_{1}=I({\mathbf{g}}\rightarrow x/L).$ (40)

Second, $I_{2}$ is the entropy of the distribution of concentrations at a
particular value of scaled position $x_{s}$, averaged over this position,

$I_{2}=\langle S[P({\mathbf{g}}|x_{s})]\rangle_{x_{s}},$ (41)

where $S[Q]$ denotes the entropy of the distribution $Q$. Finally, $I_{3}$ is
the negative of the entropy of the same distribution but restricted to embryos
of length $L$, and then averaged also over $L$,

$I_{3}=-\langle S[P_{L}({\mathbf{g}}|x_{s})]\rangle_{x_{s},L}.$ (42)

Comparing with Eq. (10) we see that the deviation from scaling can be written
as the difference between two entropies, suitably averaged:

$\Delta I=\langle S[P({\mathbf{g}}|x_{s})]\rangle_{x_{s}}-\langle
S[P_{L}({\mathbf{g}}|x_{s})]\rangle_{x_{s},L}.$ (43)

This has a very simple interpretation: There is a deviation from scaling if
specifying the length of the embryo reduces the entropy of fluctuations in
morphogen concentration at a given scaled position.

These general expressions simplify enormously in the case where we have only a
single morphogen and the conditional distributions are Gaussian. In this case

$\displaystyle P(g|x_{s})$ $\displaystyle=$
$\displaystyle{1\over{Z(x_{s})}}\exp\left[-{1\over 2}\chi^{2}(g;x_{s})\right]$
(44) $\displaystyle\chi^{2}(g;x_{s})$ $\displaystyle=$
$\displaystyle{{[g-\langle g(x_{s})\rangle]^{2}}\over{\sigma_{g}^{2}(x_{s})}}$
(45) $\displaystyle Z(x_{s})$ $\displaystyle=$
$\displaystyle\sqrt{2\pi\sigma_{g}^{2}(x_{s})},$ (46)

where $\langle g(x_{s})\rangle$ is the mean and $\sigma_{g}^{2}(x_{s})$ is the
variance of $g$ at scaled positions $x_{s}$. Importantly the entropy of a
Gaussian distribution does not depend on the mean, and we have [37, 38]

$S\left[P(g|x_{s})\right]={1\over 2}\log_{2}\left[2\pi
e\sigma_{g}^{2}(x_{s})\right].$ (47)

Thus we find the deviation from scaling is

$\Delta I={1\over
2}\langle\log_{2}[\sigma_{g}^{2}(x_{s})]\rangle_{x_{s}}-{1\over
2}\langle\log_{2}[\sigma_{g}^{2}(x_{s}|L)]\rangle_{x_{s},L},$ (48)

where $\sigma_{g}^{2}(x_{s}|L)$ is the variance in concentration at scaled
position $x_{s}$ across embryos of length $L$. In other words, there is a
deviation from scaling if the variance in morphogen concentration is reduced
by knowing the length of the embryo.

This result can be generalized to multiple morphogens if we stay in the
Gaussian approximation. Now with $d$ genes, we have

$\displaystyle P({\mathbf{g}}|x_{s})$ $\displaystyle=$
$\displaystyle{1\over{Z(x_{s})}}\exp\left[-{1\over 2}\chi^{2}(\\{g_{\rm
i}\\};x_{s})\right]$ (49) $\displaystyle\chi^{2}(\\{g_{\rm i}\\};x_{s})$
$\displaystyle=$ $\displaystyle\sum_{{\rm i}=1}^{d}\sum_{{\rm j}=1}^{d}[g_{\rm
i}-\langle g_{\rm i}(x_{s})\rangle][\Sigma^{-1}(x_{s})]_{\rm ij}[g_{\rm
j}-\langle g_{\rm j}(x_{s})\rangle]$ (50) $\displaystyle Z(x_{s})$
$\displaystyle=$
$\displaystyle\left[(2\pi)^{d}||\Sigma(x_{s})||\right]^{1/2},$ (51)

where $\Sigma(x_{s})$ is the covariance matrix of fluctuations in
concentration at scaled position $x_{s}$ and $||\Sigma(x_{s})||$ is the
determinant of this matrix. Following the same logic as in the case of one
gene we have

$\Delta I={1\over
2}\langle\log_{2}\left[||\Sigma(x_{s})||\right]\rangle_{x_{s}}-{1\over
2}\langle\log_{2}\left[||\Sigma(x_{s}|L)||\right]\rangle_{x_{s},L}.$ (52)

Even if $P_{L}({\mathbf{g}}|x_{s})$ is perfectly Gaussian, averaging over $L$
could generate non-Gaussian behavior in $P({\mathbf{g}}|x_{s})$. We are
neglecting this here, but since we find that $\Delta I$ is very small, and for
the gap genes consistent with $\Delta I=0$, both the conditional and averaged
distributions are very nearly Gaussian, as seen in Fig. 4B.

##  Appendix C Cephalic furrow and scale invariance

Upon the onset of gastrulation (i.e., three hours after fertilization), the
cephalic furrow (CF) emerges as the first macroscopic morphological feature
along the anterior–posterior axis in the developing fly embryo. It results
from collective cell movement that can be seen using bright-field microscopy.
There are hints in early experiments that this marker is positioned very
precisely [25]. Modern experiments show that, as a fraction of the embryo
length $L$, CF position $x_{\rm CF}$ is reproducible to nearly $1\%$ accuracy
[49].

When we plot CF position in absolute units as a function of embryo length we
observe a linear relationship with zero intercept (Fig. A2A). The slope of
this fit is well within 1% of the mean scaled position $\langle f_{\rm
CF}\rangle$. More generally, all the discrete positional markers that we track
(CF, pair-rule stripes, gap boundaries) have absolute positions that vary
linearly with embryo length; the intercepts of the best-fit linear relations
are zero; the slopes match the mean scaled positions of the markers (Fig. A2B)
as predicted by scale invariance [Eq. (3)]; and the precision of this match is
better than $1\%$.

Figure A2: Cephalic furrow and the proportionality of scaling. (A) The
absolute position of the cephalic furrow measured in live emrbyos [49] is
proportional to the embryo length. Red line is the best fit with 95%
confidence intervals shown as dashed lines. The entire fit is shown in the
inset, emphasizing that the intercept is zero. (B) Slopes of absolute position
vs. embryo length for multiple positional markers, each plotted vs. its mean
scaled position. (C) Replotting of data in (B) to show that slopes and scaled
positions are equal within $1\%$, as predicted for perfect scaling [Eq. (3)].

##  Appendix D Aspects of the gene expression data

We analyze gene expression patterns for the pair-rule genes, the gap genes,
and the maternal input Bicoid. In each case, the concentration of the protein
is inferred from the intensity of a fluorescence signal. In each case images
are collected by focusing on the midsagittal plane, the extent of the embryo
is defined by thresholding the fluorescence intensity, and to avoid geometric
distortions we avoid the $5\%$ of the embryo at both the anterior and
posterior poles. Fluorescence intensities are averaged over a small area, as
shown in the inset to Fig. 1B, sliding along the dorsal rim of the embryo. In
live embryos, we can keep track of time during nc14 directly, while in fixed
embryos we use the progress of the cellularization as a clock with precision
$\sim 1\,{\rm min}$.

For each gene $\rm i$ we measure an intensity $I_{\rm i}^{\alpha}(x_{s})$ as a
function of scaled position in embryo $\alpha$. In each experiment, we
normalize by assuming that the minimum mean concentration is zero and we
choose units such that the maximum mean concentration is one. This defines

$g_{\rm i}^{\alpha}(x_{s})={1\over{S_{\rm i}}}\left[I_{\rm
i}^{\alpha}(x_{s})-B_{\rm i}\right],$ (53)

where the background is

$B_{\rm i}=\min_{x_{s}}\langle I_{\rm i}^{\alpha}(x_{s})\rangle$ (54)

and the scale factor is

$S_{\rm i}=\max_{x_{s}}\langle I_{\rm
i}^{\alpha}(x_{s})\rangle-\min_{x_{s}}\langle I_{\rm
i}^{\alpha}(x_{s})\rangle.$ (55)

Importantly there is no freedom to normalize the profiles measured in
individual embryos, which would distort our estimates of noise and variability
[59].

Data on three pair-rule genes—Eve, Prd, and Run—are taken from Ref. [33]. We
fit the sum of seven Gaussians to each profile, identifying stripe positions
with the centers of the Gaussians. We have also used the average peak profiles
as templates [23] and made more restricted fits to small segments of the peaks
[27]; results are the same with all three methods. Corrections for the drift
of peak position vs. time in nc14 were made as described in the main text.

Figure A3: Apparent variance of gap gene expression across multiple
experimental sessions. Results from eight sessions are largely reproducible
for each of the four genes. In regions where mean expression levels are near
zero, variances typically are $\sigma_{g}^{2}\ll 10^{-3}$, except in a handful
of sessions with highly variable backgrounds; these are excluded from further
analysis.

Simultaneous measurements on all four gap genes also were drawn from
experiments described in Ref. [33]. Because the analyses done here are so
demanding of data, we tried to merge data from as many independent
experimental sessions as possible. Most quantities are extremely reproducible
from session to session, but in a handful of sessions, we found anomalously
large variations in background fluorescence across the individual embryos.
Concretely, if we measure the variance of expression levels for the individual
genes, we typically find that $\sigma_{g}^{2}(x_{s})\ll 10^{-3}$ in regions of
$x_{s}$ with near zero mean (Fig. A3). In a few sessions, these fluctuations
in background are much larger, and these sessions are excluded; more
precisely, since all genes are measured simultaneously, excess background
variance in one gene is sufficient to exclude those data. This leaves five
independent sessions with a total of $N_{\rm em}=301$ embryos which we pool
into one data set for all further analyses.

For the analysis of gap gene expression boundaries, we mark the points that
are half-maximal along the sharp slopes, as indicated in Fig. 2. For the weak
peak of Kni expression near $x_{s}=0.33$ we fitted a Gaussian profile and took
the marker as the center of the Gaussian.

Gap gene profiles vary slowly but significantly throughout nc14. If we don’t
treat this variation explicitly it can be confused for noise, resulting in a
substantial overestimate of the variances and entropies. To separate the
temporal dynamics of the gap genes from their noise level, we follow Ref. [32]
and detrend the variations at each position, generalizing the treatment of the
stripe positions in Fig. 1D. The alternative is to focus only on a small
window of time [33], but this limits the size of the data set we can use.

Another systematic source of variation is the dependence of gap gene profiles
on the dorsal-ventral coordinate [32]. Previous work thus has been very strict
in analyzing embryos with narrowly defined orientations. To expand our data
set we are less strict, but this is problematic for the Kni profiles in the
range $0.15<x_{s}<0.45$, which contains a small peak. When analyzing Kni
alone, or all four gap genes together, we exclude this region. The alternative
is to analyze the other three genes together across the full length of the
anterior–posterior axis; results for $\Delta I$ are the same.

Figure A4: Fluctuations in gap gene expression are approximately Gaussian.
Distributions of z-scored fluctuations, as in Eq. (56), are estimated for each
individual gap gene, pooled across positions and embryos; error bars are
standard deviations. Black curves are Gaussians with zero mean and unit
variance.

An important assumption for our analysis is that the distribution of gene
expression at a given anterior–posterior position is Gaussian, as shown
previously [47]. For completeness, we reproduce this result for our larger
data set. In a single embryo $\alpha$ we observe a gene expression level
$g^{\alpha}(x_{s})$ at scaled position $x_{s}$. We compute the mean and
standard deviation across all the embryos in our ensemble and define
normalized deviations or z-scores

$\Delta^{\alpha}(x_{s})={{g^{\alpha}(x_{s})-\langle
g(x_{s})\rangle}\over{\sigma_{g}(x_{s})}}.$ (56)

We pool across all $\alpha=1,\,2,\,\cdots,\,N_{\rm em}$ embryos and across all
positions $x_{s}$ to estimate the probability density $P(\Delta)$. Results are
in Fig. A4 for each of the four gap genes.

Finally, measurements of Bicoid concentration are taken from fluorescent
imaging of live embryos expressions a Bicoid-GFP fusion [49]; we consider only
strain 2XA, which has the Bcd dosage closest to that of wild-type flies. With
live measurements we can choose a time window, in this case $t=16\pm
2\,\rm{min}$ after the start of nc14, avoiding any temporal detrending while
still including $N_{\rm em}=582$ embryos. Some measurements with missing data
points along the length of the embryo were excluded from this set.

##  Appendix E Estimating $\Delta I$ from limited data

Entropy and information depend on probability distributions, not just on their
moments, and thus are especially difficult to estimate from a finite sample of
data. Further, the entropy is a nonlinear function of the probabilities, and
so random errors in probability become systematic errors in our estimates of
information. This problem was appreciated in the very first attempts to use
information theory in the analysis of experiments on biological systems [56].
In the subsequent decades, many approaches have been developed, driven
especially by the analysis of neural coding. The approach we take here follows
the discussion in Appendix A.8 of Ref. [38]. Rather than just saying that we
follow established methods, we repeat some details that can be found in other
contexts in hopes that our presentation thus will be more accessible.

If we estimate an information-theoretic quantity such as $\Delta I$ in Eq.
(43) based on data from measurements in $N_{\rm em}$ independent embryos, then
with any simple estimation procedure our estimates will be biased:

$\Delta I=\Delta I_{\infty}+{{A(N_{\rm bins})}\over{N_{\rm em}}}+{{B(N_{\rm
bins})}\over{N_{\rm em}^{2}}}+\cdots.$ (57)

Here $\Delta I_{\infty}$ is the true value of $\Delta I$ which we would
observe if we could collect an infinite number of samples. The notation
reminds us that if we make bins along some continuous axis, then the size of
the corrections at finite $N_{\rm em}$ depend on the number of bins $N_{\rm
bins}$. With more bins the corrections are larger, which means that a naive
estimate with a fixed number of embryos will depend on the bin size. The hope
is that we can find a regime in which the extrapolated $\Delta I_{\infty}$ is
independent of $N_{\rm bins}$.

It is important that Eq. (57) is not just a guess, but a prediction that can
be derived theoretically. Theory also gives values for the coefficients $A$
and $B$, but these depend on details such as the independence of samples; the
form is more general. This suggests a strategy in which we vary the number of
embryos that we include in our analysis and look for the predicted systematic
dependence on $N_{\rm em}$. If we can see this behavior then we can feel
confident in fitting to Eq. (57) and extracting an estimate $\Delta
I_{\infty}$ [38].

This estimation procedure is illustrated by Fig. 4 in the main text and by
Fig. A5. When we vary the number of embryos that we include in our analysis,
we can choose at random from the total number available, so we have a path
also to estimating error bars (below). In Fig. 4A we analyze $\Delta I$ for
the spatial profiles of Hb expression using the Gaussian approximation of Eq.
(48). We have to make estimates of the variance as a function of the scaled
coordinate $x_{s}$ and the embryo length $L$. As explained in the main text,
we choose bins of $\Delta x_{s}=0.01$, consistent with the observed precision
of the pair-rule stripes and with earlier work [47, 33]. Along the $L$ axis we
use adaptive bins, that is bins with boundaries chosen so that the number of
embryos in each bin is as nearly equal as possible; these bins are chosen
based on the full experimental ensembles, and not readjusted as we choose
smaller samples at random.

In Figure 4A we have chosen $N_{\rm bins}=5$ adaptive bins along the $L$ axis,
and we see a clean depending of $\Delta I$ on $N_{\rm em}$ as predicted in Eq.
(57). The dependence on $N_{\rm em}$ is dominated by the linear term $A$,
although some curvature is detectable. The quality of the fit to Eq. (57) is
very good, and the extrapolation is to $\Delta I_{\infty}=0$ with a precision
of better than $0.01\,{\rm bits}$.

In Figure A5 we see how the plot of $\Delta I$ vs. $N_{\rm em}$ changes as we
vary $N_{\rm bins}=5,\,10,\,15,\,20$. With more bins, there are fewer embryos
in each bin, which means that the minimum number of embryos that we need for a
meaningful analysis is larger. Increasing the number of bins might reveal
otherwise hidden information, but also increases the size of systematic
errors. Comparing across the panels in Fig. A5 we see that at fixed $N_{\rm
em}$ the apparent $\Delta I$ increases with $N_{\rm bins}$, and if we didn’t
explore the full dependence on $N_{\rm em}$ we might be tempted to conclude
that we are indeed revealing extra information. But this is not the case,
since the plots at all values of $N_{\rm bins}$ extrapolate to zero within
error bars.

Figure A5: Consistent estimates of $\Delta I$ with varying $N_{\rm bins}$. (A)
Repeats the results of Fig 4A on $\Delta I$ vs. $N_{\rm em}$ for Hb, analyzed
with $N_{\rm bins}=5$ adaptive bins along the $L$ axis. (B) As in (A) with
$N_{\rm bins}=10$; (C) with $N_{\rm bins}=15$; and (D) with $N_{\rm bins}=20$.
Systematic errors are large at fixed $N_{\rm em}$ and increasing $N_{\rm
bins}$, but the extrapolation $N_{\rm em}\rightarrow\infty$ is within error
bars of $\Delta I=0$ in each case.

One useful test of these extrapolation methods is to be sure that we arrive at
zero information in those cases where we know that the answer must be zero. As
an example, if we randomly permute or shuffle the data we can break
correlations that lead to nonzero information. In this case, if we randomly
reassign lengths $L$ to the embryos, then we must have $\Delta I=0$. This is
illustrated in Fig. A6, using the example of Bcd. Here the real data
extrapolate to a nonzero value of $\Delta I$ (Fig. 5), and when we shuffle we
still see significantly nonzero values at $N_{\rm em}\sim 100$. But using the
whole $N_{\rm em}$ dependence we can see this extrapolates smoothly to zero,
as it should.

Figure A6: Recovering $\Delta I=0$ in shuffled data. We permute the lengths
$L$ of the embryos at random and repeat the analysis of the Bcd profiles.
While the real data extrapolate to nonzero $\Delta I$ (Fig 5C), here we
recover $\Delta I=0$ as expected.

An essential part of this analysis is the estimation of error bars. For
reasonably large $N_{\rm em}$ the systematic and random errors are additive,
and the variance of random errors scales $\propto 1/N_{\rm em}$ as usual. This
means that if we compute the variance of $\Delta I$ across random halves of
the data, and divide by two, we should have a good estimate of the variance in
$\Delta I$ based on all of the data. If this error bar $\sigma_{\Delta I}$ is
correct, and our extrapolation is consistent, then when the true $\Delta I$ is
zero, as with shuffled data, we can form a z-score,
$z=\Delta_{\infty}/\sigma_{\Delta I}$, and $z$ should be Gaussian with zero
mean and unit variance. We can test this because the extrapolation procedure
involves taking random subsamples of the data, each of which generates a
slightly different value of $\Delta I_{\infty}$. Figure A7 shows the
distribution $P(z)$ obtained from a shuffled version of the Hb data,
illustrating a good match to the expected Gaussian distribution; the deviation
is a bias toward smaller $z$, suggesting that our estimates of the error bars
may be a bit conservative.

Figure A7: Distribution of errors in $\Delta I$. The distribution of errors of
$\Delta I$ was estimated by repeating 5000 times the entire procedure leading
to Fig. 3A, on shuffled versions of the Hb data, where we know the true value
of $\Delta I$ is 0 bits. We calculate the z-score based on our estimates of
$\Delta I$ and error bars $\sigma_{\Delta I}$. Red histogram shows the
probability distribution $P(z)$ in bins of size $\Delta z=0.1$. Shaded area is
the uncertainty estimated through bootstrapping; black line is a Gaussian
distribution with $\langle z\rangle=0$ and $\langle z^{2}\rangle=1$.

Now that we have control over both the systematic and random errors in
estimating $\Delta I$, we summarize the results. We have done separate
analyses for each of the gap genes, for all four gap genes together, and for
the maternal input Bicoid. As we see in Fig. A8A, all results for the gap
genes are consistent with $\Delta I=0\,{\rm bits}$ within error bars, while
for Bicoid we find a significantly nonzero value. These quantities are perhaps
best expressed as fractions of the information conveyed about scaled position,
as shown in Fig. A8B; estimates of $I({\mathbf{g}}\rightarrow x_{s})$ are from
Refs. [47, 55].

Figure A8: Summary of deviations from scale invariance. (A) Extrapolated
$\Delta I$ with error bars for Bicoid, each gap gene individually, and the
combination of all gap genes. (B) Deviation from scale invariance as a
percentage of the information about scaled position [47, 55].

## References

  * Ishimatsu  _et al._ [2018] K. Ishimatsu, T. W. Hiscock, C. Z. M., D. W. K. Sari, K. Lischer, D. L. Richmond, Y. Bessho, T. Matsui, and S. G. Megason, Size-reduced embryos reveal a gradient scaling based mechanism for zebrafish somite formation, Development 145, dev161257 (2018).
  * Almuedo-Castillo  _et al._ [2018] M. Almuedo-Castillo, A. Bläßle, D. Mörsdorf, L. Marcon, G. H. Soh, K. W. Rogers, A. F. Schier, and P. Müller, Nature Cell Biology 20, 1032 (2018).
  * Leibovich  _et al._ [2020] A. Leibovich, T. Edri, S. L. Klein, S. A. Moody, and A. Fainsod, Natural size variation among embryos leads to the corresponding scaling in gene expression, Developmental Biology 462, 165 (2020).
  * Huxley [1932] J. S. Huxley,  _Problems of Relative Growth_ (Meuthen, London, 1932).
  * McMahon and Bonner [1983] T. A. McMahon and J. T. Bonner,  _On Size and Life_ (Scientific American Library, New York, 1983).
  * West  _et al._ [1997] G. B. West, J. H. Brown, and B. J. Enquist, A general model for the origin of allometric scaling laws in biology, Science 276, 122 (1997).
  * Kepler [2010] J. Kepler,  _The Six–Cornered Snowflake: A New Year’s Gift_ (Paul Dry Books, Philadelphia, PA, 2010) translated from the 1611 Edition by J Bromberg.
  * Bénard [1900] H. Bénard, Les tourbillions cellulaires dans une nappe liquide transportent de la chaleur par convection en regime permanent, Ann Chim Phys 7, 62 (1900).
  * Lord Rayleigh [1916] Lord Rayleigh, On convection currents in a horizontal layer of fluid, when the higher temperature is on the under side, Phil Mag 32, 529 (1916).
  * Flesselles  _et al._ [1991] J.-M. Flesselles, A. J. Simon, and A. J. Libchaber, Dynamics of one-dimensional interfaces: An experimentalist’s view, Advances in Physics 40, 1 (1991).
  * Cross and Hohenberg [1993] M. C. Cross and P. C. Hohenberg, Pattern formation outside of equilibrium, Reviews of Modern Physics 65, 851 (1993).
  * Lappa [2009] M. Lappa,  _Thermal Convection: Patterns, Evolution and Stability_ (John Wiley and Sons, New York, 2009).
  * Langer [1989] J. S. Langer, Dendrites, viscous fingers, and the theory of pattern formation, Science 243, 1150 (1989).
  * Turing [1952] A. M. Turing, The chemical basis of morphogenesis, Philos Trans R Soc Lond B 237, 37 (1952).
  * Scott and Carroll [1987] M. P. Scott and S. B. Carroll, The segmentation and homeotic gene network in early Drosophila development, Cell 51, P689 (1987).
  * Jaeger and Verd [2020] J. Jaeger and B. Verd, Dynamic positional information: Patterning mechanism versus precision in gradient-driven systems, Current Topics in Developmental Biology 137, 219 (2020).
  * Tkačik and Gregor [2021] G. Tkačik and T. Gregor, The many bits of positional information, Development 148, dev176065 (2021).
  * Houchmandzadeh  _et al._ [2002] B. Houchmandzadeh, E. Wieschaus, and S. Leibler, Establishment of developmental precision and proportions in the early Drosophila embryo, Nature 415, 798 (2002).
  * Surkova  _et al._ [2008] S. Surkova, D. Kosman, K. Kozlov, E. Myasnikova, A. A. Samsonova, A. Spirov, C. E. Vanario-Alonso, M. Samsonova, J. Reinitz, _et al._ , Characterization of the Drosophila segment determination morphome, Developmental Biology 313, 844 (2008).
  * Holloway  _et al._ [2006] D. M. Holloway, L. G. Harrison, D. Kosman, C. E. Vanario-Alonso, and A. V. Spirov, Analysis of pattern precision shows that drosophila segmentation develops substantial independence from gradients of maternal gene products, Developmental Dynamics 235, 2949 (2006).
  * Lott  _et al._ [2007] S. E. Lott, M. Kreitman, A. Palsson, E. Alekseeva, and M. Z. Ludwig, Canalization of segmentation and its evolution in Drosophila, Proceedings of the National Academy of Sciences (USA) 104, 10926 (2007).
  * Miles  _et al._ [2011] C. M. Miles, S. E. Lott, C. L. Luengo Hendriks, M. Z. Ludwig, Manu, C. L. Williams, and M. Kreitman, Artificial selection on egg size perturbs early pattern formation in Drosophila melanogaster, Evolution 65, 33 (2011).
  * Antonetti  _et al._ [2018] V. Antonetti, W. Bialek, T. Gregor, G. Muhaxheri, M. Petkova, and M. Scheeler, Precise spatial scaling in the early fly embryo, arXiv:1812.11384 (2018).
  * Driever and Nüsslein-Volhard [1988a] W. Driever and C. Nüsslein-Volhard, A gradient of bicoid protein in Drosophila embryos, Cell 54, 83 (1988a).
  * Driever and Nüsslein-Volhard [1988b] W. Driever and C. Nüsslein-Volhard, The bicoid protein determines position in the Drosophila embryo in a concentration-dependent manner, Cell 54, 95 (1988b).
  * Krotov  _et al._ [2014] D. Krotov, J. O. Dubuis, T. Gregor, and W. Bialek, Morphogenesis at criticality, Proceedings of the National Academy of Sciences (USA) 111, 3683 (2014).
  * McGough  _et al._ [2023] L. McGough, H. Casademunt, M. Nikolić, M. Petkova, T. Gergor, and W. Bialek, Finding the last bits of positional information, arXiv:2312.05963 (2023).
  * Manu  _et al._ [2009] Manu, S. Surkova, A. V. Spirov, V. V. Gurksy, H. Janssens, K. Ah-Ram, O. Radulescu, C. E. Vanario-Alonso, D. H. Sharp, M. Samsonova, and J. Reinitz, Canalization of gene expression and domain shifts in the Drosophila blastoderm by dynamical attractors, PLoS Computational Biology 5, e1000303 (2009).
  * Vakulenko  _et al._ [2009] S. Vakulenko, Manu, J. Reinitz, and O. Radulescu, Size regulation in the segmentation of Drosophila: Interacting interfaces between localized domains of gene expression ensure robust spatial patterning, Physical Review Letters 103, 168102 (2009).
  * Wolpert [1969] L. Wolpert, Positional information and the spatial pattern of cellular differentiation, Journal of Theoretical Biology 25, 1 (1969).
  * Lawrence [1992] P. A. Lawrence,  _The Making of a Fly: The Genetics of Animal Design_ (Blackwell Scientific Publications Ltd, Oxford, 1992).
  * Dubuis  _et al._ [2013a] J. O. Dubuis, R. Samanta, and T. Gregor, Accurate measurements of dynamics and reproducibility in small genetic networks, Molecular Systems Biology 9, 639 (2013a).
  * Petkova  _et al._ [2019] M. D. Petkova, G. Tkačik, W. Bialek, E. F. Wieschaus, and T. Gregor, Optimal decoding of cellular identities in a genetic network, Cell 176, 844 (2019).
  * Gregor  _et al._ [2007a] T. Gregor, E. F. Wieschaus, A. P. McGregor, W. Bialek, and D. W. Tank, Stability and nuclear dynamics of the bicoid morphogen gradient, Cell 130, 141 (2007a).
  * Smith [2015] E. M. Smith,  _Scaling and Regulation of Gene Expression in the Developing Fly Embryo_ , Ph.D. thesis, Princeton University (2015).
  * Shannon [1948] C. E. Shannon, A mathematical theory of communication, Bell System Technical Journal 27, 379 (1948).
  * Cover and Thomas [1991] T. M. Cover and J. A. Thomas,  _Elements of Information Theory_ (Wiley and Sons, New York, 1991).
  * Bialek [2012] W. Bialek,  _Biophysics: Searching for Principles_ (Princeton University Press, 2012).
  * Howard and ten Wolde [2005] M. Howard and P. R. ten Wolde, Finding the center reliably: Robust patterns of developmental gene expression, Physical review letters 95, 208103 (2005).
  * Houchmandzadeh  _et al._ [2005] B. Houchmandzadeh, E. Wieschaus, and S. Leibler, Precise domain specification in the developing Drosophila embryo, Physical Review E 72, 061920 (2005).
  * McHale  _et al._ [2006] P. McHale, W.-J. Rappel, and H. Levine, Embryonic pattern scaling achieved by oppositely directed morphogen gradients, Physical Biology 3, 107 (2006).
  * Čapek and Müller [2019] D. Čapek and P. Müller, Positional information and tissue scaling during development and regeneration, Development 146, dev177709 (2019).
  * Frasch  _et al._ [1988] M. Frasch, R. Warrior, J. Tugwood, and M. Levine, Molecular analysis of even-skipped mutants in Drosophila development, Genes and Development 2, 1824 (1988).
  * Jaeger  _et al._ [2004a] J. Jaeger, S. Surkova, M. Blagov, H. Janssens, D. Kosman, K. N. Kozlov, E. Myasnikova, C. E. Vanario-Alonso, M. Samsonova, D. H. Sharp, and J. Reinitz, Dynamic control of positional information in the early Drosophila embryo, Nature 430, 368 (2004a).
  * Jaeger  _et al._ [2004b] J. Jaeger, M. Blagov, D. Kosman, K. N. Kozlov, Manu, E. Myasnikova, S. Surkova, C. E. Vanario-Alonso, M. Samsonova, D. H. Sharp, and J. Reinitz, Dynamical analysis of regulatory interactions in the gap gene system of Drosophila melanogaster, Genetics 167, 1721 (2004b).
  * Bothma  _et al._ [2014] J. P. Bothma, H. G. Garcia, E. Esposito, G. Schlissel, T. Gregor, and M. Levine, Dynamic regulation of eve stripe 2 expression reveals transcriptional bursts in living Drosophila embryos, Proceedings of the National Academy of Sciences (USA) 111, 10598 (2014).
  * Dubuis  _et al._ [2013b] J. O. Dubuis, G. Tkačik, E. F. Wieschaus, T. Gregor, and W. Bialek, Positional information, in bits, Proceedings of the National Academy of Sciences (USA) 110, 16301 (2013b).
  * Vincent  _et al._ [1997] A. Vincent, J. T. Blankenship, and E. Wieschaus, Integration of the head and trunk segmentation systems controls cephalic furrow formation in Drosophila, Development 124, 3747 (1997).
  * Liu  _et al._ [2013] F. Liu, A. H. Morrison, and T. Gregor, Dynamic interpretation of maternal inputs by the Drosophila segmentation gene network, Proceedings of the National Academy of Sciences (USA) 110, 6724 (2013).
  * Jaeger [2011] J. Jaeger, The gap gene network, Cellular and Molecular Life Sciences 68, 243 (2011).
  * Kauffman  _et al._ [1978] S. A. Kauffman, R. M. Shymko, and K. Trabert, Control of sequential compartment formation in Drosophila: A uniform mechanism may control the locations of successive binary developmental commitments, Science 199, 259 (1978).
  * Meinhardt [1986] H. Meinhardt, Hierarchical inductions of cell states: A model for segmentation in Drosophila, Journal of Cell Science Supplement 4, 357 (1986).
  * Albert and Othmer [2003] R. Albert and H. G. Othmer, The topology of the regulatory interactions predicts the expression pattern of the segment polarity genes in Drosophila melanogaster, Journal of Theoretical Biology 223, 1 (2003).
  * Spirov and Holloway [2003] A. V. Spirov and D. M. Holloway, Making the body plan: Precision in the genetic hierarchy of Drosophila embryo segmentation, In Silico Biology 3, 89 (2003).
  * Tkačik  _et al._ [2015] G. Tkačik, J. O. Dubuis, M. D. Petkova, and T. Gregor, Positional information, positional error, and readout precision in morphogenesis: A mathematical framework, Genetics 199, 39 (2015).
  * Miller [1955] G. Miller, Note on the bias of information estimates, in  _Information Theory in Psychology II–B: Problems and Methods_ , edited by H. Quastler (Free Press, Glencoe IL, 1955) pp. 95–100.
  * Gregor  _et al._ [2005] T. Gregor, W. Bialek, R. R. de Ruyter van Steveninck, D. W. Tank, and E. F. Wieschaus, Diffusion and scaling during early embryonic pattern formation, Proceedings of the National Academy of Sciences 102, 18403 (2005).
  * Gregor  _et al._ [2008] T. Gregor, A. P. McGregor, and E. F. Wieschaus, Shape and function of the bicoid morphogen gradient in dipteran species with different sized embryos, Developmental Biology 316, 350 (2008).
  * Gregor  _et al._ [2007b] T. Gregor, D. W. Tank, E. F. Wieschaus, and W. Bialek, Probing the limits to positional information, Cell 130, 153 (2007b).
  * Petkova  _et al._ [2014] M. D. Petkova, S. C. Little, F. Liu, and T. Gregor, Maternal origins of developmental reproducibility, Current Biology 24, 1283 (2014).
  * Gierer and Meinhardt [1972] A. Gierer and H. Meinhardt, A theory of biological pattern formation, Kybernetik 12, 30 (1972).
  * Meinhardt [2008] H. Meinhardt, Models of biological pattern formation: From elementary steps to the organization of embryonic axes, Current Topics in Developmental Biology 81, 1 (2008).
  * Seyboldt  _et al._ [2022] R. Seyboldt, J. Lavoie, A. Henry, and P. François, Latent space of a small genetic network: Geometry of dynamics and information, Proceedings of the National Academy of Sciences (USA) 119, e2113651119 (2022).
  * Seung [1996] H. S. Seung, How the brain keeps the eyes still, Proceedings of the National Academy of Sciences (USA) 93, 13339 (1996).
  * Merle  _et al._ [2023] M. Merle, L. Friedman, C. Chureau, A. Shoushtarizadeh, and T. Gregor, Precise and scalable self-organization in mammalian pseudo-embryos, arXiv:2303.17522 (2023).


# Software Architecture Challenges in Integrating Hybrid Classical-Quantum
Systems

Vlad Stirbu  University of Jyväskylä
Jyväskylä, Finland
<EMAIL_ADDRESS>Tommi Mikkonen  University of Jyväskylä
Jyväskylä, Finland
<EMAIL_ADDRESS>
###### Abstract

The emergence of quantum computing proposes a revolutionary paradigm that can
radically transform numerous scientific and industrial application domains.
The ability of quantum computers to scale computations exponentially imply
better performance and efficiency for certain algorithmic tasks than current
computers provide. However, to gain benefit from such improvement, quantum
computers must be integrated with existing software systems, a process that is
not straightforward. In this paper, we investigate challenges that emerge from
building larger hybrid classical-quantum computers, and discuss some
approaches that could be employed to overcome these challenges.

###### Index Terms:

Quantum software, software architecture, classic-quantum systems

##  I Introduction

Quantum computers have demonstrated the potential to revolutionize various
fields, including cryptography, drug discovery, materials science, and machine
learning, by leveraging the principles of quantum mechanics. However, the
current generation of quantum computers, known as noisy intermediate-scale
quantum (NISQ) computers, suffer from noise and errors, making them
challenging to operate. Additionally, the development of quantum algorithms
requires specialized knowledge not readily available to the majority of
software professionals. These factors pose a significant entry barrier for
leveraging the unique capabilities of quantum systems.

For the existing base of business applications, classical computing has
already proven its capabilities across a diverse range of solutions. However,
some of the computations they must perform can be accelerated with quantum
computing, much like GPUs are used today. Therefore, quantum systems should
not function in isolation, but they must coexist and interoperate with
classical systems. To this end, software architects play a crucial role in
achieving seamless integration, while simultaneously designing systems that
effectively meet the unique requirements of businesses.

To address the challenges associated with this integration, this paper focuses
on designing hybrid systems that integrate quantum and classical computing,
aiming to overcome architectural, design, and operational hurdles. In doing
so, we look at the software development lifestyle, the technology stack of
hybrid classic-quantum systems, and deployment techniques used today.

##  II Background

The software development lifecycle (SDLC) of hybrid classic-quantum
applications consist of a multi-faceted approach, as depicted in Fig. 1. At
the top level, the classical software development process starts by
identifying user needs and deriving them into system requirements. These
requirements are transformed into a design and implemented. The result is
verified against the requirements and validated against user needs. Once the
software system enters the operational phase, any detected anomalies are used
to identify potential new system requirements, if necessary. A dedicated track
for quantum components is followed within the SDLC [1], specific to the
implementation of quantum technology. The requirements for these components
are converted into a design, which is subsequently implemented on classic
computers, verified on simulators or real quantum hardware, and integrated
into the larger software system. During the operational phase, the quantum
software components are executed on real hardware. Scheduling ensures
efficient utilization of scarce quantum hardware, while monitoring
capabilities enable the detection of anomalies throughout the process.

Figure 1: Software development lifecycle of a hybrid classical-quantum system

A typical hybrid classic-quantum software system is understood as a classical
program that has one or more software components that are implemented using
quantum technology, as depicted in Fig. 2. A quantum component utilises
quantum algorithms [2], that are transformed into quantum circuits using a
toolkit like Cirq111https://quantumai.google/cirq or
Qiskit222https://qiskit.org. The quantum circuit describes quantum
computations in a machine-independent language using quantum assembly (QASM)
[3]. This circuit is translated by a computer that controls the quantum
computer in a machine specific circuit and a sequence of pulses that control
the operation of individual hardware qubits [4]. Due to the scarcity or
quantum hardware and the process of preparing the individual runs, the quantum
task execution process is lengthy, having the characteristics of batch
processing in classical computing. In fact, techniques used in batch
processing, such as Slurm [5], can be used to implement this step, which adds
indirection to the underlying software architecture.

Figure 2: Quantum computing model: components and interfaces

##  III Architectural concerns

###  III-A Design – Algorithms, data structures and APIs

Quantum algorithms are designed specifically to take advantage of quantum
mechanics properties such as quantum superposition and entanglement. They
provide advantages over classical equivalents for specific areas, such as
factoring or linear search. Software architects should evaluate the
feasibility to achieve quantum advantage during the component requirements
phase of the SDLC. They must ensure that the needed computational resources
are available and that data can be mapped from the classic to quantum domains.
For example, TensorFlow Quantum333https://www.tensorflow.org/quantum is a
library for rapid prototyping of hybrid quantum-classical ML models that
focuses on quantum data and hybrid quantum-classical machine learning models.

The batch nature of the quantum task execution has a profound impact on the
software architecture of a hybrid classic-quantum system. The jobs submitted
for execution are queued and scheduled using fair-share policies. As the task
execution results are not available immediately, the software system should
favour asynchronous communication. Further, the system designers must consider
the security and privacy aspects of executing tasks on quantum hardware
infrastructure shared by several organizations.

###  III-B Operations – Implementation leak and resource allocation

Quantum application written using popular libraries like Cisq and Qiskit have
a monolith nature. They combine into a single imperative program the
application logic components, the general purpose quantum circuit design, the
quantum hardware selection (e.g. backend configuration), and the
transformation of the machine specific circuit that is actually executed. To
make the software architecture modular, the general purpose part needs to be
separated from the quantum backend. Essential backend information, such as the
quantum volume (the qubit connectedness), needs to be accessed at runtime so
that the actual hardware selection can be done dynamically based on dynamic
factors, like hardware availability (if there are multiple providers) and cost
estimates. For example, Kubernetes serves as an extensible orchestration
platform that enables efficient scheduling of classic computing jobs. The
capabilities of quantum computers can be exposed in this computing
environment, while the scheduler can be enhanced to efficiently handle quantum
jobs.

##  IV Conclusions and future steps

The fundamental differences in programming models and the varying levels of
maturity in tools and practices between the classical and quantum domains
makes their seamless integration difficult. To gain insights and firsthand
experience, we intend to collaborate with the users of
HELMI444https://docs.csc.fi/computing/quantum-computing/overview/ quantum
computer, in an effort to overcome the integration barriers between classical
and quantum computing.

## Acknowledgement

This work has been supported by the Academy of Finland (project DEQSE 349945)
and Business Finland (project TORQS 8582/31/2022).

## References

  * [1] B. Weder, J. Barzen, F. Leymann, and D. Vietz, Quantum Software Development Lifecycle, pp. 61–83.  Cham: Springer International Publishing, 2022.
  * [2] A. Montanaro, “Quantum algorithms: an overview,” npj Quantum Information, vol. 2, p. 15023, Jan 2016.
  * [3] A. Cross, A. Javadi-Abhari, T. Alexander, N. De Beaudrap, L. S. Bishop, S. Heidel, C. A. Ryan, P. Sivarajah, J. Smolin, J. M. Gambetta, and B. R. Johnson, “Openqasm 3: A broader and deeper quantum assembly language,” ACM Transactions on Quantum Computing, vol. 3, sep 2022.
  * [4] T. Alexander, N. Kanazawa, D. J. Egger, L. Capelluto, C. J. Wood, A. Javadi-Abhari, and D. C. McKay, “Qiskit pulse: programming quantum computers through the cloud with pulses,” Quantum Science and Technology, vol. 5, p. 044006, aug 2020.
  * [5] A. B. Yoo, M. A. Jette, and M. Grondona, “Slurm: Simple Linux utility for resource management,” in Job Scheduling Strategies for Parallel Processing: 9th International Workshop, JSSPP 2003, Seattle, WA, USA, June 24, 2003. Revised Paper 9, pp. 44–60, Springer, 2003.