Bridging scales in a multiscale pattern-forming system

Significance Biological processes operate in a spatially and temporally ordered manner to reliably fulfill their function. This is achieved by pattern formation, which generally involves many different spatial and temporal scales. The resulting multiscale patterns exhibit complex dynamics for which it is difficult to find a simplified description at large scales while preserving information about the patterns at small scales. Here, we introduce an approach for mass-conserving reaction–diffusion systems that is based on a linear theory and therefore conceptually simple to apply. We investigate multiscale patterns of the Min protein system and show that our approach enables us to explain and predict the intricate dynamics from the large-scale mass redistribution of the total protein densities.


Introduction
Pattern formation is fundamental for the spatiotemporal organization of biological processes, such as cell division, chemotaxis, and morphogenesis.More than half a century ago, Turing showed theoretically how local interactions (chemical reactions) and diffusion of chemical species can lead to spontaneous spatial patterns [1].Such reaction-diffusion systems have been successfully used to explain pattern formation phenomena in nature that arise self-organized from a stable homogeneous steady state [2,3,4,5].The analysis proposed by Turing allows to predict the emergence of patterns with a characteristic length scale as long as the entire dynamics remains in the vicinity of the homogeneous steady state [6].The validity of Turing's approach has been also tested experimentally for coupled chemical oscillators, and was found to reliably predict the experimental observations, provided that the model parameters are spatially and temporally uniform [7].Pattern-forming systems, however, are generally heterogeneous and therefore far from homogeneity, and involve multiple spatial and temporal scales.An intriguing example of biological pattern formation is morphogenesis, in which the spatiotemporal patterns of morphogens dictate the future shape of an organism that is orders of magnitude larger than its constituents [4].On a smaller scale, protein concentration patterns in cells are essential for the spatiotemporal control of cellular processes such as cell division and motility [8,5,9].Protein patterns can exhibit fascinating multiscale characteristics [10] and form in hierarchies of patterns on several scales that affect one another [11].
Such complex multiscale biological processes involve many degrees of freedom at multiple scales, rendering it difficult to analyze them and gain insight into the underlying principles.To make progress on this issue, one needs to use systematic coarse-graining schemes that allow the dynamics to be reduced to the essential degrees of freedom at the relevant time and length scales.For instance, a well-known and powerful method is the renormalization group theory [12].Unfortunately, this method is restricted to the vicinity of critical points.The Mori-Zwanzig formalism [13] is another important approach which allows to decompose the dynamics of a system into 'fast' and 'slow' variables by means of projection operators.One arrives at a closed set of equations for the slow variables, while the fast variables are treated as noise.One property that these methods have in common is that the scales that have been integrated out or eliminated are not resolved, and cannot be recovered from the coarsegrained level of description.This is most apparent in the Mori-Zwanzig formalism, where the eliminated degrees of freedom appear effectively as noise terms on the resolved scales.For pattern-forming systems, one is however interested in the patterns on the unresolved scales 1 as they usually have a specific function in biological systems.This raises the question of whether it is possible to reconstruct information about the unresolved scales from the dynamics at the resolved scales?Indeed, amplitude equations describe the long-wavelength amplitude modulations of an underlying short-wavelength base pattern and therefore resolve both the small and the large scales.Unfortunately, however, they are limited to the vicinity of the supercritical onset of pattern formation [6] (including weakly subcritical cases) and only feasible in simple geometries where the orthonormal basis functions of the diffusion operator can be found in closed analytical form.Hence, to fill these gaps, one relies on new concepts to deal with multiscale systems.
Here, we propose a semi-phenomenological approach to overcome these mathematical limitations in the concrete context of mass-conserving reaction-diffusion (MCRD) systems.Recently, a new theoretical framework for MCRD systems has been introduced [14,15] that allows one to characterize their dynamics in the highly nonlinear regime.The basic idea is to consider reaction-diffusion system as decomposed into a set of reactive compartments which are spatially coupled by diffusion.For an isolated compartment, one can determine the steady state (local equilibrium) and its stability properties which both depend on the total densities within that compartment.Since diffusion causes the lateral redistribution of these total densities, these local equilibria will change over time.This concept of moving local equilibria enables one to study the physical mechanisms underlying pattern formation and characterize the dynamics far away from the homogeneous steady state.The fact that one is able to characterize the dynamics far from homogeneity suggests that the local equilibria theory may be a promising approach to study heterogeneous systems.We therefore asked whether the ideas from local equilibria theory would be applicable to investigate multiscale patterns?
To pursue this question, we use the Min protein system of E. coli which has emerged as a paradigmatic model system for the study of pattern formation in cell biology [16,17,18,19,20].Its dynamics is driven by two proteins, MinD and MinE, which cycle between cytosolic and membrane-bound states and interact nonlinearly on the membrane (Fig. 1A).In E. coli, these proteins oscillate from cell pole to cell pole and thereby position the cell division machinery to midcell [16,17].Studying the Min dynamics in various reconstituted systems has led to the discovery of a rich set of patterns including traveling waves and spirals [18], chaotic patterns [21,22,23,10], "homogeneous pulsing" [24,25,26], as well as quasi-stationary labyrinths, spots, and mesh-like patterns [10,27].Theoretical analysis of mathematical models has lead to the key insight -and experimentally confirmed prediction -that the average total densities of MinD and MinE and the bulk height are key control parameters for pattern formation in the reconstituted Min system [5,28].The rich set of patterns, experimental accessibility in vitro and theoretical understanding make the Min-system an ideal candidate to investigate the role of spatial heterogeneity on pattern formation.
Since varying the bulk height affects the local equilibrium state and is a key control parameter for pattern formation [5,28], we study the Min dynamics in a wedge-shaped geometry with a membrane placed on the bottom surface (Fig. 1B).While there are many distinct ways to introduce large-scale spatial heterogeneities into the system, e.g. by introducing space-dependent kinetic rates, we chose to use a wedge geometry because it is relatively easy to implement experimentally.In numerical simulations, we find that the system exhibits a striking range of transient patterns, that coexist in different spatial regions along the membrane (Movie S1 and Fig. 1C).As time progresses, patterns in different regions change and transition to other patterns.
To characterize these complex dynamics that play out on multiple spatial and temporal scales, we generalize the concept of dispersion relations (obtained from a linear stability analysis) by applying it to sections of the domain, which we term regional dispersion relations.Combining this approach with the local equilibria theory [14,15,8], we show that one can reconstruct the type and characteristics of patterns on small scales from the local protein mass densities, which we identify as the essential degrees of freedom on large spatial and temporal scales, i.e. the "hydrodynamic variables" of the system.The key to this reconstruction are correlations between the regional pattern characteristics and instantaneous, regional dispersion relations, calculated from the instantaneous regional mass densities.Over time, these masses change due to diffusive redistribution, resulting in qual-  itatively different regional dispersion relations that indicate the local pattern type in the system.This reconstruction of small-scale features (on unresolved scales), together with a coarse-grained description for the mass-redistribution dynamics on large scales allows us to understand and predict the long-term temporal evolution of the system.A major advantage of our approach is that it is based on a linear theory and therefore conceptually and technically simple to apply.
A key prediction from our numerical simulations and theoretical analysis is that different pattern types form at different positions along the wedge shaped geometry.To test this prediction experimentally, we performed experiments with a reconstituted Min system in wedge-shaped microfluidic cells.In agreement with the theoretical prediction, we find a range of transient patterns coexisting in different spatial regions along the membrane.

Results
The Min protein system in wedge geometry Mathematically, the Min-protein dynamics is described by bulk-surface coupled reactiondiffusion equations, which describe the concentrations of cytosolic proteins MinD-ATP, MinD-ADP, and MinE, c = (c DD , c DT , c E ), in the bulk volume V, and the concentrations of membrane-bound MinD and MinDE complexes, m = (m d , m de ), on the surface S. For the wedge geometry, in spatial coordinates x = (x, y, z), we place the membrane surface (with lateral dimensions L × L) in the x−y plane at z = 0 and let the bulk height vary as a linear ramp from H 0 to H 1 along the x-direction (see Fig. 1B).
The dynamics of bulk components c(x, t) is governed by the equation where D c denotes the bulk diffusion constant and the matrix Λ = diag(−λ, λ, 0) describes nucleotide exchange of MinD in the bulk.The dynamics of membrane components m(x, y, t) is constrained to the membrane surface and takes the form: where D m is the membrane diffusion constant and ∇ 2 S = ∂ 2 x + ∂ 2 y is the surface Laplacian.The membrane reactions r, which comprise attachment, detachment, and recruitment processes of Min proteins, are specified in the Materials and Methods section.
The dynamics in the bulk and on the surface are coupled by reactive boundary conditions, that describe the bulk fluxes induced by attachment and detachment of proteins at the membrane (see Materials and Methods).At the remaining boundaries, no-flux boundary conditions are imposed such that the system is closed.Together, the above dynamics conserve the average mass densities of MinD and MinE: where c D = c DD + c DT is the total cytosolic MinD concentration; • S and • V denote the mean on the surface and in the bulk respectively; |S| and |V| are the total surface area and bulk volume (see Materials and Methods).
Using finite element (FEM) simulations we investigated the spatiotemporal dynamics of the Min system in wedge geometry.Our simulations show a broad range of different patterns -including traveling waves, standing waves and chaotic patterns -coexisting in different spatial regions of the membrane (see Movie S1 and Fig. 1C).Interestingly, the regions where these patterns are found change over time as the patterns transition from one type to another.For long simulation times, we observe that patterns transition to standing waves, such that the entire domain is covered by a single pattern type in the final steady state.The pattern in steady state depends on the specific choice of parameters, and therefore can be altered by changing the model parameters (Fig. S1 and Movie S2).

Experimental implementation
We tested our theoretical prediction on this multi-scale dynamics in an experimental system consisting of a wedge-shaped microfluidic flow chamber (Fig. 2A).The bottom and top surface of the wedge were covered with a supported lipid bilayer consisting of DOPG:DOPC (30:70 %) which mimics the natural membrane composition of E. coli [29].The length of the wedge was typically about 8 − 14 mm and the width about 3 − 4 mm.The bulk height range was approximately 2 − 50 µm (Fig. 2B).Min proteins were distributed in the chamber by rapid injection of a solution containing 1 µM MinD and 1 µM MinE (including 10 % fluorescently labelled MinD and MinE proteins for visualization), together with 5 mM ATP and an ATP-regeneration system [28].
Figure 2C shows a snapshot of Min protein patterns along the bottom surface of the wedge geometry 30 minutes after injection.The experiments exhibit the same essential hallmarks of multiscale Min protein patterns that we observed in our numerical simulations.In particular, consistent with our simulations, we observe a sequence of distinct spatiotemporal patterns coexisting in different spatial regions of the membrane (Fig. 2C and Movie S3): At regions of low bulk height (approximately between 2 − 10 µm), one typically observes chaotic patterns and standing waves, whereas traveling wave patterns emerge at regions of large bulk height (> 10 µm).Furthermore, as in the simulation, we observe a sharp boundary between regions that contain traveling wave patterns and regions that contain rather chaotic and standing wave patterns, and this boundary establishes quickly within a few minutes (Fig. S2 and Movie S4).Overall, the observations provide a striking verification of the height-dependent patterns predicted in the simulations.
There are also some differences between the patterns in the experiment and in our numerical simulations.First, while we observed occasional transitions from one pattern into another in our experiments (Fig. S3 and Movie S5), these transitions occurred frequently and were more pronounced in the simulations.This is explained by the lateral length of the experimental setup, that is about an order of magnitude larger as compared to the simulation setup, which is the main reason why we observe more frequent transitions between different patterns in the simulations, as will become clear later.Second, in contrast to the simulations, we noticed some homogeneous oscillations in the experiments, which are characterized by large (homogeneous) density patches on the membrane (typically few hundred micrometers in size) that oscillate with time (Figs.S3-S4 and Movies S5-S7).We attribute this difference to the following: Due to the fabrication method of the microfluidic flow chamber, both the bottom and top surface of the wedge were covered with a supported lipid bilayer.In recent work, it has been shown that membrane-to-membrane crosstalk (i.e., between top and bottom surface) is responsible for the emergence of homogeneous oscillations [28].In our simulations, however, we assume that Min proteins can only bind to the bottom membrane, which explains why we do not observe homogeneous oscillations.
Taken together, we have a system that exhibits a fascinatingly rich transient dynamics and involves patterns and transitions between them on multiple spatial and temporal scales.We are therefore left with the key question: Can we explain the cause why different patterns form in different spatial regions and how they transition from one to another over time?Moreover, is it possible to identify and reduce the system to its essential degrees of freedom?A standard way to address these questions mathematically would be to perform a multiscale analysis and to derive amplitude equations that describe the large-scale spatiotemporal evolution of the pattern amplitudes [6].This would greatly simplify the problem as it allows to obtain a quantitative relationship between the small-scale patterns and the large-scale dynamics (slowly varying pattern amplitudes), thus ultimately enabling one to reconstruct the patterns from the reduced dynamics at large length and time scales [30,31,32,33].Carrying out this analysis requires determining the set of orthogonal eigenmodes for the diffusion operator that satisfy the boundary conditions.In a one-dimensional domain, these eigenmodes are simply Fourier modes.Unfortunately, in the wedge geometry with bulksurface coupling, the eigenmodes can not be found analytically, thus precluding the use of the amplitude equation framework.Moreover, amplitude equations are restricted to the vicinity of supercritical and weakly subcritical bifurcations [6,34].The Min patterns we observe here, however, are generically subcritical [15] and exhibit large amplitudes [14,28].We therefore aim to develop a new approach that overcomes these restrictions.

Instantaneous, regional dispersion relations predict patterns
The analysis of pattern-forming systems usually starts with calculating the homogeneous steady state (HSS) solutions and performing a linear stability analysis around these states.This yields a dispersion relation that informs about the growth rate σ(q) of small spatial perturbations with a certain wavenumber q.However, the dispersion relation is generally only informative in the vicinity of the homogeneous steady state [1,6], and thus unreliable for large amplitude patterns.Moreover, the spatial variation of parameters even precludes the existence of a global HSS, so that a global dispersion relation can no longer be determined.To overcome these limitations, we adopt a semi-phenomenological approach where we generalize the concept of dispersion relations.
Let us consider the wedge as dissected into a collection of two-dimensional slices along the direction of constant bulk height.Each slice corresponds to a rectangular geometry with a bulk height that depends on the position of the slice in the wedge (see Fig. 3A,B).Next, for each slice and at each point in time, we calculate instantaneous total densities of MinD and MinE, averaged over the slice length ñD,E y (t, x) (Materials and Methods).The average total densities, together with the local bulk height H(x), then serve as parameters for the regional dispersion relation in each slice which is straightforward to determine because the slice represents a rectangular geometry [14,23,28] (see Fig. 3A,B and Supplementary Information).While the bulk height H(x) varies linearly in space, the average total densities ñD,E y (t, x) are dynamic quantities and depend on the slice position x as well as on time t, since the diffusive coupling between the slices redistributes mass.It follows that the regional dispersion relation depends on the spatial position and is dynamic: σ(q; x, t).This generalizes classical dispersion relations, which are by definition independent of space and time.How does this spatially and temporally varying dispersion relation inform about the system's dynamics?As in uniform systems that exhibit homogeneous steady states, it serves as a criterion for the onset of pattern formation and for estimating the characteristic wavelength of the initial pattern that is formed.While these insights are generally limited to the linear regime [1,6], recent theoretical findings for the Min system in a two-dimensional rectangular geometry (representing a slice geometry) have shown that the dispersion relation reliably predicts the pattern type in the fully nonlinear regime [5].In particular, it was shown that depending on the total densities of Min proteins, nD and nE , and the bulk height H, the system exhibits a variety of different patterns on the membrane, such as chaos, standing waves, and traveling waves [14,28].Moreover, a careful analysis of nu- For each such slice, at a given instance in time, we calculate the instantaneous total densities, averaged along its length ñD,E y (t, x), from the numerical simulation data.From these slice-averaged total densities, we can then calculate the corresponding local homogeneous steady state and its dispersion relation.(C ) Dispersion relation with fastest growing mode q * and right edge of the band of unstable modes q max indicated.The ratio q max /q * has been empirically found to correlate with the type of the fully developed pattern, with a sharp transition from chaotic patterns for q max /q * < 2 to ordered patterns for q max /q * > 2. Close to the transition, standing waves are found, while travelling waves form for larger ratios q max /q * [14].(D) Mode ratio q max /q * as a function of the slice position x for a given instance in time.The background shading indicates the type of pattern expected from the "commensurability criterion."(E ) Representative snapshots of the three distinct pattern types: spatiotemporal chaos, standing waves (SW) and traveling waves (TW).merical simulations has interestingly revealed a strong one-to-one correlation between the dispersion relation and the fully developed patterns in the highly nonlinear regime [14]: A commensurability criterion between the unstable mode with the shortest wavelength q max and the fastest growing mode q * has been found that determines the pattern type (Fig. 3C-E).In short, it has been shown that q max /q * < 2 coincides with the regime of chemical turbulence (spatiotemporal chaos), whereas for q max /q * > 2 the system exhibits ordered patterns (standing/traveling waves).Standing wave patterns are found close to the commensurability transition q max /q * 2, while traveling waves are found further away from the threshold.In the following, we use this observed one-to-one correspondence between the dispersion relation and the fully developed patterns to reconstruct the small scale pattern types from coarse grained densities.
To that end, we extracted the average total densities in each slice as a function of time from the numerical simulation.Based on these densities we then calculated the instantaneous regional dispersion relation in each slice and extracted the ratio q max /q * as a function of slice position x and time t (Fig. 3C-E).The resulting pattern-type prediction is shown in the space-time plot (kymograph) in Fig. 4A. Figure 4B shows the ratio q max /q * as a function of slice position x for a set of representative times (cf.Fig. 3D).The pattern-type prediction Fig. 4A is then obtained from these ratios via the mapping shown in Fig. 3D,E.
We find that this prediction correlates well with the patterns observed in the full numerical simulation (Fig. 4C,D and Movie S8).In particular, the temporally changing position x crit (t), marking regions where q max /q * = 2 (indicated by the green arrows and dashed lines in Fig. 4B and C), agrees with the position along the wedge where traveling wave patterns transition to chaotic patterns.In the vicinity of x crit (t) we observe a band of standing waves as expected from the "commensurability criterion" [14].While the transition from chaos to order at q max /q * = 2 is sharp, we were not able to identify a sharp criterion for the transition from standing waves to traveling waves.Since the ratio q max /q * and with it x crit (t) are entirely determined by the slice-averaged masses ñD,E y (x, t), we conclude that these masses are the essential degrees of freedom of the system at large scales.
Notably, we find that there are slight differences between the predictions and the actual patterns for large times (see Fig. 4A-C).The reason for these deviations lies in the model parameters, which were chosen such that the entire domain is near the critical mode ratio q max /q * = 2 for large times.This renders the dynamics, and the prediction from the regional dispersion relation highly sensitive to slight variations of the regional total masses.Hence, the fact that our method is still able to qualitatively predict the dynamics in this case underscores the robustness of our approach.In the Supplemental Information, we provide additional results where the parameters were chosen such that the mode ratio is deep in the traveling wave regime (q max /q * > 2) for late times.In this case, we obtain an excellent agreement between our predictions and the patterns observed in the numerical simulations (see Fig. S1).
Next, we ask whether one can find an approximate coarse-grained dynamics for these redistributed masses.Such a description would enable us to predict the time evolution of the redistributed masses independently from the full numerical simulations.One can then use the commensurability criterion to predict the pattern types that will form in different spatial regions as a function of the redistributed masses.In the next section we will show how one can find such a description.

Large-scale dynamics is driven by redistribution of mass
In general, mass redistribution between different spatial regions of the wedge is caused by diffusive fluxes due to concentration gradients.Similar as in the previous section, we consider here the redistribution of mass between slices along the wedge (Fig. 3B).Since membrane diffusion is by two orders of magnitude slower than bulk diffusion it may be neglected, such that redistribution of protein mass between slices is governed by bulk diffusion  3D).The green line shows x crit (t) where q max /q * = 2, indicating the transition from chaotic to ordered patterns.Green arrows mark the position x crit (t) for the times indicated by dashed white lines.(B) Plots of the mode ratio q max /q * , determined from the local dispersion relation, as a function of spatial position x for several representative times (dashed white lines in (A)).In the second to last row, the entire domain is near the critical ratio q max /q * = 2, predicting the global emergence of standing waves (see last row).(C ) Snapshots of the membrane patterns (MinD density, cf.Fig. 1) from the full numerical simulation.The green dashed line indicates x crit (t).Note the standing wave patterns found near x crit (t).Their fronts are aligned along the bulk height gradient such that the sequence of wavenodes lies on lines of constant bulk height.(D) Machine-learning based pattern classification using ilastik [35] (see Materials and Methods).

alone (Materials and Methods)
for i = D, E. Here, the second term accounts for the spatial variation of the bulk height, and thus the different volumes of neighboring slices between which the diffusive flux D c ∂ x c i y,z redistributes mass.This can be seen by rewriting Eq. ( 6) in the form of a continuity equation with the diffusive fluxes given by J diff i := −D c ∂ x c i y,z .Since the area of slices increases along the positive x − direction, the diffusive fluxes J diff i on the right-hand side of Eq. ( 7) are rescaled by the bulk height H(x).These equations seem to be simple, but unfortunately they are not closed, since the slice-averaged cytosolic densities c i y,z (x, t) appear on the right hand side.
We are interested in the dynamics of n i y,z on timescales much longer than typical oscillation periods of the patterns.Therefore, following the intuition gained from previous works on MCRD systems [15,36], we assume that one can approximate the slice-averaged cytosol concentrations by the homogeneous steady-state concentration in each slice This assumes that the spatial average over many wavelengths in y-direction is well approximated by the instantaneous homogeneous steady state in a slice.These steady state concentrations only depend on the slices bulk height H(x) and the slice-averaged total densities n i y,z (x, t).Thus, the above approximation yields a closed set of equations for the mass-densities We will call this the reduced dynamics in the following.Since the homogeneous steady states may also undergo a saddle-node bifurcation, characterized by the emergence of three steady states (two stable, one unstable), this may lead to discontinuities in c * i .To regularize the dynamics, c i is not set identical to c * i but relaxes towards it on a fast timescale (see SI for details).
Given the initial densities n i y,z (x, 0), one can numerically solve the reduced dynamics Eq. ( 9) to predict the entire time evolution of the slice-averaged masses and hence the dispersion relation at each point along the x − direction.Figure 5C shows the regional pattern types predicted from the reduced dynamics.We find good qualitative agreement for the distribution and transition of patterns as observed in the numerical simulations (cf.Fig. 4A).The main difference to the full numerical simulations is a slight quantitative deviation in the timescale, where the dynamics predicted by Eq. ( 9) is slightly slower compared to the full numerical simulation.We also note that the reduced dynamics predicts a larger region of no instabilities as compared to the numerical simulations (cf.Figs.4A and 5C).This is because the chaotic regime is rather narrow and close to the regime for which the dispersion relation predicts no instability (cf.Figs.3D and 4B).In addition, since the patterns emerge from a subcritical bifurcation [14] (a generic property of mass-conserving systems [15]), large amplitude patterns can be excited and maintained even below the instability threshold.
Figure 5A,B compare the time evolution of the slice-averaged total densities from the full numerical simulation and the solution obtained from the reduced dynamics.The colors in the kymographs indicate the total density ratio of MinE and MinD (short, E:D ratio), which is a key control parameter in the Min-protein dynamics [14].

Discussion
Multiscale patterns in biological systems often emerge from hierarchical systems, which are organized in a modular fashion.Each level of the hierarchy instructs dynamics on the next level which operates on a smaller spatial scale.For instance, along developmental trajectories of many organisms, upstream patterns such as maternal gradients instruct downstream gene-expression patterns on increasingly smaller scales [37,11].Importantly, on each level of the hierarchy, there is a clean separation between (spatially varying) control parameters and dynamical variables.
In contrast, in the system we have studied here, there is no such separation as the globally conserved total densities play a dual role: they are both dynamical variables and act as control parameters [14,15].Building on this key feature has allowed us to explain and predict the intriguingly complex patterns found in large-scale numerical simulations.The values of the total densities of MinD and MinE locally control the pattern type: we showed that a "regional dispersion relation" calculated from the regional average densities reliably predicts the pattern type.At the same time, concentration gradients in the bulk drive mass redistribution of MinD and MinE.Therefore, the total densities are hydrodynamic variables on large scales which control pattern formation on small scales.This separation of scales enabled us to derive a reduced dynamics for the total densities on large spatial and temporal scales which predicts the long-term dynamics of the system.
Notably, the dual role of total densities as dynamic variables and control parameters also plays out at the small scale of the patterns themselves [14,15].Here, instantaneous local total densities control local equilibria and their stability, which serve as proxies for the local dynamics.The local dynamics cause gradients, which drive diffusive redistribution of the total densities-in turn causing changes in the local dynamics.In the Min system, this point of view has led to a detailed understanding of the emergence of chaos near onset and of the transition to standing and traveling waves [14].From a general perspective, the concept of local equilibria controlled by total local densities is at the core of a number of recent theoretical advances in the field of mass-conserving, pattern-forming systems [15,36,38,8].
In addition to the dynamically changing total densities, the bulk height is also a (fixed) heterogeneous control parameter in our system.The bulk height (or more generally volumeto-surface ratio) is an important control parameter for bulk-surface coupled pattern-forming systems [14,28].Here, the bulk height gradient of the wedge serves to induce spatiotemporal heterogeneities in the total densities.Alternatively, one could induce heterogeneities in the total densities via spatial gradients of the kinetic rates or by imposing a heterogeneous initial condition in the total densities.However, these alternatives are difficult to realize experimentally in a reproducible and controlled manner, which is the main reason why we chose the wedge setup in this work.In a third scenario, large-scale gradients in the densities may also emerge spontaneously and be maintained in the absence of "external" heterogeneities.
An example for this third scenario is the Aranson-Tsimring model for pattern formation in vibrated granular media [39] (see Materials and Methods for details).In the following, we briefly discuss this model to put our approach into a broader context.In particular, this model has been extensively studied using amplitude equations allowing us to connect this mathematical approach to the regional dispersion relations introduced here.The Aranson-Tsimring model considers a system with a complex order parameter ψ (describing the surface modulation of a vibrated granular layer) which is coupled to a conservation law for the grain density ρ (see Eq. ( 24) in Materials and Methods).Near the onset of pattern formation, this coupling gives rise to localized patterns that have been studied using amplitude equations [32,33,30].Figure 6 and Movie S9 illustrate how these patterns can be understood in terms of regional dispersion relations.For high densities, there are no unstable modes and no patterns form.Below a critical density ρ c , a band of unstable modes appears, giving rise to patterns through a supercritical bifurcation.Indeed, localized patterns appear only where the average regional density is below ρ c (see Fig. 6B).This demonstrates the idea of regional dispersion relations in a nutshell.Moreover, it shows that this approach gives rise to qualitatively similar insights as the technically much more involved amplitude equation formalism.The conceptual and technical simplicity of regional dispersion relations make this approach readily applicable.The caveat is that this approach lacks the mathematical rigor of the amplitude equation formalism and requires numerical solutions of the dynamics as a basis.Because the bifurcation in the Aranson-Tsimring model is supercritical [39], we can immediately read off from the regional dispersion relation where small-scale patterns will form.There is only one pattern type (stripes) and therefore, no additional information is needed to reconstruct the small-scale patterns.For the Min system considered here, on the other hand, the onset is subcritical (i.e.patterns have large amplitude at onset), and the emergence of a band of unstable modes alone does not inform about the pattern type.To overcome this problem, we used the "commensurability criterion" that enabled us to predict the small-scale patterns from regional dispersion relations.However, this criterion has so far only been shown to hold for the Min-protein system.Whether it also applies to other reaction-diffusion systems remains an open question.In general, reconstructing subcritical small-scale patterns from large-scale quantities will require that adequate criteria are first identified in simplified settings (such as the "slice geometry" used here).Supercriticality also guarantees that there is no multistability of different pattern types near onset.Multistability would lead to hysteresis in the transitions and therefore introduce memory in the system.As a result, the one-to-one correspondence between the regional dispersion relation and the pattern type that we have used here to reconstruct patterns would be lost.Handling memory effects in pattern-forming systems remains an open issue that will require the development of new methods, providing an interesting task for future research.
Since conservation laws are ubiquitous in many physical systems, we believe that our approach can be generalized to a broad class of multiscale pattern-forming systems.For instance, mass conservation is inherent to particle-based active matter systems.The local particle density controls emergent orientational order, i.e. local symmetry breaking [40,41,42].In turn, orientational order controls mass redistribution due to the particles' self-propulsion.Thus, the particle density again plays a dual role as a control parameter and a dynamic variable [42,43,44].The dynamic interplay of mass redistribution and orientational order has been shown to give rise to coexistence of different macroscopic order (polar flocks, nematic lanes) and the interconversion between them [42], not unlike the coexistence and interconversion of different patterns we found for the reaction-diffusion system studied in this work.One way to induce spatial heterogeneities in these systems is to introduce a gradient of signaling chemicals (chemoattractants) that affect the local velocity of active particles.This would dynamically lead to redistribution of the particle densities on large scales.Since the particle densities, in turn, are themselves control parameters locally, non-trivial multiscale dynamics may emerge in such a setup.Exploring the effects of such gradients in active matter systems could be therefore an exciting task for future research.
On a broader perspective, our work shows how a linear analysis on small scales, combined with a reduced description for non-linear large scale dynamics (mass redistribution) can be employed to study complex multiscale phenomena.We believe that our approach can be generalized and applied to other multiscale systems with an underlying conservation law, such as transport processes in porous media, combustion, and cell migration, to name a few examples.

Mathematical model
We adopt the Min "skeleton model" introduced in Refs.[45,46,5] which is known to qualitatively reproduce Min patterns in vivo and in vitro [46,5,28].The governing equations are given in the main text, Eqs.[1]- [3].The membrane reactions are with The reaction terms account for MinD attachment and self-recruitment to the membrane, MinE recruitment by MinD, and dissociation of MinDE complexes with subsequent detachment of both proteins into the cytosol, respectively.Coupling between cytosol and membrane is established by reactive boundary conditions at the membrane [cf.Eq. ( 3)].
The boundary fluxes are given by The model parameters used in this study are summarized in Table 1.

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

Homogeneous steady state and dispersion relation
The homogeneous steady state concentrations, (c * | z=0 (H, nD , nE ), m * (H, nD , nE )) are obtained from the stationary solutions of Eqs.[1]- [3] together with the mass conservation condition Eq. ( 4): where Φ denotes the steady state fluxes at the membrane, given by: A concise derivation of these equations and how they can be solved is provided in the Supplementary Information.For a thorough presentation of the linear stability analysis of the Min system in a 2D rectangular geometry we refer to the Supplementary Informations of Refs.[14] and [28].

Operators for spatial averaging
The operators used throughout this study to calculate mean values of densities on the membrane and in the cytosol are defined as follows: • y,z := 1 H(x) where the membrane surface area and the bulk volume for the wedge geometry are explicitly given by |S| = L 2 and |V| = L 2 (H 0 + H 1 )/2.

Instantaneous total densities at the membrane
Since only cytosolic proteins in close proximity to the membrane participate in the nonlinear dynamics at the membrane, we define instantaneous total densities at the membrane: We further averaged these densities along the y-direction to obtain the the slice-averaged total densities ñD,E y (x, t).Note that the length of a slice is much larger than the typical pattern wavelength, which also permits to approximate the slice-averaged mass at the membrane by the vertically averaged mass: ñi y (x, t) ≈ n i y,z (x, t) (see Ref. [14]).This is because the local deviations ñi − n i z largely cancel when averaging over the pattern wavelength.

Mass redistribution dynamics
Here, we provide more details on the derivation of the mass redistribution dynamics Eq. ( 7).
For specificity, we present the calculation for MinD.The calculation for MinE works along the same lines.Our starting point is the slice averaged total MinD density: The time evolution of this quantity then follows from Eq. ( 1) and Eq. ( 2): where we used the reactive boundary condition Eq. ( 3) to rewrite the integral: Note that due to mass-conservation the reaction terms at the membrane cancel.
Since the system is closed, the boundary condition at the inclined top surface of the wedge reads n • ∇c D | z=H(x) = 0, where n ∝ (−∂ x H, 0, 1) is the outward normal vector at the top surface.Writing out the boundary condition explicitly, we find that: To proceed, we substitute the relation above into Eq.( 19) and slightly rewrite the resulting equation by applying the chain rule: Here, the first term describes diffusion of the averaged membrane concentrations.The integral on the right describes diffusion of the averaged cytosolic densities, where we defined the diffusive flux J D = −D c ∂ x c D y,z .The factor H(x) in the cytosolic diffusion term accounts for the increasing area of the slice along the positive x-direction.Since protein diffusion on the membrane is much smaller than cytosolic diffusion D m D c [47,48], one can neglect membrane diffusion to arrive at the result shown in the main text (Eq.( 7)).For completeness, note that Eq. ( 22) (without membrane diffusion) can be recast as which is the form given in Eq. ( 6) in the main text.

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

Aranson-Tsimring model
As a second example we briefly discuss a phenomenological model for pattern formation in vibrated granular media introduced in [39].This model, which we call Aranson-Tsimring model in the following, couples a Ginzburg-Landau-type equation [34] for the complex order parameter ψ to a conservation law for the density ρ: where ψ denotes the complex conjugate of ψ.The coupling is such that increasing the density ρ suppresses the instability in Eq. (24a) while gradients in the amplitude |ψ| drive mass redistribution away from high amplitude regions (second term in Eq. ( 24b)).This feedback loop amplifies heterogeneities in the density and gives rise to localized patterns.These patterns have been studied in detail using amplitude equations in [32,33].Moreover, in Ref. [30] it was shown that the system Eq.( 24) appears as the amplitude equation for a mass-conserving version of the classical Swift-Hohenberg-Turing equation [6,49].The reason for this is that the conserved density appears as a second hydrodynamic variable in addition to the pattern amplitude.

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

Observation of Min patterns
We purified the Min proteins based on the method proposed in Ref. [50].Injection of Min proteins into the flow cells was performed through a syringe pump containing a solution of 0.8 M MinD, 0.2 mM MinD-Cy3, 0.8 mM MinE, 0.2 mM MinE-Cy5, 5 mM ATP, 4 mM phosphoenolpyruvate, 0.01 mg/ml pyruvate kinase, 25 mM Tris-HCl (pH 7.5), 150 mM KCl and 5 mM MgCl2.To ensure that all of the buffer solution in the microdevice is replaced by the protein solution, we chose a volume of the protein solution that was 50 times larger than the volume in the microdevice.During the filling process of the microdevice, the enire solution was rapidly injected (in 5 s) to prevent protein accumulation on the membrane.
For the generation of the fluorescence images, we used the following equipment: Olympus IX-81 inverted microscope equipped with an Andor Revolution XD spinning disk system with FRAPPA, illumination and detection system Andor Revolution and Yokogawa CSU X1, EM-CCD Andor iXon X3 DU897 camera, motorized x-y stage and a z-piezo stage, using a 20x objective (UPlansApo, NA 0.85, oil immersion).Imaging of MinD-Cy3 and MinE-Cy5 was performed with laser spectral lines at 561 nm and 640 nm, respectively, and we further used a 617/73 band-pass filter as well as a 690 long-pass filter.We imaged several uniformly sized regions at intervals of 30 s or 60 s along the lateral length of the wedge setup.To exclude membrane imperfections that may have arisen during preparation, we also imaged the membrane using the spectral line at 491 nm and a 525/50 band-pass filter.

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

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

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

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

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

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

Fig. 6 .
Fig. 6.Regional dispersion relations predict localized patterns in the Aranson-Tsimring model Eq.(24).(A) Snapshot of the order parameter magnitude ψ showing localized patterns.Dashed white line indicates the stability threshold determined from regional dispersion relations.(B) Coarse-grained density (Gaussian filter with standard deviation 10).(C ) Representative dispersion relations in the stable and unstable regimes.Domain size: 100 × 50; see Materials and Methods for model details and remaining parmeters.
which follows from mass conservation.For analytical caclulations we adapt the following change of variables as it is more convenient: We describe the bulk dynamics of MinD in terms of the variables c D = c DD + c DT and c DD , i.e. in this case one defines the bulk concentration vector c = (c D , c DD , c E ).The membrane reaction in Eq. (11a) is then slightly modified by substituting c DT = c D − c DD , and the boundary fluxes are given by f = −r on D , r off DE , r off DE − r on E .

Table 1 .
Min model parameters