Abstract
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.
1 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 (GDP-bound) 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.
2 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 membrane-bound component modeling the described transport-mediated feedback loop.
In the following, we consider a stationary bulk domain and its compact hypersurface . We denote by the outer normal on the smooth, closed surface , by the tangential gradient on and the Laplace–Beltrami operator .
Let be smooth functions denoting the chemical concentrations or species that react and diffuse on in a fixed time interval . For substances that diffuse or move by advection in the volume , we consider smooth functions . To proceed, we denote by , a smooth function representing a transport control factor, in our case, the density of actin cable ends on the surface . Furthermore, 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 . 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 reaction–diffusion–advection systemwith coupling conditionsand initial conditions at time Here, the nonlinear functions f and g, respectively, represent activation and inactivation of the species, h describes adsorption and desorption of molecules, and is the divergence-free bulk velocity field. The parameters and denote the nondimensional bulk diffusion coefficients and the surface diffusion coefficients, which are assumed to be constant. The nondimensional parameter 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 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 , we construct a divergence-free velocity field 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.
3 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, , .
We assume that the internal pool is sufficiently large and that . 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 , and we can substitute in Eq. 1c.
The system Eq. 1 then readswith Robin-type coupling conditions.and the initial conditions (2). In the following, we will denote by the vector of concentrations and by the spatially homogeneous steady state, such thatFollowing 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 and hFor brevity, we introduce the notationassuming that at , the functions satisfy the strict inequalitiesAs in Ref. 23, to determine stability conditions for the system Eq. 3, we use an expansion in spherical harmonics:with scalar functions , and the orthonormal basis of . Then, the Laplace operator can be represented asAs a result, the scalar product with leads to the linearized systemwhere the last two equations correspond to the coupling conditions. We use the following ansatz:which also guarantees that either in the whole domain or are identical to zero.
We first consider the case . Then, using and , we obtain from Eqs. 5d and 5eIn the case , it is easy to recalculate that we havewith . By contrast, for , 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 , we derivewith . With this, we finally deduce the ODE systemcoupled to two algebraic equations given byWe introduce the notationand the jacobian matrix is the system is given aswithWritingthe stability analysis reduces to an analysis of the eigenvalues of the matrix . To determine stability conditions, we compute eigenvalues λ of via the characteristic polynomialwithThe eigenvalues are now given by the zeros of polynomial (8). Hence, from Eqs. 7a–7e, as long as , we acquire that an eigenvalue λ with exists if and only if first and additionally with fulfills the conditionProposition 3.1. In , the system (3) is stable against spatially homogeneous perturbations in the variables , and w if the following condition is satisfied:in which caseholds.
If either U or , we distinguish two cases and conclude that in
Case 1 (),Case 2 ()holds.
Proof. We first consider the case . Furthermore, we assume that . Note that in this case, w is always constant and . This also implies . Then, the characteristic polynomial Eq. 8 reduces to
For the system to be asymptotically stable in , it is necessary that all eigenvalues are negative. This means that has no zeros in . We rewriteFor , equation is equivalent to
For , it holds . Since w is in this case simply a constant and , the linearized system reduces towhere are constants. Summation of Eqs. 13a and 13b yieldWith the stationary equations for U and V, we obtainThus, we getand hence, since , and , it holds that . Furthermore, together with Eq. 14b, we haverso thatBy substituting these relations into Eq. 14 and straightforward calculations, as the first condition, we obtain that this system has a nontrivial solution ifWith Eq. 10 and the relation , we further deduce thatTogether with Eq. 4, this yields .
Let us now consider the case . From Ref. 23, we know thatSince we suppose , together with Eq. 16, we obtain that . Furthermore, it holds thatIn other words, (15), (17), and imply that for , 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 , we recheckfor , , , and . Thus, we have to distinguish two cases. First, considerSince is decreasing and , on , we have the downward estimationFor the caseWe directly conclude that 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 .
Case 1 (): The system is overdetermined, and it holds that . Furthermore, we obtainand henceMoreover, for , the characteristic polynomial reduces toThe stability conditions for this case have already been discussed, and the proof can be found in Ref. 23.
Case 1 (): The system is again overdetermined, and the matrix has an eigenvalue . Moreover, it holds thatso that for (7b), we obtainThis implies that any eigenvalue λ corresponding to the linearized system is given byFor , we require that all eigenvalues have negative real parts. We claim thatFurthermore, the characteristic polynomial reduces toSince , and , we deduce that . Using Eqs. 9a and 7e, we further calculateTo prevent that has only negative eigenvalues, meaning it does not change its sign for , considerThis is fulfilled even if
Summarized, the derived conditions ensure that has only negative eigenvalues.
We next determine under which conditions small spatial perturbations from the homogeneous steady state induce instabilities.
Theorem 3.2. Assume that the system (3) satisfies the condition (10). If in further holds thatthen, the system is linearly asymptotically unstable in .
Remark 3.3. In the case of , the system is asymptotically unstable in if it fulfills the condition (11), and furthermoreIf and (12) is satisfied, then the instability condition is given by
Proof. Again, we first consider the case . Since we claim that the system becomes unstable in the presence of diffusion, we require that the characteristic polynomial has a positive zero . As already mentioned, from (7a)–(7e), as long as , we acquire that an eigenvalue λ with exists if and only if first , and additionally with (9a) and (9b), fulfills the conditionFrom Ref. 23, we know thatThis implies that
For in order to change its sign, we finally examine and get the condition (18), which is sufficient to ensure a positive zero of .
Similarly, considerandThen, (19) and (20) follow directly with the same argumentation.
Corollary 3.4. Assume that the system (3) satisfies the condition (10) and either or are chosen sufficiently large. Then, the instability condition (18) is satisfied if the following conditions hold:
Case 1:and for exists an such that .
Case 2:and with as defined above exists an with .
Remark 3.5. If and the system fulfills condition (12), then the instability condition (20) holds for sufficiently large , if the following conditions are satisfied:
Case 1:and forexists an withCase 2:and with as defined above exists an with
Remark 3.6. If and the system fulfills condition (12), then the instability condition (20) holds for sufficiently large if the following condition is satisfied.
Case 1:and forexists an with
Case 2:and with as defined above exists an withProof. We first restrict ourselves to and the case as well as . In order to achieve instability, we consider (18) and narrow down to the coefficient of which is given byWe definewhose roots are given byIn order to satisfy condition (18) and to obtain an instability, we now require . First, assume and . Then, e represents a right displaced upward open parabola which intersects the positive axis at points . Thus, with to ensure , we have to satisfy the conditions of case 1. By contrast, if is negative, then the parabola is shifted to the left, and we directly prove case 2 to obtain .
We further consider as well as the case . Since we suppose that , as before, we observe that either or becomes dominant in (18). This implies that an instability exists for sufficiently large or .
Finally, with the same argumentation as before, the analysis of the coefficient in (19) as well as in (19) deduces Remarks 3.5 and 3.6 (for the case , see also Ref. 23).
In contrast to the model of Ref. 23, the conditions in our model depend on as well as the capacity function 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.
4 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 first-order 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 dune-pdelab package [4]. Coupled bulk-surface problems are solved by the DUNE modules multidomain and multidomaingrid [20].
4.1 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 byThe particular choice of parameters for the numerical simulation of the nondimensionalized system (1) is given in Table 1. The model parameters are estimated from values reported in the literature. We based our choice on published results in Refs. 11–14, 25, and 32. A full list of these parameters is given in the Supplementary Material. The effective parameters were then obtained by rescaling with respect to the membrane diffusion coefficient , which relates to the parameter in the nondimensional formulation.
TABLE 1
| Param. | Value | Description |
|---|---|---|
| 1.0 | Diffusion coefficient of the inactive membrane-bound species | |
| 27.78 | Diffusion coefficient of the actin cables | |
| 305.5 | Diffusion coefficient of the cytosolic component | |
| 0.278 | Diffusion coefficient of the internal component | |
| Γ | 15.6 | Spatial scale factor |
| 0.056 | Basal activation rate | |
| 27.78 | Basal inactivation rate | |
| 0.025 | Basal membrane attachment rate | |
| 13.89 | Basal membrane detachment rate | |
| 10.77 | Feedback mediated membrane attachment rate | |
| 19.41 | GEF mediated activation rate | |
| 59.42 | GEF mediated feedback activation | |
| 2.23 | Negative feedback/inactivation rate | |
| 0.703 | Rate of exocytosis | |
| 8.33 | Rate of endocytosis | |
| 0.0163 | Rate controlling local endocytosis | |
| B | 27.78 | Transport gradient control rate |
| V | 1.53 | Capacity function control rate |
| Α | 419.94 | Potential flow control rate |
| 0.004 | Initial concentration of the active membrane-bound species | |
| 0.196 | Initial concentration of the inactive membrane-bound species | |
| 0.002 | Initial concentration of the actin cable density | |
| 1.504 | Initial concentration of the internal component | |
| 0.902 | Initial concentration of the cytosolic component |
4.2 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 well-mixed 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 defineIn 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 and decrease 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 and .
In all computations, we use functions 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 . This implies that the only requirement on is that it must be nonzero. In this case, the capacity function determines whether the system is stable against small perturbations or not. We further call this capacity-driven instability.
Figure 1 shows the development of u in time for distinct values of . We observe that even for large changes of , provided that , the system is always unstable and tends to form a polarized patch. It becomes clear that the capacity function as well as determines the stability behavior. The constant 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 , the system is still able to form a polarized patch, albeit after a very long time ().
FIGURE 1
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 as well as the total mass is visualized in Figure 2.
FIGURE 2
Another result of the stability analysis is the fact that we may observe polarization, even if or . 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.
FIGURE 3
The requirement yields that must not be very large to ensure instability. We have seen that even in the case , the instability conditions may be satisfied. Our numerical simulations confirm these results. Figure 4 illustrates capacity-driven polarization for the system (3), where .
FIGURE 4
4.3 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 three-dimensional 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.
FIGURE 5
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 and of the same intensity. The signals are imposed on opposite sides of the cell surface, one 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 diffusion-driven 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.
4.4 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 and of the same intensity. Whereas one signal is located near the organelles, the other is placed at the opposite side.
FIGURE 6
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 of strength .
FIGURE 7
5 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 diffusion-driven 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 diffusion-driven 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–diffusion–transport 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 and LatA-treated cells. Treating wild-type yeast cells with latrunculin A (LatA) removes the actin-dependent recycling pathway, while 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, , the capacity function 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 capacity-driven instability also generates pattern when the cytosolic and membrane diffusion rates are equal. For that reason, and since the diffusion constant has no essential impact on the stability of the system, we assert that this instability mechanism distinguishes from the Turing-type instability.
Funding
The content of this manuscript has been published as part of the thesis of [6]. NE was supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy EXC 1003, Cells in Motion, CiM, Münster. CE was supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy EXC 2044-390685587, Mathematics Münster: Dynamics–Geometry–Structure.
Statements
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.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fams.2020.570036/full#supplementary-material
References
1.
AyscoughKRStrykerJPokalaNSandersMCrewsPDrubinDG. High rates of actin filament turnover in budding yeast and roles for actin in establishment and maintenance of cell polarity revealed using the actin inhibitor latrunculin-A. J Cell Biol (1997). 137:399–416. 10.1083/jcb.137.2.399.
2.
BastianPBlattMDednerAEngwerCKlöfkornRKornhuberRet alA generic grid interface for parallel and adaptive scientific computing. part II: implementation and tests in DUNE. Computing (2008). 82:121–38. 10.1007/s00607-008-0004-9.
3.
BastianPBlattMDednerAEngwerCKlöfkornROhlbergerMet alA generic grid interface for parallel and adaptive scientific computing. part I: abstract framework. Computing (2008b). 82:103–19. 10.1007/s00607-008-0003-x
4.
BastianPHeimannFMarnachS. Generic implementation of finite element methods in the distributed and unified numerics environment. Kybernetik (2010). 46, 294–315.
5.
DagaRRChangF. Dynamic positioning of the fission yeast cell division plane. Proc Natl Acad Sci USA (2005). 102, 8228–32. 10.1073/pnas.0409021102
6.
EmkenN. A coupled bulk-surface reaction-diffusion-advection model for cell polarization. [PhD thesis]. Münster University of Münster (2016).
7.
EmkenNPüschelABurgerM. Mathematical modelling of polarizing gtpases in developing axons. arXiv preprint arXiv:1204.5725 (2012).
8.
Etienne-MannevilleSHallA. Rho gtpases in cell biology. Nature (2002). 420, 629–35. 10.1038/nature01148.
9.
FreisingerTKlünderBJohnsonJMüllerNPichlerGBeckGet alEstablishment of a robust single axis of cell polarity by coupling multiple positive feedback loops. Nat Commun (2013). 4, 1–11. 10.1038/ncomms2795
10.
GieseWEigelMWesterheideSEngwerCKlippE. Influence of cell shape, inhomogeneities and diffusion barriers in cell polarization models. Phys Biol (2015). 12, 066014. 10.1088/1478-3975/12/6/066014
11.
GoryachevABPokhilkoAV. Dynamics of cdc42 network embodies a turing-type mechanism of yeast cell polarity. FEBS Lett (2008). 582, 1437–43. 10.1016/j.febslet.2008.03.029
12.
HowellASSavageNSJohnsonSABoseIWagnerAWZylaTRet alSingularity in polarization: rewiring yeast cells to make two buds. Cell (2009). 139, 731–43. 10.1016/j.cell.2009.10.024
13.
JohnsonJMJinMLewDJ. Symmetry breaking and the establishment of cell polarity in budding yeast. Curr Opin Genet Dev (2011). 21, 740–6. 10.1016/j.gde.2011.09.007
14.
KlünderBFreisingerTWedlich-SöldnerRFreyE. Gdi-mediated cell polarization in yeast provides precise spatial and temporal control of cdc42 signaling. PLoS Comput Biol (2013). 9:e1003396. 10.1371/journal.pcbi.1003396.
15.
KnabnerPAngermanL. Numerical methods for elliptic and parabolic partial differential equations.New York, NY: Springer (2003).
16.
KraussGSchönbrunnerNCooperJ. Biochemistry of signal transduction and regulation. (Bayreuth: Wiley Online Library), Vol. 3 (2003).
17.
LaytonATSavageNSHowellASCarrollSYDrubinDGLewDJ. Modeling vesicle traffic reveals unexpected consequences for cdc42p-mediated polarity establishment. Curr Biol (2011). 21, 184–94. 10.1016/j.cub.2011.01.012
18.
MadzvamuseAChungAHWVenkataramanC. Stability analysis and simulations of coupled bulk-surface reaction-diffusion systems. Proc. R. Soc. A (2015). 471, 20140546. 10.1098/rspa.2014.0546
19.
MataJNurseP. tea1 and the microtubular cytoskeleton are important for generating global spatial order within the fission yeast cell. Cell (1997). 89, 939–49. 10.1016/s0092-8674(00)80279-2
20.
MüthingS. A flexible framework for multi physics and multi domain PDE simulations. [PhD thesis]. Institut für Parallele und Verteilte Systeme der Stuttgart: Universität Stuttgart (2015). 10.18419/opus-3620
21.
RätzA. Turing-type instabilities in bulk-surface reaction-diffusion systems. J Comput Appl Math (2015). 289:142–52. 10.1016/j.cam.2015.02.050
22.
RätzARögerM. Turing instabilities in a mathematical model for signaling networks. J Math Biol (2012). 65:1215–44. 10.1007/s00285-011-0495-4
23.
RätzARögerM. Symmetry breaking in a bulk-surface reaction-diffusion model for signalling networks. Nonlinearity (2014). 27:1805–27. 10.1088/0951-7715/27/8/1805
24.
SagotIKleeSKPellmanD. Yeast formins regulate cell polarity by controlling the assembly of actin cables. Nat Cell Biol (2002). 4:42–50. 10.1038/ncb719
25.
SavageNSLaytonATLewDJ. Mechanistic mathematical model of polarity in yeast. Mol Biol Cell (2012). 23:1998–2013. 10.1091/mbc.e11-10-0837
26.
SiegristSEDoeCQ. Microtubule-induced cortical cell polarity. Genes Dev (2007). 21:483–96. 10.1101/gad.1511207
27.
SlaughterBDDasASchwartzJWRubinsteinBLiR. Dual modes of cdc42 recycling fine-tune polarized morphogenesis. Dev Cell (2009). 17:823–35. 10.1016/j.devcel.2009.10.022
28.
SlaughterBDSmithSELiR. Symmetry breaking in the life cycle of the budding yeast. Cold Spring Harb Perspect Biol (2009). 1:a003384. 10.1101/cshperspect.a003384
29.
SymonsMH. Rho GTPases. Molecular biology intelligence unit. Springer (2004).
30.
TahirovicSBradkeF. Neuronal polarity. Cold Spring Harb Perspect Biol (2009). 1:a001644. 10.1101/cshperspect.a001644
31.
TanosBRodriguez-BoulanE. The epithelial polarity program: machineries involved and their hijacking by cancer. Oncogene (2008). 27:6939–57. 10.1038/onc.2008.345
32.
Wedlich-SoldnerRAltschulerSWuLLiR. Spontaneous cell polarization through actomyosin-based delivery of the cdc42 gtpase. Science (2003). 299:1231–5. 10.1126/science.1080944
33.
Wedlich-SöldnerRWaiSCSchmidtTLiR. Robust cell polarity is a dynamic state established by coupling transport and GTPase signaling. J Cell Biol (2004). 166:889–900. 10.1083/jcb.200405061
Summary
Keywords
polarization models, spatial simulation, spatial inhomogeneities, Cdc42, yeast, surface PDEs, advection diffucions reaction systems, pattern formation
Citation
Emken N and Engwer C (2020) A Reaction–Diffusion–Advection Model for the Establishment and Maintenance of Transport-Mediated Polarity and Symmetry Breaking. Front. Appl. Math. Stat. 6:570036. doi: 10.3389/fams.2020.570036
Received
05 June 2020
Accepted
18 September 2020
Published
16 November 2020
Volume
6 - 2020
Edited by
Sebastian Aland, Hochschule für Technik und Wirtschaft Dresden, Germany
Reviewed by
Anotida Madzvamuse, University of Sussex, United Kingdom
Andreas Rätz, Technical University Dortmund, Germany
Updates
Copyright
© 2020 Emken and Engwer.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Christian Engwer, christian.engwer@uni-muenster.de
This article was submitted to Mathematical and Statistical Physics, a section of the journal Frontiers in Applied Mathematics and Statistics
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.