Bridging scales in a multiscale pattern-forming system.

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

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 (1) showed, theoretically, how local interactions (chemical reactions) and diffusion of chemical species can lead to spontaneous spatial patterns. 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 (HSS) (2–5). The analysis proposed by Turing allows prediction of the emergence of patterns with a characteristic length scale, as long as the entire dynamics remains in the vicinity of the HSS (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 (5, 8, 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 decomposition of the dynamics of a system into “fast” and “slow” variables by means of projection operators. One arrives at a closed set of equations for the slow variables, while the fast variables are treated as noise. One property that these methods have in common is that the scales that have been integrated out or eliminated are not resolved, and cannot be recovered from the coarse-grained level of description. This is most apparent in the Mori–Zwanzig formalism, where the eliminated degrees of freedom appear effectively as noise terms on the resolved scales. For pattern-forming systems, one is, however, interested in the patterns on the unresolved scales, 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 semiphenomenological 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 the 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 HSS. 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 Escherichia coli which has emerged as a paradigmatic model system for the study of pattern formation in cell biology (16–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. 1). 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 (10, 21–23), and “homogeneous pulsing” (24–26), as well as quasi-stationary labyrinths, spots, and mesh-like patterns (10, 27). Theoretical analysis of mathematical models has led 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 in pattern formation.

Fig. 1.

(A) Schematic illustration of the Min protein reaction network. (B) Wedge geometry with a membrane surface at the bottom plane (z = 0) and bulk height H(x) increasing linearly along the x direction. (C) Snapshot of the membrane density of MinD, obtained by numerically simulating the Min dynamics Eqs. – 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.

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. 1). While there are many distinct ways to introduce large-scale spatial heterogeneities into the system, for example, 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. 1). 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 (8, 14, 15), 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, that is, the “hydrodynamic variables” of the system. The key to this reconstruction is correlations between the regional pattern characteristics and instantaneous, regional dispersion relations, calculated from the instantaneous regional mass densities. Over time, these masses change due to diffusive redistribution, resulting in qualitatively different regional dispersion relations that indicate the local pattern type in the system. This reconstruction of small-scale features (on unresolved scales), together with a coarse-grained description for the mass redistribution dynamics on large scales allows us to understand and predict the long-term temporal evolution of the system. A major advantage of our approach is that it is based on a linear theory and therefore conceptually and technically simple to apply.

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

Results

The Min Protein System in Wedge Geometry.

Mathematically, the Min protein dynamics is described by bulk surface coupled reaction–diffusion equations, which describe the concentrations of cytosolic proteins MinD-ATP, MinD-ADP, and MinE, , in the bulk volume , and the concentrations of membrane-bound MinD and MinDE complexes, , on the surface . For the wedge geometry, in spatial coordinates , we place the membrane surface (with lateral dimensions ) in the plane at and let the bulk height vary as a linear ramp from H0 to H1 along the x direction (Fig. 1).

The dynamics of bulk components is governed by the equationwhere D denotes the bulk diffusion constant, and the matrix describes nucleotide exchange of MinD in the bulk. The dynamics of membrane components is constrained to the membrane surface and takes the formwhere D is the membrane diffusion constant and is the surface Laplacian. The membrane reactions r, which comprise attachment, detachment, and recruitment processes of Min proteins, are specified in .

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 ). 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 is the total cytosolic MinD concentration; and denote the mean on the surface and in the bulk, respectively; and and are the total surface area and bulk volume (see ).

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 (Movie S1 and Fig. 1). 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 ().

Experimental Implementation.

We tested our theoretical prediction on this multiscale dynamics in an experimental system consisting of a wedge-shaped microfluidic flow chamber (Fig. 2). The bottom and top surface of the wedge were covered with a supported lipid bilayer consisting of 1,2-dioleoy- lphosphatidylglycerol: 1,2-dioleoylphosphatidylcholine (DOPG: DOPC) (30:70%) which mimics the natural membrane composition of E. coli (29). The length of the wedge was typically about 8 mm to 14 mm, and the width was about 3 mm to 4 mm. The bulk height range was approximately 2 µm to 50 µm (Fig. 2). Min proteins were distributed in the chamber by rapid injection of a solution containing 1 µM MinD and 1 µM MinE (including 10% fluorescently labeled MinD and MinE proteins for visualization), together with 5 mM adenosine 5-triphosphate (ATP) and an ATP regeneration system (28).

Fig. 2.

Experimentally observed Min patterns in a wedge flow cell. (A) Schematic presentation of the experimental setup. Both the bottom and the top surfaces (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. (Top) Shown is a merge of MinD (green) and MinE (red) channels. Bottom shows a kymograph for intensities taken along the center line in Top.

Fig. 2 shows a snapshot of Min protein patterns along the bottom surface of the wedge geometry 30 min 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. 2 and Movie S3): At regions of low bulk height (approximately between 2 µm and 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 (). 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 (), 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 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 a few hundred micrometers in size) that oscillate with time ( 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 surfaces of the wedge were covered with a supported lipid bilayer. In recent work, it has been shown that membrane-to-membrane cross-talk (i.e., between top and bottom surfaces) 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 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 obtaining 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–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 (1D) domain, these eigenmodes are simply Fourier modes. Unfortunately, in the wedge geometry with bulk–surface coupling, the eigenmodes cannot 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 an approach that overcomes these restrictions.

Instantaneous, Regional Dispersion Relations Predict Patterns.

The analysis of pattern-forming systems usually starts with calculating the HSS solutions and performing a linear stability analysis around these states. This yields a dispersion relation that informs about the growth rate of small spatial perturbations with a certain wavenumber q. However, the dispersion relation is generally only informative in the vicinity of the HSS (1, 6), and is 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 semiphenomenological approach where we generalize the concept of dispersion relations.

Let us consider the wedge as dissected into a collection of 2D 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 (Fig. 3 ). Next, for each slice and at each point in time, we calculate instantaneous total densities of MinD and MinE, averaged over the slice length (). The average total densities, together with the local bulk height H(x), then serve as parameters for the regional dispersion relation in each slice,

Fig. 3.

(A) Rectangular geometry with membrane at the bottom edge representing a slice through the 3D 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 , from the numerical simulation data. From these slice-averaged total densities, we can then calculate the corresponding local HSS and its dispersion relation. (C) Dispersion relation with fastest-growing mode and right edge of the band of unstable modes indicated. The ratio has been empirically found to correlate with the type of fully developed pattern, with a sharp transition from chaotic patterns for to ordered patterns for . Close to the transition, standing waves are found, while traveling waves form for larger ratios (14). (D) Mode ratio 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, SW, and TW.

which is straightforward to determine because the slice represents a rectangular geometry (14, 23, 28) (Fig. 3 and ). While the bulk height H(x) varies linearly in space, the average total densities 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: . 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 HSSs, 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 2D 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, and , and the bulk height H, the system exhibits a variety of different patterns on the membrane, such as chaos, standing waves, and traveling waves (14, 28). Moreover, a careful analysis of numerical simulations has interestingly revealed a strong one-to-one correlation between the dispersion relation and the fully developed patterns in the highly nonlinear regime (14): A commensurability criterion between the unstable mode with the shortest wavelength and the fastest-growing mode has been found that determines the pattern type (Fig. 3 ). In short, it has been shown that coincides with the regime of chemical turbulence (spatiotemporal chaos), whereas, for , the system exhibits ordered patterns (standing/traveling waves). Standing wave patterns are found close to the commensurability transition , while traveling waves are found farther 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 as a function of slice position x and time t (Fig. 3 ). The resulting pattern type prediction is shown in the space–time plot (kymograph) in Fig. 4. Fig. 4 shows the ratio as a function of slice position x for a set of representative times (Fig. 3). The pattern type prediction (Fig. 4) is then obtained from these ratios via the mapping shown in Fig. 3 .

Fig. 4.

(A) Kymograph showing the pattern type prediction from the commensurability criterion (Fig. 3). The green line shows where , indicating the transition from chaotic to ordered patterns. Green arrows mark the position for the times indicated by dashed white lines. (B) Plots of the mode ratio , 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 , predicting the global emergence of standing waves (see last row). (C) Snapshots of the membrane patterns (MinD density; Fig. 1) from the full numerical simulation. The green dashed line indicates . Note the SW patterns found near . Their fronts are aligned along the bulk height gradient such that the sequence of wave nodes lies on lines of constant bulk height. (D) Machine learning–based pattern classification using ilastik (35) (see ).

We find that this prediction correlates well with the patterns observed in the full numerical simulation (Fig. 4 and Movie S8). In particular, the temporally changing position , marking regions where (indicated by the green arrows and dashed lines in Fig. 4 ), agrees with the position along the wedge where traveling wave patterns transition to chaotic patterns. In the vicinity of , we observe a band of standing waves as expected from the “commensurability criterion” (14). While the transition from chaos to order at is sharp, we were not able to identify a sharp criterion for the transition from standing waves to traveling waves. Accordingly, we use a smooth color gradient to indicate the pattern type prediction near (Figs. 3 and 4). Since the ratio and, with it, are entirely determined by the slice-averaged masses , 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 (Fig. 4 ). The reason for these deviations lies in the model parameters, which were chosen such that the entire domain is near the critical mode ratio 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 , we provide additional results where the parameters were chosen such that the mode ratio is deep in the traveling wave regime ( for late times. In this case, we obtain an excellent agreement between our predictions and the patterns observed in the numerical simulations ().

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. Similarly to the previous section, we consider here the redistribution of mass between slices along the wedge (Fig. 3). Since membrane diffusion is two orders of magnitude slower than bulk diffusion, it may be neglected, such that redistribution of protein mass between slices is governed by bulk diffusion alone (),for . 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 redistributes mass. This can be seen by rewriting Eq. in the form of a continuity equationwith the diffusive fluxes given by . Since the area of slices increases along the positive x direction, the diffusive fluxes on the right-hand side of Eq. 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 appear on the right-hand side.

We are interested in the dynamics of 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 the y direction is well approximated by the instantaneous HSS in a slice. These steady-state concentrations only depend on the slices’ bulk height H(x) and the slice-averaged total densities . 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 HSSs 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 . To regularize the dynamics, c is not set identical to but relaxes toward it on a fast timescale (see for details).

Given the initial densities , one can numerically solve the reduced dynamics Eq. to predict the entire time evolution of the slice-averaged masses and hence the dispersion relation at each point along the x direction. Fig. 5 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 (Fig. 4). The main difference from the full numerical simulations is a slight quantitative deviation in the timescale, where the dynamics predicted by Eq. 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 (Figs. 4 and 5). This is because the chaotic regime is rather narrow and close to the regime for which the dispersion relation predicts no instability (Figs. 3 and 4). 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.

Fig. 5.

(A and B) Kymographs showing the 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. 4.

Fig. 5 compares 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 (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 (11, 37). 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 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 (8, 15, 36, 38).

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

An example for this third scenario is the Aranson–Tsimring model for pattern formation in vibrated granular media (39) (see 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. in ). Near the onset of pattern formation, this coupling gives rise to localized patterns that have been studied using amplitude equations (30, 32, 33). Fig. 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 ρ, 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 ρ (Fig. 6). This demonstrates the idea of regional dispersion relations in a nutshell. Moreover, it shows that this approach gives rise to qualitatively similar insights such 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.

Fig. 6.

Regional dispersion relations predict localized patterns in the Aranson–Tsimring model Eq. . (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 SD 10). (C) Representative dispersion relations in the stable and unstable regimes. Domain size: 100 × 50; see for model details and remaining parameters.

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, that is, local symmetry breaking (40–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–44). The dynamic interplay of mass redistribution and orientational order has been shown to give rise to the coexistence of different macroscopic orders (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, nontrivial 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 nonlinear 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.

Materials and Methods

Mathematical Model.

We adopt the Min “skeleton model” introduced in refs. 5, 45, and 46. which is known to qualitatively reproduce Min patterns in vivo and in vitro (5, 28, 46). The governing equations are given in Eqs. –. The membrane reactions arewith

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. Coupling between cytosol and membrane is established by reactive boundary conditions at the membrane (Eq. ). The boundary fluxes are given bywhich follows from mass conservation. For analytical calculations, we adapt the following change of variables, as it is more convenient: We describe the bulk dynamics of MinD in terms of the variables and ; that is, in this case, one defines the bulk concentration vector . The membrane reaction in Eq. is then slightly modified by substituting , and the boundary fluxes are given by

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

Table 1.

Min model parameters

Parameter	Symbol	Value
Bulk diffusion	Dc	60  μm^2·s^−1
Membrane diffusion	Dm	0.013  μm^2·s^−1
Average total MinD density	n¯_D	665 μm^−3
Average total MinE density	n¯_E	410 μm^−3
Attachment rate	k_D	0.065 μm·s^−1
MinD recruitment rate	k_dD	0.098  μm^3·s^−1
MinE recruitment rate	k_dE	0.126  μm^3·s^−1
MinDE dissociation rate	k_de	0.34 s^−1
Nucleotide exchange	λ	6 s^−1

Numerical Simulations.

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

HSS and Dispersion Relation.

The HSS concentrations are obtained from the stationary solutions of Eqs. – together with the mass conservation condition Eq. ,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 . For a thorough presentation of the linear stability analysis of the Min system in a 2D rectangular geometry, we refer to the SI Appendix 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:where the membrane surface area and the bulk volume for the wedge geometry are explicitly given by and .

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 average these densities along the y direction to obtain the slice-averaged total densities . Note that the length of a slice is much larger than the typical pattern wavelength, which also permits approximation of the slice-averaged mass at the membrane by the vertically averaged mass: (see ref. 14). This is because the local deviations 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. . 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 Eqs. and ,where we used the reactive boundary condition Eq. 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 , where 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. 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 . 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 (47, 48), one can neglect membrane diffusion to arrive at the result shown in Eq. . For completeness, note that Eq. (without membrane diffusion) can be recast aswhich is the form given in Eq. .

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. These data were then down-sampled and averaged over slices to yield an x–t space time map of pattern probabilities. To render the kymograph in Fig. 4, 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 ref. 39. This model, which we call the 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. , while gradients in the amplitude drive mass redistribution away from high-amplitude regions (second term in Eq. ). 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 refs. 32 and 33. Moreover, in ref. 30, it was shown that the system Eq. appears as the amplitude equation for a mass-conserving version of the classical Swift–Hohenberg–Turing equation (6, 49). The reason for this is that the conserved density appears as a second hydrodynamic variable in addition to the pattern amplitude.

A linear stability analysis shows that the system Eq. has a short-wavelength instability when and , where ρ0 denotes the average density. Following ref. 33, we set parameters . Localized patterns are found near the instability threshold, so we set for the simulation shown in Fig. 6 and Movie S9.

Preparation of the Wedge Flow Cell.

The microfluidic wedge chambers were prepared using two rectangular coverslips (bottom one of dimensions 22 mm × 50 mm, and top one of dimensions 5 mm × 30 mm). Close to one of the short edges of the top glass, a tiny inlet hole was drilled using a sandblaster. Coverslips were cleaned in 1 M potassium hydroxide (KOH) for 1 h followed by a methanol bath for 10 min in a sonicator bath. Surfaces of the coverslips were activated with oxygen plasma for 20 s, using oxygen plasma PREEN I (Plasmatic System, Inc.) with a O2 flow rate of 1 standard cubic feet per minute (SCFH). Furthermore, a small polydimethylsiloxane (PDMS) slab with a 0.3-mm hole was attached onto 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 prevented a collapse of the top and bottom slides (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 (). 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 TrisHCl (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 entire 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 fluorescence recovery after photobleaching and photoactivation (FRAPPA), illumination and detection system Andor Revolution and Yokogawa CSU X1, electron multiplying charge coupled device (EM-CCD) Andor iXon X3 DU897 camera, motorized x–y stage and a z-piezo stage, using a 20× objective (UPlansApo, numerical aperture 0.85, oil immersion). Imaging of MinD-Cy3 and MinE-Cy5 was performed with laser spectral lines at 561 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 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 (maximum 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 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 become comparable with the values at the center of the image.