Diffusive coupling of two well-mixed compartments elucidates elementary principles of protein-based pattern formation

Spatial organization of proteins in cells is important for many biological functions. In general, the nonlinear, spatially coupled models for protein-pattern formation are only accessible to numerical simulations, which has limited insight into the general underlying principles. To overcome this limitation, we adopt the setting of two diffusively coupled, well-mixed compartments that represents the elementary feature of any pattern---an interface. For intracellular systems, the total numbers of proteins are conserved on the relevant timescale of pattern formation. Thus, the essential dynamics is the redistribution of the globally conserved mass densities between the two compartments. We present a phase-portrait analysis in the phase-space of the redistributed masses that provides insights on the physical mechanisms underlying pattern formation. We demonstrate this approach for several paradigmatic model systems. In particular, we show that the pole-to-pole Min oscillations in Escherichia coli are relaxation oscillations of the MinD polarity orientation. This reveals a close relation between cell polarity oscillatory patterns in cells. Critically, our findings suggest that the design principles of intracellular pattern formation are found in characteristic features in these phase portraits (nullclines and fixed points) rather than the topology of the protein-interaction networks.


I. INTRODUCTION
The spatial intracellular organization of proteins by reactions (protein-protein interactions) and diffusion has received growing attention in recent years; for recent reviews see Refs. [1][2][3][4][5][6][7][8]. Gaining intuition and theoretical insight into the spatiotemporal protein dynamics remains challenging owing to the complexity arising from the spatial coupling and nonlinear reaction terms. Therefore, insights often remain restricted to specific mathematical models. A systematic understanding is hard to achieve, in particular if there are multiple protein species with several conformational states involved (complex interaction network). Thus, finding the elementary principles underpinning protein-based pattern formation still remains a largely open question.
To simplify the analysis on a technical level, systems of two diffusively coupled, well-mixed compartments (also called 'boxes', 'reactors', 'cells', or 'patches') have been widely used in earlier literature. In fact Turing himself used the setting of diffusively coupled compartments (called "cells") in his pioneering work to show that diffusion can destabilize otherwise stable reactions, thus leading to spatial pattern formation [9]. Physically, the twocompartment setting represents the elementary feature of any pattern-an interface connecting a low density region to a high density region. In the context of intracellular pattern formation, the two compartments typically represent the polar zones of rod-shaped cells, such as E. coli bacteria (see Fig. 1a), M. xanthus bacteria [10,11], and fission yeast (S. pombe) [12,13]. * frey@lmu.de In a broader context, two-compartment systems also have been realized in experiments, using diffusively coupled CSTRs (continuously stirred tank reactors) [14] and recently using nanometer scale microfluidic devices [15,16]. Furthermore, in population dynamics, they are known as "two-patch systems" and have been used to study toe role of spatial coupling and patterning in ecology, see e.g. [17][18][19].
In this manuscript, we focus on protein-based pattern formation in cells. A key property of such intracellular pattern formation is that the total number of proteins is conserved on the relevant time scale of pattern formation [6,[20][21][22][23]. Recent works [24,25] suggest that (diffusive) mass redistribution is the key physical process driving pattern formation in mass-conserving reaction-diffusion systems. Based on this insight, a framework termed local equilibria theory has been developed [25]. The basic idea of this framework is to consider the system as decomposed into (notional) compartments, small enough to be effectively well-mixed. Within each compartment, the reactive dynamics conserves the mass(es). The reactive equilibria (steady states) of the reactions within an isolated compartment, controlled by these local masses, serve as proxies for the local dynamics. Diffusive coupling of the compartments redistributes masses between them. In turn, the changing local masses shift the local reactive equilibria and potentially change their stability. Thinking about reaction-diffusion systems in terms of this interplay between mass-redistribution and shifting local equilibria has proven a powerful approach to study their complex nonlinear dynamics [24][25][26][27][28].
Here we adopt the two-compartment setting and show how this way of thinking can be made explicit in the form of simple graphical constructions and a phase portrait analysis in the phase space of the redistributed masses. This will enable us to gain insights on the physical mechanisms underlying pattern formation that would otherwise remain hidden. Importantly, and in contrast to previous works [10,11,[29][30][31], we do not assume the fast diffusing (cytosolic) components to be well mixed. In other words, we explicitly allow for cytosolic gradients between the two compartments. As we will see later, this is important understand the physical mechanisms underlying pattern formation. In particular, it is key to explain the pole-to-pole oscillations of Min proteins in E. coli. a. Motivation. Let us present the main motivation for this work using the pole-to-pole oscillations of Min proteins in E. coli as a concrete example without going into technical details (which will be presented below). Put briefly, the pole-to-pole oscillations are driven by two types of proteins, MinD and MinE, which cycle between membrane-bound and cytosolic states and interact with each other on the membrane (Fig. 1b), while the total masses of MinD and MinE (n D and n E ) remain conserved. A key insight from previous works is that spatial redistribution of such globally conserved masses constitutes the essential degrees of freedom of mass-conserving reaction diffusion systems [24]. Indeed, mapping the Min system to the two-compartment setting and tuning the diffusive exchange rates to a slow time scale retains the qualitative features of the pole-to-pole oscillations (Fig. 1c,d). On the slow time scale, only the masses in the two compartments n (1,2) D,E remain as dynamic variables. Because of mass conservation, the average masses n D,E remain constant. Defining the redistributed masses, ∆n D and ∆n E , via n (1,2) D,E =n D,E ± ∆n D,E , we can visualize the dynamics in the two-dimensional (∆n D , ∆n E )phase space, Fig. 1e, where we plot the flow field and its nullclines. Along the nullclines the rate of mass-exchange between the compartments vanishes. We hence refer to them as mass-redistribution nullclines. The phase portrait shows the characteristics of relaxation oscillations. In this paper, we show that the Min pole-to-pole oscillations are indeed spatial relaxation oscillations of the MinD polarity orientation.
This example shows how important qualitative features of mass-conserving reaction-diffusion (MCRD) systems can be obtained from a phase portrait analysis in the phase space of the redistributed masses. In the following, we show how this phase portrait can be constructed systematically, starting from the reactiondiffusion equations. We show what determines the structure of the phase space flow and derive a simple geometric relation between the mass-redistribution nullclines and the reactive nullclines of the local reaction kinetics.
b. Structure of the paper. To introduce the basic elements of our analysis, we first study MCRD systems with two components, e.g. the membrane-bound and cytosolic state of a single protein species (see Sec. II). We then generalize the nullcline-based approach to systematically derive the phase portrait of the Min system of E. coli shown in Fig. 1e. This construction then allows us to study the role of diffusive mass redistribution of MinD and MinE for the formation of Min-protein patterns. Finally, we apply the same approach to two other paradigmatic model systems: PAR polarity of C. elegans and Cdc42 polarity of budding yeast. Comparing the different nullcline geometries of these systems allows one to classify their pattern-forming mechanisms (see Sec. IV). Such a nullcline-based classification provides intuition for the role of various elements in the biochemical network. Moreover, it might guide model building and serve as a first step of analysis for systems that are biochemically not as well characterized as the aforementioned examples. In the Conclusions, Sec. V, we discuss important implications of our work, both specific to the Min system and in a broader context, and give an outlook to promising future research directions.

II. TWO-COMPONENT MCRD SYSTEMS
Two-component MCRD systems have been previously used as conceptual models for cell polarity [20,29,32,33]. In this section, we apply local equilibria theory [24,25] to these systems in the two-compartment setting. In this simplified setting, the formulation of local equilibria theory is technically simpler than in spatially continuous systems [34]. Importantly, the approach developed below for two-component MCRD systems can be generalized to systems with more components and more conserved masses such as those studied in Sections III and IV, where the new approach yields novel insights.
Let us denote the concentrations of the two components in compartment i ∈ {1, 2} by u i = (m i , c i ), where m i and c i are the concentration of membrane-bound and cytosolic proteins, respectively. The reaction kinetics f = (f, −f ) within each compartment account for the attachment and detachment to and from the membrane. Importantly, they conserve the local total density (mass) n i = m i + c i in each of the two compartments individually. Mass is transferred between the compartments by a diffusive exchange process that acts to even out concentration differences. Denoting the diffusive exchange rates in the matrix D = diag(D m , D c ), we have the coupled compartment dynamics in vector notation (1) Since both the local reactions and the diffusive exchange are mass conserving, the average total densitȳ n = (n 1 + n 2 )/2 is a constant of motion. In Appendix A, it is shown how the (diffusive) exchange rates D m,c can be related to the diffusion constants D m,c in a spatially continuous system, in such a way that the linearized dynamics of Eq. (1) near a homogeneous steady state is identical to the linearized dynamics of a single Fourier mode ∼ cos(πx/L) in the spatially continuous system on the interval [0, L] with no-flux boundary conditions. For patterns with large amplitudes, nonlinearities lead to mode coupling in a spatially continuous system. This is not captured by the two-component system which only describes the dynamics at a single length scale. Nonetheless, one can gain a good qualitative understanding of the full nonlinear pattern formation process, including the termination of the pattern-forming instability in a stationary pattern.
A. Setting the stage: phase-space geometry of two-component MCRD systems In the following, we present the key concepts of local equilibria theory in the two-compartment setting. Because of mass conservation, only the mass density difference with respect to the mean ∆n := (n 1 − n 2 )/2 is a dynamic variable, whilen is a control parameter. Thus, we can rewrite the local masses as n 1,2 (t) =n ± ∆n(t). Adding the equations for ∂ t m 1 and ∂ t c 1 (Eq. (1) yields ∂ t n 1 , and thus where ∆u = (∆m, ∆c) := u 1 − u 2 and we used that ∂ tn = 0. Observe that the reaction terms cancel because they conserve the mass in each compartment individually. Thus, the dynamics of the total density is solely determined by concentration differences in m and c between the two compartments. These concentration differences approximate the gradients in the spatially continuous system. To understand how these concentration differences are governed by the reaction kinetics, consider the (m, c)phase plane of the reaction kinetics (see Fig. 2a). While this phase plane is two-dimensional, mass conservation also implies that reactive flow (f, −f ) in each compartment i is constrained to a respective linear subspace m i + c i = n i . We term these subspaces the local phase spaces of each compartment [22,25]. Here, and in the following, the term local always refers to the properties of a single (notionally isolated) compartment. Correspondingly, we define as local reactive equilibrium the point within the local phase space where the reaction kinetics are balanced, i.e. where the reactive flow vanishes (f = 0): Geometrically, the local equilibria are the intersection points between the local phase spaces and the reactive nullcline (see Fig. 2) [35]. These local equilibria determine the steady state (reactive equilibrium) in each compartment that would be reached if, given the local masses n 1 and n 2 , the two compartments were isolated, i.e. if the diffusive exchange between the compartments was shut off. Thus, the local equilibria serve as proxies for the local reactive flow within each of the compartments (red arrows in Fig. 2a). Diffusive coupling between the compartments redistributes mass between the compartments. This is reflected in the shifting of the local phase spaces in the (m, c)-phase plane, as indicated by the purple arrows in Fig. 2a,b). As a result, the local reaction kinetics change since the local equilibria move in the (m, c)-phase plane. In the following we will elucidate this interplay between diffusive mass-redistribution and shifting local equilibria in the most elementary form.

B. Limit of slow mass exchange
To separate the roles of local reactions and diffusive mass redistribution, we consider a situation where the latter occurs on a much slower time scale than the former . The concentrations (mi, ci) in the two compartments are marked by blue dots, labelled 1 and 2, respectively. The local phase spaces corresponding to the masses in the two compartments n1,2 =n ± ∆n are shown as gray lines. Gray arrows indicate the reactive flow towards the reactive nullcline f = 0 (solid black line). Black dots mark the local equilibria (intersection points between reactive nullcline and local phase spaces) and red arrows indicate the reactive flow towards these local equilibria. (b) When diffusion is set to a slower time scale, the local concentrations adiabatically follow the reactive nullcline. Thus, the only remaining degree of freedom is the mass difference ∆n, whose dynamics is governed by the concentration differences ∆m * (∆n) and ∆c * (∆n) (see Eq. (2)). [36]. In this limit, the cytosolic and membrane concentrations in each compartment adiabatically follow the local equilibria that depend on the local masses n i , as encoded by the shape of the reactive nullcline in the (m, c)-phase plane (see Fig. 2b). We can, therefore, approximate the densities by their respective equilibrium values We term this the local quasi-steady state approximation (LQSSA). The dynamics of the mass difference ∆n is then governed by a closed equation with the shorthand notation for the concentration differences between the two compartments: ∆u * (∆n) := u * (n + ∆n) − u * (n − ∆n).
In this approximation, the roles of local reactive dynamics and diffusive mass exchange are clearly separated. The concentrations only change if the local phase spaces shift due to mass redistribution. In turn, the mass fluxes from one compartment to the other are determined by the concentration gradients ∆u * (∆n), weighted by the respective exchange rates D m,c . This nonlinear feedback between shifting equilibria and mass redistribution is the basic mechanism underlying pattern formation in massconserving reaction diffusion systems [24,25]. Importantly, the role of the reaction kinetics is fully encoded in the shape of the reactive nullcline, i.e. the functional dependence of the reactive equilibrium concentrations u * (n) on the total density n. The local masses n i within each compartment play the role of control variables [24] that determine the position of the local phase spaces (and thus the position of the reactive equilibria) within the (m, c)-phase plane. At the same time the local masses are also dynamic variables that change by means of diffusive mass redistribution between the compartments. Accordingly, we refer to the phase space of the redistributed masses as control space. In the two-component MCRD system, the only control variable is the mass difference ∆n, such that the control space is one-dimensional.
Typically, diffusion on the membrane is orders of magnitude slower than in the cytosol, D m D c such that its contribution to mass redistribution can be neglected; see e.g. Refs. [37,38]. Hence, to simplify the following analysis, we neglect the slow membrane diffusion (i.e. we set D m = 0), such that ∂ t ∆n(t) = −D c ∆c * (∆n) = −D c c * (n + ∆n) − c * (n − ∆n) .
Generalization to account for the effect of membrane diffusion is straightforward by changing variables from c to the 'mass-redistribution potential' η := c + (D m /D c )m [39]. Equation (7), has a simple geometric interpretation as shown in Fig. 3b,c. The term in the brackets in Eq. (7) expresses the difference between the nullcline (solid, black line) and its mirror image (dashed gray line) reflected at the pointn. Depending on the nullcline slope atn, the resulting dynamics ∂ t ∆n, indicated by the blue arrows, is qualitatively different. For a positive slope, ∂ n c * (n)|n > 0, following a small perturbation from the "homogeneous" state ∆n = 0 the system returns to the ∆n = 0; see Fig. 3b. In contrast, for a negative slope, ∂ n c * (n)|n < 0, the homogeneous state is unstable; see Fig. 3c. This criterion for a lateral instabil- . (c) Forn > ncrit, there are two additional intersection points between c * (n+∆n) and c * (n − ∆n), corresponding polarized steady states. The flow in control space is directed away from the homogeneous state ∆n = 0, which is therefore unstable, and drives the system towards one of the stable polarized states.
ity (instability against spatially inhomogeneous perturbations) was previously derived in Ref. [25] for spatially continuous systems. The physical mechanism for this mass-redistribution instability is that the reactive equilibrium shifts to lower concentration of the fast diffusing (cytosolic) component, c * (n), when the total density n is increased, and vice versa. Hence, a small perturbation δn results in a gradient ∆c that transports mass from the compartment with lower mass to the compartment with higher mass. This amplification mechanism drives the instability. The growth of the mass difference ∆n will stop once the cytosolic gradient ∆c * (∆n) vanishes, i.e. when the cytosolic concentration is the same in both compartments, c * (n + ∆n) = c * (n − ∆n). Thus, stationary states can be determined graphically as the intersection points of the nullcline c * (n) with its own mirror image, mirrored atn, as illustrated in Fig. 3c. The intersection point at ∆n = 0 always exists by construction, and corresponds to the homogeneous steady state. The two intersection points at ∆n = 0 represent polarized steady states.
In summary, we have shown how one can graphically construct the mass-redistribution dynamics of two-compartment systems with one conserved mass simply based on the reactive nullcline u * (n). In the next section, we will generalize this construction to systems with two conserved masses.

III. TWO-CONSERVED MASSES: THE EXAMPLE OF MIN-PROTEIN OSCILLATIONS
The Min-protein system is a paradigmatic model system for intracellular pattern formation. It was discovered in E. coli, where the pole-to-pole oscillations of the proteins MinD and MinE allow the cell to position its division machinery at midcell [40,41]. This spatial oscillation, i.e. the alternating accumulation of the proteins at the two cell poles is driven by cycling of MinD and MinE between cytosolic and membrane bound states, fuelled by ATP (details described below). Subsequent to its reconstitution in vitro [42], the Min system has been studied in great detail, both experimentally [26,[42][43][44][45][46][47][48][49][50][51][52][53] and theoretically [22,24,26,49,52,[54][55][56][57]. This research has revealed a bewildering zoo of patterns, including traveling waves, standing waves, spatiotemporal chaos, and defect mediated turbulence, observed in different experimental setups (including microfluidic devices [26,46] and vesicles [50,51]). Recent works employing local-equilibria theory to interpret data from numerical simulations and experiments have provided insights on the mechanisms underlying these patterns and their relationships among each other [6,26].
Here, we revisit the comparatively simple pole-to-pole oscillation employing the local-equilibria theory in the two-compartment setting. This offers a fresh perspective on the Min-protein dynamics as it allows us to understand this elementary dynamic pattern in terms of phase space geometry, independently of numerical simulations. In future work, this could serve as a starting point to systematically understand more complex patterns, like "stripe oscillations" (standing waves) in filamentous cells [22,41] and the zoo of patterns found in vitro [26,53,58].
Intuitively, the two-compartment system represents the two cell poles (or cell halves) of the rod-shaped E. coli bacterium, as shown in Fig. 1a (see Appendix B 1 for a systematic reduction starting from the full threedimensional cell geometry). Figure 1c,d shows that the key qualitative features of Min pole-to-pole oscillations are still captured by the two-compartment model (see also Fig. 7 in Appendix B). While this two-compartment model cannot be expected to give a detailed quantitative description of Min oscillations, it has the advantage of informing about the basic underlying mechanism. This complements earlier quantitative studies of the in vivo dynamics [22,57]. Moreover, the two-compartment model serves as a minimal system for an oscillation mode recently reported for an in vitro reconstitution of the Min system in microfluidic devices [26]. There, the oscillations go back and forth between two membrane surfaces through the bulk solution in-between them (see Fig. 8 in Fig. B). The analogy between this in vitro oscillation mode and pole-to-pole oscillations in vivo is further discussed in the conclusions, Sec. III D.
We use a well-established minimal model for the Minprotein interactions that has been shown to successfully reproduce and predict a large range of experimental findings, quantitatively in vivo and qualitatively in vitro [6,22,27,49,56]. For a detailed description of the model, we refer the reader to Refs. [22,27]. In short, the minimal model employs mass-action law kinetics to account for the attachment and detachment of MinD and MinE to and from the membrane and for their interactions there: Membrane-bound MinD amplifies the attachment of further MinD from the cytosol with rate k dD and also recruits MinE from the cytosol with rate k dE to form MinDE complexes on the membrane. In these complexes, MinE stimulates MinD hydrolysis with rate k de , leading to the dissociation of the complex and detachment of both proteins to the cytosol. In the cytosol, MinD undergoes nucleotide exchange from the ADP-bound form to the ATP-bound form, which can then attach to the membrane again.
Mathematically, the above reaction kinetics are described by the system of equations of the form Eq.
where the reaction terms account, respectively, for MinD attachment and selfrecruitment to the membrane, MinE recruitment by MinD, and dissociation of MinDE complexes with subsequent detachment of both proteins to the cytosol. The term λc DD accounts for nucleotide exchange, i.e. conversion from c DD to c DT , in the cytosol. Importantly, these reaction kinetics conserve the total number of MinD and MinE proteins,n D andn E , individually, i.e. there are two globally conserved masses that are redistributed between the two compartments (cell halves) [59]. Numerically integrating the above set of ordinary differential equations using the parameters from Ref. [22] yields pole-to-pole oscillations in good qualitative agreement with the oscillations found in the full threedimensional geometry (see Fig. 7a,b). Importantly, these oscillations persist if diffusive exchange between the compartments is set to a slow time scale compared to the reaction kinetics (see Fig. 7c) In this limit, the concentrations in the two compartments adiabatically follow the equilibrium concentrations that depend on the local masses n D,i , n E,i in the two compartments. Hence, we can again apply the LQSSA, Eq. (4), substituting the concentrations u by the reactive equilibria u * . A discussion of the validity of this approximation and potential generalizations is deferred to the Conclusions, Sec. V.
The reactive equilibria as a function of the masses n D and n E are (for each compartment) determined by (cf. Eq. (3)) where we introduced the total cytosolic MinD concentration c D = c DD + c DT . For each component in the concentration vector u * this defines a surface parametrized by n D and n E , as shown in Figure 4a,b for c * D and c * E (the respective surfaces for the membrane concentrations m * d and m * de are shown in Fig. 10 in Appendix B). We will term these reactive nullcline surfaces. In the following, we show how the dynamics of the local masses n D,i , n E,i can be inferred from these surfaces, analogously to the construction shown in Fig. 3 for two-component MCRD systems.
Because the total number of MinD and MinE proteins are conserved, only the protein masses redistributed between the two polar zones, ∆n D,E (t), are time dependent and the mass densities of MinD and MinE in the right and left polar zone are given by Analogously to the two-component system, we call the redistributed masses ∆n D,E (t) the control variables and the (∆n D , ∆n E )-phase plane the control space. The dynamics in control space are governed by where the concentration gradients (differences between the two polar zones) of the local equilibria are defined as (cf. Eq. (6)) A. From reactive nullcline surfaces to mass-redistribution nullclines To understand the qualitative structure of the controlspace dynamics Eq. (12), we first consider the lines along which there is no mass-redistribution of MinD/E respectively, ∂ t ∆n D,E = 0. We term these mass-redistribution nullclines. Importantly, these are not to be confused with the reactive nullcline (line of reactive equilibria) along which the reactive flow vanishes within a single compartment. As we shall see in Sec. III B, one can neglect the slow membrane diffusion to understand the basic oscillation mechanism of the Min system. We therefore consider this simpler case, D d = D de = 0, first. Equation (12) then reduces to describing how mass redistribution is driven by the gradients in the cytosolic protein densities, which are slaved to the local equilibria. Thus, the mass-redistribution nullclines are simply given by ∆c * D = 0 and ∆c * E = 0. Geometrically, this corresponds to the intersection lines between the reactive nullcline surfaces c * D,E (n D , n E ), and their respective point reflections, reflected at the point (n D ,n E ); see gray surfaces with dashed outlines in Fig. 4a,b. In other words, the shape of reactive nullcline surfaces encodes the essential information about the nonlinear reaction kinetics for the dynamics of the spatially coupled system. This construction is the analog to the construction for the two-component system shown in Fig. 3. In fact, in slices with n E = const, the line c * D (n D ) has the same shape as the nullcline shown in Fig. 3. This "hump shape" gives rise to the N-shaped MinD-redistribution nullcline (∂ t ∆n D = 0, see blue line in Fig. Fig. 4a,c). The two outer branches of this N-shaped nullcline represent polarized MinD states, corresponding to the two outer fixed points in the analogous two-component system construction Fig. 3c. We will make this more concrete below in Sec. III C. Ifn D lies to the left of the crest of c * D (n D , n E ), the resulting MinD-redistribution nullcline is monotonic, analogous to the single fixed point in Fig. 3b. The crest of the c * D surface defined by ∂ nD c * D (n D , n E ) = 0 approximately follows the line n E /n D ≈ k dD /k dE for sufficiently large n E (specifically in the limit n 2 . This relation is found by applying the implicit function theorem to evaluate the derivative ∂ nD c * D using the definition Eq. (10) for the reactive equilibria.
In contrast to the non-trivial MinD-redistribution nullcline, the monotonicity of the surface c * E (n D , n E ) gives rise to a monotonic MinE-redistribution nullcline (red line in Fig. 4b,c) for alln D ,n E . Mass-redistribution potentials.
In passing, let us introduce an alternative formulation of the massredistribution dynamics Eq. (12) that allows one to generalize the graphical construction presented in Fig. 4 to arbitrary values of all diffusion constants (including D d,de > 0). Using the mass-redistribution potentials (cf. Ref. [25] (12) can be written as Since these equations have the same form as Eq. (14), the construction shown in Fig. 4 can be generalized by exchanging c * D,E for η * D,E . The surfaces η * D and η * E can be interpreted as "superpositions" of the local-equilibrium surfaces of the individual components weighted by the respective exchange rates D i . The effect of reaction rates or diffusion constants on the spatial dynamics is encoded in the deformation of these surfaces under variation of these parameters (see Supplementary Movies S1 and S2).

B. Min pole-to-pole oscillations are relaxation oscillations
The nullclines enable one to read off the qualitative structure of the dynamics in the (∆n D , ∆n E )-phase plane [60,61]. Specifically, one immediately recognizes the characteristic scenario of a relaxation oscillator [62]. Recalling that the two outer branches of the N-shaped MinD-redistribution nullcline correspond to polarized MinD states, this shows that Min pole-to-pole oscillations are relaxation oscillations of the MinD-polarity direction driven by mass-redistribution of MinE between the two cell halves.
The limit cycle of relaxation oscillators can be graphically constructed in the limit where the variable with the N-shaped nullcline evolves on a much faster time scale compared to the other variable [61]. In the Min system, this corresponds to setting MinE redistribution to a much slower time scale than MinD redistribution (D D D E ). In this limit, the limit cycle deforms into a "trapezoidal" trajectory; see Fig. 4d. The dynamics slowly follows the N-shaped MinD-redistribution nullcline (polarized MinD states), driven by MinE redistribution, and rapidly switches between the left and right branches at the extrema of this nullcline, driven by MinD redistribution.
In a broader context, the above analysis demonstrates how the reactive nullcline surfaces and their intersection lines-which are the mass-redistribution nullclines-are helpful tools to explore the ability of systems to show nontrivial spatial dynamics without the need to perform a full scale finite element simulation. The shape of the reactive nullcline surfaces and thus the mass-redistribution nullclines are ultimately a consequence of the nonlinear feedback in the reaction kinetics. In the specific case of the Min system, these are the recruitment terms k dD m d c D and k dE m d c D . It is important to recall that the shape of the nullclines resulting from the reaction kinetics, and not the specific reaction kinetics per se, determines the spatial (mass-redistribution) dynamics. Hence, different reaction terms can give rise to same nullcline geometry, and thus the same spatial dynamics. Rather than classifying dynamics based on their reaction network topology, this suggests that a classification might be possible in terms of the shapes of their reactive nullcline surfaces and the resulting mass-redistribution nullclines. We demonstrate this in Sec. IV, where we analyze two further paradigmatic models for intracellular pattern formation.
C. The role of diffusive MinE transport So far, we have neglected membrane diffusion to elucidate the basic Min-oscillation mechanism. We now relax that approximation and first consider the role of MinE membrane diffusion. Using the conservation law m de + c E = n E , Eq. (12b) can be recast as This shows that diffusive transport on the membrane always counteracts cytosolic transport. In particular, if one were to set D E = D de , there would be no MinE massredistribution since Eq. 16 would reduce to ∂ t ∆n E = −D de ∆n E , such that ∆n E would simply relax to the homogeneous state ∆n E = 0. Thus, in control space, the MinE-redistribution nullcline would simply be given by ∆n E = 0, which intersects the N-shaped MinD nullcline at three points, representing the unstable homogeneous steady state in the center and two stable polarized states on the left and right, respectively. Hence, the dynamics would reduce to the one-dimensional control space for MinD redistribution which, corresponding to the scenario shown in Fig. 3. From this Gedankenexperiment, we conclude that the elementary pattern-forming mechanism of the Min system is MinD polarization and does not require spatial redistribution of MinE. The specific function of MinE in MinD polarization is that of a "local catalyst" that provides nonlinear feedback essential in shaping the non-monotonic reactive MinD nullcline c * D (n D ). Thus, while redistribution of MinE is not required for the formation of a polarized MinD pattern, it causes the emergence of oscillations by periodically inducing switching of the MinD polarization direction as we showed in the previous section. Physiologically, D E = D de would actually correspond to a scenario where free MinE remains membrane bound, i.e. MinE would cycle between the MinD-bound and the free conformation on the membrane and c E would then denote the concentration of free MinE. The stationary patterns resulting in this case provide a potential hint for the possible biomolecular features of MinE responsible for the (quasi-)stationary patterns reported in recent experiments using MinE purified with a His-tag at the C-terminus instead of the N-terminus [53]. Compared to his-MinE, MinE-his might have a strong membrane affinity causing free MinE to remain membrane-bound after the dissociation of MinDE complexes. Free MinE on the membrane diffuses much slower than in the cytosol thus suppresses the MinE redistribution that gives rise to dynamic patterns (waves and oscillations). This hypothesis suggests that increasing the MTS strength of MinE might cause a transition from dynamic to quasi-stationary patterns.
To elucidate the role of MinE transport more quantitatively, we now study the transition from stationary to oscillatory patterns as a function of the diffusion constants D E and D de . Varying these diffusion constants results in a deformation of the MinE-redistribution nullcline in control space. Specifically, the shape of the MinD-redistribution nullcline only depends on the difference D E − D de , i.e. the balance of co-polarizing diffusion of free MinE compared to the contra-polarizing diffusion of MinDE complexes. In the relaxation-oscillation limit where MinE-redistribution is much slower than MinD redistribution (D de , D E D D ), the locations of the intersection points between the MinD's and MinE's massredistribution nullclines determine whether the system is oscillator or exhibits stationary polarity (see Fig. 5b). The transition case separating these two regimes is when the MinE-redistribution nullcline intersects the MinDredistribution nullcline at its extrema. In addition, the stability of the homogeneous steady state can be obtained by a linear stability analysis in LQSSA (see Appendix C 2). The resulting "phase diagram" is shown in Fig. 5a.
This phase diagram obtained using LQSSA can be compared to the phase diagram of the full model obtained by numerical simulations (see Fig. 5b). The fact that the topology of the two phase diagrams agrees shows that the reduced dynamics, Eq. (12), accounts for the relevant physics of the in vivo Min system.

D. Concluding remarks on the Min system.
We have shown that dynamics underlying Min poleto-pole oscillations can be reduced to the redistribution of MinD and MinE mass between the two cell poles. A simple geometrical construction yields the qualitative phase space structure of the mass-redistribution dynamics. Specifically, we recovered the paradigmatic Nshaped nullcline that underlies general relaxation oscillations. This systematic reduction immediately allowed us to transfer the knowledge on relaxation oscillations to the Min pole-to-pole oscillations. The outer legs of the N-shaped MinD-redistribution nullcline correspond to oppositely polarized MinD states. MinE redistribution drives cyclic switching between these two states, giv-ing rise to the pole-to-pole oscillations. In the absence of MinE redistribution (achieved by setting D E = D de ), MinD forms stationary polarized patterns instead. This shows that the elementary pattern underlying for poleto-pole oscillations in the in vivo Min system is not oscillatory but generic cell polarity. We conclude that the oscillatory dynamics are not a direct property of the kinetic interaction network, which is the same for the oscillatory and non-oscillatory regime. Instead, oscillations arise as consequence of MinE redistribution "downstream" of MinD polarization. MinE redistribution is not necessary for MinD polarization while MinD redistribution is strictly required. This links pole-to-pole oscillation in the Min system and generic cell polarity and suggests a hierarchy of species in large multi-species multi-component systems. Notably, this also shows that the functional role of MinE for pattern formation cannot be considered to be that of an inhibitor in the sense of the "activatorinhibitor" mechanism [63,64].
The above analysis of the mass-redistribution dynamics elucidates the different roles of MinD and MinE redistribution for Min-protein pattern formation. In Sec. IV, we apply the same reduction approach to two other intracellular systems. This will allow us to compare the underlying pattern-forming mechanisms on the level of their mass-redistribution nullcline geometries.
a. Min oscillations in vitro. Let us emphasize again that the pole-to-pole oscillations emerge due to the diffusive coupling of two compartments, representing the two cell halves. An isolated compartment exhibits only stable, stationary states. In other words, the in vivo Min system is not an "oscillatory medium" of coupled oscillators. Remarkably, this is in stark contrast to the Min-protein pattern dynamics observed in classical in vitro setups with a large cytosolic bulk volume on top of a flat membrane surface. Here, a single (laterally isolated) membrane patch is coupled to an extended cytosolic reservoir, and it is this coupling that gives rise to local oscillations [6,26]. This shows that on a mechanistic level, Min protein patterns in cells are distinct from patterns in reconstituted systems with a large bulk.
In a recent work, the Min system was studied in microfluidic chambers with two flat membrane surfaces separated by a bulk solution [26] (see Fig. 8 in Appendix B). This limits the bulk volume above each membrane patch and thus suppresses the local oscillations for sufficiently low bulk height. Interestingly, for intermediate bulk height, experiments and a theoretical analysis have revealed an oscillation mode that transports mass between the two opposite membrane surfaces through the bulk in-between them. This oscillation is analogous to the in vivo pole-to-pole oscillation where the two opposite membrane patches play the role of the cell poles in vivo. Correspondingly, with regard to the in vitro geometry, the two-compartment system serves as a minimal system to represent single vertical bulk column and the membrane patches at its top and bottom; see   b. Historic note: Oscillations driven by diffusive coupling of two "dead" cells. Intriguingly, the Minoscillation mechanism described above has some parallels to a conceptual model for diffusion-driven oscillations studied by Smale already in 1974 [65]. Smale's motivation, inspired by Turing's pioneering work [9], was to show how two identical reactors that exhibit only a stable stationary state when isolated, start oscillating (in anti-phase) when coupled diffusively. Or, as Smale put it: "One has two dead (mathematically dead) cells interacting by a diffusion process which has a tendency in itself to equalize the concentrations. Yet in interaction, a state continues to pulse indefinitely." As we showed above, the in vivo Min system also has that property.
Remarkably, Smale used a relaxation oscillator as starting point to construct the diffusion driven twocompartment oscillator. In a broader view, this demonstrates how structures in phase space, like fixed points and nullclines, are powerful tools to understand and design nonlinear systems. For instance, they have been used to great success in the study of neuronal dynamics [66] and biochemical oscillators [67,68].

IV. CONTROL SPACE FLOW OF THE PAR AND CDC42 SYSTEMS
The above investigation of the Min system demonstrates that the key characteristics of the spatio-temporal protein dynamics, and the underlying pattern-forming mechanisms, can be inferred from the shapes of the reactive nullcline surfaces. In the following, we use the approach introduced above to two paradigmatic model systems for intracellular self-organization: the PAR system of C. elegans and the Cdc42 system of budding yeast (S. cerevisiae). Starting from previously established mathe- matical models on spatially continuous domains, we follow the same reduction procedure as for the Min system; details on the models, parameter choices and the reduction procedure are described in Appendix D. Put briefly, the spatially continuous dynamics is mapped to the twocompartment setting, and the LQSSA is applied such that only the redistributed masses remain as dynamic variables. The mass-redistribution dynamics can then be analyzed in terms of the reactive nullcline surfaces and the resulting phase portraits as shown in (Fig. 6). This allows us to compare the pattern forming mechanisms underlying these different systems.
a. PAR system. The first division of C. elengans embryos is asymmetric, where the fate of the daughter cells is defined by proteins called aPARs and pPARs that segregate along the long axis of the ellipsoidal cells [4]. A model for the formation of these segregated domains was introduce in Ref. [21], based on the mutual antagonism between cortex-bound A-and pPARs (see Fig. 6a). Here, we adopt this model here, to illustrate the phase portrait structure that is characteristic of the mutualantagonism mechanism. Model details and the param-eters are given in Appendix D. Since, the reaction network (and the parameters used) are symmetric, so are the reactive nullcline surfaces (Fig. 6b). From the resulting mass-redistribution nullclines (Fig. 6c), we can immediately see that the patterns form by segregation into domains where pPAR concentration is high while aPAR concentration is low and vice versa. Notably, the mass-redistribution nullclines do not intersect the lines ∆n A = 0 and ∆n P = 0 away from the origin, indicating that PAR-pattern formation the requires the redistribution of both protein species. Moreover, the topology of the phase portrait is such that oscillations cannot occur. We expect that these qualitative insights generalize to more detailed models for PAR-protein polarity, see e.g. Ref. [69]. b. Cdc42 system. Budding yeast cells divide asymmetrically by budding and growing a daughter cell. The division site is determined by the polarization of GTPbound Cdc42 to a "spot" on the membrane [70]. In wild-type cells, Cdc42 polarization is driven by a mutualrecruitment mechanism that is facilitated by the scaffold protein Bem1. Bem1 is recruited to the membrane by Cdc42-GTP. Membrane-bound Bem1 then recruits Cdc42's GEF, Cdc24, forming Bem1-GEF complexes. In turn, Bem1-GEF complexes recruit Cdc42-GDP from the cytosol and catalyze its conversion to Cdc42-GTP, thus closing the feedback loop. To illustrate the phase portrait structure that is characteristic of this mutual recruitment mechanism, we adopt a simplified form of the detailed, quantitative model introduced in Ref. [71]; see Appendix D. In the simplified model, Bem1-GEF complexes are described as a single species with a membranebound and a cytosolic state (see Fig. 6d). Figure 6e,d shows the reactive nullcline surfaces and the resulting phase portrait of this model. The location of the massredistribution nullcline intersection points, corresponding to stationary polarized states, indicates that Cdc42 and Bem1-GEF complexes co-polarize as expected. Moreover, the N-shaped Cdc42-redistribution nullcline that intersects the line ∆n B = 0 three times, indicating that polarization does not require spatial redistribution of Bem1-GEF complexes. Still, the enzymatic action of Bem1-GEF complexes in the local reaction kinetics is essential for Cdc42 polarization as they provide the nonlinear feedback that shapes the Cdc42-redistribution nullcline. In this sense, Bem1-GEF complexes play an analogous but inverse role to MinE in the Min system. In the physiological case, Bem1-GEF complexes they stabilize polarity by co-polarizing with Cdc42. In the unphysiological case that free Bem1-GEF complexes diffuse slower that membrane-bound ones (D b > D B ), contra-polarization of Bem1-GEF complexes drives cycling switching of Cdc42 polarity. Thus, the Cdc42 system and the Min system can be regarded as two complementary versions of the same mechanism in which the enzymatic function of the "secondary protein" (Bem1-GEF/MinE) is reversed such that its spatial redistribution has opposite effects in the two systems.
The above analysis has a striking implication: On the level of the pattern forming mechanisms, the Cdc42 system is closely related to the Min system, while the PAR system operates based on a fundamentally different mechanism. From the perspective of the phenomenology exhibited for physiological parameters, this is highly surprising since the Cdc42 system and the PAR system exhibit stationary polarity patterns, while the Min system exhibits pole-to-pole oscillations.

V. CONCLUSIONS
Quantitative models of biological systems are typically multi-component multi-species systems with a highdimensional parameter space. It is therefore particularly challenging to find a unifying level of description where the mechanisms underlying different models can be compared.
Here, we presented a reduction method to obtain a phase-portrait representation of mass-conserving pattern forming systems which crystallizes their key qualitative features. This reduction is based on two steps. First, a reduction of the spatially continuous domain to two well-mixed compartments coupled by diffusion. This approximation assumes that the pattern of interest is a single "interface" connecting a high density region to a low density region. This is rather the rule than the exception for protein patterns observed in cells, especially bacterial cells due to their small size. Moreover, such an interface can also be interpreted as the elementary building block of more complex patterns with many interfaces. Second, the local quasi-steady-state approximation (LQSSA), which assumes that the relaxation of the concentrations in the compartments to a reactive equilibrium (local quasi-steady state) is fast compared to slow diffusive mass exchange between the compartments. This approximation is motivated by the insights that the essential degree of freedom is the spatial redistribution of the conserved masses and that the key information about the reaction kinetics is encoded in the dependence of the reactive equilibria on these masses. Limitations and potential extensions of the LQSSA are discussed in the Outlook, Sec. V A.
After these two reduction steps, the only remaining degrees of freedom are the differences in globally conserved masses between the two compartments. In this reduced system, the dynamics of these mass differences can simply be inferred from the reactive nullcline (hyper-)surfaces. Specifically, the intersection lines of reactive nullcline surfaces act as mass-redistribution nullclines in the phase space of the redistributed masses. The mass-redistribution nullclines depend on the diffusion constants and thus inform about the role of massredistribution in the observed phenomena. Thus, they allow a classification of pattern-forming systems, as we demonstrated by comparing the phase portraits of three different protein-pattern forming systems. Attempts to classify pattern-forming systems based on the topology of the protein interaction network face the difficulty that many networks can give rise to similar phenomena, and the same network can produce different phenomena depending on parameters (e.g. stationary and oscillatory patterns in the Min system). In contrast, here we have demonstrated that the geometry of the reactive nullcline surfaces informs on the key qualitative features of the observed dynamics. This suggests that one can identify geometric design principles based on the shape of the reactive nullcline surfaces and the resulting massredistribution nullclines. Such design principles might guide future model building efforts in a similar way as the design principles that have been identified for neural excitability [66] and well-mixed biochemical oscillators [67,68].
The phase-portrait analysis in terms of massredistribution nullclines also shows that not all species need to be redistributed for patterns to form in the first place. One can construct a "core" pattern-forming system, where these species are considered non-diffusible and their kinetics absorbed into effective kinetics of the redistributed species. In the Min system and the Cdc42 system, the (local) enzymatic action of MinE / Bem1-GEF complexes is part of the core pattern-forming mechanism, whereas their cytosolic redistribution is not. Redistribution of MinD / Cdc42 is sufficient for the formation of (stationary) MinD / Cdc42 patterns. Thus, the elementary polarization mechanism is equivalent in the Min system and the Cdc42 system. The difference of these system lies in the effect of the mass redistribution of the "secondary proteins" MinE and Bem1-GEF respectively. In the Min system, redistribution of MinE by cytosolic diffusion system drives cyclic switching of the MinD polarity axis and thus gives rise to pole-to-pole oscillations. In contrast, redistribution of Bem1-GEF complexes stabilizes stationary Cdc42 polarization.
Taken together, the shape of the reactive nullcline surfaces and the resulting mass-redistribution nullclines inform about important qualitative features of a model and thus bridge the gap between nonlinear reaction kinetics and the observed phenomena. In particular, they allows one to disentangle the functional roles of each protein species in the pattern-forming mechanisms.
a. Assuming a well-mixed cytosol misses important physics. The assumption of a well-mixed cytosol is often made a priori, justified by the observation that diffusive transport on cellular scales is fast compared to membrane diffusion (and reaction kinetics); see e.g. [10,11,13,[29][30][31]. This reasoning overlooks that the relative rates of transport can be important if there is more than one protein species diffusing in the cytosol. Or put differently, setting the cytosol concentrations well-mixed neglects that the cytosol gradients of different species can have different amplitudes, which may be mechanistically relevant, even if the cytosol gradients are shallow compared to membrane gradients.
In fact, for the Min system, we find that increasing the diffusion of free MinE eventually always suppresses pattern formation in the Min system (see phase diagram Fig. 5a and Appendix B 4). This shows that the relative rate of cytosolic transport of MinD vs MinE (and, the relative amplitude of the cytosolic gradients, respectively) is important for the dynamics. This shows that one misses important physics if one assumes a well-mixed cytosol a priori.
In general, the time scales of cytosol diffusion-even if fast-and, correspondingly, the relative amplitudes of cytosolic gradients-even if shallow-can be important if there is more than one cytosolic (fast-diffusing) species. Approaches, such as the so called "local-perturbation analysis" [30], that rely on the a priori approximation to treat fast diffusing components as well-mixed, may therefore miss important features of the dynamics.
In passing, we note that explicit cytosol diffusion is also important to account for effects due to cell geometry. This is relevant for the axis selection of polarity patterns in rod-shaped or ellipsoidal cells [52,69,72]. Compartment-based models-although requiring more than two compartments-have also been employed suc-cessfully to study such geometry effects [72].
A. Outlook a. Future applications and generalizations. Going forward, it will be interesting to apply the reduction method and phase-portrait analysis presented here to other model systems, e.g. the oscillatory Cdc42polarization in fission yeast [12,13]. The phase portrait analysis might be particularly helpful to study genuinely nonlinear phenomena like stimulus induced pattern formation and stimulus-induced polarity switching [11] which are not accessible to linear stability analysis.
Potential obvious generalizations of the twocompartment setting are systems with asymmetric exchange rates, and those with heterogeneous compartments (reaction kinetics, bulk-surface ratio, size). Indeed, setting the redistribution of one species to a slow time scale in the models with two conserved masses (e.g. Min system), makes the system heterogeneous from the point of view of the fast species. The heterogeneity is determined by distribution of the slow species between the two compartments and changes on the slow time scale. Concrete application for heterogeneous two-compartment models might be Ran-GTPase driven nuclear transport, where the two compartments represent the cytoplasm and nucleoplasm, with transport between them through pores in nuclear envelope [73][74][75][76]. More broadly, two-compartment systems with asymmetric exchange rates and heterogeneous compartments have been studied in ecology [17,19], where interesting new effects compared to the symmetric case were found.
Another route of generalization is to study more than two coupled compartments. In this case, the phase space of the mass differences becomes high-dimensional and thus impractical for a phase-portrait analysis [72]. Instead, one can plot all local masses into one graph, as was done in Ref. [6]. This way, the spatial information is lost, but one can still gain insight about the role of the control space structure (surface of local equilibria and their stability) for the dynamics of the spatially coupled system.
b. Relation to parameter-space topology. A previous work on reaction-diffusion models for cell polarity has identified generic topological features of their parameter spaces [23]. In the specific case of two-component systems, the origin of these features was recently traced back to the phase space geometry, specifically the shape of the reactive nullcline of pattern forming systems (see Sec. VII in Ref. [25]). Two-compartment systems are a promising setting to generalize these findings to systems with more components and phenomena like pole-to-pole oscillations. Indeed, the way the mass-redistribution nullclines deform due to the variation of parameters (kinetic rates, diffusion constants, average masses) determines the bifurcations in parameter space. Thus, we expect a close relation between the geometry of mass-redistribution nullclines and phase space topology.
c. Relaxing the local quasi-steady state assumption. The analysis presented here relied crucially on the stability of the local equilibria and a time scale separation between reactive relaxation to the local equilibria and diffusive mass redistribution. This allowed us to make the LQSSA Eq. (4). In the absence of this time scale separation, the concentrations will deviate from the local equilibria due to the diffusive flows in the individual components. For two-component systems, this deviation from the local equilibria has only a quantitative effect but does not change the dynamics qualitatively. This is because the local phase spaces are one-dimensional such that the reactive flow is always directed straight towards the local equilibrium (see Fig. 2a). In contrast, in systems with more components, explicitly accounting for the relaxation towards local equilibria may be important to capture the salient features of the full dynamics. One potential approach is to allow for small deviations from the local equilibria along the direction of the slowest decaying eigenvector. Moreover, local equilibria may become unstable, driving the concentrations away from them [24,77]. This qualitative change of the local reaction dynamics can have profound consequences on the dynamics of the spatially extended system, as was studied in detail in Refs. [24,26]. There, it was found that destabilization of the local equilibria gives rise to chaos near the onset of pattern formation.
Even if a systematic reduction in terms of a (generalized) LQSSA is not possible, visualizing the trajectories from full numerical simulations in control space can be a powerful tool to gain insight into the underlying mechanisms [6,26]. The diffusive exchange rates D α can be related to the diffusion constants D α in a spatially continuous system in two alternative ways. First, a finite volume approximation of the Laplace operator on a line with reflective boundary conditions yields Second, we can choose the exchange rates such that the linearization of Eq. (1) for an antisymmetric perturbation u 1,2 = u * ± δu is identical to the linearization of a spatially continuous MCRD system for a Fourier mode ∼ cos qx with q = π/L: The factor 2 in the denominator originates from the linearization of the exchange terms in Eq. (1) for the antisymmetric mode where any perturbation in compartment 1 is balanced by an equal and opposite perturbation in compartment 2. For symmetric perturbations u 1,2 = u * + δu, corresponding to homogeneous perturbations of the continuous system, the exchange term in Eq. (1) cancels. For the exchange rates Eq. (A2), the small amplitude dynamics of antisymmetric perturbations of the two-compartment system exactly represent the linearized dynamics of a single mode q = π/L in the spatially continuous system, and one can use the system size L as a bifurcation parameter to sample the whole dispersion relation σ(q = π/L). The two options above differ by a factor D (LSA) α /D (FV) α = π 2 /8 ≈ 1.23. This can be interpreted as an effective rescaling of the system size L by a factor π/(2 √ 2) ≈ 1.11 due to the finite difference discretization. Throughout this study, we used the exchange rate defined by Eq. (A2).
Appendix B: Min system: Geometry reduction, parameter choice, numerical simulations and phase diagram 1. Reduction from three-dimensional spherocylinder to two-comparmetment system We model the three-dimensional cell geometry as a spherocylinder of length L = 3 µm and radius R = 0.5 µm. The surface of this spherocylinder represents the cell membrane and is the domain for the protein densities m d , m de , while its three-dimensional bulk is the domain of the cytosolic components c DT , c DD and c E . Reactive boundary conditions at the surface account for attachment and detachment of proteins at the membrane. The mathematical implementation of the Min-skeleton model in this three-dimensional bulk-surface coupled setting can be found in [26].
To reduce this geometry to the two-compartment system, we cut the spherocylinder at midplane and assume that the cytosol and membrane in both halves are well mixed. That is, we only account for concentration differences between the two cell halves which serve as a proxy for the concentration gradients along the cell. Moreover, we express the cytosol concentrations in units of surface density,ĉ = ζc, where ζ is the ratio of cytosolic bulk volume to membrane area (short bulk-surface ratio). This allows us to collect all concentrations in a vector that does not mix units. Substituting c →ĉ/ζ, all reaction rates for reaction terms involving a cytosol concentration  [22]. ζ is the bulk-surface ratio that appears because we express cytosol concentrations in units of surface density µm −2 , as explained in the text (Appendix. B 1).

Parameter
Value are rescaled by the bulk-surface ratio:k = k/ζ. In the following, we drop the hats. The bulk-surface ratio of a spherocylinder is given by which, with L ≈ 3 µm and R ≈ 3 µm for E. coli, gives ζ ≈ 0.23. For the in vitro setup using flat microchambers whose top and bottom surfaces are covered by lipid bilayers that mimic the cell membrane, the bulk-surface ratio is simply H/2, where H is the height of the microchamber; see Fig. 8. With respect to this microchamber geometry, the two-compartment system represents a single, laterally isolated cytosol column with two membrane patches at its top and bottom. Only vertical gradients in the cytosol on the scale of the microchamber height are accounted for by the two compartments.

Parameter choice
For the physiological parameters from [22], the densities enter a regime where the reaction kinetics is bistable (i.e. where there are two stable reactive equilibria for given local total densities, see Fig. 9a). This "local bistability" does not change the dynamics of the spatially coupled system qualitatively. However, it complicates the analysis in terms of LQSSA to deal with the branch switching that happens when the dynamics leaves the locally bistable region: Upon passing the saddle-node bifurcations that delimit the bistable region, the concentrations jump to the remaining branch of stable equilibria. To avoid these technical subtleties, we reduce the total densities to values where the local system no longer becomes bistable (inset in Fig. 9a). Because this also increases the minimal domain size for instability, we increase the domain length by a factor 8. The oscillation period increases due to this, but the limit cycle in control space does not change qualitatively (see Fig. 9b,c).

Simulations on 1D domain
In Figure. 5, we compare simulations of the twocompartment system to simulations in a spatially continuous domain (1D line) with no-flux boundary conditions. The dynamics in this domain is given by is the diffusion matrix. (As in the twocompartment setting, the concentrations are measured in units of surface density, µm −2 . To convert the bulk concentrations to volume densities, they must be divided by the bulk-surface ratio ζ.) The reason that we do not perform the simulations in the three-dimensional cell geometry is that we are interested in the role of lateral MinE transport, which we study by tuning the diffusion constants D E and D de . Bulk-surface coupling induces bulk-concentration gradients in the direction normal to the membrane. Those gradients control the flux onto and off the membrane (attachment-detachment dynamics). Hence, changing the cytosol diffusion constants in the bulk-surface coupled 3D system affects both transport and the reaction kinetics. Reducing the system to a 1D line geometry, which effectively amounts to neglecting vertical gradients, allows us to tune the cytosol diffusion constants to study the role of lateral mass transport alone.

No instability for well-mixed cytosol
In the limit D D , D E → ∞, the cytosol is well mixed, i.e. c (1) with f given by Eq. (8), the dynamics is governed by We now perform a linear stability analysis of the homogeneous steady states (m (1) = m (2) = m * , f (u * ) = 0) of these equations and show that they never exhibit a symmetry breaking instability. Because of the parity symmetry, 1 ↔ 2, of the homogeneous steady state, even and Nonlinear reactions (f , red arrows) account for cycling between membrane-bound and cytosolic states (concentrations m and c). Diffusive exchange is indicated by purple arrows. Time traces (center) and phase-space trajectories (right) of the redistributed masses ∆n D,E between the two cell-halves/compartments show good qualitative agreement between the full 3D simulation and the two-compartment system. Importantly, setting the diffusive exchange rates to a much slower time scale (D → εD, here ε = 10 −2 ) does not qualitatively alter the pole-to-pole oscillation (c). Illustration of an in vitro setup using a flat microchamber with two membrane surfaces (gray planes) on top and bottom of a bulk volume [26]. An individual column of that system, comprising two membrane patches and the bulk volume in-between them can be pictured as an analog to the cell geometry, where the two membrane patches correspond to the cell poles. The analogous approximation by two compartments, as shown on the right, is valid as long as the vertical bulk gradient is approximately linear. Comparing to Fig. 1a, the analogy between pole-to-pole oscillations in cells and vertical membrane-to-membrane oscillations in microchambers becomes immediately evident.
odd perturbations are decoupled. Even perturbations correspond to the stability against homogeneous perturbations. Odd perturbations correspond lateral stability, i.e. stability against spatially inhomogeneous perturbations. These are the relevant perturbations for pattern formation. For odd perturbations, mass conservation of MinD and MinE enforces δc DD = −δc DT and δc E = 0. Thus, we obtain the eigenvalue problem with the Jacobian The eigenvalues of J odd are (B7) One immediately sees that only the first eigenvalue, σ 1 , The Jacobian is evaluated at the homogeneous steady state where f (u * ) = 0. In particular, which implies k dD c * DT < k dE c * E for all steady states. Therefore, the necessary condition for instability, Eq. (B8), is never fulfilled. In conclusion, the Min skeleton model with well-mixed cytosol cannot exhibit a lateral instability (instability against spatially inhomogeneous perturbations). This result, derived in the twocompartment setting also holds in spatially continuous domains thanks to the correspondence between these two setting; see Sec. A.
Appendix C: LQSSA 1. General setup and notation Consider a system with N components, u = {u i } i=1..N , governed by local reactions f (u) that conserve M masses, n = {n α } α=1..M . The conserved masses are given in terms of the component vector via n α = s T α u where s α are "stochiometric" vectors fulfilling s T α f = 0. Denoting the diffusive exchange rates by the matrix D = diag{D i }, the dynamics in LQSSA is given by where the slaved concentration gradients ∆u * are defined in terms of the reactive equilibria as ∆u * = u * (n+∆n)− u * (n − ∆n); cf. Eq. (6). The reactive equilibria u * (n) are defined by The factor s T α D determines the diffusive mass flux of species α that results from slaved concentration gradients. We now define the "mass-redistribution potentials" [25] η α := s T α Du, which allows us to write the massredistribution dynamics as ∂ t ∆n α = −∆η * α (∆n). (C3)

Linear stability analysis
For small perturbations δn around the homogeneous steady state ∆n = 0, the dynamics is given by where, in the second line, we introduced the massredistribution Jacobian J αβ := −2∂ n β η * α |n = −2s T α D (∂ n β u * |n).
The eigenvalues of this Jacobian determine the stability of the homogeneous steady state in LQSSA.
Before we continue to calculate the derivatives ∂ n β η * α in terms of the linearized reaction kinetics, let us take moment to interpret the Jacobian J . In the case of one conserved mass n, we have the 1×1 Jacobian J = −2∂ n η * | n . Hence, we recover the nullcline-slope criterion for lateral instability ∂ n η * | n < 0 (cf. Eq. (27) in Ref. [25]). For more than one conserved mass, the entries of J are the slopes of the nullcline (hyper-)surfaces η * (n) along the directions of the conserved masses in parameter space. The eigenvalue problem for J can therefore be interpreted as a generalized slope criterion. To find the nullcline slopes ∂ n β η * α , we take the derivative of the defining equation for the reactive equilibria Eq. (C2) with respect to n α which gives (implicit function theorem) where e α is the unit vector with entry 1 in the αth component. Substituting this in Jacobian yields This can easily be implemented numerically to obtain the Jacobian and calculate its eigenvalues. a. Equivalence to perturbation theory in longwavelength limit. The Jacobian derived above for the two-compartement system in LQSSA can also be obtained by a long-wavelength perturbation theory for linear stability analysis on a continuous domain. To see why this is, consider the Jacobian on a spatially continuous domain where D = diag({D i }) is the diffusion matrix, and q is the wavenumber (i.e. −q 2 are the eigenvalues of the Laplace operator). In the long wavelength limit q → 0, we can find the eigenvalues of J q by solving the degenerate perturbation problem with q 2 as perturbation parameter. We are interested in the eigenvalue branches that emanate from 0 at q = 0, corresponding to the conservation laws. The associated left eigenvectors, (spanning the left nullspace of D u f ) are the "stochiometric vectors" s T α . The right eigenvectors of D u f associated to the eigenvalue 0 are ∂ nα u * . This follows immediately from the defining equation f (u * ) = 0 by taking the derivative w.r.to n α and using that f does not explicitly depend on n α . The first order perturbation of the degenerate 0 eigenvalues is given by eigenvalues of the matrix where we used s T α ∂ n β u * = ∂ n β (s T α u * ) = ∂ n β n α = δ αβ . Substituting the diffusion matrix D by the exchange rate matrix via Eq. (A2) yields the desired result (C11) where f is given in Eq. (8). Note that we eliminated to components from f because the system would otherwise be overdetermined owing to the two conserved masses. This gives the derivative matrix Note that the first to rows are simply s T D and s T E , which follows immediately from the definition of F. c. Inhomogeneous (asymmetric) steady states. The derivation presented above for homogeneous steady sates can be generalized to inhomogeneous steady states ∆ñ defined by ∆u * (∆ñ) = 0. The resulting Jacobian reads J αβ = s T α D · (D u F|n +∆ñ ) −1 + (D u F|n −∆ñ ) −1 e β . (m a , m p , c A , c P ) with the reaction kinetics These reactions conserve the total densities of aPARs n A = m a + c A and pPARs n P = m p + c P , respectively. Since the reaction network is symmetric under the exchange A↔P, we use reaction rates that also respect this symmetry for simplicity [23]. The diffusion matrix reads D = 4/L 2 diag(D m , D m , D c , D c ), where L ≈ 15 µm is the long half-axis of the ellipsoidal cells. Note that in LQSSA, this length only contributes to the overall time scale but does not affect the phase portrait structure.

Cdc42 polarity model
We use a simplified form of the quantitative model proposed in [71]. This model describes the dynamics of the GTPase Cdc42, its guanine nucleotide exchange factor (GEF) Cdc24 and the scaffold protein Bem1. The critical feedback loop is constituted by mutual recruitment between membrane-bound Cdc42-GTP and Bem1-GEF complexes. While the full model accounts for Bem1 and GEF separately, we lump these species into a complex species here. This retains the salient features of the model, in particular the mutual recruitment mechanism.
The variables of this simplified model are u = (m t , m d , m b , c D , c B ), accounting, respectively, for membrane-bound Cdc42-GTP, Cdc42-GTP and Bem1-GEF complexes as well as cytosolic Cdc42-GDP and Bem1-GEF complexes.
The reaction kinetics, describing attachment and detachment of Cdc42 at the membrane, hydrolysis and nucleotide exchange of Cdc42 of membrane-bound Cdc42, as well recruitment of Bem1-GEF complexes to the membrane by Cdc42-GTP are given by These reactions conserve the total densities of Cdc42 n D = m t + m d + c D and Bem1-GEF complexes n B = m b + c B , respectively.
The parameter values, given in Table III are adapted from [71]. The values of k b , k B and k tB are chosen to account for the lumped Bem1-GEF complexes species.