Phase-space geometry of mass-conserving reaction-diffusion dynamics

TL;DR

This paper extends phase space analysis to mass-conserving reaction-diffusion systems, providing geometric insights into pattern formation, bifurcations, and the role of nonlinearities, with implications for understanding protein-pattern dynamics.

Contribution

It introduces a comprehensive phase-space framework for two-component mass-conserving reaction-diffusion systems, linking local equilibria, geometric criteria, and pattern bifurcations.

Findings

Pattern formation as a mass-redistribution instability

Reactive nullcline shape influences global dynamics

Bifurcation analysis based on regional dispersion relations

Abstract

Experimental studies of protein-pattern formation have stimulated new interest in the dynamics of reaction-diffusion systems. However, a comprehensive theoretical understanding of the dynamics of such highly nonlinear, spatially extended systems is still missing. Here we show how a description in phase space, which has proven invaluable in shaping our intuition about the dynamics of nonlinear ordinary differential equations, can be generalized to mass-conserving reaction-diffusion (McRD) systems. We present a comprehensive analysis of two-component McRD systems, which serve as paradigmatic minimal systems that encapsulate the core principles and concepts of the local equilibria theory introduced in the paper. The key insight underlying this theory is that shifting local (reactive) equilibria -- controlled by the local total density -- give rise to concentration gradients that drive diffusive redistribution of total density. We show how this dynamic interplay can be embedded in the phase plane of the reaction kinetics in terms of simple geometric objects: the reactive nullcline and the diffusive flux-balance subspace. On this phase-space level, physical insight can be gained from geometric criteria and graphical constructions. The effects of nonlinearities on the global dynamics are simply encoded in the curved shape of the reactive nullcline. In particular, we show that the pattern-forming `Turing instability' in McRD systems is a mass-redistribution instability, and that the features and bifurcations of patterns can be characterized based on regional dispersion relations, associated to distinct spatial regions (plateaus and interfaces) of the patterns. In an extensive outlook section, we detail concrete approaches to generalize local equilibria theory in several directions, including systems with more than two-components, weakly-broken mass conservation, and active matter systems.

Equations

$${\partial }_{t}m(x,t)$$

$${\partial }_{t}c(x,t)$$

$$\bar{n}=\frac{1}{L}\int_{0}^{L}\kern-1.99997pt\mathrm{d}x\,\big{(}m(x,t)+c(x,t)\big{)}.$$

$$f_{\text{attach–detach}}(m,c)=a(m)c-d(m)m,$$

$${\partial }_{t}m=f(m,c),{\partial }_{t}c=-f(m,c).$$

$$\mathbf{u}^{*}(n):\left\{\begin{array}[]{rl}f(\mathbf{u}^{*})=&0,\\[0.86108pt] m^{*}+c^{*}=&n,\end{array}\right.$$

$$\begin{split}f(m,c)&\approx\sigma_{\text{loc}}(n)\cdot\big{[}m-m^{*}(n)\big{]}\\ &=-\sigma_{\text{loc}}(n)\cdot\big{[}c-c^{*}(n)\big{]}.\end{split}$$

$$s_{\text{nc}}(n):=\partial_{m}c^{*}(m)\Big{|}_{n}=-\frac{f_{m}}{f_{c}}\Big{|}_{n}.$$

$$s_{\text{nc}}(n)=-1.$$

$$D_m{\partial }_x^2\tilde{m}+f(\tilde{m},\tilde{c})$$

$$D_c{\partial }_x^2\tilde{c}-f(\tilde{m},\tilde{c})$$

$$D_m{\partial }_x\tilde{m}(x)=-D_c{\partial }_x\tilde{c}(x).$$

$$\frac{D_m}{D_c}\tilde{m}(x)+\tilde{c}(x)={\eta }_0,$$

$$\eta :=\frac{D_m}{D_c}m+c,$$

$${\partial }_{t}n(x,t)=D_c{\partial }_x^2\eta (x,t).$$

$$\tilde{f}(m,\eta ):=f\left(m,\eta -m{D}_m/D_c\right),$$

$$D_m{\partial }_x^2\tilde{m}+\tilde{f}(\tilde{m},\tilde{\eta })$$

$$D_c{\partial }_x^2\tilde{\eta }$$

$$\begin{split}f\big{(}\mathbf{u}^{*}(n(x,t))\big{)}&=0,\\ m^{*}(n(x,t))+c^{*}(n(x,t))&=n(x,t),\end{split}$$

$${\int }_0^L\mathrm{d}x\tilde{f}(\tilde{m}(x),{\eta }_0)=0.$$

$${\int }_{\tilde{m}(0)}^{\tilde{m}(L)}\mathrm{d}m\tilde{f}(m,{\eta }_0)=0.$$

$$\tilde{f}\big{(}m,\eta_{0}\big{)}\approx\sigma_{\text{loc}}(\widetilde{n}(m))\cdot\big{[}m-m^{*}(\widetilde{n}(m))\big{]},$$

$$\int_{\widetilde{m}(0)}^{\widetilde{m}(L)}\kern-1.99997pt\mathrm{d}m\,\big{[}m-m^{*}\big{(}\widetilde{n}(m)\big{)}\big{]}=0.$$

$$\begin{split}\frac{1}{2}D_{m}\big{(}\partial_{x}\widetilde{m}\big{|}_{x_{0}}\big{)}^{2}&=\mathop{\hphantom{-}}\int_{\widetilde{m}(0)}^{m_{0}(\eta_{0})}\kern-1.99997pt\mathrm{d}m\,\tilde{f}(m,\eta_{0})\\ &=-\int_{m_{0}(\eta_{0})}^{\widetilde{m}(L)}\kern-1.99997pt\mathrm{d}m\,\tilde{f}(m,\eta_{0}).\end{split}$$

$${\partial }_{t}n(x,t)=D_c{\partial }_x^{2}c(x,t)+D_m{\partial }_x^{2}m(x,t).$$

$$\mathbf{u}(x,t)\to\mathbf{u}^{*}\big{(}n(x,t)\big{)},$$

$${\partial }_{t}n(x,t)\approx D_c{\partial }_x^2{c}^{\ast }(n)+D_m{\partial }_x^2{m}^{\ast }(n).$$

$$\partial_{t}n(x,t)\approx\partial_{x}\big{[}\big{(}D_{c}\,\partial_{n}c^{*}+D_{m}\,\partial_{n}m^{*}\big{)}\partial_{x}n\big{]}$$

$$\frac{{\partial }_n{c}^{\ast }}{{\partial }_n{m}^{\ast }}=s_{\text{nc}}(n)=-\frac{f_m}{f_c}<-\frac{D_m}{D_c},$$

$$s_{\text{nc}}(\bar{n})<-\frac{D_m}{D_c}.$$

$$\partial_{t}n(x,t)\approx D_{c}\partial_{x}\big{[}\partial_{n}\eta^{*}(n)\,\partial_{x}n\big{]}.$$

$$q_{\max }^2:=\frac{f_m}{D_m}-\frac{f_c}{D_c}=\frac{{\tilde{f}}_m}{D_m}$$

Full text

mass conservation; reaction–diffusion; pattern formation; phase space; local equilibria; flux balance; self-organization

I Introduction

I.1 Motivation and background

Nonlinear systems are as prevalent in nature as they are difficult to deal with conceptually and mathematically Strogatz (2014); Izhikevich (2007); Mikhailov (1994); Murray (2002); Jackson (1989); Cross and Greenside (2009); Ott (2002). Cases in which the equations describing such systems can be solved in closed analytical form are rare, making nonlinear problems appear inaccessible to mathematical analysis at first sight. A key insight, going back to the work of Poincaré Poincaré (1899), was that geometric structures in the phase space of a system can provide qualitative information about the global dynamics (trajectories in phase space) without an explicit solution of the differential equations. The essence of this geometric reasoning can be understood by considering simple systems with only two independent variables; see e.g. Refs. Izhikevich (2007); Strogatz (2014). In this case, the key geometric objects are nullclines, defined as curves in phase space along which one of the system’s two variables is in equilibrium. The points at which nullclines intersect mark equilibria (fixed points) of the system. These geometric objects organize phase-space flow, and thereby allow one to infer the qualitative dynamics from the shapes and intersections of the nullclines. Key concepts like linear stability, excitability, multi-stability, and limit cycles can be understood in the context of such a geometric analysis Izhikevich (2007); Strogatz (2014). Transitions (bifurcations) between qualitatively different regimes are revealed by structural changes of the flow in phase space as the control parameters are varied. One key advantage of such a geometric approach to nonlinear dynamical systems is that it yields systematic physical insights into the processes driving dynamics without requiring the explicit solution of the full set of equations.

Generalizing these methods, developed for ordinary differential equations (ODEs), to phenomena that explicitly require a description on a spatially extended domain—and therefore involve partial differential equations (PDEs)—poses a huge, ongoing challenge. Instances where this has been successfully achieved are rare. One classical approach for nonlinear systems in one spatial dimension is the construction of steady-state patterns (including both stationary patterns and traveling-wave solutions in a co-moving frame). Mathematically, the steady state in this case is described by a set of ODEs, which can be analyzed based on their phase-space geometry (see e.g. Refs. Jones (1994); Kerner and Osipov (1994)). An elementary example of this is the phase-plane analysis of traveling waves of the Fisher-KPP equation Fisher (1937); Kolmogorov et al. (1937) as described in Ref. Jones et al. (2009). Here, we go beyond this approach, and gain physical insight into the global dynamics of spatially extended systems from the analysis of geometric objects in a low-dimensional phase space. Crucially, such a theory should be able to explain both, the dynamic process of pattern formation—initiated, for instance, by a lateral (Turing) instability—as well as the final stationary patterns in terms of the same concepts and principles.

In its full generality, this is likely a futile task. Here, we restrict ourselves to mass-conserving reaction–diffusion (MCRD) dynamics. A broad class of systems that can be described by MCRD dynamics are models for intracellular protein-pattern formation Halatek et al. (2018), which is essential for the spatiotemporal organization of many cellular processes, including cell division, motility and differentiation. Moreover, as we discuss in the Outlook, many pattern-forming systems are governed by a combination of mass-conserving dynamics and source terms that break mass conservation. Studying such systems in the nearly mass-conserving limit may help to tackle long-standing questions like pattern selection (wavelength selection) in the highly nonlinear regime Brauns et al. (2020a).

Recent results have indicated ways of making progress towards a general theory rooted in mass-conservation Halatek and Frey (2018). Based on numerical simulations, this study suggests a new way of thinking about pattern formation, namely in terms of mass redistribution that gives rise to moving local equilibria: A dissection of space into (notional) local compartments allows the spatiotemporal dynamics to be characterized on the basis of the ODE phase space of local reactions. As (globally conserved) masses are spatially redistributed, the local masses in the compartments act as parameters for the reactive phase-space flow within them. The properties (position and stability) of the local reactive equilibria in the compartments are shown to depend on local masses and thus act as proxies for the local phase-space flow Note:Onsager .Onsager (1931)

Diffusion acts to redistribute the conserved quantities between neighboring compartments and thereby induces changes in the local phase-space structure. This level of description proved to be very powerful in explaining chemical turbulence and transitions from chemical turbulence to long range order (standing and traveling waves) far from onset of the (subcritical) lateral instability. The prediction of chemical turbulence at onset, based on numerical simulations in Ref. Halatek and Frey (2018), was recently confirmed experimentally Denk et al. (2018). Hence, the advances in this “proof-of-principle” study Halatek and Frey (2018) suggest that a comprehensive theory of pattern formation in reaction–diffusion systems with conserved total densities (masses) can be developed based on the concept of mass-redistribution.

Here, we put this overarching idea on a general theoretical foundation. To this end, we develop a number of new theoretical concepts, exemplified by two-component mass-conserving reaction-diffusion (2C-MCRD) systems and based on simple geometric structures in the phase space of the reaction kinetics. From these concepts, general geometric criteria for lateral (Turing) instability and stimulus-induced pattern formation emerge and allow us to obtain the features and bifurcations of patterns from graphical constructions. Moreover, these advances reveal connections to other pattern forming phenomena like liquid-liquid phase separation and shear banding in complex fluids.

In our work, general two-component systems serve two purposes. First, as a paradigmatic and didactic example that encapsulates the core concepts and principles of our theory in a pedagogic and broadly accessible way. Second, they provide an elementary basis for further generalizations. Taken together, the role we envision for 2C-MCRD systems is similar to the role of two-variable systems in dynamical systems theory of ODEs.

The framework we present here has recently been employed and generalized to study the principles underlying coarsening and wavelength selection Brauns et al. (2020a), as well as systems with spatially heterogeneous reaction rates Wigbers et al. (2020a) and the role of advective flow in the cytosol Wigbers et al. (2020b)—all in the context of two-component systems. Moreover, the concepts of local equilibria and regional instabilities have recently proven useful to disentangle the interplay of several distinct instabilities that give rise to Min-protein patterns in vivo and in vitro Brauns et al. (2020b).

Potential future generalizations range from adding more components and more conserved quantities, to going beyond strictly mass-conserving systems (see Outlook). Mass conservation and, more generally, conserved quantities are inherent to the elementary processes underlying many pattern forming systems. We therefore believe that the local equilibria theory we present here offers a new perspective on a broad class of pattern-forming systems—including intracellular pattern formation, classical chemical systems such as the BZ reaction, and even agent-based active matter systems.

I.2 Structure of the paper

Put briefly, the paper is structured as follows. Section II introduces 2C-MCRD systems, and their applications, most prominently as conceptual models for cell polarization. The concepts introduced in Sec. III and Sec. IV form the foundation of the proposed framework and the subsequent analysis. The following two sections present results that are particularly relevant in the biological context of intracellular pattern formation: A characterization of the possible pattern types exhibited by 2C-MCRD systems (Sec. V) and a simple heuristic for the threshold perturbation required for stimulus-induced pattern formation (Sec. VI). Section VII delves into a more technical analysis of the generic bifurcation structure of 2C-MCRD systems. Here, we find striking similarities to the phase diagram of phase-separation phenomena (such as liquid-liquid phase separation and motility-induced phase separation). This technical section also includes weakly nonlinear analysis in the vicinity of the lateral instability onset that corresponding to a critical point in the language of phase transitions. Finally, in Sec. VIII, we provide an extensive discussion of the implications of our results and an outlook to future research directions.

II Two-component mass-conserving reaction–diffusion systems

Our goal is to find geometric structures in phase space that allow the characterization of mass-conserving reaction–diffusion (MCRD) systems. The simplest system of this type is a two-component reaction–diffusion system with two scalar densities, $m(x,t)$ and $c(x,t)$,

[TABLE]

on a one-dimensional domain of length $L$ with reflective (no-flux) boundary conditions $\partial_{x}m|_{0,L}=\partial_{x}c|_{0,L}=0$; all results can straightforwardly be generalized to periodic boundary conditions. The global average $\bar{n}$ of total density $n(x,t)=m(x,t)+c(x,t)$ is conserved:

[TABLE]

We chose the above form for its conceptual simplicity. However, the principles that characterize pattern formation for this ‘minimal’ model can be generalized to more complex systems with more components and conserved species Halatek and Frey (2018); Brauns et al. (shed), and even beyond strictly mass-conserving systems Brauns et al. (2020a); see also Sec. VIII.4.

Two-component systems of the above form were widely studied as conceptual models for cell polarization Ishihara et al. (2007); Otsuji et al. (2007); Goryachev and Pokhilko (2008); Altschuler et al. (2008); Mori et al. (2008); Jilkine and Edelstein-Keshet (2011); Jilkine et al. (2011); Edelstein-Keshet et al. (2013); Trong et al. (2014); Seirin Lee and Shibata (2015); Chiou et al. (2018); Diegmiller et al. (2018); Hubatsch et al. (2019), where Eq. (1) describes the dynamics of a protein species that cycles between membrane (slow diffusing, concentration $m(x,t)$) and cytosol (fast diffusing, concentration $c(x,t)$). In this biological context, the nonlinear kinetics term $f(m,c)$ is of the form

[TABLE]

where the non-negative terms $a(m)c$ and $d(m)m$ characterize the attachment of proteins from the cytosol to the membrane and the detachment back into the cytosol, respectively. This functional form results from the fact that in intracellular systems chemical reactions are mainly restricted to the cell membrane. We will use kinetics of the above form for illustration purposes; for specific examples see Appendix A. Importantly however, our results hold for general kinetics $f$, and are not restricted to models for intracellular pattern formation. Moreover, 2C-MCRD systems have also been studied for slime mold aggregation Keller and Segel (1970), cancer cell migration (glioma invasion) Pham et al. (2012), precipitation patterns Scheel (2009); Hilhorst et al. (2006), and simple contact processes van Wijland et al. (1998); Kessler and Levine (1998). Finally, non-isothermal solidification models Caginalp (1986) can also be rewritten in the form Eq. (1); see e.g. Refs. Morita and Ogawa (2010); Pogan and Scheel (2012).

In the mathematical literature, 2C-MCRD systems with a specific form of the reaction kinetics, $f(m,c)=c-g(m)$, have been studied extensively Morita and Ogawa (2010); Goh et al. (2011); Jimbo and Morita (2013); Latos et al. (2018). The dynamics of these systems can be mapped to a variational form (gradient flow of an effective free-energy density). In this form, the properties of the dynamics and the stationary patterns can be analyzed analogously to the Cahn–Hilliard equation which describes phase separation near thermal equilibrium (see e.g. Ref. Pismen (2006)). In particular, one can prove that these systems always exhibit uninterrupted coarsening, i.e. the fully phase separated state is the only stable stationary state of the system Jimbo and Morita (2013); Latos et al. (2018). The theory we present here is fundamentally different from these previous mathematical approaches. Instead of an abstract mapping to a variational form, our approach is grounded in concepts with clear physical interpretation that are not restricted to specific reaction terms. Local equilibria, the overarching concept of our theory, can be generalized to systems with more than two components and more complex phenomena such as waves, oscillations and chaos Halatek and Frey (2018); Brauns et al. (2020b).

In closing this section we we would like to point out that systems which are not strictly mass conserving may have a mass-conserving subsystem (or ‘core’) that captures essential aspects of the system’s pattern formation dynamics. An example is the Brusselator system Prigogine and Lefever (1968), a widely used paradigmatic model to study pattern formation. It’s reaction kinetics has a ‘core’ of the same form as Eq. (1) with additional linear production and degradation terms that break mass conservation. This broken mass conservation can give rise to interesting new phenomena, not present in the mass-conserving core. Still, the core dynamics can be useful in understanding these new phenomena by exploiting a time scale separation between (fast) mass-conserving processes and (slow) production/degradation processes Kuwamura and Morita (2015); Kuwamura and Izuhara (2017); Brauns et al. (2020a). In the Outlook, Sec. VIII.4, we briefly discuss this example and the broader prospects of such an approach to non-conservative systems.

III Setting the stage—geometric structures in phase space

In this section we introduce the basic geometric concepts in phase space, which we will later use for a full characterization of pattern formation, including pattern types, bifurcations, and the corresponding characteristic length and time scales. To this end, we will first study the spatially homogeneous (well-mixed) case where we can use the classical geometric tools for studying ordinary differential equations Strogatz (2014); Izhikevich (2007); Jackson (1989). Subsequently, we will build on the phase-space structures obtained from the well-mixed case to also understand pattern formation in spatially extended systems in terms of flow in phase space.

III.1 Phase-space analysis of a well-mixed system

For a well-mixed system, the dynamics reduces to a set of ordinary differential equations

[TABLE]

Since the reaction kinetics conserve total density (protein mass), $n=m+c$, the reactive flow in $(m,c)$-phase space is restricted to the reactive phase space where $n$ is a constant of motion (i.e. $\partial_{t}n=0$), as illustrated in Fig. 1a. The reaction kinetics are balanced at the reactive equilibria (fixed points) $\mathbf{u}^{*}=(m^{*},c^{*})$,

[TABLE]

which are given by the intersection points between the reactive reactive nullcline (NC) (or ‘line of reactive equilibria’), $f(m^{*},c^{*})=0$, and the reactive phase space for a given mass $n$. The reactive flow in the respective phase space is organized by the location and (linear) stability of these fixed points (which are both functions of $n$); see Fig. 1a and discussion below. By varying the total density $n$, i.e. by shifting the reactive phase space, one can construct a bifurcation diagram of the (reactive) equilibrium $c^{*}(n)$ as a function of the total density $n$ (Fig. 1b). The total density is a control parameter of the reactive equilibria. When the total density changes, the local equilibria shift. These shifting (or moving) local equilibria, introduced in Ref. Halatek and Frey (2018), are the key to understanding the mass-conserving reaction–diffusion dynamics as we will see repeatedly throughout this paper.

In each reactive phase space, we can eliminate (the cytosolic density) $c(t)$ and write the dynamics in terms of (the membrane density) $m(t)$ alone: \partial_{t}m(t)=f\big{(}m(t),n-m(t)\big{)}; equally well $m(t)$ could be eliminated. In the vicinity of an equilibrium $m^{*}$ the linearized reactive flow reads $\partial_{t}m(t)\approx(f_{m}-f_{c})(m(t)-m^{*})$, with $f_{i}:=\partial_{i}f|_{\mathbf{u}^{*}},i\in\{m,c\}$. We can read off the eigenvalue $\sigma_{\text{loc}}(n):=f_{m}-f_{c}$ for the local equilibrium and obtain to linear order in the vicinity of the reactive nullcline:

[TABLE]

The sign of $\sigma_{\text{loc}}$—and thereby stability of the reactive equilibria—can be inferred from the slope of the nullcline

[TABLE]

For $f_{c}>0$, which is always the case for attachment–detachment kinetics where $f_{c}=a(m)$, local equilibria are stable, $\sigma_{\text{loc}}(n)<0$, if (and only if) the slope of the reactive nullcline is less steep than the slope of the reactive phase space:

[TABLE]

Figure 1a shows an example for reaction kinetics where the dynamics is monostable except for a window of protein masses exhibiting bistability with one unstable ($\circ$) and two stable fixed points ($\bullet$). (Note that the local eigenvalue $\sigma_{\text{loc}}$ can be rewritten as $\sigma_{\text{loc}}=f_{c}\cdot(-s_{\text{nc}}-1)$, which shows why the slope criterion, Eq. (7), for local stability is reversed for $f_{c}<0$.)

III.2 Stationary patterns are embedded in a flux-balance subspace of phase space

To generalize the above approach to spatially extended systems, one has to understand the role of diffusive coupling. We start by studying stationary patterns. The insights gained from this analysis will later prove useful for studying the dynamics (instability of the homogeneous state and stimulus-induced pattern formation).

A stationary pattern $\widetilde{\mathbf{u}}(x)=[\widetilde{m}(x)$, $\widetilde{c}\kern 0.74997pt(x)]$ is a solution to the steady-state equations,

[TABLE]

under the constraint of a given average total density $\bar{n}$ (Eq. (1c)). Figure 2a shows the sketch of a typical stationary pattern $\widetilde{m}(x)$ (solid line) and the corresponding local equilibria $m^{*}(x)$ (black disks) obtained from a numerical solution of Eq. (8) with a reaction term as, for instance, in Refs. Mori et al. (2008, 2011); Chiou et al. (2018) (see Appendix A). We study patterns with monotonic density profiles, which serve as elementary building blocks for more complex stationary patterns (see Section V).

Figure 2a shows an example for a monotonic pattern profile exhibiting two plateau regions connected by an interface region with an inflection point at $x_{0}$. (We will later see that this type of pattern, termed ‘mesa’, is one of three elementary pattern types found in two-component reaction–diffusion systems; see Sec. V.3).

Here we ask what can be learned about the stationary pattern by applying geometric concepts in phase space alone, i.e. without relying on an explicit numerical solution. Observe that Eqs. (8) imply that the diffusive fluxes of $m$ and $c$ have to balance locally at each position $x$ in the spatial domain $[0,L]$ in steady state:

[TABLE]

This local flux-balance condition is obtained by adding the two steady-state equations, Eq. (8a) and Eq. (8b), integrating over $x$ and employing no-flux boundary conditions. Integrating this relation once more from the boundary to any point $x$ in the domain yields that any stationary pattern obeys the linear relation

[TABLE]

where $\eta_{0}$ is a constant of integration. (An alternative derivation of this relation that directly generalizes to higher spatial dimensions (cf. Eq. (14b) is provided below.) Equation (10) defines a family of linear subspaces in the $(m,c)$-phase plane parametrized by $\eta_{0}$. We shall call these subspaces the flux-balance subspaces (FBS), since they represents the local balance between the diffusive fluxes on the membrane and in the cytosol. Any stationary pattern is confined to one such subspace; see Fig. 2b. We will learn later (Sec. III.3), how the value of $\eta_{0}(\bar{n},L)$ is determined by the balance of reactive processes.

Equation (10) has been previously used to mathematically simplify the construction and analysis of stationary patterns in two-component systems, by introducing the new phase-space coordinate (orthogonal to the flux-balance subspace)

[TABLE]

and describing the spatiotemporal dynamics in terms of the scalar fields $n(x,t)$ and $\eta(x,t)$ (cf. Refs. Otsuji et al. (2007); Chiou et al. (2018)). The physical origin (diffusive flux balance) and the geometric interpretation (flux-balance subspace) discussed above explains why Eq. (10) has proven to be useful before (and why it will be central in our further analysis). In particular, note that by adding Eqs. (1a) and (1b) one finds that gradients in $\eta(x,t)$ drive mass redistribution:

[TABLE]

We will therefore call $\eta(x,t)$ the mass-redistribution potential. Substituting $c$ using $\eta$, the reaction term reads

[TABLE]

and the stationarity conditions, Eqs. (8a) and (8b), are replaced by

[TABLE]

From the second equation, we recover that in steady state, the mass-redistribution potential must be constant in space, $\widetilde{\eta}=\eta_{0}$, on a domain with no-flux or periodic boundary conditions. This result also holds in higher spatial dimensions, as one can see by analogy to the electric potential in a charge free space. The mass-redistribution potential plays a role analogous to the chemical potential in Model B dynamics Hohenberg and Halperin (1977). However, it does not follow from a free energy density. Instead, it is determined by the local concentrations via Eq. (11), and its spatial gradients represent the local imbalance of diffusive fluxes. Finally, note that the equation for the mass-redistribution dynamics Eq. (12) is not closed. Later, in Sec. IV, we will introduce an approximate “closure relation” for Eq. (12).

The above analysis can be generalized to systems with $N$ components whose total mass is conserved (describing, for instance, a single protein species with $N$ conformational states). The mass-redistribution potential is the sum of the concentrations weighted by their respective diffusion constants. Respectively, the flux-balance subspace in the $N$-dimensional concentration phase space is a $N-1$ dimensional hyperplane orthogonal to the vector of diffusion constants $(D_{1},D_{2},\ldots,D_{N})$.

III.3 Stationary patterns are “scaffolded” by

local equilibria

Whenever the diffusion constants $D_{m}$ and $D_{c}$ are unequal, the flux-balance subspace cannot coincide with any reactive phase space (which has slope $-1$). Hence, a non-uniform total density profile $\widetilde{n}(x):=\widetilde{m}(x)+\widetilde{c}\kern 0.74997pt(x)$ is innate to any stationary pattern (non-uniform $\widetilde{m}(x)$) whenever $D_{m}\neq D_{c}$.

As we will see next, this non-uniform total density profile is key understand the relationship between the stationary pattern in real space and the reactive nullcline in the phase plane.

Consider the system as being spatially dissected into notional local compartments Halatek and Frey (2018). Within each such compartment, local reaction kinetics induce a reactive flow $f(m,c)$ that lies in the local phase space $\{(m,c):m+c=\widetilde{n}(x)\}$ which is determined by the respective local mass $\widetilde{n}(x)$. We define the local equilibria

[TABLE]

analogously to Eq. (4), where we emphasize that the total density $n(x,t)$ is a function of position $x$ and time $t$ here. The local equilibria are geometrically determined by intersection points of the local phase space with the reactive nullcline. Together with their linear stability the local equilibria serve as proxies for the local reactive flow in each notional compartment (as in the well-mixed system discussed in Sec. III.1; cf. Eq. (5)). Thus, by thinking about a system as dissected into small compartments coupled by diffusion, we can carry over the phase-space structure of the local reaction kinetics to the spatially extended system (Fig. 2b).

What is the relationship between the local equilibria $\mathbf{u}^{*}(\widetilde{n}(x))$ and the stationary pattern $\widetilde{\mathbf{u}}(x)$? To gain some intuition, suppose the compartments were isolated from each other, i.e. diffusive coupling between them were shut off. If we choose the compartments small enough to be well-mixed, then the concentrations $m$ and $c$ in each of them will simply relax to the local equilibrium (black disks in Fig. 2) determined by the total density $\widetilde{n}(x)$ that varies from compartment to compartment. In that sense, the local equilibria act as a scaffold to which the pattern is “tied” by local reactive flows. Because total density must be conserved individually in each of the (now uncoupled) compartments the approach (red arrows) to the local equilibria is confined to the reactive phase space (gray lines) given by the total density in the compartment.

Let us now consider diffusive coupling between these compartments. In essence, diffusion acts to remove spatial gradients (as indicated by the small blue arrows along the FBS in Fig. 2b) and is counteracted by reactive flows towards the local equilibria (indicated by the red arrows from the FBS to the local equilibria in Fig. 2a,b). How does this competition play out in detail? In steady state, the net diffusive flux in and out of the compartment is balanced by the deviation from local (reactive) equilibrium (it is instructive to compare Eq. (4) for reactive equilibria and Eqs. (8) for stationary pattern). If the gradient does not change across the compartment, such that the flux into and out of the compartment are identical, the net diffusive flux vanishes: $\partial_{x}[D_{m}\partial_{x}\widetilde{m}(x)]=0$. (Thanks to the flux-balance condition Eq. (9), the same holds automatically for $\widetilde{c}\kern 0.74997pt(x)$.) In turn, the stationary pattern pattern must coincide with the local equilibria: $\widetilde{\mathbf{u}}(x)=\mathbf{u}^{*}(\widetilde{n}(x))$ when the gradient does not change across a compartment. This holds exactly at inflection points of the pattern. For plateaus, the gradient is small, $\partial_{x}\widetilde{\mathbf{u}}(x)\approx 0$, in a spatially extended region, and so is the local net flux, that is, $\partial_{x}^{2}\widetilde{\mathbf{u}}(x)\approx 0$. Hence, for plateaus we have $f(\widetilde{\mathbf{u}}(x))\approx 0$, such that the pattern can be locally approximated by the respective local equilibria, that is, $\widetilde{\mathbf{u}}(x)\approx\mathbf{u}^{*}(\widetilde{n}(x))$ in plateau regions.

III.4 The flux-balance construction on the reactive nullcline

Combining these insights with the fact that the stationary pattern must be embedded in a flux-balance subspace, we can identify plateaus and the inflection point as intersection points of a flux-balance subspace and the nullcline (FBS–NC intersections). As illustrated in Fig. 2, these ‘landmark points’ in phase space enable us to graphically construct the spatial patterns in real space as two plateaus connected by an interface (flux-balance construction).

Near the interface the densities $\widetilde{\mathbf{u}}(x)$ of a stationary pattern will deviate from the corresponding local equilibria. The ensuing reactive flows (red arrows) left and right of the inflection point are of opposite sign and correspond to attachment and detachment zones for protein patterns (see Fig. 2a) Halatek et al. (2018). Linearizing the phase-space flow around the landmark points will later enable us to further quantify the spatial profile of stationary patterns, i.e. to determine the relevant length scales.

III.5 Turnover balance determines $\eta_{0}$

Integrating one of the stationarity conditions, Eq. (14a), over the whole spatial domain $[0,L]$ yields that in steady state the total reactive turnover must vanish

[TABLE]

This total turnover balance determines the position $\eta_{0}$ of the flux-balance subspace in steady state. A mathematically more convenient form of turnover balance is obtained by multiplying Eq. (8a) with $\partial_{x}\widetilde{m}(x)$ before integrating:

[TABLE]

In this form, it becomes evident that the total turnover balance does not depend on the full density profile $\widetilde{m}(x)$, but only on the densities at the boundaries, $\widetilde{m}(0)$ and $\widetilde{m}(L)$. Total turnover balance, Eq. (17), together with the stationarity condition for $\widetilde{m}(x)$, Eq. (14a), fully determine the stationary patterns.

Geometrically, total turnover balance can be interpreted as a kind of (approximate) Maxwell construction in the $(m,c)$-phase plane (balance of areas shaded in red in Fig. 2). This requires the following approximations. First, we linearize the reactive flow around the reactive nullcline (cf. Eq. (5)):

[TABLE]

where $\widetilde{n}(m):=\eta_{0}+(1-D_{m}/D_{c})m$ because the pattern is embedded in the flux-balance subspace, cf. Eq. (10). The expression in the square brackets of Eq. (18) is simply the distance of the reactive nullcline from the flux-balance subspace measured along the respective local phase space. Further, suppose for the moment that the local eigenvalue $\sigma_{\text{loc}}(n)$ is approximately constant in the range of total densities attained by the pattern. Turnover balance, Eq. (17), is then represented by a balance of the areas between nullcline and flux-balance subspace on either side of the inflection point (see areas shaded in red (light gray) in Fig. 2b):

[TABLE]

In the characterization of pattern profiles in Sec. V, we will use that for a spatial domain size $L$ much larger than the interface width, one can approximate the plateau concentrations by FBS-NC intersections: $\widetilde{m}(0)\approx m_{-}(\eta_{0})$ and $\widetilde{m}(L)\approx m_{+}(\eta_{0})$. In this case, Eq. (17) is closed and can be solved for $\eta_{0}$, either numerically or geometrically using the approximate ‘Maxwell construction,’ Eq. (19).

Multiplying the stationarity condition, Eq. (14a), by $\partial_{x}\widetilde{m}(x)$ (as we did to obtain Eq. (17) for total turnover balance) and integrating over the spatial subinterval $[0,x_{0}]$), one obtains a relation that depends only on $\eta_{0}$ and the boundary concentrations $\widetilde{m}(0)$ and $\widetilde{m}(L)$:

[TABLE]

These equations state that the net turnover on either side of the pattern inflection point $x_{0}$, has to be balanced by the net diffusive flux across that point as illustrated by the blue arrow in Fig. 2a. Because the reactive flow changes sign at the inflection point, the reactive turnover (integrated flow) is extremal there and determines the maximal slope $\widetilde{m}^{\prime}(x_{0})$ of the pattern profile. Depending on how the reactive turnover saturates on either side of the inflection point, the system exhibits, as we will learn in Sec. V.3, three distinct characteristic elementary pattern types, classified by the shape of the concentration profile $\widetilde{m}(x)$: mesas, peaks/troughs and nearly harmonic (or ‘weakly nonlinear’) patterns.

III.6 Summary of geometric structures in phase space

Let us pause for a moment and briefly summarize our findings so far. We have established three major geometric structures in $(m,c)$-phase space: First, the reactive nullcline, $f(m,c)=0$, along which the local reaction kinetics are balanced; second, the local phase spaces, $m+c=\widetilde{n}(x)$, determined by the local total densities $\widetilde{n}(x)$—local equilibria $\mathbf{u}^{*}(\widetilde{n}(x))$ are intersections of the reactive nullcline and the local phase spaces; third, the family of flux-balance subspaces, within which diffusive flows in membrane and cytosol balance each other. The position of the flux-balance subspace, $\eta_{0}$, of a stationary pattern is determined by total turnover balance, Eq. (17), which represents a balance of reactive processes.

This geometric picture underlies the key results we present in the remainder of the paper. Up to now, we only discussed the embedding of the pattern in the $(m,c)$-phase plane. To study the possible pattern profile shapes $\widetilde{m}(x)$ in real space, we need to understand the dynamic process of pattern formation, in particular the factors determining the interface region. As we will see below (Sec. V), the interface of a pattern is inherently connected to lateral instability. We will therefore first analyze lateral instability and the dynamic process of pattern formation in the following section. With these tools at hand, we will then be able to characterize the distinct pattern types exhibited by 2C-MCRD systems.

IV Lateral instability

How can the geometric structures introduced in the previous section help us to understand the physical process of pattern formation? Previous research Halatek and Frey (2018) suggests that the total densities are the essential degrees of freedom and their redistribution is the key dynamic process. Building on this insight, we systematically connect the geometric structures established above (Sec. III) to the lateral instability, i.e. instability against spatially inhomogeneous perturbations, of a homogeneous steady state.

IV.1 Mass-redistribution instability

Consider the dynamics of the local total density $n(x,t)=c(x,t)+m(x,t)$. Because the kinetics conserve local total density, the time evolution of $n(x,t)$ is driven only by diffusion due to spatial gradients in the concentrations $c(x,t)$ and $m(x,t)$:

[TABLE]

As a result of mass redistribution, the local equilibrium concentrations $\mathbf{u}^{*}(n(x,t))$ change. In turn, these locally shifted equilibrium concentrations induce changes in the local reactive flows and thereby result in altered spatial gradients in $\mathbf{u}(x,t)$. This intricate coupling between redistribution of total mass, reactive flows, and diffusive flows drives pattern formation.

Qualitatively, this coupling between reactive and diffusive flow can be understood by observing that the dynamics depends mainly on the direction in which the local equilibria shift due to increasing or decreasing local total density. Let us therefore posit that the relevant diffusive gradients can be (qualitatively) estimated by replacing the local concentrations by the (locally stable) local equilibrium

[TABLE]

such that the local mass $n(x,t)$ is the only remaining degree of freedom:

[TABLE]

We term this the local quasi-steady state approximation. Note that this approximation becomes exact in the long wavelength limit where diffusive redistribution is much slower than chemical relaxation; see Sec. IV.2 and Appendices C and D for a detailed discussion. Applying the chain rule once, we can rewrite the mass redistribution dynamics as

[TABLE]

which is simply a diffusion equation for the total density $n(x,t)$. For locally stable equilibria, the effective diffusion constant will become negative (which entails anti-diffusion) if

[TABLE]

where $s_{\text{nc}}(n)=\partial_{m}c^{*}(m)|_{n}$ is the slope of the reactive nullcline $c^{*}(m)$ (cf. Eq. (6) in Sec. III.1; note that local stability ensures $\partial_{n}m^{*}>0$ when $f_{c}>0$, for $f_{c}<0$ the inequality Eq. (25) is reversed). Hence, starting from a homogeneous steady state $\mathbf{u}^{*}(\bar{n})$, a lateral instability due to effective anti-diffusion takes place if (and only if)

[TABLE]

This condition for lateral instability has a simple geometric interpretation in the $(m,c)$-phase plane: A spatially homogeneous steady state with total density $\bar{n}$ is laterally unstable if the slope of the nullcline is steeper than the slope of the flux-balance subspace. We term the mechanism a mass-redistribution instability to emphasize the underlying physical process and to contrast this mechanism to the “activator–inhibitor mechanism” (see Sec. VIII.2.3 in the Discussion). Importantly, the mass-redistribution instability in reaction–diffusion systems is a Turing instability Turing (1952), in the sense that it is a diffusion-driven instability of a system that is stable in a well-mixed situation (i.e. stable against spatially homogeneous perturbations) Note:Turing . The bifurcation where a homogeneous steady state becomes laterally unstable, i.e. Turing unstable, will be referred to as a Turing bifurcation.

The mass-redistribution dynamics, Eq. (24), can be rewritten most compactly using the mass-redistribution potential, $\eta$, (cf. Eq. (12)):

[TABLE]

This implies that a mass-redistribution instability occurs if an increase in total density entails a decrease of the mass-redistribution potential (i.e. $\partial_{n}\eta^{*}|_{\bar{n}}<0$).

Importantly the instability condition, Eq. (25), can be related directly to an underlying physical mechanism. For illustration purposes let us disregard membrane diffusion ($D_{m}=0$). Following a small modulation $\delta n$ of the mass on a large length scale, the local reactive equilibrium within each compartment shifts (1 in Fig. 3a). The instability criterion, Eq. (25), requires the slope of the reactive nullcline to be negative. In this case, the equilibrium shifts to lower cytosolic concentration $c^{*}(n)$ as total density $n$ is increased. In other words, in regions with a higher total density, there will be a reactive flow onto the membrane (red arrows) as the shifted local equilibrium is approached—thus creating cytosolic sinks (2 in Fig. 3a). Conversely, the regions with lower total density become cytosolic sources. The ensuing cytosolic gradient leads to diffusive mass-redistribution (3), resulting in a further shift of the local equilibria (4), thus sustaining and amplifying the diffusive flux—the cycle feeds itself. Taken together, this shows that the mass-redistribution instability is a self-amplifying mass redistribution cascade. In contrast, when the cytosolic equilibrium density rises due to an increase in total density (i.e. for positive nullcline slope $\partial_{n}c^{*}|_{\bar{n}}>0$) the compartment with more total density will become cytosolic source inducing mass redistribution that brings the system back to a homogenous state.

IV.2 Diffusion- and reaction-limited regimes

On sufficiently large length scales, diffusive relaxation (transfer of mass) is slow compared to chemical relaxation $D_{c}q^{2}\ll|\sigma_{\text{loc}}|$, such that the local quasi-steady state approximation Eq. (22) becomes exact—the concentrations are slaved to the local equilibria. This is the diffusion-limited regime: the growth rate of the lateral instability is limited by cytosolic redistribution via diffusion ($\sigma_{\text{lat}}\approx-\partial_{n}c^{*}\cdot D_{c}q^{2}$). In contrast, if cytosolic diffusion is much faster than chemical relaxation ($D_{c}q^{2}\gg|\sigma_{\text{loc}}|$), the lateral instability is limited by the rate at which the shifting equilibria are approached ($\sigma_{\text{lat}}\approx\partial_{n}c^{*}\,\cdot\,\sigma_{\text{loc}}=f_{m}$). This is the diffusion-limited regime. Importantly, the concept of shifting local equilibria still informs about the presence of the lateral instability in this regime. But it does no longer yield the growth rate quantitatively. A more detailed discussion of the local quasi-steady state approximation is provided in Appendix D.

The principle of shifting local equilibria provides insight into the spatial dynamics of systems with more than two components: In a five-component MCRD model for in vitro Min patterns, the concept of scaffolding allowed to predict the transition to chaos (qualitative change of the local attractors from stable fixed point to limit cycle) Halatek and Frey (2018). Notably, in this system, the onset of lateral instability is not a long wavelength instability but takes place for a band of unstable modes bounded away from $q=0$, corresponding to ‘type I’ instability in the Cross–Hohenberg classification scheme. Thus, the principle of shifting local equilibria is not restricted to systems with a long wavelength (‘type I’) instability.

IV.3 The marginal mode $q_{\mathrm{max}}$ reveals the role of membrane diffusion

Let us compare the instability criterion Eq. (25) to ‘classical’ linear stability analysis of Eq. (1) around the homogenous steady state (see Appendix C). There one obtains the dispersion relation $\sigma(q)$ for the growth rate $\sigma$ of a mode with wavenumber $q$ (see Fig. 3b). It exhibits a band of unstable modes, $\sigma(q)>0$ for $0<q<q_{\mathrm{max}}$, with

[TABLE]

if and only if $f_{m}/D_{m}>f_{c}/D_{c}$, i.e. exactly when the slope criterion, Eq. (25), is fulfilled. Equation (28) can be rewritten as $q_{\mathrm{max}}^{2}=f_{c}/D_{m}(-s_{\text{nc}}-D_{m}/D_{c})$, which shows why the slope criterion Eq. (25) is reversed for $f_{c}<0$.

The instability condition Eq. (25) and the expression for the edge of the edge of the band of unstable modes Eq. (28) inform about the role of membrane diffusion as counteracting the cytosolic mass-redistribution that drives the instability. This is because the membrane gradient will always be opposed to the cytosolic gradient whenever the nullcline slope is negative $(-1<\partial_{n}c^{*}<0)$ (because $\delta m=\delta n-\delta c$ and $\delta c^{*}=\delta n\,\partial_{n}c^{*}$).

The expression for $q_{\mathrm{max}}$ can be found quite easily by utilizing phase-space geometry. We start from the observation that $q_{\mathrm{max}}$ is a non-oscillatory marginal mode; it cannot be oscillatory for a locally stable fixed point, $\sigma_{\text{loc}}<0$, as shown in Appendix C. Hence, the eigenvalue $\sigma(q_{\mathrm{max}})=0$, so the mode $\sim\cos(q_{\mathrm{max}}x)$ must fulfill the steady-state condition, Eq. (8), and the corresponding eigenvector must point along a flux-balance subspace $\propto(1,-D_{m}/D_{c})^{\text{T}}$ in phase space. The steady state condition in flux-balance subspace is given by Eq. (14a), which, in linearization around a homogeneous steady state reads

[TABLE]

This equation is solved by the mode $\delta m(x)\propto\cos(q_{\mathrm{max}}x)$ with $q_{\mathrm{max}}^{2}=f_{m}/D_{m}-f_{c}/D_{c}$ (cf. Eq. (28)). To conclude, the two effects of membrane diffusion are interlinked in the expression for $q_{\mathrm{max}}$: (i) The condition $q_{\mathrm{max}}=0$ determines the critical NC-slope ($s_{\text{nc}}^{\text{crit}}=-D_{m}/D_{c}$) for the (long wavelength) onset of lateral instability. (ii) In the laterally unstable regime, $q_{\mathrm{max}}$ determines the smallest unstable length scale $\ell=q_{\mathrm{max}}^{-1}$. In the limit of large $D_{c}$, this length scale is given by $\ell^{2}\approx D_{m}/f_{m}$, i.e. a balance of membrane diffusion and reactive flows. In the next section, it will be shown that the marginal mode $q_{\mathrm{max}}$ at the pattern inflection point determines (to leading order) the interface width of a stationary pattern.

V Characterization of stationary patterns

With an intuitive picture of the principles underlying pattern formation in 2C-MCRD systems in hand, we now return to the spatially continuous system. We first study the characteristic types of stationary patterns exhibited by 2C-MCRD systems, focusing on elementary stationary patterns with monotonic concentration profiles on a domain with no-flux boundaries. More complex, non-monotonic stationary patterns (also in domains with periodic boundary conditions) can always be dissected into such elementary patterns at their extrema (recall that due to the diffusive flux-balance condition, Eq. (9), extrema in $\widetilde{m}(x)$ and $\widetilde{c}\kern 0.74997pt(x)$ must coincide). Previous studies have observed that 2C-MCRD systems typically exhibit coarsening Kang et al. (2007); Ishihara et al. (2007); Otsuji et al. (2007); Chiou et al. (2018). In a follow-up work building on the concepts presented here, we show that coarsening is indeed generic in all 2C-MCRD systems, independently of the specific form of the reaction kinetics Brauns et al. (2020a).

V.1 Interface width

The width of the interfacial region, $\ell_{\text{int}}$, is the only intrinsic length scale of the elementary patterns. Recall that the pattern inflection point, which defines the position of the interface region, is in local reactive equilibrium—geometrically determined by an FBS–NC intersection $(m_{0},c_{0})$ (cf. Fig. 2 in Sec. III.3). The interface is maintained by a balance of diffusion and the reactive flow in the vicinity of the inflection point. Therefore,the interface $\widetilde{m}(x-x_{0})\approx m_{0}+\delta\widetilde{m}(x)$ is to leading order determined by linearizing the steady-state equation, Eq. (14a), around the inflection point,

[TABLE]

where $n_{0}=m_{0}+c_{0}$, and we used the flux-balance subspace constraint, Eq. (10), to substitute the cytosol concentration $\delta\widetilde{c}\kern 0.74997pt(x)$. Equation (30) exactly resembles the equation that determines the marginal mode $\sin(q_{\mathrm{max}}x)$ in the dispersion relation (right-hand edge of the band of unstable modes). Hence, the interface length scale is determined by the marginal mode of the dispersion relation at the inflection point:

[TABLE]

The interface shape is approximated by the corresponding eigenfunction \delta\widetilde{m}(x)\propto\sin\big{(}q_{\mathrm{max}}(n_{0})x\big{)}.

Let us pause for a moment to look at the interface region from the perspective of mass redistribution: the total density $n_{0}$ at the inflection point is such that the corresponding reactive equilibrium is laterally unstable, because the nullcline slope is steeper than the FBS-slope there; see Fig. 4 below, where we elaborate on this point in terms of spatial regions. From the spectrum of modes, only the marginally stable one, $q_{\mathrm{max}}(n_{0})$, fulfills the (linearized) stationarity condition. Thus, intuitively it must be the $q_{\mathrm{max}}$-mode that determines the interface length scale. Importantly, because the pattern is formed by mass redistribution, the total density $n_{0}$ at the inflection point does not coincide in general with the average total density $\bar{n}$. The interface width is determined $q_{\mathrm{max}}(n_{0})$, not by $q_{\mathrm{max}}(\bar{n})$. This also implies that the interface width depends on the FBS-position $\eta_{0}$ because the inflection point $(m_{0},c_{0})$, and hence $n_{0}=m_{0}+c_{0}$, is determined geometrically as FBS-NC intersection point. We explicitly denote the interface width by $\ell_{\text{int}}(\eta_{0})$ when we use this relationship in the following.

Finally, to approximate the stationary concentration profile of the interface, we use that its maximal slope $\widetilde{m}^{\prime}(x_{0})$ is attained at the pattern inflection point $x_{0}$ and can be calculated by flux-turnover balance (20). Together with the harmonic mode $\delta\widetilde{m}(x)\propto\sin(\pi x/\ell_{\text{int}})$ obtained by linearizing phase space flow, we find

[TABLE]

in the vicinity of the inflection point. To go beyond this leading order approximations, one can perform a perturbative expansion of $\tilde{f}(m_{0}+\delta m,\eta_{0})$ and $\widetilde{m}(x_{0}+\delta x)$ in Eq. (14a) around the pattern inflection point $(m_{0},\eta_{0})$. The solution of this expansion can then be matched to the plateaus to obtain an approximation of the interface profile shape. Linearization around the plateaus yields exponential decay towards the plateaus (“exponential tails”) ${\sim}\exp(-x/\ell_{\pm})$, where the decay lengths are given by $\ell_{\pm}^{2}=D_{m}/\tilde{f}_{m}(n_{\pm})$.

V.2 Regions generalize the concept of local compartments

Not only the interfaces ($n_{0}$) but also the plateaus ($n_{\pm}$) of patterns correspond to FBS-NC intersection points in phase space (Fig. 4). However, in contrast to the interface, the plateaus lie on laterally stable sections of the nullcline where $s_{\text{nc}}(n_{\pm})>-D_{m}/D_{c}$; recall the slope criterion for lateral instability Eq. (25). Put more precisely, the pattern profile becomes flat in the vicinity of $n_{\pm}$ because of regional lateral stability Note:Splitting . Ivanov et al. (2010); Crampin et al. (2002)

Thus, the FBS-NC intersection points act as “landmark points” that enable us to (notionally) dissect the pattern profile into spatial regions (plateaus and interface) in such a way that these spatial regions can be associated with regional phase spaces. The (linearized) properties of the reaction–diffusion dynamics—encoded in the regional dispersion relations (Fig. 4d)—in the vicinity of these landmark points can be used to determine the regional properties in real space. Within each of these regions, we can ask what would happen there if we were to isolate it from the rest of the system, akin to the question we asked in the context of local equilibria. Just as the local equilibria scaffold the interface, these regional “attractors” serve as scaffolds for the global pattern. The regional properties depend on the average regional mass which is redistributed between regions by diffusion, and the properties of the full pattern can be pieced together by (characteristically distinct) isolated regions (plateaus and interfaces).

Taken together, the nonlinear kinetics is encoded by the nullcline shape. The internal properties of the spatial regions are determined by regionally linear properties of phase space flow, encoded in the regional dispersion relations. We will therefore refer to this as the method of regional phase-spaces and regional attractors. This method bridges the gap between the linear and the highly nonlinear regime.

Finally, let us note that based on the region decomposition (cf. Fig. 4), the interface position—determining the global spatial structure—can be pictured as a collective degree of freedom. A conceptually similar, but technically more involved approach to study interfaces (also called ‘kinks’ or ‘internal layers’) and their dynamics is singular perturbation theory (specifically matched asymptotic expansion) where one uses an asymptotic separation of spatial scales, see e.g. Ref. Ward (2006) and references therein. Such methods also facilitate a phase-space geometric analysis Jones (1994).

V.3 Pattern classification

Employing the concept of regions we will now turn to the classification of patterns. We distinguish two generic pattern types: mesas and peaks. Mesa patterns are composed of plateaus (low density and high density) connected by an interface (Fig. 5a, c and d), while the term peak refers to an interface concatenated to a plateau only on the low density site (Fig. 5b) Note:Trough . For small systems, close to onset, there is an additional pattern type comprising only an interface that spans the whole system; see Sec. VII.4.

What are the conditions for the formation of a peak pattern versus the formation of a mesa pattern? A mesa pattern requires two plateau regions, each characterized by an FBS-NC intersection point, one at low density and one at high density. The low density plateau is generically present because the densities must be positive and thus are bounded from below. In contrast, the position of the FBS-NC intersection point at high density depends sensitively on the shape of the nullcline and the slope of the diffusive flux-balance subspace $-D_{m}/D_{c}$. Let us first consider the case of fast cytosol diffusion $D_{c}\gg D_{m}$. For an N-shaped nullcline, i.e. one that has an “upwards-pointing” tail (see Fig. 5a), the flux-balance construction presented in Sec. III.3 yields a mesa pattern. The situation is different for a “$\Lambda$-shaped” nullcline that has an asymptotically flat tail for large $m$ (e.g. approaching $c^{*}(m)\to 0$ for $m\to\infty$); see Fig. 5b. In that case the third FBS-NC intersection point generically is far away from the first two; in Fig. 5b it lies out of frame. The requirement of total turnover balance (approximated by a balance of the areas shaded in light red in the Fig. 5) limits the maximum membrane concentration $\hat{m}$, such that there is no high density plateau and the pattern assumes a peak profile instead (Fig. 5b bottom). In the detailed analysis of peak patterns below, we will show that the FBS-position $\eta_{0}$, and thus the peak amplitude, is determined by the total mass in the system

Let us now consider what happens if the FBS is made steeper by lowering the ratio of diffusion constants $D_{c}/D_{m}$. As the FBS becomes steeper, the third FBS-NC intersection point moves towards lower membrane concentration. Eventually, this limits the total turnover on the right hand side of the inflection point, such that the peak amplitude saturates in a plateau, i.e. a mesa pattern forms (Fig. 5d). In the case of an N-shaped nullcline, lowering $D_{c}$ reduces the concentration difference between the two plateaus, because the FBS-NC intersection points move closer together (Fig. 5c).

Taken together, the phase-plane analysis reveals how the interplay of nonlinear reactions (encoded in the nullcline shape) and diffusion (encoded in the FBS-slope) determine the pattern type and pattern amplitude.

V.3.1 Mesa patterns

To characterize mesa patterns in the limit $L\gg\ell_{\text{int}}$, we first determine the FBS-position, $\eta_{0}$, using total turnover balance, Eq. (17) (cf. Sec. III.3). Since the plateaus are scaffolded by (laterally stable) local equilibria we can approximate the boundary concentrations

[TABLE]

where the plateau scaffolds $m_{\pm}(\eta_{0})$ are geometrically determined in phase space as intersection points $m_{\pm}(\eta_{0})$ of FBS and NC:

[TABLE]

With the approximation, Eq. (33), the total reactive turnover balance condition, Eq. (17), becomes:

[TABLE]

where $m_{\pm}^{\infty}=m_{\pm}(\eta_{0}^{\infty})$ and $\eta_{0}^{\infty}$ denotes the FBS-position in the large system size limit. Equation (35) is closed and can be solved for $\eta_{0}^{\infty}$. Once one has determined $\eta_{0}^{\infty}$, the interface width $\ell_{\text{int}}(\eta_{0}^{\infty})$ can be estimated with Eq. (31).

This total turnover balance condition implicitly determines the FBS-offset $\eta_{0}=\eta_{0}^{\infty}$. Note that equation Eq. (35), and hence $\eta_{0}^{\infty}$ depends only on the function $f$ and the ratio of the diffusion constants, but not on the average mass $\bar{n}$ or the domain size $L$ in the limit $L\gg\ell_{\text{int}}$. Instead, the average total density $\bar{n}$ determines the position $x_{0}$ of the interface. Again assuming an interface much narrower than the domain size, the contribution of the interface region can be neglected, $L\bar{n}\approx n_{-}^{\infty}x_{0}+n_{+}^{\infty}(L-x_{0})$, which yields

[TABLE]

where $n_{+}^{\infty}$ and $n_{-}^{\infty}$ are the average total densities in the plateau regions:

[TABLE]

This shows that the amplitude of mesa patterns is geometrically determined by the reactive nullcline alone and does not sensitively depend on average mass $\bar{n}$ or system size $L\gg\ell_{\text{int}}$. Moreover, far away from the critical point $D_{c}^{\mathrm{min}}$, cf. Sec. VII.1 the mesa-pattern amplitude becomes independent of the ratio of the diffusion constants. Adding mass to a mesa pattern shifts the interface position $x_{0}$ as the additional mass is redistributed between the two plateau regions.

Notably, a geometric argument shows that mesa patterns are the generic pattern for $L\to\infty$ when the ratio of the diffusion constants is nonzero $D_{m}/D_{c}>0$, and $m\geq 0$, $c\geq 0$ (as must be the case for concentrations): The FBS intersects the $m$-axis ($c=0$) at $(D_{c}/D_{m})\eta_{0}$, and hence must intersect the nullcline at some finite value $m<(D_{c}/D_{m})\eta_{0}$. For $L\to\infty$ keeping the average mass $\bar{n}$ constant, the pattern profile will eventually reach this third FBS–NC intersection point, and thus become a mesa pattern. Next, we will discuss the conditions under which peak/trough patterns occur.

The approximation Eq. (33) for the plateau densities, and in turn also Eq. (36) for the interface position, will break down when the distance of the interface to one of the system boundaries becomes smaller than the interface width $\ell_{\text{int}}(\eta_{0}^{\infty})$. Then the stationary pattern no longer exhibits a plateau on that side and instead becomes a plateau–interface pattern, forming either a peak when $\bar{n}$ is close to $n_{-}^{\infty}$, or a trough (‘anti-peak’) when $\bar{n}$ is close to $n_{+}^{\infty}$. An estimate for these transition from mesa to peak/trough patterns can be obtained based on the approximated interface position, Eq. (36):

[TABLE]

V.3.2 Peak patterns

Let us now study these peak/trough patterns. Their defining characteristic is that a plateau, corresponding to laterally stable FBS-NC intersection point, forms only only on one side of the interface. Correspondingly, the reactive turnover saturates on the side where the plateau forms, while it depends on the variable pattern amplitude on the other side. For specificity, we focus on peak patterns here. As explained above, such a peak pattern forms when the nullcline is $\Lambda$-shaped, flux-balance subspace is very shallow ($D_{c}\gg D_{m}$); see Fig. 5b.

Suppose for a moment that we can freely choose the FBS-position, $\eta_{0}$, and that the total mass is not fixed. Given some $\eta_{0}$, the FBS-NC intersection point (1) determines the low-density plateau at the foot of the interface and thus the total turnover on this side, corresponding to the enclosed area in the interval between (1) and (2). This turnover must be balanced by an equal and opposite turnover on the right. Using again the (approximate) correspondence to the enclosed phase-plane area, it becomes obvious that this balance of areas determines the peak amplitude $\hat{m}$ (3). Using that the interface profile can be approximated as $\propto\sin((x-x_{0})/\ell_{\text{int}})$, where $\ell_{\text{int}}$ is determined by the steady state equation linearized around the pattern inflection point (see Sec. V.1), we can roughly estimate the total mass in a peak as Note:RollingBallPeak

[TABLE]

The sinusoidal shape of the interface furthermore mandates that the inflection point $m_{0}$ is approximately half-way between the plateau $m_{-}$ and the maximum $\hat{m}$, such that we can eliminate $\hat{m}\approx m_{-}+2\,(m_{0}-m_{-})$ in Eq. (39). The remaining unknowns $m_{0}$ and $m_{-}$ are determined geometrically (FBS–NC intersections) as functions of $\eta_{0}$. Thus we obtain a relation for the average total density $\bar{n}(\eta_{0})\approx n_{-}(\eta_{0})+N_{\text{peak}}(\eta_{0})/L$ as a function of $\eta_{0}$. The inverse of this relation yields the FBS-position $\eta_{0}(\bar{n})$ as a function the control parameter $\bar{n}$. This estimate will hold until the peak density $\hat{m}$ reaches the third FBS-NC intersection point $m_{+}(\eta_{0})$, where a second plateau will start to form, such that the peak pattern transitions to a mesa pattern. In Appendix G we present the details of the peak approximation, and a comparison to numerical solutions.

Our estimate for the peak mass Eq. (39) and the resulting relation $\eta_{0}(\bar{n})$ show that, in contrast to mesa patterns, the amplitude of peak patterns sensitively depends on the total mass $N=L\,\bar{n}$ and the membrane diffusion constant (via $\ell_{\text{int}}^{2}\sim D_{m}/f_{m}$, cf. Eq. (31)). In addition, the position of the third FBS-NC intersection point $m_{+}(\eta_{0})$ that limits the maximum peak density, sensitively depends on the FBS slope $-D_{m}/D_{c}$. In the limit $D_{m}/D_{c}\to 0$, the third FBS-NC intersection point $m_{+}$ moves to infinity. Hence, in this limit, a system with an asymptotically flat nullcline tail never exhibits mesa patterns.

V.3.3 More general nullcline shapes

Here we considered two types of nullcline shapes—N- and $\Lambda$-shaped—that both have a single maximum in the $(m,c)$-phase plane, but differ in their ‘tail behavior’. Beyond these two prototypical nullcline shapes, more general nullcline shapes are possible. For instance, reaction kinetics of the attachment–detachment form, Eq. (2), with higher order nonlinearities (e.g. 5th-order polynomials) may exhibit nullclines with multiple maxima in the $(m,c)$-phase plane. For general reaction kinetics $f(m,c)$, more exotic shapes of the nullclines, e.g. with multiple disconnected branches, are possible. Our findings equally apply to such nullclines since our phase-space analysis is based on simple geometric properties such as slopes and intersection points with the FBS. Additional care is required if $f_{c}$ changes sign along the nullcline, since the slope criteria for local and lateral stability, Eqs. (7) and (25), are reversed for $f_{c}<0$. Conveniently, for reaction kinetics of the attachment–detachment form, Eq. (2), one has $f_{c}=a(m)$ which is generally positive for systems of biochemical origin.

V.4 Generic bifurcation structure under variation of the average mass $\bar{n}$

Now that we have classified the different types of stationary patterns exhibited by 2C-MCRD systems, we turn to study bifurcations where the patterns change structurally or in stability. The bifurcation parameter we study first is the average total density $\bar{n}$. This parameter does not affect the phase-space geometry (NC and FBS), which makes it particularly easy to study. Later, in Sec. VII, we generalize our findings to bifurcation parameters that change the phase-space geometry: diffusion constants change the FBS-slope, whereas kinetic rates affect the nullcline shape. For biological systems, the average total density $\bar{n}$ is a natural parameter as it can be tuned by up- or down-regulating the production of a protein.

Let us begin with the bifurcations where the homogenous steady state becomes laterally unstable. We already learned in Sec. IV that there is a band of unstable modes, $[0,q_{\mathrm{max}}]$, if the NC-slope $s_{\text{nc}}(\bar{n})$ is negative and steeper than the FBS-slope, $-D_{m}/D_{c}$ (cf. Eq. (25) and Fig. 3a). Hence, a band of unstable modes exists if $\bar{n}$ is in the range $(n_{\text{lat}}^{-},n_{\text{lat}}^{+})$, bounded by the points $n_{\text{lat}}^{\pm}$ where the flux-balance subspace is tangential to the reactive nullcline (dash-dotted green lines in Fig. 6a. (Note that a system of finite size $L$, is unstable if the longest wavelength mode lies in the band of unstable modes $\pi/L<q_{\mathrm{max}}(\bar{n})$, where $q_{\mathrm{max}}^{2}=\tilde{f}_{m}/D_{m}$, as defined in Eq. (28) and $\tilde{f}_{m}=f_{m}-f_{c}\,D_{m}/D_{c}=(-D_{m}/D_{c}-s_{\text{nc}})f_{c}$.)

What about the range where stationary patterns exist? The plateau scaffolds $m_{\pm}(\eta_{0})$ are geometrically determined by the reactive nullcline via the FBS–NC intersection points. The position $\eta_{0}$ of the flux-balance subspace generally depends on $\bar{n}$ and $L$ via total turnover balance, Eq. (17). However, in the large system size limit ($L\to\infty$), the FBS position $\eta_{0}^{\infty}$ is independent of $\bar{n}$ and $L$ (cf. Eq. (35)). For patterns to exist, the average total density $\bar{n}$ must lie in-between the plateau densities $n_{\pm}^{\infty}$; see Fig. 6a, cf. Eq. (36). Hence, in the limit $L\to\infty$, stationary patterns exist in the range $n_{-}^{\infty}<\bar{n}<n_{+}^{\infty}$.

Importantly the range of pattern existence generically extends beyond the range of lateral instability ($n_{\text{lat}}^{-}>n_{-}^{\infty}$ and $n_{\text{lat}}^{+}<n_{+}^{\infty}$) by geometric necessity for N-shaped nullclines; see Fig. 6a. This implies, that generically there are regions of multistability in parameter space, where stable stationary patterns exist, and the homogeneous steady state is stable (regions shaded in blue in Fig. 6).

To gain some intuition on the steady states in the multistable regimes, we performed numerical continuation (see Appendix F for details) of the stationary patterns for an example 2C-MCRD system using the attachment–detachment kinetics Eq. (54*) from Ref. Mori et al. (2008) which exhibit an N-shaped nullcline. Figure 6b shows the numerically obtained bifurcation structure where we plot the pattern amplitude $|\widetilde{m}(L)-\widetilde{m}(0)|$ against the bifurcation parameter $\bar{n}$. The star marks a typical stable mesa pattern (see inset in Fig. 6a) in the central region of the branch of stable patterns (solid blue line). As the plateaus are scaffolded by the FBS-NC intersections $m_{\pm}(\eta_{0}^{\infty})$, the pattern amplitude stays approximately constant ($m_{+}{-}m_{-}$, dotted blue line) across the whole range of $\bar{n}$ where patterns exist. Changing total average density simply shifts the interface position (cf. panels $(ii)$ and $(iii)$ in (c)). When the interface position is in the vicinity of a boundary, mesa pattern transitions to peak patterns (see. panels $(i)$ and $(iv)$ in (c)) as we learned in the previous section (Eq. (38) in Sec. V.3). The numerical continuation shows, that the peak/trough patterns are then annihilated in saddle-node bifurcations (SN), where the branch of stable patterns meets a branch of unstable patterns (dashed, red line). Due to the finite system size, the exact positions of the SN-bifurcation points are sightly offset from $n_{\pm}^{\infty}$ (by an amount ${\sim}\,\ell_{\text{int}}/L$). The branches of unstable patterns emerge from homogenous steady state in subcritical pitchfork bifurcations (P) at the Turing bifurcations ($n_{\text{lat}}^{\pm}$). (In a finite sized system, the onset of lateral instability is offset by an amount ${\sim}\,L^{2}$ from the geometrically defined points $n_{\text{lat}}^{\pm}$, because the system is unstable only if the longest wavelength mode lies within the band of unstable modes; see Sec. VII.4). In the multistable regions (shaded in blue), patterns can be triggered by a finite amplitude perturbation. The unstable patterns are “transition states” (or “critical nuclei”) that lie on the separatrix separating the basins of attraction of the stable patterns and the stable homogeneous steady state. The actual separatrix is a complicated object in the high-dimensional PDE phase space. In the next section (VI), we will show that a heuristic can be inferred from the nullcline shape to estimate the threshold for stimulus-induced pattern formation for a prototypical class of spatial perturbation profiles.

Because the unstable patterns are peak/trough patterns, they can be approximated by the ‘peak approximation’ introduced in Sec. V.3 (see also Appendix G). Thus the qualitative structure of the branches of unstable patterns is determined by $(m,c)$-phase-space geometry independently of the details of the reaction term $f(m,c)$, as long as the reactive nullcline $f(m,c)=0$ is N-shaped.

In summary, we conclude that the qualitative form of the bifurcation structure shown in Fig. 6b is determined by geometric relations in $(m,c)$-phase space. In particular, we find that 2C-MCRD systems generically have regions of multistability and that the onset of lateral instability is generically subcritical for large system size $L\gg\ell_{\text{int}}$. Further, this implies that such systems exhibit stimulus-induced pattern formation, and that there is hysteresis of stationary patterns when the total average density is varied.

VI Perturbation threshold for stimulus-induced pattern formation

Before we delve into the more technical analysis of bifurcation structures, we would like to discuss one more important aspect of pattern formation: stimulus-induced pattern formation, i.e. the ability to induce the transition from one stable attractor (homogeneous steady state) to another (stationary pattern) by a large enough perturbation (stimulus). (In the context of phase separation, this is called nucleation and growth). Stimulus-induced pattern formation is a particularly important aspect of 2C-MCRD systems, because, as we have shown above, these systems generically have regions of multistability. Furthermore, biologically it is often desirable to be able to form a pattern following an external or internal stimulus that exceeds a certain threshold (“nucleation threshold”) . As of yet, this threshold could only be determined numerically Jilkine and Edelstein-Keshet (2011). In the following, we will show how simple heuristic reasoning—based on regional lateral instability—yields a geometric criterion for the perturbation threshold in the $(m,c)$-phase plane.

As we have shown in the previous section, the hallmark of a stationary pattern is a laterally unstable region surrounding the pattern inflection point $x_{0}$ (even if the homogeneous state of the system is laterally stable). In the proposed framework, the phase-space dynamics are simply represented by the expansion of the system in the $(m,c)$-phase plane due to mass redistribution. Hence, to lead to a stationary pattern, a trajectory in the (high-dimensional) phase space of a partial differential equation (PDE) must enter and remain in a (linearly) laterally unstable region in the $(m,c)$-phase plane (shaded in green in Fig. 7). The laterally unstable region in $(m,c)$-space corresponds to a respective region in real space. If the homogenous state is laterally stable then a finite perturbation (stimulus) is required to create a laterally unstable region. Let us study a prototypical perturbation able to induce a laterally unstable region: a step function that represents moving a ‘block’ of protein mass (total density) from one end of the system to the other; for an illustration of the spatial perturbation and the resulting flows in phase space see Fig. 7. Generalization to other perturbations is straightforward and based on analogous arguments. Such perturbations can be created by various means of ‘active’ mass redistribution, e.g. active transport in the cell cortex, along microtubules, and hydrodynamic cytosolic flows; see for instance Ref. Goehring et al. (2011).

Following a (large amplitude) perturbation, there are two distinct processes that are triggered in phase space as shown in Fig. 7. On the one hand, in the laterally unstable region (green shaded area), a mass-redistribution instability will start to form a pattern, thus further amplifying the perturbation. On the other hand, because the perturbation shifts the regional reactive equilibria (black disks), there will be reactive flows (red arrows) in the regions that induce a cytosolic gradient which leads to mass redistribution between the regions by cytosolic diffusion (large blue arrows). If the cytosolic density in the laterally stable region is lower than in the laterally unstable one, the regional instability may not be sustained and the system returns to homogenous steady state (Fig. 7b). Conversely, if the cytosolic density is lower in the laterally unstable region than in the laterally stable region, then the cytosolic flow between the regions (blue arrow) will sustain the regional instability (Fig. 7c). Because the mass-redistribution instability creates a self-organized and self-sustaining cytosolic sink, the laterally unstable region can be self-sustained. The heuristic criterion for a (self-)sustained laterally unstable region is that the perturbation must cross the nullcline (see Fig. 7c). Then, the overall cytosolic concentration in the laterally unstable region is decreased by reactive flows (red arrows) such that cytosolic diffusion (blue arrow) between the regions will sustain the laterally unstable region.

In Appendix H we show that this simple criterion already provides a very good approximation for the threshold in comparison to full numerical simulation. We conclude that the reactive nullcline provides the key information for understanding pattern formation dynamics, in a similar way as for the characterization of stationary patterns (Sec. V and the analysis of the linear mass-redistribution instability (Sec. IV). Specifically it enables one to estimate the basins of attraction of the uniform steady state and the polarized pattern. We further learned that regional lateral instability underlies stimulus-induced pattern formation from laterally stable homogeneous steady states.

The threshold estimate provided here might help to understand this “nucleation” of patterns from laterally stable homogeneous steady states. The unstable peak/trough patterns (dashed red lines in Fig. 6 are part of the separatrix between to the basin of attraction of stable stationary patterns, and can be pictured as canonical critical nucleus Bates and Fife (1993). The peak approximation described in Sec. V.3 and compared to numerical continuation in Appendix G provides a simple estimate for this critical nucleus.

VII Complete bifurcation structure

Bifurcation diagrams of 2C-MCRD systems were previously studied for specific choices of the reaction kinetics $f(m,c)$ using numerical methods Mori et al. (2011). Furthermore, based on numerical studies of various models, it was hypothesized that there might be a general bifurcation scenario underlying cell polarity systems Trong et al. (2014). Here we use the insight gained on phase-space geometry to systematically build the complete general bifurcation structure of 2C-MCRD systems. Our findings generalize previous results and unify them in the context of phase-space geometry. For large system size, the bifurcation structures are fully determined geometrically. We illustrate the effect of finite system size using numerically computed bifurcation diagrams shown in Appendix F.

Above, we studied the bifurcation diagram of stationary patterns for the bifurcation parameter $\bar{n}$ in a system with monostable kinetics; see Sec. V.4 and Fig. 6 therein. Recall, that in large systems ($L\to\infty$), the bifurcation points in $\bar{n}$ can be found based on geometric reasoning in phase space: (i) Lateral instability is identified by a criterion on the nullcline slope: $s_{\text{nc}}(\bar{n})<-D_{m}/D_{c}$. Hence, the range of lateral instability is bounded by points $n_{\text{lat}}^{\pm}$ where the FBS is tangential to the NC: $s_{\text{nc}}(n_{\text{lat}}^{\pm})=-D_{m}/D_{c}$. (ii) FBS–NC intersection points $m_{\pm}(\eta_{0}^{\infty})$ provide the scaffold for the plateaus of mesa patterns, where the FBS-position, $\eta_{0}^{\infty}$, is determined by total turnover balance, Eq. (35). Mesa patterns exist as long as the average total density can be distributed between two plateaus $n_{\pm}^{\infty}$, i.e. in the range $n_{-}^{\infty}<\bar{n}<n_{+}^{\infty}$; recall that $n_{\pm}^{\infty}=n_{\pm}(\eta_{0}^{\infty})$ depend on the position $\eta_{0}^{\infty}$ and slope $-D_{m}/D_{c}$ of the FBS (cf. Eq. (37)).

Both of these geometric bifurcation criteria depend on the diffusion constants via the slope of the flux-balance subspace $-D_{m}/D_{c}$. We keep $D_{m}$ fixed—thus fixing the smallest characteristic length scale $\ell=\sqrt{D_{m}/f_{m}}$ where spatial structures can be maintained against membrane diffusion—and vary $D_{c}$ to rotate the FBS in $(m,c)$-phase space.

VII.1 Generic bifurcation structure of stationary patterns for monostable reaction kinetics

We construct the $(\bar{n},D_{c})$-bifurcation diagram by inferring $n_{\text{lat}}^{\pm}$ and $n_{\pm}^{\infty}$, as described above, as functions of $D_{c}$. Qualitatively, this can even be done manually with pen and paper in the spirit of a graphical construction (see e.g. Ref. Strogatz (2014)) based on the geometric criteria (i), (ii) above, as shown in Fig. 6. Figure 8a shows the qualitative structure obtained by this graphical construction. A quantitative construction of the bifurcation diagram can be performed with simple numerical implementation of the bifurcation criteria described above, e.g. in Mathematica (see SM File: “flux-balance-construction.nb”, and Fig. 20 for figures of quantitative bifurcation structures). As we will see in the following, the structure of the bifurcation diagram is qualitatively the same for all monostable, N-shaped nullclines, independently of the details (nonlinearities and kinetic rates) of the reaction term $f(m,c)$. The bifurcation diagram is qualitatively different when the nullcline has a segment of bistability (where $s_{\text{nc}}<-1$, cf. Fig. 1). We will analyze this case and in particular the role of bistability further below in Sec. VII.2.

As $D_{c}$ is decreased, the flux balance subspace becomes steeper, and thus the bifurcation points $n_{\text{lat}}^{\pm}$ and $n_{\pm}^{\infty}$ start to converge (see Fig. 8; cf. Fig. 6). They meet in the inflection point of the reactive nullcline, $n_{\text{inf}}$, where the nullcline slope, $s_{\text{nc}}^{\text{inf}}$, is extremal ($\partial_{n}s_{\text{nc}}|_{n_{\text{inf}}}=0$). The extremal nullcline slope at the nullcline inflection point determines the minimal cytosolic diffusion constant,

[TABLE]

above which there are three FBS–NC intersection points. When the ‘critical’ point $(n_{\text{inf}},D_{c}^{\text{min}})$ is traversed in $D_{\kern-0.52498ptc}\kern 0.74997pt$-direction, the FBS–NC intersections bifurcate in a (supercritical) pitchfork bifurcation; see Fig. 8b. Since the FBS–NC intersection points $m_{\pm}$ are the scaffolds for the plateaus (in short: plateau scaffolds; cf. Fig. 2), this bifurcation at the critical point $(n_{\text{inf}},D_{c}^{\text{min}})$ is a bifurcation of the scaffold itself. Importantly, the actual pattern is bounded by the plateau scaffolds. Thus, if there are no plateau scaffolds (i.e. only one FBS-NC intersection point), there cannot be stationary patterns. For $L\to\infty$, patterns emerge slaved to the plateau scaffold, such that the pattern bifurcation is supercritical at the nullcline inflection point ($\bar{n}=n_{\text{inf}}$). Away from the nullcline inflection point ($n\neq n_{\text{inf}}$), the lateral instability bifurcation is always subcritical for $L\to\infty$ because the range $(n_{-}^{\infty},n_{+}^{\infty})$ where patterns exist always exceeds the range $(n_{\text{lat}}^{-},n_{\text{lat}}^{+})$ of lateral instability, as we learned above in Sec. V.4 (cf. Fig. 6).

As we will see below in Sec. VII.4, for finite $L$, the bifurcation is supercritical in the vicinity of the nullcline inflection point. The transition from super- to subcriticality depends on a subtle interplay of diffusive and reactive flow together with geometric factors like nullcline curvature.

Interestingly, the regimes and their interrelation in the $(\bar{n},D_{c})$-bifurcation diagram, as shown in Fig. 8a, are phenomenologically similar to the phase diagram of (near equilibrium) phase separation kinetics for binary mixtures, described by Cahn–Hilliard equation Cates and Tjhung (2017). In a previous study using the amplitude equation formalism, Ref. Bergmann et al. (2018), a mapping from 2C-MCRD models to Model B has been found for the vicinity of the critical point, where the pattern emerges from the Turing bifurcation in a supercritical or weakly subcritical pitchfork bifurcation (see Sec. VII.4).

Strikingly, our geometric reasoning shows that the physics implied by the bifurcation diagram is the same as in phase separation kinetics (binodal and spinodal regimes) for all N-shaped nullclines, and far away from the critical point. We discuss this finding in Sec. VIII.3.

VII.2 Locally bistable kinetics

Changing the kinetic rates deforms the nullcline shape. When the nullcline slope becomes smaller than $-1$, a regime of locally bistable reaction kinetics emerges (cf. Fig. 1).

VII.2.1 Fronts in bistable media

To elucidate the role of bistability, let us first consider the case of equal diffusion constants $D_{c}=D_{m}=:D$. (Although this does not make sense in the intracellular context anymore, where typically $D_{c}>D_{m}$, we stick to the notation with concentrations $m$ and $c$.) Then mass redistribution decouples from the kinetics $\partial_{t}n=D\partial_{x}^{2}n$, i.e. the total density becomes uniform by diffusion (see the mass-redistribution dynamics, Eq. (12), and note that $\eta(x,t)=n(x,t)$ for equal diffusion constants). As a consequence, the system can be reduced to one component, for instance the membrane density

[TABLE]

where the local kinetics is bistable at every point in space (Fig. 9). This corresponds to a (classical) one-component model for bistable media which generically exhibits propagating fronts Mikhailov (1994). A standard calculation, commonly performed by ‘Newton mapping’ (briefly described in Sec. III.3) or by phase-space analysis (in $(m,\partial_{x}m)$-phase space), shows that the propagation velocity $v$ of a front is proportional to the imbalance in reactive turnover Mikhailov (1994): $v\propto\int_{m_{-}}^{m_{+}}\kern-1.99997pt\mathrm{d}m\,f(m,\bar{n}-m)$. Hence, a stationary front can only be realized by fine-tuning of parameters, e.g. the average total density $\bar{n}$, such that total turnover is balanced:

[TABLE]

The balanced (fine-tuned) case corresponds to the Allan–Cahn equation  Allen and Cahn (1972) (also called ‘model A’ dynamics Hohenberg and Halperin (1977)). (In a finite size system, there will be logarithmic coarsening, see Ref. Ward (2006) and references therein.)

With respect to the concept of local equilibria as scaffolds for patterns (cf. Sec. III.3), the bistable local equilibria (fixed points $m_{\pm}$) can be regarded as a static scaffold for front solutions; see Fig. 9. Because there is no mass-redistribution, the scaffold must remain static and can not adapt to balance the total reactive turnover. Instead, fine-tuning of parameters (e.g. $\bar{n}$), is required to obtain a balance of total turnover and thus a stationary front. In the $(\bar{n},D_{c})$-bifurcation diagram, Fig. 10a, the stationary bistable front with a static scaffold appears only at a singular point $(\bar{n},D_{c})=(n_{\text{stat}},D_{m})$.

What happens when the diffusion constants are unequal $D_{c}\neq D_{m}$? Then, mass will be redistributed, leading to shifting of the local equilibria that scaffold the pattern. As we know from our analysis for monostable reaction terms, this dynamic scaffold is able to self-balance the total reactive turnover—fine-tuning of $\bar{n}$ is no longer required to obtain stationary patterns. Interestingly, for a bistable reaction term, stationary patterns can be constructed both for $D_{m}>D_{c}$ and for $D_{m}<D_{c}$, as we discuss next; see Fig. 10b,c. To determine the stability of these patterns, we will examine below (after the description of the bifurcation diagram) how the scaffold self-balances via mass redistribution.

VII.2.2 Bifurcation diagram for locally bistable reaction kinetics

The bifurcation diagram, Fig. 10a, for the large system size limit ($L\to\infty$) is obtained using the same geometric criteria as for the case of locally monostable reaction kinetics (see Fig. 10b, cf. Fig. 6a). The presence of a bistable nullcline segment does not affect the feasibility of the geometric construction itself. However, the bifurcation diagram one obtains is qualitatively different from the monostable case, as we will see next. We discuss the regimes of stationary patterns (shaded in blue and red) first, before we analyze the regions of lateral instability (shaded in green and orange).

The region where stable stationary patterns exist (shaded in blue) is delimited by lines $n_{\pm}^{\infty}(D_{c})$. These lines converge in singular point $(n_{\text{stat}},D_{c}=D_{m})$, where a marginally stable front exists in a bistable medium without mass redistribution (Sec. VII.2.1). Along the entire line, $D_{c}=D_{m}$, in the phase diagram Fig. 10a the dynamics can be reduced to a classical one-component system (Eq. (41)). Such a system exhibits propagating waves within the region of bistability located between the two saddle-node (SN) bifurcations of local equilibria (gray area in Fig. 10a); compare the bistable (gray) area in Fig. 9 and Fig. 10b,c, where the nullcline slope is more negative than $-1$. Only for the fine-tuned value right at $\bar{n}=n_{\text{stat}}$, the front velocity is zero (purple dot). At $(n_{\text{stat}},D_{m})$, the dynamic scaffold that self-adapts via mass redistribution for $D_{c}\gtrless D_{m}$ bifurcates from the static scaffold $m_{\pm}(n_{\text{stat}})$ of the marginally stable front.

For bistable kinetics the slope at the inflection point of the nullcline, $s_{\text{nc}}^{\text{inf}}$, is necessarily more negative than $-1$ so stationary patterns may also exist for $D_{c}^{\text{min}}<D_{m}$ (cf. Eq. (40)), since they can be constructed from FBS-NC intersection points as shown in Fig. 10c. (We stick to the notation with concentrations $m$ and $c$, although they are not meaningful as “membrane” and “cytosolic” concentrations, in the case $D_{c}<D_{m}$. Instead they should be understood as abstract concentrations Note:Correlations ). The corresponding region where such stationary patterns exist is shaded in red in the bifurcation diagram shown in Fig. 10a. In the bottom half of this “balloon”-shaped region the equilibria that form the plateau of the constructed pattern are locally unstable. Hence these patterns cannot be stable. As we will see in the next subsection, all stationary patterns for $D_{c}<D_{m}$ are unstable (even if their plateaus are locally stable) since they are destabilized by the imbalance of reaction turnover induced by any (infinitesimal) perturbation. In contrast, stationary patterns for $D_{c}>D_{m}$ are stable because the self-adapting scaffold re-balances the reactive turnover.

The regions with a laterally unstable homogeneous steady state (NC-slope steeper than FBS-slope, cf. Eq. 25) are shaded in green and orange to distinguish in the bistable region which of the two locally stable reactive equilibria is laterally unstable; see Fig. 10b Note:SaddleNode .

In conclusion, we found the generic $(\bar{n},D_{c})$-bifurcation diagram of stationary patterns for an N-shaped nullcline with a bistable segment using the same geometric arguments as for the case of a monostable nullcline. Since our analysis crucially relies on the flux-balance subspace, it is limited to stationary patterns. For the special case $D_{c}=D_{m}$ the existence of non-stationary patterns (traveling fronts) in the locally bistable regime is well known. By continuity, we expect that there will be traveling fronts also $D_{c}\neq D_{m}$. Furthermore, we speculate that there will be non-stationary patterns outside the regime of local bistability, because mass redistribution may dynamically create a region of bistability that travels through the system.

VII.2.3 The dynamic scaffold self-balances by shifting the flux-balance subspace

In the following, we will assess the stability of the stationary (mesa) patterns found by the geometric construction above. (We only discuss the stability of mesa patterns, which are generic in the limit of large system size, $L\to\infty$.) The phase-portrait analysis in the phase space of chemical reactions facilitates a simple heuristic approach to study pattern stability: Instead of a full stability analysis of the stationary pattern, we consider only the stability of the FBS-position (mass-redistribution potential $\eta(x,t)$), as a proxy for the pattern stability. Intuitively, in the direction along the FBS, the pattern is quickly stabilized by due to scaffolding by local equilibria. In the following we present a simple stability criterion for stationary patterns, derived from this intuition. Details of the (ad hoc) derivation and a comparison to numerical analysis is presented in Appendix I. A mathematically rigorous stability analysis of stationary patterns (using, for instance, the “Singular Limit Eigenvalue Problem” introduced by Nishiura and Fujii (1987), see Ward (2006) for a survey), is outside the scope of this paper.

Recall that in steady state $\widetilde{\eta}(x)=\eta_{0}^{\infty}$ is spatially uniform and determined by total turnover balance, Eq. (17), that can be geometrically interpreted as a Maxwell construction (balance of the red-shaded areas in Figs. 2b and 11). How does the system evolve following a perturbation of the FBS-position $\eta$? Consider a spatially uniform shift, $\delta\bar{\eta}>0$. Then the area between NC and FBS to the left of $m_{0}$ (inflection point of the pattern) will decrease, while the area on the right will increase; see Figs. 2b and 11. The net reactive flow (“sum of the arrows in the two areas”) leads to a change of the average concentrations $m$ and $c$ (in the interface region) which amounts to a shift of the FBS-position, as shown in the insets in Fig. 11. Because the net reactive flow points along a reactive phase space (slope $-1$), the direction in which the FBS shifts due to turnover imbalance depends on its slope. For $D_{c}>D_{m}$, i.e. an FBS-slope larger than $-1$, it will move down and, hence, relax back to $\eta_{0}^{\infty}$ (recall that we considered an upwards shift as perturbation, the arguments work analogously for a downwards shift). In other words, the scaffold will adapt until total reactive turnover is balanced again. We conclude that the scaffold is self-balancing when $D_{c}>D_{m}$. Conversely, when $D_{c}<D_{m}$, the FBS will move further in the direction of the perturbation, thus destabilizing the pattern.

This qualitative stability argument can be expressed mathematically, to obtain a quantitative approximation for the growth rate of perturbations $\delta\eta(x_{0},t):=\eta(x_{0},t)-\eta_{0}^{\infty}$ in the vicinity of a stationary mesa pattern $(\widetilde{m}(x),\eta_{0}^{\infty})$ (see Appendix I):

[TABLE]

where $\tilde{f}_{\eta}=\partial_{\eta}\tilde{f}$. Comparison to numerically computed linear stability (dominant eigenvalue) of stationary patterns confirms shows that Eq. (43) is a good lowest order approximation for pattern stability; see Appendix I and Fig. 25. The integral over $\tilde{f}_{\eta}$ is the turnover imbalance due to a shift of the FBS, and thereby captures the geometric intuition based on the ‘Maxwell construction’ (area balance) we outlined above. The prefactor $(D_{m}/D_{c}-1)$ determines the direction in which FBS will shift because the integrand is always positive for $f_{c}>0$ ($\tilde{f}_{\eta}(m,\eta_{0}^{\infty})=f_{c}(m,\eta_{0}^{\infty}-m\,D_{m}/D_{c})>0$). We hence recover the stability criterion from our geometric argument above.

VII.3 The cusp scenario is generic

In some previous literature, it was argued that bistability of the reaction kinetics is an essential prerequisite for polarization to emerge in 2C-MCRD systems Jilkine and Edelstein-Keshet (2011); Mori et al. (2008). This claim was questioned recently Goryachev and Leda (2017).

Above, we conclusively showed that bistability is not necessary for pattern formation. Instead, in systems with conserved quantities, a (non-homogeneous) pattern scaffold can generically self-organize by due to shifting local equilibria when there is mass-redistribution ($D_{c}\neq D_{m}$). However, there is an interesting and more subtle connection between bistability and the ability to form patterns. This connection is revealed by studying the transition from monostable to bistable kinetics due to variation of kinetic parameters.

Variations in the kinetic parameters will change the shape of the reactive nullcline. This may not only lead to quantitative but also qualitative changes in the $(\bar{n},D_{c})$-bifurcation diagram, namely if there is a transition from a monostable to a bistable reaction kinetics: Imagine that variation of some rate $k$ in the reaction kinetics generates nullcline deformations as shown in Fig. 12a. Let us start with a nullcline that is strictly monotonically increasing with $m$. Then according to the geometric criterion, Eq. (25), there is no lateral instability and hence no stationary patterns; recall that we have shown in Sec. V.1 that an interface, the elementary element of a pattern, must be a laterally unstable region. Upon further changing the kinetic rate $k$ there may eventually be a threshold value $k_{\text{saddle}}$ beyond which the nullcline shows a region with a negative slope. A regime of lateral instability, and with it a regime where patterns can exist, emerges once the nullcline slope first becomes steeper than the slope of the FBS: $s_{\text{nc}}^{\text{crit}}=-D_{m}/D_{c}$. Eventually, at $k=k_{\text{cusp}}$, a further deformation of the nullcline may create a segment with slope $s_{\text{nc}}<-1$ where there is bistability of local equilibria. The bistable regime (shaded in gray) emanates in a cusp bifurcation where the two saddle-node (SN) bifurcations of the reactive equilibria meet in a single point, and the nullcline inflection point has slope $s_{\text{nc}}^{\text{inf}}=-1$ such that it is tangential to the reactive phase space $m+c=n_{\text{inf}}=n_{\text{cusp}}$; see Fig. 12b. The surface of reactive equilibria in the two-parameter bifurcation structure starts to fold over itself at the cusp point. Because the Turing bifurcations lie on different sheets of the folded surface of local equilibria, The Turing-bifurcation line (dash-dotted green line) crosses over itself in the bistable region of the bifurcation diagram.

At the cusp point, the topology of the $(\bar{n},D_{c})$-bifurcation diagram changes from the topology characteristic for monostable reaction kinetics (leftmost panel in Fig. 12c, cf. Fig. 8) to the one characteristic for bistable reaction kinetics (rightmost panel; cf. Fig. 10). At $k=k_{\text{cusp}}$, the bifurcation diagram has the singular topology shown in the center panel, with the cusp critical point at $(n_{\text{cusp}},D_{m})$.

Interestingly, by changing kinetic rates, the reactive nullcline of any system that is able to form stationary patterns, can be deformed to undergo such a cusp bifurcation of (local) reactive equilibria. To see how this works, first recall that to obtain stationary patterns, there must be a segment where the nullcline slope is $s_{\text{nc}}<-D_{m}/D_{c}<0$, i.e. necessarily negative. This segment can be deformed into a bistable region (slope $s_{\text{nc}}<-1$) by a dilation of the nullcline in the $c$-direction in phase space; see Fig. 13. Let us denote the scale factor by $\zeta$, then the dilated nullcline is determined by $f(m,\zeta^{-1}c)=0$. A slope $s_{\text{nc}}$ of the original nullcline $f(m,c)=0$ will become a slope $\zeta\cdot s_{\text{nc}}$ for the dilated nullcline. In particular, dilation by a factor $\zeta=|s_{\text{nc}}^{\text{inf}}|^{-1}$ leads to a nullcline with slope $-1$ at its inflection point, hallmark of a cusp bifurcation of the reactive equilibria (cf. Fig. 13). Any further dilatation leads to a bistable region (nullcline segment with slope more negative than $-1$). We conclude that the generic bifurcation scenario underlying all pattern-forming 2C-MCRD systems is a cusp bifurcation of reactive equilibria. Our geometric reasoning thus explains previous numerical observations Trong et al. (2014).

VII.4 Sub- and supercriticality of lateral instability in finite sized systems

So far, we have focused on the large system size limit ($L\to\infty$) where bifurcation diagrams can be constructed from phase-space geometry. In particular, we found that the onset of a mass-redistribution instability is generically subcritical. To analyze sub- vs. super-criticality in the case of finite system size, we use a perturbative approach (weakly nonlinear analysis, see e.g. Manneville (1990)) for the pattern close to onset. In the vicinity of the homogenous steady state $(m^{*},c^{*})$, we expand a stationary state $(\widetilde{m}(x),\eta_{0})$ in harmonic functions (eigenmodes of the Laplace operator under no-flux boundary conditions):

[TABLE]

As we have learned in Sec. IV, the band of unstable modes always extends to long wavelengths ($q\to 0$) in a 2C-MCRD system (‘type II’ instability, cf. Fig. 3b). Therefore, in a finite size system, the mode that becomes unstable first at the onset of the lateral instability is the longest wavelength mode: $\cos(\pi x/L)$. We want to study the steady-state amplitude of this mode in the vicinity of the onset bifurcation. It is not sufficient, though, to keep only this first harmonic in the mode expansion, Eq. (44), since higher harmonics couple to it through the nonlinear terms. For an expansion to third order in the first mode amplitude $\delta m_{1}$, one needs to include only the first and second harmonics in Eq. (44). Higher harmonics are not needed because they couple to $\delta m_{1}$ through higher nonlinearities $\mathcal{O}(\delta m_{1}^{4})$. To leading order the ansatz thus reads

[TABLE]

where $(m^{*},\eta^{*})$ denotes the homogenous steady state: $\tilde{f}(m^{*},\eta^{*})=0$. Using this ansatz in Eq. (14a) and keeping terms up to third order yields for the steady state pattern amplitude $\delta m_{1}$ (see Appendix J for details):

[TABLE]

The first order and third order coefficients read

[TABLE]

where $\tilde{\partial}_{m}=\partial_{m}-\frac{D_{m}}{D_{c}}\partial_{c}$ is the derivative along the direction of the flux-balance subspace.

The first harmonic amplitude $\delta m_{1}$, solution to Eq. (46) undergoes a pitchfork bifurcation at $F_{1}=0$. This bifurcation is simply the Turing bifurcation as it coincides with the onset of lateral instability: the homogenous steady state, $\delta m_{1}=0$, is laterally unstable only if the longest wavelength mode, $q_{1}=\pi/L$, is within the band of unstable modes: $\pi^{2}/L^{2}<q_{\mathrm{max}}^{2}=\tilde{f}_{m}/D_{m}$, i.e. if $F_{1}>0$ (cf. Eq. (47a)). Hence, the sign of the third order coefficient $F_{3}$, evaluated at the bifurcation point $F_{1}=0$, determines wether the bifurcation is supercritical ($F_{3}<0$) or subcritical ($F_{3}>0$); see Fig. 14b.

On the basis of this weakly nonlinear analysis, we can study the bifurcation at $F_{1}=0$ in any control parameter $\mu$ (for instance the average total density $\bar{n}$, the system size $L$, kinetic rates, and diffusion constants). With the critical value $\mu_{\text{c}}$, defined by the condition $F_{1}(\mu_{c})=0$, we introduce the reduced control parameter $\delta\mu=\mu-\mu_{\text{c}}$, and linearize Eq. (46) to lowest order in $\delta\mu$:

[TABLE]

To leading order, the branch of the solution that bifurcates at $F_{1}=0$ then reads

[TABLE]

At singular points where local stability and lateral stability change simultaneously, i.e. $\sigma_{\text{loc}}=0$ and $F_{1}=0$, the last term in $F_{3}$ diverges. Such codimension-two points require a more technically involved analysis (unfolding) that is outside the scope of this study (see e.g. Refs. Arnold (1988); Murdock (2003)).

From the third order coefficient $F_{3}$ (see Eq. (47b)), we can analyze the type of bifurcation in terms of geometric features (nullcline curvature $\kappa(n)\sim-\tilde{f}_{mm}|_{n}$, see Appendix K) together with the quantity that characterizes the relaxation rate to local equilibrium ($\sigma_{\text{loc}}$, cf. Eq. (5)). The first two terms in $F_{3}$ encode geometric properties of the reactive nullcline (curvature and its rate of change), i.e. how the local equilibria shift as mass redistribution shifts the local phase spaces. The last term in $F_{3}$ represents the rate of change of the timescale for relaxation $\sigma_{\text{loc}}$ to the local equilibria. In the following, we discuss the various regimes that arise due to the interplay of the three terms in $F_{3}$. (To simplify notation, we will implicitly assume that all coefficients ($\tilde{f}_{m}$, etc.) in $F_{3}$ are evaluated at the bifurcation point $\mu_{\text{c}}$.)

1. (i)

At the nullcline inflection point $\tilde{f}_{mm}=0$, only the first summand of $F_{3}$ remains: $F_{3}=\tilde{f}_{mmm}/8$. The third derivative $\tilde{f}_{mmm}$ is proportional to the rate of change of the curvature ($\tilde{f}_{mmm}\sim-\partial_{m}\kappa$, see Appendix K). The curvature of a typical N-shaped nullcline must be positive (bent upwards) to the right of the inflection point and negative (bent downwards) to the left of it, otherwise it is impossible to smoothly connect laterally unstable regions ($\tilde{f}_{m}<0$) to laterally stable regions ($\tilde{f}_{m}<0$). This implies that $\tilde{f}_{mmm}<0$, and we conclude that the bifurcation at the nullcline inflection point is supercritical. This confirms our geometric argument above (cf. Fig. 8, and thin gray lines in Fig. 14). 2. (ii)

In the large system size limit ($L\to\infty$) away from the inflection point ($\tilde{f}_{mm}\neq 0$), the second term in Eq. (47b) dominates. Because this term is always positive, the Turing bifurcation in large systems is generically subcritical (cf. Fig. 14) as we have already concluded from geometric arguments in Sec. V.4 above. 3. (iii)

In a finite-sized but still large system, the last two terms of $F_{3}$ are negligible sufficiently close to the nullcline inflection point where $\tilde{f}_{mm}$ vanishes. Hence, finite sized systems are supercritical in the vicinity of the nullcline inflection point, because $\tilde{f}_{mmm}$ must be negative as we argued above in point (i); (cf. Fig. 14). Solving the condition $F_{3}=0$ to to leading order in $L^{-1}$ yields the estimate

[TABLE]

for the range of supercriticality. 4. (iv)

In small systems, the transition from supercriticality to subcriticality will depend also on the last term in $F_{3}$. It can contribute either positively or negatively to $F_{3}$, depending on the details of the reaction kinetics.

Importantly, the statements (i)–(iii), regarding large systems, follow from purely geometric arguments as they are determined by the first two terms in $F_{3}$. The reason for this is that at large wavelength, chemical relaxation is fast compared to diffusion, so the pattern is slaved to the scaffold, i.e. the reactive nullcline. For the same reason, the long wavelength onset of lateral instability is determined by a geometric criterion (slope of the nullcline) as we have shown in Sec. IV.

In conclusion, we comprehensively characterized the Turing bifurcation (sub- vs. supercritical) and the bifurcations of stationary patterns, using the $(m,c)$-phase-space geometry. Because of the inherent link between geometry and the physical concepts of mass-redistribution and shifting equilibria, we are able understand the physics underlying the patterns and their bifurcations. The bifurcations in the large system size limit ($L\to\infty$)—determined by geometry—provide a good starting point to study the bifurcations in a finite sized system, e.g. by numerical continuation (see Appendix F).

VIII Conclusions and Discussion

We have introduced a local equilibria theory that allows to analyze and characterize both the initial growth and the eventual stabilization (saturation) of patterns on the basis of geometric objects in phase space: the shape of the reactive nullcline and its intersections with the flux-balance subspace. Within this framework, many properties of the nonlinear dynamics of mass-conserving reaction-diffusion systems can now be directly addressed in terms of phase space geometry, which would otherwise only be accessible by numerical analysis. In the following we will summarize the key concepts of local equilibria theory and important findings for 2C-MCRD obtained from this theory. Subsequently, we address new perspectives for the investigation of reaction–diffusion systems with conserved masses. Protein-pattern forming systems in vivo (intracellular) and in vitro (reconstituted systems) will serve as the specific context for this part of the discussion. Finally, we will give a brief outlook on upcoming work, on questions that are currently under investigation, and on future research directions. In particular, we outline how local equilibria theory might provide a unifying geometric perspective on pattern formation in mass-conserving non-equilibrium systems and how it can be generalized to systems that are not strictly mass-conserving but contain a mass-conserving ‘core’.

VIII.1 Summary of key concepts and results

Phase space analysis of (two-component) MCRD systems has shown that spatial variations in protein density give rise to spatially heterogeneous local equilibria. The relationship between mass and reactive equilibria is represented geometrically as a line of reactive equilibria (reactive nullcline) in phase space. Along the reactive nullcline, the reaction kinetics are balanced. We have shown that this nullcline is a central geometric object in phase space that organizes the spatiotemporal dynamics. Furthermore, we have identified a one-parameter family of (one-dimensional) manifolds in phase space on which diffusive fluxes balance. Any stationary pattern is embedded in one of these so-called flux-balance subspaces. The variable that parameterizes the family of flux-balance subspaces acts as a mass-redistribution potential: Its spatial gradients represent a local imbalance of diffusive fluxes that drives mass-redistribution. Thus, the mass-redistribution potential encapsulates the interplay between reactions and diffusion processes.

In this way, the spatiotemporal dynamics of the reaction–diffusion system is fully determined by geometric structures in phase space: (i) We introduced a flux-balance construction, based on intersection points between flux-balance subspace and reactive nullcline that enabled us to graphically construct stationary patterns and their complete bifurcation structure in the limit of large system size. Underlying this construction is the general insight that patterns are scaffolded by local equilibria which are, in turn, encoded by the reactive nullcline. This principle can be generalized to dynamics: (ii) The interplay between diffusive redistribution of mass and shifting local equilibria drives the pattern-forming instability that we termed mass-redistribution instability. Our analysis has shown that the slope of the nullcline provides a simple criterion for the occurrence of this instability. Importantly, we find that the onset of instability is generically subcritical in 2C-MCRD systems on a large domain. (iii) Generalizing the local equilibria concept, we have introduced a decomposition of the spatial domain into regions, such as plateaus and interfaces, which are characterized in terms of regional dispersion relations. This has enabled us to find simple heuristics for many properties of the pattern formation dynamics and stationary patterns. For instance, the width of the interface region can be approximated by the marginal mode of the regional dispersion relation at the interface. Building on this regional decomposition, one can characterize different pattern types and the transitions between them as control parameters are changed (Sec. V.3). Furthermore, based on the concept of regional (in)stability we found an inherent connection between lateral (Turing) instability and stimulus-induced pattern formation (“nucleation and growth”), which enabled us to estimate the basin of attraction (“nucleation threshold”) for stationary patterns by a simple heuristic using the reactive nullcline. As an additional advantage of such a characterization, we note that the reactive nullcline could in principle be determined experimentally for any given system in which the average total density can be controlled in a well-mixed “reactor”.

Importantly, the concepts of scaffolding by local equilibria (i), mass-redistribution instability (ii) and regional dispersion relations (iii) are directly generalizable to systems with more components and more conserved quantities. We extensively discuss these future directions in the Outlook below.

VIII.2 Reaction–diffusion systems

For the sake of specificity, we will discuss the implications and application of our findings on mass-conserving reaction–diffusion systems with respect to protein pattern formation, which operates far from thermal equilibrium, and has received growing interest over the past two-decades. Intracellular, i.e. in vivo, protein pattern formation and self-organization have been subject to a large body of research, both experimentally (see Refs. Campanale et al. (2017); Chiou et al. (2017) for recent reviews) and theoretically (see Ref. Halatek et al. (2018) and references therein.) Furthermore, the in vitro reconstitution of the MinDE system Loose et al. (2008) has made it possible to study protein pattern formation experimentally under a wide range of externally controllable conditions; see Refs. Ivanov and Mizuuchi (2010); Zieske and Schwille (2014); Vecchiarelli et al. (2014); Caspi and Dekker (2016); Kretschmer et al. (2017); Denk et al. (2018); Glock et al. (2018); Litschel et al. (2018); Ramm et al. (2018); Kohyama et al. (2019); Glock et al. (2019) and Ref. Vendel et al. (2019) for a recent review.

Taken together, these studies of both in vivo and in vitro systems have led to many important insights. However, many intriguing questions that are relevant to all reaction–diffusion systems far from equilibrium, remain open: What is the role of the (biomolecular) interaction network, and how can complex models be reduced to their essential components? What are the physical mechanisms underlying the pattern-forming instabilities and under which conditions do these instabilities arise? How can the dynamics of patterns far from the homogeneous steady state be studied systematically, i.e. how can we bridge the gap between the linear and the highly nonlinear regime? A particular question in this context is how different patterns and their characteristic length scales (interface width and wavelength) are selected in the highly nonlinear regime. In what follows, we discuss the implications of our work to these questions.

VIII.2.1 Model classification, network motifs and experimental accessibility

In recent years, several studies have employed high-throughput computational analyses of reaction–diffusion systems and graph theoretical analysis with the goal to infer the pattern-forming capabilities from the topology of the underlying reaction networks Chau et al. (2012); Marcon et al. (2016); Sugai et al. (2017); Diego et al. (2018). Our results offer an entirely new and distinct perspective on model classification and the role of the interaction-network topology for mass-conserving systems. We found a simple condition for the pattern-forming (mass-redistribution) instability in 2C-MCRD systems: The slope of the line of reactive equilibria (reactive nullcline) must be (sufficiently) negative. Broadly speaking, the reactive equilibrium of the faster diffusing (i.e. cytosolic) component has to decrease with increasing total density (cf. Sec. IV). Importantly, our approach goes beyond the classification based on linear (in)stability. It shows that the effect of nonlinearities on the dynamics is encoded in the curved shape of the nullcline. In particular, there is a direct connection between the nullcline shape and the characteristic spatial density profile of the pattern (cf. Sec. V.3). The reactive equilibria, as represented by the reactive nullcline, might therefore provide an alternative approach to model classification. Hence, a key challenge for future research will be to study how specific reaction kinetics and model parameters affect the shape of the reactive nullcline.

Moreover, a major advantage of reactive equilibria as the essential criteria for model classification is their experimental accessibility. In principle, any line of reactive equilibria can be measured directly in experiments by using a single, isolated, and well-mixed reactor and externally controlling the available conserved quantity (e.g. particle number). Such experiments would allow one to probe and classify the core mechanism quantitatively without any knowledge of the molecular details (which are irrelevant for such a classification).

The concepts that underlie local equilibria theory—mass redistribution and moving local equilibria—are not restricted to two-component systems with a single conservation law. They have previously been applied to the model of the (in vitro) MinDE system, which has two conserved protein species, MinD and MinE, and five components Halatek and Frey (2018); Brauns et al. (2020b). We believe that the results presented here are foundational for the development of a more general theory. Our analysis constitutes the first step in a long-term project to find a geometric representation of the nonlinear dynamics of spatially extended systems. In the outlook, Section VIII.4, we briefly describe various forthcoming works and future projects that build upon the present study and generalize its results.

VIII.2.2 Polarity patterns, bistability and the necessary condition for a Turing instability

Bistability is a generic feature of nonlinear systems, and its putative relation to polarity patterns has been controversially discussed in the literature Trong et al. (2014); Goryachev and Leda (2017); Mori et al. (2008); Cusseddu et al. (2018); Gamba et al. (2009). In systems without mass conservation, bistable reaction kinetics facilitate traveling fronts that connect the two stable reactive equilibria (homogeneous steady states). (This scenario appears in mass-conserving systems as the special case of equal diffusion constants, which entails the lack of lateral mass redistribution; see Sec. VII.2.) The stable equilibria of such a bistable medium can be pictured as scaffold for the traveling front. Because this scaffold is static in systems without mass redistribution, fine-tuning is required to achieve a stationary front (cf. Sec. VII.2.1). Our results show that, in mass-redistributing systems (i.e. systems with unequal diffusion constants), the scaffold becomes dynamic and thereby supports stable, stationary polarity patterns in an extended parameter regime (cf. Sec. VII.2). Most importantly, we found that, even in a monostable system, mass redistribution can facilitate the formation of a scaffold for stationary polarity patterns. Hence, bistability of the reaction kinetics is neither required nor sufficient for the formation of such patterns.

A central finding of our work is the physical mechanism by which the scaffold of a pattern emerges dynamically from a homogeneous steady state: the mass-redistribution instability (Sec. IV). Diffusive mass-redistribution requires that the diffusivities of the two components are different (“differential diffusion”); already a ratio $D_{m}/D_{c}$ slightly different than unity is sufficient. This mass-redistribution drives an instability under the condition that the reactive nullcline in the phase space of the reaction kinetics includes a segment of negative slope. More precisely, the negative slope must be steeper than the flux-balance subspace, whose slope is determined by the negative ratio of the diffusivities $-D_{m}/D_{c}$. As membrane–bound proteins are significantly less mobile than cytosolic proteins, i.e. $D_{m}\ll D_{c}$, a small negative slope of the reactive nullcline is already sufficient for a Turing instability in the intracellular context. Recall also that the criterion for bistability is a NC-slope more negative (“steeper”) than $-1$ (cf. Sec. III.1), which is obviously a more restrictive condition than that for lateral instability. Hence, a bistable region is generally surrounded by a larger region of lateral instability in parameter space (cf. Fig. 12).

One might wonder how generic nullclines with a segment of negative slope are. In fact, they are frequently encountered as N-shaped nullclines in a broad range of classical nonlinear systems Strogatz (2014); Izhikevich (2007); Novák and Tyson (2008); Murray (2002). Typically, these nullclines encode some nonlinear feedback mechanisms that give rise to widespread phenomena such as relaxation oscillations, excitability and bistability. Feedback loops in intracellular signaling networks (e.g. GTPase and phosphorylation cycles) generically lead to these phenomena Kholodenko (2006). From this perspective, it appears that N-shaped nullclines should be rather common in biological systems.

Finally, the condition for mass-redistribution instability shows that, contrary to popular belief, Turing instability in a mass-conserving system can already arise for a ratio $D_{m}/D_{c}$ only slightly below unity when the reactive nullcline has sufficiently steep segment (i.e. in the vicinity of saddle-node bifurcations of the reactive equilibria, cf. Fig. 12).

VIII.2.3 Mass-redistribution instability is mechanistically distinct from the activator–inhibitor paradigm

A key finding of our work is the physical mechanism of mass-redistribution instability: Shifting local equilibria induce gradients and hence (diffusive) mass-redistribution which, in turn, leads to further shifting of the local equilibria. Notably, mass-redistribution instability is a Turing instability in the original sense: a lateral instability due to diffusive coupling of a system that is locally stable in the absence of this coupling.

Importantly, the mechanism underpinning the mass-redistribution instability is distinct from the paradigm of short range activation and long range inhibition Gierer and Meinhardt (1972); Segel and Jackson (1972). The latter, termed “activator–inhibitor” mechanism, innately relies on some form of “production” and “degradation” processes, either of “morphogens” via gene-regulatory networks (e.g. in cell-cell signaling Karig et al. (2018) and tissue patterning Manukyan et al. (2017); Watanabe and Kondo (2015); Kondo and Asai (1995)), or of electric signals coupling excitory and inhibitory neurons (e.g. in the visual cortex where lateral instability may underlie hallucinations due to long range coupling of inhibitory neurons Ermentrout and Cowan (1979)). In contrast, intracellular pattern formation is often driven by (nearly) mass-conserving dynamics Halatek et al. (2018) and, hence, is necessarily based on a mass-redistribution instability. Because, as we emphasized above, mass-redistribution instability is a Turing instability, self-organized intracellular protein-patterns (see Ref. Halatek et al. (2018) for a review) are examples for Turing patterns (in the sense of patterns that arise from a Turing instability Note:Turing ) in biological systems.

In general, it might not always be obvious whether the instability underlying pattern formation in a given system is either a mass-redistribution instability or one that essentially requires production and degradation. In the outlook Sec. VIII.4.3, we propose a general definition to answer this question.

VIII.2.4 Subcriticality and stimulus-induced pattern formation

An important results of our analysis is that the onset of pattern formation in 2C-MCRD systems is generically subcritical (cf. Secs. V.4 and VII.4). Subcriticality may be beneficial in a biological context, as it confers robustness: Once a pattern is established, it is robust towards parameter variations due to hysteresis Eldar et al. (2004); Goehring et al. (2011); Trong et al. (2014).

Moreover, subcriticality implies the existence of parameter regimes where pattern formation can be triggered by sufficiently large perturbations (akin to “nucleation and growth” in the binodal regime of phase-separating systems). Such stimulus-induced pattern formation has been suggested as a new mechanism for pattern formation (under the term “wave-pinning”), which—it was argued—is fundamentally distinct from a lateral (Turing) instability Mori et al. (2008); Jilkine and Edelstein-Keshet (2011); Edelstein-Keshet et al. (2013). This claim has been disputed in recent works Goryachev and Leda (2017); Verschueren and Champneys (2017). Indeed, our results show that linear instability and stimulus-induced pattern formation are inherently connected: the latter is possible only where there is an adjacent regime of linear lateral instability and the underlying mechanism is a regional (linear) lateral instability; see Sec. VI. Concretely, our results show that an interface—the elementary building block of a pattern—must necessarily be a laterally unstable region (cf. Sec. V.1). Hence, the creation of a laterally unstable region is a necessary condition for the formation of a stationary pattern. This implies that, any two-component system that has a regime of stimulus-induced pattern formation, must also exhibit a regime where the homogenous steady state is laterally unstable, and this regime can always be reached by simply changing the average total density. Conversely, this suggests that subcriticality may be a generic feature of mass conserving systems since regional instability will facilitate stimulus-induced pattern formation adjacent to regimes where the homogeneous steady state is laterally unstable.

Finally, building on the concept of regional instability insights, we provided a simple geometric argument for the perturbation threshold for stimulus-induced pattern formation (“nucleation threshold”) based on the reactive nullcline.

VIII.2.5 Length scale selection

An important consequence of subcriticality is that well-known mathematical results for systems near a supercritical instability may not apply anymore. Potentially the most prominent of such results is the existence of a characteristic wavelength, determined by the fastest growing mode in the dispersion relation. This is often considered as a defining property of “Turing patterns”. However, for subcritical systems, the wavelength of the pattern cannot be inferred from a linear stability analysis of the uniform steady state in general—not even at onset. Indeed, 2C-MCRD systems always exhibit uninterrupted coarsening Brauns et al. (2020a), i.e. the wavelength selected by the fastest growing mode at onset is observed only transiently and the final steady state is fully “phase separated.”

Note, however, that coarsening is not a generic feature of all mass-conserving systems. Multicomponent MCRD models show (multi-)stable patterns with finite wavelengths Halatek and Frey (2012); Murray and Sourjik (2017); Halatek and Frey (2018). We believe that identifying and understanding mechanisms of nonlinear wavelength selection that bring the coarsening processes to a halt and stabilize patterns with finite wavelength are among the most important tasks for future research on multicomponent models. While some rather general criteria have been found for one-component systems Politi and Misbah (2004, 2006), a comprehensive understanding of multicomponent systems remains out of reach for now.

In addition to the wavelength, patterns have a second characteristic length scale—the width of interfaces. While the wavelength of patterns far from the homogeneous steady state is not determined by the dispersion relation, we have shown in Section V.1 that the interface width is determined by the marginal mode ($q_{\mathrm{max}}$) of the interface’s regional dispersion relation. The more general insight underlying this finding is that highly nonlinear patterns can be decomposed into spatial regions and studied in terms of regional phase-spaces and regional attractors (cf. Sec. V.2). In a follow-up work Brauns et al. (2020b), we use this region decomposition to classify the different types of instabilities that govern Min-protein pattern formation in vivo and in vitro. In particular, we find that the interface width of standing wave patterns is also determined by the marginal mode of the regional dispersion relation.

VIII.2.6 Pattern types are determined by the nullcline shape

Two-component MCRD systems are found to exhibit two generic pattern types, referred to as mesas and peaks. In a recent numerical study of 2C-MCRD models Chiou et al. (2018), a phenomenological “saturation point” was found to mark the peak- to mesa-pattern transition. Our geometric analysis has now revealed the phase-space structure that underlies this “saturation point”: It corresponds to an intersection point between reactive nullcline and diffusive flux-balance subspace; hence, its position dictated by the nullcline shape and the ratio of the diffusion constants. Based on this insight, have shown that the nullcline shape serves as a simple criterion that to predicts the type of patterns formed by a given 2C-MCRD system; see Fig. 5.

Put briefly, we distinguish $\Lambda$- and N-shaped nullclines, based on their tail behavior for large densities. In the case of a $\Lambda$-shaped nullcline (e.g. nullclines that asymptotically approach the $m$-axis for large $m$; see Figs. 5b and 17) and a shallow flux-balance subspace ($D_{m}/D_{c}\ll 1$), a mesa pattern forms only when the average mass is sufficiently large. For lower densities, the system exhibits peak patterns instead. Their amplitude depends sensitively on the total mass and the membrane diffusion constant. The approximation of the peak pattern amplitude provided in Sec. V.3 (details in Appendix G), provides a simple estimate of the total (protein) mass at which peak patterns transition to mesa patterns. Notably, the amplitude of peak patterns and the transition to mesa patterns depends sensitively on the ratio of the diffusion constants. In the case of an N-shaped nullcline (see Figs. 2 and 5a), mesa patterns are generic, because the high-density plateau is already formed at low average mass (even if $D_{m}/D_{c}\ll 1$). (Peak patterns are also possible for such nullclines, but require fine tuning of the average mass to the vicinity of bifurcation points (see Fig. 6).)

These findings have important biological consequences because the characteristic features of a pattern dictate the positional information it can convey, and therefore the biological function it can facilitate Wolpert (1969); Alberts et al. (2007); Lutkenhaus (2007); Bendezu and Martin (2012); Halatek et al. (2018). For example, cell division in budding yeast (S. cerevisiae) requires the formation of a single, narrow polarity site marked by a high density of the protein Cdc42, to uniquely determine the future bud site Park and Bi (2007); Martin and Arkowitz (2014). This requirement is met by peak patterns. Mesa-like patterns, in contrast, sharply separate two spatial domains. This is, for instance, a feature of PAR-protein patterns in C. elegans Goehring et al. (2011); Hoege and Hyman (2013); Gross et al. (2018); Geßele et al. (2020). Interestingly, the stem cells that polarize via PAR-protein pattern formation become smaller after each cell division during morphogenesis of the C. elegans embryo Hubatsch et al. (2019). As the cell size approaches the interface width, the patterns transition from the mesa to the ‘weakly nonlinear’ type (cf. Sec. VII.4). Finally the system size becomes so small that the Turing instability is suppressed (see discussion of the marginal mode $q_{\mathrm{max}}$ in Sec. IV). This size-dependent loss of cell polarization has been shown to be important for stem-cell fate decision Hubatsch et al. (2019).

An important difference between peak and mesa patterns is how they respond to changes in average total protein density, e.g. owing to up- or down-regulation of gene expression, or system size, e.g. due to growth. An increase in average mass makes peak patterns grow in amplitude, while for mesa patterns it leads to a shift in the interface position. Upon an increase in system size, peak patterns grow in amplitude, because the total number of proteins in the system increases. In contrast, amplitude and the relative interface position of mesa patterns remain unchanged. Hence, mesa patterns inherently scale with the system size—a property that is desirable in developmental systems. Note, however, that the interface width of mesa patterns does not scale with the system size, and is independent of the average total density (Sec. V.3).

Another difference between patterns composed of peaks vs patterns composed of mesas is their rate of coarsening by competition for mass. Peak patterns coarsen significantly faster than mesa patterns Brauns et al. (2020a); Chiou et al. (2018). Fast coarsening is important if the biological function requires selection of a single polarity site, as for instance cell division of budding yeast Park and Bi (2007); Martin and Arkowitz (2014).

VIII.3 Non-equilibrium phase separation

On the phenomenological level, the dynamics of 2C-MCRD systems closely resembles a phase separation process, as exhibited by binary mixtures that undergo liquid-liquid phase separation near thermal equilibrium. We showed that the bifurcation diagram of 2C-MCRD systems (Fig. 8), obtained by the flux-balance construction on the reactive nullcline, resembles the phase diagram of phase separation with spinodal and binodal lines that meet in a critical point. The analogous process to spinodal decomposition is the mass-redistribution instability (Sec. IV). In both cases, the condition for instability is that a potential decreases as a function of the total density. In the former case, it is the chemical potential, derived from a free energy functional, while in the latter case it is the mass-redistribution potential $\eta^{*}(n)$ derived from the reactive nullcline. Moreover, nucleation in the binodal regime of binary-mixture phase separation corresponds to stimulus-induced pattern formation in 2C-MCRD systems (Sec. VI). In addition to their equivalent phase diagrams, binary-mixture phase separation and 2C-MCRD systems and both exhibit uninterrupted coarsening Brauns et al. (2020a).

This phenomenological equivalence between binary-mixture phase separation and 2C-MCRD dynamics is quite remarkable, since the former describes systems close to thermal equilibrium, while the latter is inherently far from thermal equilibrium (driven by a chemical fuel like ATP/GTP in the case of proteins). In fact, segregation into domains of high and low density is observed in many other non-equilibrium systems, most prominently self-propelled particles that exhibit formation of polar waves and nematic lanes Chaté (2020); Huber et al. (2018) and motility-induced phase separation (MIPS) Cates and Tailleur (2015). Further examples include active contractility Bois et al. (2011), motile bacteria Liu et al. (2019), and shear banding Olmsted (2008). Some authors have used the term “active phase separation” John and Bär (2005); Bergmann et al. (2018); Agudo-Canalejo and Golestanian (2019) for such phenomena.

One interesting feature shared by many phase-separating systems is that they exhibit interrupted coarsening (‘micro-phase separation’) once the mass conservation is weakly broken by additional production and degradation terms. The strength of these terms determines the length scale where coarsening arrests. This holds true both near thermal equilibrium Glotzer et al. (1995); Oono and Shiwa (1987) and for MIPS of self-propelled particles subject to a birth–death process Cates et al. (2010); Curatolo et al. (2019). In a followup work to the present manuscript, we show that the same holds true for two-component reaction–diffusion systems with weakly-broken mass conservation Brauns et al. (2020a). Interestingly, interrupted coarsening can also occur as a consequence of non-reciprocal coupling that destroys the variational nature of the Cahn–Hilliard model Saha et al. (2020).

Given the phenomenological equivalence to phase separation near thermal equilibrium, one might be tempted to search for a mapping to an effective thermodynamic language or develop an entirely phenomenological thermodynamic description. In fact, for 2C-MCRD systems with a specific form of reaction kinetics, $f\,{=}\,c\,{-}\,g(m)$, an effective free-energy functional can be constructed Morita and Ogawa (2010); Goh et al. (2011); Jimbo and Morita (2013), giving a gradient-flow structure to the 2C-MCRD dynamics with such a reaction term. In this specific case, some of the results presented here, like the phase diagram with binodals and spinodals, can be inferred directly from the mapping to equilibrium phase separation (see e.g. Ref. Goh et al. (2011)). Moreover, the mapping implies uninterrupted coarsening for this specific form of reaction kinetics.

However, such an approach disregards the actual underlying non-equilibrium physics. Our analysis of 2C-MCRD systems shows how a framework can be developed that is rooted in the underlying physics (here chemical reactions and diffusion) and not subject to the restrictions of a mapping to an effective thermodynamic description. Concretely, we show in a followup work that uninterrupted coarsening is generic in 2C-MCRD systems independently of the specific reaction kinetics Brauns et al. (2020a). Moreover, a thermodynamic description cannot account for the rich phenomenology, including and oscillatory patterns and waves Bernitt et al. (2017); Kuwamura and Morita (2015); Kuwamura and Izuhara (2017); Radszuweit et al. (2013); Weber and Frey (2017); Halatek and Frey (2018); Brauns et al. (2020b), that arises once coupling to additional components or source terms breaking mass-conservation are included. In contrast, local equilibria theory has proven useful also in these more complex scenarios Halatek and Frey (2018); Brauns et al. (2020a, b). We further discuss this perspective in the final two subsections of the Outlook.

Instead of describing attachment and detachment of proteins at a membrane, Eq. (1) can be interpreted as the mean field equation for particles undergoing Brownian motion and switching between two states with different velocities/tumbling rates. Then the reaction term $f=f(n)$ describes the switching dynamics that depends on the local density of particles (e.g. by some quorum sensing mechanism), and $m$ and $c$ are the concentrations of slow and fast diffusing particles, respectively. This system constitutes a minimal example for MIPS and shows that MIPS and Turing instability are analogous on the mean field level. MIPS is typically studied for particles with continually varying velocity that depends on the local particle density. For particles switching between two states with different velocities, the mean velocity varies continuously as a function of the total density.

Hence, we define the average diffusion constant

[TABLE]

which is directly connected to the mass-redistribution potential via the identity $\eta=n\bar{D}/D_{c}$. With this, the nullcline-slope condition for lateral instability, $\partial_{n}\eta^{*}<0$, can be recast as

[TABLE]

(As always, the star denotes evaluation at the reactive equilibrium; the function $\bar{D}^{*}(n)$ can be obtained from the reactive nullcline, $\eta^{*}(n)$, via the relation $\bar{D}^{*}(n)=D_{c}\eta^{*}(n)/n$.) Notably, the condition for MIPS has exactly the same form as Eq. (52), where the density dependent particle velocity $v(n)$ takes the place of the average diffusivity $\bar{D}^{*}(n)$.

This reveals the common underlying principle of Turing instability in MCRD systems and the instability driving MIPS of self-propelled particles: slowing down of particles in regions of high density. For an instability to occur, particles have to slow down enough in response to an increase in density, cf. Eq. (52).

Another phenomenon that can be pictured as a phase separation process is shear banding in complex fluids. In Appendix L, show how one can establish an analogy between our flux-balance construction on the reactive nullcline for 2C-MCRD systems and the “common total stress” construction on the constitutive relation employed to analyze shear banding  Olmsted (2008); Divoux et al. (2016). This analogy shows that these physically distinct systems are topologically equivalent, i.e. share the same phase space geometry.

VIII.4 Outlook

Based on the theoretical concepts presented in this work for 2C-MCRD systems there are several promising directions for future research. First of all, for 2C-MCRD systems several follow-up works elucidate the principles at work during wavelength selection Brauns et al. (2020a) and address how self-organization may be controlled by spatiotemporal “templates” Wigbers et al. (2020a) and by advective flows Wigbers et al. (2020b). Going forward, it will be interesting to explore various avenues toward generalization to MCRD systems with more components and conserved species, including applications to various specific physical and biological systems as outlined shortly below. The long-term perspective is a generalization toward a geometric theory of MCRD systems. Below, we discuss some possible routes towards such a generalization in more detail.

In a broader context, reaction–diffusion systems are part of a large class of non-equilibrium systems that are able to form self-organized patterns. This includes models for living matter where molecular motors generate active flows Bois et al. (2011); Kumar et al. (2014) and active visco-/poroelastic deformations of the medium Radszuweit et al. (2013); Weber et al. (2018), as well as particle-based active matter Marchetti et al. (2013); Cates and Tailleur (2015); Chaté (2020) and granular media Aranson and Tsimring (2006). Ultimately, it may be fruitful to apply the concepts presented here to such systems.

VIII.4.1 Generalization to more complex phenomena

The 2C-MCRD system studied here has a comparatively simple phenomenology, exhibiting only steady states that are stationary (no oscillations) and no wavelength selection (uninterrupted coarsening). In a previous study Halatek and Frey (2018), and recently in Ref. Brauns et al. (2020b), the concepts of mass redistribution and local equilibria have already been successfully employed to analyze a multicomponent multispecies MCRD model (for MinDE in vitro pattern formation) that exhibits much more complex phenomena, like spatiotemporal chaos (chemical turbulence) at onset and a transition to order (standing and traveling waves). These phenomena, observed in numerical simulations, can be understood in terms of the changing local stability of local equilibria due to mass redistribution.

In particular, this study revealed intriguing, highly non-trivial, connections between the nonlinear pattern dynamics far from the homogeneous steady state and the dispersion relation that characterizes the linearized dynamics in the vicinity of the homogeneous steady state. However, these findings are model specific, and rely on numerical simulations. In contrast, the characterization of the 2C-MCRD systems presented here is independent of the specific model (reaction term) and enables us to predict the pattern formation dynamics from a simple graphical analysis, without the need to perform numerical simulations. The obvious next step is to generalize this level of understanding to more complex phenomena by studying three-component MCRD systems. Such models have recently been employed to model various (bio-)physical phenomena in numerical studies Murray and Sourjik (2017); Alonso and Bär (2010). Studying these models using local equilibria theory might reveal the principles underlying their dynamics and provide a good starting point for a generalization of the theory presented here.

For ODE dynamics, the seemingly small step from two to three variables increases the diversity of phenomena dramatically. We expect a similar situation for spatially extended systems. Note in particular that, with three components and one conserved quantity, the reactive phase space will be two-dimensional. This allows for more complex local dynamics and attractors, such as limit cycle oscillations. A fully general study of three-component MCRD systems will therefore probably not be possible from the outset. Instead we propose to focus on cases where a time-scale separation enables one to build on the results for two-component systems. A good starting point is to study cases where the coupling to the additional third component is slow (the model investigated in Ref. Murray and Sourjik (2017) is of that form) using, for instance, a singular perturbation analysis. Closely related to such three-component MCRD systems are nearly mass-conserving two-component systems that contain production/degradation terms on a slow timescale (see next subsection). Making use of such timescale separations is has proven to be a powerful strategy to study complex phenomena in dynamical systems; see for instance Ref. Izhikevich (2000) for a comprehensive overview of three-variable ODE systems in the context of neural excitability.

The benefit of such an approach is that the effects of strong nonlinearities may be captured by a geometric phase-portrait analysis of the fast two-variable subsystem. More complex behavior arises as the slow dynamics modifies the fast subsystem, driving it through bifurcations. Analogously, in reaction–diffusion systems, the effects of the strong nonlinearities are encoded in the shape of the nullcline, which enables one to construct the elementary patterns (mesas, peaks). Much richer phenomenology can arise by the modification of these elementary patterns due to additional linear (or weakly nonlinear) terms, whose effects could be studied using the method of regional phase spaces (Sec. V.2). This can be viewed as a dual (“strong coupling”) approach to the amplitude equation formalism (weakly nonlinear analysis), where the elementary pattern originates in the narrow band of unstable modes, and the modification of the pattern due nonlinearities that couple these modes perturbatively.

Another promising direction of study is the role of noise in stochastic MCRD systems. Noise-induced phenomena in reaction–diffusion systems have previously been studied for models with (e.g. Altschuler et al. (2008); Jilkine et al. (2011)) and without conserved quantities (e.g. Sagués et al. (2007); Schumacher et al. (2013); Biancalani et al. (2017); Karig et al. (2018); Adamer et al. (2018)). Importantly, noise in reaction–diffusion systems is not constrained by the fluctuation-dissipation theorem, but must be inferred from the stochasticity of the chemical processes by explicit coarse graining, using, for instance, path integral approaches Weber and Frey (2017).

VIII.4.2 Model reduction and classification

While multicomponent, multispecies models can exhibit complex phenomenology, as discussed above, this is not necessarily the case. For instance, cell polarization of eukaryotic cells is often brought about by a large number of interacting protein species. Well-studied examples are the Cdc42 system of S. cerevisiae Wedlich-Soldner (2003); Goryachev and Pokhilko (2008); Klünder et al. (2013); Chiou et al. (2018) and the PAR system of C. elegans Goehring et al. (2011); Trong et al. (2014); Geßele et al. (2020); Gross et al. (2018). Despite the complexity of these systems, their phenomenology—cell polarization—is simple and can already be captured by two-component systems as studied here. This raises the question whether complex models can be reduced to an underlying minimal ‘core’ that captures their essential phenomenology. The finding that redistribution of conserved quantities is the essential driver of pattern formation in MCRD systems suggests that such a reduction might be possible in the phase space of conserved quantities—the control space. Note that such reduced ‘core’ systems may comprise more than two components and more than a single conserved mass. As we detail in the next subsection, even non-conservative models can have a mass-conserving ‘core’.

The central theme of local equilibria theory is that conserved quantities are control parameters for reactive equilibria. Hence, the bifurcation scenario of reactive equilibria in the control space may serve as a criterion for classifying models. One example of such a class is the cusp bifurcation scenario, found by numerical analysis of various cell-polarization models in Ref. Trong et al. (2014), and identified here as the general bifurcation scenario underlying pattern formation in 2C-MCRD systems. In a forthcoming work, we will use our theory to elucidate the control-space geometry underlying the pole-to-pole oscillations of the in vivo MinDE system Brauns et al. (shed)). As there are two conserved quantities, the total densities of MinD and MinE respectively, the control space is two-dimensional and there are surfaces (instead of lines) of reactive equilibria. Using the local quasi-steady state approximation, these nullcline surfaces allow for a geometric analysis of the in vivo MinDE dynamics. We propose this geometrically motivated approach as an alternative to algebraically motivated model reduction methods, such as the quasi-steady-state approximation of slowly diffusing components Smith and Dalchau (2018) or (extended) center-manifold reduction Sewalt et al. (2015); Doelman et al. (2015). Such an approach offers the advantage that it does not require abstract mathematical calculations and instead enables one to gain physical intuition from elementary geometric objects and graphical constructions. Furthermore, as pointed out earlier, reactive equilibria in principle allow one to assess phase-space geometry experimentally. This could ultimately make it possible to infer theoretical models from experimental data at a mesoscopic level, especially in situations where access to more molecular information (at the protein level) is not available yet.

VIII.4.3 Beyond strict mass conservation

What can we learn from local equilibria theory about systems without strict mass conservation?

Non-equilibrium systems are dissipative, that is, they consume some sort of chemical fuel (e.g. ATP in biological systems and malonic acid in the Belousov–Zhabotinsky (BZ) reaction Zhabotinsky (1964)) that drives them far from thermal equilibrium. The chemical fuel often drives cycling of components between different states, as illustrated in Fig. 15. Examples are cycling of NTPase proteins between active and inactive states Neves (2002); phosphorylation–dephosphorylation Hansen et al. (2019); Alonso and Bär (2014) and membrane attachment–detachment of proteins Halatek et al. (2018); and cycling between molecular bromine and bromide ions and cycling of a metal catalyst (e.g. cerium) in the BZ reaction Belousov (1959); Zhabotinsky (1964); Field et al. (1972). While such cycles consume a fuel (and produce a waste), they conserve the total density of the cycling components ($\mathrm{C}_{1}$ and $\mathrm{C}_{2}$ in Fig. 15). The chemical fuels are quite generally assumed to be supplied from a large reservoir and hence are not explicitly modeled. Several recent works study explicitly the role of a finite fuel supply in the framework of stochastic thermodynamics Rao and Esposito (2018); Avanzini et al. (2019); Ehrmann et al. (2019). In particular, in Ref. Ehrmann et al. (2019), the bifurcations of bistable and oscillatory well-mixed systems are studied as they approach thermal equilibrium where detailed balances holds (vanishing chemical potential difference between fuel and waste).

In addition, there may be production and (irreversible) degradation processes that break the (strict) conservation of the cycling components. To dissect the roles of conserving and non-conserving processes, any reaction kinetics can be split up into the respective functional terms. For a general two-component reaction–diffusion system, such a splitting would take the form

[TABLE]

where $\gamma$ describes the ratio of the timescale of mass-conserving reactions encoded in $f$ and the timescale of production/degradation processes encoded in the source terms $s_{1,2}$ Note:GammaInterpolation . If $\gamma$ is small, one can study the mass-conserving subsystem (limit $\gamma\to 0$) using the theory presented here. In particular the capacity for pattern formation of the mass-conserving subsystem can be directly read off from the shape of the $f$-nullcline. The remaining effect of the source terms in this limit is to set the total average mass $\bar{n}$, which is no longer a control parameter when $\gamma\neq 0$, but is controlled by the parameters in $s_{1,2}$ instead. Geometrically, this corresponds to setting the center of mass of a pattern in the $(m,c)$-phase plane.

A decomposition into mass-conserving and non-conserving terms of the form Eq. (53) can be carried out for systems with any number of components. (Note that this should be done before model reductions like non-dimensionalization and variable eliminations that affect the mass-conserving terms.) Linear stability analysis can then be used to determine whether the mass-conserving subsystem is able to form patterns by setting $\gamma=0$. This defines two classes of systems, (i) those where pattern-forming instability is preserved in the mass-conserving case $\gamma=0$, and (ii) those that inherently depend on production and degradation (i.e. don’t exhibit lateral instability for $\gamma=0$). Put differently, class (i) systems do have a pattern-forming, mass-conserving core, while class (ii) systems don’t. Of course, the distinction between these two classes will be limited if $\gamma$ is too large such that production-degradation dynamics dominate over the mass-conserving core.

The mass-conserving ‘core’ of systems (or regimes) in class (i) can be analyzed using the local equilibria theory. An immediate conclusion is that their lateral instability is driven by a mass-redistribution cascade. This insight may inform model reduction of many-component models as we have discussed in the previous subsection. Furthermore, one can then study how the ‘core patterns’ are qualitatively and quantitatively modified by $0\neq\gamma\ll 1$. Importantly, $\gamma\neq 0$ might be a singular perturbation, i.e. weakly broken mass conservation might cause a qualitative change in the dynamics, such as interrupted coarsening and oscillations; see example ($i$) below. Importantly, for $\gamma\ll 1$ these qualitatively new phenomena will play out at large length- and timescales, whereas the behavior at short scales is still determined by the mass conserving core; see e.g. Ref. Brauns et al. (2020a). In the language of bifurcation theory, the mass-conserving ‘core’ would take the role of an “organizing center” from which the various dynamical regimes of the system “unfold” (see e.g. section 30 in Ref. Arnold (1988) and chapter 3 in Ref. Murdock (2003)). Let us provide examples for the two classes of systems defined above.

Example for class (i): Brusselator model. — Let us exemplify the reduction to a mass-conserving core for a classical model, the Brusselator Prigogine and Lefever (1968). This early, qualitative model for the BZ reaction can be written in the form Eq. (53) where $f(m,c)=m^{2}c-m$, $s_{1}=k_{\mathrm{in}}-k_{\mathrm{out}}m$, $s_{2}=0$. The nullcline of $f$ has two segments, $m=0$ and $c^{*}(m)=1/m$, which is a singular case of an N-shaped nullcline. Thus, it is immediately clear that pattern formation of the Brusselator is driven by a mass-conserving core and that its elementary stationary patterns are mesa patterns when $D_{c}\gtrsim D_{m}$ and peak patterns when $D_{c}\gg D_{m}$. Slow production and degradation lead to interrupted coarsening and splitting of the mesa patterns at length scales that depend on $\gamma$ as has been studied using singular perturbation methods in the limit $\gamma\ll 1$ Kolokolnikov et al. (2006, 2007); McKay and Kolokolnikov (2012); Tzou et al. (2013). Local equilibria theory, in particular the phase-portrait analysis of spatially extended systems that it facilitates, provides a new, intuitive approach that explains the physics underlying interrupted coarsening and mesa splitting Brauns et al. (2020a). Moreover, in the oscillatory regime of the Brusselator, the limit cycle oscillations can be constructed as relaxation oscillations on the basis of the nullcline of the mass-conserving core ($f$-nullcline) in the limit $\gamma\ll 1$ Kuwamura and Izuhara (2017). The oscillation period depends on $\gamma$ and diverges in the limit $\gamma\rightarrow 0$.

Example for class (ii): Gierer–Meinhardt model. — An example for a two-component system without a mass-consering core (class (ii)) is the “Gierer–Meinhardt” model Gierer and Meinhardt (1972). This model describes an “activator” which enhances its own production and the production of an “inhibitor” which impedes the “activator’s” production. Both “activator” and “inhibitor” are degraded at constant rates. Written in the form Eq. (53), using $\{a,h\}$ for the “activator” and “inhibitor” concentrations instead of $\{m,c\}$, one has $f=0$, $s_{1}=k_{a}^{+}+k_{\mathrm{fb}}a^{2}/h-k_{a}^{-}a$, $s_{2}=k_{h}^{+}+k_{\mathrm{fb}}a^{2}-k_{h}^{-}h$. Clearly, this system does not possess a mass-conserving core capable of pattern formation. Generally, systems in class (ii) are those where production (from a reservoir or substrate) and irreversible degradation are the dominant processes; for instance, during tissue patterning (morphogenesis) Manukyan et al. (2017); Watanabe and Kondo (2015); Kondo and Asai (1995) or cell-cell signaling in bacterial colonies Karig et al. (2018).

The role of (nearly) conserved quantities in classical pattern-forming systems. — Many chemical systems contain a mass-conserving ‘core’ of components that rapidly cycle between different states and are produced and degraded only on a much slower timescale. An example is the cycling of bromine between a molecular form and a bromide ion in the BZ reaction. It is an interesting question whether this core alone is able to produce some non-trivial behavior like oscillations or patterns (assuming that the chemical fuel driving the cycling of the core components is abundant, as e.g. ATP in protein-based pattern formation). Systems where this is the case can then be analyzed using local equilibria theory (potentially extended to account for slow production and degradation, see e.g. Brauns et al. (2020a)). Above, we showed that this is true for the Brusselator, which is a conceptual model for the BZ reaction. We hypothesize that a similar approach might/will also work for more detailed models such as the so-called FKN mechanism Field et al. (1972), whose mass-conserving core presumably contains more than two components.

Such a program might eventually lead to a perspective that unifies “classical” pattern-forming systems such as the BZ reaction and the more recently discovered biological systems, including the Min system Loose et al. (2008); Halatek and Frey (2018); Brauns et al. (2020b), intracellular actin waves Khamviwath et al. (2013); Bernitt et al. (2017), and Rho excitability Bement et al. (2015); Graessl et al. (2017); Michaux et al. (2018); Tan et al. (2020).

VIII.4.4 Beyond reaction–diffusion systems

In the systems discussed so far (excluding Sec. VIII.3), energy is fed into the system via the reaction kinetics alone, while the spatial transport process, diffusion, is passive. Such systems are part of a broader class of so-called active systems where energy is fed in on the microscopic scale. This comprises systems where the transport processes are active, driven e.g. by molecular motors such as myosin. Examples include systems with active flows generated by cortical contractions Bois et al. (2011); Kumar et al. (2014), as well as actively deforming visco- and poroelastic media Weber et al. (2018); Radszuweit et al. (2013). On a more conceptual level, several non-equilibrium generalization of the Cahn–Hilliard equation have been studied recently Zwicker et al. (2015); Tjhung et al. (2018); Weber et al. (2019); Saha et al. (2020). The phenomena exhibited by the above systems include (micro-)phase separation and more complex phenomena like waves and turbulence. Another broad sub-class of active systems are self-propelled particles that exhibit a huge variety of collective phenomena, including MIPS (Cates and Tailleur (2015), see Sec. VIII.3 above) and flocking Chaté (2020); Huber et al. (2018); Denk and Frey (2020).

What all these systems have in common with MCRD systems is the presence of conserved quantities that serve both as macroscopic variables and as local parameters of the microscopic dynamics. We therefore expect that the ideas on mass redistribution and local equilibria put forward in this article can be broadly applied to understand emergent behavior in these systems. As a concrete example, the self-organized interplay between total density and emergent orientational order (polar or nematic) was investigated recently for self-propelled particles with microscopic alignment interactions that are continuously tunable between polar and nematic symmetry Denk and Frey (2020). Here, the local particle density controls the emergent orientational order, i.e. induces local symmetry breaking. In turn, the orientational order leads to particle currents that redistribute the particle density. Strikingly, this interplay explains the coexistence of different macroscopic structures, such as polar flocks and nematic lanes, and the continual interconversion between them as recently observed in experiments and agent-based numerical simulations Huber et al. (2018).

Acknowledgements.

We would like to thank B. Lohner who was involved at early, preliminary stages of this project, and M. Wigbers and B. Polovnikov for critical reading of the manuscript and for providing valuable feedback. E.F. acknowledges financial support from the Deutsche Forschungsgemeinschaft (DFG) via project P03 within TRR174 (‘Spatiotemporal Dynamics of Bacterial Cells’). F.B. acknowledges financial support by the DFG Research Training Group GRK2062 (‘Molecular Principles of Synthetic Biology’).

F.B. and J.H. contributed equally to this work.

Appendix A Models used for illustration and numerical studies

To visualize the our findings on pattern formation in 2C-MCRD systems we use two prototypical reaction terms, $f(m,c)$, that exhibit two distinct nullcline shapes. Both variants effectively model attachment–detachment dynamics as used to describe cell-polarization systems:

(i) The reaction kinetics used in Ref. Mori et al. (2008), to conceptually describe cell-polarization based on autocatalytic recruitment (Michaelis-Menten kinetics with Hill coefficient 2) and linear detachment:

[TABLE]

We can non-dimensionalize by expressing time in units of $k_{\text{off}}^{-1}$ and densities in units of $K_{\text{d}}$. Furthermore, for specificity we set the (non-dimensional) feedback rate $k_{\text{fb}}/k_{\text{off}}=:\hat{k}_{\text{fb}}=1$, leaving only $k:=k_{\text{on}}/k_{\text{off}}$ as free parameter in the non-dimensional reaction term:

[TABLE]

Fig. 16 shows the nullcline in the $(m,c)$-phase plane of the reaction kinetics Eq. (54*) for $k=0.07$, together with a numerically determined stationary pattern (steady-state solution to Eq. (1)) and the local equilibria (spatial profile on the right).

(ii) Dynamics due to attachment together with linear self-recruitment and enzyme driven detachment (described by first order Michaelis–Menten kinetics):

[TABLE]

We non-dimensionalize by expressing time in units of the attachment rate $k_{\text{on}}$ and densities in units of the dissociation constant $K_{\text{d}}$ of the detachment kinetics. The two remaining parameters are the (non-dimensional) feedback rate $\hat{k}_{\text{fb}}:=K_{\text{d}}k_{\text{fb}}/k_{\text{on}}$ and detachment rate $\hat{k}_{\text{off}}:=k_{\text{off}}/(k_{\text{on}}K_{\text{d}})$:

[TABLE]

where we suppressed the hats. Fig. 17 shows a typical reactive nullcline for the reaction term Eq. (55*) together with a stationary peak pattern and the local equilibria that scaffold it.

Length scale.—

For convenience, we do not specify a unit of length in the domain size $L$ and the diffusion constants $D_{m,c}$. In an intracellular context a typical size would be $L\sim$10\text{,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}, and typical diffusion constants $D_{m}\sim 0.01{-}$0.1\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}^{2}\,\mathrm{s}^{-1} on the membrane and $D_{c}\sim$10\text{,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}^{2},\mathrm{s}^{-1}$$ in the cytosol. Rescaling to different spatial dimensions is straightforward.

Nullcline shape predicts the pattern type.—

Important for the distinction between peak-forming vs. mesa-forming systems is the behavior of the reactive nullcline for large $n$. If it approaches the $m$-axis monotonically for large $n$, then peak patterns are will form in a large range of $\bar{n}$ for $D_{m}/D_{c}\ll 1$ (see Fig. 23). Otherwise, mesa patterns are typical while peak/trough patterns only form in narrow regimes at the edges of the range of pattern existence $[n^{\infty}_{-},n^{\infty}_{+}]$. For attachment–detachment kinetics, one can study the nullcline behavior for large $n$ by comparing the largest powers in the denominator and numerator of the functional form $c^{*}(m)=m\,d(m)/a(m)$ of the nullcline. For the reaction term Eq. (54*), one obtains $c^{*}(m\to\infty)\to m$, i.e. typically mesa patterns. For the reaction kinetics Eq. (55*), one has $c^{*}(m\to\infty)\to 0$, favoring peak patterns for $D_{c}\gg D_{m}$.

Appendix B Numerical simulations

The reaction–diffusion dynamics were simulated on a domain with no-flux boundaries using the numerical PDE-solver routine NDSolve[] provided by Mathematica (see Supp. File “PDE-solver_minimal-setup.nb” for an example setup). In the videos we show the density distribution in phase space and the real space profile $m(x,t)$ together with local equilibria.

Appendix C Linear stability analysis

This section provides the technical details of linear stability analysis.

C.1 Canonical linear stability analysis

Linear stability analysis of a reaction–diffusion system is performed by expanding a spatial perturbation in the eigenbasis of the diffusion operator (Laplacian) in the geometry of the system. In a line geometry with reflective boundary conditions at $x=0,L$, the eigenfunctions of the Laplacian are the discrete Fourier modes $\cos(k\pi x/L)$ with $k\in\mathbb{N}$. Linearization of the dynamics of a mass-conserving two-component system around a homogeneous steady state $(m^{*},c^{*})$ yields the linear dynamics

[TABLE]

with the Jacobian

[TABLE]

where we use the abbreviations $q_{k}=k\pi/L$ for the (discrete) wavenumbers and f_{m,c}=\partial_{m,c}f\big{|}_{(m^{*},c^{*})} for the linearized kinetics at the homogeneous steady state. The eigenvalues of the Jacobian yield the growth rates $\sigma^{(i)}_{q}$ of the respective eigenmodes such that a perturbation in the spatial eigenfunction $\cos(qx)$ evolves in time as

[TABLE]

with the eigenvectors $\mathbf{e}^{(i)}_{q}$. For a given initial condition (perturbation), the coefficients $A^{(i)}_{q}$ are determined by projecting initial condition the onto the eigenbasis $\mathbf{e}^{(i)}_{q}\cos(qx)$.

To calculate the eigenvalues of the Jacobian, we use that the eigenvalues of a $2\times 2$-matrix can be expressed in terms of its trace $\tau$ and determinant $\delta$:

[TABLE]

where the indices $\{1,2\}$ correspond to $\{-,+\}$ on the RHS.

The trace and determinant of the Jacobian $J(q)$ can be written as

[TABLE]

where we used the expression Eq. (28) for $q_{\mathrm{max}}$.

From Eq. (60a), it follows that for a locally stable hom. steady state ($\sigma_{\text{loc}}<0$) the trace $\tau_{q}$ is negative for all $q$. Hence, the only way to get lateral instability is a negative determinant $\delta_{q}<0$. This implies that eigenvalues with positive real part must be purely real, since the term under the square root in Eq. (59) is positive. This means that there cannot be oscillatory lateral instability for a locally stable hom. steady state. Moreover, the instability condition $\delta_{q}<0$ immediately yields the band of unstable modes $[0,q_{\mathrm{max}}]$.

Figure 18 shows the two branches of eigenvalues, $\sigma^{(1,2)}_{q}$, for a laterally unstable case. In the limit $q\to 0$, the first branch connects to the eigenvalue of a well-mixed system $\sigma^{(1)}_{0}=\sigma_{\text{loc}}$ (cf. Sec. III.1). The corresponding eigenvector lies in the reactive phase space for $q=0$, and hence fulfills mass conservation.

The second branch, $\sigma^{(2)}_{q}$, smoothly approaches zero in the limit $q\,{\to}\,0$. The eigenvector $\mathbf{e}^{(2)}_{0}$ corresponding to the marginal eigenvalue $\sigma^{(2)}_{0}=0$ points along the reactive nullcline. It represents a perturbation that changes the total density and thus shifts the reactive equilibrium. For the stability of a closed, well-mixed system (i.e. the stability against homogeneous perturbations) such a perturbation is not relevant since it breaks mass conservation. For $q\neq 0$, the perturbation is spatially inhomogeneous, and therefore redistributes mass in the system. This mass redistribution shifts the local equilibria. That the eigenvector points along the reactive nullcline reflects the fact that , the concentrations are slaved to the local equilibria in the long wavelength limit. As one goes towards shorter wavelengths (i.e. larger $q$) the eigenvector begins to deviate from being tangential to the nullcline (see Supplementary Material of Ref. Halatek and Frey (2018)). In particular, the eigenvalue of the marginal mode $q_{\mathrm{max}}$ points along the flux-balance subspace, $\mathbf{e}^{(2)}_{q_{\mathrm{max}}}\propto(1,-D_{m}/D_{c})^{\text{T}}$.

We found that the band of unstable modes for the 2C-MCRD system always extends down to long wavelength ($q\to 0$), a situation called “type II” instability in the Cross–Hohenberg classification scheme Cross and Hohenberg (1993). In systems with more components and/or multiple conserved species, this is no longer true—the band of unstable modes can be bound away from zero (“type I” in Cross–Hohenberg scheme), see e.g. Ref. Halatek and Frey (2018).

C.2 Approximation close to the onset of lateral instability

The eigenvalues of a $2\times 2$-matrix, Eq. 59, to leading order in $\delta/\tau\ll 1$ are given by

[TABLE]

We will now use this approximation for the 2C-MCRD dynamics where the trace and determinant of the Jacobian for a mode $q$ are given by Eqs. (60). A straightforward calculation shows that $\delta_{q}$ reaches its extremum at $q^{*}=q_{\mathrm{max}}/\sqrt{2}$, with a minimal value of $\delta_{q^{*}}=-\tilde{f}_{m}^{2}D_{c}/(4D_{m})$. Furthermore, the trace $\tau_{q}\approx\sigma_{\text{loc}}$ is approximated by the local eigenvalue for $\tilde{f}_{m}\ll\sigma_{\text{loc}}$. Hence, the above approximation is valid in the vicinity of lateral instability onset ($\tilde{f}_{m}\approx 0$), far away from the local instability onset (which takes place at $\sigma_{\text{loc}}=0$). The dispersion relation then reads

[TABLE]

for $\tilde{f}_{m}^{2}D_{c}/D_{m}\ll|\sigma_{\text{loc}}|$. The first branch, $\sigma_{q}^{(1)}$, simply represents relaxation to local equilibrium. The laterally unstable branch, $\sigma_{q}^{(2)}$, shown in Fig. 19a, is the identical to the dispersion relation of the Cahn–Hilliard equation Cahn (1961) (and the more general class of Model B dynamics). We can rewrite Eq. (62b) as

[TABLE]

The two terms of this expression reflect the shifting of local equilibria (dominating at large wavelengths) driving the instability for $\partial_{n}\eta^{*}|_{\bar{n}}<0$ and the competition of local reactive flow towards equilibrium and membrane diffusion that re-stabilizes the system on short length scales. (To recast Eq. (62b) as (62b*), one uses the relation $\partial_{n}c^{*}(n)=\partial_{m}c^{*}(m)/[1+\partial_{m}c^{*}(m)]$.)

Figure 19 gives an overview of the various regimes of mass-redistribution instability and their interrelation. A typical dispersion relation (solid green line) deep in the laterally unstable regime (i.e. far from onset) is shown in Fig. 19b, together with various limiting cases (see caption for details). One can distinguish diffusion- and reaction-limited regimes. The former (i.e. long wavelength limit) was studied in detail in the main text. In the next subsection, we briefly analyze the two limits of fast cytosol diffusion and vanishing membrane diffusion.

C.3 Limits in the diffusion constants

Fast cytosol diffusion.—

In the limit $D_{c}\to\infty$, the dispersion relation approaches the function $f_{m}-D_{m}q^{2}$ for wavenumbers $q\gg|\sigma_{\text{loc}}/D_{c}|$; see dashed gray line in Fig. 19b). This shape of the dispersion relation reflects the fact that the growth rate of instability is limited by the rate of chemical relaxation to the shifting local equilibria ($\partial_{n}c^{*}\cdot\sigma_{\text{loc}}=f_{m}$) and counteracted by membrane diffusion on short scales ($-D_{m}q^{2}$).

Somewhat deceptively, this shape of the dispersion relation resembles that of a locally unstable system for $q>0$, even though the system is locally stable. In the strict limit, $D_{c}\to\infty$, the dispersion relation becomes discontinuous at $q=0$, because the zero eigenvalue at $q=0$ corresponding to the homogeneous perturbation breaking mass conservation is always present. (Recall that for the stability of a closed system against homogeneous perturbations, this mode is not relevant.)

It is important to keep in mind that the lateral instability is always driven by cytosolic mass redistribution. From this perspective the strict asymptotic limit $D_{c}\to\infty$ (and the common approximation to treat the cytosol as well mixed) is pathological as it masks the core dynamics underlying lateral instability in MCRD systems: the formation of gradients in the fast diffusing component(s) and the ensuing diffusive fluxes that redistribute total density. Furthermore, in systems with more than two components, the subtle interplay of multiple fast diffusing components might play an important role for pattern formation Denk et al. (2018). Such aspects would be missed if fast diffusing components are assumed to be well-mixed at all times (which corresponds to setting their diffusion constants to infinity).

Vanishing membrane diffusion.—

In the limit $D_{m}\to 0$, the band of unstable modes extends to arbitrarily small wavelengths (i.e. $q_{\mathrm{max}}\to\infty$); see purple, dash-dotted line in Fig. 19. For short wavelengths, the growth rate is reaction limited ($\sigma_{q}^{(2)}\approx f_{m}$) because cytosol diffusion is fast. For long wavelengths, the growth rate is diffusion limited, $\sigma_{q}^{(2)}\approx D_{c}q^{2}\,\partial_{n}c^{*}$.

Appendix D Remarks on the local quasi-steady state approximation

The local quasi-steady state approximation (LQSSA), i.e. slaving of the chemical concentrations to the local (reactive) equilibria, becomes exact when the timescales of diffusion and local reactions are separated. Specifically, let us scale the reaction terms in the reaction–diffusion dynamics Eq. (1) by $\varepsilon^{-1}$:

[TABLE]

In the limit $\varepsilon\to 0$, relaxation to local equilibria becomes arbitrarily fast compared to diffusive redistribution—the concentrations will be at a local quasi-steady state (a stable local equilibrium). The characteristic spatial scale(s) of the dynamics and of stationary patterns are given by a balance of reaction and diffusion. In particular, in Sec. V.1, we learned that the interface width is determined by $\pi(D_{m}/\tilde{f}_{m})^{1/2}$. Under the scaling of reaction rates by $\varepsilon^{-1}$, this length scale will go to zero as $\propto\varepsilon^{1/2}$. Hence, in the LQSSA, there is no “microscopic” length scale. This behavior is characteristic for a singular perturbation problem where some physics is lost when the small parameter $\varepsilon$ is set to zero Jackson (1989); Ward (2006). A rigorous analysis of Eqs. (63) could be performed in terms of singular perturbation theory. To lowest order in $\varepsilon$, any series of jumps (sharp interfaces) between two plateaus $n_{\pm}$ that fulfill $\eta^{*}(n_{+})=\eta^{*}(n_{-})$, such that $\eta^{*}(\widetilde{n}(x))=\eta_{0}$ is constant in space, is a valid steady state of Eq. (23). However, to be consistent with Eq. (63) in the limit $\varepsilon\to 0$, the FBS-position must be $\eta_{0}=\eta_{0}^{\infty}$ as determined by total-turnover balance Eq. (35). In addition, the given total density $\bar{n}$ constrains the spatial average $\langle\widetilde{n}(x)\rangle_{[0,L]}=L_{+}n_{+}+L_{-}n_{-}=\bar{n}$, where $L_{\pm}$ are the aggregate lengths of the high and low density regions.

One may compare the LQSSA to the approach used to analyze limit cycle attractors of relaxation oscillators. There the N-shaped nullcline allows an analytic construction of the limit cycle in the asymptotic timescale separation limit. Moreover, even without the timescale separation the qualitative phase-space structure that underlies the oscillations can be deduced from the nullcline shapes. Note that treating such a timescale separation in a mathematically rigorous way requires singular perturbation theory (see Appendix D).

The LQSSA can also be understood as a closure relation. In that picture, $n(x,t)$ corresponds to a “coarse-grained order parameter” with the microscopically correct dynamics given by Eq. (21). This equation is not closed, because $m(x,t)$ and $c(x,t)$ are not known. Equation. (22) is a closure for Eq. (21) at the price of losing the “microscopic” length scale. One could try construct higher order closures that also take into account deviations from the local equilibria owing to diffusion on short length scales Weyer et al. (shed).

In a forthcoming publication, we choose a less technical way to qualitatively illustrate the elementary pattern formation dynamics: we consider the dynamics not on a continuous domain but in two diffusively coupled compartments. In this ‘coarse grained’ setting, the LQSSA is well-posed because a “microscopic” length scale as been imposed externally by the discretization into two-compartments.

Appendix E Geometric construction of bifurcations

In Sec. V.4 and Sec. VII we describe how the bifurcations diagrams of stationary pattern can be constructed geometrically using the reactive nullcline and the flux-balance subspace. We implemented this procedure in Mathematica (see Supp. File “flux-balance-construction.nb” to find quantitative bifurcation structures. As an illustrative example, we present the results for the reaction kinetics Eq. (54*); see Figs. 20 and 22. Figure 20a shows the shape of the reactive nullcline for a range of the kinetic rate parameter $k$ (non-dimensional attachment rate). For $k>k_{\text{saddle}}$, the nullcline is monotonic, such that pattern formation is impossible. (Recall that we set the non-dimensional feedback rate to 1, so $k$ effectively describes the relative strength of basal attachment vs. feedback due to recruitment.) At $k=k_{\text{saddle}}$ a section of negative slope emerges on the nullcline, giving rise to lateral instability for $D_{c}/D_{m}\to\infty$. Further lowering $k$ increases the range of negative nullcline slope and increases the maximal negative nullcline slope (thus decreasing $D_{c}^{\text{min}}$, cf. Eq. (40)). Figure 20b shows the regimes of lateral instability and pattern existence for $D_{c}=10D_{m}$. At $k=k_{\text{cusp}}$ the maximal negative nullcline slope becomes $-1$, indicating a cusp bifurcation of the reactive equilibria. From this cusp point (black dot on the red nullcline in (a)), a regime of bistability emerges, section of unstable equilibria shown as dashed line in (a)). The locally bistable regime (shaded in gray in (b) is delimited by two saddle-node bifurcations (SN) which emerge from the cusp point, shown as black dot in (b). In the locally bistable regime, there exist unstable stationary patterns for $D_{c}<D_{m}$. These patterns can be constructed in the same way as stable stationary patterns for $D_{c}>D_{m}$. Their range of existence for $D_{c}=0.9D_{m}$ is shaded in red in (b), delimited by a dashed red line.

Figure 20c shows the geometrically constructed $(\bar{n},D_{c})$-bifurcation diagram for $k=0.045$, i.e. for a bistable nullcline (corresponding to the schematic bifurcation diagram shown in Fig. 10 in the main text). The $(\bar{n},D_{c})$-bifurcation diagram for a monostable nullcline is shown in Fig. 22 in Appendix F, where we also show the bifurcation structure for finite domain size, $L$, obtained by numerical continuation.

Appendix F Numerical continuation of stationary patterns

To calculate steady states and their bifurcation structures for systems with finite size, we use a standard numerical continuation scheme (pseudo-arclength continuation, see e.g. chapter 4 in Ref. Fiddler (2002)). The stationarity condition Eq. (14a) was discretized using finite differences, yielding a set of equations for the concentrations at the grid points. These equations, together with the flux-balance subspace Eq. (10) and the constraint of average total density Eq. (1c), are used to numerically determine the stationary patterns and their bifurcations (in the Mathematica software). To continue the fold bifurcations of stationary patterns, we use a bordered matrix method Fiddler (2002).

To determine the stability of the stationary patterns, we use a finite difference discretization of the reaction–diffusion dynamics (1) linearized around the steady state. The resulting eigenvalue problem is solved with Mathematica. The eigenvalue with the largest real part (“dominant eigenvalue”) determines the pattern stability (see Fig. 25, which is discussed in Appendix I).

F.0.1 Bifurcation structure for $\bar{n}$

Figure 6 in the main text shows the $\bar{n}$-bifurcation structure of stationary patterns determined by numerical continuation for the reaction kinetics (54*). In Fig 6, the pattern amplitude is plotted against $\bar{n}$. For the same bifurcation structure, Fig. 21, shows additional plots of the maximum and minimum concentrations in (a) and the FBS-position, $\eta_{0}$, in (b). For mesa patterns $\widetilde{m}(0)$ and $\widetilde{m}(L)$ are the plateau concentrations and therefore slaved to $m_{\pm}^{\infty}=m_{\pm}(\eta_{0}^{\infty})$, while the FBS-position $\eta_{0}^{\infty}$ is almost constant. At the boundaries of the range where patterns exist (limited by $[n_{-}^{\infty},n_{+}^{\infty}]$, for $L\to\infty$), the mesa patterns undergo fold bifurcations where they meet the branches of unstable peak/trough patterns (dashed lines) that emanate from the homogeneous steady state (black line, dash-dotted in the regime of lateral instability). In both plots, the prediction from the analytic approximation of (unstable) peak/trough patterns (see Appendix G) is shown as red solid lines.

F.0.2 Two-parameter $(\bar{n},D_{c})$-bifurcation diagram

Figure 22 shows the two-parameter $(\bar{n},D_{c})$-bifurcation diagram for a monostable reactive nullcline (corresponding to the schematic diagram in Fig. 14a; the fixed parameters are the same as in Fig. 6). The fold-bifurcation lines (solid blue lines) of stationary patterns at finite domain size were obtained by numerical continuation. The respective bifurcation lines in the infinite system size limit, $n_{\pm}^{\infty}$ (solid gray lines), were geometrically constructed (cf. Appendix E. The laterally unstable regime, bounded by the dash-dotted green line, was determined by linear stability analysis. On the right, a blow-up of the region around the critical point $(\bar{n}_{\text{inf}},D_{c}^{min})$. The finite sized system,

The tip of the laterally unstable regime in a finite-sized system is shifted upwards by an amount $\sim L^{-2}$ because of the stability condition $\tilde{f}_{m}=\pi^{2}D_{m}/L^{2}$. Close to critical point, the patterns emerge in a supercritical pitchfork bifurcation (dashed green line). The points where the onset becomes sub-critical are marked by in red disks. At these points the two lines of fold bifurcations of stationary patterns originate. The sub-critical lateral instability bifurcation is shown as green dash-dotted line.

Appendix G Approximation of peak/trough patterns

In the following, we detail the construction of approximate peak/trough patterns that was briefly introduced in Sec. V.3. For specificity, we present the construction for peak patterns—generalization to trough patterns is straight forward.

In our analysis of pattern types (Sec. V.3), we characterized peak patterns as composed of an interface region at a system boundary connected to a plateau region (cf. Fig. 5b). We have also characterized interfaces by linearization around the inflection point in Sec. V.1 which yielded a sinusoidal interface shape with a width $\ell_{\text{int}}(\eta_{0})$ (cf.  Eq. (31)). We now construct a peak pattern by piecing together such an (approximate) interface at the left domain boundary and a plateau at $m_{-}(\eta_{0})$ in the reminder of the system:

[TABLE]

Within this approximation, the pattern inflection point is always at $x_{0}=\ell_{\text{int}}/2$. To match the interface to the plateau continuously at $x=\ell_{\text{int}}$, the amplitude must be $A=m_{-}-m_{0}$. The plateau $m_{-}(\eta_{0})$ and inflection point $m_{0}(\eta_{0})$ are geometrically determined. To close the approximation, one has to find the FBS-position consistent with the given average total density $\bar{n}$ (to fulfill the constraint Eq. (1c)). From Eq. (64), one obtains the approximate total density average

[TABLE]

This relation can be inverted to obtain a relation $\eta_{0}(\bar{n},L)$ for peak patterns with a density-profile approximated by Eq. (64).

Peak/trough patterns are encountered in two contexts. First, stable peak patterns are typical for reaction kinetics that exhibit a strongly asymmetric nullcline shape, e.g. Eq. (55*), when $D_{m}\ll D_{c}$ (see Fig. 17 for a typical peak pattern). Secondly, unstable peak/trough patterns form the unstable branches that connect the subcritical Turing bifurcation with the stable pattern branch (see Fig. 6). These unstable peak/trough patterns play the role of ‘transition’ states since they lie on the separatrix that separates the basins of attraction of the homogeneous steady state and the stationary pattern in the multistable regimes.

For both scenarios, we compared the analytic approximation Eq. (64), where $\eta_{0}$ is determined via Eq. (65), with numerical numerical continuation of the stationary patterns.

The approximation of unstable peak/trough patterns for the reaction kinetics Eq. (54*) is shown in the bifurcation structure Fig. 21.

Figure 23 shows the $\bar{n}$-bifurcation diagram of stationary patterns for the reaction kinetics Eq. (55*). There is a large regime of peak patterns where the pattern amplitude keeps increasing with average total density $\bar{n}$. Peak patterns transition to mesa patterns when peak saturates in the third FBS-NC intersection point $m_{+}$. The amplitude of mesa patterns is almost independent of $\bar{n}$, because a change of total density merely shifts the interface position of mesa patterns (compare Fig. 6). Ultimately, mesa patterns transition to trough patterns which then undergo a fold bifurcation where they meet with the unstable branch of trough patterns that emerges from the homogeneous steady state. The asymmetry of the reactive nullcline (Fig. 17a) is reflected by the asymmetry of the bifurcation structure. The dot-dashed red lines in Fig. 23 show the analytic approximation for pattern amplitude (a), the FBS-position (b) and pattern profiles (c). The approximation of peak patterns becomes less accurate as the average total density increases and ultimately breaks down at the transition to mesa patterns (around $\bar{n}\approx 60$ in Fig. 23). The approximation of trough patterns is less accurate because the trough saturates more “abruptly” in $m_{-}$: Recall that the approximation underlying Eq. (65) is a linearization of the interface region around the inflection point (cf. Sec. V.1). This approximation breaks down in regions of high nullcline curvature, indicative of high nonlinearities. Interestingly, even though the pattern profile is not well approximated for troughs (see Fig. 23c), the estimate for the FBS-position $\eta_{0}$ is close to the true value (see inset in Fig. 23c), indicating that the relevant physics (total turnover balance) is still captured.

Appendix H Stimulus-induced pattern formation

In Sec. VI, we argued that to trigger pattern formation from a laterally stable homogenous steady state, a perturbation (stimulus) must induce a (self-)sustained laterally unstable region. Based on this intuition we provided a simple geometric heuristic for the perturbation threshold: a perturbation of the membrane concentration profile must be such that concentrations in a spatial region are pushed beyond the laterally unstable part of the nullcline in phase space. Hence, the intersection point of the line $c=c^{*}(\bar{n})$ with the laterally unstable section of the nullcline, provides an estimate $m_{\text{th}}(\bar{n})$ for the threshold that the membrane perturbation has to exceed in a spatial region (see Fig. 24b, cf. Fig. 7). This criterion does not take into account the spatial shape of a perturbation, but only its characteristics in phase space. We tested how robust the estimate is against different spatial profiles using numerical simulations of the reaction–diffusion dynamics Eqs. (1) with the reaction term Eq. (54*). We consider prototypical perturbations with a ‘step-like’ profile (see Fig. 24a):

[TABLE]

For the perturbation to conserve the global average total density, we must set $b=a\frac{L-w}{w}$. The ‘step-like’ profiles therefore form a two-parameter family of perturbations with the shape parameters $w$ (width of the region where density is increased) and $a$ (density removed uniformly from the rest of the system). Because the concentration may not drop below zero, only perturbations with $a<m^{*}(\bar{n})$ are physically sensible. The (heuristic) threshold in phase space $m_{\text{th}}(\bar{n})$ (see Fig. 24b) is exceeded in the (high density) region $x>L-w$ when

[TABLE]

Note that the threshold $m_{\text{th}}(\bar{n})$ in phase space is a function of the average total density. We tested various total average densities $\bar{n}$ across the multistable regime $n_{-}^{\infty}<\bar{n}<n_{\text{lat}}^{-}$, and varied the ‘shape parameters’ of the perturbation—amplitude $a$ and width $w$—throughout their respective maximal ranges: $0<a<m^{*}(\bar{n})$ and $0<w<L$. Figure 24 shows that there is good agreement between the geometrically estimated threshold and the actual basins of attraction determined by numerical simulation.

Appendix I Stability of stationary patterns

In our analysis of stationary patterns, we have touched the question of stability of these patterns only peripherally in Sec. VII.2. Coarsening, i.e. the instability of multi-peak / multi-mesa patterns, in 2C-MCRD systems has been studied before both numerically Otsuji et al. (2007); Chiou et al. (2018) and semi-analytically Ishihara et al. (2007) for specific choices for the reaction kinetics $f(m,c)$. For specific reaction terms $f(m,c)$ that allow a mapping of the reaction–diffusion dynamics to an effective gradient dynamics, stability of patterns can be analyzed with the help of the effective free energy that is minimized by the stationary pattern Goh et al. (2011); Morita and Ogawa (2010). In the broader class of two-component systems without conserved total density, stability of stationary patterns has been subject to numerous mathematical studies, see e.g. Refs. Nishiura and Fujii (1987); Kolokolnikov et al. (2005); McKay and Kolokolnikov (2012).

Instead of the technical tools typically employed there, we choose a more heuristic approach here, building on the physical intuition we have gained throughout this work. We restrict our analysis to the case of mesa patterns with a small interface width compared to the system size ($\ell_{\text{int}}\ll L$).

Our starting point is the insight that the stationary pattern $(\widetilde{m}(x),\widetilde{c}\kern 0.74997pt(x))$ is embedded in a flux-balance subspace, Eq. (10), whose position $\eta_{0}^{\infty}$ is determined by total turnover balance, Eq. (17). We hence write the stationary pattern as a pair $(\widetilde{m}(x),\eta_{0}^{\infty})$, where only $\widetilde{m}(x)$ depends on $x$. Next, recall that the pattern itself is scaffolded by local equilibria. In particular, the plateaus are slaved to the plateau scaffolds $m_{\pm}(\eta_{0}^{\infty})$, and the pattern inflection point is determined by $m_{0}(\eta_{0}^{\infty})$. Following a perturbation $\delta m(x,t)$ of the stationary pattern profile, the plateaus will quickly return to their stable local equilibria, the plateau scaffolds $m_{\pm}$. On the other hand, a perturbation of the mass-redistribution potential (“FBS-position”) $\delta\eta$, will not only shift the plateau scaffolds, but also cause an imbalance of total reactive turnover. This imbalance drives the dynamics of $\eta(x,t)$ and thus determines the stability of the pattern, as we will see in the following.

The dynamics of $\eta(x,t)=c(x,t)+m(x,t)D_{m}/D_{c}$ follows straightforwardly from the reaction–diffusion dynamics Eq. (1) and read:

[TABLE]

In linearization around a stationary pattern $(\widetilde{m}(x),\eta_{0}^{\infty})$, we have

[TABLE]

where the membrane perturbation $\delta m=\delta m(x,t)$ is governed by the linearization of the reaction–diffusion dynamics Eq. (1a). The intuition is that $\delta m(x,t)$ quickly relaxes to the scaffold of local equilibria. We therefore focus on the dynamics of $\delta\eta(x,t)$, which will affect the scaffold itself by shifting local equilibria, in particular the plateau scaffolds $m_{\pm}(\eta)$.

Reactive turnover balance is primarily determined in the interfacial region (cf. Fig. 2) around the pattern inflection point $x_{0}$. We therefore focus on the interface region to learn how an imbalance of total reactive turnover affects the perturbation of the mass-redistribution potential (“FBS-shift”) $\delta\eta(x,t)$. To that end, we use that the gradient of the membrane profile, $\partial_{x}\widetilde{m}(x)$, is negligible in the plateaus and whereas it peaks at the inflection point $x_{0}$. We hence multiply Eq. (I) by $\partial_{x}\widetilde{m}(x)$ and integrate over the whole domain $[0,L]$ to obtain

[TABLE]

We neglect contributions $\mathcal{O}(\delta m(x,t))$ that correspond to perturbation along the direction of the FBS and quickly relax to the scaffold on the timescale $|\sigma_{\text{loc}}|^{-1}$, fast compared to the contribution by the first term in Eq. (I). Furthermore, because mass-redistribution $\eta(x,t)$ quickly homogenizes in the (small) interface region, we can neglect the second term $~{}\partial_{x}^{2}\delta\eta(x_{0},t)$.

Fig. 25 shows a comparison of the heuristic estimates of perturbation growth/decay rate based on Eq. (43) to numerically determined dominant eigenvalues $\sigma_{\text{max}}(D_{c})$ (linear stability analysis of stationary patterns determined by numerical continuation). The dominant eigenvalue $\sigma_{\text{max}}(D_{c})$ crosses over from system size independent growth rate of perturbations for $D_{c}<D_{m}$ to system size dependent decay of perturbations for $D_{c}>D_{m}$. The instability is well estimated by turnover imbalance term in Eq. (43) (dashed red line in the inset in (b)), while the rate at which perturbations decay in the stable regime ($D_{c}>D_{m}$) is limited by diffusive transport $\sigma_{\text{diff}}\sim D_{c}/L^{2}$ between the far ends of the system. Reactive timescales will become limiting in the stable regime only when $\sigma_{\text{diff}}\approx\sigma_{\text{loc}}$, that is, for fast enough cytosolic diffusion or a small system.

Appendix J Weakly nonlinear analysis

Our goal is to find the stationary pattern in the vicinity of the onset of lateral instability (Turing bifurcation). To that end, we expand the stationary state $(\widetilde{m}(x),\eta_{0})$ in harmonic functions (eigenmodes of the Laplace operator under no-flux boundary conditions):

[TABLE]

where $\eta^{*}=D_{m}/D_{c}m^{*}+c^{*}$ is the FBS position of the homogenous steady state. Mass conservation necessitates $\delta m_{0}+\delta\eta_{0}-D_{m}/D_{c}\delta m_{0}=0$, hence

[TABLE]

We plug the ansatz Eq. (71a) into the stationarity condition, Eq. (14a), Taylor expand $\tilde{f}$, and project onto the zeroth harmonic (this is most conveniently done in CAS software, e.g Wolfram Mathematica)

[TABLE]

and and onto the second harmonic

[TABLE]

Solving for $\delta m_{0}$ and $\delta m_{2}$ we get

[TABLE]

These equations describe how asymmetry of the nullcline shape (and thereby reactive turnover) influences the pattern profile. For $\tilde{f}_{mm}=0$ the turnovers on either side of the inflection point grow symmetrically (with opposite sign) as the amplitude $\delta m_{1}$ of the pattern increases. This symmetry of the turnovers only occurs at the inflection point of the nullcline where its shape is point-symmetric. Away from the inflection point of the nullcline, its shape is no longer symmetric around the steady state $(m^{*},c^{*})$, and, respectively, the turnovers grow asymmetrically in the two halves of the system. This asymmetry creates an imbalance of the total turnover that is compensated by two effects: (i) the flux-balance subspace shifts (i.e. $\eta_{0}$ will deviate from $\eta^{*}$), and (ii) the pattern itself becomes asymmetric as the second harmonic will have a non-zero amplitude $\delta m_{2}$. Together, these two effects compensate the asymmetry of the turnovers, such that total turnover balance is reached.

The amplitude $\delta m_{1}$ of the stationary pattern is obtained by projecting (14a) onto the first harmonic $\cos(\pi x/L)$

[TABLE]

where

[TABLE]

and

[TABLE]

Since there is no second order term in Eq. (75), patterns always originate in a pitchfork bifurcation. At the bifurcation point, the first order coefficient vanishes ($F_{1}=0$). The system is laterally unstable if $F_{1}$ is positive (cf. Eq. (28)). Hence, only if $F_{3}<0$, the third order coefficient can saturate the pattern amplitude (supercritical bifurcation). For $F_{3}>0$ the bifurcation is subcritical.

The third order coefficient can be simplified further: A simple calculation shows that

[TABLE]

and therefore

[TABLE]

where $\tilde{\partial}_{m}=\partial_{m}-(D_{m}/D_{c})\partial_{c}$ is the derivative along the flux-balance subspace. With that, the second summand in the brackets in $F_{3}$ (cf. Eq. (77)) can be rewritten and we obtain

[TABLE]

We further rewrite the denominator of the first summand in the brackets as $F_{1}-3D_{m}\pi^{2}/L^{2}$, and use that $F_{1}$ vanishes at the bifurcation point (i.e. is small in its vicinity). We thus have

[TABLE]

Note that in weakly nonlinear approximation Eq. (46) in the main text, we have pulled the $\mathcal{O}(F_{1})$ out from $F_{3}$ to simplify notation. The role and physical interpretation of the three terms in $F_{3}$, as written in the form Eq. (81), are discussed in the main text in Sec. VII.4.

Appendix K Nullcline curvature approximation

In the following, we show that the nullcline curvature $\kappa$ can be approximated by

[TABLE]

in the vicinity of the Turing bifurcation (onset of lateral instability; recall the slope criterion Eq. (25)). Start by rewriting the second derivative $\tilde{f}_{mm}$ in terms of derivatives of $f$:

[TABLE]

In the vicinity of the Turing bifurcation, we have $-D_{m}/D_{c}\approx-f_{m}/f_{c}$ (from the slope criterion for lateral instability, Eq. (25)), so we obtain

[TABLE]

Comparing to the formula for the curvature $\kappa$ of an implicitly determined curve $f(m,c)=0$

[TABLE]

one sees that the numerator of the curvature formula is identical to the term in the square brackets in Eq. (84). Thus, by combining Eqs. (84) and (85), we obtain the approximation Eq. (82).

Appendix L Topological equivalence of 2C-MCRD systems shear banding in complex fluids

Analogies between the shear banding in complex fluids, phase separation near thermal equilibrium, and reaction–diffusion systems have been drawn before on mathematical grounds, i.e. using similarity of the equations used describe these phenomena Radulescu et al. (1999); Fardin et al. (2012). The phase space analysis of MCRD pattern formation presented in the present work establishes a connection to shear banding on the basis of phase space geometry. More specifically, these two phenomena can be regarded as topologically equivalent, i.e. they can be understood in terms of equivalent geometric objects in phase space. We shall briefly outline this connection in the following.

Complex fluids can exhibit a non-monotonic constitutive relationship $\Sigma(\dot{\gamma})$ between the total stress $\Sigma$ and the (homogeneous) strain rate $\dot{\gamma}$ Olmsted (2008); Divoux et al. (2016). When the total stress decreases upon an increase in strain rate, $\partial_{\dot{\gamma}}\Sigma<0$, a mechanical instability results, which leads to a separation of a sheared fluid into “shear bands” with different viscosities and strain rates, which coexist at a common total stress. The “common total stress” construction on the constitutive curve $\Sigma(\dot{\gamma})$ employed to analyze this phenomenon (in a one-dimensional system) is analogous to our flux-balance construction on the reactive nullcline for 2C-MCRD systems (see Sec. III.3), by means of a mapping $(n,\eta)\leftrightarrow(\dot{\gamma},\Sigma)$ between the phase space variables. The average strain rate $\bar{\dot{\gamma}}=L^{-1}\int_{0}^{L}\mathrm{d}x\,\dot{\gamma}$ is analogous to the average total density $\bar{n}$. The constitutive $\Sigma(\dot{\gamma})$ curve is analogous to the reactive nullcline. The steady state conditions are spatially uniform total stress and flux-balance (spatially uniform mass-redistribution potential; see III.2) respectively. The selection of the common total stress generally depends on the details of the model Dhont and Briels (2008), in particular on stress diffusion Fardin et al. (2015). For simple models, it can be pictured similarly to a Maxwell-construction (cf. total turnover balance in a 2C-MCRD system illustrated in Fig. 2). Furthermore, momentum propagation due to stress gradients in complex fluids is equivalent to mass redistribution due to concentration gradients in MCRD systems. Accordingly, the low Reynolds number limit is analogous to the $D_{c}\to\infty$ limit in the 2C-MCRD system.

Taking these analogies together, we conclude that these physically distinct phenomena are topologically equivalent and can be studied with similar phase-space geometric tools. Such a connection might benefit both fields as more involved scenarios are investigated, for instance coupling to additional degrees of freedom: For models of complex fluids, additional spatial dimensions Fielding (2010), and coupling to internal structure of the fluid Fielding (2007); Turcio et al. (2018), can lead to a variety of intricate spatiotemporal patterns; for mass-conserving reaction–diffusion systems, additional components or additional conserved species can equally lead to a broad range of phenomena Murray and Sourjik (2017); Halatek and Frey (2018); Alonso and Bär (2010). Studying such systems in terms of local equilibria theory offers an exciting new perspective for future research.