Influence of cell shape, inhomogeneities and diffusion barriers in cell polarization models

In silico experiments bear the potential for further understanding of biological transport processes by allowing a systematic modification of any spatial property and providing immediate simulation results. Cell polarization and spatial reorganization of membrane proteins are fundamental for cell division, chemotaxis and morphogenesis. We chose the yeast Saccharomyces cerevisiae as an exemplary model system which entails the shuttling of small Rho GTPases such as Cdc42 and Rho, between an active membrane-bound form and an inactive cytosolic form. We used partial differential equations to describe the membrane-cytosol shuttling of proteins. In this study, a consistent extension of a class of 1D reaction-diffusion systems into higher space dimensions is suggested. The membrane is modeled as a thin layer to allow for lateral diffusion and the cytosol is modeled as an enclosed volume. Two well-known polarization mechanisms were considered. One shows the classical Turing-instability patterns, the other exhibits wave-pinning dynamics. For both models, we investigated how cell shape and diffusion barriers like septin structures or bud scars influence the formation of signaling molecule clusters and subsequent polarization. An extensive set of in silico experiments with different modeling hypotheses illustrated the dependence of cell polarization models on local membrane curvature, cell size and inhomogeneities on the membrane and in the cytosol. In particular, the results of our computer simulations suggested that for both mechanisms, local diffusion barriers on the membrane facilitate Rho GTPase aggregation, while diffusion barriers in the cytosol and cell protrusions limit spontaneous molecule aggregations of active Rho GTPase locally.


Introduction
Fundamental processes of living cells such as cell division, chemotaxis and morphogenesis depend on prior polarization and breaking of spatial symmetry. Spatial reorganization of membrane-bound and cytosolic proteins is required to establish an axis of polarity with a distinct direction (i.e. 'front and back') to guide directed processes. In these processes, cells have to adapt and react according to multiple and often conflicting cues of the environment. Rapid technical development makes it possible to observe these processes in great detail in vivo. Advances in imaging techniques such as total internal reflection fluorescence (TIRF) and confocal and electron microscopy [1][2][3][4] provide quantitative data for spatial cellular structures such as diffusion barriers on the membrane, organelles and membranes in the cytosol. To enhance the experiment-theory feedback loop, new computational tools have to be developed, to reconstruct the observed cell shapes, spatial inhomogeneities and structures at different levels of biological detail and mathematical complexity [5]. Quantitative modeling of intracellular asymmetries and inhomogeneities is essential for testing hypotheses, for obtaining an intuition of the process as well as for motivating new experimental studies.
The theoretical treatment of cell polarization has been devoted to biological model organisms such as yeast, fish keratocytes, Dictyostelium discoideum or neutrophils. Much recent work is concerned with reaction-diffusion (RD) models which are often described by partial differential equations (PDEs) in time and space as in [6][7][8][9][10][11][12][13][14][15]. Such models may vary greatly in complexity. Some of them are based on biochemical networks with many components [11,12,14,15] gathered from signaling pathway databases. Other models focus on only a few components to accurately describe cell behavior. Detailed reviews of models for yeast cell polarity can be found in [5,16].
The above sets of models pursue different goals. Since biochemical complexity can become overwhelmingly large, the experimental validation of generated models of high complexity is basically impossible [17,18] and an accurate spatial simulation becomes prohibitively costly. The alternative is to reduce signaling pathways to very few interacting constituents which still exhibit the specific behavior one is interested in. This simplified system can then be simulated subject to a spatio-temporal mathematical model which also includes (some) geometrical properties encountered in living cells. Ideally, in vivo or in vitro phenomena can be reconstructed, rates and concentrations can be compared and matched.
Since Meyers et al [19] showed the potential influence of cell shape and size on signaling, a growing number of models take into account the interaction of biochemical signaling and cell geometry [20,21]. More recently, some studies investigated the interplay of cell motility and signaling [22][23][24][25] in a 2D or 3D environment. There have been a lot of advances in computational frameworks for complex simulations in 2D and 3D. In particular, we want to mention VCell [26] dedicated to cell biology and DUNE [27] as a general purpose framework; the latter was used in this work. An interesting review on multiphysics models at different levels of complexity can be found in [28]. Examples where microscopic imaging and spatial modeling of cell shape and signaling are connected are shown in [29][30][31]. However, there are still a lot of aspects that need physical and computational analysis, among them inhomogeneous diffusion caused by organelles and other cellular inclusions.
Spatial structures such as organelles, diffusion barriers on the membrane and cell shape are expected to have a considerable influence on cell polarization. In [19,32,33], it was shown that high intracellular cytosolic gradients can be generated. Even though diffusion in the cytosol might be fast in theory, it can differ a lot from the effective diffusion rate due to obstacles such as large organelles, membrane stuctures or a crowded environment [34]. Therefore, it is important to incorporate intracellular gradients and diffusion as a spatial process in the cytosol [35][36][37]. Furthermore, membrane inhomogeneities are likely to have an effect on signaling as well [38].
To include the mentioned spatial properties in our modeling, we present an advanced computational approach to simulate RD models which take into account important aspects of spatial complexity. In our modeling approach spatial inhomogeneities such as organelles or bud scars are represented by an inhomogeneous diffusion coefficient or as inner boundaries in the computational domain. The RD model equations connect membrane and cytosol by coupling terms which represent reaction kinetics on the membrane and depend on the concentration of membranebound and cytosolic species. In comparison to the cytosol, the cell membrane is modeled as a thin layer to allow for diffusion in lateral direction only [39]. Furthermore, we are not restricted to a circular or elliptical geometry but can choose any shape.
Our computational approach relates to recent advances in surface finite element methods (SFEM) [29,40] which also take into account surface curvature of the membrane. With these methods, the RD model equations can be solved with high accuracy and great flexibility regarding the spatial geometry.

Modeling yeast cell polarity
A prominent example for cell polarization is the yeast budding process. The location of bud formation is specified by cortical landmark proteins Rsr1/Bud1, which are inherited from the previous cell division.
One key molecule during this process is the Rho GTPase Cdc42. In the absence of Cdc42, cells fail to polarize, remain circular and are, therefore, unable to form a bud [41]. As known from all Rho GTPases, Cdc42 cycles between an active GTP-bound and an inactive GDP-bound state. A prenyl group anchors Cdc42 into the membrane where it can diffuse slowly. Cdc42 can be pulled out of the membrane by Rdi1, its guanosine nucleotide dissociation inhibitor (GDI). Active Cdc42-GTP levels can be upregulated by guanine nucleotide exchange factors (GEFs) or reduced by GTPase-activating proteins (GAPs). In the course of these interactions, a unique polarization cap of active Cdc42 is generated on the cell membrane which initiates the emerging bud on the site of the polarization cap [14]. A sketch of this process is shown in figure 1.
Interestingly, the normal bud site selection machinery, in which the bud scar from the previous cell division determines the direction of budding, could be perturbed systematically in a number of experiments. First, it could be shown that overexpression of the active form of Cdc42 leads to phenotypes that polarize randomly at multiple spots and also grow more than one bud [42]. Second, deleting the landmark protein Rsr1 still yields a polarization cap which randomly explores the cell periphery. Additional treatment with latranculin A, an actin-depolymerizing drug, causes cells to spontaneously polarize and establish a random but stable axis of polarity [7]. This indicates that polarization is possible without actin-driven rearrangements of Cdc42.
In S. cerevisiae, an alternative pathway that leads to Cdc42 activation is induced by the heterotrimeric G protein signaling during the mating process [43]. Here, shallow pheromone gradients cause the yeast cell to grow a projection (called shmoo) in the direction of a mating partner. In this case, Far1p associates with the Gβγsubunits which leads to the accumulation of Cdc42p at the site of receptor stimulation. Even when stimulated with conflicting or noisy pheromone signals, cells are still able to form a single shmoo in an arbitrary direction [30].
The cell-type specific details of the spatial molecule cycling can be very complex. For instance, signaling molecules can associate to the plasma membrane through membrane anchoring, binding to phospholipid headgroups or docking with transmembrane receptors. To be able to find generic properties of cell polarization, we aim to investigate abstracted mechanisms that represent generic features of eukaryotic cell polarization. A detailed review of polarization models for various cell types can be found in [16]. We are going to investigate basic mechnisms of Turing-type [8,11] and wave-pining type [7,10], which have both been suggested for yeast.
As the first mechanism of Turing-type, we tested the model from Goryachev and Pokhilko [11] and refer to it as GOR model. This model was derived in a bottom-up approach that takes all actin-independent components of the yeast polarization machinery into account. Analysis and reduction of the model led to a Turing-type mechanism for yeast polarization with just two components. The system starts in a homogeneous steady state which is unstable with respect to minute spatial perturbations and runs into a spatially inhomogeneous polarized state. Such perturbations may be small localized cues or spontaneous association of active molecules to the membrane.
As the second mechanism, we tested the wave-pinning (WP) model which was first suggested in Mori et al [10]. This model is an abstracted caricature of the shuttling of Rho GTPases. It has been used to characterize the polarization behavior of various cell types from neutrophils to Dictyostelium discoideum and also yeast. Experimental evidence for traveling waves of active Cdc42 were observed for instance in Ozbudak et al [7]. In the WP model, the corresponding homogeneous system has two stable and one unstable transient steady state. The steady state at lower concentration of active membrane-bound molecules corresponds to an unexcited state, while the steady state at higher surface concentration corresponds to an excited state. Stimulation of the system causes traveling waves on the cell surface that are pinned into a polar stable state due to global mass conservation. However, the stimulation needs to exceed a certain threshold to drive the system from a homogeneous unpolarized state into a polar state.
In this work, we examined and illustrated the importance of spatial properties such as cell shape, cell size, inhomogeneities on the membrane and the location of (large) organelles in the cytosol for the two models described above. There is still a debate about which kind of model is most appropriate to capture the main features of yeast polarization and which is best suited to develop an intuition of the polarization process in various situations [5]. Both models, WP and GOR, were tested in setups relevant to the polarization behavior of mating and budding yeast Saccharomyces cerevisiae.

Modeling approach for a class of mass conservative reaction-diffusion mechanisms
A traditional way to understand basic mechanisms of cell polarity is the investigation of RD models where the geometry of the cell is reduced to one space dimension. With models of this kind, basic features of the reaction-diffusion dynamics can be captured. The cell is represented by an interval L 0, [ ]which can be either interpreted as the circumference [44] of the cell or as a cell diameter transect [16].
In the first case, the ends of the interval are glued together and periodic boundary conditions are employed. In the second case, the boundary of the interval represents front and back of the cell. In the following model, the concentration of active Rho GTPase (GTP, membrane-bound) is denoted by u(x, t) and the concentration of inactive Rho GTPase (GDP, cytosolic) is denoted by v x t , . ( ) Using the interpretation as a cell diameter transect, the set of RD equations to describe the spatial dynamics reads with some initial values u , 0 ( · ) and v , 0 ( · ) in L 0, [ ] and no-flux boundary conditions A comparison and analysis of different reduced models can be found in [16,45]. A drawback of these 1D models is that they neglect potentially important details of the cell geometry, whereas models in 2D or 3D can also account for the spatial structure of a cell. We thus investigate higher dimensional models and focus on the 2D case in this article.
In the geometric representation of the cell, we incorporate the cell membrane, the cytosol and inner membranes that enclose organelles like the vacuole or the nucleus. We assume free diffusion of molecules in the cytosol which is limited by the outer cell membrane, M cell , and the inner membranes, M org , of the organelles. The cytosolic volume is denoted by V cyt and its boundary surface is V M M cyt cell org È ¶ = (see figure S1).
Here, ∂ D denotes the boundary of some domain D. The active membrane-bound form of the signaling molecule is represented by its concentration u x t , ( )  and the inactive cytosolic form is represented by its where T 0 > is the end time. In brief, in the two-dimensional case, the plasma membrane is represented as a one-dimensional curve enveloping the two-dimensional cytosolic domain.
With this geometric representation of the cell, the shuttling between membrane-bound and cytosolic signaling molecules is naturally described by a flux at the membrane-cytosolic interface which follows specific reaction kinetics depending on concentrations u and v. The model is formulated as the system of partial differential equations For f u v , 0 ( ) > we obtain a flux from the cytosol onto the membrane whereas for f u v , 0 ( ) < the active form u dissociates into the cytosol. We use the term )represents the reaction kinetics of the model's basic reaction mechanism. The function k S is positive and depends on space and time, i.e.
It accounts for signals that excite the system and cause a flux from cytosol to membrane. The signal settings in this paper vary from conflicting localized stimuli, graded signals to external noise (see Methods section).
Since the source term in equation (4) and boundary condition (6) balance each other, this system of equations is mass conservative. This means that the total mass of the considered signaling molecule is con- A derivation can be found in supplementary material S1. Note that in the two-dimensional case, which is considered in this work, the surface integral is an integral along the onedimensional curve enveloping the cell and the volume integral is an integral over a two-dimensional domain, which represents a slice of the cell. Therefore, the unit for M becomes molecules per height. To obtain actual molecule numbers in the cell slice, M needs to be integrated along the height of the cell slice.
In the general framework for cell polarization models, the reaction kinetics f u v , cyc ( )incorporates a self-amplifying feedback on the active membranebound species u. The location where most active signaling molecules are accumulated also recruits the most inactive signaling molecules from the cytosol. This has to be understood as a competition between different polarization sites since the cytosolic pool of inactive molecules is limited. There are several different model realizations for f u v , cyc ( ) with different properties [8,10,11].
We employ the kinetics of the GOR [11] and WP [10] models in the remainder of this study. The GOR model [11,14,16] represents a Turing-type mechanism and was derived as a simplification of a detailed biochemical signaling pathway involved in the yeast budding process. It employs the reaction kinetics The terms α E c u 2 v and β E c u v account for the positive feedback loop of Rho GTPase activation mediated by the corresponding GEFs. The first nonlinear term additionally assumes a cooperative effect of the active form of the Rho GTPases onto its own production. The dissociation of the membrane-bound Rho GTPase mediated by its GAPs and GDIs is described by the degradation term γu. The factors α, β and γ are constants. E c accounts for the Bem1-Cdc24 complex occuring in yeast and it is assumed constant for simplicity, see also [16]. For a cell of 5 μm diameter, Cdc42 molecule numbers in the order of 10 5 are assumed in [11]. However, recent measurements [9,15,46] have reported lower Cdc42 molecule numbers in the order of 10 3 to 10 4 . Therefore, we rescaled the parameter values from [11] (see table 1) .
) is plotted in figure S2 (see Table 1.

Model
Parameter/entity Value Unit Description and u 0 ≡ 0.054 μm μM as initial conditions for the GOR model throughout this paper. This corresponds to the non-zero homogeneous steady state of the system.
The WP model was thoroughly investigated in [10,16,47,48]. Its reaction kinetics reads The nonlinear Hill-Term u K u accounts for cooperative effects of membrane-bound Rho GTPase activation which is mediated by its GEFs. Here, the constants γ and K represent the maximal activation rate and the concentration of half occupation. The basal GEF activation rate is expressed by k 0 . The inactivation and dissociation of the membrane-bound Rho GTPase mediated by its GAPs and GDIs is described by the degradation term δu. As with the GOR model, the reaction kinetics is plotted in figure  S2 for fixed values of v and the kinetic parameters are shown in table 1. In the original work of Mori et al [10], both u and v were treated as effective concentrations with the same units anlong a 1D cell transect. Here, u represents a surface concentration on the membrane and v a bulk concentration in the cytosol. Therefore, we adapted the units of the parameters k 0 , γ and K accordingly. Furthermore, we rescaled the parameter K to be in a comparable regime of concentrations as with the GOR model. In all experiments we used v 0 ≡ 0.2 μM and u 0 ≡ 0.026 μm μM as initial conditions for the WP model, which corresponds to the stable unexcited homogeneous steady state of the system. It should be noted that for fixed corresponds to a stable unexcited, the second to an unstable transient and the third to a stable and excited state of the system. In the original works from Goryachev and Pokhilko [11], and Mori et al [10], different values for the diffusion coefficients were used even though both refer to measured values in the model organism yeast. Each set D GOR * and D WP * of considered diffusion coefficients can be compared using table 2. The membrane diffusion coefficients vary about a factor of 40. These discrepancies have been inherited in recent literature.
In [14], for instance, the membrane diffusion coefficient was employed for the same model organism yeast. However, in this work it is not our aim to resolve reported parameter variations, or different modeling choices, but to illustrate systematically that the behavior of two common RD-based polarization models can be dramatically influenced by cell shape and diffusion inhomogeneities. In the experiments shown in this work, we thus simulated both models with a set D cons * of consensus diffusion coefficients which are approximately the geometric mean of extreme values of diffusion coefficients found in the references (see table 2).
To get a first impression of our higher dimensional formulation of the GOR model and the WP model, it is reasonable to perform 1D simulations of the original models and comparable 2D simulations. We employed a fully circular cell without any inhomogeneities for the 2D simulations. It could be shown that both higher dimensional models are capable of qualitatively reproducing the 1D simulations which are performed in review article [16] with parameters given therein. A comparison of 1D and 2D simulations with two competing stimuli, for example, can be found in figure S3 in supplementary material S1.
In silico experiments (A) Protrusions in a cell locally limit molecule aggregations In a first setup, we strive to understand the influence of the cell shape on the polarization behavior by the introduction of a protrusion to an otherwise fully circular cell. In [19], it was shown that the cell shape potentially influences signaling within the cell for one cytosolic species, neglecting the membrane-cytosol shuttling. In a computational study [25] for models of motile cells like neutrophils or Dictyostelium, the polarization was also examined in the case of a static geometry. For an elliptical geometry, cells repolarized along the major axis when initially polarized along the minor axis without further stimulus or bias. For the MinCDE system, a similar effect was observed in a model with membrane-cytosol coupling and an elliptical cell shape [49]. There, pole-to-pole oscillations of Min proteins predominantly occured along the major  axis as well. However, for a non-oscillatory polarization model as used in our study and a more complex geometry, it was not obvious what effect to expect. In the following experiments, we incorporate a protrusion into a circular cell as shown in figure 2. The cell shape is motivated by the shape of a yeast cell during mating with the protrusion representing a mating projection. However, similar shapes can also be found in other cell types like dendritic spines [50]. First, a homogeneous signal was employed to excite the cell from its resting homogeneous steady state. Using the mechanism introduced in Since no spatial perturbations were present in the circular cell, it did not polarize. However, homogeneous excitation of the non-circular cell led to cluster formation opposite to the protrusion. (B) The GOR and the WP model were simulated for a cell with a protrusion and a circular cell. The signal comprised two competing stimuli S1 and S2. The amplitude of S1 was chosen 10% larger than the amplitude of S2 (see Methods section). Note: The non-circular cell has a length of 7 μm and a width of 5 μm; the circular cell has a diameter of 5.4 μm. For both cell shapes the cytosolic volume is the same, while the total arc length of the circumference is slightly larger for the non-circular cell. equation (9), the signal k s 0.03 m S m º was applied for a time period of Δ t = 100 s on the surface of a circular cell with a diameter of 5 μm and a cell which has a protrusion with length 2 μm (see figure 2(a)). As expected, the circular cell did not polarize since there were no asymmetries in the signal or in the cell shape. Once the signal vanished, the cell returned to its unexcited resting state. For the non-circular cell, we observed a different effect. A flux was induced from the cytosol to the membrane. This created a cytosolic gradient from the protrusion (lower concentration) to the opposite side of the cell (higher concentration) since the ratio surface to volume was larger in the protrusion than in the circular part. The concentration of the surface constituent, therefore, also grew faster on the opposite side of protrusion. Due to the positive feedback, the concentration gradient on the surface was amplified and a unique cluster was established. The positive feedback in both models also led to a dramatically reduced concentration of inactive molecules in the cytosol which reverses the cytosolic gradient (shown in figure 2(a)). The membrane-cytosolic flux was still positive in the cluster region but negative elsewhere. The geometry effect was more pronounced for the WP model where polarization was established very fast (less then 100 s) while the GOR model polarized much slower (∼1000 s).
Second, we imposed a signal comprising two stimuli S1 and S2 on the cell surface. The stronger stimulus S1 was located at the protrusion while the weaker stimulus S2 was located at the opposite side of the cell (see figure 2(b)). For the circular cell, we observed the expected outcome: the stronger stimulus S1 induces a cluster of active membrane-bound molecules which dominated the smaller cluster induced by S2. Eventually, the smaller cluster vanished and the system reached a polarized steady state for both models. However, we observed that the protrusion reversed the outcome. As in the former setup with a homogeneous signal, the cytosolic concentration decreased faster in the protrusion and therefore limited the cluster growth at the site of S1. Eventually, the cell was oriented toward the opposite side of the protrusion.
The effect that we observed in the simulations can be described as a kind of 'bottle neck' caused by the protrusion at the site of S1. In particular, diffusive transport into the protrusion was slightly hindered when compared to diffusion on the unperturbed membrane. Thus, the location with the stronger stimulus S1 at the protrusion was not able to 'win the competition' against the initially smaller cluster at the site of S2 due to insufficient transport of inactive signaling molecules to the site of S1. This effect might be intepreted in the context of the narrow escape problem [50]. The emerging cluster acts as an absorber of cytosolic molecules since there is a flux from the cytosol to the membrane at the emerging cluster. The mean time for a molecule in the cytosol to reach the cluster (or absorber) is greater at the protrusion than in the circular part of the cell surface. Slow cytosolic diffusion or narrowing of the protrusion increases the mean time and, hence, slows down molecule aggregation in the protrusion. For analytical results of the narrow escape time applied to dendritic spines, we refer to [51].
In [30], we stimulated yeast cells with pheromone where yeast cells formed protrusions toward a spatial gradient. When the gradients were almost uniform or the gradients were perturbed by small rapid changes, the yeast cell still polarized and grew a protrusion several times, but not in the same direction. This effect has also been observed in other experimental works for yeast [52][53][54]. Moreover, in experiments with Dictyostelium cells and neutrophils it was reported, that unless gradients were very steep, new pseudopods steered cells away from the direction of the old pseudopod [55]. In a dynamic cell shape model [29], the effect was explained by the dilution of the activator by surface expansion of the cell surface in the pseudopod region. Our study suggests that this effect can alternatively be explained by a narrowing of the volume that occurs when a pseudopod is extended.
Based on these observations, we complemented our simulations with stimuli that act as gradients on the cell surface. Again, the protrusion acted as a negative feedback and even if the protrusion was aligned toward the gradient, the polarization occured on the opposite side of the cell if the gradient was very shallow (see figure S7 in supplementary material S1). This demonstrated nicely that small perturbations in the cell shape may lead to a qualitative difference in the behavior of both models. One interpretation of this observation is that the growth of a protrusion evokes a direct feedback on the underlying biochemical process, i.e. spatial geometrical changes cause changes in the biochemical properties of signaling pathways. This effect might be a generic feature of cell polarization that damps spurious deformations and promotes reorientation in case of ambiguous signals. m m - [11,16,56], we can conclude that the root-mean square distance D t 4 c D varies between 2 μm-6.32 μm in a time interval of Δt = 1s. For a small cell with a diameter of 3 μm, which can be assumed for small yeast cells, this means that one molecule typically travels from one end of the cell to the other in less than one second. It thus is apparent that passive transport processes are much more efficient in small cells than in larger volumes. However, higher concentration differences as well as multidirectional gradients that interfere with each other are expected in larger cells.
There are several mathematical approaches for the analysis of RD models. These can be used to obtain parameter studies of crucial model parameters such as cell size, intial conditions kinetic parameters and diffusion coefficients. A classical method is linear stability analysis (LSA) which was performed in [15,57,58] for bulk-surface RD systems. In this approach, the equations are linearized around a homogeneous steady state and the polarization behavior of the model is studied in the linear regime. There is also a more recent approach called local perturbation analysis [59,60] which characterizes the nonlinear regime and reduces complex reaction-diffusion equations to a small set of ordinary differential equations. However, the analysis is based on the assumption of very fast cytosolic diffusion coefficients (e.g. D c  ¥) and very slow diffusion coefficients on the membrane (e.g. D 0 m  ) and, hence, the interdependence of fast and slow diffusion is not investigated. Therefore, we performed the classical LSA which was complemented with numerical simulations in the nonlinear regime. We point out that there is also a possibility to extend the LSA into the nonlinear regime [45], however, this analysis goes beyond the scope of this study and could be addressed in future studies.
In the following analysis, we consider circular cells of radius R. The surface constituent u was decomposed into a Fourier series with summands k exp cos k ( ) ( ) l f and k exp sin k ( ) ( ) l f, where k is the wave number and λ k the corresponding growth mode. The volume constituent was expanded into a series using the fundamental solutions for the radial part, i.e. modified Bessel functions of the first kind (see supplementary material S1 for details). The fundamental solutions for both the surface and the volume constituent are shown in figure 3(a) for different wave numbers k. Note that the wave numbers k 1  correspond to the number of polarization sites of the cell. We derived an analytical expression for the corresponding growth modes λ k . These are depicted in figure 3(b) in dependence on the cell size. For larger cell sizes, multiple peaks become unstable, while for small cells all λ k are negative and, therefore, no polarization occurs in the linear regime. A phase portrait for the influence of cytosolic and membrane diffusion coefficients is shown in figure 3(d). To get an understanding of the interplay of postive feedback, cell size and diffusion coefficients, the following estimate for λ k was derived: In experiments [42], it was shown that overexpression of Cdc42 led to multiple polarization sites and the dependency on the GEF was reduced. Motivated by these experiments, we varied the total molecule number of signaling molecules with respect to the positive feedback mediated by its GEF. In the case of the GOR model, we varied the crucial parameter E c and for the WP model the parameter γ as shown in figure 3(c). The steady state u v , 0 0 (¯¯) for the stability analysis depended on changing E c and γ for the GOR and WP model, respectively. The number of molecules was deduced from the mass of the steady state as given in equation (10) by assuming a slice of the cell with height 1 μm. For the GOR model, we observed the expected effect that higher molecule numbers led to multiple clusters. Furthermore, polarization occured already for a slow positive feedback E c if the molecule number was high. For the WP model, there was an upper bound for the molecule number that led to polarization. Here, a strong positive feedback led to a wider range of molecule numbers. For further analysis of the WP model, we refer to [47].
Simulations for circular cells with cell diameters in the range of one order of magnitude from 1.5 μm up to 15 μm were performed to demonstrate the behavior of both models in the nonlinear regime. For all cell sizes, we first imposed one single stimulus and second two conflicting stimuli S1 and S2 on the cell membrane. The amplitudes of the stimuli were fixed while the membrane areas of the stimuli were scaled proportionally to the cell size to make simulations with different diameters comparable. In general, two clusters emerged and competed with each other. It is of interest which influence the cell's diameter has on this competition.
In order to evaluate the influence of cell size, we introduced two measures to compare the degree of polarization on the membrane (see Methods section). The first was the mean polarization (POL) which related the highest concentration to the total mean concentration on the surface normalized by the surface area. The second was the polarization factor (PF) which was calculated from the area of the smallest surface patch that comprised half of the mass on the surface. For the steady states of the simulations, we calculated PF and POL for varying cell sizes, see figure 4.
The GOR and WP models exhibited essentially different polarization behaviors. For the GOR model, the clusters grew mainly in 'height' (i.e. locally high concentrations), while for the WP model a traveling wave could be observed where the cluster grew in width but not in height.
For the GOR model, we observed that the POL measure slightly decreased, since the cluster height did not grow proportionally to the cell diameter (compare figure 4). On the other hand, the relative cluster width was decreasing and therefore the PF measure was increasing with cell size. The development of the relative cluster width in relation to the cell size is shown in figure S5 in supplementary material S1. It should be noted that for a cell with a diameter smaller than 4.5 μm, polarization was not achieved for the employed parameters.
For the WP model, the POL measure decreased with increasing cell size since the cluster height was the same for smaller and larger cells. The PF measure decreased as well. The surface to volume ratio scales with 2/R. Hence, if the maximum concentration of the cluster is fixed, the cluster width grows with cell size. An explicit formula for the cluster width is derived in supplementary material S1. The formula was in almost exact agreement with the numerical 2D simulation and predicts that the relative cluster grows linearly with R. For a cell with a diameter of 15 μm, more than 80% of the cell surface was polarized (see figure S6 in supplementary material S1). For a cell with a diameter smaller than 3 μm, polarization was not achieved in the case of one and two stimuli. In the case of two stimuli at opposite sites, no steady state was reached in a time frame of less than t = 2000 s for a cell with a diameter larger than 6 μm. We also compared the duration for the build-up of polarization by measuring the time until 90% of the final PF value was reached. The increase for larger cells was especially significant. In case of one single stimulus, we observed an increase from 21 s for a cell with diameter 4.5 μm to 55 s for a cell with diameter 15 μm for the GOR model. For the WP model, we observed an increase from 17 s (cell diameter 3 μm) to 80 s (cell diameter 15 μm). Much more pronounced was the increase in case of two stimuli. In fact, for the WP model the polarization time increased from 43 s (cell diameter 3 μm) to 953 s (cell diameter 6 μm) and a single site of polarity was not established within 2000 s for cells with a diameter larger than 6 μm. For the GOR model, the polarization time for a cell with diameter 15 μm more than doubled compared to a cell with diameter 6 μm.
For both models, we observed that polarization is either not possible for very small cells or takes very long for larger cells. Regarding the chosen parameter values, we could demonstrate that there is indeed an optimal cell size for both mechanisms.

(C) Membrane barriers can amplify cluster formation
Interior subdomains such as organelles and diffusion barriers on the membrane potentially play an important role in the signaling process. In budding yeast, bud and birth scars occur after cell division and may influence subsequent cell polarization. Furthermore, it is known that during cell division diffusion barriers are established to separate material of mother and daughter cells while they still share a contiguous membrane [61][62][63]. In a recent study [14], it was shown that septin structures in yeast are formed in the early phase of polarization and that Cdc42 clusters can be trapped in these regions. Furthermore, septins also play an important role in steering the direction during chemotropic growth [52]. We introduced membrane diffusion barriers into the model to examine the influence of these inhomogeneities on polarization behavior, see figure 5.
As in the previous setups, a signal comprising two stimuli S1 and S2 was applied on the cell surface, but this time with the same magnitude to compare the effect of cell shape with the effect of slow membrane diffusion (see figure 5(a)). The cell geometry was exactly the same as in section (A). First, we simulated the effect of diffusion barriers surrounding the location of stimulus S1 (see left column in figure 5). Interestingly, the cluster induced by S1 grew steadily while the cluster induced by S2 vanished. Hence, the diffusion barriers compensated for the negative feedback caused by the protrusion as examined in (A). A similar effect was observed when applying a gradient on the cell surface. Again polarization was enhanced for both models. For the WP model the steady state was almost the same but was achieved in a shorter time frame. The effect was much more pronounced for the GOR model. Here, polarization was only achieved with added diffusion barriers.
This effect can be attributed to an accumulation of signaling molecules at the site of S1 since the transport away from S1 is blocked by the introduced barriers.
We repeated the same setup with the same parameters and cell shapes. However, this time we employed only one single diffusion barrier which was placed next to stimulus S1 at the protrusion, see center column in figure 5. In yeast cells, such impermeable regions on the membrane corresponded to septin structures or bud scars which emerged after cell division. It is a topic of current research how septin structures influence the budding process and chemotropic growth. For the GOR model, we observed that due to restricted diffusion caused by the barrier, the cluster at the weaker stimulus S1 grew much stronger than the cluster at S2. We thus had the same effect as for surrounding diffusion barriers. However, the influence of a single diffusion barrier on one side was weaker. In the case of the WP model, it could not compensate for the negative feedback induced by the protrusion. For the GOR model, we observed that the direction of polarity moved toward the diffusion barrier, an effect that was used to explain the occurance of bending shmoos during chemotropic growth of yeast cells [52]. The effects of barriers surrounding the protrusion and impermeable regions in the vicinity of the signal were qualitatively similar. During cluster formation on the membrane, active membrane-bound molecules diffuse laterally from high concentrations to low concentrations with diffusion coefficient D m . The effect of diffusion barriers introduced in a 2D model as shown in this work can be different from those in a 3D model, since in 2D diffusion barriers block diffusion while in 3D there are more possible shapes of barriers. For example in 3D a ring shaped diffusion barrier around the neck of the shmoo is more constraining than a spot-shaped or circular diffusion barrier at the side of the shmoo. However, in both cases 2D and 3D, reducing the diffusion diminishes the diffusive flux which counteracts the growth of the cluster. Therefore, diffusion barriers on the membrane have the potential to accelerate, stabilize and steer cell polarization.
(D) Organelles in the cytosol can alter polarization preferences Cells accommodate many structures of different sizes, for instance large membrane structures like the endoplasmatic reticulum, the nucleus and other organelles. Thanks to recent advances in imaging technologies, very detailed microscopic images of the cytoplasma can be produced, see for example [2,3,64]. Intracellular structures certainly influence diffusive but also vesicular transport in the cell [35,65,66]. During cell division, the spatial position of the organelles has to be organized and is most likely controlled by signaling, but effectively also influences signaling itself [67]. Even for small cells, the calculation of effective diffusion coefficients usually is inferred from molecule properties and viscosity of the medium. However, these derivations do not incorporate intracellular diffusion barriers such as organelles or impermeable membranes in general. Therefore, we examined the effects of obstacles in the medium on intracellular gradients.
To address this goal, we performed computational experiments with organelles of different size and shape. In a first setup, we placed a large organelle with elliptical shape in the vicinity of the membrane. In a second setup, we added a smaller circular obstacle opposite to the large organelle near the membrane. In each case, a vast number of simulations with different noisy signals on the cell membrane was carried out. To assess the influence of the organelles, the cell surface was devided into twelve equal parts. For each simulation, the observed steady state polarization direction was assigned to the respective part in order to obtain a likelihood distribution of preferred polarization sites, see figure 6.
The effect of the large organelle is very clear and polarization close to the organelle turned out to be very unlikely for both models WP and GOR. In 100 simulations for each setup and model, at most 3 polarization events could be observed in one of the segments next to the organelle for the GOR model and 0 for the WP model. However, in a circular cell without diffusion barriers, the expectation value of polarization events for each area segment is 100/12 ≈ 8.3. A similar but much less pronounced effect was induced by the small circular organelle. In particular, an interesting pattern could be observed for the WP model. Most polarizations occured in the neighborhood of (but not behind) the large organelle or at the opposite side. This can be explained by the complex cytosolic gradients that form due to the cytosolic barriers. In case of polarization, the growing cluster acts as an absorber of cytosolic molecules. If aggregation occurs behind the organelle, an initial cluster is quite likely to vanish due to the limited transport. If polarization takes place further away from the organelle, molecule transport is hindered by the organelle and a cytosolic gradient toward the organelle is formed which leads to a movement of the cluster toward the organelle. This effect was much more pronounced for the WP model but it could also be observed for the GOR model. To complement the experiments, a deterministic investigation of the effect for organelles of different sizes is shown in figure S4 (in supplementary material S1).
We conclude that the GOR model is more gradually influenced by the shape of the organelles in the cytosol, while the steady state of the WP model behaves rather like a switch with respect to the introduction of different cytosolic diffusion barriers.

Discussion
Cell polarization is a fundamental process during cell division, cell differentiation and directed growth. In this paper, we focused on the crucial initial phase of polarization when a polarization site is established. Since spatial aspects are often neglected in computational as well as experimental investigations, our main question was whether and how spatial parameters like cell size, cell shape and spatial inhomogeneities influence the polarization process.
Our study is based on a class of conceptual two component models which describe the shuttling of molecules from the cytosol to the membrane and vice versa. Various models of this type can be found in the literature. Since many of these models have been formulated and simulated only in a 1D environment, we suggest a consistent expansion of this class of models to 2D and 3D. In the framework presented in this study, different cell shapes, curved membranes, diffusion inhomogeneities and barriers were integrated. This allowed us to go beyond the variation of kinetic parameters to the introduction of several crucial spatial properties and their influence on the polarization process. For this, it is essential to perform simulations at least in 2D so that differences in surface to volume ratio, cell shape and organelles placed in the cytosol can be represented and taken into account.
In our study, we illustrated that the influence of cellular inhomogeneities can be quite dramatic and should certainly be considered in modeling and simulation of spatial intracellular processes. For the reaction kinetics, we chose the GOR model with a Turingtype mechanism on the one hand, and the WP model with a wave-pinning mechanism on the other hand. With these two kinetics, a vast number of computational experiments were carried out in order to examine the influence of (A) the cell shape, (B) the cell size, (C) inhomogeneities on the membrane and (D) organelles in the cytoplasm.

Importance for cell polarization in yeast
For the first setup (A), we illustrated that protrusions act as a negative feedback for both mechanisms, which led to a lower sensitivity to external signals in this region. In the presence of high pheromone concentration (≈1μM) for a long time period (t > 2h), yeast cells form multiple shmoos in distinct directions [30,54,68,69]. The periodic formation of mating projections has been explained by a downregulation and termination of polarized growth meditated by Sst2 [68,70], a regulator of G protein signaling (RGS). After termination of the growth of the first mating projection, the system polarizes again in the presence of high pheromone concentration. The second mating projection forms in a distinct direction than the first one. It would be interesting to investigate experimentally whether the reorientation is mediated by a time delayed memory e.g. localization of Sst2 in the first mating projection, or by the geometry dependence of the Cdc42 polarization module as shown in this study.
Furthermore, we demonstrated in setup (B) that using fixed kinetic parameters there exists an optimal cell size for the GOR and the WP model. For both mechanisms, the polarization measured by the average polarization (POL) decreases for larger cells, while polarization times increase dramatically. Opposite to that, small cells were unable to polarize at all. We derived estimates for theoretical upper and lower bounds for the cell size depending on the kinetic parameters that control the postive feedback of Rho GTPase activation (see setup (B) or Text S1 for a derivation). It is an interesting question whether the size limits that hold for the theoretical models can be matched with experimental investigations in yeast. An experimental study on size dependent mating partner selection in yeast can be found in [71].
Diffusion inhomogeneities on the membrane are investigated in setup (C). We showed that diffusion barriers can act as a positive feedback on cluster formation. Diffusion on the membrane usually counteracts cluster formation. Hence, local obstruction of membrane diffusion yields an amplification of the signal. Diffusion barriers in yeast are mediated for instance by cytoskeletal scaffolding proteins known as septins. A recent experimental study [52] investigated the interplay of septins and the turning of the polarization cap in life cell imaging upon stimulation with low pheromone concentration. The authors could show that the propability of turning the Cdc42 polarization cap in a short time scale of a few minutes is limited for wild type cells, where the septins form at the edge of the polarization cap. However, deleting the Sst2 resulted in cells, where the polarization cap wandered along the cell periphery and cells failed to grow persistently.
The behavior of the WP model for gradient tracking on short time scales, which is not a feature of the GOR model, might be more suitable to explain the mutant behavior of sst2Δ cells. However, with diffusion barriers both models seem to behave roughly the same, e.g. cluster formation is promoted and stabilized by septins. Both setups (A) and (C) encourage further investigation of the temporal and spatial localization of septins and RGS proteins like Sst2. We point out that diffusion barriers induced by the cytoskeleton can also be modeled as local energy barriers that modulate the obstacle strength dynamically. A model describing the interplay of slow diffusion induced by the cytoskeleton and active transport was investigated in [72], which could a be useful a approach for mating and budding yeast as well.
In experimental setup (D), we introduced different organelles into the cytosol. We illustrated that reduced transport due to obstacles led to a change of the polarization behavior despite fast cytosolic diffusion. In particular, placing an organelle in the vicinity of a cluster limited the membrane-cytosolic shuttling and, therefore, also limited the possible growth of the cluster. We observed that the WP model was more sensitive to obstacles and organelles than the GOR model. We suggest that cell polarity and the organization of organelles should be investigated in two directions; experimentally, by observing the orchestration of molecule aggregations and the postioning of organelles [64]; computationally, by incorporating active transport and the cytoskeleton as suggested in [59,73]. However, the incorporation of the effect of inclusions and organelles in a comprehensive model with active transport and cytoskeleton has to be addressed in future research.
Outlook and relevance for other organisms The results of this paper are based on static geometries which can be assumed for the initial establishment of polarization or for cells that grow slowly in comparison to the diffusion and polarization processes. These assumptions are valid for many fungi, but also plant cells or neurons. Apart from this, the observation that spatial inhomogeneities may significantly influence intracellular transport processes is generically applicable to most cells and requires further research. For instance, we believe that our approach could be beneficially embedded in forward-looking multiphysics frameworks as presented in [74] or [59], where deformable cell geometries are combined with signaling on the cell membrane.
The application of our approach for spatial modeling and simulation in a 3D setting is obvious since it only needs straightforward reformulation. We expect the influence of the cell shape to be even stronger in 3D. As in the 2D results, the local ratio of membrane and cytosol volume will be a key to characterize the influence of the cell shape. Therefore, we expect cluster emergence to preferentially occur in regions of minimal mean curvature in case of a uniform or noisy signal as also hypothesized in [10].
In the case of migratory cells such as neutrophils, fibroblasts or Dictyostelium discoideum, the underlying mechanism differs from yeast. Here, it is well-known that phosphatidylinositol-3,4,5-trisphosphate (PIP3), a membrane lipid that promotes Rho GTPase activation, increases sensitivity at the pseudopods through active transport. However, in experiments with Dictyostelium discoideum and neutrophils it was shown that the leading edge of the pseudopod splits into two and the new pseudopod steers the cell away from the direction of the old pseudopod [55]. In a theoretical study [29] this effect was explained by the dilution of signaling molecules on the cell surface due to surface expansion. Our study motivates two other possibly explanations. First, we clearly expect that narrowing and extending the leading edge increases the surface to volume ratio locally and limits the recruitment of molecules onto the cell surface. Second, we suggest that extending the volume in a dynamic geometry causes a dilution of signaling molecules in the cytosol.
Recent advances in imaging techniques such as TIRF as well as confocal and electron microscopy have the potential to provide quantitative data for the spatio-temporal organization of cellular structures and spatial inhomogeneities. Computer tomographic images of cells can be translated into computational meshes and could serve as a basis for spatial modeling and simulation using real cell data in the future.
Our results suggest many spatial effects that could be investigated with the aid of such methods. Therefore, we hope that our approach leads to experimental investigations that give more insight into fundamental processes such as cell differentiation, directed growth and cell division.

Signal repertoire
To qualitatively compare the influence of different spatial effects, the GOR and the WP models were excited and probed with different kinds of signals which are described in more detail in the following subsection.

1) Conflicting localized stimuli
In this setting, two different regions of the cell were excited simultaneously. The imposed signal comprises two localized stimuli and is given by

3) External noise
The external noise used in setup (D) was assumed to be spatially uncorrelated. Thus, the noise at each location on the cell surface was generated by selecting independent and identically distributed stationary random variables at each time period t t , .
[ ] This realization of the noise was then assigned to k x t , .