MATHEMATICAL ANALYSIS OF STEADY-STATE SOLUTIONS IN COMPARTMENT AND CONTINUUM MODELS OF CELL POLARIZATION

Cell polarization, in which substances previously uniformly distributed become asymmetric due to external or/and internal stimulation, is a fundamental process underlying cell mobility, cell division, and other polarized functions. The yeast cell S. cerevisiae has been a model system to study cell polarization. During mating, yeast cells sense shallow external spatial gradients and respond by creating steeper internal gradients of protein aligned with the external cue. The complex spatial dynamics during yeast mating polarization consists of positive feedback, degradation, global negative feedback control, and cooperative effects in protein synthesis. Understanding such complex regulations and interactions is critical to studying many important characteristics in cell polarization including signal amplification, tracking dynamic signals, and potential trade-off between achieving both objectives in a robust fashion. In this paper, we study some of these questions by analyzing several models with different spatial complexity: two compartments, three compartments, and continuum in space. The step-wise approach allows detailed characterization of properties of the steady state of the system, providing more insights for biological regulations during cell polarization. For cases without membrane diffusion, our study reveals that increasing the number of spatial compartments results in an increase in the number of steady-state solutions, in particular, the number of stable steady-state solutions, with the continuum models possessing infinitely many steady-state solutions. Through both analysis and simulations, we find that stronger positive feedback, reduced diffusion, and a shallower ligand gradient all result in more steady-state solutions, although most of these are not optimally aligned with the gradient. We explore in the different settings the relationship between the number of steady-state solutions and the extent and accuracy of the polarization. Taken together these results furnish a detailed description of the factors that influence the tradeoff between a single correctly aligned but poorly polarized stable steady-state solution versus multiple more highly polarized stable steady-state solutions that may be incorrectly aligned with the external gradient.


1.
Introduction.Breaking symmetry is fundamental to all biological systems [29].Components that were previously uniformly distributed become asymmetrically localized.This anisotropy or polarization creates specialized structures that produce complex behaviors.Cell polarization or anisotropy is involved in the differentiation, proliferation, morphogenesis of organisms and activation of the immune response [20,2,32].A key challenge during these processes is robust polarization: polarizing in the right direction, at the right time, and to the proper extent under uncertain conditions.
Internal and external cues direct cells to localize components to specific cellular locations leading to morphological changes.For example, haploid cells of the yeast Saccharomyces cerevisiae typically form a new bud at the site of the previous bud, which acts as an internal cue.In addition, haploid yeast cells can sense an external gradient of mating pheromone and form a mating projection toward the source.In both cases, a large number of proteins adopt a polarized distribution, being concentrated at the site of the morphological change [7,25].
Cell polarization can be thought of as a type of pattern formation.Turing originally proposed that complex spatial patterns could arise from simple reactiondiffusion systems [31].In particular, Meinhardt demonstrated that local positive feedback balanced by global negative feedback could give rise to cell polarization [18].More recently, researchers have constructed detailed mechanistic models in which specific molecular species and reactions are represented.One popular class of models employs a local excitation, global inhibition (LEGI) mechanism [10,12].
From a biology perspective, it is known that the cell polarity behavior is quite robust [3], in the sense that the polarization can be established even under very shallow gradient.In the literature, the focus has been on understanding how a shallow external gradient can be amplified to create a steep internal gradient of cellular components.High amplification can result in an all-or-none localization of the internal component to a narrow region.Detailed biochemical models are proposed in [16,11,24,26,13,9,15,27] and reviewed in [5,12,4], while some models aim to account for the symmetry breaking [19,28,22,24,8].
In addition to the establishment of polarity, the tracking of a moving signal source has also been acknowledged to be important.Devreotes and colleagues [5] made the distinction between directional sensing (low amplification, good tracking) and polarization (high amplification, poor tracking).Meinhardt first highlighted the potential tradeoff between amplification and tracking [19].Dawes et al. categorized some models according to gradient sensing, amplification, polarization, tracking of directional change, persistence when the stimulus is removed (i.e.multi-stability) ( [4] and references therein).These models varied in the degree of amplification (polarization), presence of multiple steady states, response to a rotating gradient, etc.
While mathematical modeling provides great insight into how this robustness is achieved and sheds light on the tradeoff between polarization and tracking, simple models are particularly favorable because it permits more rigorous theoretical investigations.Most of the literature and work on mathematical analysis of the models of cell polarization have mainly focused on the establishment and maintenance of polarity, without emphasis on the tracking of the stimuli.
Compartmental analysis has also been frequently used for analyzing models [1].The substance with a spatial distribution can be considered as distributed among a number of separate and connected compartments.The dynamics of the substance within the system is then described by ordinary differential equations in each compartment, allowing to obtain more quantitative information of the entire system.To explain both adaptation to uniform increases in chemoattractant and persistent signaling in response to gradients, Levchenko et al. [14] put forth a set of differential equations and analyzed the steady-state solutions by investigating the algebraic equations of the associated steady-state system.Recently, Mori et al. studied a simple system composing of a single active/inactive Rho protein pair with cooperative positive feedback and conservation requirement [21], based on a single unified system of actin, Rho GTPases and PIs in [4].Through analysis, the authors [21] elucidated the phenomenon of wave-pinning and demonstrated how it could account for spatial amplification and maintenance of polarity as well as sensitivity to new stimuli typical in polarization of eukaryotic cells.
Previously, we constructed both generic and mechanistic models of yeast cell polarization in response to mating pheromone spatial gradients [3], in which we used only numerical simulations of a system of reaction-diffusion equations in continuum version to demonstrate the tradeoff between the amplification necessary to tightly localize proteins at the front of the cell and the tracking necessary to follow a change in the gradient direction.In this paper, we analyze several models with similar structures but in different and simple spatial description: two-compartment,threecompartment, and continuum space.Using this approach, we are able to investigate both analytically and numerically on several important characteristics and properties of the system that have not been carefully explored in [3].One of the key features in the model is the existence of multiple steady state solutions for the system in the absence of membrane diffusion, showing the important role of membrane diffusion in regulating polarity.We systematically explore and compare three types of models with and without diffusions in this paper to examine the important factors that influence the tradeoff between a single correctly aligned but poorly polarized stable steady-state solution versus multiple more highly polarized stable steady-state solutions that may be incorrectly aligned with the external gradient.
This paper is organized as follows.In Section 2, we introduce models with two spatial compartments having different forms of positive feedback; in Section 3, we study the three-spatial-compartment version of some of the models in Section 2; in Section 4, we investigate the continuum model with the same mechanisms as in the compartment models; in Section 5, we discuss and summarize. .Schematic descriptions of two-compartment and continuous models of cell polarization.In both, the polarized species x 1 (red) becomes localized to the front of the cell through cooperative interactions (red arrow, P f ) in response to the input and through positive feedback (red +).There is a global negative feedback mediated by the species x 2 (blue).(A) Two-compartment model of cell polarization.The cell is divided into compartment 1 (front where the ligand concentration is higher) and compartment 2 (back).The species x 1 polarizes to the front (x f 1 ) and is less abundant at the back (x b 1 ).(B) Continuous model of cell polarization.The input gradient and spatial localization of x 1 is represented in a continuous fashion.
2. Two-compartment models.Throughout this paper, we will consider the interaction of two intracellular species x 1 and x 2 .Species x 1 is a membrane bound protein that polarizes when exposed to a ligand gradient input, and x 2 is a global inhibitor of x 1 which is homogeneous in space.The basic dynamics in the system include membrane diffusion of x 1 , cooperative production of x 1 , positive feedback of x 1 , degradation of x 1 and x 2 , and global inhibition of x 1 by x 2 (see Fig. 1 for illustration).The simplest model that accounts for the basic dynamics of x 1 and x 2 is a two-compartment model, in which the cell is divided into two compartments: front and back.Herein "front" refers to the end exposed to higher ligand concentration, with the back to lower ligand concentration.In this case, the front end will be where x 1 localizes.In the two-compartment setting, the membrane diffusion corresponds to a linear transport of x 1 between the compartments, and the gradient dependent cooperative production of x 1 corresponds to a higher production of x 1 at the front and lower at the back compartment.
In the following subsections, four models with different forms of positive feedback will be discussed.In those models, superscripts f and b will be used to distinguish x 1 at the front and back compartments, respectively.Since x 2 is homogeneous in space, in the two-compartment model, we have totally three variables x f 1 , x b 1 and x 2 .The rate of transport of x 1 from the back to the front is denoted by D b , and D f for the rate from the front to the back.This transport between compartments corresponds to the surface lateral diffusion in a continuum setting, so we will refer to this inter-compartment transport by the term "diffusion" in this paper.The cooperative production of x 1 induced by the ligand gradient at the front and back compartments are denoted by P f and P b , respectively.Therefore, the case of P f higher than P b corresponds to that in which the ligand gradient is from the back to the front.The degradation rate of x 1 is denoted by k 2 , and the negative feedback term is k 3 x 1 x 2 , which depends on the global inhibitor x 2 as well as x 1 .Parameter k 1 modulates the rate of positive feedback.The rate of change of the global inhibitor x 2 is k 4 , and it is proportional to the difference between a constant k ss and the averaged amount of x 1 throughout the cell.This control of x 2 will be referred to as "integral constraint", which regulates x 2 according to the total amount of x 1 .
2.1.Model 1A.The first model we consider has the positive feedback of x 1 (normalized by a constant k ss ) taking place in an exponential fashion, with the power h.The dynamics of x f 1 , x b 1 and x 2 can be described by the following ordinary differential equations: It is expected that the extent of polarization highly depends on, beside all other parameters, the power h in the positive feedback term, which is the major spatial amplification mechanism in our model that localizes x 1 to a narrow region in the front.As will be seen in the following analysis, h also dictates the number of steady states of the system: the higher h is, the more steady-state solutions there are.Since we are interested in the polarized solution with nonzero x 2 , we assume x 2 , k 4 > 0; in addition, we fix k ss to be 1 to focus on the role of diffusion, cooperativity and positive feedback on polarization.We will use h = 1, 2 to study the steady states of the system.
2.1.1.Linear case: h = 1.The steady state is unique, and the solutions of x f 1 and x b 1 are To examine the polarization of x 1 , we take the difference of x f 1 and x b 1 It can be seen that whether x 1 is higher at the front or at the back of the cell depends on the balance between the difference of production and the difference of diffusion between the two compartments.Since the cooperative production is induced by the pheromone gradient, we have P f > P b when the gradient is along the back-to-front direction.In particular, if the diffusion is spatially homogeneous, namely 1 , a solution polarizing at the front.Generally, as long as (P f − P b ) > 2(D f − D b ), x 1 will polarize at the front of the cell.This expression also reveals that diffusion counteracts the input-dependent cooperative production term.
2.1.2.Quadratic case: h = 2.The steady-state equations of Eqs. ( 1)-( 3) can be reduced to a cubic equation of x f 1 , whose three roots are where and This system has multiple real steady states if and only if s 1 = s 2 , namely, d = 0. Due to the complexity of the formula, the explicit forms of the conditions under which the system has a unique real steady state are difficult to obtain.However, knowing that this model has at most three steady states, one observes that if a 1 > 3, i.e, D f +D b +(P f +P b )/2 > k 1 , d will be positive, and consequently this system has only one real solution.In other words, if the diffusion and the cooperative production P f or P b is strong enough compared to the positive feedback, the system has a unique steady state.This expression illustrates how the number of steady states in the quadratic case depends on diffusion, input-dependent polarized production, and the positive feedback.
In order to understand the behavior of steady-state solutions of this system, we numerically solve the steady-state equations, and evaluate the local stability of each steady state by computing the eigenvalues of its Jacobian.Fig. 2 shows some examples when varying parameters P f , P b , D f , D b and k 1 , with other parameters fixed.According to Fig. 2 and other extensive numerical simulations not shown here, if we further define c 1 = P f k1 , c 2 = P b k1 , several numerical observations were obtained: (1) three real steady-state solutions appear only when c 1 < 1, c 2 < 1; (2) there is at most one stable steady-state solution, which occurs only when c 1 > 1 or c 2 > 1; (3) as k 1 increases (i.e.c 1 decreases), there are more steady-state solutions; (4) small diffusion rates D f , D b result in multiple steady-state solutions.Thus, the ratio of input-dependent cooperative production to the positive feedback is a key parameter.), where k T is a constant.The system of equations is shown in Eqs. ( 4)- (6).The whole k 1 term is positive when < k T and negative when > k T .Therefore, k T acts as a threshold at which the regulation is positive when the average amount of x 1 is lower than k T , and it is negative otherwise.In this manner, the additional term prevents the positive feedback from growing in an unstable fashion.
), ( 4) If we only consider steady-state solutions with nonzero x 2 , then has to be k ss , and the term k ) will be a constant k1 ≡ k 1 (k T − k ss ).In other words, at steady states, Model 1B is essentially the same as Model 1A except that the rate of positive feedback term is scaled by the constant (k T − k ss ).Therefore, all the conclusions for the number of steady states of Model 1A can be applied to this model by replacing k 1 in Model 1A with k1 ≡ k 1 (k T − k ss ).
When h = 1 and without diffusion, Model 1B has a unique solution which preserves the monotonicity of the ligand concentration at the front and back compartments.When h = 2, multiple steady states may arise, with at most 3 steady states, and the steady state is unique while D f + D b + (P f + P b )/2 > k1 .We performed numerical simulations with varying parameters P f , P b , D f , D b , k 1 , and the steady states are shown in Fig. 3.All the parameters used are the same as in Fig. 2, except k = 1 is used for Fig. 3.The steady states in Fig. 2 and 3 are the same, and we also have the same observations as for Model 1A: (1) three real steady-state solutions appear only when c1 < 1, c2 < 1 (c 1 = P f k1 , c2 = P b k1 ); (2) a unique stable steady-state solution occurs only when c1 > 1 or c2 > 1; (3) as k1 increases (i.e.c1 decreases), there are more steady-state solutions; (4) small diffusion rates D f , D b result in multiple steady-state solutions.
Although the behavior of steady states is similar to Model 1A, Model 1B is very different from Model 1A in the local stability of the steady states.It can be observed in Fig. 3 that the number of stable steady states (blue symbols) is much more than in Fig. 2 for Model 1A, which implies that this model has more admissible solutions.Extensive numerical simulations also reveal that there are at most 2 stable steady-state solutions, which occurs only when c1 < 1, c2 < 1. 2.3.Model 2A.In this subsection, we consider a model with the positive feedback term in a Hill function form, which possesses the Hill exponent h and the Hill halfmaximal constant 1/γ.In this manner, we replaced the exponential form of the positive feedback term with a Hill form that is more common to biological reaction descriptions.This positive feedback achieves its maximal value k 1 as x 1 approaches infinity, and it assumes its minimal value 0 as x 1 approaches 0. When h is large, the feedback response becomes switch-like.The system is as follows: We consider two cases when h = 1 or h = 2.
2.3.1.Linear case: h = 1.The steady-state equations can be reduced to a cubic equation and the relation where Since a 1 will always be positive, the leading coefficient of the polynomial will be nonzero.Hence, there are at most three steady states.
The leading coefficient of the polynomial is , which is always negative; therefore, the above equation has at most 5 steady-state solutions.
Due to the complexity of the coefficient of Eq. ( 10), we directly solve the steadystate system (7)-( 9) numerically with MATLAB, and analyze the local stability of the steady states.The results of varying parameters P f , P b , D f , D b , k 1 is displayed in Fig. 4. With the definitions c 1 = P f k1 , c 2 = P b k1 , we made four numerical observations: (1) under some parameter sets, we did find five real steady states, all observed when c 1 < 0.1, c 2 < 0.1.(2) at most three stable steady state are observed, which occurs only when c 1 < 0.1, c 2 < 0.1; therefore, Model 2A not only has more steady states, but also more stable steady states than Models 1A and 1B; (3) as k 1 increases (i.e.c 1 decreases), there are more steady-state solutions, which is also observed in Models 1A and 1B; (4) small diffusion rates D f , D b result in multiple steady-state solutions, which is also observed in Models 1A and 1B.
2.4.1.Linear case: h = 1.The steady-state equation of The above equation is a cubic equation, so there are at most 3 steady-state solutions.
2.4.2.Quadratic case: h = 2.The steady-state equation is, 1 , so there are at most 5 steady-state solutions.By solving the steady state system (11)-( 13) directly, numerical simulations with different P f , P b , D f , D b , k 1 are displayed in Fig. 5. Five real steady-state solutions were found under some parameter sets as shown in Fig. 5C.Different from Model 2A, this model requires larger k 1 to obtain 5 steady-state solutions and 3 stable steady-state solutions.In general, for high k 1 , Model 2B exhibited strong polarization (i.e.x f 1 ≈ 2) but a reduced region of multi-stability than Model 2A.Similar to all Models 1A, 1B and 2A, it is observed that small diffusion rates result in multiple steady-state solutions, and a larger k 1 results in more steady states.3. Three-compartment models.In this section, we study three-compartment models in which the cell is divided into three segments: front, middle and back.The mechanisms included in our models are the same as in Section 2, but the increase in the number of spatial compartments provide greater spatial detail while still being analytically approachable.The concentration of x 1 at the front, middle and back compartments are denoted by x f 1 , x m 1 and x b 1 , respectively.The diffusion rate of x 1 from the front/back to the middle compartment is D f /D b ; the rates at which x 1 is transported from the middle to the front and back compartment are denoted by D mf and D mb , respectively.The cooperative production of x 1 at the three compartments are P f , P m and P b .All the other parameters have the same definitions as in Section 2.
In the rest of this section, we study two models that assume the same form of positive feedback, as in Model 1A and 2B in Section 2. We still call those threecompartment versions Model 1A and Model 2B without confusion.
3.1.Model 1A.With the positive feedback in an exponential form, a threecompartment Model 1A can be described by the following system: In the following subsections, we will study the cases of h = 1 and h = 2, which correspond to different strengths of positive feedback.
3.1.1.Linear Cases: h = 1.For the simplicity of notation, we denote the sum of the cooperative production to be a positive constant c, namely, If one assumes that the diffusion rates between the compartments are uniform, namely, D mf = D mb = D f = D b = D, then the above solutions can be simplified as x m By Eqs. ( 18)- (20), we are able to conclude that the monotonicity and linearity of the cooperative production P f , P m , P b are correlated with those of x f 1 , x m 1 , x b 1 , with the proof in Proposition 1.In other words, a graded external signal would result in a graded response of x 1 , which guarantees polarization in the correct direction, and that is a desirable property for the gradient sensing model because the polarization with the input gradient is in the correct direction.Proposition 1. Suppose the diffusion rates between each compartment are uniform, namely, Proof.Using Eqs. ( 18)-( 20), one gets Therefore, the monotonicity of the input to the system is preserved at the steady state.
In particular, if P f − P m = P m − P b = α, then it can be easily seen that Next, we would like to examine how the diffusion rate D affects the polarization.Intuitively, faster diffusion of the molecules inhibits the accumulation of the molecules, and one would expect a decrease in polarization when the diffusion rate is enhanced.Proposition 2 proves the above statement, in which we define the extent of the polarization by the difference of x 1 at the front versus the back compartment.Moreover, if the cooperative production is linear, increasing D does not change the response in the middle compartment (x m 1 ), but only decreases the polarization (x f 1 − x b 1 ).Proposition 2. Suppose the diffusion rate between each compartment is uniform, i.e.
decreases.In particular, if P f , P m , P b are linear with P f − P m = P m − P b , x m 1 will be invariant with respect to D.
Proof.By Eqs. ( 18) and ( 20), we have 3D+c , and therefore if As D decreases, the difference of x 1 at the front and the back will increase.
By Eq. ( 19), x m 1 = 3(3D+Pm) 9D+c .Using the relation P f + P m + P b = c, one can easily verify the following conclusions: 1) if 1 is invariant with respect to D. Case (3) corresponds to the situation when the cooperative production P f , P m , P b are linear.The analysis tells us that in that case, changing D does not affect x m 1 but only affects the polarization of x 1 at the front and back.
3.1.2.Quadratic Cases: h = 2. Throughout this subsection, we will assume D f = D mf = D mb = D b for the simplicity of analysis.We first investigate a case when there is no diffusion, in which the communication between the compartments is merely through the global integral control.The no-diffusion case renders a relatively simple system to analyze, and we prove in the following proposition that there will be at most 9 steady states with D = 0.
x 2 by x, y, z, w respectively to avoid super-and subscripts in the equations, we have the following properties for the steady-state equations with D = 0.
• If 2y 2 − 3y + c 2 = 0, the steady-state system is equivalent to the following system: • If 2y 2 − 3y + c 2 = 0 in system (22), the system is consistent only if c 1 = c 3 and it has real solutions only if c 2 ≤ 9 8 .• There are at most 9 steady-state solutions.
Proof.The steady-state equations of Eqs. ( 14)-( 17) when h = 2 with D = 0 are: x • Summing Eqs. ( 23)-( 25), one gets Using ( 23) and ( 24), we have With Eqs. ( 26) and ( 27), Eq. ( 24) becomes 2y 3 + 2x 2 y + 2xy 2 − 9y 2 − 6xy + 9y + ( If we further eliminate x 2 y in (29) and use (28), we get Substitute x in Eq. ( 28) with Eq. ( 30), one gets the following equation for y: After expansion, it becomes Therefore, the solutions of the steady-state equations ( 23)-( 25) are solutions of system (22).Conversely, it is easy to verify that if 2y 2 −3y +c 2 = 0 and if (x, y, z) is a solution of system (22), then it is a solution to the steady-state equations ( 23)-( 26).Furthermore, if system (22) has real solutions under the condition 2y 2 − 3y + c 2 = 0, then 9 − 8c 2 ≥ 0, i.e. c 2 ≤ 9  8 .• All the solutions of system ( 23)-( 26) must satisfy Eq. (32).So system ( 23)- (26) has at most 7 solutions for y.If 2y 2 − 3y + c 2 = 0, one y corresponds to one x based on Eq. ( 30).If 2y 2 − 3y + c 2 = 0, each y corresponds to at most 2 real solutions for x according to Eq. ( 28), so system ( 23)-( 26) has also at most 9 real solutions (x, y, z, w). ).Other parameters used are as follows: The solutions are evaluated within the range 0.2 ≤ c 1 , c 2 , c 3 ≤ 2, with discretized space 0.2.Colors of red, magenta, yellow, green, cyan, blue and black stand for 1, 2, 3, 4, 5, 6, 7 real positive solutions respectively.Although Proposition 3 provides clues about the number of steady states when D = 0, due to the complexity of the system, one still relies on numerical simulations to obtain more details about how the number of steady states changes with respect to different parameters.In Fig. 6, the number of steady states is evaluated within a range of parameters c 1 , c 2 , c 3 , with different colors indicating different number of solutions.It can be observed from Fig. 6 that without diffusion, when c 1 , c 2 , c 3 are all very small, there are more steady-state solutions (up to 7 steady states for the parameters we used); as one of the c j 's is increased, the number of steady states decreases.This implies that either increasing the cooperative production or decreasing the positive feedback can reduce the number of steady states.
Next, we studied how the solution of x f 1 changes with parameters (P f , P b , k 1 and D f , D mf , D mb , D b ) in Fig. 7, without assuming diffusion rates zero.According to Fig. 7 and other numerical simulations not shown here, we made the following observations: (1) Model 1A has 7 steady-state solutions only when c 1 , c 2 , c 3 are less than 1; (2) steady-state solutions are stable only when c 1 > 1; no more than 2 stable steady-state solutions are found, and 2 stable steady states occur only when P f = P m ; (3) as k 1 increases, the number of steady-state solutions increases but the number of stable solutions decreases; (4) small diffusion results in more steady-state solutions.In particular, Fig. 7D shows that an increase in the diffusion rates results in a decrease of the number of steady states, and also results in more stable steady states.This implies that the diffusion improves the system by both reducing the number of steady states and increasing their stability.Comparing three-compartment and two-compartment models with the exponential form of feedback (Model 1A), it is found from our numerical tests that they are mainly different in the total number of steady-state solutions and the number of stable steady states: (1) there are at most 7 steady-state solutions for the three-compartment model and at most 3 steady-state solutions for the twocompartment model; (2) there are at most 2 stable steady-state solutions for the three-compartment model and at most 1 stable steady-state solution for the twocompartment model.
However, these two models have more in common: (1) the number of steadystate solutions is reduced as P b or P m increases ; (2) the maximal number of steady states occurs only when c 1 , c 2 , c 3 < 1; (3) there are stable steady states only when c 1 > 1; (4) as k 1 increases, there are more steady-state solutions and fewer stable steady-state solutions; (5) as the diffusion increases, the number of steady-state solutions is reduced, and there are more stable steady states.
Next, we consider if the monotonicity of the input P f , P m , P b can be preserved by the model.The simulations with all diffusion being zero is shown in Fig. 8A, in which there are 7 steady states for the set of parameters simulated.Only one steady state out of these seven is polarized toward the right direction, while the others either have the maximum at the middle compartment or polarize at the back.However, when the diffusion rates D f , D mf , D mb , D b are increased, the total number of steady states is reduced (Fig. 8B-8D).At a high diffusion rate such as only the front-polarizing solution is found.In summary, in the quadratic case h = 2, the monotonicity of the input may not be preserved as in the linear case, and increasing the diffusion could reduce the number of steady states, enhance the stability, and at the same time select the correctly polarized solution.
, and the non-polarizing solution (x

Model 2B.
With the positive feedback loop implemented in the Hill form and having the feedback/feedforward coincidence detection mechanism, the threecompartment model can be described by the following equations: We numerically investigated the steady states of the system with h = 2, shown in Fig. 9.By comparing Fig. 5 and Fig. 9, we can compare Model 2B with twocompartment and three-compartment.It can be observed that the general behavior of these two models are quite similar except that the three-compartment model has more steady states and more stable solutions than the two-compartment model.More precisely, in the two-compartment Model 2B, we found up to 5 steady states and 3 stable steady states, and in the three-compartment Model 2B, we found up to 17 steady states and 6 stable steady states.
If one compares the three-compartment Models 1A and 2B by comparing Fig. 7 and Fig. 9, it can be seen that there are more steady states found in Model 2B than 1A under the same sets of parameters.

Diffusion barrier at the front compartment enhances polarization.
So far, all the analysis and numerical simulations conducted are based on the "uniform diffusion" scenario, even though our model does not require the diffusion rates to be the same.However, there is abundant evidence that the plasma membrane is quite heterogeneous and that lateral diffusion can be restricted by the cytoskeleton [17].In yeast, it is known that the septins can act as a diffusion boundary during the polarization that accompanies budding [6].Here we use this model to analytically explore how differential diffusion rates could affect cell polarization.
Two sets of parameters are compared in Fig. 10.One is with uniform diffusion As expected, one observes in Fig. 10A that as diffusion increases, the polarization decreases.Another set of parameters is with a diffusion barrier in which we set D f = 0 and assume that the barrier is between the front and middle compartment that prevents the substance in the front compartment from diffusing to the middle compartment, but does not impede the transport from the middle to the front compartment.In Fig. 10B, it is observed that with D f = 0, as other diffusion rates increased, the polarization is not decreased but rather enhanced.The "unidirectional diffusion" from middle to front can be thought of as polarized transport, and according to our result, this unidirectional transport could serve as a mechanism for establishing polarization.

3.3.
Comparison between the two-compartment and three-compartment models.We close this section by comparing the two-compartment models with the three-compartment models, and summarizing some of the results described in Sections 2 and 3.
The two-compartment Models 1A, 1B, 2A, and 2B and the three-compartment Models 1A and 2B share the following common properties: k1 are critical parameters which dictate the number of steady states.The systems tend to have more steady states when c 1 , c 2 , c 3 are small.In other words, enhancing the level of gradient input or decreasing the positive feedback can reduce the number of steady states.
• As the diffusion increases, the number of steady states decreases.
• For the two-compartment and three-compartment Model 1A, stable steady states are found only when c 1 > 1, c 2 > 1 or c 3 > 1.For the two-compartment and three-compartment Model 2B, stable steady states are found only when c 1 , c 2 or c 3 are small (for example, less than 1).
The differences of those models are: • Generally, three-compartment models have more steady-state solutions than the two-compartment models.For example, the three-compartment Model 2B has up to 17 solutions while the two-compartment Model 2B has at most 5 solutions.• With a fixed number of compartments, models with positive feedback in the Hill form have more steady states and more stable steady states than those with the exponential positive feedback term.For example, the twocompartment Model 2A and 2B have up to 5 solutions while the twocompartment Model 1A and 1B can have at most 3 solutions.• For Model 1A, stable steady states are found only when c 1 > 1, c 2 > 1 or c 3 > 1, but for Model 2B, stable solutions are found only when c 1 , c 2 or c 3 are small (for example, less than 1).
4. A continuum model.In this section, we will extend our cell polarization model from the compartmentalized setting to a continuum spatial setting.A continuum model on the geometry of a cell membrane will be considered.In order to simplify the analysis, we assume that the cell membrane is a sphere embedded in a spatial gradient of ligand.With the symmetry of the geometric setup, we further assume that the distributions of the polarized membrane proteins are axisymmetric with respect to the axis aligned with the ligand gradient.Thus, the geometry of this problem can be simplified as a one-dimensional curve, from the back to the front of the cell, being parameterized by a parameter α.We denote the Cartesian coordinate of each point along this curve by (z(α), r(α)).While the cell is set to be of radius 1 µm, we choose the parameterization to be z = − cos α, r = sin α, with 0 ≤ α ≤ π.
We consider a model equipped with the same mechanisms included in the 2-and three-compartment models: This continuum model has been proposed and discussed in our previous work [3].In Eq. ( 37), the D m term is the lateral surface diffusion with a constant diffusion rate D m .This diffusion mechanism was implemented in the compartment models by the transport terms (D f , D b , D mb , D mf ).The k 0 parameter represents the cooperative production which depends on the input gradient u, where u = L mid + L slope z is a linear function of z.The form of this term is a Hill expression possessing a Hill cooperativity parameter q and a Hill half-maximal constant 1/δ.This cooperative production term corresponds to the P f , P m , P b terms in the compartment models.
The k 1 term is the positive feedback in which x 1 stimulates its own production.This autocatalytic reaction is also a cooperative reaction possessing a Hill cooperativity parameter h and a Hill half-maximal constant 1/γ.Note that there is a space dependent function p(α) in the positive feedback term.If p = 1, then the whole k 1 term is a regular positive feedback term.Here, we mostly consider p = 1 1+[δu(α)] −q 2 , which is a type of feedforward/feedback coincidence detection [23] in the positive feedback loop.As a result, the positive feedback term has a dependence on both x 1 and u, and the input-dependence is modulated by the cooperativity parameter q 2 in the Hill term.When q 1 = q 2 , p takes the same form as the k 0 term, and the implementation of the positive feedback loop will be the same as Model 2B in the compartment models.
In Eq. (38), s x 1 ds represents the surface integral of x 1 over the cell membrane, a unit sphere, and SA is the total surface area, which is 4π in our model.This global regulation mechanism was also implemented in the compartment models by the averaged x 1 term in the equations for x 2 .
Having the basic mechanisms in common, we would like to ask: What is in common to the continuum and compartment models?What are the differences?In the following, we perform analysis on the continuum model, followed by a comparison with the compartment models.

4.1.
Without lateral surface diffusion.We start the analysis of system (37)-(38) with D m = 0, i.e. there is no membrane diffusion.This happens when some membrane proteins after synthesis become anchored to the cytoskeleton and hence move locally but not globally, resulting in an effective macroscopic diffusion constant of D m = 0 [30].The steady-state equations of the system thus become a system of algebraic equations involving a global integral constraint: To simplify the notations, we let Eq. ( 39) then becomes Note that since Eq. ( 42) is a polynomial equation with the power h + 1, and hence, as h increases, the numbers of steady states will increase.With no general analytic solutions available, one way to solve the system (37)-( 38) is to solve x 1 as a function of α from Eq. ( 42), and then check if x 1 (α) satisfies the integral constraint Eq. (40).Since x 1 could have multiple roots at each α, there could be multiple x 1 (α) that globally satisfy the integral constraint.Due to the difficulties of describing the solutions analytically, in the rest of this section, we will solve the steady-state system numerically.
In the meanwhile, we would like to consider how to determine the local stability of a steady state, especially for our system which is not typical in the sense that it involves an integral constraint.First, we consider Eq. ( 37) with where w(x 1 , α) ) is a solution to Eqs. ( 39)-(40), taking ȳ as a constant, the solution is stable with respect to Eq. ( 43) alone if ) ≤ 0 for all α, and unstable if ∂ ∂x1 f (x 1 , ȳ, α) > 0 for some α.In the following proposition, we prove that the stability of Eq. ( 43) alone, taking ȳ as a fixed constant, is equivalent to the stability of the full system Eqs.(37)-(38).Therefore, to check the stability of a steady state, one only needs to consider the sign of ∂ ∂x1 f (x 1 , ȳ, α).Proposition 4. If (x 1 , ȳ) is a solution to Eqs. (39)-(40), and it is stable with respect to Eq. ( 43), this solution is stable with respect to the system Eqs.(37)-(38).
Proof.Using the notation in Eq. (43), Eqs.(37)-(38) can be written as: Let (x 1 (α), x2 ) be the steady-state solution.To analyze the stability of the steady states of Eqs.(37)-(38), we perturb the steady state by (e −λt φ 1 (α), e −λt φ 2 ): where φ 1 , φ 2 are negligible compared to x1 , x2 .Substitute x 1 and x 2 into equation (44) we get After linearization and using the fact that x1 and x2 are steady-state solutions, one obtains By substituting the Eq. ( 46) into Eq.( 45), we get Thus Integrate equation (47) on both sides over the surface ) < 0 and k 3 x1 > 0, λ has to be positive, otherwise the integrand would be negative for all α and the integral cannot be 1.
Having the simple criterion for local stability of the solution, we numerically analyze the steady states and their stability for h = 1 and h = 2. 4.1.1.Linear case: h=1.The steady-state system (39)-(40) becomes Eq. ( 48) has the roots: The "hat" is to emphasize that x1 satisfies Eq. (48) pointwise in space, but not necessarily a solution to Eq. ( 49).Since only the root with the positive sign makes x1 non-negative, there is at most one solution.
Note that in solving the roots of the polynomial (48), one views y as a constant, but since the integral constrain Eq. ( 49) also has to be satisfied, y cannot be arbitrary.It will be natural to ask how many y's there are such that solution of Eq. (48), x1 (α; y), satisfies the integral constraint?The following proposition ensures that there is at most one steady-state solution for system (48)-(49) in the case of h = 1.
Proposition 5.If h = 1, system (37)-(38) has at most one steady-state solution. Proof.Define where x1 (α; y) is defined as a solution of Eq. (48).To prove that there is at most one solution, it suffices to show that F (y) is monotonic.Here, we denote and it suffices to examine F (y): where J(α) = 2π sin(α) cos(α) 2 + sin(α) 2 is the Jacobian and is non-negative.
Note that here we assume x 1 is axisymmetric, and the geometry is a sphere.We further denote θ ≡ y − a(α), ω ≡ 2b(α)g(α)y and simplify the above equation to Because −a(α)( √ θ 2 + 2ω − θ) − ω < 0 and other terms are non-negative, we can conclude that F (y) < 0 and hence F (y) is monotonically decreasing.Therefore, there is at most one y satisfying F (y) = constant, in particular for our case F (y) = k ss .
Although h = 1 is a good model in the sense that it gives rise to a unique solution, the polarization is usually poor (Fig. 11A).When we examine how the polarization of the solution, measured by the solution at the front z = 1, changes with respect to the slopes of ligand concentration L slope , it is found that as L slope increases, the polarization increases, but peaks around L slope = 4 and slightly decreases for higher L slope .It can be seen from Fig. 11B that even with very steep slope, the polarization is still poor.4.1.2.Quadratic case: h=2.One normally expects better polarization with stronger positive feedback, which in our model corresponds to a larger h.Here we explore the case h = 2, for which the steady-state system (39)-(40) becomes In Eq. ( 51), for a given y and α, there are at most three roots, and Fig. 12 depicts how the real roots change as y changes.For a given y, we denote the set of stable real roots by A y (α), which may contain up to three elements for a fixed α.To look for solutions of Eqs.(51)-(52), we first consider a solution of Eq. (51), named x1 (α), which is a single-valued function in α.As in the h = 1 case, the notation "hat" here is to emphasize the fact that x1 (α) is a root of Eq. ( 51), therefore a subset of A y (α), but not necessarily globally satisfies the integral constraint.Since we are only interested in the polarized solutions, namely, solutions non-decreasing from the back z = −1 to the front z = 1, a natural way to form x1 (α) through the stable solution set A y (α) is: where α s is in the range [α min , α max ], in which more than one real stable roots exist.This choice of xs does not exhaust all the possible solutions, but will only pick the most polarized solution while keeping the minimal number of jumps in the function, which is a desirable property for our model.Fig. 13A shows that for a fixed y, two solutions xs corresponding to α s = α min (black) and α s = α max (green), which we later refer to as xmin and xmax , respectively.It can be seen from Fig. 13A that in between xmin and xmax , there are infinitely many functions xs satisfying Eq. ( 51) pointwise and has a jump between [α min , α max ].However, the actual solutions of the system (51)-(52) will then be determined by the integral constraint in Eq. ( 52) as explained in the following.
In order to evaluate the integral in Eq. ( 52), we define the quantity as in the case of h = 1, F xs (y) = s xs (α; y) ds SA .
Obviously, we have F xmin (y) ≤ F xs (y) ≤ F xmax (y).Hence, the solutions to the system (51)-( 52) exist if and only if F xmin (y) ≤ k ss ≤ F xmax (y).In order words, by evaluating F xmin (y) and F xmax (y), we can know if for the given y, there exists a solutions to the system (51)-(52).Moreover, for each given y, there is at most one solution x 1 , in the form of xs , to the system (51)-(52).To summarize, the steps to find the solution of system (51)-( 52) are (51)-( 52).
In Fig. 14A, we display solutions for a set of fixed parameters.Theoretically, there could be infinitely many solutions for this set of parameters, forming an "envelope of solutions" bounded by the red and green colored solutions.Only 11 solutions out of the envelope are shown in Fig. 14A.This "envelope of steady states" can be found for any h greater than 2. In Fig. 14B-D, one can observe how the solution envelope gets wider when h is increased.For a very large h, one can even find a solution polarizing at the wrong direction (magenta curve in Fig. 14D).This multiple steady state property of the model contributes to the difficulty of the system to track when the ligand gradient is reversed to the front-to-back direction, as mentioned in the introduction.The widening of the solution envelope as h becomes larger is also consistent with the argument in our previous work [3] that there is a tradeoff between tracking and amplification.The symbols with blue colors are stable roots while those with red color are unstable roots.
Next we examine systematically how the solution envelope changes with respect to different parameters or model variations.Typically, we use the value of x 1 at z = 1 to indicate the extent of polarization.Here, we additionally define an indicator called "polarization factor (P F )" to measure the extent polarization based on the width of the global distribution of the polarized component: , Ω is the surface of the sphere .
S p (x 1 ) is the surface area at the front of the cell that encompasses 50% of the polarized component x 1 and SA is the total surface area of the cell.An unpolarized cell would have a P F of 0 and an infinitely polarized cell would have a P F of 1.We concluded from our simulations that in most cases, both measures (P F and x 1 | z=1 ) conveyed the same information.In Fig. 15, the ligand slope L slope is plotted against the two measures of polarization, P F and x 1 | z=1 , for two variations of models.Fig. 15A corresponds to the model with a feedforward/feedback coincidence detection (p(α) = 1 1+[δu(α)] −q 2 ), while Fig. 15B corresponds to a normal positive feedback (p(α) = 1).Although our attention has always been on the feedforward/feedback coincidence detection model, we use this comparison to show the reason for favoring this model than the one with usual positive feedback.As analyzed above, the steady-state solutions generally forms an envelope, and in Fig. 15, we plotted the measures of the upper bound (red) and the lower bound (green) of the envelope.It is observed that the polarization of the feedforward/feedback model is much better than the normal feedback model.For both models, as the ligand slope increases, the envelope becomes narrower, and the trend is more obvious in the feedforward/feedback model with the envelope almost collapses to a single solution when the gradient is very steep.
We also investigate how the strength of the positive feedback, k 1 , affects the width of the envelope.In Fig. 16, k 1 is plotted against P F and x 1 at z = 1.It can be seen that when k 1 is small enough, the steady-state solution is unique, but as k 1 is increased, multiple solutions appear and the width of the envelope becomes larger.That implies with high gain of positive feedback, the range of steady states becomes wider, which will make the tracking harder.

4.2.
With lateral surface diffusion.As shown in the last section, when D m = 0, h ≥ 2 in the model (37)-(38), there may be an envelope of solutions which contains infinitely many solutions.It was also shown that the larger the feedback parameter h is, the wider the solution envelope is.In this section, we discuss how the envelope of solutions of model ( 37)-(38) may change when the diffusion D m is changed from zero to non-zero.4.2.1.Adding surface lateral diffusion decreases the number of solutions.As indicated by the analysis and numerical simulations for 2-and three-compartment models, the surface lateral diffusion improves the stability of the the solutions, and also reduces the number of steady states.Here we will investigate if the statement is still true for the continuum model.While the study of number of steady states is motivated by the concept of "tracking of directional change of the ligand gradient", we first compare how the tracking of the solution is different under different h and D m .In   observed that as D m = 0, for all h = 2, 4, 8, the tracking is not perfect, meaning that the tracked polarization is not of the symmetric shape as the original polarization.
However, the difference between those test cases is that for h = 2, 4, the solution gets to track to the reversed side, with the shape slightly different than the original polarization, while for h = 8, the solution does not track but rather is stuck at the front.According to the previous section, for h = 2, 4, 8, the steady-state solutions are within an envelope, but the width of the envelope increases as h increases, and h = 8 even has a reversely polarized solution.That explains why for h = 2 and h = 4, the solutions tracks when the gradient is flipped but not for h = 8. Surprisingly, when the diffusion coefficient is increased a little, for example D m = 0.001, the results for the same simulations are very different from the nodiffusion case (lower panel of Fig 17): the tracking is perfect for all h.Although these simulations do not exhaust all the possible initial conditions and parameters, they somehow shed light on the functions of surface lateral diffusion, in terms of "collapsing the envelope of solutions with D m = 0".
To examine the role of diffusion, we explore many different initial conditions, trying to find out as many steady states as possible.In our simulations, numerous initial conditions which are analytic solutions with D m = 0 are used.The choice of those initial condition is based on the expectation that by adding a small amount of diffusion, the solutions are likely to converge to functions close to one of the steady states with D m = 0.The same simulations and choice of initial conditions are performed for both h = 2, 4, 8.For h = 2, 4, we only found one steady-state solution for D m = 0.001, 0.01, 0.1, with polarization at the front (data not shown).In Fig. 18, results for h = 8 are shown.The initial conditions are displayed in Fig. 18A.For small diffusion D m = 0.001, 0.01, two steady states are found, but as the diffusion increases to D m = 0.1, only one steady state is found.
Although we did not exhaust all the initial conditions, this extensive exploration implies that for h = 2 and h = 4, the diffusion seems to "collapse" the solution envelope (D m = 0) to a single solution, and in the case of h = 8, the number of solutions is also greatly reduced, but multiple solutions may still exist.Table 1 summarizes the results, and it also contains the measure of polarization of each case.It can be seen from Table 1 that as h increases, the polarization is improved, but the number of steady states increases.This again supports our claim in [3] that there is a tradeoff between polarization and tracking.On the other hand, increased diffusion helps with picking a single solution, but at the cost of loss of polarization.• As the number of compartments increases, the number of steady states is dramatically increased.For example, the two-compartment models have up to 5 steady-state solutions, and three-compartment models have up to 17 steady states while the continuum model could have infinitely many solutions.• As the membrane lateral diffusion rate increases, the number of the steady states decreases.• As the rate constant of the positive feedback k 1 increases, the number of the steady states increases.
• When the gradient is steep, for example, in the compartment model when P f is large while keeping P b , P m fixed, and in the continuum model when the slope of the concentration is large, the number of steady state becomes small.

5.
Conclusions.In this paper, mathematical models for the polarization of proteins in a cell induced by external chemical gradients were developed and analyzed.Two different species -a polarized membrane bound protein, and a spatially homogeneous protein that is an inhibitor of that polarized protein -were considered, and amplification mechanisms such as cooperativity and positive feedback were included in the models.
We investigated these models in three different spatial settings: two-compartment, three-compartment and continuum space.The purpose was to study in a progressive fashion the number of steady states and their stability properties.It was found that as the number of compartments increased, the number of steady-state solutions also increased: in the two-compartment models, there are up to 5 steady states; in the three-compartment models, there are up to 17 steady states; and in the continuum models, one could have infinitely many steady states.
We also examined different types of models in which the functional form of the positive feedback was either exponential (Models 1A and 1B) or described by Hill terms (Models 2A and 2B).In the exponential models, the steady states tended to be less stable.
The compartmental models facilitated more detailed analysis including explicit relationships among the various parameters with respect to the steady states.These results were supported by numerical explorations over a range of parameter values.From this work, we found an interesting relationship between the number of stable steady states and the ratio of the cooperative production term to the positive feedback term.In general, stronger positive feedback relative to input-dependent cooperative amplification resulted in more stable steady states.
Interestingly, diffusion exerted a strong effect on the number of steady states in all of the models.This effect was most pronounced in the continuum models in which an envelope of infinitely many steady-state solutions in the absence of diffusion collapsed to a single solution by the presence of even low diffusion.
Finally, the simulations of the various models revealed that scenarios in which many stable steady-state solutions arose often correlated with the existence of one strongly polarized solution in the correct direction, but also with many solutions that were in the wrong direction.These incorrect solutions reduced the ability of the system to track a change in gradient direction because the simulated cell would get trapped in one of the misaligned steady states.
This analysis of the steady states and their stabilities in the different cell polarization scenarios highlights the challenges faced by the cell in its real environment with respect to performance tradeoffs.Importantly, this work also points out potential strategies for overcoming these challenges by modulating the ratio between positive feedback and cooperative production or by regulating diffusion.
In yeast cells, we believe that under wild-type conditions, the system may attempt to be monostable in order to track the direction of the gradient.However, there are circumstances in which wild-type cells will commit to a certain direction and ignore the gradient information by increasing the positive feedback resulting in the multiple steady-states depicted in the models in which the direction of the projection does not align with the direction of the gradient.In addition, mutant yeast cells may exhibit partial polarization or irregular polarization (Yi, data not shown) that is indicative of the more unusual steady-states observed in the modeling.
In the future, we plan to explore other challenges to robust cell polarization including broad dynamic range and noise attenuation.Ultimately, we are interested in how the complexity of the biological network and the regulation of spatial dynamics through positive and negative feedback helps to balance the various constraints that arise as the cell localizes its components in an asymmetric fashion according to uncertain internal or external cues.

b 1
Figure1.Schematic descriptions of two-compartment and continuous models of cell polarization.In both, the polarized species x 1 (red) becomes localized to the front of the cell through cooperative interactions (red arrow, P f ) in response to the input and through positive feedback (red +).There is a global negative feedback mediated by the species x 2 (blue).(A) Two-compartment model of cell polarization.The cell is divided into compartment 1 (front where the ligand concentration is higher) and compartment 2 (back).The species x 1 polarizes to the front (x f 1 ) and is less abundant at the back (x b 1 ).(B) Continuous model of cell polarization.The input gradient and spatial localization of x 1 is represented in a continuous fashion.

Figure 6 .
Figure 6.Number of steady-state solutions (x f 1 , x m 1 , x b 1 ) of system (14)-(17) with h = 2 under a range of parameters (c 1 , c 2 , c 3 ) ≡ ( P f k1 , Pm k1 , P b k1).Other parameters used are as follows:D mf = D f = D mb = D b = 0, k 2 = k 3 = k 4 = k ss = 1.The solutions are evaluated within the range 0.2 ≤ c 1 , c 2 , c 3 ≤ 2, with discretized space 0.2.Colors of red, magenta, yellow, green, cyan, blue and black stand for 1, 2, 3, 4, 5, 6, 7 real positive solutions respectively.Although Proposition 3 provides clues about the number of steady states when D = 0, due to the complexity of the system, one still relies on numerical simulations to obtain more details about how the number of steady states changes with respect to different parameters.In Fig.6, the number of steady states is evaluated within a range of parameters c 1 , c 2 , c 3 , with different colors indicating different number of solutions.It can be observed from Fig.6that without diffusion, when c 1 , c 2 , c 3 are all very small, there are more steady-state solutions (up to 7 steady states for the parameters we used); as one of the c j 's is increased, the number of steady states decreases.This implies that either increasing the cooperative production or decreasing the positive feedback can reduce the number of steady states.Next, we studied how the solution of x f 1 changes with parameters (P f , P b , k 1 and D f , D mf , D mb , D b ) in Fig.7, without assuming diffusion rates zero.According to Fig.7and other numerical simulations not shown here, we made the following observations: (1) Model 1A has 7 steady-state solutions only when c 1 , c 2 , c 3 are less than 1; (2) steady-state solutions are stable only when c 1 > 1; no more than 2

Figure 7 .
Figure 7. Steady-state solution of x f 1 versus P f for threecompartment Model 1A with h = 2.Other parameters used for these figures are k 1 = k 2 = k 3 = k 4 = k ss = 1.Red represents unstable steady-state solutions, and blue represents stable steady states.(A) D f = D mf = D mb = D b = 0 and P b = 0.1; the solutions marked by '.', 'o', ' ' correspond to P m = 0.2, 1, 1.8 respectively; (B) D f = D mf = D mb = D b = 0 and P m = 1.5; the solutions marked by '.', 'o', ' ' correspond to P b = 0.2, 1, 1.8 respectively; (C) D f = D mf = D mb = D b = 0 and P m = 1.5, P b = 1; the solutions marked by '.', 'o', ' ' correspond to k 1 = 0.1, 1, 10 respectively; (D) P m = 1.5, P b = 0.8, the solutions marked by '.', 'o', ' ' correspond to D f = D mf = D mb = D b = 0.001, 0.1, 1 respectively.Comparing three-compartment and two-compartment models with the exponential form of feedback (Model 1A), it is found from our numerical tests that they are mainly different in the total number of steady-state solutions and the number of stable steady states: (1) there are at most 7 steady-state solutions for the three-compartment model and at most 3 steady-state solutions for the twocompartment model; (2) there are at most 2 stable steady-state solutions for the

Figure 8 ., x m 1 , x b 1
Figure 8. Steady-state solutions of x f 1 , x m 1 , x b 1 of threecompartment Model 1A with h = 2. x f 1 , x m 1 , x b 1 are connected by straight lines to distinguish different sets of solutions.Red, blue, and yellow respectively represent the front-polarizing solution(x f 1 > x m 1 > x b 1 ), back-polarizing solution (x f 1 < x m 1 < x b 1 ), and the non-polarizing solution (xm 1 > x f 1 , x m 1 > x b 1 or x m 1 < x f 1 , x m 1 < x b 1 ); The parameters used are k 1 = k 2 = k 3 = k 4 = k ss =1, and P f = 0.3, P m = 0.2, P b = 0.1.(A) D f = D mf = D mb = D b = 0; (B) D f = D mf = D mb = D b = 0.1; (C) D f = D mf = D mb = D b = 0.5; (D) D f = D mf = D mb = D b = 1.

Fig. 13B
Fig. 13B displays F xmin (y) and F xmax (y) as functions of y.As shown in Fig. 13B, by drawing a horizontal line at F (y) = k ss , one can identify the interval for y, [y min , y max ], that admits polarized solutions satisfying both Eqs.(51)-(52).To summarize, the steps to find the solution of system (51)-(52) are

Figure 13 .
Figure 13.(A) All the real roots are plotted in circles, with blue circles representing the stable solutions and red for unstable ones; the black solid line is xmin and the green solid line is xmax ; (B) y versus F xmin (y) (green) and y versus F xmax (y) (blue); the horizontal gray dashed line is k ss (taken as 1 in this figure), and it crosses the two curves at y min and y max , respectively.[y min , y max ] is the interval that will give rise to solutions satisfying both Eqs.(51)-(52).The parameters used for the figures are k 0 = k 2 = k 3 = 1, k 1 = 10, q 1 = q 2 = 10, γ = 1, δ = 0.1, L mid = 10, L slope = 2 and k ss = 1.

Figure 14 .
Figure 14.(A) h = 2; 11 solutions out of infinitely many solutions in the envelope, bounded by the red and green colored solutions, are plotted.(B)-(D) The root curves displaying the steadystate solutions for increasing values of h; each curve represents the roots for a particular value of x 2 that satisfies the integral constraint; both stable roots (green circles) and unstable roots (red circles) are present.The highest polarized solution for each root curve is traced in blue.For h = 8, a reversed polarization solution is shown in magenta.The parameters used for the figures are k 0 = k 2 = k 3 = 1, k 1 = 10, q 1 = q 2 = 10, γ = 1, δ = 0.1, L mid = 10, L slope = 2 and k ss = 1.

4. 3 .Figure 17 .
Figure 17.A 180 • directional change is shown with and without diffusion for h = 2, 4, 8. Blue line represents initial polarization, red line represents reversed polarization, and dashed blue line represents the reversed polarization solution symmetric to initial polarization.Upper panel: D m = 0; lower panel D m = 0.001.(A) h = 2; (B) h = 4; (C) h = 8.The overlap between the dashed blue and red lines in the lower panel suggests perfect tracking.

Figure 18 .
Figure 18.Number of steady-state solutions for D m = 0.001, 0.01, 0.1 and h = 8. (A) 21 initial conditions for simulations in (B)-(D), and they are analytic steady states for D m = 0; (B) two steady-state solutions found for D m = 0.001; (C) two steady-state solutions found for D m = 0.01; (D) one steady-state solution found for D m = 0.1;

•
Scanning a wide range of y; for each y, find xmax and xmin .• Evaluating F xmax (y) and F xmin (y), and determine if k ss is in the interval [F xmin (y), F xmax (y)]; if it is, then y ∈ [y min , y max ]. • Identifying y min and y max with the previous step.• For each y ∈ [y min , y max ], one can get exactly one solution satisfying Eqs.

Table 1 .
Number of steady states Polarization (x 1 at z = 1) Number of steady states and the polarization with different D m and h.M: multiple solution; S: single solution.