Stable developmental patterns of gene expression

without morphogen gradients

Maciej Majka1,2 , Nils B. Becker3,4 , Pieter Rein ten Wolde3 , Marcin Zagorski1 , and
Thomas R. Sokolowski3,5*
1
Institute of Theoretical Physics and Mark Kac Center for Complex Systems Research,
Jagiellonian University, ul. prof. Stanislawa Lojasiewicza 11, 30-348 Kraków, Poland
arXiv:2306.00537v3 [[Link]-ph] 27 Aug 2024

2
Present address: Department of Physics, East Carolina University, Greenville NC, 27858,
USA
3
FOM Institute AMOLF, Science Park 104, 1098 XG Amsterdam, The Netherlands
4
Present address: Theoretical Systems Biology, German Cancer Research Center, 69120
Heidelberg, Germany
5
Present address: Frankfurt Institute for Advanced Studies (FIAS), Ruth-Moufang-Straße
1, 60348 Frankfurt am Main, Germany

August 29, 2024

* Corresponding author, e-mail: sokolowski@[Link]

1
 Abstract
Gene expression patterns in developing organisms are established by groups of cross-regulating
target genes that are driven by morphogen gradients. As development progresses, morphogen activity
is reduced, leaving the emergent pattern without stabilizing positional cues and at risk of rapid
deterioration due to the inherently noisy biochemical processes at the cellular level. But remarkably,
gene expression patterns remain spatially stable and reproducible over long developmental time spans
in many biological systems. Here we combine spatial-stochastic simulations with an enhanced sampling
method (Non-Stationary Forward Flux Sampling) and a recently developed stability theory to address
how spatiotemporal integrity of a gene expression pattern is maintained in developing tissue lacking
morphogen gradients. Using a minimal embryo model consisting of spatially coupled biochemical
reactor volumes, we study a prototypical stripe pattern in which weak cross-repression between nearest
neighbor expression domains alternates with strong repression between next-nearest neighbor domains,
inspired by the gap gene system in the Drosophila embryo. We find that tuning of the weak repressive
interactions to an optimal level can prolong stability of the expression patterns by orders of magnitude,
enabling stable patterns over developmentally relevant times in the absence of morphogen gradients.
The optimal parameter regime found in simulations of the embryo model closely agrees with the
predictions of our coarse-grained stability theory. To elucidate the origin of stability, we analyze a
reduced phase space defined by two measures of pattern asymmetry. We find that in the optimal
regime, intact patterns are protected via restoring forces that counteract random perturbations and
give rise to a metastable basin. Together, our results demonstrate that metastable attractors can
emerge as a property of stochastic gene expression patterns even without system-wide positional
cues, provided that the gene regulatory interactions shaping the pattern are optimally tuned.

2
1 Introduction
Maintaining the integrity of spatial gene expression patterns over time is essential in embryonic devel-
opment. In early embryo development locally expressed morphogenetic proteins spread through the tis-
sue to form gradients of chemical signals [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]. Inside developing cells,
these chemical signals are interpreted by gene regulatory networks to form remarkably precise and repro-
ducible spatial patterns of gene expression that subsequently give rise to different body parts and organs
[12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23]. However, as spatial patterns are established by read-
ing out upstream morphogen gradients, their stability is constantly subject to inherently noisy cellular
and extracellular processes [24, 25, 26, 27, 28, 29, 30]. Moreover, the activity of morphogenetic gradi-
ents that is interpreted by target cells can decrease over developmental time. This decrease in activity
can take different forms, including reduction of the relative signaling range as the embryo grows in size
[16, 31, 32], signalling pathway desensitization [33], or complete disappearance of the gradients at later
developmental stages [34, 35]. Together, the inherent cellular stochasticity and reduced role of morphogen
gradients at later stages raise the question whether stable patterns can be maintained over sufficiently
long developmental times in the absence of morphogen gradients, and, if so, how.
Focusing on the cellular stochasticity, biological cells are facing two types of noise, namely intrinsic
and extrinsic noise, with different notions of robustness against the respective noise types [24, 25, 26,
27, 28, 29, 30]. Intrinsic noise originates from the processes of gene regulation, protein production, and
intracellular transport. Thus, robustness of spatial patterns to intrinsic noise amounts to buffering random
fluctuations in the copy numbers of patterning proteins. Extrinsic noise, on the other hand, terms the
variations originating from different external conditions including cell size variability [36, 37], cell-to-cell
variation in ribosome abundance [27] or fluctuations in the cellular environment [38, 39]. Therefore, the
robustness of spatial patterns to extrinsic noise refers to the capability of producing precise patterns in spite
of imperfect initial conditions, classically termed “canalization” in Waddington’s picture of development
[40, 41]. Several gene regulatory strategies providing either type of robustness have been studied [11, 21,
29], but our understanding of how nature orchestrates them in the fully interacting wild-type organism is
still incomplete.
Among the regulatory mechanisms that drive developmental pattern formation, the regulatory motif in
which two genes mutually repress each other is particularly prevalent [4, 42, 43, 44, 45, 15, 46, 47, 20, 30].
Intriguingly, mutual repression can have a dual role in the establishment of spatial patterns. On the one
hand, in systems driven by threshold-dependent activation of patterning genes via morphogen gradients,
mutual repression is crucial for shaping out expression domains that are bounded from two sides, thus
increasing the positional information carried by the expression pattern [14, 48, 16, 20, 22]. On the other
hand, mutual repression can induce bistability leading to stochastic switching between cell fates. Hence,
it is a priori unclear to which extent mutual repression supports or counter-acts the formation of stable
spatial patterns [14, 16, 20]. This issue is particularly relevant to systems that lack external cues for
symmetry breaking, such as morphogen or maternal gradients, that could force bistable cells into one of
their opposing fates.
Here we ask whether a system of mutually interacting genes can maintain an initially arranged expres-
sion pattern in the absence of upstream input gradients. To address this question we study a spatially
resolved gene regulatory network, conceptually motivated by the gap gene system in the early embryo
of the fruit fly Drosophila melanogaster [49, 50, 51, 52, 53, 54]. This system implements a particular
regulatory architecture, in which weak and strong mutual repressive interactions between expression do-
mains of different genes alternate depending on whether the domains are adjacent or not. This motif,
termed “alternating cushions”, was earlier investigated in terms of stability and robustness against extrinsic
noise in initial conditions [43]. That study employed a reaction-diffusion model with step-like activation
functions for representing the underlying gene expression dynamics. Using the so-called “moving kink
approximation”, the study predicted an extensive basin of pattern stability in the parameter space of the
model, where the stability could be attributed to repulsive forces between mutually repressing gene ex-
pression domains (“cushions”). More recently, an exact solution was obtained in an analogous model for
the dynamics of the contact zones between two gene expression domains and for arbitrary combination of
activating or repressing interactions between the involved genes [55]. This work provided exact conditions
for stability, leading to a better quantitative understanding of the conditions under which gene expression
patterns can survive arbitrary long time. Importantly, it was shown that perfect pattern stability (i.e., a
pattern surviving infinitely long) can only be achieved for a very specific combination of system parameters;
nevertheless, in the vicinity of these states, there exists a continuity of well-defined but slowly changing

3
gene expression patterns, which can fulfill their biological role for a finite but typically sufficiently long
period of time. However, since the reaction-diffusion model considered in [55] is only a continuous and
deterministic limit of the genuinely stochastic microscopic dynamics of gene expression, it remained unclear
whether the derived stability conditions provide useful insight into the regime of strong fluctuations.
In this work, we assess the temporal stability of gene expression patterns interacting via the “alternating
cushions” motif by numerical simulations of a minimal spatial-stochastic model that features a full micro-
scopic representation of stochastic gene expression and protein diffusion, thus incorporating the relevant
intrinsic noise sources. Using Non-Stationary Forward Flux Sampling (NS-FFS) [56, 57], an enhanced
biased sampling technique for stochastic systems changing in time, we quantify for how long patterns
shaped and maintained only by mutual repression can self-sustain. Contrasting with previous approaches
[43, 58, 49], NS-FFS allows us to go beyond a local, linear stability analysis of the studied system, and
to assess the depth and the width of the emerging basin of stability from large ensembles of stochastic
trajectories of the full spatially interacting gene expression pattern. Moreover, we derive the effective,
deterministic model of simulated system that expands the stability theory from [55] to the case of multiple
interfaces and allows us to determine the parameter regime within which the distances between boundaries
of adjacent gene expression domains are predicted to remain stable. Eventually, we employ this model to
identify the mechanism enhancing the pattern survival time.
Our results show that the stability of patterns arranged in the alternating cushions scheme strongly
varies with the strength of mutual repression between adjacent gene expression domains. We find that
pattern stability time is significantly longer when spatially adjacent genes repress each other with inter-
mediate strength and the next-nearest neighbor genes repress each other strongly. This results in a broad
peak of pattern survival time for a range of interaction strength ratios, with a single maximum at the
optimal choice. In this enhanced regime, we confirm the existence of robust restoring forces and find
signatures of a metastable basin that stabilizes well-ordered patterns (dynamical attractor), in accordance
with the previous findings of [43]. Away from the optimum, forces induced by strong nearest neighbor mu-
tual repression destroy the stripe patterns rapidly, while for weaker nearest neighbor repression the forces
are imperceptible when compared to stochastic fluctuations. We manage to explain these observations
employing our deterministic, effective model and the recent exact stability theory. We determine the the-
oretical optimal interaction strength ratio, situated in the vicinity of the numerically predicted optimum.
Further analysis reveals a nuanced interplay between fluctuations and a few stabilizing mechanisms present
in the deterministic, effective model, leading to enhanced survival time in the vicinity of optimal choice
and qualitatively in agreement with numerical observations. In result, we highlight the connection between
effective restoration forces seen in simulations, moving-kink approximation model [43] and exact stability
theory [55]. Going beyond the setting studied in [43], we also show that pinning of the pattern at the
embryo boundaries, which could be achieved by very short-ranging, peripherally acting maternal inputs,
can significantly further enhance the optimal pattern stability.
Taken together, we demonstrate that forces generated in the alternating cushions scheme can maintain
the gene expression pattern subject to stochastic production and diffusion of proteins for extremely long
times, thanks to the interplay between fluctuations and deterministic dynamics, constituting emergent
noise-control mechanism for the close-to-optimal choice of mutual repression parameters.

4
 A C A

B D

STANDARD NUCLEAR VOLUME NUCLEAR VOLUME AT BOUNDARY

A2 C2 A2 C2

D2 B2 C2 D2 B2 A2 D2 B2 A2
A C A C

B2 D2 B2 D2

A2 C2 D2 A2 C2 B2 A2 C2 D2 A2 C2 B2
B D B D

Figure 1: Schematic of the spatial gene-regulatory model. We use a cylindrical lattice of reaction
volumes to mimic the arrangement of cortical nuclei in the posterior Drosophila embryo at developmental
cycle 14. In each nuclear volume (shaded squares) we simulate production, degradation, dimerization and
mutual repression of the four genes A, B, C and D via the Gillespie algorithm. Each gene is subject to
repression by the protein dimers of the other genes, as indicated by the schematic promoters. Neighboring
nuclei can exchange monomers and dimers via diffusive hopping. The system is initialized in a five-stripe
pattern of expression domains in the order A–B–C–D–A, corresponding to the experimentally observed
order in the fly embryo. The strength of mutual repression varies among gap gene pairs: genes associated
with nearest neighbor (NN) domains repress each other weakly (dashed arrows), while next-nearest neigh-
bors (NNN) domains exhibit strong mutual repression (thick arrows). By default, the concentration of A
is pinned at the system boundary where the set of modeled reactions differs from the rest of the system
by the fact that the A promoter can not be repressed. Details in the Methods, Sec. 4.

Modeling framework

In order to investigate stability of gene expression patterns without external input gradients, we performed
stochastic simulations of a spatial pattern of four mutually repressing genes, using NS-FFS. Here we opted
for a minimal spatially resolved stochastic model, shown in schematic Fig. 1, inspired by the posterior gap
gene pattern in Drosophila development [49]. The model considers four mutually interacting genes A, B,
C and D, arranged in a five-stripe pattern (with order A-B-C-D-A) along a cylindrical spatial lattice. The
four genes are analogous to the arrangement of the expression domains of gap genes hb, kr, kni and gt
in nuclear cycle 14 in the posterior half of the early fly embryo, where hb is expressed in two stripes, in
the first (anteriormost) and last (posteriormost) stripe [49, 50, 51, 52, 53, 54] ; in the following, we use
the term “expression domain” or just “domain” of gene A for referring to both of the A stripes together.
The spatial lattice consists of Nz × Nϕ equally spaced and well-stirred reaction volumes with periodic
boundary conditions in the circumferential (ϕ-) direction motivated by the arrangement of cortical nuclei
in the developing fly embryo. Protein diffusion and nuclear exchange are modeled via hopping between
neighboring reaction volumes, with a rate proportional to the diffusion coefficient. In each nucleus,
proteins of the genes A, B, C and D are produced from their corresponding promoters, dimerize and
mutually repress each other by promoter binding. Each gene can repress the promoter of each other gene.
Repression is non-competitive, i.e., each promoter has binding sites for each of the three other genes’ dimers
and is inactivated when at least one dimer is bound (“OR”-logics). The model combines transcription
and translation into one production step, neglecting some features of eukaryotic gene expression such as
transcriptional bursts and enhancer dynamics, but previous work has shown that this does not alter the
results qualitatively [59, 14]. We provide a list of the biochemical equations governing the dynamics of

5
our system in Sec. S1.1 and and a graphical summary in Fig. S12 of the Supporting Information.
In the anterior-posterior arrangement A-B-C-D-A, the genes repress each other mutually via the char-
acteristic pattern of strong next-nearest neighbor (NNN) and weaker nearest-neighbor (NN) repression
(alternating cushions), as observed in the Drosophila embryo [43, 49, 54, 60, 61, 62, 63, 64]. Specifically,
there are two pairs of strongly repressing genes, (A,C) and (B,D), and four pairs of genes that repress each
other weakly, (A,B), (B,C), (C,D) and (D,A). In our model, the difference in repression strength is tuned
via the unbinding rates of the repressors from the repressed promoter. The strong-repressor unbinding rate
ksoff is set to a fixed value such that the NNN gene pairs (A,C) and (B,D) are in the bistable regime, while
off
the weak-repressor unbinding rate kw that tunes the repression between NN gene pairs is varied (ranging
from very high unbinding rates, corresponding to very weak repression, towards rates as low as ksoff , which
also brings the NN repressive interactions into the bistable regime). Here being in the bistable regime
means that in individual nuclei only one of the two strongly repressing genes can be expressed at a high
level, while its counterpart is expressed at very low level, e.g. the state in which A is expressed at high and
C at very low levels, or vice versa. Since the repression between the two genes is assumed to be symmetric,
i.e. A and C unbind from each other’s regulatory region with the same rate ksoff (and likewise for B and
D), these two mutually exclusive stable expression states are equally probable without any further inputs
that could break their symmetry, and therefore form a perfectly symmetric ”genetic switch”. Thus, in the
absence of external cues capable of forcing the bistable systems into a preferred state, stochastic switching
is expected to eventually result in one of the domains to dominate over the respective other domain in
the NNN pair, causing its elimination and simultaneous expansion of the dominating gene’s domain. This
partial breakdown of the initial pattern can happen independently for both strongly repressing NNN pairs
and thus in random temporal order; however, ultimately one of the strong repression partners is eliminated
in each of the NNN pairs and the system settles in a new, effectively irreversible state in which only the
remaining two genes are co-expressed.
On the one hand, we expect that the presence of the third expression domain in between the NNN
pair domains can impede elimination of (one of) the NNN pair domains when additional NN repression
is present, because it can spatially move apart the strongly repressing (NNN) expression domains and
form a ”cushion” domain between them, effectively replacing one interface of strong competition by two
interfaces of weak competition that allow for local coexistence of the competitors. On the other hand,
overly strong NN repression is expected to enhance pattern breakdown because then even the overlapping
NN expression domains are brought towards the bistable regime. We therefore study the pattern stability
as a function of the repression strength ratio κ, defined as
off
κ ≡ kw /ksoff , (1)
off
where kw and ksoff are the repressor unbinding rates for weakly repressing NN pairs and strongly repressing
off
NNN pairs, respectively. κ is varied through the weak repression unbinding rate kw . For κ = 1, i.e.
off off
kw = ks , both the NNN and the NN gene pairs are deeply in the bistable regime and repress each other
strongly (because both ksoff and kw off
are low), while in the opposite limit κ → ∞ (kw off
→ ∞) only the
NNN gene pairs form bistable switches whereas the NN pairs do not affect each other at all.

2 Results

Pattern stability is quantified by asymmetry factors

In order to quantify pattern stability, here we define how we understand pattern collapse and construct
order parameters that track pattern destruction by mapping the pattern dynamics onto a low-dimensional
phase space. A typical “intact” spatial pattern of gene expression with (roughly) equally-sized domains is
shown in Fig. 2. We consider patterns in which the expression domain of one gene is lost completely as
being “destroyed” (note that in our terminology the expression domain of gene A refers to both stripes
at the system boundaries). In our system, the strong mutual NNN repression and resulting bistability
effectively prohibit coexistence of the strongly repressing genes at one location. Hence, an increase in
the size of one domain is always accompanied by a reduction in the size of the domain belonging to the
strongly interacting partner. This lead us to introduce the following two order parameters, λAC and λBD ,
here termed asymmetry factors, that measure the asymmetry of the expression domain sizes for each of

6
the two strongly antagonistic NNN pairs:

λAC ≡ max([A]tot , [C]tot )/N , λBD ≡ max([B]tot , [D]tot )/N . (2)

Here [P ]tot is the total copy number of P proteins (counting dimers twice), and N = [A]tot + [B]tot +
[C]tot + [D]tot is the total protein number in the system across all species.

Figure 2: Spatial pattern of gene expression. Snapshots of the total copy numbers of all considered
patterning proteins as a function of the axial coordinate z of the cylinder, averaged over its circumference.
Colors correspond to Fig. 1 (green = A, blue = B, red = C, black = D). Snapshots were taken every
60 min over a total simulated time of 20 h after an initial relaxation phase of 30 min, starting from
rectangular domain profiles of equal length. No-flux boundary condition at either end.

In the spatially well-ordered pattern each protein domain occupies roughly the same fraction of the
system, such that λAC ≃ λBD ≃ 0.25. As expansion of a domain progresses at the expense of its
strong antagonist, λAC (or λBD ) is enlarged and reaches values around 0.5 when the shrinking domain is
eventually lost. In order to track progress of complete pattern losing one of its domains, we use the sum
λ = λAC + λBD , with values around 0.5 for five-stripe patterns and values above 0.75 indicating pattern
breakdown.
Our initial simulations revealed that even for very low protein copy numbers (≲ 20) the waiting times
until one domain is lost are long compared to the duration of the actual breakdown event, and therefore
difficult to sample by direct simulation. This lead us to ask whether the pattern breakdown is merely a
slow random process akin to an unbiased random walk in configuration space, or whether the patterns
have intrinsic restoring capabilities which would counteract the breakdown process; in such scenario, many
more (counteracted) random attempts would be required for concluding the pattern breakdown, effectively
rendering it a barrier crossing problem in which the transition states towards destroyed patterns form the
barrier. In order to resolve which of these alternative mechanisms is responsible for the stabilization of the
expression pattern, we combined our stochastic simulations with Non-Stationary Forward Flux Sampling
(NS-FFS), which is particularly suited for enhanced sampling of non-equilibrium rare events. We used
λ as the progress coordinate for NS-FFS, which aims at generating a branched and weighted trajectory
ensemble that, in the most favorable cases, samples the relevant λ-range uniformly. This allowed us to
generate sufficient statistics of rare breakdown events even in the most stable regions of parameter space
(see Methods, Sec. 4) .
The initial simulations also showed that, in the regime of significant NN repression (small κ), the
expression domains of gene A at the boundaries of the system (green expression domains in Fig. 2) are
particularly prone to destruction by their opponent domains, as we hypothesized for two reasons: Firstly,
they can expand only in one direction, towards the interior of the system; thus, unlike all the other domains,
they cannot compensate shrinkage of the domain at one interface by an expansion at the other interface,

7
and therefore take less effort to get completely destroyed. Secondly, their NN domains (B and D) are
not counteracted by any strongly repressing partner at their interfaces with the A domain, which makes
it easier for them to invade the A domain. We therefore decided to study a setup in which we ”pin”
(keep constant) the level of A proteins at the system boundaries by locally disallowing repression in the
first and last rings of reaction volumes along the z-axis. This setup is not merely an ad hoc modeling
assumption but motivated by experimental findings in the Drosophila gap gene system, which arguably
is the most widely studied example of the alternating cushions arrangement. There, the anterior stripe
of the gap gene hunchback (hb), which corresponds to gene A in our system, not only is under stringent
control by Bicoid (Bcd) but additionally translated from maternal mRNA localized towards the anterior
pole [65, 66]; conversely, in the posterior, zygotic hb expression is driven by a second enhancer under the
control of Tailless (Tll) [67], which in turn is directly controlled by the maternal terminal system and thus
tightly localized [68, 69]. Our “pinning” prescription mimics this biological situation. To assess how the
assumed pinning influences our results, we later compare to simulations in which expression of A can be
repressed at the system boundaries, finding our main results hold up also in this less restricted system. In
particular, pinning is not necessary for enhanced stability but can further increase the maximal stability
time by more than an order of magnitude compared to the system without pinning, as we present further
below.

Long-term pattern stability requires optimal repression strengths

In order to see how varied repression strength affects pattern stability, we reweighted histograms of sim-
ulated trajectories over the reduced phase space spanned by order parameters λAC and λBD at different
times, for different values of κ ranging from strong NN repression (κ ≃ 3) to the limit of non-interacting
nearest neighbors (κ = ∞), see Fig. 3. We found that there exists a region of stable expression patterns
in phase space which is populated rapidly and then remains quasi-stationary, indicating that the system
can remain in a metastable state if NN repression is moderate. In particular, the velocity with which the
system escapes from the quasi-stationary region strongly depends on κ, with low and very high κ resulting
in quick pattern deterioration, and intermediate κ values resulting in the most long-lived quasi-stationary
states.

8
Figure 3: Pattern breakdown in the phase space spanned by asymmetry factors. Probability density
snapshots of the phase space spanned by asymmetry factors λAC and λBC , defined in (2), at different
times t for varied repression strength ratio κ. The conditions are following: strong NN repression, κ = 3.16
(top row), optimal NN repression for pattern stability, κ = 31.6 (middle row), and lack of NN repression,
κ = ∞ (bottom row). The simulation was started with the initial rectangular five-stripe A-B-C-D-A pattern
(λAC , λBD )=(0.4, 0.2) (white circle) in the pinned system. All snapshots are normalized histograms of
reweighted (λAC , λBD )-points within t ± 5 min. In the middle and bottom rows we identify three densely
populated regions: a broad region centered around (0.30, 0.30), RS , which contains five-stripe patterns,
† †
and two smaller regions close to (0.55, 0.30), RAC , and (0.30, 0.55), RBD , representing patterns with one
domain lost (region boundaries (dashed white), details in Methods, Sec. 4). Ultimately trajectories will
converge towards region centered around (0.55, 0.55), R‡ , where two domains are lost.

Stochastic fluctuations can lead to two different events corresponding to partial pattern destruction:
one in which either the A or C domain is lost first and one in which either the B or D domain is lost
first. Motivated by these observations we defined a region of stable patterns in terms of the asymmetry
factors as RS ≡ {(λAC , λBD )|λAC ≤ 0.45 and λBD ≤ 0.43}, see Fig. 3. States that lie outside of RS
† †
are considered deteriorated patterns, and accordingly we also defined two regions RAC and RBD and a
‡
region R accumulating patterns with one expression domain lost and patterns with two domains lost,

9
Figure 4: An optimal strength of nearest neighbor repression maximizes pattern stability. The
mean time until pattern destruction τD as a function of κ, the ratio between the weak and strong repressor
off-rate, for the system in which expression of gene A is fixed at the boundaries. We observe a pronounced
maximum of the stability time when the weak repression is about 30 times weaker than the strong repression
(κopt = 31.6) in the system with pinning at the boundaries (blue line). When pinning of the pattern at
the system boundaries is relaxed (red line), the maximum of stability time moves to κ = 100. The
dashed horizontal lines indicate the values for the completely uncoupled systems with κ = ∞. The dashed
vertical line (orange) shows the optimal repression strength ratio predicted analytically by our stability
theory, κtheor ≃ 76 (see last part of Results section).

RR
respectively. The pattern survival probability S(t) = RS p(λAC , λBD , t)dλAC dλBD is the probability for
the system to remain in the region of stable patterns until time t. We have never observed re-entry into
RS . We found that S(t) is well-described by an exponential decay, S(t) ∝ e−kD t , for times t larger than a
certain lag-time tlag . kD then defines a deterioration rate, corresponding to average pattern stability time
or the mean time until pattern has lost one of its domains, τD ≡ 1/kD (see Methods, Sec. 4).
By quantifying pattern stability time, we found that τD depends strongly on the repression strength
ratio, with a maximum of τD as a function of κ at κopt ≃ 30, see Fig. 4 (blue curve). For κ values close to
κopt pattern stability is still markedly enhanced. While significantly less stable than in the region around
the optimum, patterns with stability time on the order of several hours remain possible in the absence of
NN repression (κ → ∞). In contrast, when NN and NNN repression have close to equal strength (κ → 1)
patterns collapse almost immediately.
Examples of individual trajectories leading to (partial) pattern destruction as they proceed in biased
simulation time (with increasing λ) are described in Section S1.5 of the Supporting Information and shown
in SI-Figs. S7–S9, for the most stable regime, i.e. for κ = κopt in the system with pinning. These examples
demonstrate that multiple destruction pathways are possible in which the individual domains are destroyed
in different order, and that destruction of one domain of a strongly competing gene pair can (but does
not need to) facilitate subsequent destruction of a domain in the other strongly competing pair.

In the maximally stable regime restoring forces reconstitute perturbed patterns

The observation of a phase space region in which system trajectories persist for long times raises the
question whether this region constitutes a true metastable basin of attraction. We first addressed this
question next by analyzing transient behavior of the perturbed patterns. If enhanced phase space density
in certain regions of the (λAC , λBD )-space were indeed due to the presence of a metastable basin, pertur-
bations that transiently drive the system away from the stable pattern should be counteracted by restoring
forces. To test this hypothesis, we perturbed relaxed five-stripe patterns from the hypothetical basin by

10
artificially enlarging domains in which one gap gene is dominant. Using these perturbed states as initial
conditions, we then ran the spatial-stochastic simulator with higher time resolution, and checked whether
the perturbed systems relax back into the presumed basin. We investigated two types of asymmetric
perturbations: “C expansion”, in which the central C domain is unidirectionally expanded at the expense
of the posterior A domain, and the converse “A expansion”, in which the anterior A domain is enlarged
at the expense of the C domain. The perturbation experiments are described in detail in the Methods,
Sec. 4, and Sec. S1.2 and Fig. S1 of the Supporting Information.
We find that at κ = κopt , for both perturbations the perturbed pattern ensembles relax back to their
original positions on a timescale ∼ 10 h (see Supporting Fig. S1). This demonstrates that for optimal
repression strength ratio an effective restoring force counteracts deviations from the five-stripe pattern for
varied λAC . Moreover, this suggests that the probability-enriched region within RS is a real metastable
state confined by an underlying force field. In accordance, the timescale of relaxation is orders of magnitude
shorter than the timescale of pattern collapse. Thus, for κ = κopt pattern destruction is a Markovian
transition between metastable basins with transition waiting times much longer than the timescales of
intra-basin dynamics. In contrast, we could not observe clear restoring behavior in the systems with very
weak or no nearest neighbor interaction. Here perturbations of similar strength tend to result in almost
immediate pattern destruction.
In summary, for the repression strengths ratio κopt ≃ 30 that maximizes stability, pattern breakdown
appears to be an activated process characterized by a restoring force towards the initial state.

Statistical analysis of phase-space dynamics reveals a metastable basin

We further figured that the existence of a true metastable basin should manifest itself also in the statistics of
transient dynamics in phase space. Here the local velocities in the (λAC , λBD ) phase space are particularly
informative: forces that drive trajectories back into basins of attraction should translate into local mean
phase space velocities with a clear bias towards the bottom of the basin.
To extract the velocity field for our system we modeled the coarse-grained pattern dynamics as over-
damped diffusive motion in the ⃗λ ≡ (λAC , λBD ) plane, assuming that these degrees of freedom capture
the slowest time scales of the system and making a Markov approximation for the fast dynamics [70, 71].
This technique has been successfully applied in protein folding [72, 73, 74, 75]. The corresponding model
equation is
q
d⃗
λ = ⟨⃗vλ ⟩ (⃗λ) + 2Dλ (⃗λ)dW
⃗ (3)
dt

where W ⃗ is uncorrelated (2D) white noise with unit covariance. We estimated the local drift ⟨⃗vλ ⟩ (⃗λ) and
diffusion coefficient Dλ (⃗λ) from our reweighed simulated trajectories by averaging local displacements (see
Methods, Sec. 4, and Sec. S1.3 in the Supporting Information). Furthermore, ⟨⃗vλ ⟩ (⃗λ) is proportional to
the effective force acting at the reduced phase space point ⃗λ in the overdamped Langevin model. The local
mean velocity field ⃗vλ (⃗λ) is determined by the conditional transition probabilities π(⃗λ, ⃗λ′ ) between states
⃗λ and ⃗λ′ , and thus can be extracted from our transient simulation data. The resulting average velocity
field in the reduced phase space of (λAC , λBD ) for the optimal repression strength ratio (κ = 31.6) is in
Fig. 5A, and for suboptimal (κ = 1000) is shown in Fig. 5B.
Interestingly, in Fig. 5A one can identify two regions of (λAC , λBD )-space with low average velocities:
†
one within the region of stable states RS , the other within the region RAC of states in which the C
†
expression domain is lost. The region RBD in which either B or D are lost, has no clear boundaries for
optimal κ = 31.6, and only for much larger κ = 1000 a low-velocity plateau is clearly seen in this region
(Fig. 3B). Notably, in the lower-left corner of the RS plateau we notice a small region in which average
velocities are significantly higher and all pointing inwards. We refer to this region as RMB , and identify
it as the metastable basin of intact, relaxed five-stripe patterns. In accordance, the two shown exemplary
perturbed trajectories relax into RMB after randomly exploring the RS plateau, and remain confined to
the RMB for later times (Fig. 5A). However, if the system drifts far away from RMB , in the direction of
† †
RAC , the trajectories are quickly absorbed into RAC once they reach the edge of RS characterized by high
†
velocity components towards RAC .
In order to further investigate the low-velocity attraction basins and the high-velocity ridges that
separate these basins, we use a different representation of pattern asymmetry, defining the shifted difference

11
Figure 5: Phase space velocity fields and passage statistics reveal metastable basins. (A, B)
Average phase space velocity fields for the system with optimal (κ = κopt , A) and suboptimal (κ = 1000,
B) repression strength ratio, in the phase space spanned by asymmetry factors λAC and λBD as defined in
(2). The small subregion with concentrically inwards-pointing velocities towards which perturbed trajec-
tories relax, corresponding to metastable basin of five-stripe patterns, is indicated (RMB , dashed circle).
Velocity fields were obtained by averaging displacements of all trajectories that exit the local bin (see
Methods, Sec. 4). Two examples of trajectories relaxing after perturbations are shown (blue lines = pert.
from boundary, turquoise lines = pert. from center) with their starting points (circles). The boundaries of
phase space regions (thin dashed lines) are as in Fig. 3. Velocity magnitude is indicated with colors. (C)
The velocity field corresponding to κ = κopt for the alternative asymmetry factors (“differences”) δAC and
δBD , as in (4). In C the metastable basin RMB is localized around the “center point” (δAC , δBD ) = ( 12 , 12 ),
corresponding to an intact pattern with equal proportions of strongly competing genes. For additional
clarity, the inset in the upper left corner shows this region without the relaxing trajectories. Note the
almost concentric pattern of velocity vectors pointing towards the center, highlighting the presence of the
metastable basin. The magnitude unit “phase space unit per hour” (PSU/h) is specific to the chosen
asymmetry factors. (D, E) The landscapes of the “pseudopotential” − log p̃ computed from the total
number of phase space trajectories registered in the respective bin of the phase space. The contour plots
to the right of the 3D views show a projected view of the same landscapes. (F, G) Comparison of sections
in λAC and λBD directions, respectively, at λ⊥ = 0.28 between the optimal and suboptimal choice of the
repression strength ratio κ. Here the − log p̃ profile is almost identical in the metastable basin RMB , but
transitions towards the destroyed pattern states face a higher barrier in the system with optimal κ = κopt ,
in both phase space directions.

12
coordinates
1 1 1 1
δAC ≡ ([A]tot − [C]tot )/N + , δBD ≡ ([B]tot − [D]tot )/N + . (4)
2 2 2 2
These coordinates measure the deviation from an intact pattern with equal proportions of the strongly
competing gene pairs, corresponding to the point P0 = (δAC , δBD ) = ( 12 , 12 ) in the phase space of the
shifted difference coordinates, in a way that retains information about which of the antagonistic genes
becomes dominant. Similar latent-space projections have recently proven instrumental in analyzing the
temporal dynamics of the emerging gap gene expression pattern [76]. The corresponding average velocities
in (δAC , δBD )-space are shown in Fig. 5C. The low-velocity basin, corresponding to RMB , occupies the
central part of (δAC , δBD )-space in Fig. 5C. Accordingly, perturbed trajectories relax towards the region
enclosed by concentric velocity vectors pointing towards P0 (the concentric vector pattern is best seen in
the inset), again highlighting the presence of the metastable basin and restoring forces that tend to drive
back fluctuations that perturb the intact pattern. Overall, the finding of a metastable basin is in line with
the Waddington picture of canalization [40, 41], in which developmental stages are seen as successive
attractors of the underlying dynamics with the intact five-stripe pattern considered here representing such
an attractor.
In Fig. 5B we show the average velocity field for the case with weaker NN repression (κ = 1000). Here
the velocity fields are even more plateau-like in the region corresponding to weakly asymmetric patterns,
and the characteristic concentric velocity pattern indicative of the basin in the optimal case cannot be
clearly discerned any more in this case. In accordance, trajectories starting from perturbed patterns do not
relax back and progress towards patterns with at least one domain lost. See also Sec. S1.4 and Figs. S3–S6
for the corresponding velocity fields in shifted difference coordinates and additional alternative projections.
In addition to the average velocity fields of the registered phase space trajectories, the signatures of
the metastable basins are also visible in the local phase space density sampled over many trajectories that
explored the phase space during the whole sampled time interval, p(⃗λ). A suitable quantity for visualizing
the corresponding phase space “landscape” is the negative logarithm of p(⃗λ); note that in an equilibrated,
stationary system this quantity would be proportional to the energy (landscape) defining the stationary
probability distribution of the system. Since our system is genuinely non-stationary, this relationship does
not hold. Nevertheless we can consider our most stable systems transiently equilibrated in the metastable
basins or origin and akin to stationary systems until they irreversibly cross the barrier towards one of the
basins corresponding to destroyed patterns. Note that the depth of these destroyed pattern basins grows
with the amount of simulated time after the destruction events, because then the basins continue to be
explored by phase space trajectories corresponding to fluctuations of the destroyed patterns; their apparent
depth therefore depends on the prescribed maximal duration of the biased simulation trajectories in the
NS-FFS scheme, which is a technical simulation parameter. In contrast, the height difference between
the metastable basin of intact patterns and the barrier separating it from the destroyed patterns basin
is entirely determined by the biophysical parameters that set the average time scale for the stochastic
destruction process. Therefore, the barrier height, as seen from the metastable basin of intact patterns, is
a biophysical property that does not depend on technical choices (but obviously the maximal simulation
time has to be chosen larger than the fastest barrier crossings in order to make the barrier visible).
In Fig. 5D and E we plot the “pseudopotential landscape” defined as − log(p̃(⃗λ)) for optimal κ = 31.6
(D) and suboptimal κ = 1000 (E), where p̃(⃗λ) is a locally smoothened version of p(⃗λ) which equalizes
out small local spikes in p(⃗λ) but preserves the overall structure of the resulting landscape (see Methods
for details). The small plots right of the landscape visualizations show sections through the landscapes in
direction of the asymmetry factors λAC and λBD at chosen constant values of the respective orthogonal
factor (see Fig. 5 caption). In both cases we can clearly identify the metastable basin of undestroyed
patterns and a barrier separating it from the basins of (half-) destroyed patterns. The basin corresponding
to the states in which either the B or D domain is lost is less pronounced for the optimal choice of κ due
to its lower accessibility, and—more importantly—separated by a higher barrier. This is best seen in a
more detailed explicit comparison of the sections through the landscapes, shown in Fig. 5F and G. The
comparison clearly reveals that the barrier separating the metastable basin of intact patterns from the
basin in which the C domain is lost is both higher and wider for the optimal choice of κ, overall leading
to a markedly lower rate of pattern destruction.
Taken together, the analysis of both the velocity fields and the empirically sampled phase space density
demonstrate that the long-time confinement of phase space trajectories close to the five-stripe pattern at
optimal NN repression is due to the existence of a metastable basin which impedes progress towards losing

13
one of the domains by restraining the system from leaving the metastable basin. With decreasing strength
of NN repression the basin gradually disappears, thus enhancing the probability of pattern deterioration.
The finding that pattern stability is enhanced by the emergence of a metastable basin is further
supported by the quantification of the diffusion coefficient in ⃗λ space, Dλ (see Eq. 3 above). As discussed
in detail in Sec. S1.3 and Fig. S2 of the Supporting Information, we observe that the average diffusion
constant in the metastable basin of intact patterns (RS ), ⟨Dλ (⃗λ)⟩RS , monotonically decreases with growing
κ (i.e., with decreasing NN repression strength), meaning that the estimated average time for leaving the
metastable basin by random, “diffusive” motion (shown in Fig. S2B and C) monotonically increases with
increasing κ. Pattern stabilization around the optimal κ therefore cannot be explained by a decrease of
a diffusive escape rate, but rather by the emergence of restoring forces that drive deteriorating patterns
back into the metastable basin.

14
Stability enhancement does not require pinning

To assess whether pinning of the A-domains at the system boundaries is necessary for the observed stability
enhancement at intermediate NN repression, we repeated our simulations and analysis for a system without
pinning. In contrast to the system with pinning, here the promoters of gene A in the nuclei at the system
boundaries can be inhibited by the repressors of A. We found that also in the system without pinning,
pattern stability is markedly enhanced by the presence of weak interaction partners between two strongly
repressing gene domains. In Figure 4 the red curve shows the mean destruction time τD against the ratio
of the repressor off-rates κ for the system without pinning. We again find the highest pattern stability at
(np)
an optimal repression strength ratio κopt = 100 (red curve), which is close to the optimum in the system
(p)
with pinning (κopt ≃ 30, blue curve), albeit with about 10 times lower overall stability times; yet, these
stability times are still about an order of magnitude larger than without fine-tuning of NN interactions.
Overall this demonstrates that enhancement of pattern stability by at least one order of magnitude is
possible both with and without pinning of expression at the system boundaries. However, pinning alters
the proportion of destruction pathways that the collapsing patterns pursue; in particular, it prevents the
destruction of the peripherical A domain, which is the dominant destruction pathway at low and optimal
κ ≲ 100 in the system without pinning. We present this effect in more detail in Sec. S1.6 and Figs. S10
and S11 of the Supporting Information.

An analytical model of expression domain competition predicts optimal pattern stability

The problem of pattern stability has been recently addressed analytically in [55], where general and exact
stability conditions for a pattern of two interacting domains were derived. In that work “stability” refers
not only to the robustness against perturbations, but to the ability of a pattern to survive for infinitely long
time. In this section, we show that these stability conditions can be successfully applied to the multi-gene
system studied in this work, in order to obtain a coarse-grained prediction of the parameter values leading
to pattern stabilization.
The central result reported in [55] is the description of the dynamics of a contact zone between two
gene-expression domains for various levels of mutual repression between the two expressed genes. A single
expression domain can form either by overcoming the “activation threshold” in the nearby undifferentiated
tissue, resulting in asymptotically constant-velocity expansion, or emerge instantaneously in the entire
available tissue, when expression is constitutive (active by default). For two genes in the system (and
two respective domains) the scenario depends on the strength of mutual repression. If one gene cannot
prevent the expression of the other gene in the bulk of its own domain, the dominating gene overtakes
the system exponentially fast, expressing in the entire volume and without forming a meaningful contact
zone between domains. For stronger repression, which prevents gene expression deeper in the bulk of
its adversary domain, a contact zone emerges, within which both domains of active expression overlap.
However, this region of overlap grows indefinitely, albeit with asymptotically constant velocity. When the
interaction strength surpasses a critical value, an asymptotically finite-size contact zone is formed. In this
regime one domain can still shrink and the other grow, but in a coordinated manner, preserving the width
of the contact zone. Asymptotically, the contact zone drifts with a constant velocity that is determined by
the system parameters. This gives rise to a “travelling” gene expression pattern. The width and velocity
of the contact zone are stable against perturbations in this phase, acting as an attractor of the system
dynamics. However, while the travelling pattern is well-organized into two domains, it is not stable in
the sense that in finite-size systems it survives only for a limited time, until one domain “pushes out”
the other. Finally, perfectly stable patterns arise as a special case of travelling patterns with zero-velocity
drift. As such, they can survive arbitrarily long.
The simulations in this work are stochastic, tracking the chemical reactions at single-molecule reso-
lution across the set of reaction-volumes constituting the system. However, in the limit of large particle
number and small reaction volumes, this type of spatially discrete and stochastic dynamics approaches the
continuous and ultimately deterministic reaction-diffusion dynamics of the type considered in [55]. The
existence of this deterministic limit can be also seen as the manifestation of the emergent noise-control
mechanism that overtakes the system. Therefore, we compare the numerically found optimal κopt with
the theoretically predicted κtheor to assess how well the deterministic theory approximates the dynamics
in the highly stochastic regime, and to explain the nature of the emergent noise-control mechanism.
To this end, we mapped the microscopic model used in our stochastic simulations onto the effective
reaction-diffusion model analysed in [55] (see Methods). We obtained the following continuous model for

15
the expression of each gene X in {A,B,C,D}:
 
 
 X 
∂t X2 (x, t) = D ∂xx X2 (x, t) − γX2 (x, t) + H θ 
1 − ϵXY Y2 (x, t)
 , (5)
 Y∈ 
{A,B,C,D},
Y ̸=X

where X2 (x, t), Y2 (x, t) are the concentration profiles of the protein dimers (indicated by the subscript
2), D is the diffusion constant, γ the degradation constant, H a production constant, and ϵXY are gene-
gene interaction strengths. θ(. . . ) denotes the Heaviside step function, corresponding to steep Hill-type
regulatory
p kinetics. Note that the derivations in [55] only apply to systems with size L ≫ λ, where
λ ≡ D/γ is the characteristic length of gene interaction. For the systems studied here, λ ≈ 8.62 µm,
which is much smaller than the system size L ≃ 340 µm, warranting application of the theory.
While the original theory in [55] describes only the contact zone involving exactly two domain bound-
aries, we can adapt it to the four-gene system studied here. Fig. 2 shows that in the alternating cushions
system there are only two types of contact zones: (i) between two strongly interacting genes (NNN do-
mains) with the third, weakly interacting gene (NN domains), expressed in the background or (ii) between
two weakly interacting genes (NN domains), with all other genes having close-to-zero expression level.
Thus, we will consider stability of both contact zone types separately.
In the type-(i) contact zone, the dynamics of gene expression is described by the effective equations
  
 −1 H −1
 ∂
 t 2X (x, t) = D∂ X
xx 2 (x, t) − γX 2 (x, t) + Hθ 1 − K w − K s Y2 (x, t)
γ
  (6)

 H
 ∂t Y2 (x, t) = D∂xx Y2 (x, t) − γY2 (x, t) + Hθ 1 − Kw−1 − Ks−1 X2 (x, t)
γ

where we approximate that the third “background gene”, has a constant expression level over the contact
zone. The equilibrium value of this expression level is H/γ. Kw and Ks are the weak and strong repression
constants, respectively and they satisfy (cf. Eq. 1 and Methods, Sec. 4):

Kw k off
κ= = w . (7)
Ks ksoff

Type-(i) contact zones are established between genes A and C (with B or D in the background) as well as
between B and D (with C in the background). In the type-(ii) contact zone, the equations take the form:
(

∂t X2 (x, t) = D∂xx X2 (x, t) − γX2 (x, t) + Hθ 1 − Kw−1 Y2 (x, t)

(8)
∂t Y2 (x, t) = D∂xx Y2 (x, t) − γY2 (x, t) + Hθ 1 − Kw−1 X2 (x, t)

This contact zone emerges between gene pairs (A,B), (B,C), (C,D), and (D,A).
Having defined these contact zones, in Method section we adapt the more general derivation of stability
conditions from [55] for the current case. The main idea of this derivation is that the shape of expression
profile X2 (x, t), defined by Eq. (5), can be found without knowing where the domain boundaries are
located. Then, the positions of domain boundaries are sought from a separate set of equations. In [55],
it is shown that the boundaries asymptotically travel with the common constant velocity v, preserving the
distance ∆r between them. This ansatz leads to the algebraic set of equations defining v and ∆r, which
can be solved. Eventually, the conditions for pattern stability are equivalent to ensuring that v = 0 is the
correct solution. Let us define two sets of constants, one for each type of contact zone:
H
(i) C̃X = C̃Y = 1 − Kw−1 , ϵXY = ϵYX = Ks−1 , ϵXX = ϵYY = 0 , (9)
γ
(ii) C̃X = C̃Y = 1 , ϵXY = ϵYX = Kw−1 , ϵXX = ϵYY = 0 . (10)

and two auxiliary variables:

2γ C̃X 2γ C̃Y
RX = −1 , RY = −1 , (11)
ϵXY H ϵYX H

16
Then, the stability conditions, as derived in the Methods section, come down to:

RX = RY = R , −1 ≤ R ≤ 1 , (12)

Additionally, the width of the stable contact zone reads:

∆r = −sgn(R)λ ln(1 − |R|) (13)

In order to determine the range of parameters ensuring the global stability of pattern, we apply the
stability conditions (12) separately to type-(i) contact zone (Eqs. (9)) and type-(ii) contact zone (Eqs.
(10)). One can notice that in each type of contact zone the equality RX = RY is automatically satisfied,
due to the common choice of parameters γ, H and D for both genes, as well as the symmetry in gene
interactions. Sharing the same parameters is also the main reason why stability conditions (12) are much
simpler than their general counterpart reported in [55]. RX = RY = R means that each type of contact
zone is characterized by one variable:

2(1 − Kw−1 H
γ ) 2γ
R(i) = −1 R(ii) = −1 (14)
Ks−1 H
γ Kw−1 H

The remaining stability condition, −1 < R < 1, applied to R(i) and R(ii) , results in the following
inequalities:
H 1
(i) Kw ≥ , Ks ≤ H −1 ,
γ ( γ ) − Kw−1
(15)
H
(ii) Kw ≤ .
γ
These conditions show that the addition of weak interactions is instrumental for increasing system stability.
On the one hand, the type-(i) contact zone is stable (i.e., immobile) provided that the weak interaction
−1
strength Kw−1 does not exceed (H/γ) ; otherwise it would prevent the expression of strongly interacting
−1
genes in this region. On the other hand, for the type-(ii) contact zone it is necessary that Kw−1 > (H/γ) ,
as this minimal strength of repression is required to prevent co-expression of both weakly interacting genes
in the same region. In order to simultaneously stabilize both types of contact zones, one needs to negotiate
between these two largely opposite goals. This trade-off can be achieved only for the most marginal value
in both parameter ranges, Kw = H/γ, which highlights why in the alternating cushions architecture the
weak interactions have to be fine-tuned for pattern stability. In contrast, but in line with the numerical
findings, the strong interactions characterized by Ks can be arbitrarily large, Ks ≤ +∞.
The simulations in this work were performed for Ks ≃ 0.003 µm−3 with Kw varied to obtain different
values of κ, see Methods, Sec. 4. Calculated from these microscopic parameters, H/γ ≃ 0.23 µm−3 . The
resulting theoretical value of κ that ensures stability is then κtheor ≃ 76. This number is of the same order
of magnitude as the optimal κ in the simulated stochastic systems, showing slightly better agreement with
(np) (p)
the no-pinning case (κopt ≃ 100) than with the case with pinning at the boundaries (κopt ≃ 30), see
Fig. 4.

The emergent noise-control mechanism can be understood via the analytical model

The analytical deterministic model can be employed to obtain further insights into the mechanism of
increased pattern robustness against noise in the vicinity of optimal κ. For this, we must first consider
the width of type-(i) and type-(ii) contact zones in their stability regions predicted by the theory from
[55]. Inserting R(i) and R(ii) into Eq. (13) with Kw = κKs , we obtain ∆r(i) and ∆r(ii) as functions of
κ, shown in Fig. 6A. Here, ∆r > 0 indicates a no-expression region between the domains (a gap), while
∆r < 0 means that active expression regions overlap. One can instantly notice that ∆r(i) → +∞ and
∆r(ii) → −∞ at κ = κtheor . Tending to infinite values is an artefact of our analysis, in which we treat
each contact zone as a separate region, disconnected from the others. However, this behaviour conveys
an important message. At κ = κtheor the system attempts to maximize the size of each contact zone,
forming five contact zones tightly filling the entire system. In this state, any pattern perturbation distorts
at least two contact zones. Since each contact zone is stable, their maximized widths are attractors for
the deterministic dynamics [55], and consequently the system tends to remove the perturbation. This is
the origin of increased survival time of patterns at optimal κ. This restoration behavior is qualitatively

17
similar to the model of repulsive forces between domain boundaries (kinks), discussed in [43]. Although the
analytical stability theory [55] does not rely on the concept of explicit restoration forces, these forces arise
effectively, leading to the occurrence of the pseudopotential in our phase space analysis, in Fig. 5. Thus,
the effective restoration forces form a link between the exact stability theory [55] and the approximation
of interacting kinks in [43].

Figure 6: Analysis of theoretical contact zone widths uncovers a deterministic stabilization mech-
anism. (A) Plot of theoretical contact zone widths in the approximation of separate interfaces, ∆r(i) and
∆r(ii) , for type-(i) contact zones (between two NNN domains with third interacting gene expressed in the
background) and type-(ii) contact zones (between NN domains), calculated from Eq. (13) and Eqs. (14)
with κ = Kw /Ks and Ks kept constant. In their respective regimes of stability, the widths are restored by
deterministic dynamics if perturbed. Vertical lines: κ0 ≃ 38 (red, dashed) at which ∆r(ii) changes sign;
critical κtheor ≃ 76 (black, dashed) ensuring simultaneous stability of type-(i) and type-(ii) contact zones.
Horizontal line (gray, dashed): limit of ∆r(i) ≃ 31.3 [µm] without any weak interactions (κ → +∞).
(B-E) Schematic representations of system states in various regimes of κ, predicted by the deterministic
model. Solid arrows: stable contact zones (restorable width); dashed arrows: unstable contact zones
(non-restorable width); inward arrowheads indicate ∆r(ii) < 0; contact zones of type-(ii) (red), and type-
(i) (blue). (B) 0 < κ < κ0 : type-(ii) contact zone stable, ∆r(ii) > 0, no type-(i) contact zones, domain
widths lack stabilization against fluctuations. (C) κ0 < κ < κtheor : type-(ii) contact zone stable, partial
overlap of domains, ∆r(ii) < 0, provides minimal domain width stabilization against fluctuations, but
fluctuations can shift entire contact zones. (D) κ = κtheor : type-(i) and type-(ii) contact zones stable,
maximizing their widths (∆r(i) and ∆r(ii) tend to ±∞ in the approximation of separate contact zones).
Pattern is restored after any perturbation. (E) κ > κtheor type-(i) contact zones stable, but ∆r(i) ≪ L,
fluctuations can shift entire contact zones.

The existence of a rigorously sharp stability condition, Kw = H/γ, raises the question about the
deterministic dynamics for suboptimal choice of κ and its influence on the stochastic system. Let us first
consider the case κ < κtheor , in which type-(ii) contact zones are stable. In this regime, the system
forms a pattern of domains A-B-C-D-A, but the NNN domains are so distant from each other that strong
interactions are not yet important. There are four type-(ii) contact zones in this system. For κ ≃ 0, the
weak interactions are extremely repressive and ∆r(ii) → +∞. Thus, the pattern collapses. For somewhat
larger κ, a finite-size gap (∆r(ii) > 0) between NN domains emerges (see Fig. 6B) and is reduced to
zero width (∆r(ii) = 0) at κ0 ≃ 38. In this regime, the pattern can survive arbitrarily long in the absence
of fluctuations, but the domain widths are not stabilized in any way. Thus, in the presence of noise, the
survival time of the domain depends on its size (which grows with κ), as larger domains take longer to
be destroyed. For κ > κ0 , NN domains begin to overlap, as ∆r(ii) becomes negative (see Fig. 6C). This
marks the first emergence of the additional stabilizing mechanism, as the deterministic dynamics will tend
to restore ∆r(ii) in each contact zone if perturbed. This means that, in the presence of noise, ∆r(ii)
would keep returning to its deterministic value, but fluctuations can still shift a stable contact zone as
one entity. If, as a result, two contact zones meet or one is pushed to the system boundary, this causes
the collapse of a domain and partial desintegration of the pattern. As κ further approaches κtheor , the
overlap becomes large enough such that NNN domains begin to interact and type-(i) contact zones are
formed. These contact zones have a certain minimal width, but they are not stable, in the sense that this
width would not be restored if increased. The occurrence of type-(i) contact zones imposes a barrier for
the further growth of ∆r(ii) with κ (see Fig. 6D). At this stage, a major enhancement of pattern stability
occurs, as stable type-(ii) contact zones and type-(i) contact zones tightly fill the system. At κ = κtheor

18
also the type-(i) contact zones gain stability, due to the further increase of overlap between NN domains.
This results in the maximal global robustness of the pattern against noise. For κ > κtheor , the width of
the now-stable type-(i) contact zones quickly decreases and eventually saturates at ∆r(i) = 31.3 [µm] in
the limit of completely absent weak interactions (κ → +∞). In this regime, with type-(i) contact zones
width ∆r(i) < L/3 and no mechanism restoring the width of type-(ii) contact zones (see Fig. 6E), the
pattern gradually loses stability against noise. Fluctuations can shift each type-(i) contact zone as one
entity (analogous behaviour to the type-(ii) interfaces in κ0 < κ < κtheor regime), eventually leading to
pattern destruction.
In summary, the emergent deterministic dynamics, described in [55], is crucial for stabilizing the highly
stochastic system simulated in this work. The increase in survival time τD , towards κopt , illustrated
in Fig. 4, is directly associated with the gradual activation of deterministic stabilization mechanisms,
described in the paragraphs above. General principles of pattern stabilization, outlined in [55] for two
genes, apply also to the four-gene system studied here, but many-gene competition and stochasticity
results in a more nuanced picture of stabilization. A more detailed investigation would require considering
the full spatial variability of all expression profiles together, but the approximated effective model proves
useful in predicting optimal parameters.

3 Discussion
In many developing organisms, morphogen gradients provide a long-range positioning system by activating
downstream patterning genes in a concentration-dependent manner. Prominent examples are the gap
gene system in Drosophila, whose main maternal regulators are the morphogen gradients of Bcd and Cad
spreading along the embryo axis [53, 49, 77, 78, 79, 80, 81, 82, 83], and the vertebrate neural tube with
Shh and BMP/Wnt secreted from the opposite sides of the neural tube [4, 7, 8, 18, 20, 23]. For the
Drosophila embryo, multiple studies have shown that mutual interactions between gap genes play a crucial
role in abdominal segmentation [14, 49, 51, 52, 53, 48, 84, 85, 58], leading to the formation of stable
domains with slow effective dynamics [43]. However, it remains unclear how such a system could be
robust given the stochastic nature of gene expression and regulation if the emergent interactions are not
fine-tuned to mitigate the resulting noise. Moreover, it is observed that maternal regulators such as the
Bcd gradient disappear while the expression patterns invoked downstream persist [86, 87]. In support of
the view that self-coordination properties emerge in the gap gene system after maternal activation, a more
recent study which found that the gap gene expression pattern scales with the size of the embryo with
high precision, while—surprisingly—the Bcd gradient does not display any scaling properties [88]. Similar
emerging self-organizing properties have been observed in other developmental systems [15, 32, 89].
Here we asked whether a system of mutually repressing developmental patterning genes arranged in
successive expression domains can indeed be stable over developmentally relevant time intervals without
upstream morphogen gradients while facing unavoidable fluctuations in the expressed gene products. Such
copy number fluctuations can induce bistable switching at the domain boundaries, resulting in stochastic
movement of the boundary which ultimately can lead to destruction of one of the gene expression domains.
We quantified the mean stability time of a five-stripe expression pattern formed by four interacting genes
in a stochastic model conceptually inspired by the posterior Drosophila embryo in cycle 14 as a function of
the repression strength between neighboring stripes. To be able to simulate the breakdown of very stable
patterns we employed Non-Stationary Forward Flux Sampling (NS-FFS), an enhanced sampling scheme
for simulating rare events in non-stationary systems with transient dynamics [57]. We find that for an
optimal value of the repression strength between adjacent expression domains the stability of the pattern
is increased by about an order of magnitude. This stability optimum can be traced back to the fact that
bistable switching at the boundary between domains of strongly mutually repressing genes is inhibited
by an intervening cushion domain of a gene that weakly represses both strong partners. This stabilizing
mechanism works best if the spacer gene represses its nearest neighbors (NN) with moderate strength: very
weak NN repression has no effect while strong NN repression globally destabilizes overlapping domains.
At the optimal repression strength (κ = κopt ) the cushion thus slows down the random motion of the
domain boundary and subsequent pattern destruction.
Stability is enhanced even more, by one more order of magnitude, if expression of the outermost
gene is pinned at the system boundaries, which effectively anchors the whole expression pattern. Such a
situation may emerge when the outermost gene remains under control of maternal cues, such as maternally
deposited mRNA, while the other gene stripes form only by zygotic interactions. Furthermore, it resembles

19
the late stages of neural tube development in which the Shh and BMP morphogen gradients are acting
close to system boundaries [31, 16]. In the system considered here we find that five-stripe patterns form
a metastable attractor of the dynamics with a restoring force that counteracts perturbations, such as
non-perfect initial conditions. In the optimal stability regime, our observations are consistent with the
Waddington picture [40, 41] of development as canalization into successive metastable states, with the
ordered initial gap gene pattern representing one of the metastable states in this succession. Earlier
work already demonstrated that developmental attractors may emerge as an intrinsic property of the
gene expression pattern established through mutual interactions [52, 58]. Here, we demonstrate that
even without morphogen gradients metastable basins can arise and protect expression patterns against
stochastic fluctuations.
Further insight comes from the application of the stability theory derived in [55] to the model of
four genes interacting in the alternating cushions scheme. In agreement with the simulations, these
analytical calculations reveal that the presence of weak interactions is necessary for stabilizing the system
and establishing long-surviving patterns. More specifically, theoretical analysis shows that requirements
for stability of type-(i) contact zones (i.e. two strongly interacting genes with the third weakly interacting
in the background) and type-(ii) contact zones (i.e. two weakly interacting genes with other genes at very
low expression level) are to certain degree incompatible, and agreement between them can be achieved
only for the most marginal value of κ = κtheor in the respective stability range for each type. As a
consequence, simultaneously ensuring perfect stability of both contact zone types requires fine-tuning of
the weak repression strength, quantified by the corresponding dissociation constant Kw . This analytical
prediction of one optimal value of κ is in qualitative agreement with the numerical simulations, which
show a very sharp rise in the survival time of expression pattern near one particular value of κ = κopt , see
Fig. 4.
Quantitatively, the numerical κopt and theoretical κtheor agree particularly well in the no-pinning case
(p) (p) (np) (np)
(for pinning: |κopt − κtheor |/κopt = 139%, for no-pinning: |κopt − κtheor |/κopt = 24%). This is in
line with the assumptions of [55], where an open system was considered and system boundary effects,
such as pinning, were neglected. Differences between κtheor and κnum are expected due to the nature of
approximations employed in the mapping of microscopic model on its effective representation (5). It is
plausible that this discrepancy could be resolved by constructing an even higher-level stability theory that
takes into account the spatial variability of all four genes in each contact zone.
Further, using the division into type-(i) and type-(ii) contact zones, we investigated the behaviour
of effective deterministic model Eq. (5), in the entire range of κ. We found that the preference of
the system to form possibly large contact zones, combined with the stability of at least one type of
interfaces between domains, results in the increased robustness of the pattern against fluctuations, in the
vicinity of optimal κ. These observations are in agreement with our highly stochastic and microscopically
detailed simulations, for which the deterministic model is only the continuous-limit approximation. Yet,
the approximate agreement between κopt and κtheor as well as the broad peak of increased survival time
(Fig. 4), suggest that the deterministic dynamics of model (5) is still remarkably important for this system.
The interplay between deterministic and stochastic component of dynamics in simulations results in the
emergent noise-control mechanism, significantly increasing survival time of patterns. We also found that
the shifting of the stable contact zones by fluctuations is the major reason of pattern destruction for κ
away from the optimal value.
The observed stability times appear sufficient for early fly embryogenesis (≃ 2h until cycle 14) for all NN
repression strengths weaker than the optimal value, with or without pinning, even for the reduced system
size considered here for computational feasibility. We expect that the stability times will systematically
increase when a larger system size is chosen in a more realistic description. Note that the system size
can increase in two ways, either by increasing the considered lattice of reaction volumes (nuclei) or by
allowing for a larger maximal copy number per reaction volume. In the first case, stability is enhanced
because the expression states of more reaction volumes need to be switched in order to destroy the
now larger expression domains, while the local molecular noise level (which is a key determinant for the
speed of this process) remains the same. In the second case, the molecular noise is reduced, such that
detrimental cell switching events are impeded, leading to longer average destruction times for the individual
expression domains and consequently overall longer stability. Nevertheless, based on the theoretical and
numerical evidence we believe that the stability enhancing mechanisms uncovered in this study will also
apply to biologically relevant system sizes. Other factors potentially affecting stability are autoactivation
interactions and interactions with other genes not included in the simplified regulatory network studied
here, which will likely affect the dynamics of the gene expression pattern. Moreover, in our system the

20
parameters do not differ between the four interacting genes (with the exception of the repressor unbinding
rate from the A promoter at the boundaries in the system with pinning); as expected, we therefore did not
observe any directional preference for domain shifts and destruction events, unlike in Drosophila, where
the gap genes exhibit a systematic anterior shift in early developmental cycle 14 [90, 50, 46] (which appear
to require the action of shadow enhancers not considered here [91]). The possible effect of asymmetric
regulatory interactions on pattern stability is an interesting open question that could be assessed in future
iterations of the model and theory presented here. Note, however, that the adapted stability theory clearly
identifies the weak nearest-neighbor repression strength as the key parameter for enhancing stability, while
the strong repressive interactions are found not to affect stability as long as they are chosen strong enough.
Furthermore, since in the theory the predicted value of the optimal repression strength ratio is entirely
determined from the properties of the contact zones between the expression domains, this prediction does
not depend on the spatial system size, provided that it is large enough as to accommodate all the contact
zones. We therefore expect that increasing the size and realism of our spatial-stochastic model in the
described ways would alter the recorded stability times, but at the same time retain the observed key
property of strong stability enhancement at an optimal repression strength ratio.
Our work puts an interesting perspective on the role of maternal gradients in establishing and main-
taining developmental patterns. We show that sufficiently stable patterns can exist without morphogen
gradients, but at the same time that their stability is significantly enhanced by pinning the patterns at
the embryo boundaries. Taken together, this suggests that morphogens do not act deep inside the embryo
interior, which could explain why the patterns remain stable even when the morphogen inputs disappear
[86, 87]. Instead, they may predominantly act at the embryo boundaries as to break symmetry, by selecting
the desired pattern from the larger set of patterns that, by permutation, would also be stable. By acting
only at the periphery, the morphogens, which themselves do not exhibit scaling, still would allow scaling
of the downstream pattern with embryo length, in line with recent findings [88].

21
 The stabilizing mechanism arising from fine-tuning nearest-neighbor interactions in the alternating
cushions scheme can be also considered in the broader class of regulatory mechanisms providing pattern
stability against intrinsic and extrinsic noise [11, 21, 29]. In future studies, it may be instrumental to
further numerically and analytically explore the proposed model by including other biologically relevant
features. Possible extensions include growth of the tissue by cell divisions, self-correcting mechanisms
through cell-to-cell communication other than diffusive exchange of proteins, or inclusion of more specific
noise types. Another interesting scenario would be studying the alternating cushions system considered
here but with periodic boundary conditions, possibly relevant to sea urchins and sea stars that develop
pentaradial symmetry in later stages of development [92]; this likely would allow for further increase of
pattern stability in the optimal repression strength regime without any pinning, as it would provide an
alternative way of stabilizing the A domain. These extensions could further test the validity of our stability
theory under more realistic biological conditions. However, due to the remarkable agreement between our
adapted stability theory and the numerical simulations of the minimal model studied in this work, we
believe that more realistic variants of it will result in quantitative but not qualitative changes in our
predictions.

4 Methods
Details of the model
Our model is inspired by arguably the most paradigmatic developmental system in which development
of distinct cell fates is determined by local protein expression patterns driven by external morphogen
gradients, the early embryo of the fruit fly Drosophila melanogaster. We model the egg-shaped embryo
with its cortical layer of nuclei as a cylindrical array of reaction volumes coupled by diffusion of proteins.
Every volume (nucleus) contains four individual promoters for each of the genes A, B, C and D. Each
promoter can be repressed by the products of the three others with different affinities; this system of four
mutually inhibiting genes represents the gap gene system in the early fly embryo, formed by the four genes
hb, kr, kni and gt, and comprises its essential regulatory interactions. For combined repressive interactions,
we employ OR-logic, i.e. whenever one of the three repressor sites is occupied expression of the gene is
completely blocked. There is no competition for repressor sites on the promoters. In the unrepressed state
the promoters exhibit constitutive protein production, i.e. no external activator signal is required. This
deliberately mimics a situation in which activation of the genes is not provided by external morphogen
gradients but by either an omnipresent master activator or auto-activation with a low activation threshold.
Consequently, our model explicitly does not include morphogen gradients. As a simplifying assumption,
we treat the whole production process, i.e. transcription, elongation and translation, as one step governed
by a single rate β. Proteins however can form (homo)dimers and dedimerize again [93, 94], and only in
their dimeric form they act as repressors. This is to ensure that antagonistic genes form bistable pairs
for sufficiently strong mutual repression. Initially, all simulations are set up in a stripe pattern similar to
the experimentally observed order in the embryo posterior, i.e. A-B-C-D-A [53, 49, 54]. This implies a
fixed definition of “gene neighborhood” to which we refer throughout this paper: by nearest neighbors
(NN) we mean the pairs (A, B), (B, C), etc., while the pairs (A, C) and (B, D) are considered next-
nearest neighbors (NNN). A key ingredient of our model is that nearest-neighbor repression is weaker than
repression between next-nearest neighbor domains (see “Parameter choice” in Methods). By default we
pin the expression of A at the system boundaries, i.e. in nuclei on the two outermost rings of the cylinder
the A promoter is irrepressible, and therefore constitutively produces A proteins. This is motivated by the
fact that in the real Drosophila embryo the gene Hb is under strict control by the maternal morphogen
Bcd throughout the anterior half [65], while in the posterior a second enhancer exposes Hb to positive
regulation by the maternal terminal system [67, 68, 69]. We compare this system to a system in which
there is no pinning and all nuclei are identical.

Simulations
To perform rare-event sampling of the spatially resolved system we integrate our “Gap Gene Gillespie”
(GGG) simulator used in previous work [14, 59] with the NS-FFS scheme [57]. NS-FFS is used to monitor
and process a progress coordinate written out by GGG at regular simulation interrupts, at which GGG
trajectories are cloned and restarted in a way that sampling is enhanced in the direction of increased
progress coordinate, i.e. towards pattern destruction.

22
Spatially resolved stochastic simulations (GGG)
In GGG, the model is implemented via the Stochastic Simulation Algorithm by Gillespie [95, 96] on a cylin-
drical 2D lattice of reaction volumes at constant distance l = 8.5 µm, with periodic boundary conditions
in the circumferential direction of the array. An abstract graph of the reaction network that displays the
set of reactions for any of the simulated promoters is shown in Figure S12 in the Supporting Information.
Diffusive chemical species (patterning gene proteins and their dimers) hop between neighboring volumes
via the next-subvolume method [97] which integrates diffusion into the Gillespie algorithm by annihilation
of a species copy in the volume of origin and instantaneous insertion of that copy in a randomly chosen
neighboring volume with a rate kdiff = 4DP /l2 , where DP corresponds to the protein diffusion coefficient.
The source code of GGG can be downloaded from [Link]

Forward flux sampling

We employ the recently developed non-stationary forward flux sampling (NS-FFS) method [57, 56, 98] to
enhance stochastic sampling of system realizations that increase a (reaction) progress coordinate λ while
retaining correct statistical weight. NS-FFS achieves this by branching off multiple child trajectories upon
crossing predefined interfaces in undersampled regions of (λ, t)-space and pruning trajectories that cross
interfaces in oversampled regions. The NS-FFS scheme aims at equilizing the flux of simulated trajectories
in the reaction coordinate direction among the time bins. The rate of branching and pruning is calculated
from the temporal trajectory crossing statistics collected during runtime. To that purpose the time domain
is subdivided into equidistant time intervals. For a detailed account of the reweighting procedure we refer
to [57].

Progress coordinates
The choice of a suitable progress coordinate is a critical step of the FFS technique. Here, we seek to
enhance progress of the simulated patterns towards their destroyed state. The destruction events are in
particular characterized by the disappearance of one of the partners within each of the strongly repressing
gene pairs. Progress towards destruction thus is accompanied by increasing pair asymmetry, which can be
quantified for each pair separately by the following two asymmetry factors:

λAC ≡ max([A]tot , [C]tot )/N (16)

λBD ≡ max([B]tot , [D]tot )/N (17)

where N = [A]tot + [B]tot + [C]tot + [D]tot is the number of all proteins in the system. Based on this we
define our progress coordinate, which increases whenever asymmetry among any of the pairs is augmented,
via

λ ≡ λAC + λBD = [max([A]tot , [B]tot ) + max([B]tot , [D]tot )] /N . (18)

Since NS-FFS features multi-dimensional reaction coordinates we compared our standard choice to a
setup in which the two components λAC , λBD of the reaction coordinate λ are treated as two separate
reaction coordinates with an own set of interfaces each. While an orthogonal pair of reaction coordinates
captures the principal reaction paths in our system more accurately, the acquisition of crossing statistics
is prolongated because of the increased number of bins in these simulations, and we did not find any
substantial advantage of this choice in terms of branching behavior. We therefore preferred the standard
definition.

Combination of simulation methods

In order to wrap NS-FFS around the GGG simulator we run GGG for a predefined simulation time tGGG =
60 s. At the end of the simulation the reaction coordinates are calculated and passed on to the NS-FFS
module, and the end state of the simulation is recorded. The NS-FFS module then determines whether an
interface crossing has occurred and, if so, decides on whether the trajectory shall be branched or pruned.
In case of branching NS-FFS will prompt nB ≥ 1 restarts of the GGG simulator with the recorded end
state as initial condition, different random seeds and with new statistical weights. At each crossing and at
measuring times spaced by a regular interval ∆t the time, branch weight and reaction coordinate values
are stored in a tree-like data structure that facilitates later analysis.

23
 Trajectory trees are started from a standardized, regular-stripe initial condition passed to the first call
of GGG. Propagation of the tree stops when all child branches have either reached the end of the time
histogram or have been pruned. Subsequently a new tree is started with a different random seed. NS-FFS
monitors the cumulative simulated time Tcum and terminates simulation when Tcum exceeds a predefined
maximal simulation time Tmax and the last trajectory tree has been propagated towards the end. Typically,
Tmax = 3 − 7 h and Tcum = 2 − 5 · 107 s, which usually results in several thousand independent starts
from the initial condition. By default we start from an artificial pattern consisting of five non-overlapping
stripes with rectangular profiles occupying an equal part of the total system length L/5 each and equal
number of monomers (no dimers) in each nucleus close to the expected total copy numbers. We find that
these initial patterns quickly relax towards typical metastable patterns, i.e. into the metastable main basin
of attraction, which justifies our approach a posteriori.
The source code of the NS-FFS path-branching algorithm (Pabra) combined with GGG can be down-
loaded from [Link]

Parameter choice
Repression
We are mainly concerned about the importance of distinct repression strength of nearest-neighbor (NN) as
compared to next-nearest neighbor (NNN) interaction. We assume repressor binding-rates to be diffusion-
R
limited via kon = 4πσR DN , where DN is the intranuclear diffusion constant and σR an effective target
radius. Repression strength therefore is varied by changing the unbinding rates of the repressing dimers.
off
The main parameter in our simulations is κ = kw /ksoff , the ratio between NN and NNN repressor off-
off off
rate. In this work only kw is varied, while ks is chosen sufficiently low to guarantee bistability between
next-nearest neighbor genes, which is a precondition for the formation of individual stripe domains in the
first place, see Table S1 in the Supporting Information. For κ = 1 NN and NNN repressive interactions
are equally strong, while for large κ values NN repression is much weaker than NNN repression. In the
“uncoupled limit” κ → ∞ the two bistable pairs coexist without sensing each other. We do not consider
cases with κ < 1.

Dimerization
D
We set the dimerization forward rate kon to be equal to two times the diffusion-limited repressor binding
rate, which is accounting for the fact that both reaction partners are diffusing. The dimerization backward
D D
rate is set via koff = kon /VN (VN = nuclear volume) as in [99, 100, 101, 14] to ensure that at any moment
most of the proteins are dimerized.

Production and degradation

In our model both monomers and dimers are degraded. This leads to a nontrivial dependence of the total
copy number on production, degradation and (de)dimerization rates, as we discuss with more detail in
[14]. Since we did not find any experimental reports of gap protein lifetimes, we chose equal monomeric
(µM ) and equal dimeric degradation rate (µD ) for all genes and set these quantities to values that lead
to a reasonable effective lifetime of the corresponding proteins of teff ≃ 100 s. The steady-state copy
number is tuned via the production rate β. By default, we consider copy numbers as low as possible
(≃ 15) to minimize computational effort. The effect of increasing the average copy number is discussed
in the “Discussion” section.

Geometry and internuclear transport

The choice of our geometric parameters, in particular of the lattice constant, is inspired by experimental
measurements in the Drosophila embryo by Gregor et al. [1]. Information on the diffusion constants of
proteins involved in early Drosophila patterning is scarce. The diffusion constant of the morphogen Bcd
has been measured by several groups, yet its true value is still under debate [12, 102]. In our model
we therefore set for all patterning proteins an effective internuclear diffusion constant DP = 1 µm2 /s,
which comprises both protein import/export and actual diffusion. This value is a reasonable cytoplasmic
diffusion coefficient and well within the bounds reported for Bcd.

24
 The simulated lattice is 40 nuclei long so that the total system length L roughly corresponds to the
posterior 2/3 of the Drosophila embryo in cycle 14. To reduce computation effort we simulate a system
with smaller circumference (8 nuclei) as compared to the living embryo. This is justified by the fact that
for our standard
p diffusion constant DP and effective protein lifetime µeff the diffusive correlation length
lcorr = DP /µeff is ≤ 2 nuclei. A larger circumference therefore is not expected to introduce new
features into the system, but might alter the timescales of expression boundary movement and domain
desintegration. We discuss the effect of reduced system size on measured stability times in the “Discussion”
section.
A complete overview of the specific numerical values of our model parameters is found in Table S1 of
the Supporting Information.

Data analysis
Quantification of pattern stability
In order to analyse pattern stability we represent each simulated pattern as a point in (λAC , λBD ) phase
space. For every pattern simulation from time t = 0 until time t = tend the temporal sequence of these
points corresponds to a trajectory in the (λAC , λBD ) space. For each parameter choice and pinning
scenario, we restarted the simulations with 6000 trajectories started from the relaxed initial patterns at
t = 0; the trajectories ensemble is then further enriched by the branching process at the NS-FFS interfaces.
Next, the trajectories are binned with the statistical weight assigned by NS-FFS, and then the histograms
are normalized. As a result, we can identify a few distinct regions that accumulate probability.
In order to formally define these regions we define rectangular boundaries that enclose accumulated
probability regions corresponding to different types of patterns:
• the metastable main basin with five-stripe pattern:
RS ≡ {(λAC , λBD )|λAC ≤ 0.45 ∧ λBD ≤ 0.43}
• the basin in which either the A or C protein domain was lost:
†
RAC ≡ {(λAC , λBD )|λAC > 0.45 ∧ λBD ≤ 0.43}
• the basin in which either the B or D protein domain was lost:
†
RBD ≡ {(λAC , λBD )|λAC ≤ 0.45 ∧ λBD > 0.43}
• the basin in which either A or C and one of B or D were lost:
R‡ ≡ {(λAC , λBD )|λAC > 0.45 ∧ λBD > 0.43}
Note that the location of the regions slightly changes for different values of κ. We found that the
above boundary definitions constitute a good compromise. For each basin we compute the fraction of
total probability as a function of time by integrating the weights of trajectories that are within the basin
at time t. We define the RR pattern survival probability to be the integrated probability in RS at time t after
initialization: S(t) = RS p(t)dλAC dλBD . As expected, S(t) displays roughly exponential decay behavior
after a certain lag phase that can be attributed to initial relaxation. To obtain the pattern destruction
rate kD we fit a function f (x) ≡ exp(−kD (t − tlag )) to S(t). This only yields satisfactory results if
the fitting range is adapted accordingly, i.e. only S(t) values for t > tlag are taken into account. Since
tlag itself is a fitting parameter we adopted the following protocol: Starting from a value of tstart that is
clearly in the relaxation regime we perform the fit on the interval [tstart , tend ] where tend is the largest
time recorded. We then choose the fitted values kD and tlag for which |tlag − tstart | is minimal. From
this we compute the pattern stability time (average time until pattern has lost one of the domains) via
τD ≡ 1/kD . In most considered cases the patterns are very stable, i.e. kD very small, and we can expand
S(t) ≃ 1 − kD (t − tlag ). As a control, we therefore also fitted g(t) ≡ kD (t − tlag ) to 1 − S(t) for a fixed
tlag clearly in the exponential regime and obtained almost identical results.

Computation of average probability fluxes

To quantify which destruction pathways are dominant we computed the average fluxes Javg into the
regions of (partly) destroyed patterns. Here the average flux is defined as the average rate of increase in
time of the fractional probability in the region and obtained by fitting a linear function h(t) ≡ Javg t + P0
RR † †
to PR (t) ≡ R p(t)dλAC dλBD for R ∈ {RAC , RBD , R‡ } over the interval [tstart , tend ] with tstart chosen
such that ∂t PR (t) ̸= 0 for t > tstart . P0 depends on the particular choice of tstart and is discarded.

25
Computation of average flux velocities
The average local drift velocity and diffusion constant of the trajectories in the (λAC , λBD ) phase space
are computed by averaging displacements ∆λAC(BD) ≡ λAC(BD) (t + ∆t) − λAC(BD) (t) and squared
displacements ∆λ2 ≡ ∆λ2AC + ∆λ2BD on a two-dimensional lattice of bins covering the whole phase
space. Displacements ∆λAC(BD) are assigned to the bin at ⃗λ ≡ (λAC , λBD ), i.e. we are averaging
outgoing displacements and the averaged vector ⟨∆λ⟩( ⃗ ⃗λ) therefore will represent the average velocity
withhwhich trajectories
 leave this bin. The localiphase space diffusion constant is calculated as Dλ (⃗λ) ≡
1
4∆t ⟨∆λ ⟩(λ) − ⟨∆λAC ⟩2 (⃗λ) + ⟨∆λAC ⟩2 (⃗λ) . This is done in the same way for other combinations
2 ⃗

of phase space coordinates. The diffusion-drift decomposition is explained in more detail in the Supporting
Information.

Computation of “pseudopotential” landscapes

The trajectory binning procedure used for computing the average flux velocities as described above was
at the same time used for computing the “pseudopotential” − log(p̃(⃗λ). Herein p̃(⃗λ) is the local density
calculated from the reweighed number of trajectories leaving the bin at ⃗λ = (λAC , λBD ), and smoothened
afterwards by 2D median filtering over nfilt neighboring bins. For the 2D median filtering we used the
medfilt2 function from the MATLAB Image Processing Toolbox. We empirically chose nfilt = 4 as we
found that this choice efficiently removes local spikes in p(⃗λ) without changing the overall shape of the
landscape.

Perturbation experiments
Simulations starting from perturbed initial conditions were performed directly via the GGG simulator.
First the systems were relaxed to representative states within the metastable basin for a simulated time
of trelax = 30 min. The final states of these runs then were post-modified according to the following two
protocols:
1. “C expansion”: starting from mid-embryo the central C protein domain was expanded as follows:
the configurations in the nuclei just posterior to mid-embryo were copied and used to overwrite
configurations in the subsequent ∆ rows in the axial (z-) direction of the cylinder. The original
configurations were stored and for each nucleus at row zi > Nz /2 + ∆ (counting from the anterior)
the configuration was overwritten by the original configuration at zi −∆. The posterior-most nucleus
was exempted from overwriting to preserve pinning.
2. “A expansion”: here the anterior A protein domain was enlarged at the expense of the C protein
domain. To this purpose we applied the same copy-paste procedure as above starting form zi = 5,
however only nuclei up to mid-embryo (zi ≤ Nz /2) were overwritten by the original configurations
at zi − ∆.
∆ quantifies the severity of perturbation. We found that ∆ < 4 results in changes to the pattern that were
hard to distinguish from noise, while for ∆ > 12 perturbations were large enough to induce immediate
pattern destruction with high probability. We therefore limited systematic tests to perturbations with
∆ ∈ {4, 8, 12}. Starting from the perturbed initial conditions simulations were continued for tsim = 20 h
and snapshots of the current configurations in all nuclei were written out with an acquisition interval of
10 min (simulated time). 10 samples starting from 10 different perturbed initial conditions were produced
for each set of parameters.
In order to overcome the difficulties of boundary detection we quantified the motion of protein domains
by tracking their center of mass (CoM) along the z-axis of the cylinder. For each considered gene G we
define the CoM zG as
R R
zGtot (r, z)drdz
zG ≡ Rz Rr (19)
G (r, z)drdz
z r tot

where Gtot = [G] + 2[G2 ] is the total copy number. Since our system features two A domains we
calculate zA separately for the anterior (Aant ) and the posterior (Apost ) part of the embryo by restricting
z-integration adequately. While the CoM remains unchanged upon symmetric changes of the domain

26
boundaries or global copy number increase, it is well-suited to indicate relaxations from the asymmetric
perturbations that we apply. To find general trends in the time-evolution of the domains CoM trajectories
were averaged over the 10 samples.

Effective reaction-diffusion dynamics of dimer expression

In this section we map the fully microscopic model defined in Fig. 1 onto the effective model from [55].
First, we postulate that the stochastic dynamics of gene expression studied in this paper corresponds to
the following effective dynamical equations:
D
∂t X1 (x, t) = −µM X1 (x, t) − kon X12 (x, t) + koff
D
X2 (x, t) + βf ({Y2 (x, t)}Y̸=X ) ,
D D
(20)
∂t X2 (x, t) = D0 ∂xx X2 (x, t) − (koff + µD )X2 (x, t) + kon X12 (x, t) ,

where X, Y ∈ {A, B, C, D} denotes the expressed protein species, X1 , Y1 are the concentrations of its
monomer, and X2 , Y2 are concentrations of its dimers. The synthesis and decay of dimers is described by
rates kon and koff . Both, monomers and dimers degrade with rates µD and µM . In this system, only dimers
are allowed to diffuse (with diffusivity D0 ) and only monomers are primarily synthesized, with maximal
production rate β and production kinetics described by function f ({Y2 (x, t)}Y̸=X ), which we specify later.
However, we assume that in the absence of other dimers Y2 ̸= X2 the production is active by default, so
f ({0}) = 1. Since the system has cylindrical symmetry, we will treat the axis x as distinguished and treat
the system as effectively one-dimensional.
The fact that the model defined by eqs. (20) involves monomers and dimers complicates its mapping
onto the model in [55]. We therefore translate it into a simplified model, tracking the effective dynamics
of dimers only. To this end, we will first determine the ratio between stationary concentrations X̃1 and X̃2
in the absence of other dimers (Y2 (x, t) = 0) and assuming system homogeneity. In this case, equations
(20) turn into:
D
0 = −µM X̃1 − kon X̃12 + koff
D
X̃2 + β ,
(21)
0 = −(koff + µD )X̃2 + kon X̃12 .
D D

Solving for X̃1 and X̃2 , we obtain:

s !
1 k D + µD D + µ )2
(koff D k D + µD
X̃1 = D − off + 2
D β off
µ2M + 4kon ,
2kon µD µD µD
(22)
β µM
X̃2 = − X̃1 .
µD µD

We will now sum both equations in (20) to obtain

∂t (X1 (x, t) + X2 (x, t)) = D∂xx X2 (x, t) − µD X2 (x, t) − µM X1 (x, t) + βf ({Y2 }Y̸=X ) , (23)

and approximate
X̃1
X1 (x, t) ≈ X2 (x, t) . (24)
X̃2
In other words, we assume that X1 (x, t) follows strictly X2 (x, t). The advantage of this approximation
is that it becomes exact in the stationary state. This procedure results in the following effective equation
for X2 (x, t):
∂t X2 (x, t) = DX ∂xx X2 (x, t) − γX X2 (x, t) + HX f ({Y2 }Y̸=X ) , (25)
where the rescaled constants are:

D0 µD + µM X̃
X̃
1
β
2
DX = , γX = , HX = . (26)
X̃1 X̃1 X̃1
1+ X̃2
1+ X̃2
1+ X̃2

We can now specify the kinetics function. The microscopic dynamics is such that each gene X is
produced, unless it is blocked by the biding of any other dimer to its repressor site on the promoter. In
the averaged-out description, we expect that a sufficiently high concentration of free repressor particles
effectively shuts down the production of X. Similarly to [55], we will assume that this transition is steep,

27
so we can choose the functional form of the regulatory Hill function in (25) to have the overall shape of
Heaviside step function:  
X
f ({Y2 }Y̸=X ) = θ  ϵXY Y2 (x, t) − CX  . (27)
Y̸=X

Finally, we relate the effective gene interaction constants ϵXY to microscopic parameters by the fol-
lowing reasoning: In the microscopic simulations, the attachment of Y2 to the repressor site is described
R
by the constant kon and detachment by ksoff or kwoff
. Assuming that the repressor production speed can be
approximated by Michealis-Menten kinetics, with the repressor site acting like a “catalyst”, we know that
off R
KY = kw,s /kon (28)
off off
where kw,s (standing for either kw or ksoff ) is the concentration of repressor dimers Y2 at which the velocity
of production of Y1 is at the half of its maximal value. We postulate that at this point Y2 effectively
switches off the production of X, and we equate this point with reaching the threshold for production in
(27). Hence, the following is satisfied:

ϵXY KY − CX = 0 . (29)

Solving for ϵXY we obtain:

CX
ϵXY = . (30)
KY
We choose CX < 0 to ensure that the production of X is active by default, in the absence of repressive
dimers (Y2 (x, t) = 0). Since CX is now present in every term in (27), we can factor it out and neglect.
Taken together, and assuming that diffusion, degradation and production constants are the same for all
genes, that is: DX = D, γX = γ and HX = H for all X ∈ {A, B, C, D}; the microscopic dynamics of
gene expression mapped onto the effective model results in eq. (5).

Derivation of stability conditions for a contact zone between two domains

This section outlines the origin of stability conditions for the effective continuous model, derived in the
previous section, which is employed for the analysis of the four-gene pattern in the approximation of
separate contact zones. We discuss the major steps leading to stability conditions in the current case,
taking advantage of specific setting of the system studied in this work. For a detailed and more general
derivation we refer the reader to our earlier work [55].
Let us consider a pair of genes (X,Y), whose expression dynamics is described by the effective equations:
 
∂t X2 (x, t) = D ∂xx X2 (x, t) − γX2 (x, t) + H θ C̃X − ϵXY Y2 (x, t)
  (31)
∂t Y2 (x, t) = D ∂xx Y2 (x, t) − γY2 (x, t) + H θ C̃Y − ϵYX X2 (x, t)

These equations describe both type-(i) and type-(ii) contact zones, albeit for different values of constants
(see Eqs. (9) and (10)). We assume that the system is open (L → +∞) and two respective expression
domains occupy the opposite ‘ends’ of the system. That is, initially, the expression profiles of respective
dimers read:

X2 (x, 0) = AX θ(qX (0) − x) Y2 (x, 0) = AY θ(x − qY (0)) (32)

Here AX and AY are initial amplitudes, sufficiently high to initiate auto-activation, while qX (t) and qY (t)
denote the positions of expression domain boundaries.
The effective equations in the form (5) can be solved analytically in this system, providing the spatio-
temporal profile of expression for both dimers:
Z +∞ Z t Z +∞
X2 (x, t) = dx′ G(x − x′ , 0)X2 (x, 0) + H dt′ dx′ G(x − x′ , t − t′ )
−∞ 0 qX (t′ )
Z Z Z (33)
+∞ t qY (t′ )
′ ′ ′ ′ ′ ′
Y2 (x, t) = dx G(x − x , 0)Y2 (x, 0) + H dt dx G(x − x , t − t )
−∞ 0 −∞

28
where the Green’s function of eq. (5) in the open system reads:
(x−x′ )2
−γt− 4D(t−t′ )
e
G(x − x′ , t − t′ ) = p (34)
4πD(t − t′ )

However, the solution (33) is known only up to the position of domain boundaries, qX (t) and qY (t). In
order to determine these positions, the ”Free Boundary Problem” must be solved, defined by the activation
conditions at each boundary:
C̃X = ϵXY Y2 (qX (t), t)
(35)
C̃Y = ϵYX X2 (qY (t), t)
Inserting the solution (33) into these equations leads to a system of coupled nonlinear integral equations.
For t large enough that the system ‘forgets’ its initial conditions, these equations simplify into:
Z t Z qY (t′ )
C̃X
= dt′ dx′ G(qX (t) − x′ , t − t′ )
ϵXY H 0 −∞
Z t Z +∞ (36)
C̃Y ′ ′ ′ ′
= dt dx G(qY (t) − x , t − t )
ϵYX H 0 qX (t′ )

In [55] we show that the asymptotic solution of these equations is provided by the constant velocity
∞ ∞ ∞ ∞
ansatz, qX (t) = vX t + qX and qY (t) = vY t + qY , where qX and qY are constants. For this choice,
the right-hand side integrals in Eqs. (36) saturate at constant values, though corresponding to domain
boundaries travelling with constant velocities. Moreover, for sufficiently strong interactions between the
genes, these velocities must be equal, vX = vY = v, meaning that the domains change their size in a
coordinated manner. For the constant velocity ansatz, the integrals can be analytically computed. This
turns the system of integral equations system into an algebraic system, defining the common velocity v
∞ ∞
and the distance between the boundaries ∆r = qX − qY . In the t → +∞ limit, this system reads:
√
!
2C̃X γ v∆r
−
|∆r| 4Dγ+v 2 v
= sgn(∆r) − e D 2D p + sgn(∆r)
ϵXY H 4Dγ + v 2
√
! (37)
2C̃Y γ v∆r
−
|∆r| 4Dγ+v 2 v
= sgn(∆r) + e D 2D p − sgn(∆r)
ϵYX H 4Dγ + v 2

Finally, the patterns for which v = 0 are stable, that is, they do not change in the long-time limit.
Substituting v = 0 turns Eqs. (37) into:
 |∆r|

RX = sgn(∆r) 1 − e− λ
 |∆r|
 (38)
RY = sgn(∆r) 1 − e− λ

where RX and RY are defined as in Eqs. (11). In order to make v = 0 the solution of Eqs. (37), Eqs.
(38) must be solvable and satisfied by the same ∆r, as the system is over-defined. This happens, provided
that:

RX = RY −1 < RX < 1 −1 < RX < 1 (39)

which constitutes the stability conditions utilized in this work. Further, we can also derive:
 
RX
|∆r| = −λ ln 1 − (40)
sgn(∆r)

However, one can notice in Eqs. (38) that as 1 − e−|∆r|/λ > 0 for any ∆r, then:

sgn(∆r) = sgn(RX ) = sgn(RY ) (41)

This allows us to obtain ∆r, as provided by the formula (13).

29
Acknowledgments
M.M. and M.Z. were supported by the Polish National Agency for Academic Exchange and by a grant from
the Priority Research Area DigiWorld under the Strategic Programme Excellence Initiative at Jagiellonian
University. M.Z. was supported by National Science Center, Poland, no. 2021/42/E/NZ2/00188. T.R.S.,
N.B.B. and P.R.t.W. were supported by the Foundation for Fundamental Research on Matter (FOM),
which is part of the Netherlands Organisation for Scientific Research (NWO). T.R.S. was supported by
the Center for Multiscale Modelling in Life Sciences (CMMS), sponsored by the Hessian Ministry of Science
and Art. The funding agencies did not play any role in the study design, data collection and analysis,
decision to publish, or preparation of the manuscript.

References
[1] Gregor T, Wieschaus EF, McGregor AP, Bialek W, Tank DW. Stability and Nuclear Dynamics of
the Bicoid Morphogen Gradient. Cell. 2007;130(1):141–152. doi:10.1016/[Link].2007.05.026.
[2] Rogers KW, Schier AF. Morphogen gradients: from generation to interpretation. Annual Review of
Cell and Developmental Biology. 2011;27:377–407. doi:10.1146/annurev-cellbio-092910-154148.

[3] Lander AD. Pattern, growth, and control. Cell. 2011;144(6):955–969.

doi:10.1016/[Link].2011.03.009.
[4] Balaskas N, Ribeiro A, Panovska J, Dessaud E, Sasai N, Page K, et al. Gene Regulatory
Logic for Reading the Sonic Hedgehog Signaling Gradient in the Vertebrate Neural Tube. Cell.
2012;148(1):273–284. doi:10.1016/[Link].2011.10.047.
[5] Kicheva A, Cohen M, Briscoe J. Developmental pattern formation: insights from physics and biology.
Science. 2012;338(6104):210–212. doi:10.1126/science.1225182.
[6] Shvartsman SY, Baker RE. Mathematical models of morphogen gradients and their effects on
gene expression. Wiley Interdisciplinary Reviews: Developmental Biology. 2012;1(5):715–730.
doi:10.1002/wdev.55.
[7] Briscoe J, Small S. Morphogen rules: design principles of gradient-mediated embryo patterning.
Development. 2015;142(23):3996–4009. doi:10.1242/dev.129452.
[8] Bier E, De Robertis EM. BMP gradients: A paradigm for morphogen-mediated developmental
patterning. Science. 2015;348(6242):aaa5838. doi:10.1126/science.aaa5838.
[9] Stapornwongkul KS, de Gennes M, Cocconi L, Salbreux G, Vincent JP. Patterning and
growth control in vivo by an engineered GFP gradient. Science. 2020;370(6514):321–327.
doi:10.1126/science.abb8205.
[10] Stapornwongkul KS, Vincent JP. Generation of extracellular morphogen gradients: the case for
diffusion. Nature Reviews Genetics. 2021;22:393–411. doi:10.1038/s41576-021-00342-y.
[11] Simsek MF, Özbudak EM. Patterning principles of morphogen gradients. Open Biology. 2022;12(10).
doi:10.1098/rsob.220224.
[12] Gregor T, Tank DW, Wieschaus EF, Bialek W. Probing the Limits to Positional Information. Cell.
2007;130(1):153–164. doi:10.1016/[Link].2007.05.025.
[13] Dubuis JO, Tkacik G, Wieschaus EF, Gregor T, Bialek W. Positional information, in bits. Proceed-
ings of the National Academy of Sciences of the United States of America. 2011;110(41):16301–
16308. doi:10.1073/pnas.1315642110.
[14] Sokolowski TR, Erdmann T, ten Wolde PR. Mutual Repression Enhances the Steepness
and Precision of Gene Expression Boundaries. PLoS Computational Biology. 2012;8:e1002654.
doi:10.1371/[Link].1002654.

30
[15] Raspopovic J, Marcon L, Russo L, Sharpe J. Digit patterning is controlled by a Bmp-Sox9-
Wnt Turing network modulated by morphogen gradients. Science. 2014;345(6196):566–570.
doi:10.1126/science.1252960.
[16] Zagorski M, Tabata Y, Brandenberg N, Lutolf MP, Tkačik G, Bollenbach T, et al. Decod-
ing of position in the developing neural tube from antiparallel morphogen gradients. Science.
2017;356(6345):1379–1383. doi:10.1126/science.aam5887.
[17] Petkova MD, Tkačik G, Bialek W, Wieschaus EF, Gregor T. Optimal Decoding of Cellular Identities
in a Genetic Network. Cell. 2019;176(4):844–855.e15. doi:10.1016/[Link].2019.01.007.
[18] Kuzmicz-Kowalska K, Kicheva A. Regulation of size and scale in vertebrate spinal cord development.
WIREs Developmental Biology. 2020;10(3):e383. doi:[Link]
[19] Tkačik G, Gregor T. The many bits of positional information. Development. 2021;148(2):dev176065.
doi:10.1242/dev.176065.
[20] Exelby K, Herrera-Delgado E, Perez LG, Perez-Carrasco R, Sagner A, Metzis V, et al. Precision of
tissue patterning is controlled by dynamical properties of gene regulatory networks. Development.
2021;148(4):dev197566. doi:10.1242/dev.197566.
[21] Iyer KS, Prabhakara C, Mayor S, Rao M. Cellular compartmentalisation and receptor promiscuity as a
strategy for accurate and robust inference of position during morphogenesis. eLife. 2023;12:e79257.
doi:10.7554/eLife.79257.
[22] Sokolowski TR, Gregor T, Bialek W, Tkačik G. Deriving a genetic regulatory network from an
optimization principle. arXiv:230205680. 2023;doi:10.48550/arxiv.2302.05680.
[23] Minchington T, Rus S, Kicheva A. Control of tissue dimensions in the develop-
ing neural tube and somites. Current Opinion in Systems Biology. 2023; p. 100459.
doi:[Link]
[24] Swain PS, Elowitz MB, Siggia ED. Intrinsic and extrinsic contributions to stochasticity in gene
expression. Proceedings of the National Academy of Sciences of the United States of America.
2002;99(20):12795–800. doi:10.1073/pnas.162041399.
[25] Raser JM, O’Shea EK. Noise in Gene Expression: Origins, Consequences, and Control. Science.
2005;309:2010–2014. doi:10.1126/science.1105891.
[26] Paulsson J. Models of stochastic gene expression. Physics of Life Reviews. 2005;2(2):157–175.
doi:10.1016/[Link].2005.03.003.
[27] Raj A, van Oudenaarden A. Nature, Nurture, or Chance: Stochastic Gene Expression and Its
Consequences. Cell. 2008;135(2):216–226. doi:10.1016/[Link].2008.09.050.
[28] Chalancon G, Ravarani CNJ, Balaji S, Martinez-Arias A, Aravind L, Jothi R, et al. Interplay between
gene expression noise and regulatory network architecture. Trends in Genetics. 2012;28(5):221–232.
doi:10.1016/[Link].2012.01.006.
[29] Averbukh I, Gavish A, Shilo BZ, Barkai N. Dealing with noise: The challenge of buffering biological
variability. Current Opinion in Systems Biology. 2017;1:69–74. doi:10.1016/[Link].2016.12.011.
[30] Perez-Carrasco R, Guerrero P, Briscoe J, Page KM. Intrinsic Noise Profoundly Alters the Dynamics
and Steady State of Morphogen-Controlled Bistable Genetic Switches. PLoS Computational Biology.
2016;12(10):e1005154. doi:10.1371/[Link].1005154.
[31] Kicheva A, Bollenbach T, Ribeiro A, Valle HP, Lovell-Badge R, Episkopou V, et al. Co-
ordination of progenitor specification and growth in mouse and chick spinal cord. Science.
2014;345(6204):1254927. doi:10.1126/science.1254927.
[32] Almuedo-Castillo M, Bläßle A, Mörsdorf D, Marcon L, Soh GH, Rogers KW, et al. Scale-invariant
patterning by size-dependent inhibition of Nodal signalling. Nature Cell Biology. 2018;20(9):1032–
1042. doi:10.1038/s41556-018-0155-7.

31
[33] Cohen M, Kicheva A, Ribeiro A, Blassberg R, Page KM, Barnes CP, et al. Ptch1 and Gli reg-
ulate Shh signalling dynamics via multiple mechanisms. Nature Communications. 2015;6:1–12.
doi:10.1038/ncomms7709.

[34] Drocco JA, Grimm O, Tank DW, Wieschaus E. Measurement and Perturbation of Mor-
phogen Lifetime: Effects on Gradient Shape. Biophyiscal Journal. 2011;101(8):1807–1815.
doi:10.1016/[Link].2011.07.025.
[35] Durrieu L, Kirrmaier D, Schneidt T, Kats I, Raghavan S, Hufnagel L, et al. Bicoid gradient for-
mation mechanism and dynamics revealed by protein lifetime analysis. Molecular Systems Biology.
2018;14(9):1–15. doi:10.15252/msb.20188355.
[36] Huh D, Paulsson J. Non-genetic heterogeneity from stochastic partitioning at cell division. Nature
Genetics. 2011;43(2):95–100. doi:10.1038/ng.729.
[37] Thomas P. Intrinsic and extrinsic noise of gene expression in lineage trees. Scientific Reports.
2019;9(1):1–16. doi:10.1038/s41598-018-35927-x.
[38] Pedraza JM, Oudenaarden AV. Noise Propagation in Gene Networks. Science. 2005;307(5717):1965–
1969. doi:10.1126/science.1109090.
[39] Stockholm D, Edom-Vovard F, Coutant S, Sanatine P, Yamagata Y, Corre G, et al. Bistable Cell
Fate Specification as a Result of Stochastic Fluctuations and Collective Spatial Cell Behaviour.
PLOS ONE. 2011;5(12):1–12. doi:10.1371/[Link].0014441.
[40] Waddington CH. Canalization of Development and the Inheritance of Acquired Characters. Nature.
1942;150:563–565. doi:doi:10.1038/150563a0.
[41] Waddington CH. Canalization of Development and Genetic Assimilation of Acquired Characters.
Nature. 1959;183(4676):1654–1655. doi:10.1038/1831654a0.
[42] Alon U. Network motifs: theory and experimental approaches. Nature Reviews Genetics.
2007;8(6):450–61. doi:10.1038/nrg2102.
[43] Vakulenko S, Manu, Reinitz J, Radulescu O. Size Regulation in the Segmentation of Drosophila:
Interacting Interfaces between Localized Domains of Gene Expression Ensure Robust Spatial Pat-
terning. Physical Review Letters. 2009;103(16):168102. doi:10.1103/PhysRevLett.103.168102.
[44] Cotterell J, Sharpe J. An atlas of gene regulatory networks reveals multiple three-gene
mechanisms for interpreting morphogen gradients. Molecular Systems Biology. 2010;6:425.
doi:10.1038/msb.2010.74.

[45] Burda Z, Krzywicki A, Martin OC, Zagorski M. Motifs emerge from function in model gene regulatory
networks. Proceedings of the National Academy of Sciences of the United States of America.
2011;108(42):17263–17268. doi:10.1073/pnas.1109435108.
[46] Verd B, Crombach A, Jaeger J. Dynamic Maternal Gradients Control Timing and Shift-
Rates for Drosophila Gap Gene Expression. PLoS Computational Biology. 2017;13(2):e1005285.
doi:10.1371/[Link].1005285.
[47] Verd B, Monk NAM, Jaeger J. Modularity, criticality, and evolvability of a developmental gene
regulatory network. eLife. 2019;8:e42832. doi:10.7554/eLife.42832.
[48] Sokolowski TR, Tkačik G. Optimizing information flow in small genetic networks. IV. Spatial cou-
pling. Physical Review E. 2015;91(6):062710. doi:10.1103/PhysRevE.91.062710.

[49] Jaeger J, Blagov M, Kosman D, Kozlov KN, Manu, Myasnikova E, et al. Dynamical Analysis of Regu-
latory Interactions in the Gap Gene System of Drosophila melanogaster. Genetics. 2004;167(4):1721–
1737. doi:10.1534/genetics.104.027334.
[50] Jaeger J. The gap gene network. Cellular and Molecular Life Sciences. 2011;68(2):243–274.
doi:10.1007/s00018-010-0536-y.

32
[51] Dubuis JO, Samanta R, Gregor T. Accurate measurements of dynamics and reproducibility in small
genetic networks. Molecular Systems Biology. 2013;9:639.
[52] Manu, Surkova S, Spirov AV, Gursky VV, Janssens H, Kim AR, et al. Canalization of Gene Expression
in the Drosophila Blastoderm by Gap Gene Cross Regulation. PLoS Biology. 2009;7(3):591–603.
doi:10.1371/[Link].1000049.
[53] Surkova S, Kosman D, Kozlov K, Manu, Myasnikova E, Samsonova AA, et al. Characterization of
the Drosophila segment determination morphome. Developmental Biology. 2008;313(2):844–862.
doi:10.1016/[Link].2007.10.037.
[54] Clyde DE, Corado MSG, Wu X, Pare A, Papatsenko D, Small S. A self-organizing system of
repressor gradients establishes segmental complexity in Drosophila. Nature. 2003;426(6968):849–
853. doi:10.1038/nature02189.
[55] Majka M, Ho RDJG, Zagorski M. Stability of pattern formation in systems with dynamic source
regions. Physical Review Letters. 2023;130:098402. doi:10.1103/PhysRevLett.130.098402.
[56] Becker NB, ten Wolde PR. Rare switching events in non-stationary systems. Journal of Chemical
Physics. 2012;136(17):174119–174119–15. doi:10.1063/1.4704812.
[57] Becker NB, Allen RJ, ten Wolde PR. Non-stationary forward flux sampling. Journal of Chemical
Physics. 2012;136(17):174118–174118–18. doi:10.1063/1.4704810.
[58] Manu, Surkova S, Spirov AV, Gursky VV, Janssens H, Kim AR, et al. Canalization of Gene Expression
and Domain Shifts in the Drosophila Blastoderm by Dynamical Attractors. PLoS Computational
Biology. 2009;5(3):e1000303. doi:10.1371/[Link].1000303.
[59] Erdmann T, Howard M, ten Wolde PR. Role of Spatial Averaging in the Precision of Gene Expression
Patterns. Physical Review Letters. 2009;103(25):258101. doi:10.1103/PhysRevLett.103.258101.
[60] Kraut R, Levine M. Mutually repressive interactions between the gap genes giant and Kruppel define
middle body regions of the Drosophila embryo. Development. 1991;111(2):611–621.
[61] Kraut R, Levine H. Spatial regulation of the gap gene giant during Drosophila development. De-
velopment. 1991;111(2):601–609.
[62] Eldon ED, Pirrotta V. Interactions of the Drosophila gap gene giant with maternal and zygotic
pattern-forming genes. Development. 1991;111(2):367–378.
[63] Hülskamp M, Pfeifle C, Tautz D. A morphogenetic gradient of hunchback protein organizes
the expression of the gap genes Kruppel and knirps in the early Drosophila embryo. Nature.
1990;346(6284):577–580. doi:10.1038/346577a0.
[64] Jäckle H, Tautz D, Schuh R, Seifert E, Lehmann R. Cross-regulatory interactions among the gap
genes of Drosophila. Nature. 1986;324(6098):668–670. doi:10.1038/324668a0.
[65] Driever W, Nüsslein-Volhard C. The bicoid protein is a positive regulator of hunchback transcription
in the early Drosophila embryo. Nature. 1989;337:138–143.
[66] Porcher A, Abu-Arish A, Huart S, Roelens B, Fradin C, Dostatni N. The time to measure po-
sitional information: maternal Hunchback is required for the synchrony of the Bicoid transcrip-
tional response at the onset of zygotic transcription. Development. 2010;137(16):2795–2804.
doi:10.1242/dev.051300.
[67] Margolis JS, Borowsky ML, Steingrı́msson E, Shim CW, Lengyel JA, Posakony JW. Posterior stripe
expression of hunchback is driven from two promoters by a common enhancer element. Development.
1995;121:3067–3077.
[68] Casanova J. Pattern formation under the control of the terminal system in the Drosophila embryo.
Development. 1990;110:621–628.
[69] Weigel D, Jürgens G, Klingler M, Jäckle H. Two gap genes mediate maternal terminal pattern
information in Drosophila. Science. 1990;27(4954):495–498.

33
[70] Zwanzig R. Memory Effects in Irreversible Thermodynamics. Physical Review. 1961;124:983–992.
doi:10.1103/PhysRev.124.983.
[71] Mori H. Transport, Collective Motion, and Brownian Motion. Progress of Theoretical Physics.
1965;33(3):423–455. doi:10.1143/PTP.33.423.
[72] Yang S, Onuchic JN, Levine H. Effective stochastic dynamics on a protein folding energy landscape.
Journal of Chemical Physics. 2006;125(5):054910. doi:10.1063/1.2229206.
[73] Kopelevich DI, Panagiotopoulos AZ, Kevrekidis IG. Coarse-grained kinetic computations for rare
events: Application to micelle formation. Journal of Chemical Physics. 2005;122:044908.
[74] Hummer G, Kevrekidis IG. Coarse molecular dynamics of a peptide fragment: Free energy, kinetics,
and long-time dynamics computations. Journal of Chemical Physics. 2003;118:10762.
[75] Plotkin SS, Wolynes PG. Non-Markovian Configurational Diffusion and Reaction Coordinates for
Protein Folding. Physical Review Letters. 1998;80:5015–5018. doi:10.1103/PhysRevLett.80.5015.

[76] Seyboldt R, Lavoie J, Henry A, Vanaret J, Petkova MD, Gregor T, et al. Latent space of a small
genetic network: Geometry of dynamics and information. Proceedings of the National Academy of
Sciences. 2022;119(26):e2113651119. doi:10.1073/pnas.2113651119.
[77] Driever W, Nüsslein-Volhard C. A gradient of bicoid protein in Drosophila embryos. Cell.
1988;54(1):83–93.
[78] Driever W, Nüsslein-Volhard C. The bicoid protein determines position in the Drosophila embryo in
a concentration-dependent manner. Cell. 1988;54(1):95–104. doi:10.1016/0092-8674(88)90183-3.
[79] Frohnhöfer HG, Nüsslein-Volhard C. Organization of anterior pattern in the Drosophila embryo by
the maternal gene bicoid. Nature. 1986;324(6093):120–125. doi:10.1038/324120a0.

[80] Struhl G, Struhl K, Macdonald PM. The gradient morphogen bicoid is a concentration-dependent
transcriptional activator. Cell. 1989;57(7):1259–1273.
[81] Macdonald PM, Struhl G. A molecular gradient in early Drosophila embryos and its role in specifying
the body pattern. Nature. 1986;324(6097):537–545. doi:10.1038/324537a0.

[82] Schulz C, Tautz D. Zygotic caudal regulation by hunchback and its role in abdominal segment
formation of the Drosophila embryo. Development. 1995;121(4):1023–1028.
[83] Mlodzik M, Gehring WJ. Expression of the caudal gene in the germ line of Drosophila: Formation of
an RNA and protein gradient during early embryogenesis. Cell. 1987;48(3):465–478. doi:16/0092-
8674(87)90197-8.

[84] Kozlov K, Surkova S, Myasnikova E, Reinitz J, Samsonova M. Modeling of Gap Gene

Expression in Drosophila Kruppel Mutants. PLoS Computational Biology. 2012;8:e1002635.
doi:10.1371/[Link].1002635.
[85] Crombach A, Wotton KR, Cicin-Sain D, Ashyraliyev M, Jaeger J. Efficient Reverse-Engineering
of a Developmental Gene Regulatory Network. PLoS Computational Biology. 2012;8(7):e1002589.
doi:doi:10.1371/[Link].1002589.
[86] Surdej P, Jacobs-Lorena M. Developmental Regulation of bicoid mRNA Stability Is Mediated by the
First 43 Nucleotides of the 3′ Untranslated Region. Molecular and Cellular Biology. 1998;18(5):2892–
2900. doi:10.1128/mcb.18.5.2892.

[87] Spirov A, Fahmy K, Schneider M, Frei E, Noll M, Baumgartner S. Formation of the bicoid morphogen
gradient: an mRNA gradient dictates the protein gradient. Development. 2009;136(4):605–614.
doi:10.1242/dev.031195.
[88] Nikolić M, Antonetti V, Liu F, Muhaxheri G, Petkova MD, Scheeler M, et al. Scale invariance in
early embryonic development. arXiv:231217684. 2023;doi:10.48550/arXiv.2312.17684.

34
 [89] Morales JS, Raspopovic J, Marcon L. From embryos to embryoids: How external signals
and self-organization drive embryonic development. Stem Cell Reports. 2021;16(5):1039–1050.
doi:10.1016/[Link].2021.03.026.

[90] Jaeger J, Irons D, Monk N. Regulative feedback in pattern formation: towards a general relativistic
theory of positional information. Development. 2008;135(19):3175–3183. doi:10.1242/dev.018697.
[91] El-Sherif E, Levine M. Shadow Enhancers Mediate Dynamic Shifts of Gap Gene Expression in the
Drosophila Embryo. Current Biology. 2016;26(9):1164–1169. doi:10.1016/[Link].2016.02.054.

[92] McClay DR. Evolutionary crossroads in developmental biology: sea urchins. Development.
2011;138(13):2639–2648. doi:10.1242/dev.048967.
[93] McCarty AS, Kleiger G, Eisenberg D, Smale ST. Selective Dimerization of a C2H2 Zinc Finger
Subfamily. Molecular Cell. 2003;11(2):459–470.
[94] Sauer F, Jäckle H. Dimerization and the control of transcription by Krüppel. Nature. 1993;364:454–
457.
[95] Gillespie D. A general method for numerically simulating the stochastic time evolution of coupled
chemical reactions. Journal of Computational Physics. 1976;22(4):403–434. doi:10.1016/0021-
9991(76)90041-3.

[96] Gillespie D. Exact stochastic simulation of coupled chemical reactions. Journal of Chemical Physics.
1977;81(25):2340–2361. doi:10.1021/j100540a008.
[97] Hattne J, Fange D, Elf J. Stochastic reaction-diffusion simulation with MesoRD. Bioinformatics.
2005;21(12):2923–2924. doi:10.1093/bioinformatics/bti431.
[98] Allen RJ, Warren PB, ten Wolde P. Sampling Rare Switching Events in Biochemical Networks.
Physical Review Letters. 2005;94(1):018104. doi:10.1103/PhysRevLett.94.018104.
[99] Warren PB, ten Wolde PR. Enhancement of the Stability of Genetic Switches by
Overlapping Upstream Regulatory Domains. Physical Review Letters. 2004;92(12):128101.
doi:10.1103/PhysRevLett.92.128101.

[100] Warren PB, ten Wolde PR. Chemical Models of Genetic Toggle Switches. Journal of Physical
Chemistry B. 2005;109(14):6812–6823. doi:10.1021/jp045523y.
[101] Morelli MJ, Tănase-Nicola S, Allen RJ, ten Wolde PR. Reaction Coordinates for the Flipping of
Genetic Switches. Biophyiscal Journal. 2008;94(9):3413–3423. doi:10.1529/biophysj.107.116699.
[102] Abu-Arish A, Porcher A, Czerwonka A, Dostatni N, Fradin C. High Mobility of Bicoid Captured
by Fluorescence Correlation Spectroscopy: Implication for the Rapid Establishment of Its Gradient.
Biophyiscal Journal. 2010;99(4):L33–L35. doi:10.1016/[Link].2010.05.031.

35
 1

SUPPLEMENTARY INFORMATION
Stable developmental patterns of gene expression without
morphogen gradients
Maciej Majka1 , Nils B. Becker2,3 , Pieter Rein ten Wolde2 , Marcin Zagorski1 , and Thomas R. Sokolowski2,4,∗
arXiv:2306.00537v3 [[Link]-ph] 27 Aug 2024

1 Institute of Theoretical Physics and Mark Kac Center for Complex Systems Research,
Jagiellonian University, ul. prof. Stanislawa Lojasiewicza 11, 30-348 Kraków, Poland
2 FOM Institute AMOLF, Science Park 104, 1098 XG Amsterdam, The Netherlands
3 Present address: Theoretical Systems Biology, German Cancer Research Center,
69120 Heidelberg, Germany
4 Present address: Frankfurt Institute for Advanced Studies (FIAS), Ruth-Moufang-Straße 1,
60348 Frankfurt am Main, Germany
∗ Corresponding author, e-mail: sokolowski@[Link]

S1.1 Chemical equations describing the gene regulatory network

The biochemical system that we model is governed by the following set of chemical equations:

Monomer production, dimerization, monomer and dimer degradation

D
β kon µM µD
PA −−→ A −−⇀
A↽ −− A2
D
A −−→ ∅ A2 −−→ ∅
koff

Promoter repression and derepression (1st repressor dimer binding)

kR [C2 ] kR [B2 ] kR [D2 ]
PA ↽ A
PA ↽ A
PA ↽ A
on on on
−−−
−−−−
−⇀− PC2 −−−
−−−−
−⇀− PB2 −−−
−−−−−⇀− PD2
off off off
ks kw kw

Promoter repression and derepression (2nd repressor dimer binding)

kR [C2 ] kR [B2 ]
PAB2 ↽ A
PAC2 ↽ A
on on
−−−
−−−−
−⇀− PB2 C2 −−−
−−−−
−⇀− PB2 C2
off off
ks kw

kR [D2 ] kR [B2 ]
PAB2 ↽ A
PAD2 ↽ A
on on
−−−
−−−−−⇀− PB2 D2 −−−
−−−−
−⇀− PB2 D2
off off
kw kw

kR [D2 ] kR [C2 ]
PAC2 ↽ A
PAD2 ↽ A
on on
−−−
−−−−−⇀− PC2 D2 −−−
−−−−
−⇀− PC2 D2
off off
kw ks

Promoter repression and derepression (3rd repressor dimer binding)

kR [C2 ] kR [D2 ] kR [B2 ]
PAB2 D2 − A
PAB2 C2 − A
PAC2 D2 − A
on on on
−−
↽ −−−−
−⇀− PB2 C2 D2 −−
↽ −−−−−⇀− PB2 C2 D2 −−
↽ −−−−
−⇀− PB2 C2 D2
off off off
ks kw kw
 2

Note that here, for brevity, we only show the complete set of reactions for the promoter of gene A, and its
products (A monomers and dimers A2 ). An analogous set of equations holds for the promoters of genes B, C
and D, with PA replaced by PB , PC or PD , respectively, all occurrences of the corresponding dimers (B2 , C2
or D2 ) replaced by A2 and corresponding monomers A replaced by B, C, or D, respectively, and the unbinding
rates adjusted to the respective strength of mutual interactions.
In addition to the shown set of equations, we provide a graph-like summary of the reaction network in
Fig. S12.

S1.2 Analysis of perturbation experiments for assessing pattern restor-

ing forces
Perturbation experiments on patterns initially relaxed into their long-persisting intact state were carried out by
expanding expression domains at the expense of the respective antagonistic gene’s domains and simulating the
subsequent dynamics of the system, as described in the main text and Methods, Sec. 4. The stochastic time
trajectories of such perturbation experiments that were repeated many times were then analysed as follows:
In order to quantify the spatial properties of a given domain G we used its center of mass, zG , where G
stands for one of the five domains {Aa , B, C, D, Ap }, with Aa and Ap marking the anterior and posterior parts
of A. Advantageously, zG is a robust measure, as it avoids ambiguity associated with determining domain
boundaries in the presence of gene expression noise. Figure S1 shows, for κ = κopt , time traces of zG for the
two types of perturbations, averaged over 10 independent samples in each case. For both perturbations the
average centers of copy number relax back to their original positions on a timescale ∼ 10 h. This demonstrates
that for optimal repression strength ratio an effective restoring force counteracts deviations from the five-stripe
pattern for varied λAC .

300 300

250 250

200 200
zG [µm]

zG [µm]

150 150

100 100

50 50

0 0
0 2 4 6 8 10 12 14 16 18 20 0 2 4 6 8 10 12 14 16 18 20
t [h] t [h]
Aa Ap B C D Aa Ap B C D

Figure S1. Perturbed trajectories are restored to their origin at optimal NN repression strength.
Shown are averaged time traces of the copy number center-of-masses zG for the five domains of the stable
pattern (Aa = anterior A domain, Ap = posterior A domain) at κ = κopt for two different perturbations:
(A) “C expansion”, i.e. prolongation of the central C domain by ∆ = 8 nuclei into the posterior and (B) “A
expansion”, i.e. prolongation of the anterior A domain by ∆ = 8 nuclei towards the center of the embryo.
The gray-dashed line marks the center of the system. In both cases we observe a restoration of the metastable
state on a timescale ≲ 10 h.
 3

S1.3 Estimation of phase space diffusion coefficient from overdamped

Langevin dynamics
Let f (X, Y, t) be a twice differentiable real function depending on a two-dimensional diffusion-drift processes
⃗ = (X, Y ) and time t (explicitly). In the overdamped Langevin limit, i.e. assuming that the displacements
X
of the random walker are governed only by forces that stem from an underlying force field and by Gaussian
noise, and that its accelerations and inertia are negligible, we can describe this random processes via
⃗ = ⃗v dt + σdW
dX ⃗ (S1)

where W ⃗ is a (two-dimensional) Wiener processes and ⃗v = (vX , vY ) a (local) drift velocity resulting from the
potential forces.
We then can calculate the differential of f with Itō’s Lemma (as a generalization of Taylor expansion) as
follows:
∂f ∂f
df (X, Y, t) = σ dWX + σ dWY
∂X
 ∂Y 
∂f ∂f ∂f σ2 ∂ 2 f σ2 ∂ 2 f 2
2 ∂ f
+ + vX + vY + + + ζσ dt
∂t ∂X ∂Y 2 ∂X 2 2 ∂Y 2 ∂X∂Y
(S2)

Here ζ measures the correlation between X and Y .

In order to apply this general formula to the specific diffusion-drift problem for the phase space coordinates
(λAC , λBD ) defined in the main text Eq. (2), we assign for brevity λx = λAC and λy = λBD , and then we set
X = λx , Y = λy and f (X, Y, t) = f (λx , λy ) = (λx − λx0 )2 + (λy − λy0 )2 ≡ ∆λ2 (the squared displacement
function).
Itō’s Lemma now reads (note that the time and mixed derivatives vanish):
   
d(∆λ2 ) = d (∆λx )2 + (∆λy )2 = d (λx − λx0 )2 + (λy − λy0 )2
≃ 2(λx − λx0 )(vλx dt + σdWx ) + 2(λy − λy0 )(vλy dt + σdWy )
σ2 σ2
+ 2dt + 2dt (S3)
2 2
To relate the above formula to the displacements sampled in our simulations with a fixed acquisition
time interval ∆t we shall integrate the infinitesimal contributions over this interval. At the same time we
take the ensemble average to account for the averaging of independent samples, which causes the Gaussian
terms σdWx and σdWy to vanish. We further assume that, to a good approximation, the drift velocities and
diffusion coefficients are constant over the time interval ∆t and diffusion isotropic in λx and λy direction, i.e.
 4

√
Dλx = Dλy = Dλ (⃗λ). Finally, using σ = 2Dλ , we obtain:
Z 
∆λ2 = d(∆λ2 )
*Z∆t + *Z +
∆t ∆t
= 2 [λx (t) − λx (0)] vλx dt + 2 [λy (t) − λy (0)] vλy dt
0 | {z } |{z} 0 | {z } |{z}
≃⟨vλx ⟩t ≃⟨vλx ⟩ ≃⟨vλy ⟩t ≃⟨vλy ⟩
*Z + Z Z
∆t
+ 4Dλ dt + ⟨2∆λx σdWx ⟩ + ⟨2∆λy σdWy ⟩
0 ∆t | {z } ∆t | {z }
0 0
* Z + * Z +
∆t ∆t
2 2
≃ ⟨vλx ⟩ 2tdt + vλy 2tdt + 4 ⟨Dλ ⟩ ∆t
0 0
2 2
≃ ⟨vλx ∆t⟩ + vλy ∆t + 4 ⟨Dλ ⟩ ∆t
2 2
= ⟨∆λx ⟩ + ⟨∆λy ⟩ + 4 ⟨Dλ ⟩ ∆t (S4)

The final result shows that, knowing the average displacements ⟨∆λx ⟩ and ⟨∆λx ⟩ and average squared
displacements ∆λ2 at ⃗λ, we can compute the average diffusion coefficient ⟨Dλ ⟩ (⃗λ) via:

1 h 
2 2
i 1
⟨Dλ ⟩ (⃗λ) = ∆λ2 (⃗λ) − ⟨∆λx ⟩ (⃗λ) + ⟨∆λy ⟩ (⃗λ) = V⟨λ⟩ (⃗λ) (S5)
4∆t 4∆t
The bracket term containing the first moments corrects the mean squared displacement for the contributions
coming from the deterministic drift and tends to zero as the process becomes purely diffusive.
The results of this analysis applied to our simulated systems are shown in Fig. S2, where we plot the diffusion
coefficients and the resulting expected diffusion times to the edge of the basin of initial patterns as a function
of our main parameter, the repression strength ratio κ. Interestingly, we find that the diffusion constant
remains almost constant for large κ, corresponding to weak nearest-neighbor repression (Fig. S2A). In this
regime, the estimated time of diffusion to the edge of the stable basin, τ0.2 is about 12 hours and comparable
to the recorded average stability times (∼ 20 hours, panels B and C). When κ is reduced, meaning that the
nearest-neighbor repression is increased, we observe an increase of the diffusion constant and correspondingly a
reduction of τ0.2 to values below 10 hours. Since in this regime the pattern stability increases dramatically, our
analysis corroborates the finding that the stability increase is not due to a slow-down of the pattern boundary
dynamics, but due to the emergence and deepening of a metastable basin that generates restoring forces to
pattern perturbations.
 5

A B

Figure S2. Average phase space diffusion coefficients as a function of κ. (A) The diffusion coefficient
of phase space trajectories in the (λx , λy ) = (λAC , λBD ) space are obtained from the overdamped Langevin
analysis, see section S1.3. The diffusion coefficients are averaged over the phase space region RP = [0.3, 0.4]2 ,
which is part of the diffusive plateau, for different repression strength ratios κ. The PSU stands for the phase
space units. (B) The resulting approximate diffusion times from the phase space region of the five-stripe
relaxed patterns, towards the edge of the diffusive plateau as a function of κ, assuming a distance of 0.2
PSU for the initial phase space distance to the edge. The edge of diffusive plateau is defined as a region of
(λAC , λBD ) from which the states are quickly absorbed into the regions with one the domains lost.
 6

S1.4 Supplementary velocity field figures

Figures S3, S4, S5 and S6 show phase space velocity fields and trajectories starting from perturbed initial
conditions for three different projections of the reaction coordinates, for the following cases:
• for the repression strength ratio κ ≃ 31.6 that optimizes pattern stability with pinning, see Fig. S3,
• for weaker NN repression strength κ = 1000 with pinning, see Fig. S4,

• for the optimal repression strength ratio κ = 100 that optimizes pattern stability without pinning, see
Fig. S5,
• for weaker NN repression strength κ = 1000 without pinning, see Fig. S6.

S1.5 Example pattern destruction processes

Figures S7, S8 and S9 summarize typical pattern destruction processes traced in biased simulated time (within
the NS-FFS scheme) as the reaction coordinate λ = λAC + λBD increases, for the system with pinning and at
optimal repression strength ratio κ ≃ 30 (maximally stable regime).
In all figures, panel (A) shows snapshots of the circumference-averaged expression patterns, initially consist-
ing of four expression domains, some of which are destroyed as λ increases. Panels (B) show the corresponding
increase of λ. Panels (C) illustrate how the four gene expression domain sizes (defined as the number of posi-
tions at which the expression level of the respective gene is above a predefined threshold of 3 copies) change as
biased simulation time and λ increase–initially starting from approximately equal domain sizes, some domains
are eliminated while their counterpart domains increase in size and ultimately can span the whole system.
More specifically, Fig. S7 shows an example destruction process in which the expression domain of gene
C (red) is destroyed first. Notice that once the C domain is eliminated, this appears to enable subsequent
expansion of the D domain (black) at the expense of the opposing B domain (blue).
Fig. S8 shows an example in which first the B domain (blue) is destroyed by the counterpart D domain
(black). This is accompanied by a steady shrinkage of the C domain (red) by the A domain (green); however,
the A domain does not manage to destroy the C domain completely by the end of the simulated time interval.
In contrast, Fig. S9 shows a case in which two genes manage to destroy their opposing counterparts’
expression domains by that time point: First the A domain grows at the expense of the C domain, while
growth of B at the expense of D follows afterwards with reduced speed. Notice that the early asymmetry
between B and D is transiently reduced before B starts to push out D ultimately completely, highlighting a
(partial) pattern restoration event.
 7

Figure S3. Average phase space velocities for the maximally stable system (κ ≃ 31.6) with pinning.
Left plots (A, C, E) show local average phase space velocities, right plots (B, D, F) additionally show example
trajectories for the two types of perturbations considered in the restoration experiments (blue = pert. from
boundary, turquoise = pert. from center). Starting points are marked by black bullets.
 8

Figure S4. Average phase space velocities for weaker NN interaction (κ = 1000) in the system with
pinning. Left plots (A, C, E) show local average phase space velocities, right plots (B, D, F) additionally
show example trajectories for the two types of perturbations considered in the restoration experiments (blue
= pert. from boundary, turquoise = pert. from center). Starting points are marked by black bullets.
 9

Figure S5. Average phase space velocities for the maximally stable system (κ = 100) without pinning.
Left plots (A, C, E) show local average phase space velocities, right plots (B, D, F) additionally show example
trajectories for the two types of perturbations considered in the restoration experiments (blue = pert. from
boundary, turquoise = pert. from center; purple = pert. from center resulting in complete pattern destruction).
Starting points are marked by black bullets.
 10

Figure S6. Average phase space velocities for weaker NN interaction (κ = 1000) in the system without
pinning. Left plots (A, C, E) show local average phase space velocities, right plots (B, D, F) additionally
show example trajectories for the two types of perturbations considered in the restoration experiments (blue
= pert. from boundary, turquoise = pert. from center). Starting points are marked by black bullets.
 11

Figure S7. Example destruction event: Destruction of the C domain at optimal NN repression
strength (κ = κopt ≃ 30) in the system with pinning. Panel (A) shows snapshots of gene expression
patterns (averaged along the circumference) at the indicated times of the biased simulation trajectory. Panel
(B) shows the progress of the reaction coordinate λ = λAC + λBD along the simulated time trajectory; red
bullets correspond to the snapshots in (A) in the same temporal order, red dotted lines show the corresponding
λ levels. Panel (C) shows the corresponding evolution of the expression domain sizes of the four individual
genes, here defined by the number of positions at which the (circumference-averaged) expression level is above
3 copies (dotted lines in (A) panels). The C domain (red) starts to disappear from 2.13 h onwards and is
completely destroyed by the A domain (green) by 3.73 h. Notice the posterior movement of the boundary
between the B and D domains that appear to precede A destruction in this case. Later the B-D boundary
moves back towards the anterior again.
 12

Figure S8. Example destruction event: Destruction of the B domain at optimal NN repression
strength (κ = κopt ≃ 30) in the system with pinning. Panel (A) shows snapshots of gene expression
patterns (averaged along the circumference) at the indicated times of the biased simulation trajectory. Panel
(B) shows the progress of the reaction coordinate λ = λAC + λBD along the simulated time trajectory; red
bullets correspond to the snapshots in (A) in the same temporal order, red dotted lines show the corresponding
λ levels. Panel (C) shows the corresponding evolution of the expression domain sizes of the four individual
genes, here defined by the number of positions at which the (circumference-averaged) expression level is above
3 copies (dotted lines in (A) panels). The B domain (blue) starts to be suppressed by the D domain (black)
already from 1.03 h onwards and is completely destroyed by 2.17 h. Notice that the destruction of the B
domain appears to be accompanied by a compression of the C (red) domain by the posterior A (green) domain
at 0.57 h.
 13

Figure S9. Example destruction event: Destruction of both the C and D domains at optimal
NN repression strength (κ = κopt ≃ 30) in the system with pinning. Panel (A) shows snapshots of
gene expression patterns (averaged along the circumference) at the indicated times of the biased simulation
trajectory. Panel (B) shows the progress of the reaction coordinate λ = λAC + λBD along the simulated time
trajectory; red bullets correspond to the snapshots in (A) in the same temporal order, red dotted lines show the
corresponding λ levels. Panel (C) shows the corresponding evolution of the expression domain sizes of the four
individual genes, here defined by the number of positions at which the (circumference-averaged) expression
level is above 3 copies (dotted lines in (A) panels). Notice that in this example the C (red) domain is under
pressure already at 0.75 h, but then rebounds again at 1.08 h, demonstrating the pattern’s restoring capability.
In parallel, the interface between B (blue) and D (black) moves posteriorly, which initiates the destruction of
the C domain at 1.83 h, completed by 3.10 h. Subsequently, the B domain also expels the D domain across
the posterior system boundary, enabled by the fact that D is not pinned there (in contrast to A).
 14

Figure S10. Pinning affects destruction pathways We plot here the average probability fluxes from the
region of stable patterns RS into the different remote basins identified in the (λAC , λBD ) phase space as a
function of the repression strength ratio κ for the systems with and without pinning. Here the flux is defined
as the average increase per time of the total probability in the basin. Basin boundaries and flux quantities are
†
described in detail in Methods, Sec. 4. Shown are the flux into the basin RAC , corresponding to destruction
†
of either the A or C domain (red lines), the flux into the basin RBD , in which either the B or D domain breaks
down (blue lines), and the total outflux from RS , which equals the pattern destruction rate kD (black bullets).
Solid lines and triangles show the data for the system with pinning, dashed lines and circles the values for the
†
system without pinning. Clearly, in both with and without pinning and for all κ considered here, RAC is the
dominant fraction of the flux, reflecting that the dominant pathway to destruction is the one that starts with
the disappearance of either the A or C domain. Pinning of A expression at the system boundaries leads to a
pronounced reduction of the flux through this pathway for κ = 10 − 100.
 15

S1.6 Pinning alters the proportion of destruction pathways

Both with and without pinning of A expression at the system boundaries, pattern stability is maximal at
an optimal strength of NN repression. Stability times, however, are significantly higher in the system with
pinning. In order to understand whether this is simply due to the fact that pinning prohibits destruction of
the A domains or due to other pinning-induced effects, we compared the different pathways to destruction
by computing probability fluxes through distinct reaction pathways (see Methods, Sec. 4 for details). The
different reaction pathways are defined by the order in which gap gene domains are destroyed. In our system
there are two major pathways: the A-C destruction pathway (either the A or the C domain vanishes first)
and the B-D destruction pathway (either the B or D domain vanishes first). The phase space histograms
in Figure 3 of the main text demonstrate that simultaneous destruction of two domains, corresponding to
trajectories that progress diagonally in (λAC , λBD ) space, is highly improbable. We find that, while in general
the A-C destruction pathway prevails, the fact that the A-destruction pathway is dominant for κ ≤ 100 in the
system without pinning accounts for the strong enhancement of pattern stability due to pinning.
In Figure S10 we plot for different repression strength ratios κ the magnitude of average fluxes from the
region of intact patterns RS in the (λAC , λBD ) space into the respective neighboring regions that correspond
to states in which one expression domain vanished. The figure reveals that for all κ the flux through the A-C
destruction pathway is approximately ten times higher than the flux through the B-D pathway, for systems both
with and without pinning. The figure also shows that pinning indeed reduces the flux through the dominant,
i.e. A-C, pathway, most significantly for κ ≃ 10 − 100, i.e. around the optimal value κopt . This gives rise to
the pronounced stability enhancement. The simultaneous reduction of the flux through the B-D pathway is
not relevant for overall stability.
We analysed further the detailed composition of fluxes through the dominant (A-C) pathway by computing
the average flux into the regions of destroyed states in ([A]tot /Ntot , [C]tot /Ntot ) space, see Fig. S11. As
expected, in the systems with pinning the entire flux through the dominant pathway goes into the C-destroyed
state. Interestingly, this is also the case for the weakly coupled systems without pinning. Here the flux into
the A-destroyed state is clearly dominant over the flux into the C-destroyed state for strong NN interaction.
This explains why pinning, which prohibits exit through the B-destruction pathway, increases stability in the
κ ≲ 100 regime. While the flux through the C-destruction pathway is minimal at κ = 31.6 with or without
pinning, in the system without pinning the accessibility of A destruction shifts the minimum of the combined
flux through both pathways towards κ = 100, see Fig. S10.
 16

Figure S11. Pinning shifts the destruction flux balance in the dominant (A-C) destruction pathway
The figure shows the contributions of the A-destruction and C-destruction pathways to the outflux from RS as
a function of the repression strength ratio κ for the systems with (A) and without (B) pinning. See Methods,
Sec. 4 for the definitions of basin boundaries and details of flux calculation. Without pinning and for strong
NN repression, the preferred pathway to destruction is the one in which the A domains are destroyed first,
while for weaker coupling (large κ) destruction begins via annihilation of the C domain. Interestingly, in both
cases the flux through the C-destruction pathway is minimal at κ ≃ 31.6. However, in the system without
pinning this value falls into the regime in which the flux through the A-pathway markedly increases. Pinning
forbids destruction via the A pathway and thus dramatically reduces the overall destruction flux for low κ in
the system with pinning, giving rise to the additional enhancement of optimal stability at κ ≃ 31.6.
 17

Quantity Symbol Value Unit

Geometry
Nuclear radius rN 2.5 µm
Nuclear volume VN 65.4 µm3
No. of nuclei in axial direction Nz 40
- resulting system length L 340 µm
No. of nuclei in circumferential direction Nϕ 8
Production / degradation
Protein production rate β 0.20 s−1
Monomer degradation rate µM 0.05 s−1
Dimer degradation rate µD 0.005 s−1
- resulting effective degr. rate µeff 0.0095 s−1
Binding / unbinding
Intranuclear diffusion const. DN 3.2 µm2 /s
Repressor target site radius σR 0.5 µm
R
- resulting (diff. ltd.) repressor on-rate kon 20.1 µm3 /s
Standard (strong) repressor off-rate off
ks 0.06 s−1
off R,s
Weak repressor off-rate kw varied ≥ koff
Monomer protein radius σM 0.05 µm
D
- resulting (diff. ltd.) dimerization forward rate kon 4.0 µm3 /s
D
Dimerization backward rate koff 0.062 µm3 /s
Internuclear diffusion
Standard internuclear diffusion const. D 1.0 µm2 /s
Internuclear lattice distance l 8.5 µm

Table S1. The standard parameters of the simulated model of four

mutually repressing genes prearranged in the “alternating cush-
ions” pattern on a cylindrical lattice of expressing nuclei.
 18

β δ
PA A A2
ϵ
µM µD

α[B2 ] α[C2 ] σ ω ∅
ω α[D2 ]

PBA2 PCA2 PDA2

α[D2 ] ω α[C2 ] σ α[B2 ] ω α[D2 ] ω α[C2 ] σ α[B2 ] ω

PBA2 D2 PBA2 C2 PCA2 D2 PBA2 D2

α[D2 ] ω α[B2 ] ω
σ α[C2 ]
α[C2 ] σ

PBA2 C2 D2

Figure S12. Reaction network. This schematic shows the set of reactions that affect production and
degradation of a single gap gene species A. The strong repressor of A is denoted by B, the weak interaction
partners by C and D. For each species, X denotes the monomer, X2 the dimer. For easy readability here
R
we abbreviate: α ≡ kon = diffusion limited repressor binding rate; σ ≡ ksoff = next-nearest neighbor /
off D
strong repressor unbinding rate; ω ≡ kw = nearest-neighbor / weak repressor unbinding rate; δ ≡ kon =
D
dimerization forward rate; ϵ ≡ koff = dimerization backward rate. The schematic shows the reactions for the
promoter of species A and its protein products, which we denote by (A||C, B, D), but holds similarly for all
other combinations of regulated and regulating species, (B||D, A, C), (C||A, B, D) and (D||B, A, C); the order
of the regulating species is not strictly alphabetical because the strongly repressing species is indicated first.