A Reaction–Diffusion–Advection Model for the Establishment and Maintenance of Transport-Mediated Polarity and Symmetry Breaking

Cell polarity is a fundamental process in many different cell types. The yeast cell Saccharomyces cerevisiae provides an exemplary model system to study the underlying mechanisms. By combining biological experiments and mathematical simulations, previous studies suggested that the clustering of the most important polarity regulator Cdc42 relies on multiple parallel acting mechanisms, including a transport-driven feedback. Up to now, many models explain symmetry breaking by a Turing-type mechanism which results from different diffusion rates between the plasma membrane and the cytosol. But active transport processes, like vesicle transport, can have significant influence on the polarization. To simulate vesicular-mediated transport, stochastic equations were commonly used. The novelty in this paper is a continuous formulation for modeling active transport, like actin-mediated vesicle transport. Another important novelty is the actin part which is simulated by an inhomogeneous diffusion controlled by a capacity function which in turn depends on the active membrane bound form. The article is based on the PhD thesis of N. Emken, where it is used to model budding yeast using a reaction–diffusion–advection system. Model reduction and nondimensionalization make it possible to study this model in terms of distinct cell types. Similar to the approach of Rätz and Röger, we present a linear stability analysis and derive conditions for a transport-mediated instability. We complement our theoretical analysis by numerical simulations that confirm our findings. Using a locally mass conservative control volume finite element method, we present simulations in 2D and 3D, and compare the results to previous ones from the literature.


INTRODUCTION
The development and maintenance of cell polarity is essential for many biological processes like cell growth, cell morphogenesis, cell migration, cell differentiation, proliferation, and signal transmission. Also known as symmetry breaking, it describes the process by which cells generate an internal, functional, structural, and molecular axis. This asymmetric arrangement often arises due to intrinsic or extrinsic cues which are amplified by transport processes or pathways of diffusing and interacting molecules. The budding yeast (Saccharomyces cerevisiae) is an exemplary model system to study the underlying mechanisms of cell polarization. Whereas in these cells, the small family GTPase Cdc42 is a key regulator of cell polarity, GTPases in general are exemplary for a complex system with symmetry breaking in many eukaryotic cells [8,13,28].
GTPase molecules are able to change between three forms: an active (GTP-bound) membrane-bound state, an inactive (GDPbound) membrane-bound state, and an inactive (GDI-bound) cytosolic state. The regulation of this cycle is controlled by certain exchange factors, GEFs (GTPase-activating proteins), GAPs (guanine nucleotide exchange factors), and GDIs (GTPase dissociation inhibitors), leading to shuttling between the cytosol and the plasma membrane [8,16,29]. Thus, the GTPase cycle is characterized by a coupled bulk (cytosol) and surface (plasma membrane) reaction-diffusion system.
Since coupled bulk-surface reaction-diffusion systems naturally arise in many biological processes, a huge number of studies concerning such systems, like, for example, Refs. 18, 21, and 23, has recently been published. All these models were always based on reaction-diffusion equations posed on the bulk and surface coupled by Robin-type boundary conditions which generate symmetry breaking by Turing-type instabilities. But many cells exhibit a transport machinery characterized by actin filaments or microtubules (see, e.g., Refs. 9,19,24,and 31), which further influence spatial patterns. For example, the budding yeast generates polarity by the coupling of reaction-diffusion to transport systems. Actin cables which are aligned along the plasma membrane transport vesicles containing key proteins required for cell polarization from the interior of the cell to the polarized site (exocytosis) [32]. Simultaneously, molecules are internalized from the plasma membrane to the interior of the cell (endocytosis). To simulate vesicular mediated transport, stochastic equations were commonly used, see e.g. [14,17]. It is observed that transport-mediated recycling of molecules plays a key role in polarity establishment and maintenance as well [33]. For that reason, here, we consider a coupled bulk-surface reaction-diffusion-advection system to investigate the contribution of transport to cell polarity. Following the approach proposed in Ref. 23, we perform a linear stability analysis and derive conditions for a transport-mediated instability, which are confirmed numerically.

Contribution of This Paper
Our main contribution is a continuous model for vesicle transport based on active transport, together with an analysis of its contribution to symmetry breaking. Previous models of actin-mediated polarization were solely based on stochastic simulations [9]. Our continuous PDE model allows for a better characterization of the conditions for polarization.
In Section 2, we introduce the model in its nondimensional form. For details on the model reduction and nondimensionalization, we refer to the Supplementary Material. In Section 3, we analyze in detail under which conditions the model can induce pattern formation and complement these results by numerical experiments in Section 4. The stability results confirm that actin-mediated Cdc42 recruitment can increase the robustness of the system. And, we show the ability of the system to polarize via two independent pathways, as it was observed in experiments for budding yeast cells [27,33].
We further investigate numerically the interplay of active transport and geometrical features like organelles, where our experiments indicate that the presence an actin-mediated pathway accelerates polarization and can even induce different patterns.
We conclude with a discussion on the biological implications of our numerical findings.

MODEL DESCRIPTION
We consider a generic reaction-diffusion-transport system that is based on a complex model for cell polarization proposed in [6]. This model was motivated by the influence of vesicle transport along actin cables on the cell polarization, as described in Refs. 9 and 33 (see Supplementary Material for model reduction and nondimensionalization). This system differentiates among one active membrane-bound, one inactive membrane-bound, and a cytosolic state. This model includes the distribution of actin cables on the membrane as an additional component. Its dynamics are described by an inhomogeneous diffusion proportional to the membranebound component modeling the described transportmediated feedback loop.
In the following, we consider a stationary bulk domain Ω and its compact hypersurface Γ : zΩ. We denote by n → the outer normal on the smooth, closed surface Γ, by ∇ Γ the tangential gradient on Γ and the Laplace-Beltrami operator Δ Γ . Let u, v : Γ × I → R be smooth functions denoting the chemical concentrations or species that react and diffuse on Γ in a fixed time interval I : [0, T] ⊂ R. For substances that diffuse or move by advection in the volume Ω ⊂ R n , we consider smooth functions U, V : Ω × I → R. To proceed, we denote by w : Γ × I → R, a smooth function representing a transport control factor, in our case, the density of actin cable ends on the surface Γ. Furthermore, c(u) > 0 describes a capacity function controlling w and hence impacts actin cable assembly. This model follows the observation that actin is essential for cell polarization [1], and vesicle transport (in yeast cells) happens along actin cables. We further include the fact that actin has a reduced dissociation rate, where Cdc42 concentration is high [33], which leads to a u-dependent actin density in the membrane. We model this by an inhomogeneous diffusion and the nonlinear capacity function c(u). Where c is large, the likelihood for actin to bind is high, whereas a small c increases the probability that the actin cable moves away. Further details can be found in the Supplementary Material.
This leads to the following nondimensional coupled reactiondiffusion-advection system with coupling conditions and initial conditions at time t 0 Here, the nonlinear functions f and g, respectively, represent activation and inactivation of the species, h describes adsorption and desorption of molecules, and v → is the divergence-free bulk velocity field. The parameters D U and D V denote the nondimensional bulk diffusion coefficients and d v , d w > 0 the surface diffusion coefficients, which are assumed to be constant. The nondimensional parameter c > 0 relates to the spatial scale of the cell. REMARK 2.1 (mass conservation). Note that this formulation implies conservation of mass. This means that with dσ denoting the integration with respect to the surface area measure and M the total mass, the system satisfies the condition REMARK 2.2 (velocity field). As vesicles only release their content when being integrated into the membrane, the velocity field v → is conservative, that is, divergence free.
The outflow rate on the membrane depends (potentially nonlinear) on the concentration w of actin cable ends on the membrane. Given an outflow function j(w), we construct a divergence-free velocity field v → ∇ϕ as the gradient of a scalar function ϕ. This potential flow from the internal to the external membrane is caused by constructing a divergence-free velocity field and is computed by solving.
where α describes a potential flow control rate, which limits the transport capacity.

LINEAR STABILITY ANALYSIS
Here, we present a stability analysis of the generic system which mainly follows the analysis shown in Ref. 23, to determine conditions required for pattern formation. We restrict ourselves to the spherical case, that is, Ω : B 1 (0), Γ : zB S 2 .
We assume that the internal pool is sufficiently large and that D u ≫ 1. This ensures a well-mixed internal pool, similar to the assumptions in Ref. 23, so that the feedback loop between the component u (active form) and w (actin) is dominating. This simplification has clear limits; in particular, it cannot handle any transport-limited cases. The validity of this assumption will be backed by the numerical simulations in Section 4.2, comparing the full model with the simplified coupling used in the following analysis.
As the rate of transport indirectly depends on the amount of w on Γ, the actin concentration on the membrane is now only governed by an inhomogeneous diffusion controlled by the capacity function c(u), and we can substitutew : w · c(u) − 1 in Eq. 1c.
The system Eq. 1 then reads with Robin-type coupling conditions.
and the initial conditions (2). In the following, we will denote by x : (u, v, w, U, V) T the vector of concentrations and by x * : (u * , v * , w * , U * , V * ) ∈ R 5 + the spatially homogeneous steady state, such that Following the approach of Ref. 23, we analyze the stability of system Eq. 1 at its stationary states. Focusing on the GTPase cycle, we can interpret f as an activation rate and g as the flux describing membrane attachment and detachment of the GTPase. The function h describes the flux induced by exocytosis and endocytosis. This interpretation corresponds to the following conditions on f , g, and h z v f ≥ 0, z v g ≤ 0, z v g ≤ z u g, and z U h ≥ 0.
For brevity, we introduce the notation assuming that at s * (u * , v * , w * , U * , V * ), the functions satisfy the strict inequalities As in Ref. 23, to determine stability conditions for the system Eq. 3, we use an expansion in spherical harmonics: with scalar functions ψ lm , χ lm : [0, 1] → R, and the orthonormal basis {φ lm } l ∈ N 0 ,m ∈ Z,|m| ≤ l of L 2 (Γ). Then, the Laplace operator can be represented as As a result, the L 2 (Γ) scalar product with φ lm leads to the linearized system where the last two equations correspond to the coupling conditions. We use the following ansatz: which also guarantees that either U lm , V lm ≠ 0 in the whole domain or are identical to zero. We first consider the case U lm , V lm ≠ 0. Then, using U ′ lm λ lm U lm and V ′ lm μ lm V lm , we obtain from Eqs. 5d and 5e In the case λ lm , μ lm 0, it is easy to recalculate that we have with α lm , β lm ∈ R. By contrast, for λ lm , μ lm > 0, Eqs. 6a and 6b are modified versions of Bessel differential equations whose solutions are defined by Bessel functions of first kind. Hence, using the respective modified Bessel function J l+ 1 2 , we derive with ξ l π 2r J l+ 1 2 (r). With this, we finally deduce the ODE system coupled to two algebraic equations given by and the jacobian matrix is the system is given as the stability analysis reduces to an analysis of the eigenvalues of the matrix J F . To determine stability conditions, we compute eigenvalues λ of J F via the characteristic polynomial The eigenvalues are now given by the zeros of polynomial (8). Hence, from Eqs. 7a-7e, as long as U lm , V lm ≠ 0, we acquire that an eigenvalue λ with Re(λ) > 0 exists if and only if first λ λ lm μ lm ∈ R + 0 and additionally with (9b) λ lm fulfills the condition (3) is stable against spatially homogeneous perturbations in the variables u, v, and w if the following condition is satisfied: in which case holds. If either U or V 0, we distinguish two cases and conclude that in Case 1 (U 0), holds. PROOF. We first consider the case l 0. Furthermore, we assume that U 00 , V 00 ≠ 0. Note that in this case, w is always constant and w w 0 . This also implies h w 0. Then, the characteristic polynomial Eq. 8 reduces to For the system to be asymptotically stable in (u * , v * , w * , U * , V * ), it is necessary that all eigenvalues are negative. This means that P 0 (λ) has no zeros in [0, ∞). We rewrite .
For λ 0, it holds P 0 (0) 0. Since w is in this case simply a constant and w w 0 , the linearized system reduces to where u, v, U, and V are constants. Summation of Eqs. 13a and 13b yield With the stationary equations for U and V, we obtain Thus, we get and hence, since u, U, and h U > 0, it holds that h u < 0. Furthermore, together with Eq. 14b, we haver Frontiers in Applied Mathematics and Statistics | www.frontiersin.org November 2020 | Volume 6 | Article 570036 By substituting these relations into Eq. 14 and straightforward calculations, as the first condition, we obtain that this system has a nontrivial solution if With Eq. 10 and the relation g v ≤ g u , we further deduce that Let us now consider the case λ ∈ (0, ∞). From Ref. 23, we know that Since we suppose g V , h U > 0, together with Eq. 16, we obtain that In other words, (15), (17), andw * w* c(u*) imply that for λ > 0, if the conditions from Proposition 3.1 are satisfied, the characteristic polynomial has no change of sign. This inequality is necessary for the stability of the homogeneous steady state.
To investigate if this term is also sufficient to exclude an eigenvalue λ with Reλ > 0, we recheck , and h u < 0. Thus, we have to distinguish two cases. First, consider Sinceκ is decreasing andκ ≤ 1 3 , on [0, ∞), we have the downward estimatioñ For the case We directly conclude thatP 0 (λ) > 0 and prove the assertion for the full system. In order to investigate the system for stability conditions in the absence of some species, we proceed with the special cases V lm , U lm 0.
Case 1 (U lm 0): The system is overdetermined, and it holds that μ lm 0. Furthermore, we obtain Moreover, for U lm 0, the characteristic polynomial reduces to The stability conditions for this case have already been discussed, and the proof can be found in Ref. 23.
Case 1 (V lm 0): The system is again overdetermined, and the matrix has an eigenvalue λ lm 0. Moreover, it holds that This implies that any eigenvalue λ corresponding to the linearized system is given by For l 0, we require that all eigenvalues have negative real parts.
We claim that Furthermore, the characteristic polynomial reduces to Since h U > 0, f v > 0, and h u < 0, we deduce that lim λ → ∞ H 0 (λ) +∞. Using Eqs. 9a and 7e, we further calculatẽ To prevent thatH 0 has only negative eigenvalues, meaning it does not change its sign for λ ∈ [0, ∞), consider This is fulfilled even if Summarized, the derived conditions ensure thatH 0 has only negative eigenvalues.
For P l (λ) in order to change its sign, we finally examine lim λ → 0 P l (λ) and get the condition (18), which is sufficient to ensure a positive zero of P l (λ).
Then, (19) and (20) follow directly with the same argumentation. COROLLARY 3.4. Assume that the system (3) satisfies the condition (10) and either D u or D v are chosen sufficiently large. Then, the instability condition (18) is satisfied if the following conditions hold: Case 1: and for r ± 1 2dv (C 2 ± Q √ ) exists an l ∈ N such that r − < l(l+1) c < r + . Case 2: and with r + as defined above exists an l ∈ N with l(l+1) c < r + .
Frontiers in Applied Mathematics and Statistics | www.frontiersin.org November 2020 | Volume 6 | Article 570036 REMARK 3.5. If U lm 0 and the system fulfills condition (12), then the instability condition (20) holds for sufficiently large D v , if the following conditions are satisfied: Case 1: Case 2: and with r + as defined above exists an l ∈ N with REMARK 3.6. If V lm 0 and the system fulfills condition (12), then the instability condition (20) holds for sufficiently large D u if the following condition is satisfied.
Case 1: Case 2: and with r + as defined above exists an l ∈ N with l(l + 1) c < r + .
PROOF. We first restrict ourselves to U lm , V lm ≠ 0 and the case D u ≫ 1 as well as D v ≫ 1. In order to achieve instability, we consider (18) and narrow down to the coefficient of D u · D v which is given by We defineε , whose roots are given by In order to satisfy condition (18) and to obtain an instability, we now require e < 0. First, assume C 1 ≥ 0 and C 2 > 0. Then, e represents a right displaced upward open parabola which intersects the positive axis at points r ± . Thus, with l ∈ N to ensure e < 0, we have to satisfy the conditions of case 1. By contrast, if C 1 is negative, then the parabola is shifted to the left, and we directly prove case 2 to obtain e < 0. We further consider D u ≫ 1 as well as the case D u ≈ 1. Since we suppose that D v ≫ 1, as before, we observe that either D u le or D v le becomes dominant in (18). This implies that an instability exists for sufficiently large D u or D v .
Finally, with the same argumentation as before, the analysis of the coefficient D u in (19) as well as D v in (19) deduces Remarks 3.5 and 3.6 (for the case U 0, see also Ref. 23).
In contrast to the model of Ref. 23, the conditions in our model depend on w * as well as the capacity function c(u * ) at steady state. This is a direct consequence of the actin part which is simulated by an inhomogeneous diffusion controlled by a capacity function which in turn depends on the active membrane bound form. As a consequence, we have shown that the actin feedback can directly contribute to system instability. This actin-mediated feedback was reported in Ref. 33, and it was suggested that it increases robustness of the polarization and even can ensure polarization in the absence of GDI.

NUMERICAL SIMULATIONS
We follow the methods of line approach to handle time derivatives independent from spatial derivatives.
Throughout this work, we employ a control volume finite element (CVFE) method using first-order trial functions and constant test functions on the dual mesh to discretize in space.
These methods are also known as vertex-centered finite volume scheme and can be formulated as a Petrov-Galerkin method. Advective terms are stabilized using upwinding. In particular, the CVFE method has the property to be locally mass conservative and thus our discrete model recovers this feature of the continuous model. For details of the methods, we refer to text books, for example, Ref. 15.
The temporal evolution is discretized using a simple firstorder implicit Euler method. We solve the arising fully coupled nonlinear system using a Newton-Krylov solver with an AMG preconditioner.
The implementation is based on the Distributed and Unified Numerics Environment (DUNE) framework [2,3] and the dunepdelab package [4]. Coupled bulk-surface problems are solved by the DUNE modules multidomain and multidomaingrid [20].

Reaction Kinetics and Parameters
The generic formulation described in Section 2 allows us to investigate cell polarization under consideration of distinct protein kinetics.
For a particular choice of the kinetics f and g, we simulate the application to different geometries. It serves as an exemplary model to study transport-mediated polarity in different cell types (see Supplementary Material for the derivation of the reaction kinetics).
The functions are given by The particular choice of parameters for the numerical simulation of the nondimensionalized system (1) is given in

Actin-Mediated Cell Polarization
In the following, we confirm the results of the linear stability analysis performed in the previous section. In particular, we compare simulations of the full system (1) and the simplified system (3), which was used in the analysis. As we assumed a wellmixed pool, the effect of exocytosis and endocytosis were assumed to dominate over the actual vesicle transport along actin cables. To simulate transport via exocytosis and endocytosis, we define In the simplified (well-mixed) case, the transport to the membrane is slower, due to the nearly homogeneous distribution of U. Thus, we had to increase e 1 and decrease e 2 to obtain similar results as for the full system, where molecules are actively transported. These rates are chosen such that we obtain similar ratios between internal and membrane components as before. We set e 1 84.3 and e 2 4.167.  In all computations, we use functions f , g as given in (22) and (23), respectively. We use initial concentrations and parameters as given in Table 1.
We numerically solve system (3) for the different cases to investigate its behavior.
The most interesting outcome of the stability analysis is the fact that the conditions determining instability are completely independent of the diffusion parameter d w . This implies that the only requirement on d w is that it must be nonzero. In this case, the capacity function c(u) determines whether the system is stable against small perturbations or not. We further call this capacitydriven instability. Figure 1 shows the development of u in time for distinct values of d w . We observe that even for large changes of d w , provided that d w ≠ 0, the system is always unstable and tends to form a polarized patch. It becomes clear that the capacity function c(u) as well as w * determines the stability behavior. The constant d w only changes the temporal dynamic of polarization (see Figure 1). For reduced rates, the maximum value of u is reached much later. It can be shown that even for d w ≪ 1, the system is still able to form a polarized patch, albeit after a very long time (t > 30).
As mentioned in Remark 2.1, the model is mass conservative. Our numerical model adequately reproduce this behavior due to the use of a locally mass conservative method. The evolution of mass of the different components u, v, U, V as well as the total mass is visualized in Figure 2.
Another result of the stability analysis is the fact that we may observe polarization, even if V 0 or U 0. From Figure 3, we see that the generic system is able to represent these cases. Even in the absence of a cytosolic exchange or a transport mechanism, the system becomes unstable and forms a polarized cluster.
The requirement D u ≫ 1 yields that D v must not be very large to ensure instability. We have seen that even in the case D v ≈ 1, the instability conditions may be satisfied. Our numerical simulations confirm these results. Figure 4 illustrates capacitydriven polarization for the system (3), where D v 1.

Cell Shape Influences Transport-Driven Cell Polarization
Active transport of molecules plays a significant role in many cell types. For example, in the fission yeast, neurons and the Caenorhabditis elegans zygote microtubules may mediate the transport of important regulators of cell polarization and in this way ensure its correct location [19,26,30]. Therefore, our modeling approach can be used to investigate polarization for a range of different cell types with distinct shapes.
In order to understand the influence of the cell shape on polarization, we simulate the system for different threedimensional model geometries. We employ a random signal to drive the cell out of its uniform state. The results are shown in Figure 5. In all cases, we obtain an enhanced peak of the nondimensional concentration u.  One observes that transport-mediated polarization is significantly accelerated in nonspherical cells. In this case, the gradient increases or decreases with the length or broadness of the shape, respectively. Regarding the polarity direction, our results show that transport can change the spatial location of the polarized patch. This becomes particularly obvious in Figure 5D which shows polarity in a cell that features a small bud. In this case, we excite the cell from its homogeneous state by a signal comprising of two stimuli S 1 and S 2 of the same intensity. The signals are imposed on opposite sides of the cell surface, one  November 2020 | Volume 6 | Article 570036 located at the protrusion. Depending on the presence of transport, different patterns are obtained. In the presence of active transport, a peak forms at the bud, without u clusters at the opposite side. The influence of protrusions on diffusiondriven polarization in a cell has already been studied in Ref. 10.
Their results have shown that protrusions locally limit molecule aggregations. Diffusive transport into the protrusion is slightly hindered so that the cytosolic concentration decreases faster in this region. As a result, the cluster emerges at another location. Interestingly, our results demonstrate that for sufficiently high rates of active transport, this kind of "bottle neck" can be compensated, and the cluster forms at the protrusion.
Depending on the particular rates, feedback strength, and the interplay between transport and reaction kinetics, transport can either enhance or disturb polarity. For some choices, it even perturbs the system so strongly that it is no longer capable of polarization.

Influence of Internal Components on Cell Polarization
Cells contain many different cell components of distinct shape and size like for instance the nucleus, the Golgi, or the endoplasmic reticulum. All these structures serve as a kind of diffusion and transport barrier within the cell. In this way, the spatial position of organelles can influence signaling pathways, including the accumulation of polarization molecules.
How internal barriers control diffusion-driven cell polarization has already been investigated in Ref. 10. The results have demonstrated that the cluster formation close to organelles is very unlikely. Diffusion-driven polarization mostly occurred in the neighborhood of large organelles, but not behind them. The local accumulation of substances at the opposite side of protrusions or in regions with low curvature is more likely [10]. In order to investigate whether active transport alters the results, we perform similar computational experiments. We consider the two-dimensional case, where the cell is characterized by a circle. Organelles are modeled by elliptic or circular shapes placed in the cell interior. The results are shown in Figure 6. Again, we excite the cell from its homogeneous state by a signal comprising two stimuli S 1 and S 2 of the same intensity. Whereas one signal is located near the organelles, the other is placed at the opposite side.
Without consideration of transport effects, we obtain similar results as presented in Ref. 10. The organelles near the surface negatively affect cluster formation at this site. Contrarily, we see that under consideration of active molecule transport, the polar cluster forms behind the internal component. In this case, organelles support a nearby spatial location of the polarity patch.
As mentioned before, protrusions positively influence transport-mediated polarization too. This raises the question of how polarity behaves in cells exhibiting both a complex shape and internal barriers. Figure 7 illustrates this interplay. It becomes clear that since protrusions as well as diffusion barriers can promote polarization, the localization of organelles next to protrusions strongly enhances polarity. Conversely, we see that an opposed position leads to a competing situation. As long as the organelle is sufficiently far away from the surface and centrally located, the cluster still forms at the bud. In contrast, when the organelle is placed near the membrane, but opposed to the protrusion, we obtain polarization behind the organelle. Only a very strong stimulus at the protrusion reverses the outcome. This is demonstrated by the last computational experiment illustrated in Figure 7, where the cell is excited at the bud tip with a signal S 1 of strength s 1 0.33.

DISCUSSION
Based on a complex bulk-surface reaction-diffusion-advection system for cell polarization proposed in Ref. 6, in this work, we have introduced a generic approach for the simulation of transport-mediated cell polarization. We performed numerical simulations with distinct cell geometries and cell types, and compared the results to those found in the literature. Since our main interest was to analyze the conditions leading to cluster formation, we further performed a linear stability analysis considering a spherical cell.
The results have shown that vesicular transport may not only influence the robustness, shape, and intensity of the polar cluster but also its spatial location. Particularly, in cells with complex shapes, we observed different patterns between simulations with and without active molecule transport. Here, protrusions and narrower domains differently affected symmetry breaking. Whereas complex shapes rather inhibit diffusion-driven symmetry breaking, transport-mediated polarization can be enhanced under these circumstances.
However, cells are able to robustly polarize at sites of complex protrusions. For example, the tip of the future axon is strongly polarized during neuronal development. These findings suggest that, especially in nonspherical cells, active transport may be required to ensure the correct location of the polarized patch, which is in line with previous finding in Ref. 7.
Based on a complex bulk-surface system for the simulation of cell polarization, we have presented a reduced generic system of bulk-surface reaction-diffusion-advection equations. Our main interest here was to analyze the conditions leading to pattern formation. Therefore, using a spherical cell, we applied a linear stability analysis to a simplified system composing three surface quantities and two bulk concentrations. Our results have demonstrated that two different main mechanisms lead to symmetry breaking. The first one is related to a classical diffusion-driven instability studied in Refs. 22 and 23. The second mechanism is controlled by a capacity-dependent inhomogeneous diffusion of the transport triggering factor. Such dependence has the capability to induce a positive feedback leading to spatial patterns.
However, we have restricted our analytical and numerical studies to stationary domains. In many cases, biological processes induce the development of cell shapes. Thus, the consideration of surfaces which evolve continuously in time would be of great interest. But this implies a more complicated modeling, analysis, and simulation of the coupled system and could be focus of further studies. The results have shown that vesicular transport may not only influence the robustness, shape, and intensity of the polar cluster but also its spatial location. Particularly, in cells with complex shapes, we observed different patterns between simulations with and without active molecule transport. Here, protrusions and narrower domains differently affected symmetry breaking. Whereas complex shapes rather inhibit diffusiondriven symmetry breaking, transport-mediated polarization can be enhanced under these circumstances.
Another outcome of the computational results is the distinct role of organelles. Whereas internal barriers inhibit diffusiondriven polarization behind them, active transport is able to overcome this negative feedback to facilitate polarity next to organelles. The influence of internal components on the direction of cluster formation has already been shown by biological experiments. To give an example, studies with the fission yeast have demonstrated that the position of the interphase nucleus dictates the future site of cell division [5]. These findings together with our results emphasize that it is of particular importance to consider spatial aspects in the mathematical study of cell polarization. As a consequence, to investigate such biological processes in greater detail, the application of more complex mathematical models, including coupling bulk-surface PDEs, must take on greater significance.
Unfortunately, with growing complexity, the analysis of mathematical models becomes increasingly challenging. To enable a linear stability analysis, we continued with a reduction of the generic approach given by reaction-diffusion-advection equations to a minimal coupled bulk-surface reaction-diffusiontransport system. The stability analysis has shown that the reduced generic system is able to generate spatial patterns under certain conditions. These conditions confirm that the transport process derived in this work can increase the robustness of the system. The reason is that two distinct mechanisms act in parallel to generate symmetry breaking. These can explain polarization in Δrdi1 and LatA-treated cells. Treating wild-type yeast cells with latrunculin A (LatA) removes the actin-dependent recycling pathway, while Δrdi1 denotes cells with removed GDI, which both can establish polarization [27].
The first one relates to a classical Turing instability which requires a large difference in the cytosolic and membrane diffusion coefficient. Even if there is no transport of molecules from and to an internal compartment, this mechanism is able to achieve polarization. Since this case has already been analyzed in detail, at this point, we refer the reader to Ref. 23.
The second mechanism is based on a capacity function that regulates the concentration of the component driving transport. Under certain conditions, this mechanism can induce symmetry breaking, even if the cytosolic exchange is blocked. Hence, this case explains symmetry breaking in cells lacking the cytosolic component. In this case, d w ≠ 0, the capacity function c(u) together with the homogeneous state of w entirely determines the stability behavior.
By the performance of numerical simulations, we finally confirmed the results of the stability analysis and demonstrated that our model is able to show the different cases derived. Furthermore, we have shown that this capacitydriven instability also generates pattern when the cytosolic and membrane diffusion rates are equal. For that reason, and since the diffusion constant d w has no essential impact on the stability of the system, we assert that this instability mechanism distinguishes from the Turing-type instability.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
NE contributed to the modeling, analysis, numerical simulations, and preparation of the manuscript. CE contributed to the modeling and numerical simulations, guided the research, and contributed to the preparation of the manuscript.

FUNDING
The content of this manuscript has been published as part of the thesis of [6]