Stochastic switching between multistable oscillation patterns of the Min-system

The spatiotemporal oscillation patterns of the proteins MinD and MinE are used by the bacterium E. coli to sense its own geometry. Strikingly, both computer simulations and experiments have recently shown that for the same geometry of the reaction volume, different oscillation patterns can be stable, with stochastic switching between them. Here we use particle-based Brownian dynamics simulations to predict the relative frequency of different oscillation patterns over a large range of three-dimensional compartment geometries, in excellent agreement with experimental results. Fourier analyses as well as pattern recognition algorithms are used to automatically identify the different oscillation patterns and the switching rates between them. We also identify novel oscillation patterns in three-dimensional compartments with membrane-covered walls and identify a linear relation between the bound Min-protein densities and the volume-to-surface ratio. In general, our work shows how geometry sensing is limited by multistability and stochastic fluctuations.


Introduction
The bacterial Min-proteins are a well studied example of a pattern-forming protein system that gives rise to rich spatiotemporal oscillations. It was discovered as a spatial regulator in bacterial cell division, where it ensures symmetric division by precise localization of the divisome to midcell [1]. The dynamic nature of this protein system was demonstrated by live cell imaging in E. coli bacteria, where these proteins oscillate along the longitudinal axis between the cell poles of the rod-shaped bacterium, forming so-called polar zones [2][3][4][5].
Most bacteria use a cytoskeletal structure, a so-called Z-ring, for the completion of bacterial cytokinesis [6]. This Z-ring self-assembles from filaments of polymerized FtsZ-proteins, the prokaryotic homolog of the eukaryotic protein tubulin [7], which serve as a scaffold structure for midcell constriction and the eventual septum formation in the midplane. If successful, this process creates two equally sized daughter cells with an identical set of genetic information [8]. A necessary prerequisite for successful symmetric cell division is hence the targeted assembly of FtsZ towards midcell. In E. coli cells this is mediated by two independent mechanisms, nucleoid occlusion and the dynamic oscillation of the MinCDE proteins [5,9,10]. While nucleoid occlusion prevents division near the chromosome, the Min-system actively keeps the divisome away from the cell poles through the MinC-protein acting as FtsZ-polymerization inhibitor [5]. The characteristic pole-to-pole oscillations create a time-averaged concentration gradient with a minimal inhibitor concentration of MinC at midcell, suppressing Z-ring assembly at the cell poles [3][4][5]. Although MinC is indispensable for correct division site placement, it acts only as a passenger molecule, passively following the oscillatory dynamics of MinD and MinE [2,3,11].
On the molecular level, the oscillations emerge from the cycling of the ATPase MinD between a freely diffusing state in the cytosolic bulk and a membrane-bound state, induced by its activator MinE under that are considered here with computer simulations. While geometry A uses a flat membrane patch, similar to flat patterned substrates [26], geometry B corresponds to microfabricated chambers with an open upper side [27,28]. Geometry C corresponds to the cell sculpting approach [29,30].
Like for other pattern-forming systems, the theory of reaction-diffusion processes offers a suitable framework to address the Min-oscillations from a theoretical point of view [31][32][33][34]. Many theoretical models have been proposed to unravel the physical principles behind this intriguing self-organizing protein system and to explain the origin of its rich spatiotemporal dynamics. While all of them rely on a reaction-diffusion mechanism similar to the Turing model, they differ severely in their details. The first class of mathematical models used an effective one-dimensional PDE-approach and relied strongly on phenomenological nonlinearities in the reaction terms [35][36][37]. Although all of them successfully gave rise to pole-to-pole oscillations, they did not allow a clear interpretation of the underlying biomolecular processes and were not in agreement with all experimental observations, such as MinE-ring formation and the dependence of the oscillation frequency on biological parameters.
The next advance in model building was the focus on the decisive role of MinD aggregation and the relevance of MinD being present in two states (ADPand ATP-bound) [37,38]. Very importantly, this highlighted the interplay between unhindered diffusion with a nucleotide exchange reaction in the bulk as a delay element for MinD reattachment [39]. Subsequent models shared a common core framework but still differed strongly in the functional form of the protein binding kinetics and the transport properties of membrane-bound molecules [40,41]. The main difference between the more recent models was the dimensionality, ranging from onedimension [35][36][37]42] to two [26] and three [39,[43][44][45][46][47]. Moreover, the models can be classified as deterministic PDE-models [24,26,30,[35][36][37][38][39][47][48][49][50] or using stochastic simulation frameworks [22,[43][44][45][46][50][51][52], and whether they neglected membrane diffusion [35,36,39,43,48] or not. While some models contained higher than second order nonlinearities in concentrations of the reaction terms [24], it is the prevailing opinion to rely on at most second order nonlinearities, allowing for a clear interpretation in terms of bimolecular reactions. Following the same line of thought, a strong effort was made to distill a minimal system that explains the oscillation mechanism without the necessity of spatial templates or prelocalized determinants [48] and neglecting secondary processes like filament formation [45,46,48].
The most influential minimal model for the Min-system has been suggested by Huang and coworkers [39]. It has been further simplified by discarding cooperative MinD recruitment by MinDE complexes on the membrane [44,49], allowing a clear view on the core mechanisms: the cycling of MinD between bulk and membrane, cooperativity of MinD-recruitment and diffusion in bulk and along the membrane. Using the minimal model, it has been shown that the canalized transfer of proteins from one polar zone to the other underlies the robustness of the Min-oscillations [44,49]. Because the deterministic variants of the minimal model [39,49] do not allow us to address the role of stochastic fluctuations, a stochastic and fully threedimensional version has been introduced to study the effect of stochastic fluctuations in patterned environments [52]. For rectangular patterns of 5 μm×10 μm, it was found that the system can be bistable, with transverse pole-to-pole oscillations along the minor and longitudinal striped oscillations along the major axis, respectively. In this early work, it was observed that the stable phase emerged depending on the initial conditions and that sometimes switching occurred, but the statistics was not sufficiently good to observe switching in quantitative detail. Indeed such multistability has been observed experimentally in sculptured cells over a large range of cell shapes [30] and the deterministic minimal model has been used to explain the relative frequency of the different oscillations patterns for a given shape using a perturbation scheme [47]. However, as a deterministic model, this approach was not able to address the rate with which one pattern stochastically switches into another.
Here we address this important subject by using particle-based Brownian dynamics computer simulations. Compared to earlier work along these lines [52], we have developed new methods to efficiently simulate and analyze the switching process. We find excellent agreement with experimental data and measure for the first time the switching time of multistable oscillation patterns. We also use our model to confirm that it contains the minimal ingredients for the emergence of Min-oscillations. In addition, we use our stochastic model to investigate the three-dimensional concentration profiles in different geometries and in particular the role of edges in membrane-covered compartments. We identify novel oscillation patterns in compartments with membrane-covered walls and find a surprisingly simple (linear) relation between the bound Min-protein densities and the volume-to-surface ratio, which might be relevant for geometry sensing by E. Coli cells.

Reaction-diffusion model and parameter choice
For the particle-based simulation, we use the reaction scheme of the minimal model for cooperative attachment [39,47,49,52]. The model uses the following interactions between Min-proteins and the inner bacterial membrane (schematically shown in figure 1(A) [30] and set C in [47,49]).

Simulation algorithm
We use custom-written code to simulate the stochastic dynamics of the Min-system with very good statistics. For all simulations we use a fixed discrete time step of D =t 10 s 4 . During every time step each particle is first propagated in space. Thereafter every particle can react according to the previously introduced Min reaction scheme (1a)-(1e). The movement of both free and membrane-bound particles is realized through Brownian dynamics. Individual molecules are treated as point-like particles without orientation. Therefore we can monitor the propagation separately for each Cartesian coordinate. During a simulation step of Dt the displacements of the diffusing particles with diffusion constant D are drawn from a Gaussian distribution with standard deviation s = D D t 2 x [53] such that where p X G is the probability distribution of X G . The same update step is used for the y and z direction. Free particles in the bulk of the simulated volume undergo three-dimensional diffusion with reflective boundary conditions at the borders of the simulation volume. The membrane-bound particles perform a two-dimensional diffusion on the membrane with a much smaller diffusion constant D bound (compare table 1). Membranebound particles are allowed to diffuse between different membrane areas that are in contact with each other. Table 1. Parameters sets used for simulations of the Min-system. Set A is the main parameter set used here following [39,52]. Parameter set B is taken from [30] and C from [47,49]. The different reactions in the Min reaction scheme (1a)-(1e) can be classified into three different types (more details on the corresponding implementations are given in the supplementary information). The first type considered here are first order reactions without explicit spatial dependence. The conversion of MinD ADP to MinD ATP and the unbinding of the MinDE bound complex from the membrane are of this type. Such reactions are treated as a simple Poisson process. For a reaction rate κ, the probability to react during a time step Dt is given by The second type is also a first order reaction, but with confinement to a reactive area at a border of the simulated volume. The membrane attachment of MinD ATP proteins is a reaction of this type. For a given reaction rate σ, we implement these reactions by allowing particles that are closer to the membrane than m = d 0.02 m to attempt membrane attachment with a Poisson rate k s = d. This results in a reaction The last reaction type is a second order reaction between free and membrane-bound particles. The cooperative recruitment of cytosolic MinD ATP and MinE to membrane-bound MinD are reactions of this third type. In our simulation we adopt the algorithm implemented in the software package Smoldyn, which has been used earlier to simulate the Min-system [52]. This algorithm is based on the Smoluchowski framework in which two particles react upon collision [54]. However, the classical treatment by Smoluchowski only considers diffusion-limited reactions and therefore assumes instantaneous reactions upon collision. In order to take finite reaction rates into account, one imposes a radiation boundary condition [55][56][57]. From the diffusion constant D, the reaction rate σ and the simulation time step Dt , a reaction radius r σ is calculated [58]. Whenever a freely diffusing particle comes within the distance of r σ to a membrane-bound particle, the free particle reacts. For intermediate values of Dt (such as the time step of -10 s 4 that we use for the Min-system) the value of is obtained numerically [58]. Those numerical values are taken from the Smoldyn software. For example, for parameter set A the reaction radius for the rate s dD is m = In our simulations we use rectangular reaction compartments. We considered three different membrane setups as illustrated in figure 1(B). To mimic in vitro experiments, where substrates or open compartments are functionalized with a membrane layer [26][27][28], we place the reactive membrane at the bottom (geometry A) or at the side walls and the bottom of the simulation compartment (geometry B). To simulate rectangular shaped E. coli cells, inspired by the cell sculpting approach from [30,47], fully membrane-covered volumes are used (geometry C). We refer to the long side of the lateral extension as the major or the x-axis, and the smaller side as minor or y-axis, and accordingly consider the compartment height to extend in the z-direction, aligning the rectangular geometry perpendicular with the coordinate frame.
In our simulations we investigate a wide range of compartment dimensions. For a simulation box with a length of 10 μm, width of 5 μm and height of 0.5 μm, we use 6003 MinD ATP particles, 6003 MinD ADP particles and 3124 MinE particles as initial condition [52]. These particle numbers amount to a total MinD concentration of m 0.797 M and a MinE concentration of m 0.207 M. For other simulation compartment sizes we scale the particle numbers linear with the volume, since in experiments E. coli bacteria typically have a constant Minprotein concentration [30].

Identification of oscillation modes
In our simulations of the Min-system different oscillation patterns emerge along the major or minor axis of the simulation compartment. In order to analyze the frequency of different modes and the stability of the oscillations in the large amount of simulation data, an oscillation mode recognition algorithm is needed. Therefore we monitor the MinD protein densities at the poles of the different axes over time. To determine the axis along which the oscillation takes place, we compare the Fourier transformation of the normalized densities over time (r t i where i denotes the discretized time resolution). If an oscillation takes place, there is a dominant peak in the Fourier spectrum and the overall maximal amplitude of the Fourier spectrum is significantly higher than the one from the non-oscillating axis. The same Fourier spectrum is also used to determine the oscillation period T. To differentiate between pole-to-pole oscillations and striped oscillations of a given axis in the system, we extract the phase difference between the density oscillation at the poles of the cell.
When identifying switching events, one has to be careful because they have to be identified in local time windows and stochastic fluctuations might lead to temporal changes that might be mistaken to be mode switches. We therefore smoothen the data. In detail, we calculate the convolution C i between the densities over time and a Gaussian time Here τ is the time between successive density measurements and we set w = 100 s as width of the time window. The current oscillation mode is now determined from the convoluted densities C i and assigned to the time t i . In this way, only switches are identified that persist for a sufficiently long time.

Oscillation patterns in geometry A
First we investigated the oscillations that emerge in geometry A with parameter set A using a rectangular simulation volume with dimensions m m ḿ10 m 5 m 0.5 m (x y z , , ). With this particular choice the width of the system approximately matches the typical length of wild-type E. coli cells and the length of the system corresponds to the length of a grown E. coli cell which can roughly double in length before septum formation and division. As shown by the kymographs in figure 2 and in agreement with the results of Hoffmann and Schwarz [52], in our simulations two different oscillation modes occur (compare also supplemental movie S1). Note from the color legend that dark and light colors correspond to low and high concentrations, respectively, as used throughout this work. In the first mode the Min-proteins oscillate along the minor y-axis from one pole to the other (pole-to-pole oscillation). In the second mode the proteins oscillate along the major x-axis between the poles and the middle of the compartment (striped oscillation). The system stochastically switches between the two modes, sometimes via a short oscillation along the diagonal of the compartment. The mode switching behavior of the Min-system in large volumes is in agreement with the experimental results of Wu et al [30] and cannot be analyzed completely with conventional PDE-models of the Min-oscillations because they do not account for the noise in the system leading to the stochastic switch.  kymographs of the free particles (middle row in figure 3(A)) are averaged over all heights and have the inverse shape of the corresponding kymographs of the bound particles (top row in figure 3(A)). Where the density of bound particles is high, the density of free particles is low and vice versa. During the simulations the spatial density differences of both MinD and MinE are higher for membrane-bound particles than for free particles in the bulk. Therefore in the two bottom kymographs in figure 3(A), which are showing total particle densities, the pattern of membrane-bound particles is dominant.
The kymographs of the striped oscillation in figure 3(B) have the same structure as the ones of the pole-topole oscillations. However, the edges of the bound Min-protein clusters in the kymographs, that indicate the traveling Min-waves, are curved, in contrast to the straight lines that we see for the pole-to-pole oscillations as shown in figure 3(A). The time-averaged density profiles of MinD proteins for the pole-to-pole and striped oscillations are shown in figures 3(C) and (D), respectively. As expected the density of the proteins is minimal between the oscillation nodes of the emerging standing wave patterns.
It is highly instructive to compare the time evolution of the MinD and MinE protein densities. In figure 4(A) we plot the time evolution of the particle densities of the transverse pole-to-pole oscillation at a fixed position m = y 4.9 m, which is at the edge of the minor axis along which the oscillation takes place. The shape of the transient density profiles is similar to experimentally observed density profiles of traveling Min-protein waves on flat membrane surfaces [24,25]. The period of both oscillations modes was =  ( ) T 33.8 0.1 s, which is in agreement with the results of Huang et al [39].
To analyze the influence of the bulk volume on the oscillations, we have also monitored how the MinD and MinE densities changes at different heights above the membrane. For this we again use geometry A, but now with a bottom surface of only m ḿ 6 m 3 m, which robustly produces longitudinal pole-to-pole oscillations along the major axis. For a compartment height of 4 μm the volume was divided in three layers, and for each layer the mean density on one side of the major axis was plotted over time as shown in figure 4(B). We see that the largest density changes take place in the layer directly above the membrane, however, the oscillations persist up to the top layer even for the highest compartment with m = z 4 m. Strikingly, the MinD-density is always much higher than the one of MinE. Furthermore we reduced the compartment height to 0.2 μm. Interestingly, the pole-to-pole oscillations were still present. This implies that the bulk of the simulation volume has only a mild effect on the oscillations in geometry A, despite the fact that there are appreciable density variations in the bulk.

Essential model elements
We next checked that our model indeed includes the essential elements for the emergence of Min-oscillations. Because the MinD-switching between bulk and membrane is clearly indispensible, here we consider the relevance of diffusion on the membrane and of cooperative recruitment. All simulations in this section are carried out using geometry A in a m m ḿ5 m 2.5 m 0.5 mcompartment geometry. We chose this geometry since we expect it to give rise to stable longitudinal pole-to-pole oscillations and hence can closely monitor any deviations from this default oscillation mode. All other parameters were kept fixed and we used parameter set A for this test. where stable pole-to-pole oscillations emerge along the major axis. By further increasing the membrane-diffusion-coefficient the pole-to-pole oscillations still robustly emerge but increasingly smear out. This behavior is most clearly illustrated by looking at the time-averaged density profiles as shown in figure 5(C). The inset in figure 5(C) shows that without membrane diffusion, no oscillations emerge at all (square symbol), while a slow diffusivity as used in parameter set A seems optimal.
In a similar fashion we analyzed the influence of the cooperative membrane recruitment of MinD. Figure 5(B) shows three kymographs for no cooperative recruitment (s = 0 dD ), the default value from parameter set A and a two-fold increase, respectively. Without cooperativity in this process no oscillations emerge at all. The second panel shows again the default oscillation mode while already a two-fold increase also leads to unstable behavior without any patterns emerging. This sensitivity with respect to the cooperative recruitment rate s dD is also summarized in figure 5(D), which illustrates that only in a vary narrow range stable oscillations emerge. Although a complete parameter scan is out of the question at the current stage for reasons of computer time, we conclude that stable Min-oscillations emerge only for certain parameter values and that parameter set A performs very well in this respect.

Edge oscillations in geometry B
As we have seen in the preceding section, the simulation compartment height has little influence on the emerging oscillation modes when the membrane only covers the bottom area of the simulation compartment (geometry A). However, when we extend the membrane to cover bottom and side walls of the compartment (geometry B), we find that the oscillation modes change with increasing wall height. First, we consider a compartment that only exhibits longitudinal pole-to-pole oscillations along the major axis (wall height of 0.5 μm and bottom area of m ḿ 5 m 2.5 m). In contrast to the flat membrane geometry A, here the MinD protein density decreases in the vicinity of the membrane edges between the side walls and the bottom area, as shown in figure 6. This is due to the decreased volume per membrane area ratio in the vicinity of the membrane corners, leading to a decreased density of bound MinD proteins in these regions. This effect is clearly visible by comparing the time-averaged density profile in figure 6 with the one shown in 3(C).
The full oscillation pattern is shown in 7(A) for the lowest height value (black label) and resembles a pole-topole oscillation for geometry A, but now between the two opposing walls along the major axis. Now we increase the wall height in increments of 0.5 μm. The oscillation changes at a wall height of m 3.5 m to a new oscillation mode (green label, compare also supplemental movie S4). There the proteins start to oscillate between the middle of the bottom area and the upper edges of the walls. This newly identified edge oscillation with one polar zone at the bottom and one top polar zone at each of the four side walls persists until the wall height reaches 9.5 μm. Thereafter the oscillation mode changes to a striped edge oscillation along the walls (yellow label).
In figure 7(B) we show results for a simulation volume that gives rise to striped oscillations along the major axis (bottom area of the length of 9 μm and width of 4 μm). At a wall height of m = z 0.5 m the simulation gives rise to longitudinal striped oscillations along the major axis or transverse pole-to-pole oscillations along the minor axis. This is the same kind of bistability that we already observed in geometry A using a compartment of For an edge oscillation the free MinD protein densities in the bulk of the compartment along one of the bottom area axis (y-axis) and along the walls (z-axis) are shown in figure 8(A). We see that spatial oscillations in the bulk only take place along the z-axis. On the y-axis kymograph (figure 8(A) top kymograph) the change of the density only takes place along the temporal axis. This corresponds to the change of total numbers of free and bound MinD proteins during the oscillation and no spatial oscillation takes place along the y-axis. On the z-axis kymograph ( figure 8(A) bottom kymograph) we observe a spatial pole-to-pole oscillation. In figure 8(B) we show the densities of the free MinD proteins in the bulk of the compartment during the large pole-to-pole oscillation for m = z 4 m. We see that the bulk proteins oscillate only along one axis (here the y-axis) and no spatial oscillations occur along the walls (z-axis).
From our study of geometry B, we conclude that in non-flat membrane geometries the oscillations of the Min-proteins do not only take place along the membrane, but are rather defined by the canalized transfer of the proteins through the bulk. This shows the importance of three-dimensional simulations for non-flat membrane geometries in order to determine the self-organized oscillation modes.

Geometrical determinants of bound particle densities
We have also monitored the amount of membrane-bound MinD and MinE proteins and compared it to the total amount of proteins in the system. Interestingly, we observed that these two quantities have a linear dependence µ N N bound total as shown in figure 9(A). Surprisingly, this relation seems to hold for all geometries and dimensions, as indicated in figure 9(A) by the different symbols. Since in our simulations we keep density constant, the total amount of proteins scales linearly with the compartment volume ( µ N V total ). Overall, we conclude the following relation for the mean bound Min-protein density where A denotes the total reactive membrane area. Thus the mean bound Min-protein density increases linearly with the volume-to-area ratio, as verified by figure 9(B). Strikingly, we again observe a dramatic difference between MinD and MinE. Although this scaling is the same for both, MinE is much closer to being completely bound (gray lines), indicating that once a stable oscillation emerges, MinE is almost completely depleted from the bulk. We next assessed the physiological relevance of these observations. In order to investigate how the volumeto-area ratio V/A changes during cell growth, we approximate the shape of an E. coli bacterium by a spherocylinder of length L and width W and evaluate V/A analytically. Since E. coli bacteria mainly grow in length and stay constant in width [59], we show the evolution of V/A as function of the cell length L for various fixed cell widths in figure 10(A). One clearly sees that V/A remains nearly constant as a function of the cell length L. Recalling that our previous observation stated that the mean membrane-bound Min-protein densities r bound    scale linearly with V/A, this would suggest that living bacteria grow in a fashion that keeps both V/A and thus consequently r bound constant, which might be advantageous for the stability and robustness of the Minoscillations and related processes, such as formation of the FtsZ-ring. A similar qualitative behavior is observed if we instead of a spherocylinder would assume an ellipsoidal shape as is shown in figure 10(B).

Oscillation mode switching
Above we have seen that multistability is a recurrent phenomenon in the Min-system, both in geometries A and B. We next turn to a systematic investigation of the stochastic switching between two different oscillation modes. An oscillation mode transition of this type occurs frequently in a m m ḿ8 m 2 m 0.5 m compartment using geometry B. In this geometry the Min-system gives rise to both pole-to-pole and striped oscillations along the same (major) axis (kymograph in figure 11(A)). We measure the lifetimes (oscillation duration before a mode switch occurs) of the two modes during a 50 000 s long simulation trajectory. The resulting histograms of the lifetimes are shown in figure 12. We can approximate the mode switching as a Poisson process by assuming that the switching probability p(t) obeys . Fitting an exponential to the switching times histograms we obtain the switching rates of these processes. Here we neglect the measurements of short lifetimes below t = 100 s c since our oscillation mode detection algorithm can miss a mode transition if its lifetime is shorter. The switching rate for the pole-to-pole oscillation is = k 0.004 22 s p 1 and for the striped oscillation we find = k 0.004 06 s s 1 , thus the two modes seem to be equally frequent in this case. We have also performed the same analysis with identical compartment geometry for the two other parameter sets (set B and C) as presented in table 1. Parameter set B gives rise to stable longitudinal pole-to-pole oscillations along the major axis (kymograph in figure 11(B)) in agreement with the results from [30], while parameter set C shows a qualitatively similar switching behavior as parameter set A (data not shown) with frequent mode transitions.  To analyze the effect of the Min-protein concentration on the oscillation mode switching, we have performed the same simulation with increased Min-protein particle densities using again parameter set A. The resulting fractions of the two oscillation modes during the 50 000 s simulation trajectory are shown in figure 13, and their corresponding switching rates k p and k s are given in table 2. We note that the fraction of striped oscillations increases monotonously with increasing particle density. In agreement with this, the rate k s decreases with increasing particle density. In contrast, the rate k p does not show a systematic change and seems to fluctuate strongly. The shift to the striped oscillations suggests that the Min-oscillations can be used not only to sense geometry, but also to sense concentrations.
3.6. Phase diagrams for geometry C Wu et al have experimentally measured the fractions of different oscillation modes of Min-proteins in rectangular cell geometries of E. coli bacteria of various sizes with constant height [30]. In order to address these observations, we now turn to geometry C, which is a rectangular and fully membrane-covered compartment as sketched in figure 1(B). Throughout this section we keep the compartment height fixed at m = z 1 m, while we vary the length (major axis) between 2 and m 10 m and the width (minor axis) of the compartment between 1 and m 5 m in steps of m 1 m, respectively. To determine the oscillation mode fractions we calculate an ensemble average of the oscillation modes after 500 s simulation time. For each compartment size a sample of 20 independent simulations is used.
In the following analysis we only focus on the three main oscillation modes, which here are longitudinal pole-to-pole and striped oscillation (along the major axis) and transverse pole-to-pole oscillation (along the minor axis). Other oscillation modes as they have been reported in [30] were also observed using our framework, however, due to their low probability, a much higher sample size would be necessary to analyze their occurrence frequency with sufficient statistics. The results for the three oscillation modes are shown in figures 14(A)-(C). In general, each of the three oscillation modes dominates in one region of the phase diagram, but the transitions are fuzzy and therefore bistabilities occur. For compartments with length below m 7 m and width below m 4 m, only pole-to-pole oscillations are observed. Most of those oscillations occur along the major axis of the system. The transverse pole-to-pole oscillations along the y-axis emerge most frequently in compartments with quadratic bottom area. Increasing the width further increases the fractions of transverse pole-to-pole oscillations at the expense of longitudinal pole-to-pole oscillations. Increasing the length for a fixed width shows a sharp transition from longitudinal pole-to-pole to longitudinal striped patterns at around 6 μm. For compartments with length larger than m 7 m, the longitudinal pole-to-pole oscillations vanish almost entirely. In the region of both large long sides (around 7-10 μm in length) and large short sides (around 4-5 μm in width), the oscillation mode fractions are rather equally shared between longitudinal striped and transverse pole-to-pole oscillations, which is also in line with the bistability that we observed between these two patterns in the m m ḿ10 m 5 m 0.5 m compartment of geometry A as shown in section 3.1. Figure 14(D) summarizes these findings in a phase diagram that considers only the dominating mode (expect in the regions of clear bistability). Overall we find the

Conclusion
Using a stochastic particle-based simulation framework, we have investigated the stochastic switching between multistable oscillation modes of the Min-system in different three-dimensional compartment geometries.
Although it is well known that geometrical constraints have a strong impact on the dynamic oscillations of the Min-proteins, multistability and mode switching have only recently been investigated in more detail [47,52].
Our stochastic framework provides a suitable approach to address the question of oscillation mode selection and stability, since it naturally incorporates fluctuations due to finite copy numbers in the system. This allowed us to address the influence of the three-dimensional shape of the compartment and the boundary conditions, with a close match to existing experimental assays (flat supported bilayers [24,26], functionalized compartments [27,28], and cell sculpting [30,47]). For example, we addressed the role of compartment height and demonstrated the emergence of new oscillation modes with increasing height, underlining the importance of an explicit three-dimensional representation of the system. We also showed that diffusion long the membrane and cooperative recruitment are essential elements for the emergence of Min-oscillations. In general, particle-based stochastic computer simulations are a great tool for explorative research and in the future could be used to explore more details of the different scenarios that have been suggested for the molecular mechanisms shaping the Min-oscillations [42,50]. Our simulations demonstrated several features that might be related to the physiological function of the Min-system in E. coli. First we observe that there is a dominating length scale of m 5 m, which happens to be the natural length of an interphase E. coli cell. Second we found a linear relation between the density of membranebound Min-proteins and the volume-to-surface ratio, which tends to be constant during growth of E. coli. Third we observed that the relative frequency of competing oscillation modes depends on concentration, suggesting that the Min-oscillations can be used not only to sense geometry, but also concentration. Fourth, we found that stable oscillations strongly deplete MinE from the bulk. For the future, it would be interesting to study possible feedback between protein production and Min-oscillations.
For our simulations, we mainly used the established parameter set A from table 1 and achieved excellent agreement with experimental results regarding the relative frequency of the three main oscillation modes in different cell geometries [30]. Again the critical length scale around 5 μm plays an important role in transitions between different regimes, which then are the regions of high bistability. However, we also note that different parameter choices lead to different outcomes and that it would be interesting to perform an exhaustive exploration of parameter space to better understand how robustness and multistability depends on kinetic rates, diffusion constants and concentrations. In the future, such a complete scan might become possible by using GPU-code rather than the CPU-code developed here.
Most importantly, however, our stochastic approach allowed us to measure for the first time the switching rates between different competing oscillation patterns of the Min-system. This was done for parameter set A from table 1. Interestingly, parameter set C gave similar results in this respect, while parameter set B results in very stable oscillations without switching, in agreement with the experimental observations for the cell sculpting experiments [30]. In the future, it would be interesting to investigate more systematically how the effective barriers between two competing oscillation patterns depend on model parameters and compartment geometry. In general, the Min-system is an excellent model system to study not only geometry sensing, but also the role of spatiotemporal fluctuations in molecular systems.