Bridging scales in a multiscale pattern-forming system

Laeschkir Würthner1,*, Fridtjof Brauns1,2,*, Grzegorz Pawlik3, Jacob
Halatek1,4, Jacob Kerssemakers3, Cees Dekker3,†, and Erwin Frey1,5,†

1Arnold Sommerfeld Center for Theoretical Physics and Center for NanoScience,
Department of Physics, Ludwig-Maximilians-Universität München, Theresienstraße
37, D-80333 München, Germany

2Present address: Kavli Institute for Theoretical Physics, University of
California Santa Barbara, Santa Barbara, CA 93106, USA

3Department of Bionanoscience, Kavli Institute of Nanoscience Delft, Delft
University of Technology, Van der Maasweg 9, 2629 HZ Delft, the Netherlands

4Research Department, Oxford BioMedica Ltd., Windrush Court, Transport Way,
Oxford, OX4 6LT, UK

5Max Planck School Matter to Life, Hofgartenstraße 8, D-80539 Munich, Germany

*L.W. and F.B. contributed equally to this work.
†Corresponding authors. Email<EMAIL_ADDRESS>or<EMAIL_ADDRESS>

###### Abstract

Self-organized pattern formation is vital for many biological processes.
Reaction-diffusion models have advanced our understanding of how biological
systems develop spatial structures, starting from homogeneity. However,
biological processes inherently involve multiple spatial and temporal scales
and transition from one pattern to another over time, rather than progressing
from homogeneity to a pattern. To deal with such multiscale systems, coarse-
graining methods are needed that allow the dynamics to be reduced to the
relevant degrees of freedom at large scales, but without losing information
about the patterns at the small scales. Here, we present a semi-
phenomenological approach which exploits mass-conservation in pattern
formation, and enables to reconstruct information about patterns from the
large-scale dynamics. The basic idea is to partition the domain into distinct
regions (coarse-grain) and determine instantaneous dispersion relations in
each region, which ultimately inform about local pattern-forming
instabilities. We illustrate our approach by studying the Min system, a
paradigmatic model for protein pattern formation. By performing simulations,
we first show that the Min system produces multiscale patterns in a spatially
heterogeneous geometry. This prediction is confirmed experimentally by in
vitro reconstitution of the Min system. Using a recently developed theoretical
framework for mass-conserving reaction-diffusion systems, we show that the
spatiotemporal evolution of the total protein densities on large-scales
reliably predicts the pattern-forming dynamics. Our approach provides an
alternative and versatile theoretical framework for complex systems where
analytical coarse-graining methods are not applicable, and can in principle be
applied to a wide range of systems with an underlying conservation law.

## Introduction

Pattern formation is fundamental for the spatiotemporal organization of
biological processes, such as cell division, chemotaxis, and morphogenesis.
More than half a century ago, Turing showed theoretically how local
interactions (chemical reactions) and diffusion of chemical species can lead
to spontaneous spatial patterns [1]. Such reaction–diffusion systems have been
successfully used to explain pattern formation phenomena in nature that arise
self-organized from a stable homogeneous steady state [2, 3, 4, 5]. The
analysis proposed by Turing allows to predict the emergence of patterns with a
characteristic length scale as long as the entire dynamics remains in the
vicinity of the homogeneous steady state [6]. The validity of Turing’s
approach has been also tested experimentally for coupled chemical oscillators,
and was found to reliably predict the experimental observations, provided that
the model parameters are spatially and temporally uniform [7]. Pattern-forming
systems, however, are generally heterogeneous and therefore far from
homogeneity, and involve multiple spatial and temporal scales. An intriguing
example of biological pattern formation is morphogenesis, in which the
spatiotemporal patterns of morphogens dictate the future shape of an organism
that is orders of magnitude larger than its constituents [4]. On a smaller
scale, protein concentration patterns in cells are essential for the
spatiotemporal control of cellular processes such as cell division and
motility [8, 5, 9]. Protein patterns can exhibit fascinating multiscale
characteristics [10] and form in hierarchies of patterns on several scales
that affect one another [11].

Such complex multiscale biological processes involve many degrees of freedom
at multiple scales, rendering it difficult to analyze them and gain insight
into the underlying principles. To make progress on this issue, one needs to
use systematic coarse-graining schemes that allow the dynamics to be reduced
to the essential degrees of freedom at the relevant time and length scales.
For instance, a well-known and powerful method is the renormalization group
theory [12]. Unfortunately, this method is restricted to the vicinity of
critical points. The Mori-Zwanzig formalism [13] is another important approach
which allows to decompose the dynamics of a system into ‘fast’ and ‘slow’
variables by means of projection operators. One arrives at a closed set of
equations for the slow variables, while the fast variables are treated as
noise. One property that these methods have in common is that the scales that
have been integrated out or eliminated are not _resolved_ , and cannot be
recovered from the coarse-grained level of description. This is most apparent
in the Mori-Zwanzig formalism, where the eliminated degrees of freedom appear
effectively as noise terms on the resolved scales. For pattern-forming
systems, one is however interested in the patterns on the unresolved
scales111We adapt the term unresolved scales from the computational fluid
dynamics literature to refer to the (small) scales that have been integrated
out in the coarse-grained description. as they usually have a specific
function in biological systems. This raises the question of whether it is
possible to reconstruct information about the unresolved scales from the
dynamics at the resolved scales? Indeed, amplitude equations describe the
long-wavelength amplitude modulations of an underlying short-wavelength base
pattern and therefore resolve both the small and the large scales.
Unfortunately, however, they are limited to the vicinity of the supercritical
onset of pattern formation [6] (including weakly subcritical cases) and only
feasible in simple geometries where the orthonormal basis functions of the
diffusion operator can be found in closed analytical form. Hence, to fill
these gaps, one relies on new concepts to deal with multiscale systems.

Here, we propose a semi-phenomenological approach to overcome these
mathematical limitations in the concrete context of mass–conserving
reaction–diffusion (MCRD) systems. Recently, a new theoretical framework for
MCRD systems has been introduced [14, 15] that allows one to characterize
their dynamics in the highly nonlinear regime. The basic idea is to consider
reaction–diffusion system as decomposed into a set of reactive compartments
which are spatially coupled by diffusion. For an isolated compartment, one can
determine the steady state (_local equilibrium_) and its stability properties
which both depend on the total densities within that compartment. Since
diffusion causes the lateral redistribution of these total densities, these
local equilibria will change over time. This concept of _moving local
equilibria_ enables one to study the physical mechanisms underlying pattern
formation and characterize the dynamics far away from the homogeneous steady
state. The fact that one is able to characterize the dynamics far from
homogeneity suggests that the _local equilibria theory_ may be a promising
approach to study heterogeneous systems. We therefore asked whether the ideas
from local equilibria theory would be applicable to investigate multiscale
patterns?

To pursue this question, we use the Min protein system of E. coli which has
emerged as a paradigmatic model system for the study of pattern formation in
cell biology [16, 17, 18, 19, 20]. Its dynamics is driven by two proteins,
MinD and MinE, which cycle between cytosolic and membrane–bound states and
interact nonlinearly on the membrane (Fig. 1A). In E. coli, these proteins
oscillate from cell pole to cell pole and thereby position the cell division
machinery to midcell [16, 17]. Studying the Min dynamics in various
reconstituted systems has led to the discovery of a rich set of patterns
including traveling waves and spirals [18], chaotic patterns [21, 22, 23, 10],
“homogeneous pulsing” [24, 25, 26], as well as quasi-stationary labyrinths,
spots, and mesh-like patterns [10, 27]. Theoretical analysis of mathematical
models has lead to the key insight — and experimentally confirmed prediction —
that the average total densities of MinD and MinE and the bulk height are key
control parameters for pattern formation in the reconstituted Min system [5,
28]. The rich set of patterns, experimental accessibility in vitro and
theoretical understanding make the Min-system an ideal candidate to
investigate the role of spatial heterogeneity on pattern formation.

Since varying the bulk height affects the local equilibrium state and is a key
control parameter for pattern formation [5, 28], we study the Min dynamics in
a wedge–shaped geometry with a membrane placed on the bottom surface (Fig.
1B). While there are many distinct ways to introduce large–scale spatial
heterogeneities into the system, e.g. by introducing space-dependent kinetic
rates, we chose to use a wedge geometry because it is relatively easy to
implement experimentally. In numerical simulations, we find that the system
exhibits a striking range of transient patterns, that coexist in different
spatial regions along the membrane (Movie S1 and Fig. 1C). As time progresses,
patterns in different regions change and transition to other patterns.

Figure 1: (A) Schematic illustration of the Min-protein reaction network. (B)
Wedge-geometry with a membrane surface at the bottom plane ($z=0$) and bulk
height $H(x)$ increasing linearly along the $x$ direction. (C) Snapshot of the
membrane-density of MinD, obtained by numerically simulating the Min dynamics
Eqs. 1–3 in the geometry shown in (B). One observes regions with chaotic
patterns, standing waves (SW, dashed green outline) and traveling waves (TW)
along the membrane and at different bulk heights; see Movie S1.

To characterize these complex dynamics that play out on multiple spatial and
temporal scales, we generalize the concept of dispersion relations (obtained
from a linear stability analysis) by applying it to sections of the domain,
which we term _regional dispersion relations_. Combining this approach with
the local equilibria theory [14, 15, 8], we show that one can reconstruct the
type and characteristics of patterns on small scales from the local protein
mass densities, which we identify as the essential degrees of freedom on large
spatial and temporal scales, i.e. the “hydrodynamic variables” of the system.
The key to this reconstruction are correlations between the regional pattern
characteristics and instantaneous, regional dispersion relations, calculated
from the instantaneous regional mass densities. Over time, these masses change
due to diffusive redistribution, resulting in qualitatively different regional
dispersion relations that indicate the local pattern type in the system. This
reconstruction of small-scale features (on unresolved scales), together with a
coarse-grained description for the mass-redistribution dynamics on large
scales allows us to understand and predict the long–term temporal evolution of
the system. A major advantage of our approach is that it is based on a linear
theory and therefore conceptually and technically simple to apply.

A key prediction from our numerical simulations and theoretical analysis is
that different pattern types form at different positions along the wedge
shaped geometry. To test this prediction experimentally, we performed
experiments with a reconstituted Min system in wedge-shaped microfluidic
cells. In agreement with the theoretical prediction, we find a range of
transient patterns coexisting in different spatial regions along the membrane.

## Results

### The Min protein system in wedge geometry

Mathematically, the Min-protein dynamics is described by bulk-surface coupled
reaction–diffusion equations, which describe the concentrations of cytosolic
proteins MinD-ATP, MinD-ADP, and MinE,
$\mathbf{c}=(c_{\mathrm{DD}},c_{\mathrm{DT}},c_{\mathrm{E}})$, in the bulk
volume $\mathcal{V}$, and the concentrations of membrane-bound MinD and MinDE
complexes, $\mathbf{m}=(m_{\mathrm{d}},m_{\mathrm{de}})$, on the surface
$\mathcal{S}$. For the wedge geometry, in spatial coordinates
$\mathbf{x}=(x,y,z)$, we place the membrane surface (with lateral dimensions
$L\times L$) in the $x{-}y$ plane at $z=0$ and let the bulk height vary as a
linear ramp from $H_{0}$ to $H_{1}$ along the $x$-direction (see Fig. 1B).

The dynamics of bulk components $\textbf{c}(\textbf{x},t)$ is governed by the
equation

$\partial_{t}\mathbf{c}(\mathbf{x},t)=D_{c}\nabla^{2}\mathbf{c}+\Lambda\mathbf{c},$
(1)

where $D_{c}$ denotes the bulk diffusion constant and the matrix
${\Lambda=\text{diag}(-\lambda,\lambda,0)}$ describes nucleotide exchange of
MinD in the bulk. The dynamics of membrane components $\mathbf{m}(x,y,t)$ is
constrained to the membrane surface and takes the form:

$\partial_{t}\mathbf{m}(x,y,t)=D_{m}\nabla^{2}_{\mathcal{S}}\mathbf{m}+\mathbf{r}(\mathbf{c}|_{z=0},\mathbf{m}),$
(2)

where $D_{m}$ is the membrane diffusion constant and
${\nabla^{2}_{\mathcal{S}}=\partial_{x}^{2}+\partial_{y}^{2}}$ is the surface
Laplacian. The membrane reactions $\mathbf{r}$, which comprise attachment,
detachment, and recruitment processes of Min proteins, are specified in the
Materials and Methods section.

The dynamics in the bulk and on the surface are coupled by reactive boundary
conditions,

${-}D_{c}\partial_{z}\mathbf{c}|_{z=0}=\boldsymbol{f}(\mathbf{c}|_{z=0},\mathbf{m}),$
(3)

that describe the bulk fluxes induced by attachment and detachment of proteins
at the membrane (see Materials and Methods). At the remaining boundaries, no-
flux boundary conditions are imposed such that the system is closed. Together,
the above dynamics conserve the average mass densities of MinD and MinE:

$\displaystyle\bar{n}_{\mathrm{D}}\,|\mathcal{V}|$ $\displaystyle=\langle
m_{\mathrm{d}}+m_{\mathrm{de}}\rangle_{\mathcal{S}}\,|\mathcal{S}|+\langle
c_{\mathrm{D}}\,\rangle_{\mathcal{V}}|\mathcal{V}|\,,$ (4a)
$\displaystyle\bar{n}_{\mathrm{E}}\,|\mathcal{V}|$ $\displaystyle=\langle
m_{\mathrm{de}}\rangle_{\mathcal{S}}\,|\mathcal{S}|+\langle
c_{\mathrm{E}}\rangle_{\mathcal{V}}\,|\mathcal{V}|\,,$ (4b)

where $c_{\mathrm{D}}=c_{\mathrm{DD}}+c_{\mathrm{DT}}$ is the total cytosolic
MinD concentration; $\langle\cdot\rangle_{\mathcal{S}}$ and
$\langle\cdot\rangle_{\mathcal{V}}$ denote the mean on the surface and in the
bulk respectively; $|\mathcal{S}|$ and $|\mathcal{V}|$ are the total surface
area and bulk volume (see Materials and Methods).

Using finite element (FEM) simulations we investigated the spatiotemporal
dynamics of the Min system in wedge geometry. Our simulations show a broad
range of different patterns — including traveling waves, standing waves and
chaotic patterns — coexisting in different spatial regions of the membrane
(see Movie S1 and Fig. 1C). Interestingly, the regions where these patterns
are found change over time as the patterns transition from one type to
another. For long simulation times, we observe that patterns transition to
standing waves, such that the entire domain is covered by a single pattern
type in the final steady state. The pattern in steady state depends on the
specific choice of parameters, and therefore can be altered by changing the
model parameters (Fig. S1 and Movie S2).

### Experimental implementation

We tested our theoretical prediction on this multi-scale dynamics in an
experimental system consisting of a wedge-shaped microfluidic flow chamber
(Fig. 2A). The bottom and top surface of the wedge were covered with a
supported lipid bilayer consisting of DOPG:DOPC (30:70 %) which mimics the
natural membrane composition of E. coli [29]. The length of the wedge was
typically about $814\text{\,}\mathrm{m}\mathrm{m}$ and the width about
$34\text{\,}\mathrm{m}\mathrm{m}$. The bulk height range was approximately
$250\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$ (Fig. 2B). Min proteins
were distributed in the chamber by rapid injection of a solution containing
$1\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{M}$ MinD and
$1\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{M}$ MinE (including 10 %
fluorescently labelled MinD and MinE proteins for visualization), together
with $5\text{\,}\mathrm{m}\mathrm{M}$ ATP and an ATP-regeneration system [28].

Figure 2: Experimentally observed Min patterns in a wedge flow cell. (A)
Schematic presentation of the experimental setup. Both, the bottom and the top
surface (glass slides) are covered with a lipid bilayer. (B) Measurement of
the bulk height profile of the flow cell versus distance along the lateral
length of the wedge. The height was measured microscopically by z-stacks at
multiple spots. (C) Snapshot of the Min pattern along the wedge, the picture
was obtained by stitching individual adjacent images. Shown is a merge of MinD
(green) and MinE (red) channels. The bottom figure shows a kymograph for
intensities taken along the center line in the top figure.

Figure 2C shows a snapshot of Min protein patterns along the bottom surface of
the wedge geometry 30 minutes after injection. The experiments exhibit the
same essential hallmarks of multiscale Min protein patterns that we observed
in our numerical simulations. In particular, consistent with our simulations,
we observe a sequence of distinct spatiotemporal patterns coexisting in
different spatial regions of the membrane (Fig. 2C and Movie S3): At regions
of low bulk height (approximately between
$210\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$), one typically observes
chaotic patterns and standing waves, whereas traveling wave patterns emerge at
regions of large bulk height
($>$10\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$$). Furthermore, as in
the simulation, we observe a sharp boundary between regions that contain
traveling wave patterns and regions that contain rather chaotic and standing
wave patterns, and this boundary establishes quickly within a few minutes
(Fig. S2 and Movie S4). Overall, the observations provide a striking
verification of the height-dependent patterns predicted in the simulations.

There are also some differences between the patterns in the experiment and in
our numerical simulations. First, while we observed occasional transitions
from one pattern into another in our experiments (Fig. S3 and Movie S5), these
transitions occurred frequently and were more pronounced in the simulations.
This is explained by the lateral length of the experimental setup, that is
about an order of magnitude larger as compared to the simulation setup, which
is the main reason why we observe more frequent transitions between different
patterns in the simulations, as will become clear later. Second, in contrast
to the simulations, we noticed some homogeneous oscillations in the
experiments, which are characterized by large (homogeneous) density patches on
the membrane (typically few hundred micrometers in size) that oscillate with
time (Figs. S3– S4 and Movies S5– S7). We attribute this difference to the
following: Due to the fabrication method of the microfluidic flow chamber,
both the bottom and top surface of the wedge were covered with a supported
lipid bilayer. In recent work, it has been shown that membrane-to-membrane
crosstalk (i.e., between top and bottom surface) is responsible for the
emergence of homogeneous oscillations [28]. In our simulations, however, we
assume that Min proteins can only bind to the bottom membrane, which explains
why we do not observe homogeneous oscillations.

Taken together, we have a system that exhibits a fascinatingly rich transient
dynamics and involves patterns and transitions between them on multiple
spatial and temporal scales. We are therefore left with the key question: Can
we explain the cause why different patterns form in different spatial regions
and how they transition from one to another over time? Moreover, is it
possible to identify and reduce the system to its essential degrees of
freedom? A standard way to address these questions mathematically would be to
perform a multiscale analysis and to derive amplitude equations that describe
the large-scale spatiotemporal evolution of the pattern amplitudes [6]. This
would greatly simplify the problem as it allows to obtain a quantitative
relationship between the small-scale patterns and the large-scale dynamics
(slowly varying pattern amplitudes), thus ultimately enabling one to
reconstruct the patterns from the reduced dynamics at large length and time
scales [30, 31, 32, 33]. Carrying out this analysis requires determining the
set of orthogonal eigenmodes for the diffusion operator that satisfy the
boundary conditions. In a one-dimensional domain, these eigenmodes are simply
Fourier modes. Unfortunately, in the wedge geometry with bulk-surface
coupling, the eigenmodes can not be found analytically, thus precluding the
use of the amplitude equation framework. Moreover, amplitude equations are
restricted to the vicinity of supercritical and weakly subcritical
bifurcations [6, 34]. The Min patterns we observe here, however, are
generically subcritical [15] and exhibit large amplitudes [14, 28]. We
therefore aim to develop a new approach that overcomes these restrictions.

### Instantaneous, regional dispersion relations predict patterns

The analysis of pattern-forming systems usually starts with calculating the
homogeneous steady state (HSS) solutions and performing a linear stability
analysis around these states. This yields a _dispersion relation_ that informs
about the growth rate $\sigma(q)$ of small spatial perturbations with a
certain wavenumber $q$. However, the dispersion relation is generally only
informative in the vicinity of the homogeneous steady state [1, 6], and thus
unreliable for large amplitude patterns. Moreover, the spatial variation of
parameters even precludes the existence of a global HSS, so that a global
dispersion relation can no longer be determined. To overcome these
limitations, we adopt a semi-phenomenological approach where we generalize the
concept of dispersion relations.

Let us consider the wedge as dissected into a collection of two-dimensional
slices along the direction of constant bulk height. Each slice corresponds to
a rectangular geometry with a bulk height that depends on the position of the
slice in the wedge (see Fig. 3A,B). Next, for each slice and at each point in
time, we calculate instantaneous total densities of MinD and MinE, averaged
over the slice length $\langle\tilde{n}_{\mathrm{D,E}}\rangle_{y}(t,x)$
(Materials and Methods). The average total densities, together with the local
bulk height $H(x)$, then serve as parameters for the _regional dispersion
relation_ in each slice

$\sigma\left(q;H(x),\langle\tilde{n}_{\mathrm{D,E}}\rangle_{y}(t,x)\right),$
(5)

which is straightforward to determine because the slice represents a
rectangular geometry [14, 23, 28] (see Fig. 3A,B and Supplementary
Information). While the bulk height $H(x)$ varies linearly in space, the
average total densities $\langle\tilde{n}_{\mathrm{D,E}}\rangle_{y}(t,x)$ are
dynamic quantities and depend on the slice position $x$ as well as on time
$t$, since the diffusive coupling between the slices redistributes mass. It
follows that the regional dispersion relation depends on the spatial position
and is dynamic: $\sigma(q;x,t)$. This generalizes classical dispersion
relations, which are by definition independent of space and time.

Figure 3: (A) Rectangular geometry with membrane at the bottom edge
representing a slice through the three-dimensional in vitro system. (B) A
slice through the wedge geometry. For each such slice, at a given instance in
time, we calculate the instantaneous total densities, averaged along its
length $\langle\tilde{n}_{\mathrm{D,E}}\rangle_{y}(t,x)$, from the numerical
simulation data. From these slice-averaged total densities, we can then
calculate the corresponding local homogeneous steady state and its dispersion
relation. (C) _Dispersion relation_ with fastest growing mode $q^{*}$ and
right edge of the band of unstable modes $q_{\text{max}}$ indicated. The ratio
$q_{\text{max}}/q^{*}$ has been empirically found to correlate with the type
of the fully developed pattern, with a sharp transition from chaotic patterns
for $q_{\text{max}}/q^{*}<2$ to ordered patterns for $q_{\text{max}}/q^{*}>2$.
Close to the transition, standing waves are found, while travelling waves form
for larger ratios $q_{\text{max}}/q^{*}$ [14]. (D) Mode ratio
$q_{\text{max}}/q^{*}$ as a function of the slice position $x$ for a given
instance in time. The background shading indicates the type of pattern
expected from the “commensurability criterion.” (E) Representative snapshots
of the three distinct pattern types: spatiotemporal chaos, standing waves (SW)
and traveling waves (TW).

How does this spatially and temporally varying dispersion relation inform
about the system’s dynamics? As in uniform systems that exhibit homogeneous
steady states, it serves as a criterion for the onset of pattern formation and
for estimating the characteristic wavelength of the initial pattern that is
formed. While these insights are generally limited to the linear regime [1,
6], recent theoretical findings for the Min system in a two-dimensional
rectangular geometry (representing a slice geometry) have shown that the
dispersion relation reliably predicts the pattern type in the fully nonlinear
regime [5]. In particular, it was shown that depending on the total densities
of Min proteins, $\bar{n}_{\text{D}}$ and $\bar{n}_{\text{E}}$, and the bulk
height $H$, the system exhibits a variety of different patterns on the
membrane, such as chaos, standing waves, and traveling waves [14, 28].
Moreover, a careful analysis of numerical simulations has interestingly
revealed a strong one-to-one correlation between the dispersion relation and
the fully developed patterns in the highly nonlinear regime [14]: A
commensurability criterion between the unstable mode with the shortest
wavelength $q_{\text{max}}$ and the fastest growing mode $q^{*}$ has been
found that determines the pattern type (Fig. 3C–E). In short, it has been
shown that $q_{\text{max}}/q^{*}<2$ coincides with the regime of chemical
turbulence (spatiotemporal chaos), whereas for $q_{\text{max}}/q^{*}>2$ the
system exhibits ordered patterns (standing/traveling waves). Standing wave
patterns are found close to the commensurability transition
$q_{\text{max}}/q^{*}\gtrsim 2$, while traveling waves are found further away
from the threshold. In the following, we use this observed one-to-one
correspondence between the dispersion relation and the fully developed
patterns to reconstruct the small scale pattern types from coarse grained
densities.

To that end, we extracted the average total densities in each slice as a
function of time from the numerical simulation. Based on these densities we
then calculated the instantaneous regional dispersion relation in each slice
and extracted the ratio $q_{\text{max}}/q^{*}$ as a function of slice position
$x$ and time $t$ (Fig. 3C–E). The resulting pattern-type prediction is shown
in the space-time plot (kymograph) in Fig. 4A. Figure 4B shows the ratio
$q_{\text{max}}/q^{*}$ as a function of slice position $x$ for a set of
representative times (cf. Fig. 3D). The pattern-type prediction Fig. 4A is
then obtained from these ratios via the mapping shown in Fig. 3D,E.

Figure 4: (A) Kymograph showing the pattern-type prediction from the
commensurability criterion (cf. Fig. 3D). The green line shows
$x_{\mathrm{crit}}(t)$ where $q_{\text{max}}/q^{*}=2$, indicating the
transition from chaotic to ordered patterns. Green arrows mark the position
$x_{\mathrm{crit}}(t)$ for the times indicated by dashed white lines. (B)
Plots of the mode ratio $q_{\text{max}}/q^{*}$, determined from the local
dispersion relation, as a function of spatial position $x$ for several
representative times (dashed white lines in (A)). In the second to last row,
the entire domain is near the critical ratio $q_{\text{max}}/q^{*}=2$,
predicting the global emergence of standing waves (see last row). (C)
Snapshots of the membrane patterns (MinD density, cf. Fig. 1) from the full
numerical simulation. The green dashed line indicates $x_{\mathrm{crit}}(t)$.
Note the standing wave patterns found near $x_{\mathrm{crit}}(t)$. Their
fronts are aligned along the bulk height gradient such that the sequence of
wavenodes lies on lines of constant bulk height. (D) Machine-learning based
pattern classification using ilastik [35] (see Materials and Methods).

We find that this prediction correlates well with the patterns observed in the
full numerical simulation (Fig. 4C,D and Movie S8). In particular, the
temporally changing position $x_{\mathrm{crit}}(t)$, marking regions where
$q_{\text{max}}/q^{*}=2$ (indicated by the green arrows and dashed lines in
Fig. 4B and C), agrees with the position along the wedge where traveling wave
patterns transition to chaotic patterns. In the vicinity of
$x_{\mathrm{crit}}(t)$ we observe a band of standing waves as expected from
the “commensurability criterion” [14]. While the transition from chaos to
order at $q_{\text{max}}/q^{*}=2$ is sharp, we were not able to identify a
sharp criterion for the transition from standing waves to traveling waves.
Since the ratio $q_{\text{max}}/q^{*}$ and with it $x_{\mathrm{crit}}(t)$ are
entirely determined by the slice-averaged masses
$\langle\tilde{n}_{\mathrm{D,E}}\rangle_{y}(x,t)$, we conclude that these
masses are the essential degrees of freedom of the system at large scales.

Notably, we find that there are slight differences between the predictions and
the actual patterns for large times (see Fig. 4A–C). The reason for these
deviations lies in the model parameters, which were chosen such that the
entire domain is near the critical mode ratio $q_{\text{max}}/q^{*}=2$ for
large times. This renders the dynamics, and the prediction from the regional
dispersion relation highly sensitive to slight variations of the regional
total masses. Hence, the fact that our method is still able to qualitatively
predict the dynamics in this case underscores the robustness of our approach.
In the Supplemental Information, we provide additional results where the
parameters were chosen such that the mode ratio is deep in the traveling wave
regime ($q_{\text{max}}/q^{*}>2)$ for late times. In this case, we obtain an
excellent agreement between our predictions and the patterns observed in the
numerical simulations (see Fig. S1).

Next, we ask whether one can find an approximate coarse-grained dynamics for
these redistributed masses. Such a description would enable us to predict the
time evolution of the redistributed masses independently from the full
numerical simulations. One can then use the commensurability criterion to
predict the pattern types that will form in different spatial regions as a
function of the redistributed masses. In the next section we will show how one
can find such a description.

### Large-scale dynamics is driven by redistribution of mass

In general, mass redistribution between different spatial regions of the wedge
is caused by diffusive fluxes due to concentration gradients. Similar as in
the previous section, we consider here the redistribution of mass between
slices along the wedge (Fig. 3B). Since membrane diffusion is by two orders of
magnitude slower than bulk diffusion it may be neglected, such that
redistribution of protein mass between slices is governed by bulk diffusion
alone (Materials and Methods)

$\partial_{t}\langle n_{i}\rangle_{y,z}(x,t)\approx
D_{c}\langle\partial_{x}^{2}c_{i}\rangle_{y,z}+D_{c}\tfrac{\partial_{x}H(x)}{H(x)}\langle\partial_{x}c_{i}\rangle_{y,z},$
(6)

$\text{for }i=\mathrm{D,E}$. Here, the second term accounts for the spatial
variation of the bulk height, and thus the different volumes of neighboring
slices between which the diffusive flux
$D_{c}\langle\partial_{x}c_{i}\rangle_{y,z}$ redistributes mass. This can be
seen by rewriting Eq. (6) in the form of a continuity equation

$\partial_{t}\big{[}H(x)\cdot\langle
n_{i}\rangle_{y,z}(x,t)\big{]}\approx-\partial_{x}\big{[}H(x)\cdot
J^{\mathrm{diff}}_{i}\big{]}$ (7)

with the diffusive fluxes given by
$J^{\mathrm{diff}}_{i}:=-D_{c}\langle\partial_{x}c_{i}\rangle_{y,z}$. Since
the area of slices increases along the positive $x-\text{direction}$, the
diffusive fluxes $J^{\mathrm{diff}}_{i}$ on the right-hand side of Eq. (7) are
rescaled by the bulk height $H(x)$. These equations seem to be simple, but
unfortunately they are not closed, since the slice-averaged cytosolic
densities $\langle c_{i}\rangle_{y,z}(x,t)$ appear on the right hand side.

Figure 5: (A,B) Kymographs showing the total-density ratio of MinE to MinD
(E:D ratio) from the full numerical simulation (A) and from local-equilibria
based reduced dynamics (B). (C) Kymograph showing the pattern-type prediction
using the commensurability criterion based on the total densities from the
reduced dynamics. Note the excellent qualitative agreement to the pattern-type
prediction based on total densities from the full numerical simulation in Fig.
4A.

We are interested in the dynamics of $\langle n_{i}\rangle_{y,z}$ on
timescales much longer than typical oscillation periods of the patterns.
Therefore, following the intuition gained from previous works on MCRD systems
[15, 36], we assume that one can approximate the slice-averaged cytosol
concentrations by the homogeneous steady-state concentration in each slice

$\langle c_{i}\rangle_{y,z}(x,t)\approx
c_{i}^{*}(x,t):=c_{i}^{*}\bigl{(}H(x),\langle
n_{\mathrm{D}}\rangle_{y,z},\langle n_{\mathrm{E}}\rangle_{y,z}\bigr{)}.$ (8)

This assumes that the spatial average over many wavelengths in $y$-direction
is well approximated by the instantaneous homogeneous steady state in a slice.
These steady state concentrations only depend on the slices bulk height $H(x)$
and the slice-averaged total densities $\langle n_{i}\rangle_{y,z}(x,t)$.
Thus, the above approximation yields a closed set of equations for the mass-
densities

$\partial_{t}\langle n_{i}\rangle_{y,z}(x,t)\approx
D_{c}\partial_{x}^{2}c^{*}_{i}(x,t)+D_{c}\tfrac{\partial_{x}H(x)}{H(x)}\partial_{x}c^{*}_{i}(x,t).$
(9)

We will call this the _reduced dynamics_ in the following. Since the
homogeneous steady states may also undergo a saddle-node bifurcation,
characterized by the emergence of three steady states (two stable, one
unstable), this may lead to discontinuities in $c^{*}_{i}$. To regularize the
dynamics, $c_{i}$ is not set identical to $c_{i}^{*}$ but relaxes towards it
on a fast timescale (see SI for details).

Given the initial densities $\langle n_{i}\rangle_{y,z}(x,0)$, one can
numerically solve the reduced dynamics Eq. (9) to predict the entire time
evolution of the slice-averaged masses and hence the dispersion relation at
each point along the $x-\text{direction}$. Figure 5C shows the regional
pattern types predicted from the reduced dynamics. We find good qualitative
agreement for the distribution and transition of patterns as observed in the
numerical simulations (cf. Fig. 4A). The main difference to the full numerical
simulations is a slight quantitative deviation in the timescale, where the
dynamics predicted by Eq. (9) is slightly slower compared to the full
numerical simulation. We also note that the reduced dynamics predicts a larger
region of no instabilities as compared to the numerical simulations (cf. Figs.
4A and 5C). This is because the chaotic regime is rather narrow and close to
the regime for which the dispersion relation predicts no instability (cf.
Figs. 3D and 4B). In addition, since the patterns emerge from a subcritical
bifurcation [14] (a generic property of mass-conserving systems [15]), large
amplitude patterns can be excited and maintained even below the instability
threshold.

Figure 5A,B compare the time evolution of the slice-averaged total densities
from the full numerical simulation and the solution obtained from the reduced
dynamics. The colors in the kymographs indicate the total density ratio of
MinE and MinD (short, E:D ratio), which is a key control parameter in the Min-
protein dynamics [14].

## Discussion

Multiscale patterns in biological systems often emerge from hierarchical
systems, which are organized in a modular fashion. Each level of the hierarchy
instructs dynamics on the next level which operates on a smaller spatial
scale. For instance, along developmental trajectories of many organisms,
upstream patterns such as maternal gradients instruct downstream gene-
expression patterns on increasingly smaller scales [37, 11]. Importantly, on
each level of the hierarchy, there is a clean separation between (spatially
varying) control parameters and dynamical variables.

In contrast, in the system we have studied here, there is no such separation
as the globally conserved total densities play a dual role: they are both
dynamical variables and act as control parameters [14, 15]. Building on this
key feature has allowed us to explain and predict the intriguingly complex
patterns found in large-scale numerical simulations. The values of the total
densities of MinD and MinE locally control the pattern type: we showed that a
“regional dispersion relation” calculated from the regional average densities
reliably predicts the pattern type. At the same time, concentration gradients
in the bulk drive mass redistribution of MinD and MinE. Therefore, the total
densities are hydrodynamic variables on large scales which control pattern
formation on small scales. This separation of scales enabled us to derive a
reduced dynamics for the total densities on large spatial and temporal scales
which predicts the long-term dynamics of the system.

Notably, the dual role of total densities as dynamic variables and control
parameters also plays out at the small scale of the patterns themselves [14,
15]. Here, instantaneous _local_ total densities control _local_ equilibria
and their stability, which serve as proxies for the local dynamics. The local
dynamics cause gradients, which drive diffusive redistribution of the total
densities—in turn causing changes in the local dynamics. In the Min system,
this point of view has led to a detailed understanding of the emergence of
chaos near onset and of the transition to standing and traveling waves [14].
From a general perspective, the concept of local equilibria controlled by
total local densities is at the core of a number of recent theoretical
advances in the field of mass-conserving, pattern-forming systems [15, 36, 38,
8].

In addition to the dynamically changing total densities, the bulk height is
also a (fixed) heterogeneous control parameter in our system. The bulk height
(or more generally volume-to-surface ratio) is an important control parameter
for bulk-surface coupled pattern-forming systems [14, 28]. Here, the bulk
height gradient of the wedge serves to induce spatiotemporal heterogeneities
in the total densities. Alternatively, one could induce heterogeneities in the
total densities via spatial gradients of the kinetic rates or by imposing a
heterogeneous initial condition in the total densities. However, these
alternatives are difficult to realize experimentally in a reproducible and
controlled manner, which is the main reason why we chose the wedge setup in
this work. In a third scenario, large-scale gradients in the densities may
also emerge spontaneously and be maintained in the absence of “external”
heterogeneities.

An example for this third scenario is the Aranson–Tsimring model for pattern
formation in vibrated granular media [39] (see Materials and Methods for
details). In the following, we briefly discuss this model to put our approach
into a broader context. In particular, this model has been extensively studied
using amplitude equations allowing us to connect this mathematical approach to
the regional dispersion relations introduced here. The Aranson–Tsimring model
considers a system with a complex order parameter $\psi$ (describing the
surface modulation of a vibrated granular layer) which is coupled to a
conservation law for the grain density $\rho$ (see Eq. (24) in Materials and
Methods). Near the onset of pattern formation, this coupling gives rise to
localized patterns that have been studied using amplitude equations [32, 33,
30]. Figure 6 and Movie S9 illustrate how these patterns can be understood in
terms of regional dispersion relations. For high densities, there are no
unstable modes and no patterns form. Below a critical density $\rho_{c}$, a
band of unstable modes appears, giving rise to patterns through a
supercritical bifurcation. Indeed, localized patterns appear only where the
average regional density is below $\rho_{c}$ (see Fig. 6B). This demonstrates
the idea of regional dispersion relations in a nutshell. Moreover, it shows
that this approach gives rise to qualitatively similar insights as the
technically much more involved amplitude equation formalism. The conceptual
and technical simplicity of regional dispersion relations make this approach
readily applicable. The caveat is that this approach lacks the mathematical
rigor of the amplitude equation formalism and requires numerical solutions of
the dynamics as a basis.

Figure 6: Regional dispersion relations predict localized patterns in the
Aranson–Tsimring model Eq. (24). (A) Snapshot of the order parameter magnitude
$\psi$ showing localized patterns. Dashed white line indicates the stability
threshold determined from regional dispersion relations. (B) Coarse-grained
density (Gaussian filter with standard deviation 10). (C) Representative
dispersion relations in the stable and unstable regimes. Domain size:
$100\times 50$; see Materials and Methods for model details and remaining
parmeters.

Because the bifurcation in the Aranson–Tsimring model is supercritical [39],
we can immediately read off from the regional dispersion relation where small-
scale patterns will form. There is only one pattern type (stripes) and
therefore, no additional information is needed to reconstruct the small-scale
patterns. For the Min system considered here, on the other hand, the onset is
subcritical (i.e. patterns have large amplitude at onset), and the emergence
of a band of unstable modes alone does not inform about the pattern type. To
overcome this problem, we used the “commensurability criterion” that enabled
us to predict the small-scale patterns from regional dispersion relations.
However, this criterion has so far only been shown to hold for the Min-protein
system. Whether it also applies to other reaction–diffusion systems remains an
open question. In general, reconstructing subcritical small-scale patterns
from large-scale quantities will require that adequate criteria are first
identified in simplified settings (such as the “slice geometry” used here).
Supercriticality also guarantees that there is no multistability of different
pattern types near onset. Multistability would lead to hysteresis in the
transitions and therefore introduce memory in the system. As a result, the
one-to-one correspondence between the regional dispersion relation and the
pattern type that we have used here to reconstruct patterns would be lost.
Handling memory effects in pattern-forming systems remains an open issue that
will require the development of new methods, providing an interesting task for
future research.

Since conservation laws are ubiquitous in many physical systems, we believe
that our approach can be generalized to a broad class of multiscale pattern-
forming systems. For instance, mass conservation is inherent to particle–based
active matter systems. The local particle density controls emergent
orientational order, i.e. local symmetry breaking [40, 41, 42]. In turn,
orientational order controls mass redistribution due to the particles’ self-
propulsion. Thus, the particle density again plays a dual role as a control
parameter and a dynamic variable [42, 43, 44]. The dynamic interplay of mass
redistribution and orientational order has been shown to give rise to
coexistence of different macroscopic order (polar flocks, nematic lanes) and
the interconversion between them [42], not unlike the coexistence and
interconversion of different patterns we found for the reaction–diffusion
system studied in this work. One way to induce spatial heterogeneities in
these systems is to introduce a gradient of signaling chemicals
(chemoattractants) that affect the local velocity of active particles. This
would dynamically lead to redistribution of the particle densities on large
scales. Since the particle densities, in turn, are themselves control
parameters locally, non-trivial multiscale dynamics may emerge in such a
setup. Exploring the effects of such gradients in active matter systems could
be therefore an exciting task for future research.

On a broader perspective, our work shows how a linear analysis on small
scales, combined with a reduced description for non-linear large scale
dynamics (mass redistribution) can be employed to study complex multiscale
phenomena. We believe that our approach can be generalized and applied to
other multiscale systems with an underlying conservation law, such as
transport processes in porous media, combustion, and cell migration, to name a
few examples.

#### Acknowledgments.

E.F. acknowledges support from the Deutsche Forschungsgemeinschaft (DFG,
German Research Foundation) through the Collaborative Research Center (SFB)
1032 – Project-ID 201269156 – and the Excellence Cluster ORIGINS under
Germany’s Excellence Strategy – EXC-2094 – 390783311. C.D. acknowledges
support from the NWO/OCW Gravitation program NanoFront and the ERC Advanced
Grant 883684.

## Materials and Methods

### Mathematical model

We adopt the Min “skeleton model” introduced in Refs. [45, 46, 5] which is
known to qualitatively reproduce Min patterns _in vivo_ and _in vitro_ [46, 5,
28]. The governing equations are given in the main text, Eqs. [1]–[3]. The
membrane reactions are

$\mathbf{r}=\left[r_{\mathrm{D}}^{\mathrm{on}}-r_{\mathrm{E}}^{\mathrm{on}},r_{\mathrm{E}}^{\mathrm{on}}-r_{\mathrm{DE}}^{\mathrm{off}}\right]^{\top},$
(10)

with

$\displaystyle r_{\mathrm{D}}^{\mathrm{on}}$
$\displaystyle=(k_{\mathrm{D}}+k_{\mathrm{dD}}m_{\mathrm{d}})c_{\mathrm{DT}}\,,$
(11a) $\displaystyle r_{\mathrm{E}}^{\mathrm{on}}$
$\displaystyle=k_{\mathrm{dE}}m_{\mathrm{d}}c_{\mathrm{E}}\,,$ (11b)
$\displaystyle r_{\mathrm{DE}}^{\mathrm{off}}$
$\displaystyle=k_{\mathrm{de}}m_{\mathrm{de}}\,.$ (11c)

The reaction terms account for MinD attachment and self-recruitment to the
membrane, MinE recruitment by MinD, and dissociation of MinDE complexes with
subsequent detachment of both proteins into the cytosol, respectively.
Coupling between cytosol and membrane is established by reactive boundary
conditions at the membrane [cf. Eq. (3)]. The boundary fluxes are given by

$\boldsymbol{f}=\left[r_{\mathrm{DE}}^{\mathrm{off}},-r_{\mathrm{D}}^{\mathrm{on}},r_{\mathrm{DE}}^{\mathrm{off}}-r_{\mathrm{E}}^{\mathrm{on}}\right]^{\top}\,,$
(12)

which follows from mass conservation. For analytical caclulations we adapt the
following change of variables as it is more convenient: We describe the bulk
dynamics of MinD in terms of the variables
$c_{\mathrm{D}}=c_{\mathrm{DD}}+c_{\mathrm{DT}}$ and $c_{\mathrm{DD}}$, i.e.
in this case one defines the bulk concentration vector
$\mathbf{c}=(c_{\mathrm{D}},c_{\mathrm{DD}},c_{\mathrm{E}})$. The membrane
reaction in Eq. (11a) is then slightly modified by substituting
$c_{\mathrm{DT}}=c_{\mathrm{D}}-c_{\mathrm{DD}}$, and the boundary fluxes are
given by

$\boldsymbol{f}=\left[-r_{\mathrm{D}}^{\mathrm{on}},r_{\mathrm{DE}}^{\mathrm{off}},r_{\mathrm{DE}}^{\mathrm{off}}-r_{\mathrm{E}}^{\mathrm{on}}\right]^{\top}\,.$
(13)

The model parameters used in this study are summarized in Table 1.

Table 1: Min model parameters Parameter |  Symbol |  Value
---|---|---
Bulk diffusion |  $D_{c}$ |  $60\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}^{2}\,\mathrm{s}^{-1}$
Membrane diffusion |  $D_{m}$ |  $0.013\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}^{2}\,\mathrm{s}^{-1}$
Average total MinD density |  $\bar{n}_{\mathrm{D}}$ |  $665\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}^{-3}$
Aveage total MinE density |  $\bar{n}_{\mathrm{E}}$ |  $410\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}^{-3}$
Attachment rate |  $k_{\mathrm{D}}$ |  $0.065\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}\,\mathrm{s}^{-1}$
MinD recruitment rate |  $k_{\mathrm{dD}}$ |  $0.098\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}^{3}\,\mathrm{s}^{-1}$
MinE recruitment rate |  $k_{\mathrm{dE}}$ |  $0.126\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}^{3}\,\mathrm{s}^{-1}$
MinDE dissociation rate |  $k_{\mathrm{de}}$ |  $0.34\text{\,}\mathrm{s}^{-1}$
Nucleotide exchange |  $\lambda$ |  $6\text{\,}\mathrm{s}^{-1}$

### Numerical simulations

To investigate the dynamics of the system, we performed 3D finite element
(FEM) simulations using the commercially available software COMSOL
Multiphysics v5.6. Numerical simulations were performed for a wedge geometry
with lateral length $L=$500\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$$
and bulk height $H(x)$ linearly increasing from
$H_{0}=$5\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$$ to
$H_{1}=$50\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$$. The simulation was
initialized with the Min proteins uniformly distributed in the bulk and a
small random spatial perturbation around this uniform state.

### Homogeneous steady state and dispersion relation

The homogeneous steady state concentrations,
($\mathbf{c}^{*}|_{z=0}(H,\bar{n}_{\mathrm{D}},\bar{n}_{\mathrm{E}})$,
$\mathbf{m}^{*}(H,\bar{n}_{\mathrm{D}},\bar{n}_{\mathrm{E}})$) are obtained
from the stationary solutions of Eqs. [1]–[3] together with the mass
conservation condition Eq. (4):

$\left\\{\begin{array}[]{rl}\mathbf{r}(\mathbf{c}^{*}|_{z=0},\mathbf{m}^{*})=&\kern-5.16663pt\mathbf{0},\\\
\boldsymbol{f}\big{(}\mathbf{c}^{*}|_{z=0},\mathbf{m}^{*})=&\kern-5.16663pt\boldsymbol{\Phi},\\\
c^{*}_{\mathrm{D}}|_{z=0}+(m^{*}_{\mathrm{d}}+m^{*}_{\mathrm{de}})/H=&\kern-5.16663pt\bar{n}_{\mathrm{D}},\\\
c^{*}_{\mathrm{E}}|_{z=0}+m^{*}_{\mathrm{de}}/H=&\kern-5.16663pt\bar{n}_{\mathrm{E}},\end{array}\right.$
(14)

where $\boldsymbol{\Phi}$ denotes the steady state fluxes at the membrane,
given by:

$\displaystyle\mathbf{\Phi}$ $\displaystyle=\left[0,\phi,0\right]^{\top},$
(15a) $\displaystyle\phi$
$\displaystyle:=\sqrt{D_{c}\lambda}\,\tanh\left(\sqrt{\lambda/D_{c}}\,H\right)\,c^{*}_{\mathrm{DD}}|_{z=0}.$
(15b)

A concise derivation of these equations and how they can be solved is provided
in the Supplementary Information. For a thorough presentation of the linear
stability analysis of the Min system in a 2D rectangular geometry we refer to
the Supplementary Informations of Refs. [14] and [28].

### Operators for spatial averaging

The operators used throughout this study to calculate mean values of densities
on the membrane and in the cytosol are defined as follows:

$\displaystyle\langle m\rangle_{\mathcal{S}}$
$\displaystyle:=|\mathcal{S}|^{-1}\int_{\mathcal{S}}\mathrm{d}x\mathrm{d}y\,m,$
(16a) $\displaystyle\langle c\rangle_{\mathcal{V}}$
$\displaystyle:=|\mathcal{V}|^{-1}\int_{\mathcal{S}}\mathrm{d}x\mathrm{d}y\int_{0}^{H(x)}\mathrm{d}z\,c,$
(16b) $\displaystyle\langle\cdot\rangle_{y}$
$\displaystyle:=\frac{1}{L}\int_{0}^{L}\mathrm{d}y\,(\cdot),$ (16c)
$\displaystyle\langle\cdot\rangle_{y,z}$
$\displaystyle:=\frac{1}{\,H(x)}\int_{0}^{H(x)}\mathrm{d}z\,\langle\cdot\rangle_{y},$
(16d)

where the membrane surface area and the bulk volume for the wedge geometry are
explicitly given by $|\mathcal{S}|=L^{2}$ and
$|\mathcal{V}|=L^{2}\,(H_{0}+H_{1})/2$.

### Instantaneous total densities _at_ the membrane

Since only cytosolic proteins in close proximity to the membrane participate
in the nonlinear dynamics at the membrane, we define instantaneous total
densities at the membrane:

$\displaystyle\tilde{n}_{\mathrm{D}}(x,y,t)$
$\displaystyle:=\frac{1}{H(x)}(m_{\mathrm{d}}+m_{\mathrm{de}})+c_{\mathrm{D}}|_{z=0}\,,$
(17a) $\displaystyle\tilde{n}_{\mathrm{E}}(x,y,t)$
$\displaystyle:=\frac{1}{H(x)}m_{\mathrm{de}}+c_{\mathrm{E}}|_{z=0}\,.$ (17b)

We further averaged these densities along the $y$–direction to obtain the the
slice-averaged total densities
$\langle\tilde{n}_{\mathrm{D,E}}\rangle_{y}(x,t)$. Note that the length of a
slice is much larger than the typical pattern wavelength, which also permits
to approximate the slice-averaged mass at the membrane by the vertically
averaged mass: $\langle\tilde{n}_{i}\rangle_{y}(x,t)\approx\langle
n_{i}\rangle_{y,z}(x,t)$ (see Ref. [14]). This is because the local deviations
$\tilde{n}_{i}-\langle n_{i}\rangle_{z}$ largely cancel when averaging over
the pattern wavelength.

### Mass redistribution dynamics

Here, we provide more details on the derivation of the mass redistribution
dynamics Eq. (7). For specificity, we present the calculation for MinD. The
calculation for MinE works along the same lines. Our starting point is the
slice averaged total MinD density:

$\langle n_{\mathrm{D}}\rangle_{y,z}(x,t):=\frac{1}{\,H(x)}\left\langle
m_{\mathrm{d}}+m_{\mathrm{de}}+\int_{0}^{H(x)}\mathrm{d}z\,c_{\mathrm{D}}\right\rangle_{y}.$
(18)

The time evolution of this quantity then follows from Eq. (1) and Eq. (2):

$H(x)\,\partial_{t}\langle
n_{\mathrm{D}}\rangle_{y,z}(x,t)=D_{m}\partial_{x}^{2}\langle
m_{\mathrm{d}}+m_{\mathrm{de}}\rangle_{y}\\\ +D_{c}\partial_{z}\langle
c_{\mathrm{D}}\rangle_{y}\big{|}_{z=H(x)}+\int_{0}^{H(x)}\mathrm{d}z\,D_{c}\partial_{x}^{2}\langle
c_{\mathrm{D}}\rangle_{y},$ (19)

where we used the reactive boundary condition Eq. (3) to rewrite the integral:

$\displaystyle\int_{0}^{H(x)}\mathrm{d}z\,D_{c}\partial_{z}^{2}c_{\mathrm{D}}$
$\displaystyle=D_{c}\partial_{z}c_{\mathrm{D}}\big{|}_{z=H(x)}-D_{c}\partial_{z}c_{\mathrm{D}}\big{|}_{z=0}$
$\displaystyle=D_{c}\partial_{z}c_{\mathrm{D}}\big{|}_{z=H(x)}+r_{\mathrm{DE}}^{\mathrm{off}}-r_{\mathrm{D}}^{\mathrm{on}}.$
(20)

Note that due to mass-conservation the reaction terms at the membrane cancel.

Since the system is closed, the boundary condition at the inclined top surface
of the wedge reads $\mathbf{n}\cdot\nabla c_{\mathrm{D}}|_{z=H(x)}=0$, where
$\mathbf{n}\propto(-\partial_{x}H,0,1)$ is the outward normal vector at the
top surface. Writing out the boundary condition explicitly, we find that:

$\partial_{z}c_{\mathrm{D}}|_{z=H(x)}=(\partial_{x}H)\,\partial_{x}c_{\mathrm{D}}|_{z=H(x)}.$
(21)

To proceed, we substitute the relation above into Eq. (19) and slightly
rewrite the resulting equation by applying the chain rule:

$H(x)\,\partial_{t}\langle
n_{\mathrm{D}}\rangle_{y,z}(x,t)=D_{m}\partial_{x}^{2}\langle
m_{\mathrm{d}}+m_{\mathrm{de}}\rangle_{y}+\partial_{x}\underbrace{\int_{0}^{H(x)}\mathrm{d}z\,D_{c}\partial_{x}\langle
c_{\mathrm{D}}\rangle_{y}}_{\displaystyle=:-H(x)J_{\mathrm{D}}(x)}.$ (22)

Here, the first term describes diffusion of the averaged membrane
concentrations. The integral on the right describes diffusion of the averaged
cytosolic densities, where we defined the diffusive flux
$J_{\mathrm{D}}=-D_{c}\langle\partial_{x}c_{D}\rangle_{y,z}$. The factor
$H(x)$ in the cytosolic diffusion term accounts for the increasing area of the
slice along the positive $x$–direction.

Since protein diffusion on the membrane is much smaller than cytosolic
diffusion $D_{m}\ll D_{c}$ [47, 48], one can neglect membrane diffusion to
arrive at the result shown in the main text (Eq. (7)). For completeness, note
that Eq. (22) (without membrane diffusion) can be recast as

$\displaystyle\partial_{t}\langle n_{\mathrm{D}}\rangle_{y,z}(x,t)$
$\displaystyle\approx\frac{1}{H(x)}\,\partial_{x}\int_{0}^{H(x)}\mathrm{d}z\,D_{c}\partial_{x}\langle
c_{\mathrm{D}}\rangle_{y},$
$\displaystyle=D_{c}\partial_{x}\langle\partial_{x}c_{\mathrm{D}}\rangle_{y,z}+D_{c}\frac{\partial_{x}H(x)}{H(x)}\langle\partial_{x}c_{\mathrm{D}}\rangle_{y,z},$
(23)

which is the form given in Eq. (6) in the main text.

### Machine-learning based pattern classification

We used the pixel classifier provided by the software _ilastik_ [35]. The
classifier was trained based on a few representative snapshots, by manually
marking areas where the pattern type (no pattern, chaos, standing wave, or
traveling wave) is easily identified by visual inspection. The trained
classifier then yields probabilities for each pattern type at each pixel. The
classifier was applied to snapshots in 20 s intervals. This data was then
downsampled and averaged over slices to yield an $x$–$t$ space time map of
pattern probabilities. To render the kymograph in Fig. 4D each pixel was
colored based on the most probable pattern.

### Aranson–Tsimring model

As a second example we briefly discuss a phenomenological model for pattern
formation in vibrated granular media introduced in [39]. This model, which we
call Aranson–Tsimring model in the following, couples a Ginzburg–Landau-type
equation [34] for the complex order parameter $\psi$ to a conservation law for
the density $\rho$:

$\displaystyle\partial_{t}\psi$
$\displaystyle=\gamma\bar{\psi}-(1-i\omega)\psi+(1+ib)\nabla^{2}\psi-|\psi|^{2}\psi-\rho\psi,$
(24a) $\displaystyle\partial_{t}\rho$
$\displaystyle=\beta\nabla^{2}\rho+\alpha\nabla\cdot(\rho\nabla|\psi|^{2}),$
(24b)

where $\bar{\psi}$ denotes the complex conjugate of $\psi$. The coupling is
such that increasing the density $\rho$ suppresses the instability in Eq.
(24a) while gradients in the amplitude $|\psi|$ drive mass redistribution away
from high amplitude regions (second term in Eq. (24b)). This feedback loop
amplifies heterogeneities in the density and gives rise to localized patterns.
These patterns have been studied in detail using amplitude equations in [32,
33]. Moreover, in Ref. [30] it was shown that the system Eq. (24) appears as
the amplitude equation for a mass-conserving version of the classical
Swift–Hohenberg–Turing equation [6, 49]. The reason for this is that the
conserved density appears as a second hydrodynamic variable in addition to the
pattern amplitude.

A linear stability analysis shows that the system Eq. (24) has a short
wavelength instability when $b\omega-1-\rho_{0}>0$ and
$\gamma>\gamma_{c}=(\omega+b(1+\rho_{0}))/\sqrt{1+b^{2}}$, where $\rho_{0}$
denotes the average density. Following Ref. [33], we set parameters
$b=1,\omega=2.5,\alpha=1.3,\beta=0.3,\rho_{0}=0.3$. Localized patterns are
found near the instability threshold, so we set $\gamma=1.001\gamma_{c}$ for
the simulation shown in Fig. 6 and Movie S9.

### Preparation of the wedge flow cell

The microfluidic wedge chambers were prepared using two rectangular cover
slips (bottom one of dimensions $22$/$50\text{\,}\mathrm{m}\mathrm{m}$, and
top one of dimensions $5$/$30\text{\,}\mathrm{m}\mathrm{m}$). Close to one of
the short edges of a top glass a tiny inlet hole was drilled using a
sandblaster. Cover slips were cleaned in $1$ M KOH for $1$ h followed by a
methanol bath for $10$ min in a sonicator bath. Surfaces of the cover slips
were activated with oxygen plasma for $20$ s, using oxygen plasma PREEN I
(Plasmatic System, Inc.) with a O2 flow rate of $1$ SCFH. Furthermore, a small
PDMS slab with a $0.3\text{\,}\mathrm{m}\mathrm{m}$ hole was attached on to
the top glass slide, such that it matches the hole in the PDMS glass slide and
a metal connector was inserted in the hole for connecting the syringe pump.
Tilt of the top glass slide was achieved by placing a piece of aluminum foil
between the top and bottom slide at the end, with the largest height between
top and bottom at the side of the inlet. At the opposite side with the
smallest distance between top and bottom slide,
$2\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$ polystyrene beads that were
deposited on the bottom slide provided an outlet and prevent a collapse of the
top and bottom slides (see Fig. 2). The lateral sides of the microchamber were
sealed with a two-component epoxy resin leaving the short edge at the low
height-side open for liquid flow (Fig. S4). The microfluidic cell was then
filled with a solution of small unilamellar vesicles (SUVs) through an
injection tube at the inlet of the PDMS slab and incubated for 30 min at
$30\text{\,}\mathrm{\SIUnitSymbolCelsius}$–yielding full lipid membrane
coverage of the bottom and top slides. SUVs were prepared as described in Ref.
[28]. Subsequently, the flow cell was thoroughly washed with a buffer to
remove excess SUVs and Min protein experiments were started.

### Observation of Min patterns

We purified the Min proteins based on the method proposed in Ref. [50].
Injection of Min proteins into the flow cells was performed through a syringe
pump containing a solution of $0.8\text{\,}\mathrm{M}$ MinD,
$0.2\text{\,}\mathrm{m}\mathrm{M}$ MinD-Cy3,
$0.8\text{\,}\mathrm{m}\mathrm{M}$ MinE, $0.2\text{\,}\mathrm{m}\mathrm{M}$
MinE-Cy5, $5\text{\,}\mathrm{m}\mathrm{M}$ ATP,
$4\text{\,}\mathrm{m}\mathrm{M}$ phosphoenolpyruvate,
$0.01\text{\,}\mathrm{m}\mathrm{g}\mathrm{/}\mathrm{m}\mathrm{l}$ pyruvate
kinase, $25\text{\,}\mathrm{m}\mathrm{M}$ Tris-HCl (pH 7.5),
$150\text{\,}\mathrm{m}\mathrm{M}$ KCl and $5\text{\,}\mathrm{m}\mathrm{M}$
MgCl2. To ensure that all of the buffer solution in the microdevice is
replaced by the protein solution, we chose a volume of the protein solution
that was 50 times larger than the volume in the microdevice. During the
filling process of the microdevice, the enire solution was rapidly injected
(in 5 s) to prevent protein accumulation on the membrane.

For the generation of the fluorescence images, we used the following
equipment: Olympus IX-81 inverted microscope equipped with an Andor Revolution
XD spinning disk system with FRAPPA, illumination and detection system Andor
Revolution and Yokogawa CSU X1, EM-CCD Andor iXon X3 DU897 camera, motorized
x-y stage and a z-piezo stage, using a 20x objective (UPlansApo, NA 0.85, oil
immersion). Imaging of MinD-Cy3 and MinE-Cy5 was performed with laser spectral
lines at 561 nm and 640 nm, respectively, and we further used a 617/73 band-
pass filter as well as a 690 long-pass filter. We imaged several uniformly
sized regions at intervals of 30 s or 60 s along the lateral length of the
wedge setup. To exclude membrane imperfections that may have arisen during
preparation, we also imaged the membrane using the spectral line at
$491\text{\,}\mathrm{n}\mathrm{m}$ and a 525/50 band-pass filter.

### Image sequence processing

We processed the fluorescence images using the following software packages:
Andor iQ3 v3.1, ImageJ 1.52j, and custom written Matlab 2016a scripts. For
better visualization, we additionally applied background correction and
filtering of artifacts. In detail, these were carried out as follows: For the
generation of the movies, each frame was first corrected for fluorescence
bleaching (max. 20 % decay of the intensity for long movies) by normalizing to
the mean intensity of the respective frame. Then, we generated two different
modifications of the images: First, we averaged out all transient features
(i.e., patterns) in the frames to obtain ‘static background’-images which we
shall call Imstat. Second, we smoothed out the images, determined the average
of all movie frames, and normalized the corresponding result with respect to
its maximum. This way, we obtained an ‘illumination correction’ image Imillum.
In the final step, each frame Immovie was corrected according to the rule
Imcorrected = (Immovie \- Imstat)/Imillum. On one hand, this ensures that
irregularities in each image are suppressed, and on the other hand, the
intensity amplitudes at the edges becomes comparable with the values at the
center of the image.

## References

  * [1] Alan M. Turing.  The chemical basis of morphogenesis.  Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences, 237(641):37–72, 1952.
  * [2] Shigeru Kondo and Rihito Asai.  A reaction–diffusion wave on the skin of the marine angelfish Pomacanthus.  Nature, 376(6543):765–768, 1995.
  * [3] Christopher A. Klausmeier.  Regular and irregular patterns in semiarid vegetation.  Science (New York, N.Y.), 284(5421):1826–1828, 1999.
  * [4] J. B. A. Green and J. Sharpe.  Positional information and reaction-diffusion: Two big ideas in developmental biology combine.  Development (Cambridge, England), 142(7):1203–1211, 2015.
  * [5] Jacob Halatek, Fridtjof Brauns, and Erwin Frey.  Self-organization principles of intracellular pattern formation.  Philosophical Transactions of the Royal Society B: Biological Sciences, 373(1747):20170107, 2018.
  * [6] M. C. Cross and P. C. Hohenberg.  Pattern formation outside of equilibrium.  Reviews of Modern Physics, 65(3):851–1112, 1993.
  * [7] N. Tompkins, N. Li, C. Girabawe, M. Heymann, G. B. Ermentrout, I. R. Epstein, and S. Fraden.  Testing Turing’s theory of morphogenesis in chemical cells.  Proceedings of the National Academy of Sciences, 111(12):4397–4402, 2014.
  * [8] Erwin Frey and Fridtjof Brauns.  Self-organisation of protein patterns.  arXiv:2012.01797 [nlin, physics:physics], 2020.
  * [9] Andrew B. Goryachev and Marcin Leda.  Many roads to symmetry breaking: Molecular mechanisms and theoretical models of yeast cell polarity.  Molecular Biology of the Cell, 28(3):370–380, 2017.
  * [10] Philipp Glock, Fridtjof Brauns, Jacob Halatek, Erwin Frey, and Petra Schwille.  Design of biochemical pattern forming systems from minimal motifs.  eLife, 8:e48646, 2019.
  * [11] Manon C. Wigbers, Tzer Han Tan, Fridtjof Brauns, Jinghui Liu, S. Zachary Swartz, Erwin Frey, and Nikta Fakhri.  A hierarchy of protein patterns robustly decodes cell shape information.  Nature Physics, pages 1–7, 2021.
  * [12] Uwe C. Täuber.  Critical Dynamics: A Field Theory Approach to Equilibrium and Non-Equilibrium Scaling Behavior.  2014\.
  * [13] Robert Zwanzig.  Nonequilibrium Statistical Mechanics.  Oxford University Press, Oxford ; New York, 2001.
  * [14] J. Halatek and E. Frey.  Rethinking pattern formation in reaction–diffusion systems.  Nature Physics, 14(5):507, 2018.
  * [15] Fridtjof Brauns, Jacob Halatek, and Erwin Frey.  Phase-space geometry of mass-conserving reaction-diffusion dynamics.  Physical Review X, 10(4):041036, 2020.
  * [16] H. I. Adler, W. D. Fisher, A. Cohen, and A. A. Hardigree.  Miniature escherichia coli cells deficient in DNA.  Proceedings of the National Academy of Sciences, 57(2):321–326, 1967.
  * [17] Piet A. J. de Boer, Robin E. Crossley, and Lawrence I. Rothfield.  A division inhibitor and a topological specificity factor coded for by the minicell locus determine proper placement of the division septum in E. coli.  Cell, 56(4):641–649, 1989.
  * [18] M. Loose, E. Fischer-Friedrich, J. Ries, K. Kruse, and P. Schwille.  Spatial regulators for bacterial cell division self-organize into surface waves in vitro.  Science (New York, N.Y.), 320(5877):789–792, 2008.
  * [19] Erwin Frey, Jacob Halatek, Simon Kretschmer, and Petra Schwille.  Protein pattern formation.  In Patricia Bassereau and Pierre Sens, editors, Physics of Biological Membranes, pages 229–260. Springer International Publishing, Cham, 2018.
  * [20] Beatrice Ramm, Tamara Heermann, and Petra Schwille.  The E. coli MinCDE system in the regulation of protein patterns and gradients.  Cellular and Molecular Life Sciences, 76(21):4245–4273, 2019.
  * [21] V. Ivanov and K. Mizuuchi.  Multiple modes of interconverting dynamic pattern formation by bacterial cell division proteins.  Proceedings of the National Academy of Sciences, 107(18):8071–8078, 2010.
  * [22] Anthony G. Vecchiarelli, Min Li, Michiyo Mizuuchi, Ling Chin Hwang, Yeonee Seol, Keir C. Neuman, and Kiyoshi Mizuuchi.  Membrane-bound MinDE complex acts as a toggle switch that drives Min oscillation coupled to cytoplasmic depletion of MinD.  Proceedings of the National Academy of Sciences, 113(11):E1479–E1488, 2016.
  * [23] Jonas Denk, Simon Kretschmer, Jacob Halatek, Caroline Hartl, Petra Schwille, and Erwin Frey.  MinE conformational switching confers robustness on self-organized Min protein patterns.  Proceedings of the National Academy of Sciences, 115(18):4553–4558, 2018.
  * [24] Thomas Litschel, Beatrice Ramm, Roel Maas, Michael Heymann, and Petra Schwille.  Beating vesicles: Encapsulated protein oscillations cause dynamic membrane deformations.  Angewandte Chemie International Edition, 57(50):16286–16290, 2018\.
  * [25] Elisa Godino, Jonás Noguera López, David Foschepoth, Céline Cleij, Anne Doerr, Clara Ferrer Castellà, and Christophe Danelon.  De novo synthesized Min proteins drive oscillatory liposome deformation and regulate FtsA-FtsZ cytoskeletal patterns.  Nature Communications, 10(1):4969, 2019.
  * [26] Shunshi Kohyama, Natsuhiko Yoshinaga, Miho Yanagisawa, Kei Fujiwara, and Nobuhide Doi.  Cell-sized confinement controls generation and stability of a protein wave for spatiotemporal regulation in cells.  eLife, 8:e44591, 2019.
  * [27] Philipp Glock, Beatrice Ramm, Tamara Heermann, Simon Kretschmer, Jakob Schweizer, Jonas Mücksch, Gökberk Alagöz, and Petra Schwille.  Stationary patterns in a two-protein reaction-diffusion system.  ACS Synthetic Biology, 8(1):148–157, 2019.
  * [28] Fridtjof Brauns, Grzegorz Pawlik, Jacob Halatek, Jacob Kerssemakers, Erwin Frey, and Cees Dekker.  Bulk-surface coupling identifies the mechanistic connection between Min-protein patterns in vivo and in vitro.  Nature Communications, 12(1), 2021.
  * [29] Anthony G. Vecchiarelli, Min Li, Michiyo Mizuuchi, and Kiyoshi Mizuuchi.  Differential affinities of MinD and MinE to anionic phospholipid influence Min patterning dynamicsin vitro.  Molecular Microbiology, 93(3):453–463, 2014.
  * [30] P. C. Matthews and S. M. Cox.  Pattern formation with a conservation law.  Nonlinearity, 13(4):1293–1320, 2000.
  * [31] S. M. Cox and P. C. Matthews.  Instability and localisation of patterns due to a conserved quantity.  Physica D: Nonlinear Phenomena, 175(3-4):196–219, 2003.
  * [32] D. M. Winterbottom, P. C. Matthews, and S. M. Cox.  Oscillatory pattern formation with a conserved quantity.  Nonlinearity, 18(3):1031–1056, 2005.
  * [33] D. M. Winterbottom, S. M. Cox, and P. C. Matthews.  Pattern formation in a model of a vibrated granular layer.  SIAM Journal on Applied Dynamical Systems, 7(1):63–78, 2008.
  * [34] Igor S. Aranson and Lorenz Kramer.  The world of the complex Ginzburg-Landau equation.  Reviews of Modern Physics, 74(1):99–143, 2002.
  * [35] Stuart Berg, Dominik Kutra, Thorben Kroeger, Christoph N. Straehle, Bernhard X. Kausler, Carsten Haubold, Martin Schiegg, Janez Ales, Thorsten Beier, Markus Rudy, Kemal Eren, Jaime I. Cervantes, Buote Xu, Fynn Beuttenmueller, Adrian Wolny, Chong Zhang, Ullrich Koethe, Fred A. Hamprecht, and Anna Kreshuk.  Ilastik: Interactive machine learning for (bio)image analysis.  Nature Methods, 16(12):1226–1232, 2019.
  * [36] Fridtjof Brauns, Jacob Halatek, and Erwin Frey.  Diffusive coupling of two well-mixed compartments elucidates elementary principles of protein-based pattern formation.  Physical Review Research, 3(1), 2021.
  * [37] Mariela D. Petkova, Gašper Tkačik, William Bialek, Eric F. Wieschaus, and Thomas Gregor.  Optimal decoding of cellular identities in a genetic network.  Cell, 176(4):844–855.e15, 2019.
  * [38] Fridtjof Brauns, Henrik Weyer, Jacob Halatek, Junghoon Yoon, and Erwin Frey.  Wavelength selection by interrupted coarsening in reaction-diffusion systems.  Physical Review Letters, 126(10), 2021.
  * [39] Lev S. Tsimring and Igor S. Aranson.  Localized and cellular patterns in a vibrated granular layer.  Physical Review Letters, 79(2):213–216, 1997.
  * [40] M. C. Marchetti, J. F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, Madan Rao, and R. Aditi Simha.  Hydrodynamics of soft active matter.  Reviews of Modern Physics, 85(3):1143–1189, 2013.
  * [41] Markus Bär, Robert Großmann, Sebastian Heidenreich, and Fernando Peruani.  Self-propelled rods: Insights and perspectives for active matter.  Annual Review of Condensed Matter Physics, 11(1):441–466, 2020\.
  * [42] Jonas Denk and Erwin Frey.  Pattern-induced local symmetry breaking in active-matter systems.  Proceedings of the National Academy of Sciences, 117(50):31623–31630, 2020.
  * [43] Olivier Dauchot and Hartmut Löwen.  Chemical physics of active matter.  The Journal of Chemical Physics, 151(11), 2019.
  * [44] Robert Großmann, Igor S. Aranson, and Fernando Peruani.  A particle-field approach bridges phase separation and collective motion in active matter.  Nature Communications, 11(1), 2020.
  * [45] K. C. Huang, Y. Meir, and N. S. Wingreen.  Dynamic structures in Escherichia coli: Spontaneous formation of MinE rings and MinD polar zones.  Proceedings of the National Academy of Sciences, 100(22):12724–12728, 2003.
  * [46] Jacob Halatek and Erwin Frey.  Highly canalized MinD transfer and MinE sequestration explain the origin of robust MinCDE-Protein dynamics.  Cell Reports, 1(6):741–752, 2012.
  * [47] Martin Loose, Karsten Kruse, and Petra Schwille.  Protein self-organization: Lessons from the min system.  Annual Review of Biophysics, 40(1):315–336, 2011.
  * [48] G. Meacci, J. Ries, E. Fischer-Friedrich, N. Kahya, P. Schwille, and K. Kruse.  Mobility of Min-proteins in ¡i¿Escherichia coli¡/i¿ measured by fluorescence correlation spectroscopy.  Physical Biology, 3(4):255–263, 2006.
  * [49] Jonathan H. P. Dawes.  After 1952: The later development of Alan Turing’s ideas on the mathematics of pattern formation.  Historia Mathematica, 43(1):49–64, 2016.
  * [50] Yaron Caspi and Cees Dekker.  Mapping out Min protein patterns in fully confined fluidic chambers.  eLife, 5, 2016.