10

††thanks: These three authors contributed equally††thanks: These three authors
contributed equally††thanks: These three authors contributed equally

# Pattern localization to a domain edge

Manon C. Wigbers  Fridtjof Brauns  Tobias Hermann  Erwin Frey<EMAIL_ADDRESS>Arnold Sommerfeld Center for Theoretical Physics (ASC) and Center for
NanoScience (CeNS), Department of Physics, Ludwig-Maximilians-Universität
München, Theresienstraße 37, D–80333 München, Germany

###### Abstract

The formation of protein patterns inside cells is generically described by
reaction–diffusion models. The study of such systems goes back to Turing, who
showed how patterns can emerge from a homogenous steady state when two
reactive components have different diffusivities (e.g. membrane-bound and
cytosolic states). However, in nature, systems typically develop in a
heterogeneous environment, where upstream protein patterns affect the
formation of protein patterns downstream. Examples for this are the
polarization of Cdc42 adjacent to the previous bud-site in budding yeast, and
the formation of an actin-recruiter ring that forms around a PIP3 domain in
macropinocytosis. This suggests that previously established protein patterns
can serve as a template for downstream proteins and that these downstream
proteins can ‘sense’ the edge of the template. A mechanism for how this edge
sensing may work remains elusive.

Here we demonstrate and analyze a generic and robust edge-sensing mechanism,
based on a two-component mass-conserving reaction-diffusion (McRD) model. Our
analysis is rooted in a recently developed theoretical framework for McRD
systems, termed local equilibria theory. We extend this framework to capture
the spatially heterogeneous reaction kinetics due to the template. This
enables us to graphically construct the stationary patterns in the phase space
of the reaction kinetics. Furthermore, we show that the protein template can
trigger a regional mass-redistribution instability near the template edge,
leading to the accumulation of protein mass, which eventually results in a
stationary peak at the template edge. We show that simple geometric criteria
on the reactive nullcline’s shape predict when this edge-sensing mechanism is
operational. Thus, our results provide guidance for future studies of
biological systems, and for the design of synthetic pattern forming systems.

###### pacs:

##  I Introduction

###  I.1 Background and motivation

Many cellular processes, such as cell division and cell motility, rely
crucially on the localization of proteins in space and time. Strikingly, these
protein localization patterns can emerge from the collective coordination of
transport and local molecular interactions of proteins. Diffusion in the
cytosol is a simple means of protein transport that accounts for many self-
organization processes Halatek  _et al._ (2018). To analyze how the interplay
of diffusive protein transport and protein-protein interactions on a nanometer
scale influences the protein patterns on the cellular scale, mass-conserving
reaction–diffusion models have proven useful Howard  _et al._ (2001); Huang
_et al._ (2003); Mori  _et al._ (2008); Goryachev and Pokhilko (2008);
Ishihara  _et al._ (2007); Otsuji  _et al._ (2007); Halatek and Frey (2012);
Klünder  _et al._ (2013); Trong  _et al._ (2014); Alonso and Bär (2014); Wu
_et al._ (2016); Goryachev and Leda (2017); Murray and Sourjik (2017); Denk
_et al._ (2018); Cusseddu  _et al._ (2018); Chiou  _et al._ (2018); Glock  _et
al._ (2019). The study of reaction–diffusion systems in general goes back to
Turing Turing (1952), who showed how patterns can emerge from a homogenous
steady state when two reactive components have different diffusivities. In
cells, differential diffusivities are generic because many proteins have
membrane-bound and cytosolic states, where diffusion on the membrane is orders
of magnitude slower than in the cytoplasm. Turing’s pioneering work Turing
(1952) has led to vast advances in the field on how protein patterns arise
from homogeneous (initial) steady states on spatially homogeneous domains.
However, as Turing already pointed out Turing (1952), “most of an organism,
most of the time, is developing from one pattern into another, rather than
from homogeneity into a pattern.”

For example, previously formed protein patterns can control pattern formation
of proteins downstream by affecting their local interactions, such that the
upstream pattern acts as a _spatial template_ for the downstream proteins. A
biological system where such “templating” has been suggested is
macropinocytosis Williams  _et al._ (2019). Here, a high density domain of
PIP3 (a charged phospholipid) and a Ras-GTPase111GTPases are hydrolase enzymes
that can bind and hydrolyze guanosine triphosphate (GTP). Ras is a subfamily
of small GTPases. have been suggested to serve as a template for a ring of
actin recruiters (SCAR complex, Arp2/3), that forms around the PIP3 domain
edge Veltman  _et al._ (2016). Recruitment of an actomyosin ring, controlled
by GTPases, is also key for single-cell wound healing. Following the rupture
of the cell wall, two GTPases– Abr and Cdc42 –are recruited to the wound edge,
where they organize into two concentric rings of high protein concentration
Vaughan  _et al._ (2011). Cdc42 in turn recruits actomyosin which contracts to
close the wound and repair the underlying cytoskeleton. Mutations of Abr,
which forms the inner ring, leads to a loss of the outside Cdc42 ring,
suggesting hierarchical interaction between Abr and Cdc42 Vaughan  _et al._
(2011). Thus, the inner Abr-ring may be pictured as a template for the outer
Cdc42-ring. Yet another example where protein patterns act as a spatial
template, can be found during cell division in budding yeast. Here, landmark
proteins direct the polarization of the GTPase Cdc42, such that the Cdc42
cluster emerges either adjacent to the previous bud-site, or at the opposite
cell pole, depending on the cell-type Chant and Pringle (1995). Various
mutations or deletions of individual landmark proteins lead to Cdc42 clusters
right on top of the previous bud-site or at a random position Tong  _et al._
(2007); Lo  _et al._ (2013); Miller  _et al._ (2017). Hence, the landmark
proteins may be pictured as a template that controls Cdc42 pattern formation.
Common to all the above examples is that the downstream proteins localize at
the edge of some template. Both the specific (molecular) mechanisms, and the
general principles underlying this ‘edge sensing’ remain elusive.

Here, we present a pattern-forming mechanism capable of robust edge sensing
and provide criteria for its operation based on simple geometric relations in
the phase space of the reaction kinetics. To find these criteria, we use a
recently developed framework, termed _local equilibria theory_ which enables
us to gain insight into the dynamics of _mass-conserving reaction–diffusion_
(McRD) systems Halatek and Frey (2018); Brauns  _et al._ (2018). We briefly
review the key elements of the local equilibria theory for a paradigmatic
model for cell polarization in Section I.2. We then introduce a step-like
template that imposes heterogeneity in the reaction kinetics, and generalize
the framework to study the dynamics of such systems. This enables us to
explain why and under which conditions a density peak forms at the edge of the
template. Thus, our results may provide guidance for the design of patterns in
synthetic systems and may help to identify molecular mechanisms underlying
edge-sensing in biological systems.

###  I.2 Local equilibria theory

Figure 1:  Illustration of the phase-space geometric analysis for two-
component McRD systems. (a) The reactive equilibria (black dots) are given by
the intersections between the reactive subspaces $m+c=n$ (grey lines) and the
reactive nullcline $f=0$ (black line) in the $(m,c)$-phase space. Hence, the
reactive nullcline encodes the qualitative structure of the reactive flow as
illustrated by the red arrows. (b) Sketch of the membrane profile of a mesa
pattern composed of a high- and low-density domains, $m^{+}$ and $m^{-}$,
connected by a diffusive interface around the inflection point $x_{0}$. (c)
Flux-balance construction of the mesa pattern in phase space (cf. (a)). The
intersections of the flux-balance subspace (FBS) (purple dashed line) and the
reactive nullcline yield the concentrations at the plateaus $m^{\pm}$ and the
inflection point $m^{0}$. The balance of net reactive flows in the system (red
arrows) determines the FBS-offset $\eta_{0}$. In the regime where the slope of
the reactive nullcline is steeper than the slope of the FBS, an homogenous
steady state is laterally unstable. (d) Linearization of the dynamics in the
vicinity of the homogeneous steady state, yields a dispersion relation for
growth rates of the eigenfunctions (Fourier modes indexed by wavenumber $q$).
The fastest growing mode, $q_{c}$, dominates the length scale of the initial
dynamics.

We consider the dynamics of one protein species on a one-dimensional domain of
length $L$, as in Ref. Brauns  _et al._ (2018). The proteins can cycle between
a membrane-bound state (concentration $m(x,t)$) and a cytosolic state
(concentration $c(x,t)$), with diffusion constants $D_{m}$ and $D_{c}$
respectively. In cells, the diffusion constants of membrane-bound proteins and
cytosolic proteins are typically widely different, such that $D_{m}\ll D_{c}$.
The reaction–diffusion equations for the membrane density $m$ and the
cytosolic density $c$ read

$\displaystyle\partial_{t}m(x,t)$ $\displaystyle=D_{m}$
$\displaystyle\partial_{x}^{2}\,m$ $\displaystyle+f(m,c),$ (1a)
$\displaystyle\partial_{t}c(x,t)$ $\displaystyle=\;D_{c}$
$\displaystyle\partial_{x}^{2}\,c$ $\displaystyle-f(m,c),$ (1b)

where the reaction term $f(m,c)$ describes the attachment–detachment dynamics
of the proteins. Specific examples of such systems exhibiting self-organized
pattern formation can be found in Ishihara  _et al._ (2007); Otsuji  _et al._
(2007); Chiou  _et al._ (2018); Mori  _et al._ (2008). At the boundaries, we
impose no-flux conditions
$D_{c}\partial_{x}c|_{0,L}=D_{m}\partial_{x}m|_{0,L}=0$. The dynamics
conserves total protein density

$\bar{n}=\frac{1}{L}\int_{0}^{L}\\!\mathrm{d}x\;n(x,t),$ (2)

with the local total density $n(x,t)=m(x,t)+c(x,t)$.

To characterize the dynamics and steady states of McRD systems, we recently
introduced a framework, termed _local equilibria theory_ Halatek and Frey
(2018); Brauns  _et al._ (2018). This theory proposes to analyze spatially
extended systems as a collection of small diffusively coupled compartments.
The local reaction kinetics inside each of the compartments, then serves a
proxy for the spatially extended dynamics, enabling a quantitative phase
portrait analysis of the spatially extended system in the phase space of
reaction kinetics Brauns  _et al._ (2018). In the following we briefly review
the key results of local equilibria theory for the two-component McRD system
and generalize this framework to analyze pattern formation in the presence of
a spatial template. For a comprehensive analysis of the two-component McRD
system on a homogeneous domain, we refer to Ref. Brauns  _et al._ (2018).

The reaction kinetics of McRD systems conserves total protein mass, which
implies that the reactive flow must point along the reactive phase spaces
$n=c+m$, indicated by the gray lines in Fig. 1(a). The reactive flow vanishes
along the _reactive nullcline_ (NC), given by $f(m,c)=0$. Intersections of the
reactive nullcline with reactive phase spaces, given by the total density
(mass) $n$, determine the reactive equilibria $(m^{*}(n),c^{*}(n))$ shown as
black dots in Fig. 1(a). Hence, the shape of the nullcline encodes how the
reactive equilibria move when total density $n$ is changed, highlighting that
the total density $n$ is a _control parameter_ for the reaction dynamics.
Within each reactive phase space, the flow is directed towards a stable
reactive equilibrium, as illustrated by the red arrows in Fig. 1(a).

In a spatially extended system, the total density $n(x,t)=m(x,t)+c(x,t)$ is
generically inhomogeneous, and its dynamics is driven by diffusion, as can be
seen by adding Eqs. (1a) and (1b)

$\partial_{t}n(x,t)=D_{c}\partial_{x}^{2}\eta(x,t),$ (3)

where we introduced the _mass-redistribution potential_ , defined as Brauns
_et al._ (2018)

$\eta(x,t):=c(x,t)+\frac{D_{m}}{D_{c}}m(x,t).$ (4)

To study the interplay of local reactions and diffusive mass-transport in
spatially extended systems, local equilibria theory proposes to analyze such
systems as a collection of diffusively coupled compartments. These notional
compartments are chosen small enough that each of them can be regarded as
well-mixed. Thus, local dynamics within each compartment can be characterized
in the ODE phase space of reactions which is determined by the density
$n(x,t)$ within that compartment. In this characterization, the _local
(reactive) equilibria_ and their stability in each local phase space serve as
proxies for the reactive dynamics in each compartment. This becomes clear when
one imagines the compartments as isolated, for a given total density profile
$n(x)$. Then each compartment will approach a stable local equilibrium,
parametrized by the local density $n(x)$. In the spatially coupled system, the
total density $n(x,t)$ is diffusively redistributed due to concentration
gradients between the compartments (cf. Eq. (3)). Consequently, the local
equilibria shift and their stability may change Halatek  _et al._ (2018);
Brauns  _et al._ (2018). This interplay between shifting local equilibria and
mass transport is at the core of _local equilibria theory_.

In the remainder of this section, we recapitulate two key results from the
phase-portrait analysis of two-component McRD systems Brauns  _et al._ (2018).
We will later generalize this analysis to systems on a spatially heterogeneous
domain.

Flux-balance construction. — From the dynamics of the total mass Eq. (3)–(4)
it follows that, for any stationary pattern (denoted by $\widetilde{m}(x)$,
$\widetilde{c}(x)$), $\widetilde{\eta}(x)$ must be constant in space on a
domain with no-flux (or periodic) boundary conditions Brauns  _et al._ (2018):

$\eta_{0}=\widetilde{c}(x)+\frac{D_{m}}{D_{c}}\widetilde{m}(x)=\mathrm{const}.$
(5)

This relationship defines a linear subspace, termed _flux-balance subspace_
(FBS), of the $(m,c)$-phase space of reaction kinetics (purple dashed line in
Fig. 1(c)). Any stationary pattern must be embedded in a single FBS. This
reflects that, in steady state, the diffusive fluxes in $m$ and $c$ are
balanced against each other such that there is no net transport of mass.

We can use this condition, Eq. (5), to geometrically construct the steady
state density profile in the $(m,c)$-phase space and from that estimate the
real space density profile. The key insight is that we can approximate the
concentrations at the plateaus, and the inflection point of the pattern by the
local equilibria at the FBS-NC intersections (see Fig. 1(b)). We denote these
intersection points by $m^{-}$, $m^{0}$ and $m^{+}$, where $m^{\pm}$
correspond to the concentrations at the plateaus and $m^{0}$ to the
concentration at the inflection point of the pattern (Fig. 1(b,c)). Thus, the
FBS-offset, $\eta_{0}$, fully determines these concentrations.

To determine the FBS-offset, $\eta_{0}$, one uses that in steady state the net
reactive flow within the whole system must be balanced222Note that before
integration, Eq. (1a) is multiplied with $\partial_{x}\widetilde{m}(x)$ as a
mathematical trick to substitute the integral over space by an integral over
$m$.:

$\int_{m^{-}(\eta_{0})}^{m^{+}(\eta_{0})}\mathrm{d}m\,f\left(m,\eta_{0}-\frac{D_{m}}{D_{c}}m\right)=0,$
(6)

where the plateau concentrations far away from the interface are approximated
by the FBS-NC intersections $m^{\pm}(\eta_{0})$. This _total turnover balance_
condition implicitly determines the FBS-offset $\eta_{0}$. Note that on a
large domain (much larger than the interface width, where the approximation
$m(0,L)\approx m_{\pm}$ holds), total turnover Eq. (6), and hence $\eta_{0}$,
depends only on the function $f$ and the ratio of the diffusion constants.
This implies that $\eta_{0}$ is not dependent on the average mass $\bar{n}$ in
this approximation.

We will show next that the average mass $\bar{n}$ determines the relative size
of the low- and high-density regions and with that the position of the
pattern’s interface. This interface is marked by the position of the
inflection point $x_{0}$ of the pattern profile. For a domain size much larger
than the interface width, we can neglect the finite width of the interface
region, such that the average mass can be approximated by

$L\,\bar{n}\approx x_{0}\,n^{-}(\eta_{0})+(L-x_{0})\,n^{+}(\eta_{0}).$ (7)

Conversely, $x_{0}$ can be determined for a given $\bar{n}$. Thus, this
geometric construction, termed _flux-balance construction_ , shows that
significant features of the steady state profile are determined by the shape
of the nullcline.

Mass-redistribution instability. — In addition to the construction of
stationary patterns, it was shown in Ref. Brauns  _et al._ (2018) that the
nullcline shape determines the stability of a homogeneous steady state, and
that the mechanism underlying lateral (“Turing”) instability is a mass-
redistribution cascade. Specifically, it was found that a homogenous steady
state is laterally unstable when the _slope_ of the nullcline
$\chi(\bar{n}):=\partial_{m}c^{*}|_{\bar{n}}$ is steeper than the slope of the
FBS (see Section II.D1 in Ref. Brauns  _et al._ (2018) for a derivation),

$\chi(\bar{n})<-\frac{D_{m}}{D_{c}},$ (8)

which, using the mass-redistribution potential, Eq. (4), is equivalent to
$\partial_{n}\eta^{*}<0$. If this condition is fulfilled, high-density regions
act as cytosolic sinks, leading to further accumulation of mass and hence a
mass-redistribution cascade. This motivates the corresponding name _mass-
redistribution instability_.

Starting from a homogeneous steady state with a small random perturbation, the
initial dynamics is dominated by the fastest growing eigenfunction of the
linearized dynamics. At the onset of this instability there is a dominant
eigenfunction that determines the initial dynamics of the system. We can find
this dominant eigenfunction by linearizing the system around its homogenous
steady state (linear stability analysis). For a homogenous steady state, these
eigenfunction are Fourier modes and their growth rates are given by the
dispersion relation, shown as inset in Fig. 1(d). The fastest growing mode,
$q_{\text{c}}$, determines the length scale of the initially growing pattern
as illustrated in Fig. 1(d)). Subsequently the pattern coarsens into a single
peak Ishihara  _et al._ (2007); Otsuji  _et al._ (2007); Chiou  _et al._
(2018); Brauns  _et al._ (2018).

###  I.3 Pattern formation with a step-like template

A common feature of the biological examples we discussed in Section I.1 is
that the templates have a sharp edge, and that a downstream protein pattern
localizes to this edge. To obtain a conceptual understanding of how such an
edge-sensing mechanism might work, we study how an idealized step-like
template affects the pattern formation of the two-component McRD model as a
paradigmatic example.

We consider a step-like template profile $\theta(x)$ with a sharp edge at
$x_{\mathrm{E}}$

$\displaystyle\theta(x)$
$\displaystyle:=\begin{cases}\theta_{\mathrm{A}}&x\leq x_{\mathrm{E}}\\\
\theta_{\mathrm{B}}&x>x_{\mathrm{E}}\end{cases}\;,$ (9)

as illustrated in Fig. 2(a). Such a template defines two spatial subdomains
(labelled A and B). Here, we consider a template that couples to the
downstream pattern forming system via the local reactions $f(m,c;\theta(x)$,
such that different reactive dynamics

$f_{\mathrm{A,B}}(m,c):=f(m,c;\theta_{\mathrm{A,B}}),$ (10)

govern the system in the two subdomains (see Fig. 2(b)).

Figure 2:  Illustration of a step-like template (a) that acts on the reaction
kinetics defines two subdomains, labelled A and B, with respective reaction
kinetics $f_{\mathrm{A}}$ and $f_{\mathrm{B}}$. (b) We consider a two-
component system describing, for example, a single protein species that cycles
between a membrane bound and a cytosolic state. (c) The reactive flow due to
the reaction kinetics can be visualized in the $(m,c)$-phase space of
concentrations. Due to mass-conservation, the flow points along the reactive
phase spaces $m+c=n$ (indicated by gray lines). Along the nullclines
$f_{\mathrm{A,B}}$, the flow vanishes. Therefore, the reactive flow in each
subdomain is qualitatively captured by the shape of the respective nullcline.
(d) With the reactive flow encoded by the nullclines, the phase portraits of
the two subdomains can be combined (‘overlaid’) into a single phase space.
This combined phase space will be used for the construction of stationary
states.

Pattern forming systems with a step-like (also called “jump-type”)
heterogeneity have a rich history in the mathematical literature (see e.g.
Maini  _et al._ (1992); Ermentrout and Rinzel (1996); Bär  _et al._ (1996);
Yuan  _et al._ (2007); Nishiura  _et al._ (2007); Miyazaki and Kinoshita
(2007); Nishi  _et al._ (2013); Doelman  _et al._ (2016); van Heijster  _et
al._ (2019); Scheel and Weinburd (2018)) where, they have been studied in the
context of front-pinning Doelman  _et al._ (2016), pulse localization van
Heijster  _et al._ (2019), and wavenumber selection Scheel and Weinburd
(2018), to name a few recent examples. These studies predominantly focussed on
excitable media and the models studied are not mass-conserving. Furthermore,
the prevalent methods employed in these studies are singular perturbation
theory (Refs. Jones (1995); Ward (2006) may serve as general introductions)
and normal form theory (see e.g. Ref. Guckenheimer and Holmes (1983)). The
former method uses matched asymptotics and is based on a separation of spatial
scales; the latter applies in the vicinity of a bifurcation.

Here, we choose a conceptually different approach building on the recently
developed local equilibria theory Halatek and Frey (2018); Brauns  _et al._
(2018), specifically the phase-portrait analysis outlined in Sec. I.2 and
introduced in detail in Brauns  _et al._ (2018). Our starting point is to use
the reactive nullclines of the two subdomains as proxies for the respective
reactive flows (Fig. 2(c)). For this purpose, we combine the phase portraits
of the subdomains into a single phase portrait as shown in Fig. 2(d). This
‘overlaying’ of the phase portraits facilitates a geometric analysis of the
system with a step-like template based on the approach presented previously
for the two-component system on a homogeneous domain Brauns  _et al._ (2018)
(see recap in Sec. I.2). Throughout the paper we will use nullcline shapes as
illustrated in Fig. 2(d) (see Appendix A for the specific equations and
parameters). A different nullcline ‘arrangement’, and the general role of the
nullcline shapes are discussed in Appendix B.

The remainder of the paper is structured as follows. In Sec. II.1, we first
discuss how the spatial template leads to steady states with a monotonic
density profile, which we call _spatially heterogeneous base states_. We
analyze these base states in the phase space of reaction kinetics and show
that they are determined by a balance of reactive turnovers. Importantly, we
find that this balance can break down, such that there are regimes where no
base state (i.e. monotonic steady state) exists. In Sec. II.2 we show that, in
these regimes, _stationary patterns_ (i.e. non-monotonic steady states)
emerge, with a density peak either at the system boundary or at the template
edge. Specifically, we show that the peak at the template edge exist if the
nullclines intersect at a point where only one of them has negative slope, as
illustrated in Fig. 2(d). In Sec. II.3, we introduce the concept of _regional
instability_ to understand the transition between the monotonic base states to
the non-monotonic patterns and the conditions under which a peak emerges at
the template edge (‘edge sensing’). Importantly, we find simple geometric
criterion for edge sensing: The reactive nullclines of the two subdomains have
to intersect at a point where only one of them has negative slope. Finally, in
Sec. III, we ask how the peak responds to a moving template edge and show that
this can lead to depinning or suppression of the peak when the template edge
moves too fast.

##  II Construction of steady states and their bifurcations

The goal of this section is to characterize the steady states of the two-
component McRD system with a step-like template. These systems generically
don’t exhibit spatially homogenous states. Non-homogeneous steady states can
be categorized based on the monotonicity of their density profiles. In Sec.
II.1, we characterize the monotonic steady states, which consist of two
plateaus connected by a monotonic interface at the template edge
$x_{\mathrm{E}}$, as shown in Fig. 3(a). We call these monotonic steady states
the _(spatially) heterogenous base-states_. In subsection B we extend our
analysis to non-monotonic steady states which are the ‘genuine’ self-organized
patterns of the system.

Note that we restrict our construction to the nullcline shapes shown in Fig.
2(d) here. Generalizations to other arrangements follow the same principles
and can be worked out analogously.

###  II.1 Monotonic steady states (base states)

Figure 3:  Illustration of the construction of monotonic steady states
(heterogeneous base states) in $(m,c)$-phase space. (a) Density profile with
the inflection point at the template edge $x_{\mathrm{E}}$. (b) Phase space
including the reactive nullclines of subdomain A (orange) and B (blue) with
the corresponding density distribution (thick orange and blue line). (c) Total
turnover as a function of the membrane density at the template edge
$m_{\mathrm{E}}$. Total turnover balance determines the steady state value
$m^{*}_{\mathrm{E}}$. (d) Density profile for a case where the inflection
point of the profile, $x_{0}$, lies within subdomain A, that is, when
$m_{\mathrm{A}}^{0}>m_{\mathrm{B}}$ as illustrated in the corresponding phase
space representation (e). (f) In this case, the total turnover becomes a non-
monotonic function of $m_{\mathrm{E}}$, such that total turnover balance may
have multiple solutions, or no solution at all. These different cases are
sketched here for different average global total densities $\bar{n}$. (g)
Bifurcation diagram of the base state, showing the saddle-node bifurcation due
to breakdown of total-turnover balance. In the region beyond the saddle-node
bifurcation (shaded in gray) no monotonic base state exists. Note that
monotonicity enforces $m^{*}_{\mathrm{E}}-m_{\mathrm{B}}<0$ here.

In steady state, the net diffusive flux that redistributes mass must vanish
and reactive flows in $m$ and $c$ must be balanced. This balance is encoded in
the constraint that a stationary pattern’s phase space distribution must be
embedded in a flux-balance subspace (cf. Eq. (5)). This constraint is
independent of the local reactions, and, hence, also holds when the local
reaction are heterogeneous due to a template (purple dotted line in Fig.
3(b)).

Analogously to the construction of mesa patterns in Sec. I.2, we can
graphically construct the steady-state density profile in real space, as
illustrated in Fig. 3(b, e). To that end, we approximate the density profile
at the plateaus—where diffusive fluxes cancel everywhere—by the concentration
at the FBS-NC intersections, such that the concentration at the plateaus is
fully determined by the FBS-offset, $\eta_{0}$. Let us denote the membrane
concentration at the FBS-NC intersections as $m_{\mathrm{A}}^{-}(\eta_{0})$,
$m_{\mathrm{A}}^{0}(\eta_{0})$ and $m_{\mathrm{A}}^{+}(\eta_{0})$ for
subdomain A (orange nullcline) and as $m_{\mathrm{B}}(\eta_{0})$ for subdomain
B (blue nullcline). For specificity, we consider in the following a base state
that approaches the plateaus $m_{\mathrm{A}}^{-}$ and $m_{\mathrm{B}}$ far
away from the template edge. An analogous construction can be made for a
monotonic state that connects $m_{\mathrm{A}}^{+}$ and $m_{\mathrm{B}}$.333
Furthermore, note that we ignore potential intersection points of FBS and
B-nullcline at higher masses. In this regime, subdomain B will also exhibit
lateral instability. Here, we restrict ourselves to the regime where only
subdomain A becomes laterally unstable.

The conservation of average total density $\bar{n}$ enforces a constraint on
the construction of the base state. This constraint can be used to estimate
$\eta_{0}$ from the average total density $\bar{n}$. For a domain much larger
than the profile interface, the interface region can be neglected and the
average total density can be approximated by the weighted average of the
plateau densities in the two subdomains

$\begin{split}&\bar{n}\approx
x_{\mathrm{E}}\,n^{-}_{\mathrm{A}}(\eta_{0})+(L-x_{\mathrm{E}})\,n_{\mathrm{B}}(\eta_{0}),\end{split}$
(11)

with $n^{-}_{\mathrm{A}}$ and $n_{\mathrm{B}}$ the total density at the
plateaus in subsystem A and B respectively. This determines an implicit,
approximate relation between the control parameter $\bar{n}$ and the FBS-
offset $\eta_{0}$.

Upon changing $\bar{n}$, the plateau concentrations of the density profile
must change, and hence, $\eta_{0}$ must shift (cf. Eq. (11)). For the
laterally stable plateaus, the nullcline slope at the corresponding FBS-NC
intersections ($m_{\mathrm{A}}^{-}$ and $m_{\mathrm{B}}$) is larger than the
FBS-slope ($\partial_{n}\eta^{*}_{\mathrm{A,B}}(n)>0$, cf. NC-slope criterion
Eq. (8)). Hence, the relationship $\eta_{0}(\bar{n})$ must be monotonically
increasing for stable base states, as one sees by taking the derivative of Eq.
(11) w.r.t $\eta_{0}$, and using that monotonicity of a function implies
monotonicity of its inverse.

Note that, even though the base state looks similar to a mesa pattern in a
system on a homogeneous domain, the relationship between $\eta_{0}$ and
$\bar{n}$ makes a key difference between the two cases. As we discussed in the
introduction (Sec. IB), in a system on a homogeneous domain, changing
$\bar{n}$ does not affect $\eta_{0}$ in a system much larger than the
interface width, but instead shifts the pattern interface (cf. Eq. (7)). In
contrast, the interface position of a heterogeneous base state is determined
by the position of the template edge $x_{\mathrm{E}}$. Hence, to accommodate a
given average mass $\bar{n}$ the plateau concentrations
$n^{-}_{\mathrm{A}}(\eta_{0})$ and $n_{\mathrm{B}}(\eta_{0})$, determined by
FBS-NC intersection points, must adapt (cf. Eq. (11)). Thus, $\eta_{0}$ of a
heterogeneous base state depends directly on $\bar{n}$.

So far we have estimated the concentration at the two plateaus of the
monotonic steady state profile. To determine how these two plateaus are
connected at the template edge position $x_{\mathrm{E}}$, we use the condition
that in steady state the total reactive turnover in the system must vanish. In
the vicinity of the template edge at $x_{\mathrm{E}}$, the concentrations
deviate from the local equilibria, such that there are reactive flows
(illustrated by red arrows in Fig. 3(b)). Since the template introduces two
subdomains with different reaction kinetics, the total reactive turnover in a
system with a template is given by the sum over the turnover in the two
subdomains,

$\begin{split}F(m_{\mathrm{E}};\eta_{0})&=F_{\mathrm{A}}(m_{\mathrm{E}};\eta_{0})+F_{\mathrm{B}}(m_{\mathrm{E}};\eta_{0})\\\
&=\int_{m^{-}_{\mathrm{A}}}^{m_{\mathrm{E}}}{dm}\,f_{\mathrm{A}}(m,\eta_{0}-\frac{D_{m}}{D_{c}}m)\\\
&+\int_{m_{\mathrm{E}}}^{m_{\mathrm{B}}}{dm}\,f_{\mathrm{B}}(m,\eta_{0}-\frac{D_{m}}{D_{c}}m),\end{split}$
(12)

where $m_{\mathrm{E}}$ is the membrane concentration at the template edge. In
steady state, the total turnover $F(m_{\mathrm{E}};\eta_{0})$ must vanish such
that all reactive flows in the system balance. Thus, the solution of
$F(m^{*}_{\mathrm{E}};\eta_{0})=0$ (see Fig. 3(c)) determines the steady state
concentration at the template edge
$\widetilde{m}(x_{\mathrm{E}})=m^{*}_{\mathrm{E}}$. Note that, due to
monotonicity, $\eta_{0}$ and $m^{*}_{\mathrm{E}}$ uniquely identify a base
state of a given system.

For small enough $\bar{n}$, the second FBS-NC intersection
$m_{\mathrm{A}}^{0}$ for the A-nullcline is larger than the FBS-NC
intersection for the B-nullcline $m_{\mathrm{B}}$ as illustrated in Fig. 3(b).
In this case, both summands of Eq. (12) are monotonic in
$m_{\mathrm{E}}\in[m_{\mathrm{A}}^{-},m_{\mathrm{B}}]$ because the reactive
flow does not change sign within either subdomain, i.e. the inflection point
of the profile coincides with the template edge. Hence, there is only a single
solution $m^{*}_{\mathrm{E}}$ that fulfills total turnover balance. For larger
$\bar{n}$ (and thus $\eta_{0}$, cf. Eq. (11)), $m_{\mathrm{A}}^{0}$ can become
smaller than $m_{\mathrm{B}}$, as illustrated in the sketch in Fig. 3(d,e).
This entails that the position where the reactive flows change sign (i.e.
inflection point) lies in subdomain A (see Fig. 3(d,e)). Thus,
$F_{\mathrm{A}}(m_{\mathrm{E}};\eta_{0})$, and thereby also the total turnover
$F$ as a function of $m_{\mathrm{E}}$ becomes non-monotonic and may thus have
multiple roots. Indeed, for increasing $\bar{n}$, our flux-balance
construction predicts three different regimes: (i) A regime where there is one
solution in the interval $[m_{\mathrm{A}}^{-},m_{\mathrm{B}}]$, (ii) a regime
with two solutions, and (iii) a regime with no solution (as illustrated in the
sketch in Fig. 3(f)). In the last regime, total turnover balance becomes
impossible for a _monotonic_ steady state (base state). In Sec. II.2, we will
see how total turnover balance can be reached in this regime by a _non-
monotonic_ steady state.

The roots of $F$ correspond to different base states which we characterize by
the amplitude of the density profile in Subdomain B,
$m^{*}_{\mathrm{E}}-m_{\mathrm{B}}$; see Fig. 3(g). For monotonic states (i.e.
base states), $m^{*}_{\mathrm{E}}-m_{\mathrm{B}}$ is negative444Note that for
high-mass base states, monotonicity enforces
$m^{*}_{\mathrm{E}}-m_{\mathrm{B}}>0$, in the case of a nullcline arrangement
as shown in Fig. 2.. At the transition from regime (ii) to (iii), the base
state undergoes a saddle-node bifurcation at $\bar{n}_{\mathrm{SN}}$. From the
flux-balance construction and total turnover balance, we can estimate the
position of this bifurcation. At the saddle-node bifurcation point, the
minimum of $F$ coincides with the root of $F$. From Eq. (12) it follows that
$F$ reaches its minimum at
$f_{\mathrm{A}}(m_{\mathrm{min}},\eta)=f_{\mathrm{B}}(m_{\mathrm{min}},\eta)$.
Thus, this condition, together with $F(m_{\text{min}};\eta_{\text{SN}})=0$
implicitly determines the value of $\eta_{\text{SN}}$ at the saddle-node
bifurcation. From this we can then estimate $\bar{n}_{\text{SN}}$ via (11).

To test this approximate construction of steady states, we use specific
reaction terms $f_{\mathrm{A}}$ and $f_{\mathrm{B}}$ as specified in Appendix
A and compare the steady state profiles obtained from the flux-construction to
the profiles obtained from numerical continuation (see Appendix E for a short
description of numerical continuation and the comparison of steady states. A
more detailed explanation of continuation methods can be found in Ref.
Krauskopf  _et al._ (2007)). We find that the flux-balance construction gives
a estimate of the steady states profiles for sufficiently large system sizes
(see Appendix E).

As we noted above, there is also a family of base states which connects a
plateau at $m_{A}^{+}$ (instead of $m_{A}^{-}$) in subdomain A to
$m_{\mathrm{B}}$ in subdomain B. These base states have a high average mass,
and we will refer to them as ‘high-mass’ base states. Following the same
arguments as above, we find that these base states undergo a saddle-node
bifurcation when the average total mass is decreased below a critical average
mass, analogously to the saddle-node bifurcation of ‘low-mass’ base states
discussed in this section.

In summary, we have shown how to find monotonic steady states (base states)
with a flux-balance construction. Notably, we found that for a range of total
mass $\bar{n}$, this flux-balance construction has no solution, and hence,
there exist no monotonic steady states. In this regime, the steady states must
be non-monotonic. We next study these non-monontonic steady states, which we
refer to as _patterns_.

###  II.2 Non-monotonic steady states (patterns)

Figure 4:  Bifurcation structure and phase space construction of stationary
patterns. (a) Bifurcation structure of stationary states for the average mass,
$\bar{n}$, obtained by numerical continuation, together with spatial profiles
for the stable steady states (montonic base states and non-monotonic
patterns). In the area shaded in gray, no monotonic base states exist. Solid
(dashed) lines indicate stable (unstable) branches. The unstable branches in
the central region, marked by numbers correspond to unstable patterns with
multiple inflection points in subdomain A (see Fig. 13 in Appendix E). In
Supplementary Movies 1 and 2 we show finite element simulations in which we
adiabatically increase and decrease the mass through the low- and high-mass
saddle-node bifurcations. (b,c) Sketches of stable, non-monotonic, stationary
patterns together with the corresponding phase space constructions. At the
density profile’s extremum, marked by the dash dotted vertical line, there are
no gradients ($\partial_{x}\widetilde{m}=0=\partial_{x}\widetilde{c}$), such
that a notional no-flux boundary can be introduced. Thus, the resulting two
segments (labelled I and II), can now be studied separately. In the phase
portraits, the density distributions of the two segments connect at the
extremal concentrations $m_{\mathrm{A}}^{-}$ and $m_{\mathrm{A}}^{+}$, in (a)
and (b) respectively. They are shown slightly offset from the FBS (dashed
purple line) for visual clarity. The true density distribution must of course
be embedded in a single FBS to fulfill diffusive flux-balance.
Parameters for the bifurcation diagram: $D_{m}=0.01$, $D_{c}=0.5$,
$\hat{k}_{\mathrm{fb}}=0.25$, $\hat{k}_{\mathrm{off}}=4$,
$\theta_{\mathrm{B}}=20$, $\theta_{\mathrm{A}}=2$, $L=10$, $x_{\mathrm{E}}=5$.

To gain some intuition about the structure of the stationary patterns, we
first calculate them numerically as a function of the average total density
$\bar{n}$ using numerical continuation555See Appendix E for a brief
description of the core idea behind numerical continuation. An excellent
overview is provided in Ref. Krauskopf  _et al._ (2007). for a specific choice
of the reaction term $f(m,c;\theta)$ specified in Appendix A. The resulting
one-parameter bifurcation structure shows that the system exhibits two stable
patterns, one with a peak (high density region) at the template edge and one
with a peak at the system boundary, respectively (see Fig. 4(a)). In the
bifurcation structure, the two branches of stable patterns are connected by a
cascade of unstable patterns exhibiting multiple peaks; see Fig. 13 in
Appendix E. Their instability can heuristically be understood as a coarsening
process due to a competition for total density, similarly as in a system on a
homogeneous domain Brauns  _et al._ (2018). Some more technical aspects of
this bifurcation structure are discussed in Appendix E.

Can we use the flux-balance construction to construct non-monotonic steady
states as well? At the extrema of any stationary pattern in a two-component
McRD system, the gradients (and hence diffusive flux) in _both_ membrane and
cytosol concentration vanish simultaneously (cf. diffusive flux-balance Eq.
(5)). This allows us to place notional reflective boundaries at extrema,
effectively splitting the non-monotonic profile into monotonic segments. Thus,
we can use the flux-balance construction as described in Sec. II.1 to
construct the steady states in the two segments separately, with the
additional constraint of continuity at the boundaries connecting the segments.

The stable patterns in the two-component McRD system with a step-like template
have only a single peak, and hence only a single extremum within the domain
that splits the system into two segments (labelled I and II; see Fig. 4(b,c)).
Segment I is fully embedded in subdomain A, i.e. it is a system on a
homogeneous (sub)domain. Hence, for sufficiently large domain size, its steady
state is a mesa pattern666In the vicinity of the saddle-node bifurcations of
these patterns, segment I exhibits a peak pattern instead of a mesa pattern.
To obtain an approximation for this case, one would need to generalize the
peak approximation as discussed in Ref. Brauns  _et al._ (2018). as introduced
in Sec. I.2. The orientation of the mesa pattern in segment I determines
whether the density peak is located at the left domain boundary or at the
template edge. segment II contains the template edge, such that the steady
state in segment II is a heterogeneous base state.

By continuity, the FBS-offset $\eta_{0}$ must be identical in both segments.
Recall that for a mesa pattern, $\eta_{0}$ is determined by total turnover
balance, and independent of the average mass and domain size in the large
domain size limit (see Sec. I.2 and Ref. Brauns  _et al._ (2018)). We can thus
find $\eta_{0}$ solely by total-turnover balance in segment I, without
specifying the position of the boundary between the segments, and without
specifying the average masses in the two segments respectively. Instead, given
$\eta_{0}$, we find the average mass in segment II, $\bar{n}_{\mathrm{II}}$,
via Eq. (11), which depends on the choice for the orientation of the mesa
pattern in segment I. In segment II, subdomain B plays the role of a mass-
reservoir that absorbs a fraction of the total average mass and thus reduces
the mass available to the mesa pattern in segment I,
$\bar{n}_{\mathrm{I}}=\bar{n}-\bar{n}_{\mathrm{II}}$. Finally,
$\bar{n}_{\mathrm{I}}$ determines the position of the mesa pattern’s interface
in segment I via Eq. (7). Similarly, we can construct the (unstable) patterns
with multiple peaks by splitting the system into more than two segments.

We conclude that the flux-balance construction fully characterizes the
stationary patterns of the system with a step-like template. These steady
states exhibit density peaks similar to a system on a homogeneous domain,
however, the position of the density peak depends on the template edge
position.

We next ask which of the two stable patterns emerges as the base state ceases
to exist. To that end, we use finite element simulations and adiabatically
increase the total average density in a system such that it passes through the
base-state bifurcation (Supplementary Movie 1). The system evolves into a
pattern with a peak at the template edge (corresponding to the upper branch in
Fig. 4(a)). Upon further increase of $\bar{n}$ the peak widens and eventually
transitions into a mesa pattern at the template edge. The right hand interface
of the mesa patter remains localized at the template edge while the left hand
interface moves into subdomain A to accommodate the additional mass (cf. Eq.
7). Eventually for even larger $\bar{n}$, the mesa pattern ceases to exist as
its interface hits the left boundary of the domain. The system then
transitions to the ‘high-mass’ base state which connects the FBS-NC
intersection points $m_{\mathrm{A}}^{+}$ and $m_{\mathrm{B}}$. Going backwards
by adiabatically decreasing $\bar{n}$, the system passes through the ‘high-
mass’ base state’s saddle-node bifurcation. The corresponding regional
instability leads to the formation of a trough, rather than a peak, at the
template edge. The resulting stationary pattern has a minimum at the template
edge and a maximum at the left boundary (corresponding to the lower branch in
Fig. 4(a)) (Supplementary Movie 2). Upon further decrease of $\bar{n}$, the
interface of this pattern will reach the left boundary of the domain such that
the system transitions back into the low-mass base state.

These transitions also take place when the average mass is changed non-
adiabatically but still so slow that mass-transport across the system by
cytosolic diffusion (${\sim}L^{2}/D_{c}$) is faster than the rate at which
mass is added or removed. Interestingly, when we increase the mass on a non-
adiabatic timescale we observe multiple transient patterns, which we
characterize in Appendix D. Furthermore, we show that we can get similar
transitions between the base state and the patterns if, instead of increasing
the average mass, the local reactions $f_{\text{A,B}}$ in the two subdomains
are varied over time (Appendix C).

Taken together, we have shown that the flux-balance construction can be used
to construct non-monotonic steady states by splitting the density profile into
monotonic segments. This is possible because, the stationary pattern profile
can be split at extrema, where all diffusive fluxes vanish. We found that the
system can exhibit two patterns, with a density peak either at the system
boundary or at the template edge. The peak at the template edge only exists
when the two nullclines intersect at a point where only one of them has
negative slope. Furthermore, our finite element simulations show that
increasing the mass, starting from the low-mass base state, leads to a peak at
the template edge. Vice versa, decreasing mass, starting from the high-mass
base state leads to a peak at the system boundary. In the next Section we
provide a heuristic argument to understand under which conditions the peak
emerges at the template edge.

###  II.3 Template-induced regional instability

We next want to understand the mechanism of pattern selection as the system
goes through the saddle-node bifurcation(s) where the base state ceases to
exist. We first show that a (numerical) linear stability analysis explains why
either a peak or a trough pattern grows at the template edge, as the system
goes through the bifurcation. We then provide a heuristic argument to explain
this edge-sensing mechanism, and formulate a geometric criterion under which
this mechanism works, based on the shape of the nullclines.

To study how the base state develops into a pattern, as the system goes
through the saddle-node bifurcation, we consider a base state
$(\widetilde{m}(x),\widetilde{c}(x))$ in the vicinity of the bifurcation point
and analyze the dynamics of a small perturbation $(\delta m(x,t),\delta
c(x,t))$. The dynamics of the perturbed state, up to linear order is given by:

$\displaystyle\partial_{t}\delta m(x,t)$
$\displaystyle=D_{m}\partial_{x}^{2}\delta m+\tilde{f}_{m}(x)\delta
m+\tilde{f}_{c}(x)\delta c\,,$ (13a) $\displaystyle\partial_{t}\delta c(x,t)$
$\displaystyle=\kern 3.00003ptD_{c}\partial_{x}^{2}\delta c\kern
4.49997pt-\tilde{f}_{m}(x)\delta m-\tilde{f}_{c}(x)\delta c\,.$ (13b)

The linearized reaction coefficients

$\tilde{f}_{m,c}(x)=\partial_{m,c}f\big{|}_{(\widetilde{m}(x),\widetilde{c}(x))}$
(14)

are not homogeneous in space and hence Eq. (13) is a set of linear PDEs with
nonconstant coefficients. We seek solutions of the form $\delta
m(x,t)=\Phi_{m}(x)e^{\sigma t}$, $\delta c(x,t)=\Phi_{c}(x)e^{\sigma t}$. With
this ansatz, Eq. (13) turns into the Sturm–Liouville eigenvalue problem

$\displaystyle\sigma\Phi_{m}(x)$
$\displaystyle=D_{m}\partial_{x}^{2}\Phi_{m}+\tilde{f}_{m}(x)\Phi_{m}+\tilde{f}_{c}(x)\Phi_{c}\,,$
(15a) $\displaystyle\sigma\Phi_{c}(x)$ $\displaystyle=\kern
3.30002ptD_{c}\partial_{x}^{2}\Phi_{c}\kern
3.99994pt-\tilde{f}_{m}(x)\Phi_{m}-\tilde{f}_{c}(x)\Phi_{c}\,,$ (15b)

for the eigenvalues $\sigma$ and the associated eigenfunctions
$(\Phi_{m},\Phi_{c})(x)$.

The defining feature of a a saddle-node bifurcation is that one eigenvalue
vanishes exactly at the bifurcation point. The associated eigenfunction
reveals the flow structure of the dynamics on the slowest timescale close to
the bifurcation point (center manifold, see e.g. Ref. Guckenheimer and Holmes
(1983)). From this we can gain intuition about the fate of the system upon
passing through the bifurcation.

Figure 5:  Regional lateral instability at the base state’s saddle-node
bifurcation. (a) The numerically calculated eigenfunction, $\Phi_{m}$,
associated to the vanishing eigenvalue at the saddle-node bifurcation is
localized at the template edge. (Parameters as in Fig. 4 at the saddle-node
bifurcation $\bar{n}\approx 2.65541$). (b) The concentration profile
$\widetilde{m}(x)$ of the base state at the saddle-node bifurcation. A spatial
region that is fully contained in subdomain A and centered around the
profile’s inflection point $x_{0}$ is marked in green. (c) This spatial region
corresponds to a phase space region where the dynamics is guided by a section
of the nullcline with a negative slope, i.e. where the system is laterally
unstable. (The phase space is shown as a sketch for visual clarity. See Fig.
8(a) in Appendix B for a plot from numerical simulation.)

At the base-state saddle-node bifurcation (marked SN in Fig. 4(a)), the
numerically calculated777The Sturm–Liouville eigenvalue problem at the
numerically calculated base-state saddle-node bifurcation is solved by
discretizing the spatial derivatives (Laplace operator) and solving the
resulting eigenvalue problem numerically in Mathematica. For further details
see e.g. Ref. Pryce (1993). eigenfunction is peaked in the vicinity of the
template edge (see Fig. 5(a)). This localized eigenfunction indicates that the
density profile will change most in the vicinity of the template edge, giving
rise to either a peak or a trough as the system goes through the saddle-node
bifurcation.

An intuition why the neutral eigenfunction at the base-state saddle-node
bifurcation is peaked at the template edge can be gained from the phase
portrait as sketched in Fig. 5(c)). Recall, that the inflection point of the
base state’s density profile $m_{\mathrm{A}}^{0}=\widetilde{m}(x_{0})$ lies
within subdomain A (cf. Fig. 3(d,e)). Consider a region centered around
$x_{0}$, fully contained within subdomain A, as indicated in green in Fig.
5(b). In phase space this point lies on a section of the A-nullcline with a
slope steeper than the FBS. Suppose for a moment that this region is isolated
from the rest of the system. Then, as the slope of the nullcline at
$m_{\mathrm{A}}^{0}$ is steeper than the slope of the FBS, the homogeneous
equilibrium in this region will be unstable due to a mass-redistribution
instability. This instability will set in when the region is large enough to
contain the shortest growing mode888More precisely, the size of the region
centered around the inflection point, $2(x_{\mathrm{E}}-x_{0}$), must be
larger than the wavelength of the shortest growing mode
$\pi/q_{\text{min}}(n_{0})$.. We call this a regional (mass-redistribution)
instability.

A necessary condition to trigger a regional instability at the template edge
is that the nullclines of the two subdomains cross at a point where the the
A-nullcline fulfills the nullcline-slope criterion for lateral instability,
Eq. (8) (see Fig. 6(a)). When this edge-sensing criterion is not fulfilled, as
shown in Fig. 6(b), the regional instability sets in at the system boundary
first, giving rise to a peak at the system boundary and not at the template
edge (as illustrated in Fig. 8(b) in Appendix B). Because the shapes of the
nullclines in the two subdomains relative to each other depend on how the
template affects the reaction kinetics, the edge-sensing criterion constrains
models that can exhibit edge sensing. In Appendix F we demonstrate that the
edge-sensing criterion precisely predicts the regime of edge sensing for a
phenomenological model of Cdc42. We furthermore show that edge sensing only
works if the template increases both the attachment and detachment rate of
Cdc42 in one of the subdomains. This may help to identify the relevant
molecular players in biological systems.

Figure 6:  Nullcline criterion for edge sensing, i.e. the emergence of a
stable density peak at the template edge. Edge sensing is possible is the
nullclines intersect in a point where one of them has negative slope (a). If
they don’t intersect (b) or intersect in a point where they have both positive
(or both negative) slope, stable peaks only exist at a system boundary, and
hence, edge sensing is not possible.

The concept of regional instability was already discussed in Ref. Brauns  _et
al._ (2018) in the context of excitability (“nucleation and growth”) for a
homogeneous domain. There, a stable homogeneous steady state is perturbed by
moving mass from the system into small region. When this region contains
sufficient mass, it can become laterally unstable and thus formation of a peak
pattern is “nucleated.” With that, the regional instability at the template
edge can also be understood in terms of lateral excitability. From the
perspective of subdomain A, the differing reaction kinetics in subdomain B
induces a perturbation at the subdomain interface, that is, the template edge
(orange line in Fig. 5(b)). In that sense, the base state is a perturbed
homogeneous steady state in each subdomain. If in subdomain A, this
perturbation becomes large enough to cross the nullcline in a section of
negative slope (Fig. 5(c)), it triggers a lateral instability and thus the
formation of density peak at the template edge. This relationship to
excitability highlights that it is the spatial _gradient_ of the reaction
kinetics due to the template that localizes the instability. Heuristically
this can be pictured as sensing the spatial derivative of the template. Here
we focussed on a sharp template edge. For future work, it will be interesting
to study also a smooth template edge. Intuitively, if the template gradient is
too shallow, it will not induce a laterally unstable region, because the
induced deviation from the local equilibria that effectively acts as the
perturbation in the analogy to excitability will be too small.

In conclusion, we showed in this section that the template localizes the
patterns and determines the position of the instability from which they
emerge. Importantly, this instability determines which of the two stable
stationary patterns forms when the system passes through a bifurcation of the
base state. Finally, we presented a simple geometric criterion, shown in Fig.
6, that determines when the regional instability is localized at the template
edge.

##  III Moving template edge

Until now we considered pattern formation with a fixed template edge position
$x_{\mathrm{E}}$. We next ask what happens when $x_{\mathrm{E}}$ moves after a
peak at the template edge has been established (Fig. 7).

When the template edge moves, the peak must adapt to the new template edge
position. In order to reach the new stationary state, mass must be transported
from one side of the peak to the other. Thus, we expect that the peak follows
the template edge position as long as the mass is transported faster than the
velocity $v_{\mathrm{E}}$ at which the template edge moves, i.e. for
$v_{\mathrm{E}}\ll D_{c}/w$, where $w$ is the width of the peak. To test the
intuition for this adiabatic case, and study what happens in the non-adiabatic
case of a fast-moving template edge, we turn to numerical simulations. To
probe a range of template velocities $v_{\mathrm{E}}$, we quadratically
increase the template edge velocity during the simulation. At a distance
$0.3\,L$ from the system boundary, the template movement is stopped.
Furthermore, we move the interface either to the right, away from the peak
Fig. 7(a,c) or to the left, towards the peak Fig. 7(b,d). In the adiabatic
case, we find—in agreement with our expectation—that the peak remains pinned
to the template edge position (see Fig. 7(a,b) and SI Movies 3 and 4).

Figure 7:  Time evolution of $m$-profile due to dynamic template interface
position. The supplementary material contains a movie for each of the four
scenarios. (a) Slow pull (the template edge moves away from the peak): This
results in pinning of the peak to the template edge (see SI Movie 3). (b) Slow
push (the template edge moves towards the peak): This results in pinning of
the peak to the template edge (see SI Movie 4). (c) Fast pull: This results in
depinning as the peak cannot follow the template edge. However, the peak stays
where it depins as it is quasi-stable in region A on the observed timescale
(see SI Movie 5). (d) Fast push: This results in suppression as the peak is
not stable in region B. As the peak dissipates the average mass at the new
edge position slowly increases and potentially (if there is enough total mass
in the system) leads to to a re-entrance of the peak (see SI Movie 6).
Parameters: $D_{m}=0.01$, $D_{c}=10$, $k_{\mathrm{on}}=1$,
$\hat{k}_{\mathrm{fb}}=1$, $k_{\mathrm{off}}=2$, $\hat{K}_{\mathrm{d}}=1$,
$\bar{n}=5$, $\theta_{\mathrm{B}}=20$, $\theta_{\mathrm{A}}=0.5$, $L=20$,
$x_{\mathrm{E}}=3/5\,L$, $v_{\mathrm{E}}(t)=\pm 2.4\times 10^{-14}\,t^{2}$ in
(a) and (b), $v_{\mathrm{E}}(t)=\pm 8.9\times 10^{-13}\,t^{2}$ in (c) and (d).

In the non-adiabatic case, when the template edge moves faster than the peak
can follow by diffusive mass transfer, the peak position will shift relative
to the template edge. We find that, when the template moves away from the peak
(‘pulling’), the peak depins at a critical velocity and stops moving (see Fig.
7(c) and SI Movie 5). In contrast, when the template moves towards the peak
(‘pushing’), the peak is suppressed at some critical edge-velocity (see Fig.
7(d) and SI Movie 6). Interestingly, the critical velocity for depinning while
pulling is much lower than the critical velocity for suppression while
pushing.

To heuristically understand this difference of critical velocities, we take a
closer look at how the peak adapts to a shifted template position. Suppose for
a moment that we keep the peak profile frozen while shifting the template edge
$x_{\mathrm{E}}$ by a small amount $\delta x_{\mathrm{E}}$. This will change
the local reactions in the vicinity of the template edge. When the peak
profile is now “released”, these reactions will lead to changes in the
concentrations $m$ and $c$, and hence in the mass-redistribution potential
$\eta(x,t)$ between the original and shifted edge positions $x_{\mathrm{E}}$
and $x_{\mathrm{E}}+\delta x_{\mathrm{E}}$. The resulting $\eta$ gradient
leads to mass transport (Eqs. (3)), which in turn causes a movement of the
peak.999A geometric analysis in phase space, similar to the one presented in
Figs. 3 and 4, shows that the $\eta$ gradient is always such that the peak
moves in the direction that the template edge was shifted, until it reaches
it’s new stationary (pinned) position at the shifted edge position. The
difference between pulling and pushing, i.e. moving the template edge away
from the peak or towards it, is the amplitude of the $\eta$ gradient that
builds up as the template edge is shifted. In the case of pulling, the
template edge moves into the flat tail of the peak, such that the $\eta$
gradient decreases with increasing distance between peak and template edge. In
contrast, while pushing, the template edge moves into the steep interface of
the peak, leading to an continually increasing $\eta$ gradient. Only when the
template edge has shifted beyond the maximum of the peak, the induced $\eta$
gradient starts decreasing again.

This effect qualitatively explains the different critical velocities for
depinning (while pulling) and suppression (while pushing). In the former case
the induced gradients in $\eta$, and hence, the rate of mass transport that
shift the peak towards the moving template edge are small and decreasing with
peak-to-edge distance. Therefore, depinning is self-enhancing. In the case of
pushing, the $\eta$ gradient keeps increasing as the peak-to-edge distance
decreases. This in turn, increases the speed of the peak due to faster mass
transport. Only when the template edge has crossed the peak maximum, the peak
will be suppressed because it then lies mostly within the laterally stable
subdomain $B$ that does not support a stable peak.

In summary, we showed that in the case of an adiabatically slow template
motion, the peak stays pinned at template edge. In the non-adiabatic case we
found a qualitatively and quantitatively different behavior between pushing
and pulling. Pulling leads to depinning, while pushing eventually leads to
suppression. Furthermore, we heuristically explained why critical velocity for
pulling (depinning) is lower than the critical velocity for pushing
(suppresion). Going forward, it would be interesting to study this behavior
more systematically, both in numerical simulations and on an analytic level.
For example, concepts like response functions Ziepke  _et al._ (2016),
projection onto slow manifolds Nishiura  _et al._ (2007) and singular
perturbation theory Doelman  _et al._ (2016); van Heijster  _et al._ (2019)
may help to estimate the critical velocities for depinning and suppression.

##  IV Discussion

We showed how protein pattern formation can be controlled by a spatial
template (e.g. an upstream protein pattern), which acts on the proteins
interaction kinetics. In particular, we demonstrated for two-component McRD
systems how a step-like template—which defines two subdomains with different
reaction kinetics—can localize the formation of a peak pattern to the template
edge. We explained this edge-sensing mechanism by a regional (mass-
redistribution) instability that emerges at the template edge position. This
is in contrast to pattern formation on a homogeneous domain, where the
instability is generically “delocalized” (Fourier modes, cf. Fig. 1(d)), such
that noisy initial conditions have a strong impact on pattern formation
process.

Our analysis is based on a recently developed theoretical framework Halatek
and Frey (2018); Brauns  _et al._ (2018), termed local equilibria theory. This
theory proposes to analyze McRD systems as dissected into diffusively coupled
compartments, so small that each of the compartments can be considered as
well-mixed. For the paradigmatic (minimal) case of two-component systems, this
framework enables one to to perform phase-portrait analysis of the interplay
between local reactions and diffusive transport in the phase space of the
reaction kinetics. Here, we have extended this phase-portrait analysis to
incorporate the effect of a step-like heterogeneity of the reaction kinetics
in the spatial domain. We were able to construct the bifurcation diagram for
the average mass $\bar{n}$, which is a natural control parameter as it can be
controlled by production or degradation/sequestration of proteins in cells
(e.g. in a cell cycle dependent manner). We found that, at a critical average
mass, the system’s base state undergoes a saddle-node bifurcation, such that
the system transitions to a stationary pattern, with either a peak at the
template edge or at the system boundary. The peak forms at the template edge
if the template triggers a regional (mass-redistribution) instability at the
template edge. The phase-portrait analysis enables us to formulate a geometric
criterion for this edge-sensing mechanism. In particular, we show that the
step-like template can trigger a regional instability at the template edge, if
the nullclines of the reaction kinetics in the two subdomains intersect at a
point where one of them has a negative slope (more precisely, a slope steeper
than the negative ratio of the diffusion constants, $-D_{m}/D_{c}$), as
illustrated in Fig. 6. Finally, we showed that the edge-localized peak is
stable when the template edge is slowly moved and demonstrated that
qualitatively different processes—depinning vs. suppression—lead to the loss
of the edge-localized peak when the template is shifted too rapidly away from
the peak (‘pulling’) or towards it (‘pushing’).

We speculate that the edge-sensing mechanism studied here might underly
formation of the actomysin ring during macropinocytosis and cellular wound
healing, as well as the direction of Cdc42 polarization in budding yeast
adjacent to the previous bud-site. In macropinocytosis a high density PIP3
domain has been suggested to act as a template for a ring of actin nucleators
that localize to its periphery Veltman  _et al._ (2016). Similarly, during
cellular wound healing, the inside Abr could act as a template for the outside
Cdc42 ring which then drives recruitment of actin via formins Vaughan  _et
al._ (2011). Finally, in budding yeast, landmark proteins that localize to the
previous bud-site can be pictured as a template that suppresses Cdc42
accumulation at the previous bud-site and simultaneously localizes Cdc42
cluster to its vicinity Chant and Pringle (1995); Tong  _et al._ (2007); Lo
_et al._ (2013); Miller  _et al._ (2017). Furthermore, the spatio-temporal
organization of intracellular membranes, like the Golgi apparatus, endosomes,
and the endoplasmatic reticulum, involves cascades of coupled GTPases Pfeffer
(2012); Novick (2016); Noack and Jaillais (2017). Hence, we speculate that
this organization may rely on similar domain-edge sensing mechanisms.

The edge-sensing criterion (cf. Fig. 6) based on the shape of the nullcline,
may provide guidance for the mathematical modeling of these systems and
thereby help to identify the key bio-molecular players and processes. The
nullcline shapes of a given model constrain the ability of this model for
edge-sensing. As an example, we showed for an phenomenological two-component
model for Cdc42 pattern formation Mori  _et al._ (2008) that edge sensing
requires a template which increases both the attachment and detachment rate of
Cdc42 in one subdomain. Indeed in single-cell wound healing, the protein Abr
could provide such a template for Cdc42 since it is both a guanine exchange
factor (GEF) and a GTPase-activation protein (GAP) for Cdc42 Vaughan  _et al._
(2011).

Beyond the understanding of living systems, our results may also advance the
field of synthetic biology. Previous studies have explored mechanisms by which
a gradient can position a sharp front pattern via a bistable reaction-
diffusion system Rulands  _et al._ (2013); Zadorin  _et al._ (2017). The edge-
sensing mechanism, presented in this paper, is a candidate for a further
building block to design spatial protein patterns.

In future studies, it would be interesting to generalize our results beyond
the paradigmatic case of a single, stationary, step-like template in one
spatial dimension. Templates with multiple steps may be dissected into
segments with single steps that can then be studied separately. Furthermore,
it has been shown that the geometry of a cell indirectly affect the
attachment–detachment kinetics via the ratio of bulk-volume to surface-area
Thalmeier  _et al._ (2016); Ziepke  _et al._ (2019), and curvature sensing
proteins Haupt and Minc (2018). Another promising direction is to incorporate
the dynamics of the template itself as a self-organized pattern-forming
system, and include a feedback from the downstream pattern to the template.
Such feedback may give rise to complex spatio-temporal behavior like
oscillatory patterns and traveling waves that can then be characterized by
building on the phase-portrait analysis presented here. Our results on the
moving template (Sec. III) and non-adiabatic upregulation of average mass
(Appendix D) indicate that the edge sensing works beyond the adiabatic regime.

Finally, even for the elementary case studied here, many important questions
remain open. We showed that a moving template will lead to a loss of the edge-
localized peak due to depinning (while pulling) or suppression (while pushing)
at different critical edge-velocities. In future work, these transitions
should be studied in more detail both numerically and analytically, e.g. using
a response function formalism Ziepke  _et al._ (2016). Furthermore, an
analytic approach employing asymptotic methods like singular perturbation
theory Ward (2006); van Heijster  _et al._ (2019); Doelman  _et al._ (2016);
Nishiura  _et al._ (2007) may help to cast our heuristic explanation of the
localized eigenfunction, based on the concept of regional instability, into a
more rigorous argument. Similarly, such an approach may elucidate the
transition from edge-localized peaks to boundary-localized peaks for too fast
mass upregulation (cf. Appendix D, Fig. 11). In general, we expect that
combining mathematical tools like singular perturbation theory with the local
equilibria framework will be a fruitful approach to systematically study
complex pattern-forming systems.

###### Acknowledgements.

We thank Jacob Halatek and Korbinian Pöppel for discussions and feedback on
the manuscript. This research was supported by the Deutsche
Forschungsgemeinschaft (DFG) via project A6 within the Collaborative Research
Center “Spatiotemporal dynamics of bacterial cells” (TRR 174). We also
gratefully acknowledge financial support by the DFG Research Training Group
GRK2062 (Molecular Principles of Synthetic Biology). M.C.W. is supported by a
DFG fellowship within the Graduate School of Quantitative Biosciences Munich.

##  Appendix A Reaction kinetics and template definition

Throughout this paper we use a two-component McRD model on a one-dimensional
domain with one protein species (cf. Eq. (1a) and (1b)). The proteins cycle
between a cytosolic state (concentration $c(x,t)$) and membrane bound state
(concentration $m(x,t)$) as specified by the reaction term $f(m,c)$.
Importantly, our results are based on the shape of the reactive nullclines
and, hence, don’t depend on the specific choice for $f(m,c)$. To illustrate
our findings, we use biochemically motivated reaction kinetics, comprising
attachment, $a(m)$, and detachment, $d(m)$, reactions

$\displaystyle f(m,c):=a(m)c-d(m)m\ .$ (16)

Specifically we use

$\displaystyle a(m)$
$\displaystyle:=\left(k_{\mathrm{on}}+k_{\mathrm{fb}}\,m\right),$ (17a)
$\displaystyle d(m)$
$\displaystyle:=\frac{k_{\mathrm{off}}}{K_{\mathrm{d}}+\,m},$ (17b)

as introduced before in Ref. Brauns  _et al._ (2018). These reaction kinetics
describe a protein species that can attach from the cytosol to the membrane
with a rate $k_{\mathrm{on}}$ and get recruited from the cytosol to the
membrane by membrane bound proteins with a rate $k_{\mathrm{fb}}$. Membrane
bound proteins detach from the membrane via an enzymatic process described by
first order Michaelis-Menten kinetics, parameterized by the rate
$k_{\mathrm{off}}$ and the dissociation constant $K_{\mathrm{d}}$.

We consider systems were the reaction rates are different in subdomains A and
B. This externally imposed heterogeneity was introduced as a step-like
template in Sec. I.3 (cf. Eq. (9)). For specificity, we choose a template that
affects the reaction rates, such that the reactive nullcline in subdomain B is
stretched along the $m$-axis with respect to the reactive nullcline in
subdomain A. To that end, we rescale the feedback rate and the dissociation
constant scale with the template value, such that these rates become space
dependent

$\displaystyle k_{\mathrm{fb}}(x)$
$\displaystyle:=\hat{k}_{\mathrm{fb}}/\theta(x),$ (18a) $\displaystyle
K_{\mathrm{d}}(x)$ $\displaystyle:=\hat{K}_{\mathrm{d}}\,\theta(x).$ (18b)

The reaction term then becomes

$\displaystyle
f(m,c;\theta)=\left(k_{\mathrm{on}}+\hat{k}_{\mathrm{fb}}\,\frac{m}{\theta(x)}\right)\,c-\frac{k_{\mathrm{off}}\,\frac{m}{\theta(x)}}{\hat{K}_{\mathrm{d}}+\frac{m}{\theta(x)}}\;.$
(19)

For convenience, we do not specify units of length and time. In an
intracellular context a typical size would be $L\sim 10\;\mu\text{m}$, and
typical diffusion constants are $D_{m}\sim
0.01{-}0.1\;\mu\text{m}^{2}\text{s}^{-1}$ on the membrane and $D_{c}\sim
10\;\mu\text{m}^{2}\text{s}^{-1}$ in the cytosol. Rescaling to different
spatial dimensions is straightforward. To fix a timescale, the kinetic rates
can be rescaled with respect to the attachment rate $k_{\mathrm{on}}$. In an
intracellular context, typical attachment timescales are on the order of
seconds, i.e. $k_{\mathrm{on}}\sim\mathrm{s}^{-1}$.

##  Appendix B Nullclines without edge-sensing

In Sec. II.3 in the main text, we analyzed the edge-sensing mechanism based on
the reactive nullclines in phase space. In this analysis, the effect of the
heterogeneous reaction kinetics, i.e. the template, is captured by the shapes
of the reactive nullclines. From our analysis, we found a criterion for the
edge-sensing mechanism, as illustrated in Fig. 5. The nullclines need to
intersect at a point where only one nullcline is steeper than FBS (cf. Eq.
(8)).

In this appendix, we discuss a case where the criterion for edge sensing is
not fulfilled, such that a localization of the regional instability to the
edge is not possible, and a peak at the template edge does not exist (cf. Fig.
6(b)). As an example, we consider a template, affecting the reaction kinetics
such that the nullcline is stretched along the $c$-axes instead of the
$m$-axis, as shown in Fig. 8(b).

For the specific attachment–detachment reaction kinetics (Eq. (17)), the
nullclines is stretched along the $c$-axis via the off-rate while keeping all
other rates constant:

$k_{\mathrm{off}}(x):=\hat{k}_{\mathrm{off}}\,\theta(x).$ (20)

The resulting reaction term then reads

$\displaystyle
f(m,c;\theta)=\left(k_{\mathrm{on}}+k_{\mathrm{fb}}\,m\right)\,c-\theta(x)\frac{\hat{k}_{\mathrm{off}}\,m}{K_{\mathrm{d}}+m}.$
(21)

Following the same arguments as presented in Sec. II.1,we can construct the
base states (monotonic steady state). Starting from the low-mass base state,
illustrated in Fig. 8(a), and increasing the average mass $\bar{n}$ results in
an upwards shift of the FBS by the same argument as presented in Sec. II.1.
When the mass is further increased the FBS moves to the level where the two
FBS-NC intersection points on the A-nullcline, $m_{A}^{-}$ and $m_{A}^{0}$,
annihilate in a saddle-node bifurcation, as shown in Fig. 8(b). If $\eta_{0}$
increases beyond that point, the base state vanishes. Note that the origin of
the saddle-node bifurcation lies in the annihilation of the two FBS-NC
intersection points. This is different from the saddle-node bifurcation that
occurs for the nullclines we analyzed in Sec. II.1 in the main text. There,
the base state vanishes due to a breakdown of total turnover balance which
becomes apparent by the ‘annihilation’ of the two solutions for
$m_{\mathrm{E}}$ of Eq. (12); cf. Fig. 3.

To study the dynamics in the vicinity of the base state bifurcation, we use
the concept of regional instability (cf. Sec. II.3). The part of the density
distribution that enters the laterally unstable region in phase space
corresponds to the concentration at the left hand system boundary ($x=0$)
(Fig. 8(b)). Hence, upon crossing the base state’s saddle-node bifurcation, a
regional instability emerges at this system boundary, and a peak forms there.

Moreover, for nullclines shown in Fig. 8, there is only one way to construct a
stationary pattern with a single interface (inflection point) within subdomain
A (in addition to the interface imposed by the template step). This pattern
always has a density peak at the system boundary ($x=0$) and decreases
monotonically in space towards $x=L$. A stable peak localized to the template
edge does not exist in this case.

Figure 8:  Base states for nullclines shapes that do not facilitate edge-
sensing. (a) The low-mass base state with a high-concentration plateau in
subdomain A and a low-concentration plateau in subdomain B. (b) At a critical
average mass, the base state undergoes a saddle-node bifurcation. The
bifurcation arises from the ‘annihilation’ of the FBS-NC intersection point
$m^{-}_{A}$ and $m^{0}_{A}$, and not from the break down of turnover balance
as discussed in Sec. II.3 the main text. In this case, a regional instability
is triggered at the system boundary. An adiabatic sweep of $\bar{n}$ through
the saddle-node bifurcation of the FBS-NC intersection points is shown in
Supplementary Movie 7.

##  Appendix C Temporal variation of the template

Figure 9:  An adiabatically changing template triggers peak formation at the
template edge. See also Supplementary Movie 8.
(a) Initial state with a homogeneous template profile
$\theta_{\mathrm{A}}=\theta_{\mathrm{B}}$. (b) Base state analogous to Fig. 3.
(c) Base state right before the saddle-node bifurcation. The region that
becomes laterally unstable in the bifurcation is highlighted in light green in
the spatial profile and in phase space (cf. Fig. 5). (d) Peak pattern state
after the system transitioned through the bifurcation. This state is
qualitatively the same as sketched in Fig. 4(c).
Parameters as in Fig. 7 but with $L=5$, $\theta_{\mathrm{B}}=20$, and
$\theta_{\mathrm{A}}^{\mathrm{{f}}}=0.5$.

In Sec. II in the main text, we considered upregulation of average mass
$\bar{n}$ while the spatial template was kept constant. One can perform a
similar analysis for a varying template, while keeping the average mass
constant. In this scenario, not $\bar{n}$ but the parametrization of the
template (e.g. $\theta_{\mathrm{A}}$) serves as a bifurcation parameter. In
the biological context, this corresponds to a dynamically varying upstream
protein pattern.

To demonstrate that varying the template induces equivalent bifurcations as
variation of the average mass, we use the template as in the main text (cf.
Eq. (19)) and perform a numeric simulation where we let the template
$\theta(x,t)$ slowly vary with time (for simplicity only in subdomain A):

$\displaystyle\theta(x,t)$
$\displaystyle=\begin{cases}\theta_{\mathrm{A}}(t)&x\leq x_{\mathrm{E}},\\\
\theta_{\mathrm{B}}&x>x_{\mathrm{E}}.\end{cases}\;$ (22)

We initialize the template value in subdomain A at $\theta_{A}(0)=\theta_{B}$
and let it change via a sinusoidal ramp to its final value
$\theta_{A}(t_{\mathrm{f}})=\theta_{A}^{\mathrm{f}}$.

An exemplary simulation is shown in Fig. 9 and Supplementary Movie 8. At the
start of the simulation, the reaction rates, and thus the reactive nullclines,
are the same in subdomain A and subdomain B, which is illustrated by the
overlapping nullclines in Fig. 9(a). The template is homogeneous, and the
corresponding steady state is a homogeneous density profile for sufficiently
low average mass. Upon decreasing $\theta_{\mathrm{A}}$, the reaction rates in
subdomain A change, leading to a change in the shape of the reactive nullcline
(Fig. 9(b)). The resulting base state is equivalent to the low-mass base
state, similar to the case for mass upregulation analyzed in Sec. II.1. When
$\theta_{\mathrm{A}}$ is further decreased, the density profile becomes
regionally unstable at the template edge(Fig. 9(c)), which triggers a peak at
the template edge (Fig. 9(d)). This shows that the pattern formation process
as discussed in the main text can also be realized with a dynamic template.

##  Appendix D Non-adiabatic mass upregulation

In Sec. II in the main text, we analyzed the steady states as a function of
average mass and found that the system undergoes a transition from base states
to patterns through a saddle-node bifurcation. To analyze how the patterns
emerge as the system goes through this bifurcation, we performed numerical
simulations where we adiabatically increased the average mass by a global
cytosolic source with rate $\kappa_{\mathrm{s}}$

$\partial_{t}c=D_{c}\,\partial_{x}^{2}c-f(m,c)+\kappa_{\mathrm{s}},$ (23)

which entails $\partial_{t}\bar{n}=\kappa_{\mathrm{s}}$. In this appendix, we
explore the emergence of patterns when the average mass is non-adiabatically
increased beyond the base-state bifurcation. We initialize the system at zero
mass and increase mass up to a value $\bar{n}_{\mathrm{f}}$ within the regime
where no base state exists. Varying the rate of mass inflow
$\kappa_{\mathrm{s}}$, we find six regimes with qualitatively different
transient dynamics (see Supplementary Movies 9-14). Below, we briefly describe
the observed dynamics of the density profile in real space in these regimes
going from slow to fast $\kappa_{\mathrm{s}}$. We then use the $(m,c)$-phase
space and the concept of regional instability (cf. Sec. II.3) to heuristically
explain the observed dynamics.

  1. 1.

_Template-edge peak._ For small $\kappa_{\mathrm{s}}$, a density peak emerges
at the template edge, even though the density profile does not relax to a
quasi-steady state. This highlights that the edge-sensing mechanism is robust
against rate of mass inflow. (Supplementary Movie 9)

  2. 2.

_Transition regime._ Here, we observe two peaks emerging simultaneously, one
at the outer boundary of subdomain A and one at the template edge. This is an
intermediate regime between the boundary peak regime (3) and the template edge
peak regime (1), as both peaks emerge simultaneously. (Supplementary Movie 10)

  3. 3.

_System-boundary peak._ Here, we observe one peak forming at the system
boundary at $x=0$. (Supplementary Movie 11)

  4. 4.

_Sequential peak formation._ First, a peak forms at the domain boundary of
subdomain A, as in the system- boundary peak regime (3). Then, after the first
peak already formed, another peak forms at the template edge as in regime
(1).(Supplementary Movie 12)

  5. 5.

_Multiple peaks._ Here, the pattern-formation process is very similar to the
system-boundary peak regime (3). However, multiple peaks form at the outer
boundary of subdomain A simultaneously. (Supplementary Movie 13)

  6. 6.

_Quenched system._ Here the mass is upregulated almost instantaneously,
$\kappa_{\mathrm{s}}^{-1}{\rightarrow}\,0$, which is equivalent to a system
initialized with the complete mass in the cytosol. This results in a sequence
of peaks forms in subdomain A, starting from the template edge in a process
akin to front invasion into an unstable state van Saarloos (2003).
(Supplementary Movie 14)

Note that all states with multiple peaks (corresponding to multiple inflection
points of the pattern profile in subdomain A) are unstable due to coarsening.
The final steady state is always a pattern with a peak either at the template
edge or at the system boundary $x=0$.

In Sec. II.3 and Appendix B, we showed that the formation of a density peak is
determined by the position of a regional instability, when the system is in
quasi-steady state. The position of the regional instability can be found from
a phase portrait analysis, since the nullcline slope criterion (cf. Eq. (8))
determines which part of the density profile becomes unstable. To use the same
heuristic for understanding the formation of these transient peaks, we analyze
the density distribution in phase space obtained from numerical simulations in
these non-adiabatic regimes. When mass is added to the system on a faster
timescale than it can relax to its steady state, the density distribution in
phase space is no longer embedded in a single FBS. Instead, the density
distribution in phase space follows a ‘zig-zag’ shape, as illustrated in Fig.
10. This indicates that the density in the vicinity of the template edge is
still embedded in a FBS, with offset $\eta_{\mathrm{int}}$, but the density
far away from the template edge deviates from this FBS. Instead, these ‘quasi-
plateaus’ are slaved to the nullcline which indicates that they are locally
close to reactive equilibrium, and that their relaxation is limited by
diffusive mass transport. Accordingly, for faster inflow of mass, the zig-zag
shape is more pronounced—that is, the quasi-plateaus deviate more from
$\eta_{\mathrm{int}}$. If inflow of mass is faster than diffusive transport
across the quasi-plateaus, a region at the system boundary in subdomain A
enters the lateral unstable region in phase space first, as illustrated in
Fig. 10(b). This leads to the emergence of a peak at the system boundary in
regime (3), and to more complex pattern formation in regimes (4)–(6). In the
transition regime (2), mass inflow and mass transport roughly balance, such
that a region at the system boundary and a region at the template edge enter
the laterally unstable region in phase space at the same time.

The $(L^{2}/D_{c},\kappa_{\mathrm{s}}^{-1})$-phase diagram shown in Fig. 11
confirms the intuition that pattern emergence depends on the competition
between the time scales of mass inflow, $\kappa_{\mathrm{s}}$, and diffusive
mass transport across the system ${\sim}L^{2}/D_{c}$. Indeed, the regime
boundaries in the phase diagram are roughly straight lines emanating from the
origin. In particular, the transition from a ‘template-edge peak’ to a
‘system-boundary peak’ corresponds to a line $\kappa_{\mathrm{s}}^{-1}\approx
L^{2}/D_{c}$. This confirms the intuition that edge sensing is only possible
when the inflow of mass into the system is slower than the timescale of
diffusive mass transport. For comparison, in an intracellular context one has
$L\approx 10\;\mu\mathrm{m}$ and $D_{c}\approx
10\;\mu\mathrm{m}^{2}\mathrm{s}^{-1}$, such that the timescale of mass
transport across the cell is on the order of seconds. This is fast compared to
changes in average protein concentrations (for instance, due to protein
expression or release from the nucleus). Hence, the edge-sense mechanism is a
realistic candidate for template guided intracellular pattern formation.

Figure 10:  Effect of non-adiabatic mass upregulation on the pattern formation
dynamics. (a) Peak formation at the template edge: The density distribution is
not embedded in a single FBS, leading to the ‘zig-zag’-shaped density
distribution in phase space. The regional instability is still triggered at
the template edge, as highlighted by the (green) shaded region. (b) Peak
formation at the system boundary: The faster mass-inflow leads to a more
pronounced ‘zig-zag’-shaped density distribution in phase space. The regional
instability is now first triggered at the system boundary. This results in a
peak forming at the system boundary as shown in Fig. 4(b). Figure 11:  Phase
diagram for the timescale of mass upregulation (global cytosolic inflow
$\kappa_{\mathrm{s}}^{-1}$) against the timescale of mass-redistribution
across the entire domain ($L^{2}/D_{c}$). Edge sensing, i.e. formation of a
single density peak a the template edge, is operational in regime (1), see
Supplementary Movie 9. In regimes (2)-(6) a peak at the system boundary
($x=0$) or multiple peaks form, see Supplementary Movies 10-14.
Parameters: $D_{m}=0.01$, $D_{c}=10$, $k_{\mathrm{on}}=1$,
$\hat{k}_{\mathrm{fb}}=1$, $k_{\mathrm{off}}=2$, $\hat{K}_{\mathrm{d}}=1$,
$\theta_{\mathrm{A}}=0.5$, $\theta_{\mathrm{B}}=20$, $x_{E}=3/5\,L$,
$\bar{n}_{\mathrm{f}}=5$.

##  Appendix E Numerical continuation, steady state construction and
bifurcation scenarios

Numerical continuation. — To numerically calculate steady state solutions of
the two-component McRD system (Eq. (1)), we choose a finite-difference
discretization of the PDEs. For steady states, this yields an algebraic system
of equations that can be solved with an iterative Newton method. The basic
idea of numerical continuation is to follow a solution branch through
parameter space (see for instance, Ref. Krauskopf  _et al._ (2007) for an
excellent overview over continuation methods). This “path-following” is often
performed by emplying a predictor-–corrector scheme: Starting from one
solution, the next solution along the branch is predicted from the tangent
space of the solution branch which can be obtained from the Jacobian.

Steady state construction and finite domain size effects. — In order to test
the geometric constructions introduced in Sec. II we characterize the steady
state of the system with the quantity $m^{*}_{\mathrm{E}}-m_{\mathrm{B}}$,
which must be negative for low-mass base states (monotonic steady states in
the low-mass regime, cf. Sec. II.1) and positive for non-monotonic steady
states (stationary peak pattern localized at the template edge, cf. Sec.
II.2). We perform numerical continuation and compare the results from the
simulation (solid lines in Fig. 12) with the approximation from the geometric
construction (red dots and dash-dotted line). The geometric construction
serves as a good approximation for the steady for sufficiently large system
sizes.

Also note that for small system sizes the base state smoothly transitions into
the pattern state (purple line corresponding to $L=5$ in Fig. 12).

Figure 12:  One-parameter bifurcation structure in $\bar{n}$ (average mass)
connecting the low-mass base state (monotonic, i.e.
$m_{\mathrm{E}}<m_{\mathrm{B}}$ and the peak pattern at the template edge
(non-monotonic, i.e. $m_{\mathrm{E}}>m_{\mathrm{B}}$). Solid (dashed) lines
indicate stable (unstable) branches from numerical continuation for different
domain sizes (increasing from dark to light lines). Solutions from the
analytic construction of base states in the large domain size limit
($L\rightarrow\infty$) are shown as red dots and (position of the saddle-node
bifurcation denoted by $n_{\mathrm{bs}}^{\infty}$). Note that for small domain
size ($L=5$), the saddle-node bifurcations vanish and the base state smoothly
transitions into a stable peak pattern upon increasing $\bar{n}$. The red,
dash dotted line indicates the analytically constructed edge-localized pattern
(limit $L\rightarrow\infty$). The lower bound in average mass for the
existence of these patterns is denoted by $n_{\mathrm{pattern}}^{\infty}$.
Fixed parameters as in Fig. 4(a).

Bifurcation scenario. — The bifurcation scenario connecting the base state and
the patterns can be understood as a series of imperfect subcritical pitchfork
bifurcations. The imperfection is caused by the template that breaks mirror
symmetry of the system. On a homogeneous domain (i.e. without a template), the
bifurcations from homogeneous steady state to patterns are subcritical
pitchfork bifurcation that become supercritical for small system sizes / large
wavenumbers Brauns  _et al._ (2018). A more detailed analysis of the
bifurcation scenario induced by the step-like template is left for future
work. One interesting starting point would be to analyze the two-parameter
bifurcation diagram in the $(\bar{n},\theta_{\mathrm{A}})$-plane, where the
line $\theta_{\mathrm{A}}=\theta_{\mathrm{B}}$ correspond to the homogeneous
domain. Alternatively, one can investigate the bifurcation scenario in the
template edge position (i.e. the $(\bar{n},x_{\mathrm{E}})$ parameter plane),
where $x_{\mathrm{E}}=0$ and $x_{\mathrm{E}}=L$ correspond to homogeneous
domains.

Unstable multi-interface patterns. — The dotted branches in the bifurcation
structure Fig. 4(a) correspond to patterns with multiple self-organized
interfaces (i.e. more than two inflection points in the spatial profile, since
the template edge enforces one inflection point at $x_{\mathrm{E}}$). Figure
13 shows spatial profiles at the numbered points in Fig. 4(a) representative
for the respective branches. The spiral structure of the bifurcation structure
reflects an increasing number of of peaks from the outside to the center of
the spiral. For the branches numbered 1-3 (4-6) the concentration difference
$m^{*}_{\mathrm{E}}-m_{\mathrm{B}}$ is positive (negative), corresponding to a
peak (trough) at the template edge.

All multi-interface patterns are unstable due to a competition for mass and
decay into one of the two stable patterns, with a peak either at the system
boundary or at the template edge, in a coarsening process.

Figure 13:  Spatial profiles representative of the unstable branches numbered
(1-6) in the $\bar{n}$-bifurcation diagram Fig. 4(a).

##  Appendix F Demonstration of the edge-sensing criterion

Figure 14:  Demonstration of the edge-sensing criterion for a phenomenological
model for Cdc42 pattern formation. (a) Nullclines of the reaction kinetics Eq.
(24) in subdomains A (orange) and B (blue), defined by the reaction rates in
Eq. (27). Multiple B-nullclines are shown for off-rate factors
$\gamma_{\mathrm{off}}=1,1.5,\ldots,5$ while the on-rate factor is fixed at
$\gamma_{\mathrm{on}}=8$. The critical values $\gamma_{\mathrm{off}}^{\pm}$
between which the A-nullcline is intersected at negative slope (i.e. the edge-
sensing criterion is fulfilled) are shown as dashed and dash-dotted lines
respectively. (b) Phase diagram of on- and off-rate factors
$(\gamma_{\mathrm{on}},\gamma_{\mathrm{off}})$. The shaded region shows the
regime where edge-sensing is observed in numerical simulations with
adiabatically increasing average mass (cf. Eq. (23)). Dashed red lines show
the curves $\gamma_{\mathrm{off}}^{\pm}(\gamma_{\mathrm{on}})$, cf. Eq. (29),
between which the edge-sensing criterion is fulfilled. (c) Quasi-stationary
density profiles obtained from numerical simulations with adiabatically
increasing average mass $\bar{n}$ illustrating the typical pattern emerging
from the base-states’s saddle-node bifurcation in the three regimes of the
$(\gamma_{\mathrm{on}},\gamma_{\mathrm{off}})$ phase diagram
($\gamma_{\mathrm{on}}=8$; top: $\gamma_{\mathrm{off}}=5.5,\bar{n}=2.9$;
$\gamma_{\mathrm{off}}=3,\bar{n}=2.35$; bottom:
$\gamma_{\mathrm{off}}=1.5,\bar{n}=2.24$.). Fixed parameters:
$k_{\mathrm{fb}}=1,k_{\mathrm{on}}=0.07,k_{\mathrm{off}}=1,D_{m}=10^{-4},D_{c}=0.1,L=1,x_{\mathrm{E}}=0.5,\kappa_{\mathrm{s}}=10^{-5}$.

To illustrate how the edge sensing criterion might help in modeling of
biological systems we consider a concrete biological example. Cdc42 pattern
formation is described by a phenomenological two-component model of the form
Eq. (1) with the reaction kinetics Mori  _et al._ (2008)

$f(m,c)=\left(k_{\mathrm{on}}+k_{\mathrm{fb}}\frac{m^{2}}{1+m^{2}}\right)c-k_{\mathrm{off}}m$
(24)

The nullcline $c^{*}(m)$ determined by $f(m,c^{*}(m))=0$ has a section of
negative slope for $k_{\mathrm{fb}}>8k_{\mathrm{on}}$ in the interval
$[m_{\mathrm{max}},m_{\mathrm{min}}]$ given by

$m_{\mathrm{min,max}}=\sqrt{\frac{1}{1+k_{\mathrm{on}}/k_{\mathrm{fb}}}\left(\frac{1\mp\nu}{2}-\frac{k_{\mathrm{on}}}{k_{\mathrm{fb}}}\right)}$
(25)

with

$\nu=\sqrt{1-8\frac{k_{\mathrm{on}}}{k_{\mathrm{fb}}}}.$ (26)

Since cytosolic diffusion is multiple orders of magnitude faster than
diffusion of membrane bound Cdc42, we consider the limit $D_{c}\gg D_{m}$
where the FBS-slope is zero. Then the criterium for lateral instability is
that the nullcline slope be negative $\partial_{m}c^{*}(m)<0$; and the edge-
sensing criterion is that nullclines must intersect at a point where one of
them has negative slope.

In cellular wound healing, Abr forms a ring of high concentration around the
wound edge, followed by the formation of a Cdc42 ring around the Abr ring
Vaughan  _et al._ (2011). Cdc42 then activates the actomyosin machinery that
drives the contraction of the cell membrane to close the wound. Protein
mutation on Abr and Cdc42 studies suggest that the high density Abr ring acts
as a template for Cdc42 Vaughan  _et al._ (2011). Furthermore, it has been
shown that Abr acts as both a nucleotide exchange factor (GEF) and a GTPase-
activating protein (GAP) for Cdc42 Vaughan  _et al._ (2011). We therefore
model the effect of Abr on the Cdc42 kinetics as a factor increasing both the
attachment rate $k_{\mathrm{on}}$ and the detachment rate $k_{\mathrm{off}}$.
The reaction terms in the two subdomains with low and high Abr density, are
thus defined as
$f_{\mathrm{A,B}}(m,c)=f(m,c;k_{\mathrm{on}}^{\mathrm{A,B}},k_{\mathrm{off}}^{\mathrm{A,B}})$
with

$\displaystyle k_{\mathrm{on}}^{\mathrm{A}}$ $\displaystyle=k_{\mathrm{on}},$
$\displaystyle\;k_{\mathrm{on}}^{\mathrm{B}}$
$\displaystyle=\gamma_{\mathrm{on}}k_{\mathrm{on}},$ (27) $\displaystyle
k_{\mathrm{off}}^{\mathrm{A}}$ $\displaystyle=k_{\mathrm{off}},$
$\displaystyle\;k_{\mathrm{off}}^{\mathrm{B}}$
$\displaystyle=\gamma_{\mathrm{off}}k_{\mathrm{off}},$ (28)

where the factors $\gamma_{\mathrm{on}},\gamma_{\mathrm{off}}>1$ encode the
relative enhancement of attachment and detachment rates in subdomain B
compared to subdomain A.

The edge-sensing criterion (nullclines intersect at a point where the
A-nullcline has negative slope) is fulfilled when the intersection point
$(m_{\mathrm{i}},c_{\mathrm{i}})$ of the two nullclines, defined by
$f_{\mathrm{A}}(m_{\mathrm{i}},c_{\mathrm{i}})=0\;\&\;f_{\mathrm{B}}(m_{\mathrm{i}},c_{\mathrm{i}})=0$,
lies in the range $[m_{\mathrm{max}},m_{\mathrm{min}}]$ (see Fig. 14(a)).
Solving these equations for $\gamma_{\mathrm{off}}$, we find that nullclines
fulfill the edge-sensing criterion when $\gamma_{\mathrm{off}}$ lies in the
range $[\gamma_{\mathrm{off}}^{-},\gamma_{\mathrm{off}}^{+}]$ as a function of
$\gamma_{\mathrm{on}}$:

$\gamma_{\mathrm{off}}^{\pm}(\gamma_{\mathrm{on}})=\frac{1}{1+k_{\mathrm{fb}}/k_{\mathrm{on}}}\left(\frac{k_{\mathrm{fb}}}{k_{\mathrm{on}}}+\frac{2}{\nu\mp
1}+\gamma_{\mathrm{on}}\frac{\nu\mp 3}{\nu\mp 1}\right).$ (29)

Note that $\gamma_{\mathrm{off}}^{\pm}(1)=1$ corresponds to the singular case
where the template has no effect, i.e. the reaction kinetics in the two
subdomains are identical. The dashed, red lines in Fig. 14(b) show the curves
$\gamma_{\mathrm{off}}^{\pm}(\gamma_{\mathrm{on}})$ that delineate the regime
where the edge-sensing criterion is fulfilled.

To test the criterion, we ran numerical simulations, for different
combinations $(\gamma_{\mathrm{on}},\gamma_{\mathrm{off}})$. In each
simulation, the average mass is adiabatically increased at a rate
$\kappa_{\mathrm{s}}=10^{-5}$, cf. Eq. (23).

We find three different types of patterns emerging from the base-states’s
saddle-node bifurcation. Examples of these are shown in Fig. 14(c). The regime
where we find edge-localized peaks agrees well with the prediction from the
geometric edge-sensing criterion (shown as dashed, red lines in (b)).

Interestingly, the $(\gamma_{\mathrm{on}},\gamma_{\mathrm{off}})$-phase
diagram shows that both attachment and detachment need to be enhanced in
subdomain B. This implies that proteins, like Abr, that have both GEF- and
GAP-catalytic domains may play a crucial role for edge sensing by GTPases.

## References

  * Halatek  _et al._ (2018) J. Halatek, F. Brauns, and E. Frey, Philosophical Transactions of the Royal Society B: Biological Sciences 373, 20170107 (2018).
  * Howard  _et al._ (2001) M. Howard, A. D. Rutenberg, and S. de Vet, Physical Review Letters 87 (2001), 10.1103/PhysRevLett.87.278102.
  * Huang  _et al._ (2003) K. C. Huang, Y. Meir, and N. S. Wingreen, Proceedings of the National Academy of Sciences 100, 12724 (2003).
  * Mori  _et al._ (2008) Y. Mori, A. Jilkine, and L. Edelstein-Keshet, Biophysical Journal 94, 3684 (2008).
  * Goryachev and Pokhilko (2008) A. B. Goryachev and A. V. Pokhilko, FEBS Letters 582, 1437 (2008).
  * Ishihara  _et al._ (2007) S. Ishihara, M. Otsuji, and A. Mochizuki, Physical Review E 75 (2007), 10.1103/PhysRevE.75.015203.
  * Otsuji  _et al._ (2007) M. Otsuji, S. Ishihara, C. Co, K. Kaibuchi, A. Mochizuki, and S. Kuroda, PLoS Computational Biology 3, e108 (2007).
  * Halatek and Frey (2012) J. Halatek and E. Frey, Cell Reports 1, 741 (2012).
  * Klünder  _et al._ (2013) B. Klünder, T. Freisinger, R. Wedlich-Söldner, and E. Frey, PLOS Computational Biology 9, e1003396 (2013).
  * Trong  _et al._ (2014) P. K. Trong, E. M. Nicola, N. W. Goehring, K. V. Kumar, and S. W. Grill, New Journal of Physics 16, 065009 (2014).
  * Alonso and Bär (2014) S. Alonso and M. Bär, EPJ Nonlinear Biomed Phys 2, 1 (2014).
  * Wu  _et al._ (2016) F. Wu, J. Halatek, M. Reiter, E. Kingma, E. Frey, and C. Dekker, Molecular Systems Biology 12, 873 (2016).
  * Goryachev and Leda (2017) A. B. Goryachev and M. Leda, Molecular Biology of the Cell 28, 370 (2017).
  * Murray and Sourjik (2017) S. M. Murray and V. Sourjik, Nature Physics 13, 1006 (2017).
  * Denk  _et al._ (2018) J. Denk, S. Kretschmer, J. Halatek, C. Hartl, P. Schwille, and E. Frey, Proc Natl Acad Sci USA 115, 4553 (2018).
  * Cusseddu  _et al._ (2018) D. Cusseddu, L. Edelstein-Keshet, J. Mackenzie, S. Portet, and A. Madzvamuse, Journal of Theoretical Biology (2018), 10.1016/j.jtbi.2018.09.008.
  * Chiou  _et al._ (2018) J.-G. Chiou, S. A. Ramirez, T. C. Elston, T. P. Witelski, D. G. Schaeffer, and D. J. Lew, PLOS Computational Biology 14, e1006095 (2018).
  * Glock  _et al._ (2019) P. Glock, F. Brauns, J. Halatek, E. Frey, and P. Schwille, bioRxiv (2019), 10.1101/666362.
  * Turing (1952) A. M. Turing, Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences 237, 37 (1952).
  * Williams  _et al._ (2019) T. D. Williams, P. I. Paschke, and R. R. Kay, Philosophical Transactions of the Royal Society B: Biological Sciences 374, 20180150 (2019).
  * Veltman  _et al._ (2016) D. M. Veltman, T. D. Williams, G. Bloomfield, B.-C. Chen, E. Betzig, R. H. Insall, and R. R. Kay, eLife 5, e20085 (2016).
  * Vaughan  _et al._ (2011) E. M. Vaughan, A. L. Miller, H.-Y. E. Yu, and W. M. Bement, Current Biology 21, 270 (2011).
  * Chant and Pringle (1995) J. Chant and J. R. Pringle, J. Cell Biol. 129, 751 (1995).
  * Tong  _et al._ (2007) Z. Tong, X.-D. Gao, A. S. Howell, I. Bose, D. J. Lew, and E. Bi, J Cell Biol 179, 1375 (2007).
  * Lo  _et al._ (2013) W.-C. Lo, M. E. Lee, M. Narayan, C.-S. Chou, and H.-O. Park, PLoS ONE 8, e56665 (2013).
  * Miller  _et al._ (2017) K. E. Miller, W.-C. Lo, M. E. Lee, P. J. Kang, and H.-O. Park, Molecular Biology of the Cell 28, 3773 (2017).
  * Halatek and Frey (2018) J. Halatek and E. Frey, Nature Physics 14, 507 (2018).
  * Brauns  _et al._ (2018) F. Brauns, J. Halatek, and E. Frey, arXiv:1812.08684 (2018), arXiv:1812.08684 .
  * Maini  _et al._ (1992) P. K. Maini, D. L. Benson, and J. A. Sherratt, Math Med Biol 9, 197 (1992).
  * Ermentrout and Rinzel (1996) G. B. Ermentrout and J. Rinzel, SIAM Journal on Applied Mathematics 56, 1107 (1996).
  * Bär  _et al._ (1996) M. Bär, A. K. Bangia, I. G. Kevrekidis, G. Haas, H.-H. Rotermund, and G. Ertl, The Journal of Physical Chemistry 100, 19106 (1996).
  * Yuan  _et al._ (2007) X. Yuan, T. Teramoto, and Y. Nishiura, Phys. Rev. E 75, 036220 (2007).
  * Nishiura  _et al._ (2007) Y. Nishiura, T. Teramoto, X. Yuan, and K.-I. Ueda, Chaos: An Interdisciplinary Journal of Nonlinear Science 17, 037104 (2007).
  * Miyazaki and Kinoshita (2007) J. Miyazaki and S. Kinoshita, Physical Review E 76 (2007), 10.1103/PhysRevE.76.066201.
  * Nishi  _et al._ (2013) K. Nishi, Y. Nishiura, and T. Teramoto, Japan Journal of Industrial and Applied Mathematics 30, 351 (2013).
  * Doelman  _et al._ (2016) A. Doelman, P. van Heijster, and F. Xie, SIAM Journal on Applied Dynamical Systems 15, 655 (2016).
  * van Heijster  _et al._ (2019) P. van Heijster, C.-N. Chen, Y. Nishiura, and T. Teramoto, J Dyn Diff Equat 31, 153 (2019).
  * Scheel and Weinburd (2018) A. Scheel and J. Weinburd, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 376, 20170191 (2018).
  * Jones (1995) C. K. R. T. Jones, in  _Dynamical Systems_, Vol. 1609, edited by R. Johnson (Springer, Berlin, Heidelberg, 1995) pp. 44–118.
  * Ward (2006) M. J. Ward, Bulletin of Mathematical Biology 68, 1151 (2006).
  * Guckenheimer and Holmes (1983) J. Guckenheimer and P. Holmes,  _Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields_, Applied Mathematical Sciences, Vol. 42 (Springer New York, New York, NY, 1983).
  * Krauskopf  _et al._ (2007) B. Krauskopf, H. M. Osinga, and J. Galan-Vioque, eds.,  _Numerical Continuation Methods for Dynamical Systems: Path Following and Boundary Value Problems_, Understanding Complex Systems (Springer Netherlands, Dordrecht, 2007).
  * Pryce (1993) J. D. Pryce,  _Numerical Solution of Sturm-Liouville Problems_ , Monographs on Numerical Analysis (Clarendon Press, Oxford; New York, 1993).
  * Ziepke  _et al._ (2016) A. Ziepke, S. Martens, and H. Engel, J. Chem. Phys. 145, 094108 (2016).
  * Pfeffer (2012) S. R. Pfeffer, Biochemical Society Transactions 40, 1373 (2012).
  * Novick (2016) P. Novick, Small GTPases 7, 252 (2016).
  * Noack and Jaillais (2017) L. C. Noack and Y. Jaillais, Current Opinion in Plant Biology 40, 22 (2017).
  * Rulands  _et al._ (2013) S. Rulands, B. Klünder, and E. Frey, Phys. Rev. Lett. 110, 038102 (2013).
  * Zadorin  _et al._ (2017) A. S. Zadorin, Y. Rondelez, G. Gines, V. Dilhas, G. Urtel, A. Zambrano, J.-C. Galas, and A. Estevez-Torres, Nature Chemistry 9, 990 (2017).
  * Thalmeier  _et al._ (2016) D. Thalmeier, J. Halatek, and E. Frey, Proc Natl Acad Sci USA 113, 548 (2016).
  * Ziepke  _et al._ (2019) A. Ziepke, S. Martens, and H. Engel, SIAM J. Appl. Dyn. Syst. 18, 1015 (2019).
  * Haupt and Minc (2018) A. Haupt and N. Minc, Journal of Cell Science 131, jcs214015 (2018).
  * van Saarloos (2003) W. van Saarloos, Physics Reports 386, 29 (2003).