# Finding the last bits of positional information

Lauren McGough,a,b Helena Casademunt,c Miloš Nikolić,a Mariela D. Petkova,d
Thomas Gregor,a,e and William Bialek,a aJoseph Henry Laboratories of Physics
and Lewis–Sigler Institute for Integrative Genomics, Princeton University,
Princeton NJ 08544 USA
bDepartment of Ecology and Evolution, The University of Chicago, Chicago IL
60637
cDepartment of Physics and dProgram in Biophysics, Harvard University,
Cambridge MA 02138
eDepartment of Developmental and Stem Cell Biology UMR3738, Institut Pasteur,
75015 Paris, France

###### Abstract

In a developing embryo, information about the position of cells is encoded in
the concentrations of “morphogen” molecules. In the fruit fly, the local
concentrations of just a handful of proteins encoded by the gap genes are
sufficient to specify position with a precision comparable to the spacing
between cells along the anterior–posterior axis. This matches the precision of
downstream events such as the striped patterns of expression in the pair-rule
genes, but is not quite sufficient to define unique identities for individual
cells. We demonstrate theoretically that this information gap can be bridged
if positional errors are spatially correlated, with relatively long
correlation lengths. We then show experimentally that these correlations are
present, with the required strength, in the fluctuating positions of the pair-
rule stripes, and this can be traced back to the gap genes. Taking account of
these correlations, the available information matches the information needed
for unique cellular specification, within error bars of $\sim 2\%$. These
observation support a precisionist view of information flow through the
underlying genetic networks, in which accurate signals are available from the
start and preserved as they are transformed into the final spatial patterns.

##  I Introduction

During the development of an embryo, cell fates are determined in part by the
concentrations of specific morphogen molecules that carry information about
position turing_52 ; wolpert_69 ; tkacik+gregor_21 . For the early stages of
fruit fly development, all of these molecules have been identified
nusslein+wieschaus_80 ; lawrence_92 ; ew+cnv_16 . For patterning along the
main body axis, spanning from anterior to posterior (AP), information flows
from primary maternal morphogens to an interacting network of gap genes to the
pair-rule genes riverapomar_96 ; jaeger_11 , whose striped patterns of
expression provide a precursor of the segmented body plan in the fully
developed organism, visible within three hours after the egg is laid (Fig. 1).
It has been known for some time that, at this stage in development,
essentially every cell “knows” it’s fate gergen+al_86 , so it is natural to
ask how this information is encoded, quantitatively, in the concentrations of
the relevant morphogens.

Figure 1: Segmented Drosophila body plan. (A) Brightfield color image of a
$5\,$mm long $3^{\rm rd}$ instar larva of the fruit fly Drosophila
melanogaster bugs  with clearly visible segments. (B) An optical section
through an embryo stained for three of the “pair-rule” proteins, $50\,$min
into nuclear cycle 14 ($\sim 3\,{\rm h}$ after oviposition), showing striped
patterns that align with the body segments; data from Ref petkova+al_19 . (C)
As in (B), from multiple embryos, illustrating the pattern reproducibility.
Time in nuclear cycle 14 indicated at bottom right of each profile.

The expression levels of the gap genes provide enough information to specify
the positions of individual cells with an accuracy $\sim 1\%$ of the embryo’s
length dubuis+al_13 . This matches the precision with which the stripes of
pair-rule expression are positioned, and the precision of macroscopic
developmental events such as the formation of the cephalic furrow liu+al_13 .
Further, the algorithm that extracts optimal estimates of position from the
expression levels of the gap genes also predicts, quantitatively, the
distortions of the striped pattern in mutant flies with deletions of the
maternal inputs petkova+al_19 . At the moment when pair-rule stripes are fully
formed, just before gastrulation, there are fewer than one hundred rows of
cells along the length of the embryo, so it is tempting to think that
positional signals with $1\%$ accuracy define unique cellular identities. In
fact, this is not quite correct dubuis+al_13 : if each cell makes independent
positional errors drawn from a Gaussian distribution, then there is a small
but significant probability that neighboring cells will get “crossed signals,”
driving errors in cell fate determination.

The small difference between $1\%$ positional errors and unique cellular
identities provides an interesting test case in the search for a more
quantitative understanding of living systems. In physics, we are used to the
idea that small quantitative discrepancies can be signs of qualitatively new
ideas or mechanisms. But in complex biological systems one might worry that
small discrepancies reflect experimental errors or over–simplifications in
interpretation. If correct, these concerns would limit our ambitions for
quantitative theory in the physics tradition. However the small discrepancies
need to be re–examined in light of dramatic improvements in experimental
precision yu+al_06 ; taniguchi+al_10 ; dubuis+al_13b .

Here we take the small quantitative discrepancy in positional information
seriously. On the theoretical side, we clarify the problem, defining an
“information gap,” and show that this gap can be closed if errors in the
positional signals are spatially correlated over relatively long distances.
Early work by Lott and colleagues lott+al_07  detected such correlations in
mRNA levels of gap and pair-rule genes; subsequent work found that noise in
different combinations of protein levels in the gap gene network are
correlated significantly over the entire length of the embryo krotov+al_14 .
On the experimental side we re–examine these correlations, measuring the
positions of stripes in the concentrations of pair-rule proteins. We find that
the extent of these correlations is what is needed to close the information
gap between positional errors and unique cellular identities, quantitatively.

##  II Defining the problem

In the early fly embryo, cells have access to the concentrations of
morphogens, and these concentrations are continuously graded. From these
concentrations, it is possible to decode an estimate of position, which we
label as $\hat{x}_{\rm n}$ in cell $\rm n$ petkova+al_19 . We expect that
these estimates are correct on average, so that $\langle\hat{x}_{\rm
n}\rangle={\rm n}L/N$, where there are $N$ cells along the length $L$ of the
embryo.111For simplicity we imagine that the problem is one–dimensional so
that cells need to know their position only along one axis. In the early fly
embryo, patterning signals along the two major axes are largely independent
nusslein-volhard_89 ; nusslein-volhard_91 , justifying this simplification.
However the signals are noisy, so decoding in one cell will have errors,

$\displaystyle\hat{x}_{\rm n}$ $\displaystyle=$ $\displaystyle{\rm
n}L/N+\delta x_{\rm n},$ (1) $\displaystyle\langle(\delta x_{\rm
n})^{2}\rangle$ $\displaystyle=$ $\displaystyle\sigma_{x}^{2}.$ (2)

For simplicity, but guided by the experimental observations dubuis+al_13 ;
tkacik+al_15 ; petkova+al_19 , we assume that $\sigma_{x}$ is the same for all
cells and that the distribution of $\delta x_{\rm n}$ is Gaussian. Here we are
interested in the question of whether cells get signals that define the
correct ordering along the axis so that $\hat{x}_{\rm n+1}>\hat{x}_{\rm n}$
for all cells, or whether they can get “crossed signals” such that
$\hat{x}_{\rm n+1}<\hat{x}_{\rm n}$.

If we look at two neighboring cells, then the probability of incorrect
ordering is

$P_{\rm cross}\equiv{\rm Pr}(\hat{x}_{\rm n+1}<\hat{x}_{\rm n}).$ (3)

To find the probability of a wrong ordering we can take a look at the
distribution of the distance to the next cell $y=\hat{x}_{\rm
n+1}-\hat{x}_{\rm n}$. But since $\hat{x}_{\rm n+1}$ and $\hat{x}_{\rm n}$
both are Gaussian, their difference $y$ is also Gaussian, with mean equal to
$\langle y\rangle=L/N$. If the noise is independent in each cell, then the
variance of this difference signal will be $\langle(\delta
y)^{2}\rangle=2\sigma_{x}^{2}$. Incorrect ordering happens when $y<0$, which
then has probability

$\displaystyle P_{\rm cross}$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{0}{{dy}\over{\sqrt{4\pi\sigma_{x}^{2}}}}e^{-(y-L/N)^{2}/4\sigma_{x}^{2}}$
(4) $\displaystyle=$
$\displaystyle{1\over\sqrt{4\pi}}\int_{1/z}^{\infty}dx\,e^{-x^{2}/4},$ (5)

with $z=\sigma_{x}(N/L)$, as shown in Fig. 2. If positional errors are
comparable to the spacing between cells, $\sigma_{x}\sim L/N$, the probability
of an error is nearly 24%; for the experimental value
$\sigma_{x}\sim(0.74)L/N$ dubuis+al_13 , crossed signals will occur in $\sim
16\%$ of cells. With $N\sim 74\pm 5$ rows of cells along the AP axis
dubuis+al_13 , the probability that all signals come in the right order would
be vanishingly small.222This uncertainty in $N$ may seem large, but what will
matter below is the information required to specify unique cellular
identities, $I_{\rm unique}=\log_{2}N$. Although $\delta N/N$ is nearly ten
percent, $\delta I_{\rm unique}/I_{\rm unique}$ is less than two percent.

Figure 2: Probability of “crossed signals” between two neighboring cells as a
function of the positional error, assuming that noise is independent in each
cell [Eq (5)]. Dashed vertical line marks the experimental value of positional
noise, $\sigma_{x}\sim 0.01L$, which corresponds to less than the mean
distance between neighboring cells $L/N$ dubuis+al_13 .

This failure to specify unique cellular identities can be given a simple
information-theoretic interpretation. To specify one cell uniquely out of $N$
requires $I_{\rm unique}=\log_{2}N\,{\rm bits}$ of information shannon_48 ;
bialek_12 . On the other hand, if we have signals that represent a continuous
position $x$ drawn uniformly from the range $0<x\leq L$, and these signals
have Gaussian noise with (small) standard deviation $\sigma_{x}$, as described
above, then the amount of information the signal conveys about position is

$I_{\rm position}=\log_{2}L-\log_{2}\left(\sqrt{2\pi e}\sigma_{x}\right),$ (6)

where the first term is the entropy of the uniform distribution of positions
and the second term is the entropy of the Gaussian noise distribution
bialek_12 . Combining these we can define an “information gap”

$I_{\rm gap}\equiv I_{\rm unique}-I_{\rm
position}=\log_{2}\left({{N\sigma_{x}}\over L}\sqrt{2\pi e}\right).$ (7)

As discussed below, we obtain a more accurate estimate of the information gap
by averaging over measurements of $\sigma_{x}$ at multiple points along the
embryo, defined by the pair rule stripes, and we find $I_{\rm gap}=1.39\pm
0.08\,{\rm bits}$ (Appendix A). Importantly this gap is measured per cell: it
is not that the embryo is missing $\sim 1.4\,{\rm bits}$ of information, but
rather that every cell is missing this information.

##  III Extra information from correlations: Theory

In order to address this information gap directly, we leverage the concept
that correlated noise facilitates enhanced information transmission. While
correlated noise is typically viewed as challenging due to its resistance to
averaging, in the context of neighboring cells making correlated errors in
position, it mitigates the probability of receiving “crossed signals,” as
previously defined. Here we develop these considerations more formally.

Information is roughly the difference in entropy between the signal and the
noise, where entropy measures the (log) volume in phase space that is occupied
by a set of points. When random variables become correlated, the volume and
hence the entropy is reduced, even if the variances of the individual
variables are unchanged. In our example, with correlations, the full pattern
of points $\\{\hat{x}_{1},\,\hat{x}_{2},\,\cdots,\,\hat{x}_{N}\\}$ fills a
smaller volume in the space $[0,L]^{N}$ of possible positions for all the
cells, and thus the embryo as a whole has access to more positional
information.

More formally, we can define the correlation matrix $C$,

$\langle\delta x_{\rm n}\delta x_{\rm m}\rangle=\sigma_{x}^{2}{C}_{\rm nm},$
(8)

with diagonal elements $C_{\rm nn}=1$. Assuming again that the noise $\delta
x_{\rm n}$ is Gaussian, the reduction in noise entropy for the entire set of
variables $\\{\delta x_{\rm n}\\}$ is given by the determinant of this matrix
bialek_12 ,

$\Delta S=-\frac{1}{2}\log_{2}\det C\,\,\,\,{\rm bits},$ (9)

and this reduction in entropy is the gain in information. Entropy is an
extensive quantity, so that when $N$ is large the information gain per cell
$I_{\rm extra}=\Delta S/N$ is finite. Can $I_{\rm extra}$ be large enough to
compensate for the information gap $I_{\rm gap}$?

We expect that the correlation between fluctuations of positional signals in
different cells depends on their spatial separation. Then $C_{\rm nm}$ is a
function of the distance between cells $\rm n$ and $\rm m$, $d_{\rm nm}=|{\rm
n}-{\rm m}|L/N$. A natural functional form is an exponential decay of
correlations,

$C_{\rm nm}=e^{-d_{\rm nm}/\xi},$ (10)

with correlation length $\xi$. This is what we would see if signals were
encoded in the gradient of a single molecular species that has a lifetime
$\tau$ and diffusion constant $D$, with $\xi=\sqrt{D\tau}$. Although this is
over–simplified, it is useful for building intuition about how the range of
correlations determines the additional information. Within this model it is
straightforward to evaluate $\Delta S$ numerically, with results shown in Fig.
3A.

Figure 3: Extra information from correlations, as a function of the
correlation length. (A) Numerical results for $N=50$ and $N=100$ from Eq (9)
with the correlation matrix in Eq (10); analytic results for
$N\rightarrow\infty$ from Eq (15). Compare with the information gap from
Appendix A (solid black line bracketed by dashed error bars). (B) Probability
of at least two signals being “crossed,” $\hat{x}_{{\rm n}+1}<\hat{x}_{\rm n}$
in a line of $N=74$ cells, with $\sigma_{x}/L=0.01$.

We can also give an analytic theory for $\Delta S$ in the large $N$ limit,
leading to Eq (15) and the red line in Fig. 3. If we define eigenvalues and
eigenvectors of the matrix $C_{\rm nm}$,

$\sum_{\rm m}C_{\rm nm}\phi_{\rm m}^{\mu}=\lambda_{\mu}\phi_{\rm n}^{\mu},$
(11)

then we have

$\Delta S=-\frac{1}{2}\sum_{\mu}\log_{2}\lambda_{\mu}\,\,\,\,{\rm bits}.$ (12)

In the limit of large $N$ at fixed $N/L$, the ends of the embryo are far away,
and there is an effective translation invariance. This means that the
eigenvectors $\phi_{\rm n}^{\mu}$ are complex exponentials, $\phi_{\rm
n}^{\mu}\propto\exp(iq_{\mu}{\rm n})$, or equivalently that the matrix $C_{\rm
nm}$ is diagonalized by a discrete Fourier transform;333The discreteness is
important. If we take a continuum limit, so that the sum in Eq (13) becomes an
integral, the calculation is a bit simpler but leads to a significant
over–estimate of $\Delta S$, even at large values of $\xi N/L$.  allowed
values of $q_{\mu}$ are in the interval $-\pi\leq q<\pi$. Then as
$N\rightarrow\infty$ we find the eigenvalues

$\lambda(q)\rightarrow\sum_{{\rm n}=-\infty}^{\infty}e^{-|{\rm
n}|L/N\xi}e^{iq{\rm n}}={{\sinh(L/N\xi)}\over{\cosh(L/N\xi)-\cos(q)}},$ (13)

and the change in entropy

$\displaystyle\Delta S/N$ $\displaystyle\rightarrow$ $\displaystyle-{1\over
2}\int_{-\pi}^{\pi}{{dq}\over{2\pi}}\log_{2}\lambda(q)$ (14) $\displaystyle=$
$\displaystyle-{1\over
2}\log_{2}\left[\frac{2\sinh(L/N\xi)}{\sinh(L/N\xi)+\cosh(L/N\xi)}\right].$
(15)

In Fig. 3A we see that this analytic result agrees with numerical results at
$N=50$ and $N=100$, which agree with one another, confirming that the fly
embryo is large enough for the entropy to be extensive. We conclude that an
information gap of $\sim$1.4 bits can be closed if correlations extend over
distances $\xi\sim 13(L/N)\sim 0.18L$. Lott and colleagues saw significant
correlations across this range of distances for all the genes that they probed
lott+al_07 , and combinations of gap gene protein levels have even longer
correlation lengths krotov+al_14 .

Beyond the perhaps abstract information theoretic measures, we can evaluate
the probability that all cells receive signals that are in the correct order,
that is $\hat{x}_{\rm n+1}>\hat{x}_{\rm n}$ for all ${\rm
n}=1,\,2\,\cdots,\,N$. If correlations extend over a distance $\xi\sim
13(L/N)$, then proper ordering will occur in more than $99\%$ of embryos, as
illustrated in Fig. 3B.

##  IV Extra information from correlations: Experiment

Taking the information gap seriously, we predict that the noise in positional
signals should be correlated over distances $\xi\sim 0.2L$. These distances
are long compared to the separation between neighboring cells. The first
indication that such correlations exist came from experiments marking the
boundaries of gene expression domains as seen through measurements of mRNA for
selected gap genes and the pair rule gene eve lott+al_07 . At the same time,
it was reported that fluctuations in the concentration of a single gap gene
product protein are correlated only over short distances gregor+al_07b .
Analyzing simultaneous measurement on protein concentrations of four gap genes
demonstrated that different combinations or modes of the network have
different correlation lengths krotov+al_14 ; the longest correlation lengths
are a significant fraction of the length of the embryo. Finally, early
analyses showed that errors in relative position are smaller than errors in
absolute position dubuis+al_13 . All of this suggests that the noise in
positional signals is spatially correlated. Can we make this statement more
quantitative?

We analyze the experiments in Ref petkova+al_19 , which used
immunofluorescence stainings to measure spatial profiles of protein
concentration for three of the pair-rule genes eve, prd, and rnt (Fig. 1). The
data include $N_{\rm em}=109$ embryos, fixed and stained in the time window
from 35 to 60 min after the start of nuclear cycle 14. This is the period of
cellularization, and as in previous work, the progress of the cellularization
membrane provides a time marker with an accuracy of up to one minute
dubuis+al_13b . For each of the three genes, the seven peaks in the striped
concentration profile can be found automatically, and their locations vary
linearly with time throughout this period antonetti+al_18 . If we don’t
correct for this systematic dynamical behavior, the variance of stripe
positions will be large and their fluctuations will be correlated,
artificially. We consider the noise in position to be the deviation from the
best fit linear relation for each individual stripe marker. The standard
deviations then are consistently slightly below $\sigma_{x}\sim 0.01L$, and
the distribution of fluctuations is well approximated by a Gaussian. These
results agree with previous work dubuis+al_13 ; petkova+al_19 ;
antonetti+al_18 , and are summarized in Appendix A.

Before analyzing correlations, we can use these data to make a more precise
estimate of the information gap. If each cell has access to a positional
signal with errors $\sigma_{x}({\rm n})$, that might vary with $\rm n$, the
average positional information available to a single cell is

$I_{\rm position}=\log_{2}L-{\bigg{\langle}}\log_{2}\left[\sqrt{2\pi
e}\sigma_{x}({\rm n})\right]{\bigg{\rangle}}_{\rm n},$ (16)

where $\langle\cdots\rangle_{\rm n}$ denotes an average over cells,
generalizing Eq (6). Rather than making inferences about single cells, we have
direct access to the signals that mark the locations of the stripes in the
expression of three pair-rule genes, for a total of 21 features spread across
half the AP axis. The mean separation between the nearest stripes is
$\Delta\bar{x}=0.023L$, just a few times larger than the spacing between
cells. Rather than introducing a model that would interpolate, we take the
stripe positions themselves as the signals $x_{\rm n}$, now with ${\rm
n}=1,\,2\,\cdots,\,21$, and the average in Eq (16) becomes an average over
stripes.

The challenge in evaluating the positional information is that random errors
in our estimates of the errors $\sigma_{x}({\rm n})$ become systematic errors
in estimates of information. This problem of systematic errors was appreciated
in the very first efforts to use information theoretic concepts to analyze
biological experiments miller_55 . The analysis of neural codes has been an
important testing ground for methods to address these errors panzeri+treves_95
; strong+al_98 ; paninski_03 ; for a review see Appendix A.8 of Ref bialek_12
. The approach we take here uses the fact that naive entropy estimates depend
systematically on the size of the sample; if we can detect this systematic
dependence then we can extrapolate to infinite data, as described in Appendix
A. The result is that $I_{\rm gap}=1.39\pm 0.08\,{\rm bits/cell}$.

The idea of positional information is that cells have access to a signal that
represents position along the axis of the embryo wolpert_69 ; tkacik+al_15 .
In the discussion above we have taken this idea very seriously, identifying
the signal in each cell as $\hat{x}_{\rm n}$. But the signals we observe are
the positions of stripes in three different pair-rule genes, and the different
stripes for each gene are controlled by different enhancers responding to
distinct combinations of transcription factors. We need to test the hypothesis
that these multidimensional molecular concentrations encode a single
positional variable.

Figure 4: Correlations between noise in peak positions of the eve, run, and
prd stripe patterns, as a function of the mean separation between stripes.
Error bars estimated from the standard deviation across random halves of the
data. With three genes, each having seven stripes, we observe $(21\times
20)/2=210$ distinct elements of the correlation matrix $C_{\rm nm}$. Solid red
line is a smooth curve to guide the eye.

We are looking at fluctuations in the positions of the stripes, $\delta x_{\rm
n}$. Fig. 4 shows the elements of the correlation matrix

$C_{\rm nm}\equiv{{\langle\delta x_{\rm n}\delta x_{\rm
m}\rangle}\over{[\langle(\delta x_{\rm n})^{2}\rangle\langle(\delta x_{\rm
m})^{2}\rangle]^{1/2}}},$ (17)

as a function of the mean separation $\Delta{\bar{x}}_{\rm nm}$ between
stripes $\rm n$ and $\rm m$. We see that, within experimental error, the
correlations really are a function of distance. There is no obvious pattern
linked to the identity of the enhancers that control these different features,
or to the identity of the transcription factors to which the enhancers
respond: nearby stripes are highly correlated, the decay of correlations with
distance is the same whether we are looking at correlations between the same
or different genes, and different pairs of stripes with same mean separation
have the same correlation. This suggests that, as in the theoretical
discussion above, we can think about an abstract positional signal that is
transmitted to each cell and controls the placement of the pair-rule stripes.
Correspondingly, there are strong indications that the correlations are
inherited from the structure of the noise in gap gene expression (Appendix C).

Qualitatively, the correlations that we see in Fig. 4 decay over distances
$\xi\sim 0.15L$, consistent with the scale needed to close the information
gap, and with early measurements lott+al_07 . Quantitatively, the decay of
correlations is not well described by a single exponential function of
distance, so we cannot simply transcribe the predictions of the theory.
Instead, we would like to make a direct estimate of the positional information
from the data. Conceptually this is simple: we estimate the correlation matrix
from the data, then compute the (log) determinant of this matrix following Eq
(9). As with the information gap itself (above), the problem is that random
errors in our estimates of individual matrix elements become systematic errors
in the entropy. We follow the same strategy of identifying the dependence of
this error on the number of embryos that we include in our analysis and
extrapolating to large data sets, as described in Appendix B.

Figure 5: Extra information in correlations per cell, $\Delta S/N$, computed
from the observed correlations in pair-rule stripe fluctuations $C_{\rm nm}$,
including different numbers of contiguous stripes. Circles and error bars
(blue) are the extrapolated estimates from Appendix B. Beyond $N\sim 14$
stripes there is a plateau $\Delta S/N=1.51\pm 0.08\,{\rm bits/cell}$,
bracketed by the dashed lines. Square and error bars (red) are the best
estimate of the information gap $I_{\rm gap}=1.39\pm 0.08\,{\rm bits/cell}$
from Appendix A.

By definition, to see the extra information hidden in correlations we have to
look at the positions of multiple stripes. We start with two neighboring
stripes, and gradually work out toward all twenty–one stripes. We see in Fig.
5 that beyond $N\sim 10$ stripes, the information per stripe reaches a plateau
at $\Delta S/N=1.51\pm 0.08\,{\rm bits/stripe}$. This agrees, within
experimental error, with our estimate of the information gap $I_{\rm
gap}=1.39\pm 0.08\,{\rm bits/cell}$.

##  V Discussion

There is strong evidence that, early in embryonic development, each cell
acquires a distinct identity gergen+al_86 ; it is less clear how this
information is encoded. In the fruit fly embryo, positional information along
the anterior–posterior axis is orchestrated through a sequential cascade
involving three primary maternal inputs, a select number of gap genes, and the
pair rule genes. The conventional perspective suggests that the information
flow through this cascade entails a gradual refinement, with noisy inputs
ultimately generating a precise and reproducible pattern arias+hayward_06 ;
lacalli_22 , in the spirit of the Waddington landscape waddington_57 .

In contrast to the picture of noisy inputs and precise outputs, at least one
maternal input itself exhibits a high level of precision, consistently
reproducible across embryos gregor+al_07b ; petkova+al_14 . Moreover, the
expression levels of gap genes within a single cell prove sufficient to
determine positions with an error smaller the distance between neighboring
cells dubuis+al_13 ; petkova+al_19 . Notably, this precision agrees with that
observed in downstream events such as the pair-rule stripes. In parallel,
crucial developmental events exhibit highly reproducible temporal trajectories
ferree+al_16 . These quantitative observations challenge the conventional view
of refinement and error correction, supporting instead a precisionist
perspective in which locally available information is processed and preserved
with near optimal efficiency. Given that all relevant molecules are present at
low copy numbers, this places significant constraints on the architecture of
the underlying networks ferree+al_16 ; tkacik+al_09 ; walczak+al_10 ;
sokolowski+al_23 .

Despite their precision, local signals in the fly embryo do not quite provide
enough information to uniquely specify all $N=74\pm 5$ cellular identities
along the AP axis, $I_{\rm unique}=\log_{2}N$: errors in the position that a
cell can infer from molecular concentrations come from a distribution, and
distributions have tails dubuis+al_13 . The result is that there is a
substantial ($\sim 22\%$) gap between the information provided by the gap
genes, or the pair-rule stripes, and $I_{\rm unique}$.

Previous measurements have characterized the noise in local estimates of
position for each cell individually. But there are many hints from previous
work that this noise is correlated lott+al_07 ; dubuis+al_13 ; krotov+al_14 .
Extra information can be hiding in these correlations, and we have seen in
§III that if correlations extend over distances $\xi\sim 0.15L$ then this
would be enough to close the information gap. This prompts a more detailed
examination of the noise correlations, which really do seem to be a function
of distance independent of gene identity (Fig. 4).

The perhaps surprising conclusion of §IV is that the extra information
contained in the correlations, $\Delta S/N$, matches the information gap
$I_{\rm gap}$ to within a few percent of $I_{\rm unique}$, with the remaining
difference essentially equal to our error bars:

$I_{\rm gap}-\Delta S/N=(-0.019\pm 0.018)I_{\rm unique}.$ (18)

This agreement supports, strongly, the precisionist view of information flow
in this system.

Historically, the lack of precise data on gene expression levels, with
uncertainties extending to factors of two, led to skepticism regarding the
relevance of more refined measurements to general mechanisms of genetic
control. These expectations stood in contrast, for example, to our
understanding of signaling in rod photoreceptors, where the quantitative
reproducibility of responses to single molecular events provides important
constraints on the underlying biochemical mechanisms doan+al_06 .

The fly embryo has provided a laboratory within which to explore the precision
vs. noisiness in the function of an intact living system. We have seen
reproducible protein and mRNA concentrations across embryos with an accuracy
of $10\%$ gregor+al_07b ; dubuis+al_13b ; petkova+al_14 , and these
concentrations encode position with an accuracy of $\sim 1\%$ of the embryo’s
length dubuis+al_13 ; petkova+al_19 ; tkacik+al_15 . The current study adds a
layer to this understanding, demonstrating that the available positional
information, including the subtle effects of correlated noise, matches the
threshold for specifying unique cellular identities, and this match itself has
an accuracy of just a few percent. Beyond the fly embryo, these results
suggest a more general conclusion: quantitative measurements in living systems
merit serious consideration, even at high precision, as in other areas of
physics.

###### Acknowledgements.

We thank Eric Wieschaus for many inspiring discussions. This work was
supported in part by the US National Science Foundation, through the Center
for the Physics of Biological Function (PHY–1734030); by National Institutes
of Health Grant R01GM097275; by the Howard Hughes Medical Institute; and by
the Simons and John Simon Guggenheim Memorial Foundations.

##  Appendix A Statistics of individual stripes

The raw data for our analyses are the profiles of fluorescence intensity vs
position along the length of the embryo, as in Fig. 1. These embryos have been
fixed and stained with antibodies against the proteins encoded by the pair
rule genes eve, prd, and rnt, and fluorescently tagged antibodies against
those antibodies petkova+al_19 . Independent experiments demonstrate that
these classical staining methods, used carefully, yield fluorescence
intensities that are linear in protein concentrations dubuis+al_13b . The data
set used here, which contains a large number of wild type embryos, comes from
Ref petkova+al_19 .

We briefly summarize the imaging protocol and describe the procedure for
localizing the stripe positions. Images are taken in the midsaggital plane
showing a row of nuclei along the dorsal and ventral side of the embryo. For
consistency and to avoid geometric distortion, we focus on the dorsal
profiles, as was done previously. In order to include the entire embryos in a
single image, large field-of-view images, with pixel size $445\,{\rm nm}$ are
acquired with a $20\times$ $0.7$NA objective on a Leica SP5 confocal
microscope. Fluorescence intensity is averaged inside a sliding window of the
size of a nucleus and the position of the window center is recorded. In a
given embryo, positions of the 7 stripes are first roughly identified by
finding local maxima in the profile of an individual embryo. To make this
quantitative, we tried several methods. First, we used an iterative procedure
in which the mean peak shape is used as a template antonetti+al_18 . Second,
we fitted a model of seven Gaussians with variable amplitudes and widths to
the entire profile. Finally, we fit individual Gaussians to each stripe, using
a window centered on the local maximum with width of 5% embryo length. These
methods give consistent results, and importantly global fits do not generate
larger correlations than local fits. In the end we use the local Gaussian
fits, as in Fig. A1A.

Figure A1: Pair rule stripe positions. (A) Concentration of Eve protein in a
single embryo. Colored circles indicate regions which were fitted with a
Gaussian function to calculate the stripe position. Each stripe is fitted
individually, with fits shown in red. Red triangles indicate centers of each
fitted peak. (B) Stripe positions as a function of time in the nuclear cycle
14. Linear fits from Eq (19) are shown as black lines. (C) Peak positions
$x_{\rm n}(t_{0})$ corrected to $t_{0}=45$ min. (D) Positional error of the
pair rule stripes. Magnitude of the error $\sigma_{x}({\rm n})$ is plotted
against the mean position $\bar{x}_{\rm n}$ for each of the eve, prd, and rnt
stripes. Errors in $\bar{x}_{\rm n}$ are standard errors of the mean; errors
in $\sigma_{x}$ are standard deviations across random halves of the data.
Dashed line marks the rough estimate $\sigma_{x}/L\sim 0.01$.

The age of embryos is estimated to 1 minute precision in nuclear cycle 14 by
measuring the length of the cellularization membrane dubuis+al_13 . At 30 min
into this cycle, the stripes of _prd_ first start to become visible and the
other two genes have a well defined stripes by that time, so we confine our
attention to $t>30\,{\rm min}$.

Stripe patterns are dynamic, with positions that depend on time. If we don’t
take account of this systematic variation, then across an ensemble of embryos
with different ages we would see artificial correlations among fluctuations in
stripe position. Stripe movement is small, however, and we can use a linear
fit (separately for each of the 21 stripes) across the population of embryos,

$x_{\rm n}(t)=x_{\rm n}(t_{0})+s_{\rm n}(t-t_{0}).$ (19)

Results are shown in Fig. A1B and C. For each embryo we find an equivalent
position of all the stripes at a reference time $t_{0}=45$ min antonetti+al_18
.

With $x_{\rm n}$ the positions of each pair rule stripe, we have the mean and
variance

$\displaystyle\bar{x}_{\rm n}$ $\displaystyle=$ $\displaystyle\langle x_{\rm
n}\rangle$ (20) $\displaystyle\sigma_{x}^{2}({\rm n})$ $\displaystyle=$
$\displaystyle\langle(x_{\rm n}-\bar{x}_{\rm n})^{2}\rangle,$ (21)

where $\langle\cdots\rangle$ denotes an average over our complete experimental
ensemble of $N_{\rm em}=109$ embryos. Results are shown in Fig. A1 D, where we
confirm that positional errors are almost all smaller than $1\%$ of the embryo
length.

Beyond measuring the variance, we can estimate the distribution of positional
errors. Since the different stripes have slightly different $\sigma_{x}$, we
normalize the positional errors for each stripe individually,

$z_{\rm n}=(x_{\rm n}-\bar{x}_{\rm n})/\sigma_{x}({\rm n}).$ (22)

With this normalization we can pool across all 21 stripes, and we estimate the
distribution of $z$ as usual by making bins and counting the number of
examples in each bin, with results shown at left in Fig. A2. Qualitatively the
distribution is close to being Gaussian, but what matters for our analysis is
the entropy of this distribution.

When we estimate a probability distribution and use this estimate to compute
the entropy, the random errors in the distribution that arise from the
finiteness of our sample become systematic errors in the entropy. The general
version of this problem goes back to the very first efforts to use information
theoretic concepts to analyze biological experiments miller_55 ; for a review
see Appendix A.8 of Ref bialek_12 . Briefly, naive entropy estimates depend
systematically on the size of the sample, and if we can detect this systematic
dependence we can extrapolate to infinite data. At right in Fig. A2 we show
the difference between the entropy of the estimated distribution $P(z)$ and
the entropy of a Gaussian. We see that when we base our estimates on $N_{\rm
em}$ embryos there is a term $\sim 1/N_{\rm em}$. Extrapolating $N_{\rm
em}\rightarrow\infty$ we see that the entropy difference goes to zero within
the small ($<0.01\,{\rm bit}$) error bars. We conclude that, for the purposes
of our discussion, it is safe to approximate the positional errors as being
Gaussian.

Finally we can use the same extrapolation methods to provide a better estimate
of the “information gap” defined in the main text. Equation (16) defines the
positional information contained in the local signals, $I_{\rm position}$, and
the information gap is the difference between this and $I_{\rm
unique}=\log_{2}N$. Fig. A2 shows the values of

$I_{\rm gap}=I_{\rm unique}-I_{\rm
position}={\bigg{\langle}}\log_{2}\left[\sqrt{2\pi e}{{N\sigma_{x}({\rm
n})}\over L}\right]{\bigg{\rangle}}_{\rm n}$ (23)

estimated from fractions of our data set and then extrapolated. The result is
$I_{\rm gap}=1.39\pm 0.08\,{\rm bits}$ (Fig. A2).

Figure A2: (A) Positional errors are well approximated as Gaussian. An
estimate of the distribution of normalized errors, Eq (22). Open circles are
means pooled across all stripes and embryos; error bars are standard
deviations across random halves of the embryos; and the line is the Gaussian
with zero mean and unit variance. (B) The entropy difference between this
estimated distribution and the Gaussian, as a function of the (inverse) number
of embryos we include in our analysis. Points (cyan) are examples from random
choices out of the full ensemble of embryos; open circles with error bars are
the mean and standard deviations of these points; and the line is a linear
extrapolation miller_55 ; panzeri+treves_95 ; strong+al_98 ; paninski_03 ;
bialek_12 . (C) Estimates of the information gap, Eq (23). Points (cyan) are
examples from random choices out of the full ensemble of embryos; open circles
(blue) with error bars are the mean and standard deviations of these points;
and the line is a linear extrapolation to $I_{\rm gap}=1.39\pm 0.08\,{\rm
bits}$.

##  Appendix B Entropy estimates

Fig. A3 shows estimates of the extra information $\Delta S/N$ [Eq (9)] based
on measurements in different numbers of embryos, for $N=10$ and $N=20$
contiguous pair rule stripes. We see the expected dependence on $1/N_{\rm
em}$, and the steepness of this dependence is twice as large at $N=20$ than at
$N=10$. This gives us confidence in the extrapolation $N_{\rm
em}\rightarrow\infty$ miller_55 ; panzeri+treves_95 ; strong+al_98 ;
paninski_03 ; bialek_12 .

Figure A3: Entropy reduction by correlations among the pair rule stripe
fluctuations, estimated from different numbers of embryos $N_{\rm em}$; $N=10$
stripes at left and $N=20$ stripes at right. Points (cyan) are examples from
random choices out of the full ensemble of embryos; open circles (blue) with
error bars are the mean and standard deviations of these points; and the line
is a linear extrapolation to the square.

##  Appendix C Origin of the correlations

The precision of pair rule stripe placement matches, quantitatively, the noise
in optimal estimates of position based on the local expression levels of the
gap genes dubuis+al_13 ; petkova+al_19 . To be consistent with this result,
the correlations should also be visible in the gap genes. As noted above, Lott
and colleagues saw correlations in expression boundaries for selected gap
genes lott+al_07 , and later measurements showed that combinations of gap gene
expression levels have correlations extending over a significant fraction of
the embryo krotov+al_14 . Here we revisit these measurements and connect
fluctuations in gap gene expression to positional noise. Notice that for the
pair rule genes we can work directly with the positions of the stripes, but
for the gap genes we have to think more carefully about how positions are
encoded in expression levels.

We start with a brief review of ideas about decoding positional information
petkova+al_19 . Measurements of gap gene expression in multiple embryos
provide samples from the conditional distribution $P(\\{g_{\rm i}\\}|x)$, at
all values of the position $x$ along the anterior–posterior axis; we focus on
the $d=4$ gap genes expressed in the middle $\sim 80\%$ of the embryo,
hunchback, giant, krüppel, and knirps. To a good approximation this
distribution is Gaussian,

$\displaystyle P(\\{g_{\rm i}\\}|x)$ $\displaystyle=$
$\displaystyle{1\over{Z(x)}}\exp\left[-{1\over 2}\chi^{2}\left(\\{g_{\rm
i}\\};x\right)\right]$ (24) $\displaystyle Z(x)$ $\displaystyle=$
$\displaystyle\left[(2\pi)^{d}\det\hat{C}(x)\right]^{1/2}$ (25)
$\displaystyle\chi^{2}\left(\\{g_{\rm i}\\};x\right)$ $\displaystyle=$
$\displaystyle\sum_{{\rm i,j}=1}^{d}\left[g_{\rm i}-\bar{g}_{\rm
i}(x)\right]\left[\hat{C}^{-1}(x)\right]_{\rm ij}\left[g_{\rm j}-\bar{g}_{\rm
j}(x)\right],$

where $\bar{g}_{\rm i}(x)$ is the mean expression level of gene $\rm i$ at
position $x$ and

$\left[\hat{C}(x)\right]_{\rm ij}=\langle\delta g_{\rm i}\delta g_{\rm
j}\rangle_{x}$ (27)

is the covariance matrix of fluctuations around these means. To decode the
position of a cell from the local expression levels we need to construct

$P\left(x|\\{g_{\rm i}\\}\right)={{P(\\{g_{\rm i}\\}|x)P(x)}\over{P(\\{g_{\rm
i}\\})}}.$ (28)

But because nuclei are arrayed uniformly along the length of the embryo,
$P(x)$ is uniform and hence the dependence on $x$ is captured in Eq (24).

Figure A4: Decoding gap gene expression levels in a single embryo and
correlations in the resulting pattern of positional errors. (A) Expression of
Hb (blue), Kr (green), Gt (cyan), and Kni (red). Thin solid lines are means
across $N_{\rm em}=38$ embryos in a small window $40\leq t\leq 44\,{\rm min}$
in nuclear cycle 14; dense points are data from a single embryo petkova+al_19
. (B) Positional errors computed from Eq (C). (C) Correlations in the
positional noise inferred from gap gene expression. For each embryo $\alpha$
we compute the correlation function in Eq (35) and then normalize to give
$\tilde{C}(\Delta x)=C(\Delta x)/C(0)$. Blue circles with error bars are mean
and standard error across $N_{\rm em}=38$ embryos; solid red line is a smooth
curve to guide the eye.

A cell at the actual position $x_{\rm true}$ has expression levels

$g_{\rm i}=\bar{g}_{\rm i}(x_{\rm true})+\delta g_{\rm i},$ (29)

and if the positional noise is small we can write

$\bar{g}_{\rm i}(x)=\bar{g}_{\rm i}(x_{\rm true})+(x-x_{\rm
true}){{d\bar{g}_{\rm i}(x)}\over{dx}}{\bigg{|}}_{x=x_{\rm true}}+\cdots.$
(30)

If the noise is small, then the the best estimate of position based on the gap
gene expression levels is the value of $x$ which minimizes $\chi^{2}$, and
this can be written as

$\displaystyle\hat{x}$ $\displaystyle=$ $\displaystyle x_{\rm true}+\delta x$
(31) $\displaystyle\delta x(x_{\rm true})$ $\displaystyle=$
$\displaystyle\left[\sigma_{x}^{2}(x)\sum_{{\rm i,j}=1}^{d}\delta g_{\rm
i}\left[\hat{C}^{-1}(x)\right]_{\rm ij}{{d\bar{g}_{\rm
j}(x)}\over{dx}}\right]_{x=x_{\rm true}},$

where the variance of positional noise is defined by

${1\over{\sigma_{x}^{2}(x)}}=\sum_{{\rm i,j}=1}^{d}{{d\bar{g}_{\rm
i}(x)}\over{dx}}\left[\hat{C}^{-1}(x)\right]_{\rm ij}{{d\bar{g}_{\rm
j}(x)}\over{dx}};$ (33)

for consistency we have

$\langle[\delta x(x)]^{2}\rangle=\sigma_{x}^{2}(x).$ (34)

Note the connection to Eqs (1) and (2) in §II.

Previous work has emphasized the scale of positional errors $\sigma_{x}$
dubuis+al_13 ; tkacik+al_15 ; petkova+al_19 . But the optimal decoding of gap
gene expression levels petkova+al_19  maps the deviation of expression levels
from the mean into a decoding error for each embryo individually, as in Eq
(C). An example is in Fig. A4, where the small fluctuations of expression
levels around the mean (A) translate into proportionally small errors $\delta
x$ (B).

For each embryo $\alpha$ we can take the positional errors $\delta
x_{\alpha}(x)$ and compute the correlation function

$C_{\alpha}(\Delta x)={1\over{L-\Delta x}}\int dx\,\delta x_{\alpha}(x)\delta
x_{\alpha}(x+\Delta x).$ (35)

Fig. A4C shows the mean and standard error of the normalized correlation
function across all $N_{\rm em}=38$ embryos in our experimental ensemble.
Qualitatively, correlations in the positional noise encoded by the gap genes
extend over distances similar to the correlation in positional noise of the
pair rule stripes (Fig. 4). Quantitatively, the gap gene correlations include
an additional component with a short correlation length. One possibility is
that this component is averaged away by interactions among neighboring cells
during expression of the pair rule stripes. Another possibility is that a
modest fraction of the noise in gap gene expression reflects local noise in
the measurements, as discussed previously dubuis+al_13b ; this measurement
noise has only a small impact on our estimates of the effective noise
$\sigma_{x}$ but a larger impact on the shape of the correlation function. It
seems likely that both effects contribute. Nonetheless, it is clear that
relatively long ranged correlations, which are crucial to closing the
information gap, are present already in the gap gene expression levels, as
suggested in earlier work lott+al_07 ; dubuis+al_13 ; krotov+al_14 . New
experiments will be needed to give a reliable estimate of the information that
is encoded in these correlations.

## References

  * (1) AM Turing, The chemical basis of morphogenesis. Philos Trans R Soc Lond B, Biol Sci 237, 37–71 (1952).
  * (2) L Wolpert, Positional information and the spatial pattern of cellular differentiation. J Theor Biol 25, 1–47 (1969).
  * (3) G Tkačik and T Gregor, The many bits of positional information. Development 148, dev176065 (2021).
  * (4) C Nusslein–Vollhard and E Wieschaus, Mutations affecting segment number and polarity in Drosophila. Nature 287, 795–801 (1980).
  * (5) PA Lawrence. The Making of a Fly: The Genetics of Animal Design (Blackwell Scientific, Oxford, 1992).
  * (6) E Wieschaus and C Nüsslein–Volhard, The Heidelberg screen for pattern mutants of Drosophila: A personal account. Annu Rev Cell Dev Biol 32, 1–46 (2016).
  * (7) R Rivera–Pomar and H Jackle, From gradients to stripes in Drosophila embryogenesis: Filling in the gaps. Trends Genet 12, 478–483 (1996).
  * (8) J Jaeger, The gap gene network. Cell Mol Life Sci 68, 243–274 (2011).
  * (9) JP Gergen, D Coulter, and EF Wieschaus, Segmental pattern and blastoderm cell identities. In Gametogenesis and The Early Embryo, JG Gall, ed, pp 195–220 (Liss, New York, 1986).
  * (10) See https://bugguide.net/node/view/1308194/bgimage.
  * (11) JO Dubuis, G Tkačik, EF Wieschaus, T Gregor, and W Bialek, Positional information, in bits. Proc Natl Acad Sci (USA) 110, 16301–16308 (2013).
  * (12) F Liu, AH Morrsison, and T Gregor, Dynamic interpretation of maternal inputs by the Drosophila segmentation gene network. Proc Natl Acad Sci (USA) 110, 6724–6729 (2013).
  * (13) MD Petkova, G Tkačik, W Bialek, EF Wieschaus, and T Gregor, Optimal decoding of cellular identities in a genetic network. Cell 176, 844–855 (2019).
  * (14) J Yu, J Xiao, X Ren, K Lao, and XS Xie, Probing gene expression in live cells, one protein molecule at a time. Science 311, 1600–1603 (2006).
  * (15) Y Taniguchi, PJ Choi, G–W Li, H Chen, M Babu, J Hearn, A Emili, and XS Xie, Quantifying E. coli proteome and transcriptome with single–molecule sensitivity in single cells. Science 329, 533–538 (2010).
  * (16) JO Dubuis, R Samanta, and T Gregor, Accurate measurements of dynamics and reproducibility in small genetic networks. Mol Sys Biol 9, 639 (2013).
  * (17) SE Lott, M Kreitman, A Palsson, E Alekseeva, and MZ Ludwig, Canalization of segmentation and its evolution in Drosophila. Proc Natl Acad Sci (USA) 104, 10926–10931 (2007).
  * (18) D Krotov, JO Dubuis, T Gregor, and W Bialek, Morphogenesis at criticality. Proc Natl Acad Sci (USA) 111, 3683–3688 (2014).
  * (19) C Nüsslein-Volhard and S Roth, Axis determination in insect embryos. Ciba Found Symp 144, 37–55 (1989).
  * (20) C Nüsslein-Volhard, Determination of the embryonic axes of Drosophila. Dev Suppl 1, 1-10 (1991).
  * (21) G Tkačik, JO Dubuis, MD Petkova, and T Gregor, Positional information, positional error, and read–out precision in morphogenesis: a mathematical framework. Genetics 199, 39–59 (2015).
  * (22) CE Shannon, A mathematical theory of communication. Bell Sys Tech J 27, 379–423 and 623–656 (1948).
  * (23) W Bialek, Biophysics: Searching for Principles (Princeton University Press, Princeton NJ, 2012).
  * (24) T Gregor, DW Tank, EF Wieschaus, and W Bialek, Probing the limits to positional information. Cell 130, 153–164 (2007).
  * (25) V Antonetti, W Bialek, T Gregor, G Muhaxheri, M Petkova, and M Scheeler, Precise spatial scaling in the early fly embryo. arXiv:1812.11384 [q–bio.MN] (2018).
  * (26) GA Miller, Note on the bias of information estimates. In Information Theory in Psychology: Problems and Methods II-B, H Quastler, ed., pp. 95–100 (Free Press, Glencoe IL, 1955).
  * (27) S Panzeri and A Treves, The upward bias in measures of information derived from limited data samples. Neural Comp 7, 399–407 (1995).
  * (28) SP Strong, R Koberle, RR de Ruyter van Steveninck, and W Bialek, Entropy and information in neural spike trains. Phys Rev Lett 80, 197–200 (1998).
  * (29) L Paninski, Estimation of entropy and mutual information. Neural Comp 15, 1191–1253 (2003).
  * (30) AM Arias and P Hayward, Filtering transcriptional noise during development: concepts and mechanisms. Nat Rev Genet 7, 34–44 (2006).
  * (31) TC Lacalli, Patterning, from conifers to consciousness: Turing’s theory and order from fluctuations. Front Cell Dev Biol 10, 871950 (2022).
  * (32) CH Waddington, The Strategy of the Genes. A Discussion of Some Aspects of Theoretical Biology. (Allen and Unwin, London, 1957).
  * (33) MD Petkova, SC Little, F Liu, and T Gregor, Maternal origins of developmental reproducibility. Curr Biol 24, 1283–1288 (2014).
  * (34) PL Ferree, VE Deneke, and S Di Talia, Measuring time during early embryonic development. Semin Cell Dev Biol 55, 80–88 (2016).
  * (35) G Tkačik, AM Walczak, and W Bialek, Optimizing information flow in small genetic networks. Phys Rev E 80, 031920 (2009).
  * (36) AM Walczak, G Tkačik, and W Bialek, Optimizing information flow in small genetic networks. II: Feed–forward interaction. Phys Rev E 81, 041905 (2010).
  * (37) TR Sokolowski, T Gregor, W Bialek, and G Tkačik, Deriving a genetic regulatory network from an optimization principle. arXiv:2302.05680 [physics.bio–ph] (2023).
  * (38) T Doan, A Mendez, PB Detwiler, J Chen, and F Rieke, Multiple phosphorylation sites confer reproducibility of the rod’s single-photon responses. Science 313, 530–533 (2006).


# Reshaping and Enzymatic Activity allow Viruses to move through the Mucus

Falko Ziebert<EMAIL_ADDRESS>Institute for Theoretical
Physics, Heidelberg University, Philosophenweg 19, 69120 Heidelberg, Germany
and
BioQuant, Heidelberg University, Im Neuenheimer Feld 267, 69120 Heidelberg,
Germany  Kenan G. Dokonon  Department of Chemical Engineering and
Biotechnology, University of Cambridge, Cambridge CB3 0AS, UK.  Igor M. Kulić
<EMAIL_ADDRESS>Institut Charles Sadron UPR22-CNRS, 67034 Strasbourg, France
and
Institute Theory of Polymers, Leibniz-Institute of Polymer Research, D-01069
Dresden, Germany

###### Abstract

Filamentous viruses like influenza and torovirus often display systematic
bends and arcs of mysterious physical origin. We propose that such viruses
undergo an instability from a cylindrically symmetric to a toroidally curved
state. This “toro-elastic” state emerges via a spontaneous symmetry breaking
under prestress, induced via short range spike protein interactions and
magnified by the filament’s surface topography. Once surface stresses become
sufficiently large, the filament buckles and the toroidal, curved state
constitutes a soft mode that can propagate through the filament’s material
frame around a ”mexican-hat” potential. In the mucus of our airways, glycan
chains are omnipresent that influenza’s spike proteins can bind to and cut. We
show that when coupled to such a non-equilibrium chemical reaction, the curved
toro-elastic state can attain a spontaneous rotation for sufficiently strong
enzymatic activity, leading to a whole body reshaping propulsion similar to –
but different from – eukaryotic flagella and spirochetes.

Thin fibers are a common structural element in biology, from biofilaments to
the body shapes of viruses and bacteria. Notably, biological filaments like
the bacterial flagellar hook [1, 2], possibly microtubules [3, 4], and
intermediate filaments [5] exhibit a highly unusual, common motif: their
straight ground state can become unstable, reshaping them into polymorphic
toroids and superhelices. Not less surprisingly, when driven out of
equilibrium, stimuli-responsive but otherwise straight fibers can acquire a
dissipative toroidal steady state of spontaneous rotation [7, 6]. Visually
inspecting the shapes of filamentous viruses like Influenza C (IV-C) [8, 9] or
Torovirus [10], see Fig. 1, the natural question arises if viral envelopes
could posses a similar symmetry-broken, toroidal ground state. And as viruses
like influenza A (IV-A) and C [11, 12, 13], as well as paramyxovirus [14] have
been shown to be able to roll on surfaces due to their spike proteins’
catalytic activity, can this toroidal state be driven catalytically and used
for virus propulsion?

In the following we consider a plausible model for the emergence of what we
call the toro-elastic state of the viral envelope. It relies on surface
switchabilty at the interfaces of spike proteins and a mismatch-stress with
layers further inside the virus. The switch corresponds to the making-and-
breaking of any reversible physical bond – including hydrogen bonds, screened-
electrostatic, van der Waals bonds – as well as protein-tail-mediated
interactions. Notably, recently resolved surface structures of influenza C
[15] have shown extensive spike-spike contacts, establishing a hexagonal
lattice on the envelope. The natural cylindrical curvature of this spike-spike
network is further modulated by the internal transmembrane and nucleoprotein
matrix proteins, indicating some form of core-shell strain mismatch for
influenza C.

Figure 1:  Emergent toroidal shapes of various biofilaments: a) torovirus
[18]; b) the bacterial flagellar hook [1]; c) Influenza C: long filamentous
viruses show both straight and oscillatory shapes [9]; inset: short, curved
filamentous IV-C actively rolling on glycan coated surfaces [12].

After demonstrating the mechanism giving rise to the toro-elastic state, we
explain how its continuous ground state looks much like a built-in “wheel
within the material” [7], mirroring the “universal joint” model of the
bacterial flagellar hook proposed by Namba and coworkers [2]. Finally, we
study the effect of subtle dynamic force imbalances, generated by polymeric
ligands being bound and cut on the viral envelope. This situation reflects
what influenza is confronted within the mucus of our airways, which is a
highly viscous environment containing polymers with sugary ends (glycans) the
virus interacts with. We show that a broken detailed-balance via enzymatic
cutting can force the toroidally curved state to spontaneously rotate once the
enzyme surpasses a certain minimal activity threshold. This dynamic shape
rotation mode should allow the virus to move at small but noticeable speed,
and could explain why catalytic activity increases mucus penetration so
strongly for influenza [16, 17].

Toro-elastic state. Consider the virus surface covered with radial spikes – in
case of influenza C called Hemagglutinin-esterase-fusion protein (HEF) – with
initial spacing $l_{0}$ (along the filament axis, see Fig. 2a). If the virus
curves, variations of strains $\varepsilon$ and thus distances between the
spikes are amplified proportional to the effective radius of the cylindrical
virus, $R$, and the centerline curvature $\kappa$ like
$l(\kappa)\approx\varepsilon l_{0}=\kappa Rl_{0}$. In addition, there is a
geometric “shortcut factor” for finite deformations: For spiky surface
elements that interact along the shortest spatial distance, this distance is
given by the secant line (see SI for details).

Introducing the vectorial curvature perpendicular to the axis
$\vec{\kappa}=(\kappa_{x},\kappa_{y})=\kappa(\cos\phi,\sin\phi)$ with
orientation angle $\phi$ and magnitude $\kappa$, and parameterizing the
material-fixed spike positions with angle $\theta$ (with $\theta=\phi$
corresponding to the orientation of the curvature vector, i.e. towards the
bend), the spacing between surface contacts is given by, see Fig. 2b,c and SI,

$l(\kappa,\theta-\phi)=l_{0}\frac{\sin(\kappa l_{0}/2)}{\kappa
l_{0}/2}(1-R\kappa\cos\left(\theta-\phi\right))\,.$ (1)

Here the first term is the shortcut factor described above and the second
shows the dependence on the angle differences. At this point we can fix
$\vec{\kappa}$ to point in the $x$-direction, setting $\phi=0$.

Figure 2:  Geometry of the virus and the toroidal mode: (a) spike proteins are
arranged on a typically staggered column lattice with spacing $\l_{0}$ and
effective radius $R$ (including the spikes). (b) Curving the lattice by
curvature $\kappa$ allows some spikes to interact at the inside of the bend.
(c) Which particular spikes can interact depends on their orientation angle
$\theta$ with respect to the angle $\phi$ of the curvature vector
$\vec{\kappa}$.

Every spike can interact with its axial neighboring spikes, in case they are
within range, letting them gain a binding free energy (per unit length),
$v(l)$. On the flip side, binding along its long axis implies a curving of the
cylindrical virus, which costs elastic bending energy. Modeling the latter
contribution as a harmonic term with bending stiffness $B$, the total energy
per cross section is hence given by

$E=E_{bend}+E_{bind}=\frac{1}{2}B\kappa^{2}+\int_{0}^{2\pi}v\left(l(\kappa,\theta)\right)d\theta.$
(2)

Apart from its magnitude and short ranged nature, the detailed form of $v$
turns out to be non-essential for the phenomenology and emergence of the toro-
elastic state. An instructive potential is the linear-exponential one,

$v(l)=-v_{0}\frac{l}{\sigma}\exp\left(1-\frac{l}{\sigma}\right),$ (3)

as sketched in Fig. 3a.

Figure 3:  (a) A spike interaction potential with a minimum at a finite
distance $\sigma$. (b) The straight ground state can become unstable through
spike-spike interactions and give rise to a Mexican hat-potential and a “toro-
elastic” mode: in angular direction, it does not cost energy to change the
orientation of the spontaneous curvature (similar to usual buckling).
(c)“Phase diagram” in the plane $v_{0}$ vs. $\sigma$. In the elastic state,
$B_{eff}>0$ and the energy has its minimum at $\kappa=0$. In the toro-elastic
state, $B_{eff}<0$ and an energy minimum at finite $\kappa$ emerges. In the
coexistence region $\kappa=0$ re-emerges as a local minimum, and the toro-
elastic state is separated from it by a barrier. Parameters: $\sigma$ in nm
and $v_{0}$ in $\frac{k_{B}T}{\rm nm}$ for fixed $l_{0}=25$ nm. Virus radius
$R=50$ nm; bending rigidity $B=3\,k_{B}T\,\mu{\rm m}$.

It has a minimum, corresponding to the bound state, at $l=\sigma$ with
$v(\sigma)=-v_{0}$ and quickly flattens for $l>\sigma$. In the SI we also
discuss a Gaussian potential, which has the width of the well as an additional
adjustable parameter. While the latter makes the scenario more complex, it
still exhibits the same states as now discussed for the potential given by Eq.
(3).

Inserting Eqs. (1) and (3) into Eq. (2), the energy can be calculated
analytically [19]. Importantly, depending on parameters, it can have the shape
of a Mexican hat, see Fig. 3b: this is the “toro-elastic” state, where a
finite curvature is preferred whose direction (i.e angle $\phi$) is arbitrary
and hence corresponds to a soft mode. Fig. 3c shows the different energy
landscapes/shapes when varying the depth, $v_{0}$, and the minimum position,
$\sigma$, of the potential for fixed spike distance $l_{0}$ and virus
stiffness $B$. This phase diagram displays: (i) an elastic state, where the
stable state is a single point $\overrightarrow{\kappa}=0$ corresponding to a
straight virus; (ii) the toro-elastic region, where the stable state space is
a circle with $\left|\overrightarrow{\kappa}\right|=\kappa_{tor}$; and (iii) a
coexistence region, where the (locally) stable state space is a circle and the
origin, with a barrier to reach the toro-elastic state.

To understand when the straight state, $\kappa=0$, becomes unstable we expand
the energy up to second order in curvature, yielding
$E=E_{0}+\frac{1}{2}B_{eff}\kappa^{2}$ with $E_{0}=-\frac{2\pi
l_{0}v_{0}}{\sigma}e^{1-\frac{l_{0}}{\sigma}}$ and an effective bending
stiffness $B_{eff}=B+B_{bind}$ with

$\frac{B_{bind}}{2\pi}=v_{0}l_{0}^{2}e^{1-l_{0}/\sigma}\left[\frac{R^{2}}{\sigma^{2}}\left(1-\frac{l_{0}}{2\sigma}\right)+\frac{l_{0}}{12\sigma}\left(1-\frac{l_{0}}{\sigma}\right)\right].\,\,$
(4)

This binding-induced part of the effective stiffness can become negative. In
fact, this condition can be always satisfied for $l_{0}\gtrapprox 2\sigma$ if
in addition the interaction strength is sufficiently large, $v_{0}\gg
e^{l_{0}/\sigma}B/l_{0}^{2}$: then $B_{bind}$ overrules $B$ and the filament
curves with higher order terms stabilizing a finite toroidal curvature
$\kappa_{tor}$. This toro-elastic state generically appears for strong surface
interaction and sufficiently large (but not too large) spike spacing vs.
interaction distance, cf. also the study of the more general Gaussian
potential in the SI.

The curvature of the toro-elastic state adjusts such that a significant
portion of the cross-section (typically half of it), benefits from the surface
interactions. For that, the elements of the surface need to shift from their
unperturbed distance $l_{0}$ to $\sigma$, i.e. by an amount
$l_{0}-\sigma\sim\kappa_{tor}Rl_{0}$, which yields the following estimate for
the curvature of the toro-elastic state

$\kappa_{tor}\sim R^{-1}\left(1-\sigma/l_{0}\right).$ (5)

The stiffness of the toroelastic state (i.e. the curvature in the rim of the
Mexican hat, cf. Fig. 3b, along the radial direction), is dominated by the
surface interactions. Assuming a curvature variation $\delta\kappa$ around
$\kappa_{tor}$, such that the binding distance changes by $\sim\sigma$, we
have $\delta\kappa\sim\frac{\sigma}{Rl_{0}}$ and the binding energy varies by
$\delta E\approx|E_{0}|$. Equating
$\frac{1}{2}K_{tor}\left(\delta\kappa\right)^{2}\approx|E_{0}|$ then allows to
estimate the stiffness of the toroelastic state to be

$K_{tor}\approx
4\pi\left(\frac{l_{0}}{\sigma}\right)^{3}e^{1-\frac{l_{0}}{\sigma}}R^{2}v_{0}$
(6)

up to subdominant contributions from the bending elasticity.

Considering numbers, typical length scales are $l_{0}=25$ nm, $\sigma=10$ nm
and virus radius $R=50$ nm; a reasonable estimate for the contact potential
strength is $v_{0}=1\frac{k_{B}T}{\rm nm}$. The bending rigidity of a virus is
hard to measure [20], but considering the virus as a membrane tube we estimate
$B=2\pi R\cdot 10k_{B}T\simeq 3\,k_{B}T\,\mu{\rm m}$. Using these values, one
finds $B_{bind}\simeq-5.7\,k_{B}T\,\mu{\rm m}$, hence indeed $B_{eff}<0$ is at
reach. Further we get $\kappa_{tor}\simeq(100\,{\rm nm})^{-1}$ estimated from
Eq. (5); the toro-elastic stiffness, $K_{tor}\simeq 100\,k_{B}T\,\mu{\rm m}$,
is much larger than the usual bending rigidity.

Toroelastic state driven by ligand dynamics. Having seen how the toro-elastic
state arises, we now show how it can be driven in an out-of-equilibrium
situation. We assume that the virus is embedded in a solvent that contains a
substantial concentration of glycan coated polymers (as in the mucus), which
we model here for simplicity as di-glycan ”dumbbells”. If these glycans attach
to two neighboring spikes, they introduce subtle force imbalances on top of
the toro-elastic state. Moreover, these imbalances are dynamic due to the
enzymatic cutting ability of the HEF spikes.

We assume the di-glycans to be harmonic polymer springs with spring constant
$S$ and with a vanishing preferred length. They can be unbound (state
$B_{0}$), bound to one spike (state $B_{1}$) or to two neighboring spikes
(along the backbone), defining the double bound state $D$, see the upper part
of the reaction scheme in Fig. 4a. In this double bound state, the glycans
stretch out to a non-zero actual length as already given by Eq. (1), with
$\phi$ the orientation of the curvature vector $\overrightarrow{\kappa}$ and
$\theta$ the angular position of the spike it attaches to in the considered
cross-section. In general, a cross-section will have an azimuthal imbalance of
the (double-) bound density $D\left(\theta\right)$. The total glycan
stretching energy per cross-section for a virus segment of length $l_{0}$ then
reads

$E_{str}\left(D\right)=\frac{S}{2}\int
D\left(\theta\right)l^{2}(\kappa,\theta).$ (7)

Assuming small curvatures $\kappa R,\,\kappa l_{0}\ll 1$ and applying a mode
ansatz $D(\theta)=D_{0}+D_{s}\sin\theta+D_{c}\cos\theta$, to linear order
$E_{str}\left(D\right)\approx\frac{S}{2}l_{0}^{2}\left(D_{0}-R\kappa\left(D_{c}\cos\phi+D_{s}\sin\phi\right)\right)$.
The corresponding torque (per $l_{0}$ section) is then given by
$M_{dr}=-\frac{\partial E_{str}}{\partial\phi}$ and in the frame co-moving
with the toro-elastic curvature evaluates to

$M_{dr}=\frac{Sl_{0}^{2}}{2}R\kappa D_{s}=m_{dr}D_{s}.$ (8)

This is the driving torque from the inhomogeneously distributed double bound
glycans, with $m_{dr}$ the characteristic torque scale. With $S\sim
0.01-0.1\frac{k_{B}T}{{\rm nm}^{2}}$, the above parameters and
$\kappa=\kappa_{tor}$, $m_{dr}$ can be estimated to be $2-20\,k_{B}T$, which
is appreciable [21].

Figure 4:  (a) (left) Interaction of tension-inducing polymeric ligands
(dumbbell springs) with the virus spikes. (right) The considered reaction
kinetics of equilibrium binding and catalytic cutting. (b) The toroidal mode’s
angular frequency (in units of $\rm{s}^{-1})$ as a function of the cutting
rate $k_{cut}$ (units $\rm{s}^{-1})$. The red curve shows the analytical
result, Eq. (11), obtained without friction, and the black curve shows the
numerical result with realistic frictional torque coefficient
$\xi_{h}/m_{drive}=0.01$ (see SI for details). Parameters: glycan
concentration $G=3\,{\rm mM}$; $k_{1}=0.1{\rm mM}^{-1}{\rm s}^{-1}$,
$k_{-1}=0.1{\rm s}^{-1}$ (for IV-A [22, 23]); $k_{2}^{0}=1{\rm mM}^{-1}{\rm
s}^{-1}$, $k_{-2}=0.1{\rm s}^{-1}$ (estimated). (c) The uniform, shape
invariant rotation of a short toroidal virus or an elongated helical virus
with frequency $\omega$ leads to a swimming velocity $v$ with the kinematics
of fluid vortex lines [25]-[27].

Beyond neighboring spikes getting linked by freely floating di-glycans, we now
also consider the enzymatic cutting activity of the HEF spike. The initial
equilibrium reaction scheme has to be extended by two more states, as sketched
in the bottom part of Fig. 4a. First, a double bound glycan can be cut into
two, to a state $C_{2}$ that still blocks both spikes for further attachments.
If one of these cut parts detaches, or if a single bound glycan is cut, this
yields state $C_{1}$. Only if all cut parts have left, one recovers state
$B_{0}$ and the spike can be “used again” for attachment. The full reaction
scheme of Fig. 4a translates into

$\displaystyle\frac{d}{dt}B_{1}$ $\displaystyle=$ $\displaystyle
k_{1}GB_{0}-\left(k_{-1}+k_{2}\left(\theta-\phi\right)+k_{cut}\right)B_{1}+k_{-2}D,$
$\displaystyle\frac{d}{dt}D$ $\displaystyle=$ $\displaystyle
k_{2}\left(\theta-\phi\right)B_{1}-k_{-2}D-k_{cut}D,$
$\displaystyle\frac{d}{dt}C_{1}$ $\displaystyle=$ $\displaystyle
k_{cut}B_{1}+k_{-1}C_{2}-k_{-1}C_{1},$ $\displaystyle\frac{d}{dt}C_{2}$
$\displaystyle=$ $\displaystyle k_{cut}D-k_{-1}C_{2},$ (9)

with an additional equation for $B_{0}$. This latter can be eliminated via
$B_{0}+B_{1}+C_{1}+C_{2}+D=1$, since all quantities are site population
fractions. Note that all considered states are functions of $\theta$. Assuming
a curved virus that is rotating with constant angular frequency $\omega$, i.e.
$\phi(t)=\omega t$, and $\kappa=const.$ as suggested by the high stiffness of
the toroidal state estimated above, we can move into the co-moving frame,
transforming $\phi\rightarrow 0$ and
$\frac{d}{dt}\rightarrow\omega\frac{\partial}{\partial\theta}$.

Importantly, the double binding ”on” constant $k_{2}\left(\theta\right)$
depends on the curvature – since di-glycans bind more easily at the inside of
bent regions as they have to stretch more on the outside of the bend. For
small curvatures $\kappa R,\kappa l_{0}\ll 1$, one can relate the reaction
rates on the two opposing sides (see SI)

$\displaystyle\frac{k_{2}\left(\theta\right)}{k_{2}\left(\theta+\pi\right)}=e^{\frac{S}{2}\frac{-l^{2}\left(\theta\right)+l^{2}\left(\theta+\pi\right)}{k_{B}T}}\approx
1+\left(\frac{2Sl_{0}^{2}}{k_{B}T}R\kappa\right)\cos\theta$

yielding

$k_{2}\left(\theta\right)=k_{2}^{0}\left(1+\alpha\cos\theta\right)\,,\,\,\alpha=\frac{Sl_{0}^{2}}{k_{B}T}R\kappa$
(10)

where $\alpha$ is a dimensionless coupling constant, in the following assumed
to be small for simplicity.

The equations for the steady state dynamics in the co-moving coordinate system
can be solved by applying a mode ansatz, where as above for the driving
torque, all the chemical species can have a constant, a $\cos$\- and a
$\sin$-mode, $X=X^{0}+X^{c}\cos\theta+X^{s}\sin\theta$ for
$X=B_{1},D,C_{1},C_{2}$. In the case of negligible fluid friction, and for
fast binding of the second bond of the dumbbell, $k_{2}^{0}\gg k_{others}$,
the solution can be given analytically, see SI. One finds that $\omega=0$ is
always a solution, while solutions with finite angular velocity, given by

$\omega^{2}=\frac{k_{1}Gk_{cut}}{2}\left(1+\sqrt{1+\frac{8k_{-1}^{2}}{k_{1}Gk_{cut}}}\right)-k_{-1}^{2},$
(11)

bifurcate from this stationary branch for sufficiently rapid glycan cutting
rate

$k_{cut}>k_{cut}^{c,0}=\frac{k_{-1}^{2}}{3k_{1}G}.$ (12)

It is also possible to solve the problem in presence of friction. Assuming a
hydrodynamic (Stokes) rotational friction torque $M_{h}=\xi_{h}\omega=4\pi\eta
R^{2}l_{0}\omega$ and equating it with the driving torque, Eq. (8), one gets
$k_{cut}=k_{cut}^{c,0}+k_{cut}^{1}\left(\eta\right)$ with a correction that is
to leading order linear in fluid viscosity $\eta$, proportional to
$\frac{R}{\alpha^{2}Sl_{0}\kappa}$ and that in addition depends on the
chemical rates, see SI for details. Fig. 4b) shows the angular velocity
$\omega$ as a function of the cutting rate $k_{cut}$: the red curve shows (the
positive branch of) Eq. (11) and the black curve are numerical results
including the frictional torque. Since IV-C has not been characterized yet
concerning its kinetic parameters, we had to use estimates for IV-A [22, 23].
Note that cutting rates are estimated very conservatively, as IV-A’s
neuraminidase can reach cutting rates of up to $k_{cut}=15\,{\rm s}^{-1}$
[24].

Helix formation. If the virus spike lattice consisted of discrete uncoupled
sections, the toro-elastic energy would only depend on the magnitude of the
local curvature. Consequently, each cross-section would behave independently
with its curvature vector pointing in random directions. However, in general
one expects multiple interdigitated and mutually staggered sub-lattices of
spike-pairs, as sketched in Fig. 2a as pairs 1-2 and 3-4, with each of the
lattices experiencing a mean effective curvature over their length scale
($\sim l$). Assuming for simplicity two sub-lattices, I and II, the effective
curvature stemming from sub-lattice I and felt by sub-lattice II is
$\kappa_{II}=\left|\frac{1}{l}\int_{-l/2}^{+l/2}\overrightarrow{\kappa}_{I}\left(s\right)ds\right|$.
For a curvature vector of constant modulus that slowly winds with a torsion
rate $\tau\ll l^{-1}$ like
$\overrightarrow{\kappa}\left(s\right)=\kappa\left(\cos\left(\tau
s\right)\mathbf{e}_{x}+\sin\left(\tau s\right)\mathbf{e}_{y}\right)$, the
total energy of the two sub-lattices can be written as, cf. SI,

$\frac{E_{tot}\left(\kappa,\tau\right)}{K_{tor}}=\frac{1}{4}\left(\kappa_{tor}-\kappa\right)^{2}+\frac{1}{4}\left[\left(\kappa_{tor}-\kappa\right)+\kappa\frac{l^{2}\,\tau^{2}}{24}\right]^{2}.\quad$
(13)

The energy is positive definite, with a curved, untwisted ground state
$\kappa=\kappa_{tor}$ and $\tau=0$. The first term is similar to a classical
bending energy, followed by an unusual torsion-bend coupling that gives rise
to a quartic torsional energy for $\kappa=\kappa_{tor}$. Notably, for
curvatures larger than optimal, $\kappa>\kappa_{tor}$, the torsion becomes
bistable with two preferred states
$\tau_{1,2}=\pm\frac{\sqrt{24}}{l}\sqrt{1-\kappa_{tor}/\kappa}$. That is, we
expect the untwisted toroidal ground state to turn into a bistable helix (of
any handedness) upon additional external bending torques. This includes in
particular the torques induced by the above-discussed glycan binding. In
practice, the symmetry of the two helical states is likely to be broken by a
geometrical chirality of the natural spike lattice – an effect neglected here
for simplicity.

Self-propulsion. For a toroidal or helical object uniformly rotating its body
along the toroidal mode, while keeping its shape invariant (see Fig.4c), we
expect the swimming kinematics to be that of a line of fixed vorticity
$\propto\omega$, that is self-advecting through the fluid [25]. In particular
the swimming kinematics becomes in leading order independent of the (small)
helical torsion and occurs at typical propulsive velocities of spinning tori
$v\sim\omega\kappa R^{2}$ [26, 27]. Note that this shape invariant swimming
contrasts usual body-reshaping swimmers (propagating bending waves), where the
swimming speed depends on the undulation amplitude and wavelength but not on
the radius [29].

In conclusion, we have seen that viruses possessing the toro-elastic state can
attain a spontaneous rotation upon binding and cutting interactions with the
glycans floating in the mucus. The predicted angular velocities associated
with this dissipative rotary shape-flexing dynamics are comparable to the
rolling angular velocities of surface attached viruses (fractions of ${\rm
s}^{-1}$) [12, 13], where the glycan cutting and binding rates both contribute
to the speed. The typical active swimming velocity of tens of $\rm{nm}/s$ is
rather modest compared to swimming speeds of organisms – bacterial flagellae
turn at 50 Hz in water, flagella-mediated body-reshaping in spirochetes has a
similar frequency [29] – but is still of the order of one virus size per
second. Hence it is preferable to being effectively immobilized in the
extremely viscous, sticky environment. The proposed mechanism hence adds
another possible motility mode to the virus’ repertoire, beyond ratchet-like
motion [30] and rolling [11, 12].

###### Acknowledgements.

Acknowledgements. IMK thanks Jens-Uwe Sommer for continuous discussion and
support.

## References

  * [1] S. Shibata, H. Matsunami, S.-I. Aizawa and M. Wolf, Nat. Struct. Mol. Biol. 26, 941 (2019).
  * [2] T. Kato, F. Makino, M. Tomoko, P. Horvath and K. Namba, Nat. Commun. 10, 5295 (2019).
  * [3] H. Mohrbach, A. Johner and I.M. Kulić, Phys. Rev. Lett., 2010, 105, 268102.
  * [4] F. Ziebert, H. Mohrbach and I.M. Kulić, Phys. Rev. Lett., 2015, 114, 147101.
  * [5] L. Bouzar, M. M. Müller, R. Messina, B. Nöding, S. Köster, H. Mohrbach and I. M. Kulić, Phys. Rev. Lett. 122, 098101 (2019)
  * [6] A. Bazir, A. Baumann, F. Ziebert and I. M. Kulić, Soft Matter 16 5210 (2020).
  * [7] A. Baumann, A. Sánchez-Ferrer, L. Jacomine, P. Martinoty, V. Le Houerou, F. Ziebert and I. M. Kulić, Nat. Mater. 17, 523 (2018).
  * [8] M. D. Badham and J. S. Rossman, Curr. Clin. Microbiol. Rep. 3, 155 (2016).
  * [9] M. L. Martin, E. L. Palmer and A. P. Kendal, J. Clin. Microbiol. 6, 84 (1977).
  * [10] A. E. Hoet and M. C. Horzinek, Torovirus, in Encyclopedia of Virology, 3rd Edition, p. 151, 2008\.
  * [11] T. Sakai, S. I. Nishimura, T. Naito and M. Saito, Sci. Rep. 7, 45043 (2017).
  * [12] T. Sakai, H. Takagi, Y. Muraki and M. Saito, J. Virol. 92, e01522-17 (2018).
  * [13] F. Ziebert and I. M. Kulić, Phys. Rev. Lett. 126, 218101 (2021).
  * [14] X. Wu, M. Goebbels, L. Chao, T. Wennekes, F. J. M. van Kuppeveld, E. de Vries and C. A. M. de Haan, PLoS Pathog. 19, e1011273 (2023).
  * [15] S. Halldorsson, K. Sader, J. Turner, L. J. Calder and P. B. Rosenthal, Nat. Commun. 12, 1694 (2021).
  * [16] M. Cohen, X.-Q. Zhang, H. P. Senaati, H.-W. Chen, N. M. Varki, R. T. Schooley and P. Gagneux, Virol. J. 10, 321 (2013).
  * [17] X. Yang, L. Steukers, K. Forier, R. Xiong, K. Braeckmans, K. Van Reeth and H. Nauwynck, PLoS ONE 9, e110026 (2014).
  * [18] D. Rodriguez Aguirre, Virus Taxonomy – Classification and Nomenclature of Viruses. Online-Edition. Chp. “Coronaviridae”. In: ICTV Reports. International Committee on Taxonomy of Viruses, 2011; Adapted from from: https://en.wikipedia.org/wiki/Torovirus.
  * [19] The result is  $E=\frac{1}{2}B\kappa^{2}-2\pi v_{0}\frac{a}{\kappa R}e^{\left(1-\frac{a}{\kappa R}\right)}\left(I_{0}(a)-\kappa R\,I_{1}(a)\right),$ (14) where $a=\frac{2R}{\sigma}\sin(\kappa l_{0}/2)$ and $I_{0},I_{1}$ are the modified Bessel functions of the first kind.
  * [20] I. A. T. Schaap and F. Eghiaian and A. des Georges and C. Veigel, J. Biol. Chem. 287, 41078 (2012).
  * [21] Note that this value is large, considering that it acts on a short section of size $l_{0}$. As a reference, the rotational friction at an angular frequency $\omega=2\pi{\rm s}^{-1}$ in a rather viscous medium of $\eta=0.1\,{\rm Pa\,s}$ (100 times of water) for $R=50\,{\rm nm}$ is approximately $M\approx 0.05k_{B}T$, i.e. largely subdominant compared to the linker-induced torque $M_{dr}$ for not too low $D_{s}$.
  * [22] N. K. Sauter, M. D. Bednarski, B. A. Wurzburg, J. E. Hanson, G. M. Whitesides, et al., Biochem. 28, 8388 (1989).
  * [23] N. K. Sauter, J. E. Hanson, G. D. Glick, J. H. Brown, R. L. Crowther, et al., Biochem. 31, 9609 (1992).
  * [24] S. E. Adams, N. Lee, V. Y.Lugovtsev, A. Kan, R. P. Donnelly and N. A. Ilyushina, Antiviral Res. 169, 104539 (2019).
  * [25] C. Pozrikidis, Fluid Dynamics: Theory, Computation and Numerical Simulation Ch.11. Second Edition, 2009 Springer.
  * [26] A.T. Chwang and W.S. Hwang, Phys. Fluids 8, 1309 (1990).
  * [27] R.M. Thaokar, H. Schiessel and I.M Kulić, Eur. Phys. J. B 60, 325-336 (2007).
  * [28] Charon et al., Annu. Rev. Microbiol. 66, 349 (2012).
  * [29] D. K. Vig and C. W. Wolgemuth, Phys. Rev. Lett. 109, 218104 (2012).
  * [30] M. D. Vahey and D. A. Fletcher, eLife 8, e43764 (2019).