Long time behavior and stable patterns in high-dimensional polarity models of asymmetric cell division

Asymmetric cell division is one of the fundamental processes to create cell diversity in the early stage of embryonic development. During this process, the polarity formation in the cell membrane has been considered as a key process by which the entire polarity formation in the cytosol is controlled, and it has been extensively studied in both experiments and mathematical models. Nonetheless, a mathematically rigorous analysis of the polarity formation in the asymmetric cell division has been little explored, particularly for bulk-surface models. In this article, we deal with polarity models proposed for describing the PAR polarity formation in the asymmetric cell division of a C. elegans embryo. Using a simpler but mathematically consistent model, we exhibit the long time behavior of the polarity formation of a bulk-surface cell. Moreover, we mathematically prove the existence of stable polarity solutions of the model equation in an arbitrary high-dimensional domain and analyse how the boundary position of polarity domain is determined. Our results propose that the existence and dynamics of the polarity in the asymmetric cell division can be understood universally in terms of basic mathematical structures.


Introduction
The polarity formation emerging during the early stage of embryonic development is a spectacular mechanism wherein a single mother cell of fertilised egg creates completely different daughter cells through asymmetric cell divisions (Campanale et al. 2017;Gönczy 2005). Because the polarity formation has initial and central roles during the entire process of asymmetric cell divisions, an elucidation of the polarity mechanism has been extensively explored in both experimental and theoretical approaches (see the recent review papers; Lang and Munro 2017;Rappel and Edelstein-Keshet 2017;Cortes et al. 2018). In particular, the PAR polarity formation in a C. elegans embryo has been well studied as one of a representative biological model. In the asymmetric cell division of C. elegans embryo, anterior PARs (aPAR) and posterior PARs (pPAR) are exclusively formed in an asymmetrical manner, which plays a key role in redistributing the cytosol substrates (Fig. 1a). PAR proteins are mostly upstream regulators that control the downstream proteins and a series of the processes of asymmetric cell division (Cuenca et al. 2002;Hoege and Hyman 2013;Wu et al. 2018). During the initial stage, aPAR and pPAR are homogeneously distributed in the membrane and cytosol, respectively. After the symmetry is broken by sperm entry, the pPAR spreads from the posterior pole, and the growth of the domain of pPAR polarity stops at approximately half the egg length with the formation of the two segregated polarity domains of aPAR and pPAR. This is consequently maintained for approximately 16 min, and a mother cell prepares for cell division (Gönczy 2005).
Owing to the key role of PAR polarity in asymmetric cell division, several mathematical models for the PAR polarity formation have been suggested and well-explored in recent years (Cortes et al. 2018;Goehring et al. 2011b;Marée et al. 2006;Seirin-Lee and Shibata 2015;Seirin-Lee 2016b, 2020Trong et al. 2014). Although the details of their model descriptions are different, a mathematical structure that generates a polarity pattern is taken into account universally, and the bi-stability structure in the kinetic terms and the property of mass conservation have been noted as the basal conditions. Nonetheless, few studies on the PAR polarity in terms of high-dimensional models reflecting the cell geometry of a bulk cytosol space and a surface cell membrane have been conducted numerically or mathematically. In particular, a rigorous mathematical analysis of the PAR polarity formation has been poorly explored, and most of the previous studies have been based on numerical observations. Although some mathematical studies have been conducted on PAR polarity models, they have dealt with one-dimensional problems or reduced systems with less variables (Kuwamura et al. 2018; Morita and Ogawa 2010;Morita 2012;Seirin-Lee et al. 2020). Moreover, the most fundamental question for the existence and stability of polarity solutions with two exclusive polarity domains has been infrequently studied owing to mathematical difficulties. Here, related to bulk cytosol-surface models, we refer to studies on Turing-type instability (Levine and Rappel 2005;Morita andSakamoto 2018, 2020;Röger 2012, 2014), the existence and stability of a polarized solution reduced on the sphere (Diegmiller et al. 2018), and polarized patterns numerically shown in 3dimensional domains with complex geometries (Cusseddu et al. 2019), where different types of reaction-diffusion systems are investigated. The readers who are interested in a basic mathematical theory may refer to Sharma and Morgan (2016) for the existence of time global solutions.
In this article, we consider the aPAR-pPAR polarity models suggested by Seirin-Lee and Shibata (2015) and Goehring et al. (2011b), and extend these models to include a high-dimensional case with a cell geometry composed of a bulk cytosol space and a cell membrane. The cell membrane consists of a lipid bilayer and contains membrane proteins which diffuse laterally through the membrane (Cooper 2000). Because a typical eukaryotic cell has a cell membrane thickness of 5-10 nm compared to a cell diameter of about 50 µm, the cell membrane can be considered as a surface domain with thin thickness (Fig. 1b(b1)). Then, the model system is described as follows: where is a bulk cytosol and ε = × D ε . (= ∂ ) is the surface of a bulk cytosol space and D ε is the thickness of cell membrane which is sufficiently small. In addition, P m and A m denote the concentrations of pPAR and aPAR proteins in the membrane, respectively, while P c and A c denote the concentrations of pPAR and aPAR proteins in the cytosol, respectively. stands for the Laplace operator defined in the bulk domain and the cell membrane domain ε . We do not consider the flux between cell membrane and extracellular space, so that zero flux boundary condition is assumed in the outer surface of cell membrane.
As |D ε | → 0, the model (1.1) can be directly reduced to a bulk-surface model ( Fig. 1b(b2)) described as follows: where is the Laplace-Beltrami operator on . The kinetic terms, F off and F off , imply the translocation of pPAR by aPAR and aPAR by pPAR from the membrane to the cytosol. In many polarity models, the offrate functions have been expressed in terms of nonlinear functions, allowing a bi-stable property. In this study, we consider two types of off-rate functions, which have been a b c Fig. 1 PARs polarity in asymmetric cell division and schematic diagram of model. a PAR polarity process in a C. elegans embryo is divided into three phases: symmetry breaking during the initial phase, the patterning(establishment) phase of an emerging pattern, and the maintenance phase of the stationary state of polarity. A, P indicate the points of anterior and posterior poles, respectively. b Schematic diagram of model reductions is shown. Gray-coloured regions imply cytoplasm and black-colored thick line region implies cell membrane. D ε and D ε are the cross-sectional areas of membrane and cytoplasmic spaces which are separated by the inner boundary region ( ) between cell membrane and cytosol. indicates bulk cytosol space in R N , and ε is the cytosol region around the cell membrane region, ε . c The domain shapes of = (0, L) × D are shown with respect to Neumann and periodic boundary conditions. D is defined by the edge points of line, the vertical line, and cross-sectional area in one, two, and three dimensional spaces, respectively (color figure online) explored as a Hill function type in Seirin-Lee and Shibata (2015) using and have been studied in Goehring et al. (2011b) through where α and α are the basal off-rates of pPAR and aPAR, respectively, and K , K 1 , K , K 1 , K 2 , and K 2 are positive constants determining the off-rate effects. In addition, γ and γ are on-rates parameters. Note that the original models for PAR polarity include advection terms (flow effects) because an acto-myosin contraction causes an advective transport (flow) in both the membrane and cytosol after the symmetry breaking (Goehring et al. 2011b;Niwayama et al. 2011). However, we neglect the advection terms in our model because the flow effect is ceased in the maintenance phase (approximately, 13 minutes later after symmetry breaking) where the polarity domains stop spreading (Munro and Nance 2004;Goehring et al. 2011b). In addition, the purpose of this study is to focus on the long time behavior of the polarity solutions.
In this paper, we explore the high-dimensional polarity models, (1.1) and (1.2), and confirm the existence of stable polarity solutions using numerical simulations. We then prove the existence and stability of such solutions with respect to a highdimensional case through a simplification of the model (1.1), called a cell membrane periphery model, and a model reduction to a shadow system of the bulk-surface model (1.2). In our study, a rigorous proof for the existence of stable polarity solutions in the aPAR and pPAR polarity models within an arbitrary high-dimensional domain is proposed. Based on our analysis, we further explore how the boundary position of the polarity domains is determined. Our results suggest that the existence and dynamics of the PAR polarity during asymmetric cell division can be understood based on a basic mathematical structure, which should be held universally without dependence on a specific choice of parameter values.
We remark that our mathematical analysis is carried out by making use of the energy functional under the assumption of a large difference of diffusion coefficients of the proteins in the membrane and the cytosol. We refer to Mori et al. (2011) for another approach to identify the polarity boundary in a simpler model called the wave pinning model of two-component reaction-diffusion equations. In the wave pinning model, the authors succeeded in obtaining a reasonable approximation of polarity boundaries in one-dimensional domain by using asymptotic analysis with the membrane diffusion coefficient as a small parameter of the membrane diffusion coefficient. In Cusseddu et al. (2019), a similar analysis is performed on the bulk-surface version of the model in Mori et al. (2011).

Nonlinearity of off-rate functions
We here show that the model of (1.4) suggested by Goehring et al. (2011b) is close to the model of (1.3) suggested by Seirin-Lee and Shibata (2015). By using an approximation through a Taylor expansion around small concentration of A m and small concentration of P m for the off-rate functions (1.3), respectively, we have (2.1) This approximation implies that the model of (1.4) can be considered as a special case of the model of (1.3) by taking K 2 = K 1 /K and K 2 = K 1 /K . Furthermore, the two models show little differences in the qualitative dynamics in the long time behavior and it may be sufficient to choose either of them in order to understand the mathematical structure of a stable polarity pattern (see Sect. 2.2). Therefore, we choose the nonlinear off-rate functions by (2.1) when we perform mathematical investigation.

A cell membrane periphery model
The diffusion coefficient of protein in the cytosol is generally much larger than that in the membrane (Goehring et al. 2011b;Kuhn et al. 2011). We thus assume that the fast diffusions of aPAR and pPAR in the cytosol lead to a well-mixed state and that the concentrations of aPAR and pPAR in the cytosol quickly approach to a uniform state in the region far from the membrane. Then, as for the polarity patterns in the model systems (1.1) with (1.3) or (1.4), it is reasonable to consider the peripheral region of the membrane (Fig. 1b(b4)). Thus, the model (1.1) can be directly considered in a cell membrane periphery domain by defining the cytoplasmic domain with ε ≡ × D ε instead of , where D ε is the sufficiently small cross-sectional area of cytosol region in the periphery of the cell membrane. Note that with the assumption of |D ε | = |D ε | 1, we may consider the model (1.1) on a domain where the cell membrane and cytoplasm are overlapped. That is, where |D| = |D ε | = |D ε |. Therefore, we can simplify the model (1.1) to the above model defined on a high-dimensional space under the Neumann boundary conditions or periodic boundary conditions. With setting in (2.1), and in terms of the new variables and parameters defined by we rewrite the system (2.2) as follows: in a bounded domain (⊂ R N ) with the Neumann boundary condition or periodic boundary conditions. The condition (2.7) is due to the fact that the diffusion coefficient of protein in the cytosol is generally larger than that in the membrane (Goehring et al. 2011a, b;Kuhn et al. 2011). is given to × D = (0, L) × D where D is a bounded domain of R N −1 . Figure 1b(b4)-(b6) and 1C show the cases of N = 1, 2, 3. L represents the perimeter length of the cell and |D| represents the thickness or crosssectional area (which does not have to be a disk) of the cell membrane or cytoplasmic domain in two and three dimensional spaces, respectively. We call the system (2.3)-(2.6) a cell membrane periphery model. The dynamics of PAR polarity in a C. elegans embryo can be understood as the case of a cell membrane periphery model with periodic boundary conditions.

A shadow system of bulk-surface model
We introduce a shadow system of the bulk-surface model (1.2) in this subsection. As seen in the later Sect. 3.2, we reduce the stationary equations of the cell periphery model (2.3)-(2.6) to a simpler system. Then, we show that the reduced equations can be obtained as the Euler-Lagrange equations of some energy functional and, interestingly, the shadow system of the bulk-surface model (1.2) with (1.4) has a similar form to the gradient flow of the functional. Thus, the mathematical results for the cell membrane periphery model (2.3)-(2.6) can be directly applied to the shadow system of the bulk-surface model (1.2) with (1.4) ( Fig. 1b(b3), (b7))(See Section 3.4). Let us return to a bulk-surface diffusion model (1.2). Namely, we consider the equation for P c and A c in a bulk domain , where the mass transport takes place on the boundary = ∂ . The bulk-surface diffusion with mass transport on the boundary in a non-dimensional form is given by This system has the following conservation of mass: We set We reduce the equations to a shadow system on by taking D c , D c → ∞ (Fig. 1b(b2),(b3)) as and the system turns to be and (2.9) Then we have equations for ξ and η aṡ which are derived by taking the spatial average of the bulk equations for P c and A c and applying the Green formula. However, these equations for ξ and η are obtained by differentiating (2.9) with respect to t and using (2.8). Hence, it suffices to handle (2.8) under the constraint (2.9) as the shadow system. Introducing the new variables for the model (2.8) with the off-rate function (1.3), we reduce the parameters as and let us set (2.10) Dropping of t in the equations for u and v leads to In a similar manner to the above derivation, we can also derive the shadow system of the cell membrane periphery model (2.3)-(2.6) in a bounded domain ⊂ R N as (2.14) In the sequel, we are concerned with the shadow system of (2.3)-(2.6) in , Note that the approximation of shadow system induces a model defined on a surface domain which does not include the thickness of cell membrane. However, we can deal with the shadow system mathematically similarly to a cell membrane periphery model assumed sufficiently small D as shown in Fig. 1b. As the simplest case, we can handle the above system in = (0, L) with periodic boundary condition. Although (2.17)-(2.18) with (2.19) are simplified model equations, it is very useful to perform numerical simulations with low numerical costs in higher-dimensional domains. We will numerically show that the solution of this shadow system exhibits similar behaviors as the other models in the following section.

Numerical simulation results
In this study, we consider the long time behavior of segregation pattern of P m and A m in the system, namely, a stable steady-state polarity patterns of pPAR and aPAR. Thus, we first compare the quantitative dynamics of the two bulk-surface models (1.2) with (1.3) and (1.4) for the long time behavior (see Appendix A for the numerical method and scheme for solving bulk-surface model). Figure 2a shows that both of models successfully generate exclusive polarity patterns and that there are no notable differences between the two models. Furthermore, as shown in Fig. 2b of the profile of the membrane PAR proteins, the two polarity domains are overlapped with a small interface region. This may allow us to assume that the off-rate effect of P m by A m is negligible in the region where the concentration of P m is low (i.e. aPAR dominated region) while the off-rate effect is likely to be dominated in the region where the concentration of P m is high (i.e. pPAR dominated region). Similarly, the off-rate effect of A m by P m is likely to be dominated in the region where the concentration of A m is high (i.e. aPAR dominated region). The approximation by a Taylor expansion in (2.1) may imply this phenomenon. Figure 2c shows stable polarity solutions of the cell membrane periphery model (2.3)-(2.6) on (⊂ R 2 ) under the Neumann boundary conditions. Fig. 2d shows nonconstant equilibrium solutions of shadow system (2.17)-(2.18) on (⊂ R 1 ) for both Neumann and periodic boundary conditions. The simulation results show, at least apparently, that there are no differences in the qualitative dynamics of polarity solutions between the bulk-surface model, the cell membrane periphery model, and the shadow system.
Note that the kinetic parameter values used in the representative simulations were chosen as arbitrary values which basically satisfy the condition that two stable equilibria exist in kinetic equations. The diffusion coefficients were chosen as the values of biologically feasible scale based on the data of C. elegans embryo (Goehring et al. 2011a;Kuhn et al. 2011;Seirin-Lee 2021). Although we showed the polarity solutions with representative parameter values, our rigorous analysis and mathematical results of this paper support that the existence of polarity pattern is less sensitive to parameter choice and holds over a wide range of the parameter space.

Basic properties of solutions to the system
Let be a bounded domain in R N with smooth boundary ∂ . In this section, we consider the system (2.3)-(2.6) in with the Neumann boundary conditions Assume nonnegative continuous initial data where u i,0 , v i,0 (i = 1, 2) are not identically zero. We aim to show the fundamental mathematical results on the positivity and global boundedness of the solution to the cell periphery model (2.3)-(2.6) with (3.1) and (3.2). The results here is easily interpreted to the case of periodic boundary condition. The first lemma assures the positivity of the solution.
for t > 0 in the maximal interval.
Next, in order to prove that the solution can be extended globally in time, we introduce the new variables z where u = (u 1 , u 2 ), z = (z 1 , z 2 ), and θ : (i = 1, 2) are given in (3.9). Then, we can prove that E plays as a Lyapunov function of (3.3)-(3.6) and obtain the next result.

Lemma 3.2 Given positive parameters
, there exists a constant C 1 , depending on the initial data, such that Using this lemma, we get to the global boundedness of the solution; The proofs of the above three lemmas are given in Appendix B.1.

Remark 3.1
In the work of Latos et al. (2018), a similar argument for the proof of the uniform boundedness can be found. They study a two-component system having the nonlinear terms with linear growth at infinity by which the uniform boundedness can be shown without the restriction of space dimensions. Another approach is also found in Jimbo and Morita (2013) for the similar system to Latos et al. (2018).

Stability of equilibrium solutions
From the system (2.3)-(2.6) with (3.1), we see so that the system allows mass conservations We set This mass conservation is expressed in the system (3.3)-(3.6) as In this subsection, we investigate the stationary problem of the system (2.3)-(2.6), that is, of the system (3.3)-(3.6), These equations turn out to be We put and let θ i (i = 1, 2) be same as in Lemma 3.2. The system (3.10) has a variational structure. As a matter of fact, consider the following functional of u = (u 1 , u 2 ). (3.12) Then, it is easy to see that (3.10) with the Neumann boundary conditions is the Euler-Lagrange equation of E s . Moreover, the following time evolution system serves as the gradient flow of the energy functional E s . We note that corresponding to a solution (u * 1 , u * 2 ) of (3.10), provides an equilibrium solution to (3.3)-(3.6). In the rest of this subsection, we show that the stability of (u * 1 , z * 1 , u * 2 , z * 2 ) is closely related to that of (u * 1 , u * 2 ) in (3.10). We first observe some property of a local minimizer of E s .

Lemma 3.4 Let
Then, there exists an ε 1 > 0 such that for any ε ∈ (0, ε 1 /4], we can take Proof By a slight modification of the proof of Lemma 7 of Latos and Suzuki (2014), we can get to the assertion of the lemma. Thus, we omit the details.
We decompose where θ i (i = 1, 2) are as in (3.9). From the next lemma, we see that the stability of a solution (u * , z * ) is inherited from that of u * .
The proof of this lemma is given in Appendix B.2.

Existence of stable nonconstant solutions
In Sect. 2.2, we confirmed numerically that there exist stable polarity solutions in the cell periphery model (2.3)-(2.6) in Fig. 2c. One can also find that the stable equilibrium solutions in the system (2.3)-(2.6) and the transformed system (3.13) are very consistent, as shown in Fig. 3a. As seen below, the system (3.13) of two variables gives a good mathematical simplification for rigorous analysis, while it is very useful to understand the dynamics of polarity patterns to a stationary state with reducing a numerical cost.
In this subsection, we prove the existence of stable nonconstant equilibrium solutions of (2.3)-(2.6) with (3.1) in a parameter regime. To this end, in view of Lemma 3.4, it suffices to show the existence of a nonconstant minimizer of the energy functional E s of (3.12). As the first step, we ensures the existence of a positive minimizer of E s .

Lemma 3.6 There is a minimizer
Proof By the direct method of calculus of variations, we have a minimizer u * = (u * 1 , u * 2 ) of E s . We easily exclude the case u * , a contradiction. Similarly, u * 2 ≤ 0, u * 2 ≡ 0 cannot happen. This implies that each u * i has a positive maximum or u i ≡ 0. Since (u * 1 , u * 2 ) satisfies (3.10), u i ≡ 0 (i = 1, 2) is excluded. By the Hopf lemma, we see that the maximum point of u * i exsits in the interior of the domain unless u * i is constant. Moreover, hold by applying the maximum principle to the first and second equations of (3.10) (use a contradiction argument for the proof). We exclude the case that u * i takes the non-positive minimum. If it happens, the minimum is achieved in the interior of by the Hopf lemma. Let x m ∈ be such a point. Then which is a contradiction. In conclusion, u * i > 0 (i = 1, 2). We next consider a positive constant equilibrium solution, which is obtained by a solution of (3.16) Given solution (ξ, η) to the system (3.16), we have where β i and θ i (i = 1, 2) are defined in (3.11) and (3.9), respectively. We easily see that for fixed m i (i = 1, 2) and k, it is impossible to realize the three conditions simultaneously. Hence, if we found a family of nonnegative functions parametrized by, say ε > 0, as E s allows a nonconstant minimizer for sufficiently small ε. Moreover, applying Lemma 3.5, we can assert the existence of a stable nonconstant positive equilibrium of the 4component system of (2.3)-(2.6). However, we are interested in the profile of the nonconstant minimizer. As a matter of fact, in numerical simulations, the spatial segregation pattern of u 1 and u 2 can be observed. Considering those simulations, we expect that the interface separating two regions of u 1 and u 2 emerges in the singular limit by varying appropriate parameters. Here we aim to construct an approximate solution with less energy than any constant solution and identify the location of the interface by minimizing the energy of the approximate solution. Our goal is to provide a reasonable approximation matching with the results of numerics.
Let us set where serves as a location of an interface in the approximate solution. We consider the parameter regime and use the following test functions: We compute the energy E s (φ, ψ). By a lengthy but straightforward computation (for the details, see Appendix C), we obtain Note that by (C.2), Put E s ( ) := E s (φ, ψ). In view of (3.18) and (3.20), we clarify the terms with up to the system of (2.3)-(2.6) in with (3.1) possesses a stable nonconstant equilibrium solution.
Consequently, we obtain This estimate might be useful for a mathematically rigorous study of the limiting behavior in the future work. Note that for the minimizer u We also expect that (φ, ψ) with = b gives a reasonable approximation of the minimizer of E s in the parameter regime stated in Proposition 3.7, though (φ, ψ) is not smooth at = b . As a matter of fact, the numerical test confirms that E s has a minimal energy at b = 0.750762 (Fig. 3b). Note that b is corresponding to the location of the interface edge of pPAR domain boundary. When we define the edge of pPAR interface by the region where pPAR concentration is included in (0.1%, 0.5%) of the minimal value of pPAR concentration in the full system (2.3)-(2.6) (Figs. 2c, 3a), the range of edge region of pPAR domain is approximated as a value in (0.745, 0.765), indicating that b = 0.750762 gives a good approximation with respect to the location of polarity boundary. Fig. 3c also confirms that the approximate solution (3.21) for = b well captures the characteristic of the shape of polarity solution of the full system (2.3)-(2.6) shown in Fig. 3a. We will leave a mathematically rigorous proof of the minimizer of E s in a future work because it is beyond the scope of this paper.

Bulk-surface model and its shadow system
In this subsection, we show the existence of nonconstant equilibrium solutions for the shadow system (2.17)-(2.18).
As in §3.2, we define the following energy functional: (3.29) Then we can directly prove the existence of nonconstant equilibrium solution for (2.17)-(2.18) with (2.19), provided that the system has the scaled parameters (2.10) with the condition in Proposition 3.7. Numerical simulations also confirm that nonconstant equilibrium solutions exist for both Neumann and periodic boundary conditions (Fig. 2d). As seen in the derivation of the shadow system from the bulk-surface model, it is reasonable to consider the shadow system in a 2-dimensional closed surface of the cell membrane. It is also interesting to study how the geometry of the surface affects the spatial pattern. Those issues are beyond the present scope and will be good as a future work. However, we expect that the above observation in a general domain will be useful for a further study in various domains. Goehring et al. (2011b), the following model equations for plolarization of PAR proteins are proposed:

Remark 3.2 In
where A and P stand for the local membrane concentrations of aPARs and pPARs, respectively. On the other hand, A cyto and P cyto stand for the uniform cytoplasmic concentrations of aPARs and pPARs, respectively, that is, those are assumed to be constants. In this model, they set a bistable character of the system, that is, the system allows two stable spatially uniform steady states by tuning parameters suitably. Goehring et al. (2011b) shows that the cortical flow velocity, ν = ν(x, t), plays a role in serving as a mechanical trigger for the pattern formation in the initial stage and the bistable character plays to create a segregation pattern of A and P after the flow effect ceases. However, when ν = 0, there is no stable equilibrium state because the system with ν = 0 belongs to a class of competition-diffusion systems. In fact, Kishimoto (1981) and Kishimoto and Weinberger (1985) show the instability of nonconstant equilibrium solutions of two component reaction-diffusion equations of competition type in an interval or a convex domain with the Neumann boundary condition or periodic boundary condition. On the other hand, by virtue of the nonlocal constraint of (2.16), the system (2.17)-(2.18) allows a stable nonconstant equilibrium in a parameter regime as seen in §3.3. Although the mathematical theorem ensures non-existence of stable segregation patterns in the system (3.30)-(3.31) with ν = 0, a very slow motion of the transition layer might happen in Goehring et al. (2011b), as seen in the scalar reaction-diffusion equation by Car and Pego (1989) and Fusco and Hale (1989).

Energy functional and biological implication
The length of the PAR polarity domain in asymmetric cell division is considered an important factor in regulating the location of the cell cleavage plane during asymmetric cell division (Coffmana et al. 2016;Morton et al. 2002;Rose and Kemphues 1998). However, the mechanism by which the polarity length is determined remains elusive, although this issue has been partially noted in both experimental and theoretical studies (Dawes and Iron 2013;Goehring et al. 2011b; Seirin-Lee and Shibata 2015; Seirin-Lee 2020). Herein, we explore how the size of polarity domains (namely, the position of We first explore how the boundary location of the pPAR polarity (namely, b ) is determined, which also defines the relative size of a polarity domain with respect to the size of the cell membrane. The boundary location of the aPAR polarity is be given to L − b . In the model systems (2.3)-(2.6) for = (0, L) × D, L can be regarded as the cell circumference of the cell membrane and |D| is the membrane thickness. We thus fix the size of |D| and vary L. We calculated b with respect to L D = L/|D| by using the energy functionals, (3.24), (3.25), and (3.27), the results of which are shown in Fig. 4a. We first found that when the size of the cell membrane is relatively small (i.e. in a small cell), the size of the pPAR domain decreases as the length of the cell membrane increases. This suggests that the length of the polarity can be critically affected by the shape of the cell (Fig. 4a, b). Because the circumference of the cell membrane in an elliptic cell is larger than that in a circular cell under the same volume, a circular cell has a longer domain of pPAR than the elliptic cell. That is, we  (Fig. 4b). By contrast, when the volume of the cell is large, namely the length of the cell membrane is large, the relative domain size of pPAR becomes constant without a dependence on L D (Fig. 4a, c). This suggests that the boundary position of the polarity is determined very robustly without the sensitivity of the cell shape. This fact is likely to be held for any cell shape, such as in the case of wrinkled cells (Fig. 4b, c).
Next, we investigate the effect of the total mass. We fix the total mass of aPAR (m 2 ) and vary the total mass of pPAR (m 1 ), which consequently leads to a variation of the total masses of pPAR and aPAR (m tot = m 1 + m 2 ) because m 1 /m 2 = (m tot /m 2 − 1). As m 1 (or m tot ) increases, the length of the polarity domain is increased (Fig. 4d). Unexpectedly, however, we found that once the total mass becomes larger than a certain amount, the length of the polarity domain is constant despite the increase in the total mass. This implies that if the total mass is sufficiently large, it may not be a dominant factor in determining the polarity length, while the other kinetic parameters or elements are likely to play an important role.
Finally, we explored both the effects of the length of the cell membrane (L D ) and the total mass of the polarity proteins on the length of the polarity domain (Fig. 4e). We found that a proper relation between the total mass and the cell membrane length is required to obtain a proper size of the polarity domain (a proper position of the polarity boundary) (Zone A). Specifically, if the length of the cell membrane is too small or the total mass is insufficient, the polarity domain either spreads out to a wide region of the cell membrane or remains in an area with a small polarity (Zone C). We also found that a wide parameter region exists where the sensitivity of L D and the ratio of the total mass to the size of the polarity domain becomes negligible (Zone B). Note that, although the reduced energy functional e in (3.27) does not show a good approximation quantitatively compared to the other two energy functionals E s in (3.24) and E s in (3.25) (Fig. 4a, d), it still captures well a qualitative characteristic of the nonconstant stationary solution. Thus, we can directly obtain similar conclusions with the Eq. (3.28).

Discussion
The pattern formation system characterised by bi-stability and mass conservation has been highlighted using a cell polarity model and explored based on several kinetic types of models (Goehring et al. 2011b;Mori et al. 2008;Otsuji et al. 2007; Seirin-Lee and Shibata 2015; Seirin-Lee 2021; Trong et al. 2014). In this study, we have explored the long time behavior of stable patterns for high-dimensional polarity models describing the PAR polarity occurring in asymmetric cell division, which are based on the models suggested by Goehring et al. (2011b) and Seirin-Lee and Shibata (2015). We showed that the long time behavior around the stationary solution in the bulk-surface model can be understood by reduced models; a cell membrane periphery model and shadow system have a mathematically simpler form but are able to capture the dynamics of the stationary polarity solutions in a high-dimensional space, and we rigorously established the existence and stability of the polarity solutions.
In this study, we also succeeded in constructing a reasonable approximate solution that provides a direct calculation of the energy functionals. Using these energy functionals, we have found the detailed effect of the cell membrane length and total mass on determining the size of polarity domain (the position of the polarity boundary). From this analysis, we also found that the size of the polarity domain can be sensitively altered depending on the cell shape when the cell size is small. By contrast, the relative location of the boundary is not changed by the cell shape when the cell size is large, which implies that a bigger cell may need to have fewer factors to maintain a robust length of the polarity domain than a smaller cell.
A previous study by Seirin-Lee and Shibata (2015) has found that the pPAR polarity length is linearly dependent on the total mass. However, their observations were restricted to a specific parameter range because the total mass was not an independent parameter but was determined based on the choice of kinetic parameters. A numerical analysis carried out on full partial-differential equations with respect to multiple parameter values often limits the computations in terms of the cost of the numerical calculation, though the recent study on the numerical method and scheme shows the development (Hao and Xue 2020;Uecker et al. 2014). By contrast, we can deal with the total mass as an independent parameter in our analysis using the energy functionals. This provides a highly effective numerical cost when we investigate the effects of the parameters on the dynamics of the solutions with a polarity profile.
In this study, our main approach to the model reductions has been based on the assumption of homogeneous state of cytosol protein which is basically resulted from a fast diffusion in the cytosol space. The approximation of shadow system requires a mathematical condition of diffusion infinity in the cytosol (D c → ∞), but it is not a biologically feasible condition. On the other hand, the cell membrane periphery model does not require such a limit condition and we can obtain mathematically rigorous results under the biologically reasonable condition that the diffusion in the cytosol is faster than that on the membrane (D c > D m ). This is generally observed in a cell (Kuhn et al. 2011). Thus, our analytical approach to the shadow system via a cell membrane periphery model supports that the mathematical assumption of fast diffusion in the cytosol is reasonable to understand the essential structure of polarity formation.
Although our mathematical analysis captures well the dynamics of the polarity in qualitative terms, it would be interesting to confirm the effect of the cell geometry on the length of the polarity domain in real biological systems. We are currently considering this, and leave it as a subject of future study. and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

A.1 Bulk-surface modeling using phase-field method
To describe the bulk-cytosol space and the membrane surface, we used a method of combining a phase-field function with the reaction-diffusion system. We apply the method proposed by Seirin-Lee (2016a, 2017 with respect to the case of a single cell. We solved the following equation in order to generate a cell shape: where α(> 0) is the intensity constant of the energy for the volume V (t), and V is the target volume of the cell. With a phase-field function satisfying (A.1) for some time t * , φ(x, t * ), we define a fixed cell domain as follows : In next, we combine the bulk-surface model (1.2) with the cell phase-field function φ(x, t * ) (Kockelkoren et al. 2003;Levine and Rappel 2005;Wang et al. 2017). The model system combined with the phase-field function which we used for the numerical simulations of bulk-surface model (1.2) is given to where B(φ) = νφ 2 (1 − φ) 2 (ν > 0), a function defining the membrane region. Because the cell shape can be easily generated by the Eq. (A.1) by adjusting the target volume V and initial condition, the phase-field-combinded modeling is very

A.2 Numerical scheme and stability
We solved the form of the bulk-surface system (A.2) using C language by a standard explicit numerical scheme on square grids. The system domain is given to [0, 1]×[0, 1] and has been spaced by 400 × 400 grids which defines x = 0.0025. t = 1/300 is chosen so that the stability condition of explicit scheme is sufficiently satisfied (Morton and Mayers 1994). We also have confirmed that the simulation results are not changed with smaller spatial and temporal grid size. Figure S5 confirms that the system (A.2) combined with phase-field functions is wellconserved numerically with respect to the total mass. On the other hand, we solved one and two dimensional reaction-diffusion systems with Neumann and periodic boundary conditions by implicit method using LU decomposition and ADI method (Morton and Mayers 1994).

A.3 Initial conditions
The initial condition for simulating the bulk-surface model (1.2) was given to locally concentrated initial conditions (Seirin-Lee and Shibata 2015) by where σ is the strength of the signal, and ω is the sufficiently small region in which the signal is imposed. δ is a positive constant and has a very small value. Typically, we set σ = 40 with |ω| = B(φ)(0.05) 2 and δ = 0.001. ψ(X) is a random function of uniform distribution, which takes values in the range of [−0.5, 0.5]. The detailed values used in the representative simulations are A 0 m = P 0 m = 1.02,

Proof of Lemma 3.1 We first consider
in with the same initial condition, where ( u(x, ·)) + := max{ u(x, ·), 0}. This system has a unique classical solution in ×[0, T ) for some T > 0 [see Rothe (1984)]. Since we use the maximum principle to assert v 1 (x, t), v 2 (x, t) ≥ 0 (0 ≤ t < T ). This leads us u 1 (x, t), u 2 (x, t) ≥ 0 (0 ≤ t < T ) by applying the maximum principle to the u 1 and u 2 equations. The strong maximum principle shows the strict positivity of the solutions for t > 0. In the sequence the solutions, (u 1 , u 2 , v 1 , v 2 ) and ( u 1 , u 2 , v 1 , v 2 ) coinside in × [0, T ) because of the uniqueness of the solution.
Before going to the proof of Lemma 3.3, we prepare some estimate for the heat equation. Let u(·, t; u 0 ) be a unique solution of ∂ t u = d u in , ∂u/∂n = 0 on ∂ , u(·, 0) = u 0 (∈ C( ; R)), and define the the semigroup {e td } t≥0 by e td u 0 := u(·, t; u 0 ). We will use the following inequality ( Rothe (1984)):

Proof of Lemma 3.3
In view of u i (·, t) + v i (·, t) = m i and the positivity of the solution, we see Then, there is a positive C 2 such that for i = 1, 2 holds.
In the rest of the proof, we omit to state the boundary condition of the equations since it is fixed. We prove the assertion of the lemma for only τ = 1 below since the same argument can be easily applied to the case τ = 1 with a slight modification.
When N = 1, since H 1 ( ) ⊂ C 1/2 ( ), we obtain the uniform boundedness for bounds v i (·, t) by the maximum principle. Hence, we easily see the uniform boundedness for v i (·, t) for i = 1, 2. For N ≥ 2, taking q = 2 and r with 1 ≤ r < 2N /(N − 2) + , we have and by (B.1) (precisely, C 0 = C 0 (d, 2, r ), but we omit the dependence on d and r ). The equations of (3.4) and (3.6) allow the integral form as for i = 1, 2. In view of (B.2), (3.8) and the uniform boundedness of the initial data, there is a positive constant C 3 such that Since the solution u i is bounded by a solution u i to with the same initial condition, we estimate the solution to (B.5). We write (B.5) in the integral form as We notice that the inequality N 2 1 r < 1 2 leads to N < r < 2N /(N − 2) + , so (N − 2) + < 2. Thus, for N = 3, by taking q = 5 and r = ∞ in (B.1), we have Applying this to (B.6) yields that there is C 5 > 0 such that In terms of (B.4), the uniform boundedness of u i (·, t) L ∞ follows.
In order to obtain the uniform boundedness of v i , we consider (B.3) and make use of the maximum principle again.
As for N = 2, we can take any q greater than 2 and r = ∞ in (B.1). A similar argument leads to the desired assertion.
In the sequence, we obtain Combining this inequality with (B.11), we obtain the desired assertion.