Balancing reaction-diffusion network for cell polarization pattern with stability and asymmetry

eLife Assessment

This manuscript makes important contributions to our understanding of cell polarization dynamics by demonstrating how compensatory regulatory and spatial mechanisms enhance the robustness of polarization patterns. By integrating a computational pipeline with comparisons to experimental data, the authors provide convincing evidence that stability and asymmetry in reaction-diffusion networks are crucial for polarization in C. elegans zygotes. Their findings offer novel insights into essential biological processes such as cell migration, division, and symmetry breaking. Future theoretical and experimental work could refine the model by addressing its acknowledged limitations.

https://doi.org/10.7554/eLife.96421.3.sa0

Significance of the findings:

Important: Findings that have theoretical or practical implications beyond a single subfield

Landmark
Fundamental
Important
Valuable
Useful

Strength of evidence:

Convincing: Appropriate and validated methodology in line with current state-of-the-art

Exceptional
Compelling
Convincing
Solid
Incomplete
Inadequate

During the peer-review process the editor and reviewers write an eLife Assessment that summarises the significance of the findings reported in the article (on a scale ranging from landmark to useful) and the strength of the evidence (on a scale ranging from exceptional to inadequate). Learn more about eLife Assessments

Abstract
Introduction
Results
Discussion
Appendix 1
Appendix 2
Data availability
References
Article and author information
Metrics

Abstract

Cell polarization is a critical process that separates molecular species into two distinct regions in prokaryotic and eukaryotic cells, guiding biological processes such as cell division and cell differentiation. Although several underlying antagonistic reaction-diffusion networks capable of setting up cell polarization have been identified experimentally and theoretically, our understanding of how to manipulate pattern stability and asymmetry remains incomplete, especially when only a subset of network components is known. Here, we present numerical results to show that the polarized pattern of an antagonistic 2-node network collapses into a homogeneous state when subjected to single-sided self-regulation, single-sided additional regulation, or unequal system parameters. However, polarity restoration can be achieved by combining two modifications with opposing effects. Additionally, spatially inhomogeneous parameters favoring respective domains stabilize their interface at designated locations. To connect our findings to cell polarity studies of the nematode Caenorhabditis elegans zygote, we reconstituted a 5-node network where a 4-node circuit with full mutual inhibitions between anterior and posterior is modified by a mutual activation in the anterior and an additional mutual inhibition between the anterior and posterior. Once again, a generic set of kinetic parameters moves the interface towards either the anterior or posterior end, yet a polarized pattern can be stabilized through tuning of one or more parameters coupled to intracellular or extracellular spatial cues. A user-friendly software, PolarSim, is constructed to facilitate the exploration of networks with alternative node numbers, parameter values, and regulatory pathways.

Introduction

Cell polarization is a biophysical and biochemical process where a cell acquires spatial anisotropy by establishing directional gradients of molecular species across its membrane or in the cytosol (Knoblich, 2001; Goodrich and Strutt, 2011). Polarity establishment is an essential step in a wide range of biological phenomena, including embryonic development, wound healing, immune activity, and so forth (Bilder et al., 2000; Etienne-Manneville and Hall, 2002). During cytokinesis, a polarized cell can allocate its molecular contents unequally to its daughter cells, leading to asymmetric cell size and cell fate (Macara, 2004). At the multicellular level, a group of polarized cells can undergo collective movements and constitute stereotypical architectures (Bilder et al., 2000; Etienne-Manneville and Hall, 2002). Thus, loss or disorder of cell polarization could severely violate normal biological processes, for example, resulting in embryonic lethality and cancerous tumors (Bilder et al., 2000; Kim et al., 2007). To this day, cell polarization has been a long-term research focus in cell and developmental biology, where more and more efforts have been dedicated to uncovering both the regulatory pathways and design principles involved (Doe and Bowerman, 2001; Tostevin and Howard, 2008; Chau et al., 2012; Koorman et al., 2016).

Much of our knowledge on cell polarization is derived from the zygote of the nematode Caenorhabditis elegans (referred to as C. elegans), which has served as a prominent model organism for studying cell polarization (Rose and Gönczy, 2014; Lang and Munro, 2017). Following fertilization, the C. elegans zygote becomes polarized, with the sperm entry site (located at one pole of the ellipsoidal egg) designating the posterior of the embryo (Goldstein and Hird, 1996; Motegi et al., 2011). Driven by cell polarization, the embryo proceeds through four successive rounds of asymmetric cell divisions and produces four somatic founder cells sequentially (Sulston et al., 1983; Hubatsch et al., 2019).

The asymmetric division in the C. elegans zygote is governed by a protein family termed partitioning-defective protein (PAR) (Kemphues et al., 1988). The initial PAR family only consisted of the PAR-3/PAR-6/PKC-3 complex and PAR-1/PAR-2 complex, which are stably accumulated in the anterior and posterior domains on the membrane of the zygote before its division, respectively (Cuenca et al., 2003). During the establishment phase triggered by the sperm entry that defines the zygote’s posterior, PAR-3/PAR-6/PKC-3, initially distributed across the entire cell membrane, are anteriorly directed and localized via the contraction of the cortical actomyosin marked by non-muscle myosin II (NMY-2), which allows the PAR-2 (along with the recruitment of PAR-1) to be enriched on the posterior membrane conversely (Rose and Gönczy, 2014; Lang and Munro, 2017). Subsequently, the maintenance phase follows as the cortical actomyosin contracts halfway across the zygote, with cortical NMY-2 being regulated downstream by PAR proteins and playing a minor role in maintaining the cell polarization pattern (Goehring et al., 2011a; Beatty et al., 2013; Rose and Gönczy, 2014). The cell polarization patterns, defined by the spatial concentration distributions of PAR proteins, are stably maintained to guide the assembly of the cytokinetic ring, which is vital for the first asymmetric cell division (Goehring et al., 2011a; Rose and Gönczy, 2014; Jankele et al., 2021).

A series of experiments (incl., RNA interference, immunoprecipitations, fluorescence recovery after photobleaching, and time-lapse single-molecule imaging) further uncovered the mutual inhibition between PAR-3/PAR-6/PKC-3 and PAR-1/PAR-2 upon their association with the membrane and demonstrated its essential role in their robust spatial separation (Cuenca et al., 2003; Motegi et al., 2011). Soon, theoretical and numerical studies proved that mutual inhibition forms the backbone (hereafter referred to as an antagonistic 2-node network) of a polarized pattern (Tostevin and Howard, 2008; Chau et al., 2012). In the C. elegans zygote, the antagonism between these two protein groups directs the uneven distribution in downstream cell fate determinants (e.g., PIE-1 and P granules), which is inherited by the two daughter cells during the following asymmetric cell division, leading to cell differentiation (Schubert et al., 2000; Wang and Seydoux, 2013).

In recent years, more proteins have been identified to significantly interact with the antagonistic 2-node network during the maintenance phase of C. elegans zygotic cell polarization, such as CDC-42, LGL-1, and CHIN-1 (Beatty et al., 2010; Kumfer et al., 2010; Sailer et al., 2015). The addition of these players increases the complexity of the cell polarization network tremendously, not only because of the explosion of the associated kinetic parameters, but also their role in related cellular processes such as cytoskeleton organization and localization at the poles that determines cell division dynamics (Ajduk and Zernicka-Goetz, 2016; Lim et al., 2021). Therefore, this growing complexity underscores the importance of developing a more generalized understanding of the stability and asymmetry of cell polarization patterns, which will help to streamline experimental data analysis and interpretation. In a separate advancement, the theoretical knowledge acquired from nature is urgently needed for de novo construction of cells with designated characteristics, where the capability of cell polarization has been a design target for over a decade ago (Tostevin and Howard, 2008; Chau et al., 2012; Lin et al., 2021; Watson et al., 2023).

In this work, we focus on the generation of stable asymmetric patterns in both the widely-used 2-node network and a more realistic C. elegans 5-node network. Starting from a symmetric antagonistic network, a polarized pattern can be stabilized at any interface location between the two antagonistic domains, thereby manifesting translational symmetry. Unbalanced modification of kinetic parameters triggers movement of the interface in favor of one of the coexisting domains. Three types of such unbalanced modification are reported: single-sided self-regulation, single-sided additional regulation, and unequal system parameters. Nevertheless, the combination of two or more unbalanced modifications can recover pattern stability through fine-tuning of kinetic parameters. Alternatively, we show that the interface can also be stabilized at a designated location with a step-like spatial profile of one or more of the kinetic parameters, with values leading to opposing velocities when the interface is displaced. Intriguingly, such a program strategy is found to be employed in the C. elegans cell polarization network to maintain pattern stability along with considerable parameter robustness, while inducing pattern asymmetry by interface localization control.

Results

A computational pipeline for simulating cell polarization

To investigate both the simple network and the realistic network consisting of various node numbers and regulatory pathways (Goehring et al., 2011b; Lang and Munro, 2017), we propose a computational pipeline for numerical exploration of a given reaction-diffusion network’s dynamics, specifically targeting the maintenance phase of stable cell polarization after its initial establishment (Motegi et al., 2011; Goehring et al., 2011b; Seirin-Lee et al., 2020b). Numerous biological experiments on different organisms have demonstrated that such a reaction-diffusion network typically comprises two groups of molecular species. Each group, once associated with a part of the cell membrane from cell cytosol, inhibits the association of the other group in the same region, thereby creating two types of distinct domains (Knoblich, 2001; Doe and Bowerman, 2001). At the interface between the two domains, the membrane association of proteins in either group is compromised due to the elevated level of antagonists. Nevertheless, one of the domains may expand at the expense of the other, leading to a finite interface velocity in general. For simplicity, three assumptions previously used in research are adopted to establish a reaction-diffusion model to describe the dynamics of each molecular species (denoted by $[X]$ ) during cell polarization:

The cellular space is reduced to a one-dimensional line of length $L = 0.5$ , where $x = - 0.25$ (anterior, where the concentration of molecular species accumulated on the cell membrane is denoted by $[A_{m}]$ ) and 0.25 (posterior, where the molecular species accumulated on the cell membrane is denoted by $[P_{m}]$ ) are its two poles (Figure 1a; Gross et al., 2019; Seirin-Lee, 2020a).
A molecular species $[X]$ can associate with the cell membrane from cell cytosol at a rate $F_{o n}^{X} (x, t)$ and dissociate from the cell membrane into cell cytosol at a rate $F_{o f f}^{X} (x, t)$ , where $t$ represents time. Both rates include leaky term and regulatory pathways affected by other molecular species (Seirin-Lee and Shibata, 2015; Seirin-Lee et al., 2020b).
Based on previous experimental measurements, it has been reported that the diffusion rate of molecular species involved in cell polarization (i.e., PAR-2 and PAR-6) is two orders of magnitude higher in the cell cytosol compared to the cell membrane (Goehring et al., 2011a; Blanchoud et al., 2015; Gross et al., 2019; Lim et al., 2021). Consequently, we adopted 2.748×10^–5 and 1.472×10^–5 from previous research (nondimensionalized based on experimental measurements) (Goehring et al., 2011a; Goehring et al., 2011b; Seirin-Lee et al., 2020b), as diffusion coefficients of anterior and posterior molecular species, while the cytoplasmic molecules are regarded as well-mixed with infinite diffusion (Kravtsova and Dawes, 2014; Gross et al., 2019; Morita and Seirin-Lee, 2021). We further assume that these molecular species are sufficiently abundant in the pool of cell cytosol, where the concentration is conserved as constant $[X_{c}]$ , although the analysis below can be extended to the more general case (Kravtsova and Dawes, 2014; Goehring et al., 2011b).

Figure 1 with 6 supplements see all

Download asset Open asset

The single modification on the antagonistic 2-node network causes the collapse of cell polarization pattern.

(a) Basic network with the transition plane (mid-plane) marked by gray dashed line. (b) Two subtypes of single-sided self-regulation. (c) Four subtypes of single-sided additional regulation. The corresponding spatial concentration distribution of $[L_{m}]$ at $t = 0, 100, 200, 300, 400, a n d 500$ , and the subtle difference between left and right panels are detailed in Figure 1—figure supplement 3. (d) Two subtypes of unequal system parameters, exemplified by unequal inhibition intensity and unequal cytoplasmic concentration. For each network, the corresponding spatial concentration distribution of $[A_{m}]$ and $[P_{m}]$ at $t = 0, 100, 200, 300, 400, a n d 500$ are shown beneath with a color scheme listed in the bottom left corner. Note that within a network, normal arrows and blunt arrows symbolize activation and inhibition respectively.

For the concentration of each molecular species $[X]$ on the cell membrane, $[X_{m}]$ , the evolution equation is described by:

\frac{\partial [X_{m}] (x, t)}{\partial t} = D_{m}^{X} \frac{\partial^{2}}{\partial x^{2}} [X_{m}] + F_{o n}^{X} (x, t) [X_{c}] - F_{o f f}^{X} (x, t) [X_{m}]

where $D_{m}^{X}$ is the diffusion coefficient of molecular species $[X]$ on the cell membrane. The first term on the right side represents the diffusion across the cell membrane; the second term represents the association with the cell membrane from cell cytosol; the third term represents the dissociation from the cell membrane into cell cytosol. These regulated association and dissociation are given by the Hill equation (Seirin-Lee and Shibata, 2015; Seirin-Lee, 2021):

F_{o n}^{X} (x, t) = γ_{X} + \sum_{Y} \frac{q_{Y 2}^{X} {[Y_{m}]}^{n_{Y q}^{X}}}{1 + q_{Y 1}^{X} {[Y_{m}]}^{n_{Y q}^{X}}}

F_{o f f}^{X} (x, t) = α_{X} + \sum_{Y} \frac{k_{Y 2}^{X} {[Y_{m}]}^{n_{Y k}^{X}}}{1 + k_{Y 1}^{X} {[Y_{m}]}^{n_{Y k}^{X}}}

where $F_{o n}^{X}$ consists of a basal on-rate (i.e., association rate) $γ_{X}$ and recruitment by other membrane-bound molecular species ( $Y \neq X$ ) and itself ( $Y = X$ , i.e., self-activation), quantified by positive parameters $q_{Y i}^{X} (i = 1, 2)$ ; $F_{o f f}^{X}$ consists of a basal off-rate (i.e., dissociation rate) $α_{X}$ and exclusion by other molecular species ( $Y \neq X$ ) and itself ( $Y = X$ , i.e., self-inhibition), quantified by positive parameters $k_{Y i}^{X} (i = 1, 2)$ ; the Hill coefficients $n_{Y q}^{X}$ and $n_{Y k}^{X}$ are set to 2 as used before (Seirin-Lee and Shibata, 2015; Seirin-Lee, 2021). To simplify the model for numerical exploration with an affordable computational cost, we set the same responsive concentration $q_{Y 1}^{X} = k_{Y 1}^{X} = k_{1}$ for both activation and inhibition pathways; the inhibition intensities $k_{Y 2}^{X}$ are set to 1 as a numerically normalized value. We also set self and mutual activation parameters $q_{Y 2}^{X} = q_{2}$ to be the same for synergistic proteins in the same group. Finally, the basal association and dissociation rates are regarded the same so that these spatially independent effects can be neutralized, in other words, $γ_{X} = α_{X} = γ$ . Consequently, there are only three independent parameters governing this system: $γ$ , $k_{1}$ , and $q_{2}$ . A detailed description of the parameters is listed in Appendix 1—table 1. This dimensionally reduced parameter configuration is sufficient to describe the temporal evolution of cell polarization pattern on the membrane under the well-mixed cytoplasmic protein concentration. All the simplified parameter value assignments above will be extensively explored by giving different values to different molecular species and pathways as well as setting a heterogeneous spatial distribution.

Simulations are performed by systematically scanning the dimensionless parameter set $Ω (γ, k_{1}, q_{2})$ on a three-dimensional (3D) grid $γ \in [0,0.05]$ in steps $∆ γ = 0.001$ , $k_{1} \in [0,5]$ in steps ${∆ k}_{1} = 0.05$ , and $q_{2} \in [0,0.05]$ in steps $∆ q_{2} = 0.001$ . The derived nondimensionalized parameter attains values on the same order of magnitude as those observed in reality, confirming the fidelity of the proposed model in representing the real system (Appendix 1). When it comes to a statement of switch from pattern stability to instability, the three-parameter simplification will be removed to free the parameter value assignments (i.e., uniform random variables $γ \in [0, 0.05], α \in [0, 0.05], k_{1} \in [0, 5], k_{2} \in [0, 5], q_{1} \in [0, 5], q_{2} \in [0, 0.05]$ , and $[X_{c}] \in [0.05,5]$ added for setting up a network with cell polarization pattern stability, searched by the Monte Carlo method) concerning all nodes and regulations, auxiliarily validating the conclusion (Fishman, 1995).

Focusing on the cell polarization pattern stability and asymmetry during the maintenance phase, we simplify and set the initial distribution of molecular species on the cell membrane, $[X_{m}] (x, 0)$ , which is polarized during the establishment phase through active actomyosin contractility and flow (Munro et al., 2004; Kravtsova and Dawes, 2014), using the sigmoid function:

{\begin{cases} [X_{m}] (x, 0) = 1 - \frac{1}{1 + e^{- 20 x}}, f o r a n t e r i o r m o l e c u l a r s p e c i e s \\ [X_{m}] (x, 0) = \frac{1}{1 + e^{- 20 x}}, f o r p o s t e r i o r m o l e c u l a r s p e c i e s \end{cases}

Then the set of partial differential Equation 1 evolves for 500 steps with a time step of 1 (Figure 1—figure supplement 1). The solution at $t = 500$ is saved for further analysis if:

(i) All the molecular species still have a polarized pattern, as defined by $[X_{m}] (- L / 2, 500) > 0.5$ and $[X_{m}] (L / 2, 500) < 0.5$ for the anterior molecular species and $[X_{m}] (- L / 2, 500) < 0.5$ and $[X_{m}] (L / 2, 500) > 0.5$ for the posterior molecular species.

(ii) The cell polarization pattern is stable or nearly intact over time for each molecular species, as defined by

| \frac{[X_{m}] (x_{0}, 500) - [X_{m}] (x_{0}, 499)}{[X_{m}] (x_{0}, 499)} | < 10^{- 4}, \forall x_{0} \in [- 0.25, 0.25]

Implemented on Matlab 2022b (The MathWork Inc, 2022), such a computational pipeline can describe the spatiotemporal dynamics of any reaction-diffusion network design, for evaluating its capability of maintaining cell polarization (Figure 1—figure supplement 1). Meanwhile, we construct an automatic software, PolarSim, for users to explore networks with arbitrarily set node numbers, parameter values, and regulatory pathways (Appendix 2, Source code 1).

An unbalanced network structure or parameter leads to the collapse of a polarized pattern

Previous experimental and theoretical discoveries have uncovered the mutual inhibition between two molecular species as a fundamental design capable of generating cell polarization (Kemphues et al., 1988; Cuenca et al., 2003; Tostevin and Howard, 2008; Chau et al., 2012). Thus, we utilize the completely symmetric antagonistic 2-node network to investigate the behavior of its cell polarization pattern when the network structure or parameter is modified. The two nodes (i.e., molecular species) placed in the anterior and posterior are denoted by $[A]$ and $[P]$ . As there is no activation here then $q_{2}$ is not considered, the computational pipeline above establishes a total of 122 viable parameter sets $(γ, k_{1})$ (122 among 5,151 sets, ~2.37%) that can achieve stable cell polarization. We use $(γ = 0.05, k_{1} = 0.05)$ as a representative to show how the corresponding cell polarization pattern behaves under alternative initial conditions and elementary modification (Figure 1, Figure 1—figure supplement 2).

Regarding the initial conditions, while uniform initial distributions maintain a homogeneous pattern (Figure 1—figure supplement 2a, b), a polarized initial distribution reliably maintains the cell polarization regardless of variations in the initial distribution, such as shifts in transition planes (Figure 1—figure supplement 2c, d), weak polarization (Figure 1—figure supplement 2e, f), or asymmetric curve shapes between the two molecular species (Figure 1—figure supplement 2f, g). In contrast, an initial distribution with only Gaussian noise produces a multipolar pattern, with the two molecular species still locally excluding each other (Figure 1—figure supplement 2h). Altogether, these findings demonstrate the necessity of additional mechanisms, such as the active actomyosin contractility and flow, for establishing initial cell polarization independently of the reaction-diffusion network during the establishment phase (Cuenca et al., 2003; Gross et al., 2019).

Regarding the elementary modification, a total of three types are exerted on the node $[A]$ :

Single-sided self-regulation (Figure 1b): a self-activation ( $q_{A 2}^{A} = 0.012$ ) or self-inhibition ( $k_{A 2}^{A} = 0.1$ ) is added on the node $[A]$ .
Single-sided additional regulation (Figure 1c, Figure 1—figure supplement 3): a new node $[L]$ is added in the anterior or posterior, with mutual activation ( $q_{L 2}^{A} = q_{A 2}^{L} = 0.012$ ) or inhibition ( $k_{L 2}^{A} = k_{A 2}^{L} = 0.025$ ) with $[A]$ . The location of the additional node $[L]$ decides the initial polarized pattern represented by the sigmoid function, where an anterior node starts with a high concentration in the anterior and a low concentration in the posterior, and vice versa.
Unequal system parameters (Figure 1d): on one hand, the inhibition intensity $k_{2}$ from $[P]$ to $[A]$ is increased from 1 to 1.6; on the other hand, the cytoplasmic concentration of $[A]$ is increased from 1 to 1.25. The other system parameters (e.g., basal on-rate $γ$ and basal off-rate $α$ ) will be explored independently in the next section.

In comparison with the stable cell polarization pattern generated by the basic network, the spatial concentration distribution dynamics of eight modified conditions are simulated as shown in Figure 1 and Figure 1—video 1, where the domain of one molecular species keeps invading the domain of the other one, ending in a homogeneous distribution. Such pattern instability caused by unbalanced network modification still holds when the network is initially set up using different initial conditions and parameters (i.e., $γ$ , $α$ , $k_{1}$ , $k_{2}$ , and $[X_{c}]$ ) assigned with various values concerning all nodes and regulations (Figure 1—figure supplement 4, Figure 1—figure supplement 5 – 1^st row, Appendix 1—table 2). This intuitively reveals that the antagonistic 2-node reaction-diffusion network is prone to become unstable with its interface keeps moving when the perturbation is introduced to the network structure or parameter.

The combination of two modifications can recover the cell polarization pattern stability

Since a single modification of reaction-diffusion network structure or parameter is enough to break the cell polarization pattern, an appealing question just comes up: how can the network be designed to maintain pattern stability considering such modifications? This is crucial for cell polarization in reality where stability is essential for cell function and survival; without stability, the spatial information defined by the pattern is inaccurate for guiding downstream biological events such as cell division and cell differentiation (Wang and Seydoux, 2013; Hubatsch et al., 2019). The simplest idea for stability recovery is to combine two kinds of modifications with opposite trends, for example, adding self-activation and additional anteroposterior mutual inhibition on $[A]$ simultaneously. For the three types of modification, we arbitrarily select one subtype within each of them from Figure 1b–d so three combinatorial networks are constituted, all of whose interfaces turn out to be stabilized finally (Figure 2, Figure 2—video 1).

Figure 2 with 1 supplement see all

Download asset Open asset

The combination of two opposite modifications recovers the cell polarization pattern stability.

The basic network and the ones added with a single modification are shown in the 1^st and 2^nd columns, respectively; the three combinatorial networks composed of any two of the three single modifications are shown in the 3^rd column. For each network, the corresponding spatial concentration distribution of $[A_{m}]$ and $[P_{m}]$ at $t = 0, 100, 200, 300, 400, a n d 500$ are shown beneath with a color scheme listed on the right. Here, the value assignments on the modifications in the 3^rd column are as follows: $q_{A 2}^{A} = 0.012$ and $k_{L 2}^{A} & k_{A 2}^{L} = 0.01$ for 1^st row, $q_{A 2}^{A} = 0.012$ and $k_{P 2}^{A} = 1.24$ for 2^nd row, and $q_{A 2}^{A} = 0.012$ and $k_{A 2}^{P} & k_{A 2}^{L} = 2$ for 3^rd row. Note that within a network, normal arrows and blunt arrows symbolize activation and inhibition respectively.

To more precisely elucidate how two opposite modifications coordinate the pattern stability together, we make use of the network with self-activation and additional inhibition on $[A]$ (shown in the top right corner of Figure 2), where the two corresponding intensities $q_{A 2}^{A}$ and $k_{L 2}^{A} & k_{A 2}^{L}$ compose a phase diagram that distinguishes the final state of the reaction-diffusion pattern. The pattern movement velocity is defined as follows:

\begin{matrix} v_{M} (t) = \frac{1}{N} \sum_{X} \frac{\int_{- \frac{L}{2}}^{\frac{L}{2}} | [X_{m}] (x, t) - [X_{m}] (x, t - 1) | d x}{L} \end{matrix}

where $N$ is the total number of the molecular species. Here, we classify the final state by calculating the pattern movement velocity at the final time $t = 500$ ; the region with $v_{M} < 10^{- 4}$ around the diagonal line is marked as the polarized state, while the regions upon and beneath it are homogeneous states dominated by $[A]$ and $[P]$ respectively (Figure 3a). Therefore, the triphase diagram and the exemplary patterns generated by the parameter assignments ①②③ shown in Figure 3a suggest the necessity of the balance between two opposite modifications for setting up a stable cell polarization pattern.

Figure 3

Download asset Open asset

The balance between system parameters is needed for maintaining pattern stability.

(a) The phase diagram between $q_{A 2}^{A}$ and $k_{L 2}^{A} & k_{A 2}^{L}$ in the network modified by self-activation (quantified by $q_{A 2}^{A}$ ) and additional inhibition (quantified by $k_{L 2}^{A} & k_{A 2}^{L}$ ) on $[A]$ . The representative parameter assignment for each phase is marked with ① (*i.e.*, $q_{A 2}^{A} = 0.015$ and $k_{L 2}^{A} & k_{A 2}^{L} = 0$ with a homogeneous state dominated by $[A]$ ), ② (*i.e.*, $q_{A 2}^{A} = 0.015$ and $k_{L 2}^{A} & k_{A 2}^{L} = 0.0135$ with a stable polarized state), and ③ (*i.e.*, $q_{A 2}^{A} = 0.02$ and $k_{L 2}^{A} & k_{A 2}^{L} = 0$ with a homogeneous state dominated by $[P]$ ). The corresponding spatial concentration distribution of $[A_{m}]$ and $[P_{m}]$ at $t = 0, 100, 200, 300, 400, a n d 500$ are shown around the phase diagram with a color scheme listed on top. (b) The phase diagram between responsive concentration $k_{1}$ , basal on-rate $γ$ , basal off-rate $α$ , cytoplasmic concentration $[X_{c}]$ , and inhibition intensity $k_{2}$ . For each phase diagram in (a) (b), the final state dominated by $[A]$ or $[P]$ or stably polarized is colored in orange, green, and gray, respectively. Note that within a network, normal arrows and blunt arrows symbolize activation and inhibition respectively.

Apart from the unequal inhibition intensity and cytoplasmic concentration mentioned in Figure 1d, we further ask if all the kinetic parameters with biophysical significance (incl., inhibition intensity $k_{2}$ , responsive concentration $k_{1}$ , basal on-rate $γ$ , basal off-rate $α$ , and cytoplasmic concentration $[X_{c}]$ ) also require a balance for achieving pattern stability. For this purpose, we generate the phase diagrams between each of them and the inhibition intensity $k_{A 2}^{P}$ . Fascinatingly, monotonic correlations exist between all those system parameters, suggesting that they can be tuned to maintain pattern stability (Figure 3b). In the gray region, the symmetric-broken parameters can be weighed against each other to realize the pattern movement velocity close to zero.

The velocity and position of the interface can be adjusted by setting up spatial cues

The required balance between parameters provides insight into controlling the pattern stability. However, in biological organisms, factors such as cellular actomyosin flow and intercellular signals could lead to changeable parameter values and molecular concentration in different parts of space (Arata et al., 2010; Beatty et al., 2013). How do cells maintain cell polarization pattern stability in the case of nonuniform initial conditions and parameter values, namely finding zero-velocity solutions at interfaces? Here, given a specific pattern, its position of the transition plane or interface $x_{T} (t)$ is defined by the mean position of all the molecular species where the absolute derivative value of the concentration curve reaches its maximum.

First, based on the symmetric 2-node network ( $γ = 0.05, k_{1} = 0.05$ , Figure 4a, the same simulation as the one in Figure 1a), we alter its initial condition by shifting the interface position or exerting a uniform or Gaussian noise distribution without any positional information (Figure 1—figure supplement 2); it is found that the initial polarized pattern is essential for constricting only one single interface with predetermined position (Figure 1—figure supplement 2c, d, e, f, g), while the uniform or Gaussian noise distribution results in no or multiple interfaces (Figure 1—figure supplement 2a, b, h). This indicates the importance of additional mechanisms (e.g., the cortical flow on the cell membrane that transports proteins related to cell polarity) establishing the initial cell polarization pattern before the pattern enters the maintenance phase, when the network structure and parameters account for its stability in the long term (Cuenca et al., 2003; Motegi et al., 2011).

Figure 4 with 1 supplement see all

Download asset Open asset

Adopting parameter sets corresponding to opposite interface velocities on two sides of the interface, the zero-velocity interface position of the polarity pattern is shiftable and regulatable.

(a) Spatially inhomogeneous parameters (2^nd column) of a symmetric 2-node network generate a symmetric pattern (1^st column), from the same simulation in Figure 1a. (b) Using the parameter combination with the posterior-shifting interface on the left and anterior-shifting interface on the right, a stable polarity pattern (1^st column) can be remained by increasing $k_{A 2}^{P}$ to 1.5 at $x < 0$ and decreasing $γ_{A}$ to 0.01 at (2^nd column). (**c–d**) The stable interface position can be optionally adjusted by setting the change position of the step-up function (2^nd column). (c) As in (b), but changing the step position to $x = - 0.1$ , the interface stabilizes around $x = - 0.1$ (1^st column). (d) As in (b), but changing the step position to $x = 0.1$ , the interface stabilizes around $x = 0.1$ (1^st column). For each parameter set, the corresponding spatial concentration distribution of $[A_{m}]$ and $[P_{m}]$ at $t = 0, 100, 200, 300, 400, a n d 500$ are shown beneath with a color scheme listed in the bottom right corner. Note that all the interface positions in (**a–d**) (1^st column) are marked by gray dashed lines. Note that all the spatially homogeneous parameter distributions in (a) (2^nd column) are marked by gray dashed lines in (**b–d**) (2^nd column).

Next, to probe into how inhomogeneous biophysical parameter values can control the interface in concert, we employ the following adjustment on the parameter set used in Figure 4a:

(i) Increasing the inhibitory intensity of $[A]$ on $[P]$ , $k_{A 2}^{P}$ , from 1 to 1.5, leads to the interface continuously shifting to the posterior, and $[A]$ finally dominates.

(ii) Decreasing the basal on-rate of $[A]$ , $γ_{A}$ , from 0.05 to 0.01, leads to the interface continuously shifting to the anterior, and $[P]$ finally dominates.

When these two systems are combined, the stable polarity pattern remains with (i) used in the region $x < 0$ while (ii) used in $x > 0$ (Figure 4b). By moving the step position of the parameter function to $x = - 0.1$ (Figure 4c) and $x = 0.1$ (Figure 4d), the interface stabilizes around the step, which means the interface position is tunable (Figure 4—video 1). This indicates that by using parameter sets with values corresponding to opposite interface velocities on two sides of the interface, the interface velocity and profile of the polarity pattern can be precisely regulated. Such inhomogeneous spatial distribution of parameters imitates signaling perception, which provides a strategy for more robust control of the adjustable interface position in response to variable cues.

Reconstruction of the molecular interaction network in C. elegans zygote and its design principle of structure and parameter trade-off

In the simple 2-node network, the cell polarization pattern stability could be broken by a single modification (Figure 1) and recovered by two combinatorial modifications (Figures 2 and 3), or with the zero-velocity interface position regulated by spatially inhomogeneous parameters (Figure 4). We further ask if such a fundamental rule is employed in the cell polarization network programming in a real system. To this end, we focus on the zygote of the nematode C. elegans, which has been a popular model for cell polarization study from both experimental and theoretical perspectives for more than three decades (Kemphues et al., 1988; Goehring et al., 2011b; Lang and Munro, 2017; Lim et al., 2021). By conducting an exhaustive literature search (a total of 19 references), we summarize the molecular interaction network in C. elegans zygote that consists of five interacting molecular species: PAR-3/PAR-6/PKC-3 complex (abbr., $[A]$ ) and CDC-42 protein (abbr., $[C]$ ) accumulated in the anterior, and PAR-1/PAR-2 complex (abbr., $[P]$ ), LGL-1 (abbr., $[L]$ ) and CHIN-1 (abbr., $[H]$ ) accumulated in the posterior (Figure 5a). The detailed description of the biochemical mechanism of each regulatory pathway as well as the corresponding supporting references are listed in Supplementary file 1. Based on the network obtained experimentally, the computational pipeline described in Section (A computational pipeline for simulating cell polarization) establishes a total of 602 viable parameter sets $Ω (γ, k_{1}, q_{2})$ (602 among 262,701 sets, ~0.229%) that can achieve stable cell polarization (Figure 5d and e).

Figure 5 with 6 supplements see all

Download asset Open asset

The molecular interaction network in C. *elegans* zygote and its natural advantages in terms of pattern stability, balanced network configuration, viable parameter sets, and parameter robustness.

(a) The schematic diagram of the network composed of five molecular species, each of which has a polarized spatial concentration distribution on the cell membrane shown beneath. The left nodes PAR-3/PAR-6/PKC-3 (*i.e.*, $[A]$ ) and CDC-42 (*i.e.*, $[C]$ ) exhibit high concentration in the anterior pole and low concentration in the posterior pole, schemed by the purple line beneath; the right nodes PAR-1/PAR-2 (*i.e.*, $[P]$ ), LGL-1 (*i.e.*, $[L]$ ), and CHIN-1 (*i.e.*, $[H]$ ) exhibit low concentration in the anterior pole and high concentration in the posterior pole, schemed by the cyan line beneath. Note that within the network, normal arrows and blunt arrows symbolize activation and inhibition respectively. (**b–d**) The structure of 4-Node, LGL-1^-, and WT networks (1^st row). The final spatial concentration distribution ( $t = 500$ ) averaged over all established viable parameter sets for each molecular species (*i.e.*, 62 among 262,701 sets,~0.024%, for the 4-Node network; 62 among 5,151 sets,~1.2%, for the LGL-1^- network; 602 among 262,701 sets,~0.229%, for the WT network), shown by a solid line (2^nd row). For each position, MEAN ±STD (*i.e.*, standard deviation) calculated with all viable parameter sets is shown by shadow. The pattern movement velocity ( $v_{M}$ ) over evolution time ( $t$ ) (3^rd row). For each subfigure in the 3^rd row, each line with a unique color represents the simulation of a unique viable parameter set, and $v_{M} = 10^{- 4}$ is marked by a gray dashed line. For (**a–d**), the legend for the relationship between molecular species and corresponding color is placed in the bottom right corner. (e) The viable parameter sets of WT (blue points; 602 among 262,701 sets, ~0.229%) and 4-Node networks (red points; 62 among 5,151 sets, ~1.2%). (f) The interplay between $[A]$ ~ $[C]$ mutual activation and $[A]$ ~ $[L]$ mutual inhibition on the interface velocity. The contour map of the interface velocity $v_{I}$ with different parameter combinations of $q_{C 2}^{A} & q_{A 2}^{C}$ and $k_{L 2}^{A} & k_{A 2}^{L}$ represents its moving trend. The darker the red is, the faster the interface is moving toward the posterior pole; the darker the blue is, the faster the interface is moving toward the anterior pole; the closer the color is to white, the stabler the interface is. (g) The averaged pattern error in a perturbed condition compared to the original 4-Node, LGL-1^-, and WT networks. The Gaussian noise is simultaneously exerted on all the parameters $γ$ , $α$ , $k_{1}$ , $k_{2}$ , $q_{1}$ and $q_{2}$ , where standard deviation $σ$ denotes the noise amplification. For a specific noise level, the pattern error is averaged over 1,000 independent simulations.

To verify whether the computational pipeline is reliable enough to simulate the protein distribution dynamics of the C. elegans network in vivo, we further reproduce a series of perturbation experiments reported previously. It’s worth noting that the reconstructed C. elegans 5-node network appears to be an integration of an antagonistic 2-node network doubled into an antagonistic 4-node network (Figure 5b) with full connections, and the $[A]$ ~ $[C]$ mutual activation (Figure 5c) and $[A]$ ~ $[L]$ mutual inhibition (Figure 5d). As the symmetric 2-node network has been revealed to be a fundamental structure capable of stable cell polarization (Tostevin and Howard, 2008; Goehring et al., 2011b; Chau et al., 2012), hereafter we focus on how the added interactions involved with $[C]$ (the protein accumulated in the anterior and with an additional mutual activation with $[A]$ in the anterior) and $[L]$ (the protein accumulated in the posterior and with an additional mutual inhibition with $[A]$ in the anterior) impact the pattern stability concordantly (Figure 5—figure supplement 1).

Taking $[L]$ as the first example for its current absence in theoretical model analysis (Seirin-Lee et al., 2020b; Lim et al., 2021), a series of experimental groups in different conditions were conducted before by knocking down or overexpressing $[L]$ , followed by a measurement of the lethality in embryo individuals (defined by death rate, a consequence of polarity loss in the zygote, which is essential for asymmetric cell division and cell differentiation) (Beatty et al., 2010; Beatty et al., 2013; Rodriguez et al., 2017; Jankele et al., 2021). Here, we utilize the pattern error in a mutant (abbr., MT) embryo compared to the wild-type (abbr., WT) one, $E r r o r = \frac{1}{M} \sum_{i = 1}^{M} \frac{\int_{- \frac{L}{2}}^{\frac{L}{2}} |{[X_{m}]}_{W T} (x, 500) - {[X_{m}]}_{M T} (x, 500)| d x}{L}$ , to represent the level of polarity loss in simulation and lethality in experiment (measured in previous literature), where $M = 602$ is the number of viable parameter sets (602 among 262,701 sets, ~0.229%). The first experimental group is that the double depletion on $[P]$ and $[L]$ leads to much more lethality than single depletion on either $[P]$ or $[L]$ , which is recapitulated by Figure 5—figure supplement 2a, b (cell polarization pattern characteristics in simulation: larger pattern error in double depletion on $[P]$ and $[L]$ , compared to single depletion on either $[P]$ or $[L]$ ) (Hoege et al., 2010). The second experimental group is that the overexpression of $[L]$ lowers the lethality induced by single depletion on $[P]$ and such effect is weakened when $[H]$ is depleted as well, which is recapitulated by Figure 5—figure supplement 2c,d (cell polarization pattern characteristics in simulation: smaller pattern error in single depletion on $[P]$ supplemented by overexpression of $[L]$ , compared to single depletion on $[P]$ only; less smaller pattern error in double depletion on $[P]$ and $[L]$ supplemented by overexpression of $[L]$ , compared to single depletion on $[P]$ only) (Beatty et al., 2010; Beatty et al., 2013).

The second example of theoretical model analysis involves $[C]$ . A series of experimental groups in different conditions were conducted before by knocking down or overexpressing $[C]$ or blocking the mutual activation between $[A]$ and $[C]$ , followed by an observation of the cell polarity in embryo individuals. The first experimental group is that the single depletion on $[C]$ results in polarity defects and mislocalization of PAR-6 and PAR-2, which is recapitulated by Figure 5—figure supplement 3a, b (cell polarization pattern characteristics in simulation: less polarized distributions of $[A]$ and $[P]$ , compared to WT) (Gotta et al., 2001; Aceto et al., 2006). The second experimental group is that the overexpression of $[C]$ increases the size of the anterior domain, which is recapitulated by Figure 5—figure supplement 3a, c (cell polarization pattern characteristics in simulation: posteriorly shifted interfaces, compared to WT) (Aceto et al., 2006). The third experimental group is that blocking the interaction between PAR-6 and CDC-42 results in polarity defects very similar to those associated with single depletion on CDC-42, which is recapitulated by Figure 5—figure supplement 3a, d (cell polarization pattern characteristics in simulation: less polarized distribution of $[A]$ in both conditions, compared to WT) (Aceto et al., 2006).

To sum up, apart from the cell polarization in a wild-type embryo, our modeling framework is further validated by reproducing two series of perturbation experiments, allowing further computational investigation of the network dynamics.

With the well-validated model of the C. elegans cell polarization network, here we study if it follows the balance design revealed by the exhaustive study on the antagonistic 2-node network. Interestingly, the central structure of the C. elegans network is shown to be a completely symmetric network composed of 2 nodes in the anterior (i.e., $[A]$ and $[C]$ ) and posterior (i.e., $[P]$ and $[H]$ ) respectively; this symmetric 4-node network is modified by a mutual activation in the anterior (i.e., $[A]$ ~ $[C]$ ) firstly and an additional mutual inhibition between the anterior and posterior (i.e., $[A]$ ~ $[L]$ ) subsequently, turning into an asymmetric 5-node network (Figure 5b-d, Figure 5—figure supplement 1). The abovementioned symmetric structure, the mutual activation, and the additional inhibition are in analogy with the basic 2-node network modified by self-activation and additional inhibition as shown in the top right corner of Figure 2, Figure 2—video 1 and Figure 3a.

Here, we seek to compare the spatial concentration distribution and movement velocity of the patterns generated by the three successive network structures termed ‘4-Node’ (the symmetric structure) (Figure 5b), ‘LGL-1^-’ (the symmetric structure added with $[A]$ ~ $[C]$ , that is the mutant network with $[A]$ ~ $[L]$ eliminated from the wild-type network) (Figure 5c), and ‘WT’ (the symmetric structure added with both $[A]$ ~ $[C]$ and $[A]$ ~ $[L]$ , i.e., the wild-type network; Figure 5d). As expected, both simulations on the 4-Node and WT networks successfully pass the computational pipeline described in Section (A computational pipeline for simulating cell polarization) with 62 (62 among 5,151 sets, ~1.2%) and 602 (602 among 262,701 sets, ~0.229%) viable parameter sets respectively (Figure 5e), all of which stabilize into a cell polarization pattern with the spatial concentration distribution highly intact and the movement velocity of each parameter set continuously declining below $10^{- 4}$ (Figure 5b and d, Figure 5—video 1). Nonetheless, the intermediate one (generated by depleting $[A]$ ~ $[L]$ from the WT network) with $[A]$ ~ $[C]$ but without $[A]$ ~ $[L]$ fails (Figure 5c, Figure 5—video 1), exhibiting a dispersed concentration distribution and stubbornly high movement velocity for its pattern. Such pattern instability caused by unbalanced network modification still holds when the network is initially set up using different initial conditions and parameters (i.e., $γ$ , $α$ , $k_{1}$ , $k_{2}$ , $q_{1}$ , $q_{2}$ , and $[X_{c}]$ ) assigned with various values concerning all nodes and regulations (Figure 1—figure supplement 5 – 2^nd row, Figure 5—figure supplement 4, Appendix 1—table 3), just like the simple 2-node network (Figure 1—figure supplement 4, Figure 1—figure supplement 5 – 1^st row, Appendix 1—table 2). Therefore, WT shows the largest variable parameter sets in the limited 3D grid and the strongest stability among the three network structures.

Remarkably, the balance between these two modifications is required to achieve stable cell polarization (Figure 5f). To evaluate how the interface moves, we quantify the transition plane or interface velocity by the rate of its position $x_{T}$ with respect to time. To avoid the influence of the initial state, the interface velocity is then calculated as the average velocity from $t = 300$ to $t = 500$ .

\begin{matrix} v_{I} = \frac{x_{T} (t = 500) - x_{T} (t = 300)}{200} \end{matrix}

Scanning the intensities of $[A]$ ~ $[C]$ mutual activation (i.e., $q_{A 2}^{C} & q_{C 2}^{A}$ ) and $[A]$ ~ $[L]$ mutual inhibition (i.e., $k_{A 2}^{L} & k_{L 2}^{A}$ ), a contour map is shown in consideration of the interface velocity averaged over all viable parameter sets (i.e., $v_{I}$ ): the system can be stably polarized when they are in balance with its interface velocity close to zero, but the overshoot of either intensity leads to a homogeneous state in the end (Figure 5f). The redder the contour color is, the faster the interface moves posteriorly, with aPARs invading; the bluer the contour color is, the faster the interface moves anteriorly, with pPARs invading.

Next, we wonder if such combinatorial modifications may be an optimal choice selected during evolution. To this end, we regard the $[A]$ ~ $[C]$ mutual activation as a primary modification that induces pattern asymmetry as proposed before (Seirin-Lee et al., 2020b; Lim et al., 2021), then the existence, as well as the form of the feedback loops between $[L]$ and a preexisting molecular species, is a supporting modification that consolidates stable cell polarization. Considering the symmetric structure of the LGL-1^- network, there is an identical role between $[A]$ and $[C]$ and between $[P]$ and $[H]$ , so $[L]$ can be effectively connected to $[A]$ or $[P]$ ; meanwhile, there are three types of directional regulation between them: activation, inhibition, and none. Thus, in theory, there are 32 possible network structures with an additional regulation, in which $[L]$ must interact with one of the existing nodes (Figure 5—figure supplement 5). Strikingly, only the WT network with mutual inhibition between $[A]$ and $[L]$ passes the computational pipeline with viable parameter sets. Thus, without any parametric asymmetry concessions, the configuration of the C. elegans network in nature is well optimized among all other alternatives for maintaining cell polarization pattern stability.

As pattern stability means how fast the pattern moves over time, does the lack of pattern stability result in a more dispersed concentration distribution when the system parameters fluctuate in time? To test this hypothesis, for each viable parameter set, we exert Gaussian noise on all the original values of system parameters $γ$ , $α$ , $k_{1}$ , $k_{2}$ , $q_{1}$ , and $q_{2}$ , and initiate 1000 independent simulations, where the noise amplification is represented by the standard deviation $σ$ of the Gaussian noise (Guan et al., 2021). To compare the variance of perturbed condition (abbr., PT) to the original pattern (abbr., OP, including, 4-Node, LGL-1^-, and WT), the pattern error, averaged over all molecular species, all viable parameter sets, and all independent simulations, is defined as follows:

\begin{matrix} E r r o r = \frac{1}{N M Q} \sum_{i = 1}^{Q} \sum_{j = 1}^{M} \sum_{X} \frac{\int_{- \frac{L}{2}}^{\frac{L}{2}} | {[X_{m}]}_{O P} (x, 500) - {[X_{m}]}_{P T} (x, 500) | d x}{L} \end{matrix}

where $N, M, Q$ represents the number of the molecular species types (i.e., five for WT, four for 4-Node and LGL-1^-), viable parameter sets (i.e., 602 among 262,701,~0.229%, for WT; 62 among 5,151,~1.20%, for 4-Node; 62 among 262,701,~0.024%, for LGL-1^-), and independent simulations (i.e., 1000). It turns out that the pattern error is always the smallest in the WT network and biggest in the LGL-1^- network, no matter how strong the noise is, indicating parameter robustness as the companion advantage of pattern stability (Figure 5g).

A protocol to identify responsive parameters for interface positioning

For a cell polarization pattern to be stable, fine-tuning of the kinetic parameters is required to reach the vanishing interface velocity (numerically identified by Equation 5). Using the C. elegans 5-node network as an example, we outline a method to delineate iso-velocity surfaces $v_{I} (Ω) = c o n s t a n t$ in the high-dimensional space of the parameter set $Ω$ . We show that the information gained can be used to quantify the role of individual molecular components in controlling the polarized pattern. Additionally, this knowledge enables us to design experiments to produce desired patterns.

The method of iso-velocity surface delineation is outlined as the following. In Section (Reconstruction of the molecular interaction network in C. elegans zygote and its design principle of structure and parameter trade-off), 602 (602 among 262,701, ~0.229%) groups of parameters capable of reaching stable polarity patterns are found by scanning through parameter space. When plotted in the 3D parameter space shown in Figure 5e, they span an iso-velocity surface (blue points in Figure 5e) at vanishing velocity, and in other words, their pattern interfaces are stabilized with little or no movement. To illustrate the local orientation of the surface, we selected 9,261 points (enumerated on a 3D grid $γ \in [0.034,0.044]$ in steps $∆ γ = 0.0005$ , $k_{1} \in [1.05,2.05]$ in steps ${∆ k}_{1} = 0.05$ , and $q_{2} \in [0.045,0.055]$ in steps ${∆ q}_{2} = 0.0005$ ) in the neighborhood of a benchmark point $Ω^{*} (γ^{*} = 0.039, k_{1}^{*} = 1.55, q_{2}^{*} = 0.05)$ , as shown by the orange box in the top of Figure 6a. The least-squares fit of computed interface velocity to a linear function yields:

f (γ, k_{1}, q_{2}) = - 0.0308 \times (γ - γ^{*}) - 6.08 \times 10^{- 4} \times (k_{1} - k_{1}^{*}) + 0.02 \times (q_{2} - q_{2}^{*})

Figure 6 with 1 supplement see all

Download asset Open asset

The control of the interface velocity and position by adjusting parameters in a multi-dimensional system.

(a) The parameter space (1^st row) and the linear relationship (2^nd row) between interface velocity and parameters. The discrete viable parameter space of the WT network is fitted by a blue curved surface to represent its null surface with little or no pattern interface moving (1^st row). The benchmark point $Ω^{*} (γ^{*} = 0.039, k_{1}^{*} = 1.55, q_{2}^{*} = 0.05)$ and its neighborhood are marked by an orange cube. Centering on the benchmark point $Ω^{*}$ , the relationship between the velocity interface and parameters in the region marked by the orange cube is shown by slice planes orthogonal to the $γ$ -axis at the values 0.034, 0.036, 0.038, 0.04, 0.042, and 0.044 (2^nd row). The darker the red is, the faster the interface is moving toward the posterior pole; the darker the blue is, the faster the interface is moving toward the anterior pole; the closer the color is to white, the stabler the interface is. (**b–e**) As in Figure 4, the control of the zero-velocity interface position of the polarity pattern by spatially inhomogeneous parameters corresponding to opposite interface velocities on two sides of the interface can be applied to the *C. elegans* 5-node network. (b) Using $Ω^{*}$ as a representative, spatially homogeneous parameters (2^nd column) generate a stable polarity pattern (1^st column). (c) On top of (b), a stable polarity pattern (1^st column) with its interface around $x = 0$ can be obtained by increasing $q_{C 2}^{A}$ & $q_{A 2}^{C}$ to 0.12 at $x < 0$ and increasing $γ_{P}, γ_{L}, a n d γ_{H}$ to 0.06 at $x > 0$ (2^nd column). (d) As in (c) but changing the step position to $x = - 0.1$ (2^nd column), the interface stabilizes around $x = - 0.1$ (1^st column). (e) As in (c), but changing the step position to $x = 0.1$ (2^nd column), the interface stabilizes around $x = 0.1$ (1^st column). For each parameter set, the corresponding spatial concentration distribution of $[A_{m}]$ , $[C_{m}]$ , $[P_{m}]$ , $[L_{m}]$ , and $[H_{m}]$ at $t = 0, 100, 200, 300, 400, a n d 500$ are shown beneath with a color scheme listed on the right. Note that all the interface positions in (**b–e**) (1^st column) are marked by gray dashed lines. Note that all the spatially homogeneous parameter distributions in (b) (2^nd column) are marked by gray dashed lines in (**c–e**) (2^nd column).

with the coefficient of determination:

\begin{matrix} R^{2} = 1 - \frac{\sum_{i = 1}^{9261} {[v_{I} (Ω_{i}) - f (Ω_{i})]}^{2}}{\sum_{i = 1}^{9261} {[v_{I} (Ω_{i}) - \frac{1}{9261} \sum_{j = 1}^{9261} v_{I} (Ω_{j})]}^{2}} = 0.9918 \end{matrix}

The coefficients in Equation 9 give the ascending gradient of the interface velocity (Figure 6a – bottom).

Within the assumed relationships adopted for parameter reduction (see Section (A computational pipeline for simulating cell polarization)), Equation 9 shows a strong dependence of the interface velocity on basal on/off-rates $γ$ and activation intensity $q_{2}$ , but weak sensitivity to responsive concentration $k_{1}$ in the inspected parameter region. Interestingly, the increase in the basal on/off-rates $γ$ and self-activation rate $q_{2}$ have opposite effects on the interface velocity, which can be attributed to the asymmetric 5-node network structure. In a more realistic scenario, one may consider basal on/off-rates for different molecular species and their pairwise interactions separately and explore the relationships in a much larger parameter space. A formula similar to Equation 9, using a linear expansion around the parameter sets balanced for a zero-velocity interface, enables quantitative prediction of the interface velocity against individual or simultaneous changes of parameters, even when knowledge of the cell polarization network is incomplete. In particular, the parameter-dependent interface velocity picture potentially enables biologists to manipulate and even synthesize the polarity pattern by rational modulation of parameters, controlling the interface location for specific physiological functions, such as designating a desired cell volume partition ratio (Hubatsch et al., 2019; Guan et al., 2021).

Experimental studies have revealed that the interface localization in the C. elegans zygotic cell polarization pattern provides an accurate spatial cue to regulate the dynamics of their downstream molecules, like the protein LIN-5 and microtubules that control cell division and volume partition (Schneider and Bowerman, 2003; Ajduk and Zernicka-Goetz, 2016), while the acquired cell volume has been reported to have a chain effect in cell cycle, cell position, and other cell behaviors (Fickentscher et al., 2018; Guan et al., 2021; Fickentscher et al., 2016; Tian et al., 2020). As in the simple 2-node network, the precise control of the stable interface position can be realized by taking the parameters corresponding to the opposite direction of the interface velocity on both sides of the interface. Based on the parameter set $Ω^{*}$ which gives out the stablest pattern (Figure 6b, Figure 6—video 1), two modifications are adopted to generate patterns with opposite interface velocities:

(i) Increasing the activation intensity between $[A]$ and $[C]$ ( $q_{C 2}^{A}$ & $q_{A 2}^{C}$ ) to 0.12 results in the interface traveling backward, and $[A]$ and $[C]$ finally dominate the whole domain of the cell membrane.

(ii) Increasing the basal on-rate of posterior proteins ( $γ_{P}, γ_{L}, γ_{H}$ ) to 0.06 results in the interface traveling forward, and $[P]$ , $[L]$ and $[H]$ finally dominate the whole domain of the cell membrane.

Combining the two sets of parameters with its step switch located at $x = 0$ ((i) on the left and (ii) on the right), the stabilized polarity pattern interface settles around $x = 0$ (Figure 6c, Figure 6—video 1). The interface position turns tunable as the step switch moves to $x = - 0.1$ (Figure 6d, Figure 6—video 1) and $x = 0.1$ (Figure 6e, Figure 6—video 1). One pre-known example in reality that matches this scheme is the asymmetric division of the P2 and P3 cells (the granddaughter and great-granddaughter cells of the C. elegans zygote P0), where the extracellular protein MES-1 and intracellular protein SRC-1 transduct signals from the EMS and E cells (the sister cell of P2 and its posterior daughter cell) and induce polarity reversal, possibly by affecting the local on-rate of PAR-2 (Arata et al., 2010; Seirin-Lee, 2016). Such a scheme depicts the introduction of intracellular or extracellular cues that break the spatial homogeneity of parameters and tune the stable localization of polarity pattern interface, serving as a theoretical basis for controlling oriented cell division with designated volume partition and fate differentiation (Arata et al., 2010; Hubatsch et al., 2019; Schubert et al., 2000).

Discussion

Cell polarization is a fundamental issue in both prokaryotes and eukaryotes, playing crucial roles in diverse biological phenomena ranging from chemotaxis to embryogenesis (Nance, 2014; Kondo et al., 2019). The reaction-diffusion network responsible for cell polarization has been a long-term research focus for both experimentalists and theorists, who have identified many interactive functional molecules, discovered underlying design principles for networks, and even synthesized new systems de novo (Tostevin and Howard, 2008; Chau et al., 2012; Koorman et al., 2016; Lang and Munro, 2017; Lin et al., 2021; Watson et al., 2023). In this paper, we focus on a basic problem – how to control the cell polarization pattern stability over time and regulate the stabilized interface localization, from the perspective of a reaction-diffusion network. First, we established a computational pipeline to search the viable parameter sets in an N-node (molecular species) network that can achieve a stable cell polarization pattern. The simple antagonistic network with only two nodes was revealed to be unstable (i.e., transit from a polarized distribution to a homogeneous distribution spontaneously) when any of the three unbalanced modifications (i.e., single-sided self-regulation, single-sided additional regulation, and unequal system parameters) were introduced. To recover stable cell polarization, two strategies are proposed: (i) the combination of two unbalanced modifications with opposite effects; (ii) the combination of two spatially inhomogeneous parameter values, either of which can lead to opposite interface velocity. Additionally, the stable interface localization can be discretionarily regulated by the step-like parameter profile and contributes to the asymmetric geometry of the cell polarization pattern, potentially providing a spatial cue for significant physiological functions like unequal cell volume partition and consequent punctuated cell movement (Fickentscher et al., 2016; Fickentscher et al., 2018; Tian et al., 2020; Guan et al., 2021). Analogous conclusions are further identified in the experimentally summarized C. elegans network which obtains large parameter space and high robustness against parameter perturbations, supporting them as strategies indeed applied in real biological scenarios. Importantly, the linear relationship between the interface velocity and biophysical parameters serves as a useful tool to characterize and predict how the cell polarization pattern, especially where the interface is located or moved toward, responds to parameter changes; this not only helps understand how the regulatory molecular species and pathways, as well as the biophysical parameters learned from experiments, play their role, but also guides the design of artificial cell polarization systems with desired features. The joint study on a simple 2-node network and C. elegans 5-node network demonstrated that the cell polarization pattern with both stability and asymmetry can be explicitly realized by a combinatorial network and spatially inhomogeneous parameters, which is expected to facilitate the interpretation of natural systems and design of artificial systems.

Although we deciphered the stability and asymmetry of the pattern generated by the C. elegans cell polarization network, additional experimental details are required to achieve a more comprehensive understanding of this system, as outlined below.

Although the model in this paper has tried to incorporate as many molecular species as possible for the maintenance phase of cell polarization, the actual network may include yet unidentified molecular species. This also extends to interactions, as well as their intensity, between molecular species. For instance, the components of the molecular complex (e.g., PAR-3/PAR-6/PKC-3), often treated as a single entity in theoretical studies (Tostevin and Howard, 2008; Gross et al., 2019; Seirin-Lee et al., 2020b; Lim et al., 2021), might not be ideally identical concerning their slightly different polarized patterns (Wang et al., 2017). The PAR-6 and PKC-3 also bind CDC-42, in addition to PAR-3, to mediate membrane association of CHIN-1 via phosphorylation, and such joint effect was previously described to be governed by a single node, CDC-42 (Lim et al., 2021). Besides, CDC-42 may inhibit CHIN-1 not through direct interaction but indirectly by recruiting the PAR-6/PKC-3 complex. Future experimental validation is needed to detail whether and how the molecular species described in Supplementary file 1 (and beyond) interact with each other.
The large number of molecular species introduces substantial parameter complexity, with myriad possible value assignments. Simplifying the model is necessary to manage computational costs unless all parameters can be measured accurately in vivo (Wang and You, 2020).
Quantitative biological details, such as the spatial concentration distributions, diffusions, productions, and degradations of molecular species on the cell membrane versus in the cell cytosol (Tostevin and Howard, 2008; Goehring et al., 2011b; Kravtsova and Dawes, 2014; Sailer et al., 2015; Seirin-Lee and Shibata, 2015; Seirin-Lee, 2016; Gross et al., 2019; Seirin-Lee et al., 2020c; Lim et al., 2021) and the precise mathematical forms of regulatory interactions (e.g., Hill equations or product equations) (Goehring et al., 2011b; Seirin-Lee et al., 2020b; Lim et al., 2021) challenge model assumptions and hinder the unification of parameter value with unit between theory and experiment (Tostevin and Howard, 2008; Kravtsova and Dawes, 2014). Nonetheless, prior studies integrating theory and experiments have demonstrated that unveiling mathematical principles for regulatory networks does not strictly depend on unit unification; numerical network structure and parameter analysis can still effectively draw the mathematical principles for guiding the synthesis of artificial networks in reality, including the one for cell polarization (Elowitz and Leibler, 2000; Gardner et al., 2000; Ma et al., 2009; Chau et al., 2012).
Incorporating biological processes across temporal stages and spatial scales could offer a more holistic understanding of cell polarization. On one hand, early stages, such as the sperm entry and establishment phase, which involve active actomyosin contractility and flow driving PAR-3/PAR-6/PKC-3 to be localized in the anterior, can provide insights into how initial polarized concentration distributions are achieved (Cuenca et al., 2003; Motegi et al., 2011). On the other hand, molecular (e.g., downstream localization of the myosin motor NMY-2 to the anterior cortex mediated by CDC-42 and PAR-2 Goehring et al., 2011a; Rose and Gönczy, 2014; Geßele et al., 2020) and cellular (asymmetric segregation of cell volume and molecular content, mediated by molecular species like MES-1 and SRC-1 Arata et al., 2010; Seirin-Lee, 2016) effects of cell polarization can further elucidate its various roles in embryogenesis.

Beyond the C. elegans zygote studied in this paper, cell polarization also exists in later stages of C. elegans embryogenesis as well as in other organisms, where diverse functional dynamics have been reported (Knoblich, 2001; Rose and Gönczy, 2014). Governed by the cell polarization, the C. elegans embryo actually proceeds through four rounds of asymmetric cell divisions to produce four somatic founder cells. Then, such a reaction-diffusion pattern will lose its sharp transition plane or interface as it is unscalable over cell size, leading to symmetric cell division at last (Hubatsch et al., 2019). Furthermore, as the reaction-diffusion network itself is not particularly constricted in the cellular scale, but also possibly in subcellular and multicellular scales, it would be interesting to investigate if the general principles proposed in this paper are also at work for the network programming and pattern modulation in other scales (Chao et al., 2014; Guan et al., 2022).

Nowadays, the de novo construction of artificial cells with designated new functions is emerging, demanding theoretical guidance about how to synthesize the molecular control circuits beneath (Zhou et al., 2023; Zhu et al., 2023). Like toggle switch, oscillation, and adaptation (Gardner et al., 2000; Sun et al., 2022; Zhou et al., 2023), the polarized concentration distribution of molecular species (accounting for cell polarization) is also another elementary behavior in a cell, and the corresponding circuits have been synthesized successfully over a decade ago (Chau et al., 2012). However, it remains unclear whether a molecular species’ distribution pattern (e.g., the transition plane or interface) can be altered while maintaining stability. Taking the C. elegans network as an example, our computational framework has been demonstrated to be capable of deciphering such an experimentally summarized reaction-diffusion network with pattern stability and asymmetry and has been packed as a user-friendly software, PolarSim (Appendix 2, Source code 1); thus we believe it can be used for not only understanding the natural molecular interaction networks reported before but also designing new molecular interaction networks in the polarized cells with more physiological functions, for example, with tunable transition plane or interface of molecular species to guide the allocation of downstream fate determinants and unequal cleavage during cytokinesis (Fickentscher and Weiss, 2017). In the future, our computational framework for cell polarization networks could be linked to the one about cytoskeleton activity, achieving a more comprehensive modeling description for cell division (Pavin et al., 2012; Ma et al., 2014).

Appendix 1

Parameter nondimensionalization and order of magnitude consistency

The parameter nondimensionalization of the reaction-diffusion model is implemented according to a previous theoretical study (Seirin-Lee et al., 2020b). The dimensional evolution equation (corresponding to the nondimensional one in Equation 1) can be written as:

\frac{\partial [\bar{X_{m}}]}{\partial \bar{t}} = \bar{D_{m}^{X}} \frac{\partial^{2}}{\partial {\bar{x}}^{2}} \bar{[X_{m}]} + {\bar{γ_{X}} + \sum_{Y} \frac{\bar{q_{Y 2}^{X}} \bar{{[Y_{m}]}^{2}}}{1 + \bar{q_{Y 1}^{X}} {\bar{[Y_{m}]}}^{2}}} \bar{[X_{c}]} - {\bar{α_{X}} + \sum_{Y} \frac{\bar{k_{Y 2}^{X}} \bar{{[Y_{m}]}^{2}}}{1 + \bar{k_{Y 1}^{X}} \bar{{[Y_{m}]}^{2}}}} \bar{[X_{m}]}

where the dimensions (or units) of these original parameters are $\bar{t} \sim s$ , $\bar{x} ~ μ m$ , $\bar{[X_{m}]} \sim \bar{[X_{c}]} \sim μ m^{- 1}$ , $\bar{D_{X}} \sim μ m^{2} \cdot s^{- 1}$ , $\bar{γ_{X}} \sim \bar{α_{X}} \sim s^{- 1}$ , $\bar{q_{Y 1}^{X}} \sim \bar{k_{Y 1}^{X}} \sim μ m^{2}$ , and $\bar{q_{Y 2}^{X}} \sim \bar{k_{Y 2}^{X}} \sim μ m^{2} \cdot s^{- 1}$ . Then we let $t = \frac{\bar{t}}{τ}$ , $x = \frac{\bar{x}}{x_{0}}$ , $[X_{m}] = \frac{\bar{[X_{m}]}}{X_{0}}$ , $[X_{c}] = \frac{\bar{[X_{c}]}}{X_{0}}$ , and $[Y_{m}] = \frac{\bar{[Y_{m}]}}{Y_{0}}$ , where $τ$ , $x_{0}$ , $X_{0}$ , and $Y_{0}$ are arbitrary constants with dimesions as $s$ , $μ m$ , ${μ m}^{- 1}$ , and ${μ m}^{- 1}$ correspondingly. In this paper, they adopt values from literature (Seirin-Lee et al., 2020b): $τ = 2 s$ and $x_{0} = 142 μ m$ with $X_{0} = Y_{0} = 1 {μ m}^{- 1}$ (as a numerically normalized value). Finally, we derive the nondimensional evolution equation:

\frac{\partial [X_{m}]}{\partial t} = \frac{\bar{D_{m}^{X}} τ}{x_{0}^{2}} \frac{\partial^{2} [X_{m}]}{\partial x^{2}} + {\bar{γ_{X}} τ + \sum_{Y} \frac{\bar{q_{Y 2}^{X}} Y_{0}^{2} τ {[Y_{m}]}^{2}}{1 + \bar{q_{Y 1}^{X}} Y_{0}^{2} {[Y_{m}]}^{2}}} [X_{c}] - {\bar{α_{X}} τ + \sum_{Y} \frac{\bar{k_{Y 2}^{X}} Y_{0}^{2} τ {[Y_{m}]}^{2}}{1 + \bar{k_{Y 1}^{X}} Y_{0}^{2} {[Y_{m}]}^{2}}} [X_{m}]

where the nondimensional parameters are expressed by the original ones as the following: $D_{m}^{X} = \frac{\bar{D_{m}^{X}} τ}{x_{0}^{2}}$ , $γ_{X} = \bar{γ_{X}} τ$ , $α_{X} = \bar{α_{X}} τ$ , $q_{Y 1}^{X} = \bar{q_{Y 1}^{X}} Y_{0}^{2}$ , $q_{Y 2}^{X} = \bar{q_{Y 2}^{X}} Y_{0}^{2} τ$ , $k_{Y 1}^{X} = \bar{k_{Y 1}^{X}} Y_{0}^{2}$ , $k_{Y 2}^{X} = \bar{k_{Y 2}^{X}} Y_{0}^{2} τ$ . Equation A1 is equivalent to Equation 1 in Section (A computational pipeline for simulating cell polarization).

Next, we reveal the consistency in order of magnitude between numerical and realistic values using four examples.

Evolution time: With a time step $τ = 2 s$ , the evolution time equals to $τ = 500 s t e p s \times 2 s p e r s t e p = 16 m i n 40 s$ , in the same order of magnitude as the period of maintenance phase of C. elegans zygotic cell polarization, ~10 min, measured in vivo (Blanchoud et al., 2015).
Membrane diffusion coefficients: This was done in a previous theoretical study (Seirin-Lee et al., 2020b), based on experimental measurements (Goehring et al., 2011a; Goehring et al., 2011b; Geßele et al., 2020). Given the lengths of the major and minor semi-axis of the ellipsoidal C. elegans zygote are roughly $a = 27.0 μ m$ and $b = 14.8 μ m$ respectively, the perimeter of an ellipse after dimension reduction is roughly $x_{0} = 2 π b + 4 (a - b) \approx 142 μ m$ . As the membrane diffusion rate of anterior PAR protein (exemplified by PAR-6) is $0.28 μ m^{2} \cdot s^{- 1}$ , while that of posterior PAR protein (exemplified by PAR-2) is $0.15 {μ m}^{2} \cdot s^{- 1}$ , their nondimensional values are determined as $D_{m}^{X} = \frac{\bar{D_{m}^{X}} τ}{x_{0}^{2}} = \frac{0.28 \times 2}{142^{2}} \approx 2.748 \times 10^{- 5}$ ( $X \in [A, C]$ ) and $D_{m}^{X} = \frac{\bar{D_{m}^{X}} τ}{x_{0}^{2}} = \frac{0.15 \times 2}{142^{2}} \approx 1.472 \times 10^{- 5}$ ( $X \in [P, L, H]$ ) respectively, which are used in this paper.
Basal off-rate: As the basal off-rate of anterior PAR protein (exemplified by PAR-6) is approximately $5.4 \times 10^{- 3} s^{- 1}$ , while that of posterior PAR protein (exemplified by PAR-2) is approxiamtely $7.3 \times 10^{- 3} s^{- 1}$ (Goehring et al., 2011a; Goehring et al., 2011b), their nondimensional values are determined as $α_{X} = \bar{α_{X}} τ = 5.4 \times 10^{- 3} \times 2 = 0.0104$ ( $X \in [A, C]$ ) and $α_{X} = \bar{α_{X}} τ = 7.3 \times 10^{- 3} \times 2 = 0.0146$ ( $X \in [P, L, H]$ ) respectively, which match their scanned region $[0, 0.05]$ in this paper.
Inhibition intensity: As PAR-3/PAR-6/PKC-3 can be phosphorylated by PAR-2 with the rate of $0.02 \sim 0.8 μ m^{2} \cdot s^{- 1}$ on the membrane, while PAR-1/PAR-2 can be released from the membrane by phosphorylation of PKC-3 with the rate of $0.06 \sim 1.6 {μ m}^{2} \cdot s^{- 1}$ (Geßele et al., 2020; Goehring et al., 2011b), their nondimensional values are determined as between $k_{Y 2 (M I N)}^{X} = \bar{k_{Y 2 (M I N)}^{X}} Y_{0}^{2} τ = 0.02 \times 1^{2} \times 2 = 0.04$ and $k_{Y 2 (M A X)}^{X} = \bar{k_{Y 2 (M A X)}^{X}} Y_{0}^{2} τ = 1.6 \times 1^{2} \times 2 = 3.2$ , which match their scanned region $[0, 5]$ in this paper.

Appendix 1—table 1

The parameter description of the reaction-diffusion model for simulating a cell polarization network.

Parameter	Biological implication	Defaulted nondimensionalvalue in simulation
$x$	Position along the anterior-posterior axis	[–0.25,0.25] in steps 0.005
$t$	Time	[0,500] in steps 1 (corresponding to 2 sec per step) ( Blanchoud et al., 2015; Seirin-Lee et al., 2020b)
$[X]$	Molecular species	/
$[X_{m}]$	Spatial concentration distribution on the cell membrane	/
$[X_{c}]$	Spatial concentration distribution in the cell cytosol	1
$D_{m}^{X}$	Diffusion coefficient across the cell membrane	2.748×10^–5 for anteriorly-located molecular species and 1.472×10^–5 for posteriorly-located molecular species (Goehring et al., 2011a; Goehring et al., 2011b; Geßele et al., 2020; Seirin-Lee et al., 2020b)
$γ_{X}$	Basal on-rate (cell membrane association)	[0,0.05] in steps 0.001
$α_{X}$	Basal off-rate (cell membrane dissociation)	[0,0.05] in steps 0.001
$q_{Y 1}^{X}$	Responsive concentration in the activation pathway from $[Y]$ to $[X]$	[0,5] in steps 0.05
$k_{Y 1}^{X}$	Responsive concentration in the inhibition pathway from $[Y]$ to $[X]$	[0,5] in steps 0.05
$q_{Y 2}^{X}$	Activation intensity from $[Y]$ to $[X]$	[0,0.05] in steps 0.001
$k_{Y 2}^{X}$	Inhibition intensity from $[Y]$ to $[X]$	1
$Ω$	A three-dimensional (3D) grid representing nondimensionalized parameter set	$(γ, k_{1}, q_{2})$
$n$	Hill coefficient	2 (Seirin-Lee and Shibata, 2015; Seirin-Lee, 2021)
$L$	Cell/embryo length	0.5
$v_{M}$	Pattern movement velocity	/
$x_{T}$	Pattern transition plane or interface	/
$v_{I}$	Interface velocity	/
$N$	Number of molecular species types	/
$M$	Number of viable parameter sets	/
$Q$	Number of independent simulations	/

Appendix 1—table 2

Viable parameter set searched by Monte Carlo method for the simple 2-node network.

For the 10 parameters in the simple 2-node network, the corresponding scanned ranges are $γ \in [0, 0.05], α \in [0, 0.05], k_{1} \in [0.05, 5], k_{2} \in [0.05, 5]$ , and . Among 100,000 independent simulations, only 20 viable parameter sets pass the computational pipeline. One representative viable parameter set used to reproduce the conclusions in Figure 1; Figure 1—figure supplement 5 -1^st row is listed below.

Parameter	NondimensionalValue in Simulation	Parameter	NondimensionalValue in Simulation
$[A_{c}]$	1.5361	$[P_{c}]$	0.6662
$γ_{A}$	0.0251	$α_{A}$	0.0437
$γ_{P}$	0.0312	$α_{P}$	0.0245
$k_{P 1}^{A}$	0.1995	$k_{P 2}^{A}$	3.8471
$k_{A 1}^{P}$	1.3638	$k_{A 2}^{P}$	2.2571

Appendix 1—table 3

Viable parameter set searched by Monte Carlo method for the C. elegans 5-node network.

For the 39 parameters in the C. elegans 5-node network, the corresponding scanned ranges are $γ \in [0, 0.05], α \in [0, 0.05], k_{1} \in [0.05, 5], k_{2} \in [0.05, 5], q_{1} \in [0.05, 5], q_{2} \in [0, 0.05]$ , and $[X_{c}] \in [0.05, 5]$ . Among 100,000 independent simulations, only 6 viable parameter sets pass the computational pipeline. One representative viable parameter set used to reproduce the conclusions in Figure 1; Figure 1—figure supplement 5 -2^nd row is listed below.

Parameter	NondimensionalValue in Simulation	Parameter	NondimensionalValue in Simulation
$[A_{c}]$	2.4420	$[P_{c}]$	3.6739
$[C_{c}]$	1.8330	$[L_{c}]$	2.3793
$[H_{c}]$	3.4683
$γ_{A}$	0.0448	$α_{A}$	0.0310
$γ_{P}$	0.0164	$α_{P}$	0.0334
$γ_{C}$	0.0158	$α_{C}$	0.0071
$γ_{L}$	0.0146	$α_{L}$	0.0221
$γ_{H}$	0.0197	$α_{H}$	0.0410
$q_{A 1}^{C}$	4.1521	$q_{A 2}^{C}$	0.0495
$q_{C 1}^{A}$	1.3029	$q_{C 2}^{A}$	0.0201
$k_{P 1}^{A}$	3.7955	$k_{P 2}^{A}$	0.2009
$k_{L 1}^{A}$	2.3597	$k_{L 2}^{A}$	3.7357
$k_{H 1}^{A}$	1.9557	$k_{H 2}^{A}$	2.2441
$k_{A 1}^{P}$	3.7180	$k_{A 2}^{P}$	3.7342
$k_{C 1}^{P}$	4.5544	$k_{C 2}^{P}$	1.5792
$k_{P 1}^{C}$	0.5495	$k_{P 2}^{C}$	2.8838
$k_{H 1}^{C}$	2.6518	$k_{H 2}^{C}$	1.2697
$k_{A 1}^{L}$	2.4436	$k_{A 2}^{L}$	4.8334
$k_{A 1}^{H}$	4.4505	$k_{A 2}^{H}$	0.3688
$k_{C 1}^{H}$	3.3712	$k_{C 2}^{H}$	2.0272

Appendix 2

Instructions for PolarSim

2.1. Introduction

PolarSim is a graphical user interface (GUI) on Matlab 2022b (The MathWork Inc, 2022) for simulating the evolution of cell polarization patterns. Based on the reaction-diffusion model, the GUI allows users to compute the behaviors of cell polarization networks in different biological scenarios. All the simulations are tested with a 12th Gen Intel(R) Core(TM) i7-1260P CPU.

2.2. Tutorials

Download the folder “PolarSim” from https://github.com/YixuanChen0726/Cell-Polarization/tree/main/PolarSim.
Open MATLAB under the ‘PolarSim’ folder path and execute script ‘GUI.m’. Click ‘Run’ and then an interactive interface pops up (Appendix 2—figure 1).
With the following parameters inputted, the GUI gives out a group of ‘Pattern_*.mat’ files containing pattern information.

Appendix 2—figure 1

Download asset Open asset

The instructions to open the *PolarSim* GUI.

① The files in the folder ‘PolarSim’. ② Open *Matlab* under the path of the folder ‘PolarSim’ and double-click to open ‘GUI.m’. ③ Click ‘Run’ to open the *PolarSim* GUI shown by ④.

Appendix 2—table 1

The instructions for the format of Network Parameter in an Excel table.

In the 1^st row and 1^st column, ‘Node N’ represents the name of the node. The 2^nd row explains the characteristic parameters of each node (i.e., location, cytoplasmic concentration, basal on-rate, and basal off-rate as listed from left to right). The interaction parameters start from the 3^rd row and 2^nd column, where the i^th row and j^th column describe the activation/inhibition effect from Node j to Node i. Note that $q_{Y 2}^{X}$ and $k_{Y 2}^{X}$ should be set to 0 when no activation or inhibition is exerted on X from Y respectively. ‘Location’ should be assigned with the string ‘a’ or ‘p’ while the description of the other parameters is detailed in Appendix 1—table 1.

	Node A	Node B	…	Node N
	Location (a or p)/ $[A_{c}]$ / $γ_{A}$ / $α_{A}$	Location (a or p)/ $[B_{c}]$ / $γ_{B}$ / $α_{B}$	…	Location (a or p)/ $[N_{c}]$ / $γ_{N}$ / $α_{N}$
Node A	$q_{A 1}^{A}$ / $q_{A 2}^{A}$ / $n_{A q}^{A}$ / $k_{A 1}^{A}$ / $k_{A 2}^{A}$ / $n_{A k}^{A}$	$q_{B 1}^{A}$ / $q_{B 2}^{A}$ / $n_{B q}^{A}$ / $k_{B 1}^{A}$ / $k_{B 2}^{A}$ / $n_{B k}^{A}$	…	$q_{N 1}^{A}$ / $q_{N 2}^{A}$ / $n_{N q}^{A}$ / $k_{N 1}^{A}$ / $k_{N 2}^{A}$ / $n_{N k}^{A}$
Node B	$q_{A 1}^{B}$ / $q_{A 2}^{B}$ / $n_{A q}^{B}$ / $k_{A 1}^{B}$ / $k_{A 2}^{B}$ / $n_{A k}^{B}$	$q_{B 1}^{B}$ / $q_{B 2}^{B}$ / $n_{B q}^{B}$ / $k_{B 1}^{B}$ / $k_{B 2}^{B}$ / $n_{B k}^{B}$	…	$q_{N 1}^{B}$ / $q_{N 2}^{B}$ / $n_{N q}^{B}$ / $k_{N 1}^{B}$ / $k_{N 2}^{B}$ / $n_{N k}^{B}$
…	…	…	…	…
Node N	$q_{A 1}^{N}$ / $q_{A 2}^{N}$ / $n_{A q}^{N}$ / $k_{A 1}^{N}$ / $k_{A 2}^{N}$ / $n_{A k}^{N}$	$q_{B 1}^{N}$ / $q_{B 2}^{N}$ / $n_{B q}^{N}$ / $k_{B 1}^{N}$ / $k_{B 2}^{N}$ / $n_{B k}^{N}$	…	$q_{N 1}^{N}$ / $q_{N 2}^{N}$ / $n_{N q}^{N}$ / $k_{N 1}^{N}$ / $k_{N 2}^{N}$ / $n_{N k}^{N}$

Appendix 2—figure 2

Download asset Open asset

The examples of network parameters.

(a) ‘Simple Antagonistic 2-Node Network.xlsx lists the parameters’ in Figure 1a. (b) ‘C. elegans Wild-Type 5-Node Network.xlsx’ lists the parameters for the benchmark point $Ω^{*} (γ^{*} = 0.039, k_{1}^{*} = 1.55, q_{2}^{*} = 0.05)$ in Figure 6b. (c) ‘LGL-1 Mutant 4-Node Network.xlsx’ lists the parameters as in (b) but with the Node $[L]$ knocked out.

(i) Import: Network Parameter.

We give out three examples ‘Simple Antagonistic 2_Node Network.xlsx’, ‘LGL-1 Mutant 4-Node Network.xlsx’ and ‘C. elegans Wild-Type 5-Node Network.xlsx’ in the folder ‘PolarSim’, respectively representing the typical networks in this paper (Appendix 2—figure 2). The Excel table for parameter value assignments should follow the format in Appendix 2—table 1.
(ii) Import: Cell Length

We take ‘0.5’ as an example. Any positive number is allowed in this box. The effects of cell length on cell polarization patterns are shown in Appendix 2—figure 6.
(iii) Import: Output Folder

Give a folder name for storing the output results (e.g., ‘Output 2-Node’)
(iv) Setting: Time for Simulation

Simulation duration ‘500’ is used in this paper. Any positive number is allowed in this box.
(v) Setting: Time Step Length

We take ‘1’ as an example. Larger values are not recommended, in consideration of possible error accumulation, and one may try a smaller step length while the error tolerance is exceeded.
(vi) Setting: Saving Interval

It must be an integer multiple of the time step length. The time point will be saved in this designated interval and can be used for pattern plotting later.

Click ‘Run’ in the interface, and then its status on the progress bar is shown (Appendix 2—figure 3a). A folder named ‘Import: Output Folder’ is generated in the current path to store the output ‘Pattern_*.mat’ data containing the molecular species or node’s name, location, and concentration distribution on the cell membrane, where ‘*’ denotes the in silico time corresponding to each file (Appendix 2—figure 4a - right).
Import the pathway of the file outputted by PolarSim into the box ‘Import: Pattern’.
Click ‘Plot’ and then a figure pops up to show the cell polarization pattern while the position of the transition plane appears in the box ‘Output: Transition Plane’.
Import the pathway of two files outputted by PolarSim into the box ‘Import: Start Time’ and ‘Import: End Time’.
Click ‘Run’ and then the mean interface velocity between two input time points appears in the box ‘Output: Interface Velocity’.

Appendix 2—figure 3

Download asset Open asset

The progress bar, files, and subfolders of *PolarSim*.

(a) The progress bar showing the running progress. (b) All files and output subfolders in the folder ‘PolarSim’.

Appendix 2—figure 4

Download asset Open asset

The results of Example 1 (‘Simple Antagonistic 2-Node Network.xlsx’) and Example 2 (‘C. elegans Wild-Type 5-Node Network.xlsx’) of *PolarSim*.

(a) The flow chart for computing Example 1: ① input parameters; ② click ‘Run’; ③ the simulation is completed with a progress bar shown and files saved in the folder ‘Output 2-Node’; the file ‘Pattern_500.mat’ is used to show the data format on the right, where the first part stores the molecular species or nodes’ names and locations while the second part stores their spatial concentration distributions on the cell membrane; ④ input the pathway of the outputted pattern file; ⑤ click ‘Plot’; ⑥ the simulation is completed with a figure shown and the position of the transition plane given in the box ‘Output: Transition Plane’; ⑦ input the pathway of the two outputted files ‘Pattern_*.mat’ at different time points into ‘Import: Start Time’ and ‘Import: End Time’, where ‘*’ denotes the start time and end time respectively; ⑧ click ‘Run’; ⑨ the simulation is completed with the interface velocity given in the box ‘Output: Interface Velocity’. (b) The same as a but for the ‘C. elegans Wild-Type 5-Node Network.xlsx’.

2.3. Examples

Example 1: ‘Simple Antagonistic 2-Node Network.xlsx’ simulation as Figure 1a.

Input the following parameters: ‘Simple Antagonistic 2-Node Network.xlsx’ into ‘Import: Network Parameter’; ‘0.5’ into ‘Import: Cell Length’; ‘Output 2-Node’ into ‘Import: Output Folder’; ‘500’ into ‘Setting: Time for Simulation’; ‘1’ into ‘Setting: Time Step Length’ and ‘100’ into ‘Setting: Saving Interval’. (Appendix 2—figure 4a– ①)
Click ‘Run’, and then six ‘Pattern_*.mat’ files at different time points are saved in the subfolder ‘Output 2-Node’ (Appendix 2—figure 4a– ②-③).
Input ‘Output 2-Node\Pattern_500.mat’ into the box ‘Import: Pattern’, and a figure at $t = 500$ appears. Then the position of the transition plane $x_{T} = 0.00252525$ is given in the box ‘Output: Transition Plane’ (Appendix 2—figure 4a – ④-⑥).
Input ‘Output 2-Node\Pattern_300.mat’ and ‘Output 2-Node\Pattern_500.mat’ into the box ‘Import: Start Time’ and ‘Import: End Time’ respectively. Then the interface velocity between and is given in the box ‘Output: Interface Velocity’ (Appendix 2—figure 4a– ⑦-⑨). The 2-node network approaches a stable polarized pattern with its interface velocity being 0.

Example 2: ‘C. elegans Wild-Type 5-Node Network.xlsx’ simulation as Figure 6b.

Input the following parameters: ‘C. elegans Wild-Type 5-Node Network.xlsx’ into ‘Import: Network Parameter’; ‘0.5’ into ‘Import: Cell Length’; ‘Output 5-Node’ into ‘Import: Output Folder’; ‘500’ into ‘Setting: Time for Simulation’; ‘1’ into ‘Setting: Time Step Length’ and ‘100’ into ‘Setting: Saving Interval’ (Appendix 2—figure 4b– ①).
Click ‘Run’”, and then six ‘Pattern_*.mat’ files at different time points are saved in the subfolder ‘Output 5-Node’ (Appendix 2—figure 4b– ②-③).
Input ‘Output 5-Node\Pattern_500.mat’ into the box ‘Import: Pattern’”, and a figure at $t = 500$ appears. Then the position of the transition plane $x_{T} = - 0.0176768$ is given in the box ‘Output: Transition Plane’ (Appendix 2—figure 4b– ④-⑥).
Input ‘Output 5-Node\Pattern_300.mat’ and ‘Output 5-Node\Pattern_500.mat’ in the box ‘Import: Start Time’ and ‘Import: End Time’, respectively. Then the interface velocity between and is given in the box “Output: Interface Velocity” (Appendix 2—figure 4b – ⑦-⑨). The 5-node network approaches a stable polarized pattern with its interface velocity being 0.

Example 3: ‘LGL-1 Mutant 4-Node Network.xlsx’ simulation, originated from ‘C. elegans Wild-Type 5-Node Network.xlsx’ but with the Node $[L]$ knocked out.

Input the following parameters: ‘LGL-1 Mutant 4-Node Network.xlsx’ into ‘Import: Network Parameter’; ‘0.5’ into ‘Import: Cell Length’; ‘Output 4-Node’ into ‘Import: Output Folder’; ‘1000’ into ‘Setting: Time for Simulation’; ‘1’ into ‘Setting: Time Step Length’ and ‘100’ in ‘Setting: Saving Interval’ (Appendix 2—figure 5a).
Click ‘Run’, and then 11 ‘Pattern_*.mat’ files at different time points are saved in the subfolder ‘Output 4-Node’.
Input ‘Output 4-Node\Pattern_500.mat’ into the box ‘Import: Pattern’”, and then a figure of LGL-1 Mutant 4-Node Network at appears with the transition plane close to the posterior pole. Then the position of transition plane is given in the box ‘Output: Transition Plane’ (Appendix 2—figure 5a).
Input ‘Output 4-Node\Pattern_300.mat’ and ‘Output 4-Node\Pattern_500.mat’ in the box ‘Import: Start Time’ and ‘Import: End Time’ respectively. Then the interface velocity between $t = 300$ and $t = 500$ is given in the box ‘Output: Interface Velocity’ (Appendix 2—figure 5a). The interface of the 4-node network keeps moving toward the posterior with $v_{I} = 0.000429293$ .
Input ‘Output 4-Node\Pattern_1000.mat’ in the box ‘Import: Pattern’, and then a figure at appears. The string ‘Transition plane doesn’t exist’ appears in the box ‘Output: Transition Plane’ as the pattern reaches a homogeneous state dominated by and (Appendix 2—figure 5b).
Input ‘Output 4-Node\Pattern_800.mat’ and ‘Output 4-Node\Pattern_1000.mat’ in the box ‘Import: Start Time’ and ‘Import: End Time’, respectively. Then the interface velocity between and doesn’t exist with the string ‘It isn’t a polarized pattern’ appearing in the box ‘Output: Interface Velocity’ to hint (Appendix 2—figure 5b). Note that the interface velocity can’t be calculated when either pattern at start time or end time is homogeneous.

Appendix 2—figure 5

Download asset Open asset

The results of Example 3 (‘LGL-1 Mutant 4-Node Network.xlsx’) of *PolarSim*.

(a) The transition plane is shown at $t = 500$ , with the value $x_{T} = 0.207071$ . The interface velocity is calculated between $t = 300$ and $t = 500$ , with the value $v_{I} = 0.000429293$ representing an unstable pattern (top). The figure is plotted at $t = 500$ (bottom). (b) The transition plane and the figure are shown at $t = 1000$ and the interface velocity is calculated between $t = 800$ and $t = 1000$ . The pattern collapses to a homogeneous state with $[A]$ and $[C]$ invading the posterior domain at $t = 1000$ , and thereby the transition plane doesn’t exist and the interface can’t be calculated.

2.4. Extensive application

Our PolarSim is extensively applicable to similar biological systems. Here, we take the cell size (length) as an exemplary research target to study how the concentration distribution on the cell membrane depends on it. Different cell lengths are applied to the Example 1 (‘Simple Antagonistic 2-Node Network.xlsx’) and Example 2 (‘C. elegans Wild-Type 5-Node Network.xlsx’) to see whether there is a cell size threshold limiting cell polarization as discovered before (Hubatsch et al., 2019). Patterns at $t = 500$ are plotted with the cell length ranging from 0.1 to 0.5 in steps 0.1 (Appendix 2—figure 6).

Appendix 2—figure 6

Download asset Open asset

The effects of cell size (length) on the cell polarization pattern.

(**a–e**) The pattern of ‘Simple Antagonistic 2-Node Network.xlsx’ at $t = 500$ . From left to right, the cell lengths are 0.1, 0.2, 0.3, 0.4 and 0.5, respectively. (**f–j**) The same as (**a–e**) but for the ‘C. elegans Wild-Type 5-Node Network.xlsx’.

The PolarSim-based simulations above indicate that the proper cell polarization demands a reasonable spatial scale, which to some extent gives an explicit constraint for the volume of a cell in reality when it needs to divide asymmetrically by reading out its cell polarization pattern interface (Hubatsch et al., 2019). It’s worth noting that our prediction of cell polarization pattern collapse (with less and less distinguishable interface) over cell size decrease aligns with previous experimental and theoretical research (Hubatsch et al., 2019). In all, PolarSim provides a user-friendly tool for more applications on the studies in cell polarization.

2.5. Contact

All the scripts of the PolarSim GUI have been uploaded onto GitHub https://github.com/YixuanChen0726/Cell-Polarization/tree/main/PolarSim. If there is any question, please contact Yixuan Chen (yixuanchen@stu.pku.edu.cn) or Guoye Guan (guanguoye@gmail.com) anytime.

Data availability

The current manuscript is a theoretical and computational study, so no data has been generated for this manuscript. The computational pipeline for simulating cell polarization is packed as a user-friendly software, PolarSim, all associated source code is available on GitHub at https://github.com/YixuanChen0726/Cell-Polarization/tree/main/PolarSim (copy archived at Chen, 2025; also see Appendix 2 and Source code 1).

References

(2006) Interaction of PAR-6 with CDC-42 is required for maintenance but not establishment of PAR asymmetry in C. elegans
Developmental Biology 299:386–397.

https://doi.org/10.1016/j.ydbio.2006.08.002
- PubMed
- Google Scholar
1. Ajduk A
2. Zernicka-Goetz M
(2016) Polarity and cell division orientation in the cleavage embryo: from worm to human
Molecular Human Reproduction 22:691–703.

https://doi.org/10.1093/molehr/gav068
- PubMed
- Google Scholar
1. Arata Y
2. Lee J-Y
3. Goldstein B
4. Sawa H
(2010) Extracellular control of PAR protein localization during asymmetric cell division in the C. elegans embryo
Development 137:3337–3345.

https://doi.org/10.1242/dev.054742
- PubMed
- Google Scholar
(2010) The C. elegans homolog of Drosophila Lethal giant larvae functions redundantly with PAR-2 to maintain polarity in the early embryo
Development 137:3995–4004.

https://doi.org/10.1242/dev.056028
- PubMed
- Google Scholar
(2013) PAR-2, LGL-1 and the CDC-42 GAP CHIN-1 act in distinct pathways to maintain polarity in the C. elegans embryo
Development 140:2005–2014.

https://doi.org/10.1242/dev.088310
- PubMed
- Google Scholar
1. Bilder D
2. Li M
3. Perrimon N
(2000) Cooperative regulation of cell polarity and growth by Drosophila tumor suppressors
Science 289:113–116.

https://doi.org/10.1126/science.289.5476.113
- PubMed
- Google Scholar
1. Blanchoud S
2. Busso C
3. Naef F
4. Gönczy P
(2015) Quantitative analysis and modeling probe polarity establishment in C. elegans embryos
Biophysical Journal 108:799–809.

https://doi.org/10.1016/j.bpj.2014.12.022
- PubMed
- Google Scholar
1. Chao JT
2. Wong AKO
3. Tavassoli S
4. Young BP
5. Chruscicki A
6. Fang NN
7. Howe LJ
8. Mayor T
9. Foster LJ
10. Loewen CJR
(2014) Polarization of the endoplasmic reticulum by ER-septin tethering
Cell 158:620–632.

https://doi.org/10.1016/j.cell.2014.06.033
- PubMed
- Google Scholar
1. Chau AH
2. Walter JM
3. Gerardin J
4. Tang C
5. Lim WA
(2012) Designing synthetic regulatory networks capable of self-organizing cell polarization
Cell 151:320–332.

https://doi.org/10.1016/j.cell.2012.08.040
- PubMed
- Google Scholar
Software
1. Chen Y
(2025) PolarSim, version swh:1:rev:220bd11e10dc3f97d10e863f4bf0be0505a12011
Software Heritage.

https://archive.softwareheritage.org/swh:1:dir:764f088785dab5f8a6f4a962997d06894dd89149;origin=https://github.com/YixuanChen0726/Cell-Polarization;visit=swh:1:snp:eaddeda940f355f8d2202c9e0638f4985bcd93e5;anchor=swh:1:rev:220bd11e10dc3f97d10e863f4bf0be0505a12011
1. Cuenca AA
2. Schetter A
3. Aceto D
4. Kemphues K
5. Seydoux G
(2003) Polarization of the C. elegans zygote proceeds via distinct establishment and maintenance phases
Development 130:1255–1265.

https://doi.org/10.1242/dev.00284
- PubMed
- Google Scholar
1. Doe CQ
2. Bowerman B
(2001) Asymmetric cell division: fly neuroblast meets worm zygote
Current Opinion in Cell Biology 13:68–75.

https://doi.org/10.1016/s0955-0674(00)00176-9
- PubMed
- Google Scholar
1. Elowitz MB
2. Leibler S
(2000) A synthetic oscillatory network of transcriptional regulators
Nature 403:335–338.

https://doi.org/10.1038/35002125
- PubMed
- Google Scholar
1. Etienne-Manneville S
2. Hall A
(2002) Rho GTPases in cell biology
Nature 420:629–635.

https://doi.org/10.1038/nature01148
- PubMed
- Google Scholar
(2016) Setting the clock for fail-safe early Embryogenesis
Physical Review Letters 117:188101.

https://doi.org/10.1103/PhysRevLett.117.188101
- PubMed
- Google Scholar
1. Fickentscher R
2. Weiss M
(2017) Physical determinants of asymmetric cell divisions in the early development of Caenorhabditis elegans
Scientific Reports 7:9369.

https://doi.org/10.1038/s41598-017-09690-4
- PubMed
- Google Scholar
(2018) Anti-correlation of cell volumes and cell-cycle times during the embryogenesis of a simple model organism
New Journal of Physics 20:113001.

https://doi.org/10.1088/1367-2630/aaea91
- Google Scholar
Book
1. Fishman GS
(1995) Monte Carlo: Concepts, Algorithms, and Applications
New York, NY: Springer.

https://doi.org/10.1007/978-1-4757-2553-7
- Google Scholar
(2000) Construction of a genetic toggle switch in Escherichia coli
Nature 403:339–342.

https://doi.org/10.1038/35002131
- PubMed
- Google Scholar
(2020) Geometric cues stabilise long-axis polarisation of PAR protein patterns in C. elegans
Nature Communications 11:539.

https://doi.org/10.1038/s41467-020-14317-w
- PubMed
- Google Scholar
1. Goehring NW
2. Hoege C
3. Grill SW
4. Hyman AA
(2011a) PAR proteins diffuse freely across the anterior-posterior boundary in polarized C. elegans embryos
The Journal of Cell Biology 193:583–594.

https://doi.org/10.1083/jcb.201011094
- PubMed
- Google Scholar
1. Goehring NW
2. Trong PK
3. Bois JS
4. Chowdhury D
5. Nicola EM
6. Hyman AA
7. Grill SW
(2011b) Polarization of PAR proteins by advective triggering of a pattern-forming system
Science 334:1137–1141.

https://doi.org/10.1126/science.1208619
- PubMed
- Google Scholar
1. Goldstein B
2. Hird SN
(1996) Specification of the anteroposterior axis in Caenorhabditis elegans
Development 122:1467–1474.

https://doi.org/10.1242/dev.122.5.1467
- PubMed
- Google Scholar
1. Goodrich LV
2. Strutt D
(2011) Principles of planar polarity in animal development
Development 138:1877–1892.

https://doi.org/10.1242/dev.054080
- PubMed
- Google Scholar
(2001) CDC-42 controls early cell polarity and spindle orientation in C. elegans
Current Biology 11:482–488.

https://doi.org/10.1016/s0960-9822(01)00142-7
- PubMed
- Google Scholar
1. Gross P
2. Kumar KV
3. Goehring NW
4. Bois JS
5. Hoege C
6. Jülicher F
7. Grill SW
(2019) Guiding self-organized pattern formation in cell polarity establishment
Nature Physics 15:293–300.

https://doi.org/10.1038/s41567-018-0358-7
- PubMed
- Google Scholar
1. Guan G
2. Wong M-K
3. Zhao Z
4. Tang L-H
5. Tang C
(2021) Volume segregation programming in a nematode’s early embryogenesis
Physical Review. E 104:054409.

https://doi.org/10.1103/PhysRevE.104.054409
- PubMed
- Google Scholar
1. Guan C
2. Qiao L
3. Xiong Y
4. Zhang L
5. Jiao Y
(2022) Coactivation of antagonistic genes stabilizes polarity patterning during shoot organogenesis
Science Advances 8:eabn0368.

https://doi.org/10.1126/sciadv.abn0368
- PubMed
- Google Scholar
(2010) LGL can partition the cortex of one-cell Caenorhabditis elegans embryos into two domains
Current Biology 20:1296–1303.

https://doi.org/10.1016/j.cub.2010.05.061
- PubMed
- Google Scholar
(2019) A cell size threshold limits cell polarity and asymmetric division potential
Nature Physics 15:1075–1085.

https://doi.org/10.1038/s41567-019-0601-x
- PubMed
- Google Scholar
(2021) Physically asymmetric division of the C. elegans zygote ensures invariably successful embryogenesis
eLife 10:e61714.

https://doi.org/10.7554/eLife.61714
- PubMed
- Google Scholar
(1988) Identification of genes required for cytoplasmic localization in early C. elegans embryos
Cell 52:311–320.

https://doi.org/10.1016/s0092-8674(88)80024-2
- PubMed
- Google Scholar
1. Kim M
2. Datta A
3. Brakeman P
4. Yu W
5. Mostov KE
(2007) Polarity proteins PAR6 and aPKC regulate cell death through GSK-3beta in 3D epithelial morphogenesis
Journal of Cell Science 120:2309–2317.

https://doi.org/10.1242/jcs.007443
- PubMed
- Google Scholar
1. Knoblich JA
(2001) Asymmetric cell division during animal development
Nature Reviews Molecular Cell Biology 2:11–20.

https://doi.org/10.1038/35048085
- PubMed
- Google Scholar
(2019) Frontline Science: Escherichia coli use LPS as decoy to impair neutrophil chemotaxis and defeat antimicrobial host defense
Journal of Leukocyte Biology 106:1211–1219.

https://doi.org/10.1002/JLB.4HI0319-109R
- Google Scholar
(2016) A combined binary interaction and phenotypic map of C. elegans cell polarity proteins
Nature Cell Biology 18:337–346.

https://doi.org/10.1038/ncb3300
- PubMed
- Google Scholar
1. Kravtsova N
2. Dawes AT
(2014) Actomyosin regulation and symmetry breaking in a model of polarization in the early Caenorhabditis elegans embryo: symmetry breaking in cell polarization
Bulletin of Mathematical Biology 76:2426–2448.

https://doi.org/10.1007/s11538-014-0016-x
- PubMed
- Google Scholar
1. Kumfer KT
2. Cook SJ
3. Squirrell JM
4. Eliceiri KW
5. Peel N
6. O’Connell KF
7. White JG
(2010) CGEF-1 and CHIN-1 regulate CDC-42 activity during asymmetric division in the Caenorhabditis elegans embryo
Molecular Biology of the Cell 21:266–277.

https://doi.org/10.1091/mbc.e09-01-0060
- PubMed
- Google Scholar
1. Lang CF
2. Munro E
(2017) The PAR proteins: from molecular circuits to dynamic self-stabilizing cell polarity
Development 144:3405–3416.

https://doi.org/10.1242/dev.139063
- PubMed
- Google Scholar
1. Lim YW
2. Wen F-L
3. Shankar P
4. Shibata T
5. Motegi F
(2021) A balance between antagonizing PAR proteins specifies the pattern of asymmetric and symmetric divisions in C. elegans embryogenesis
Cell Reports 36:109326.

https://doi.org/10.1016/j.celrep.2021.109326
- PubMed
- Google Scholar
1. Lin D-W
2. Liu Y
3. Lee Y-Q
4. Yang P-J
5. Ho C-T
6. Hong J-C
7. Hsiao J-C
8. Liao D-C
9. Liang A-J
10. Hung T-C
11. Chen Y-C
12. Tu H-L
13. Hsu C-P
14. Huang H-C
(2021) Construction of intracellular asymmetry and asymmetric division in Escherichia coli
Nature Communications 12:888.

https://doi.org/10.1038/s41467-021-21135-1
- PubMed
- Google Scholar
1. Ma W
2. Trusina A
3. El-Samad H
4. Lim WA
5. Tang C
(2009) Defining network topologies that can achieve biochemical adaptation
Cell 138:760–773.

https://doi.org/10.1016/j.cell.2009.06.013
- PubMed
- Google Scholar
1. Ma R
2. Laan L
3. Dogterom M
4. Pavin N
5. Jülicher F
(2014) General theory for the mechanics of confined microtubule asters
New Journal of Physics 16:013018.

https://doi.org/10.1088/1367-2630/16/1/013018
- Google Scholar
1. Macara IG
(2004) Parsing the polarity code
Nature Reviews Molecular Cell Biology 5:220–231.

https://doi.org/10.1038/nrm1332
- PubMed
- Google Scholar
1. Morita Y
2. Seirin-Lee S
(2021) Long time behavior and stable patterns in high-dimensional polarity models of asymmetric cell division
Journal of Mathematical Biology 82:66.

https://doi.org/10.1007/s00285-021-01619-w
- PubMed
- Google Scholar
1. Motegi F
2. Zonies S
3. Hao Y
4. Cuenca AA
5. Griffin E
6. Seydoux G
(2011) Microtubules induce self-organization of polarized PAR domains in Caenorhabditis elegans zygotes
Nature Cell Biology 13:1361–1367.

https://doi.org/10.1038/ncb2354
- PubMed
- Google Scholar
(2004) Cortical flows powered by asymmetrical contraction transport PAR proteins to establish and maintain anterior-posterior polarity in the early C. elegans embryo
Developmental Cell 7:413–424.

https://doi.org/10.1016/j.devcel.2004.08.001
- PubMed
- Google Scholar
1. Nance J
(2014) Getting to know your neighbor: cell polarization in early embryos
The Journal of Cell Biology 206:823–832.

https://doi.org/10.1083/jcb.201407064
- PubMed
- Google Scholar
1. Pavin N
2. Laan L
3. Ma R
4. Dogterom M
5. Jülicher F
(2012) Positioning of microtubule organizing centers by cortical pushing and pulling forces
New Journal of Physics 14:105025.

https://doi.org/10.1088/1367-2630/14/10/105025
- Google Scholar
1. Rodriguez J
2. Peglion F
3. Martin J
4. Hubatsch L
5. Reich J
6. Hirani N
7. Gubieda AG
8. Roffey J
9. Fernandes AR
10. St Johnston D
11. Ahringer J
12. Goehring NW
(2017) aPKC cycles between functionally distinct PAR protein assemblies to drive cell polarity
Developmental Cell 42:400–415.

https://doi.org/10.1016/j.devcel.2017.07.007
- PubMed
- Google Scholar
1. Rose L
2. Gönczy P
(2014) Polarity establishment, asymmetric division and segregation of fate determinants in early C. elegans embryos
WormBook 1:1–43.

https://doi.org/10.1895/wormbook.1.30.2
- PubMed
- Google Scholar
1. Sailer A
2. Anneken A
3. Li Y
4. Lee S
5. Munro E
(2015) Dynamic opposition of clustered proteins stabilizes cortical polarity in the C. elegans zygote
Developmental Cell 35:131–142.

https://doi.org/10.1016/j.devcel.2015.09.006
- PubMed
- Google Scholar
1. Schneider SQ
2. Bowerman B
(2003) Cell polarity and the cytoskeleton in the Caenorhabditis elegans zygote
Annual Review of Genetics 37:221–249.

https://doi.org/10.1146/annurev.genet.37.110801.142443
- PubMed
- Google Scholar
(2000) MEX-5 and MEX-6 function to establish soma/germline asymmetry in early C. elegans embryos
Molecular Cell 5:671–682.

https://doi.org/10.1016/s1097-2765(00)80246-4
- PubMed
- Google Scholar
1. Seirin-Lee S
2. Shibata T
(2015) Self-organization and advective transport in the cell polarity formation for asymmetric cell division
Journal of Theoretical Biology 382:1–14.

https://doi.org/10.1016/j.jtbi.2015.06.032
- PubMed
- Google Scholar
1. Seirin-Lee S
(2016) Positioning of polarity formation by extracellular signaling during asymmetric cell division
Journal of Theoretical Biology 400:52–64.

https://doi.org/10.1016/j.jtbi.2016.04.004
- PubMed
- Google Scholar
1. Seirin-Lee S
(2020a) Asymmetric cell division from a cell to cells: shape, length, and location of polarity domain
Development, Growth & Differentiation 62:188–195.

https://doi.org/10.1111/dgd.12652
- PubMed
- Google Scholar
(2020b) CDC-42 Interactions with Par proteins are critical for proper patterning in polarization
Cells 9:2036.

https://doi.org/10.3390/cells9092036
- PubMed
- Google Scholar
1. Seirin-Lee S
2. Sukekawa T
3. Nakahara T
4. Ishii H
5. Ei S-I
(2020c) Transitions to slow or fast diffusions provide a general property for in-phase or anti-phase polarity in a cell
Journal of Mathematical Biology 80:1885–1917.

https://doi.org/10.1007/s00285-020-01484-z
- PubMed
- Google Scholar
1. Seirin-Lee S
(2021) The role of cytoplasmic MEX-5/6 polarity in asymmetric cell division
Bulletin of Mathematical Biology 83:29.

https://doi.org/10.1007/s11538-021-00860-0
- PubMed
- Google Scholar
(1983) The embryonic cell lineage of the nematode Caenorhabditis elegans
Developmental Biology 100:64–119.

https://doi.org/10.1016/0012-1606(83)90201-4
- PubMed
- Google Scholar
1. Sun Z
2. Wei W
3. Zhang M
4. Shi W
5. Zong Y
6. Chen Y
7. Yang X
8. Yu B
9. Tang C
10. Lou C
(2022) Synthetic robust perfect adaptation achieved by negative feedback coupling with linear weak positive feedback
Nucleic Acids Research 50:2377–2386.

https://doi.org/10.1093/nar/gkac066
- PubMed
- Google Scholar
Software
1. The MathWork Inc
(2022) The mathworks inc, version 9.13.0 (R2022b)
MATLAB.

https://in.mathworks.com/products/matlab.html
1. Tian B
2. Guan G
3. Tang L-H
4. Tang C
(2020) Why and how the nematode’s early embryogenesis can be precise and robust: a mechanical perspective
Physical Biology 17:026001.

https://doi.org/10.1088/1478-3975/ab6356
- PubMed
- Google Scholar
1. Tostevin F
2. Howard M
(2008) Modeling the establishment of PAR protein polarity in the one-cell C. elegans embryo
Biophysical Journal 95:4512–4522.

https://doi.org/10.1529/biophysj.108.136416
- PubMed
- Google Scholar
1. Wang JT
2. Seydoux G
(2013) Germ cell specification
Advances in Experimental Medicine and Biology 757:17–39.

https://doi.org/10.1007/978-1-4614-4015-4_2
- PubMed
- Google Scholar
1. Wang S-C
2. Low TYF
3. Nishimura Y
4. Gole L
5. Yu W
6. Motegi F
(2017) Cortical forces and CDC-42 control clustering of PAR proteins for Caenorhabditis elegans embryonic polarization
Nature Cell Biology 19:988–995.

https://doi.org/10.1038/ncb3577
- PubMed
- Google Scholar
1. Wang T
2. You L
(2020) The persistence potential of transferable plasmids
Nature Communications 11:5589.

https://doi.org/10.1038/s41467-020-19368-7
- Google Scholar
(2023) Synthetic Par polarity induces cytoskeleton asymmetry in unpolarized mammalian cells
Cell 186:4710–4727.

https://doi.org/10.1016/j.cell.2023.08.034
- PubMed
- Google Scholar
1. Zhou Z
2. Liu Y
3. Feng Y
4. Klepin S
5. Tsimring LS
6. Pillus L
7. Hasty J
8. Hao N
(2023) Engineering longevity-design of a synthetic gene oscillator to slow cellular aging
Science 380:376–381.

https://doi.org/10.1126/science.add7631
- PubMed
- Google Scholar
(2023) Reconstitution of morphogen shuttling circuits
Science Advances 9:eadf9336.

https://doi.org/10.1126/sciadv.adf9336
- PubMed
- Google Scholar

Article and author information

Author details

Yixuan Chen
1. South Bay Interdisciplinary Science Center, Songshan Lake Materials Laboratory, Dongguan, China
2. Department of Physics, Sichuan University, Chengdu, China
3. School of Physics, Peking University, Beijing, China
Contribution
Conceptualization, Data curation, Software, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing – original draft

Contributed equally with
Guoye Guan

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0009-0004-3695-8752
Guoye Guan
1. South Bay Interdisciplinary Science Center, Songshan Lake Materials Laboratory, Dongguan, China
2. Center for Quantitative Biology, Peking University, Beijing, China
3. Department of Physics, Hong Kong Baptist University, Hong Kong, China
Present address
1. Department of Systems Biology, Harvard Medical School, Harvard Medical School, Boston, United States
2. Department of Data Science, Dana-Farber Cancer Institute, Dana-Farber Cancer Institute, Boston, United States
Contribution
Conceptualization, Data curation, Software, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing – original draft

Contributed equally with
Yixuan Chen

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0003-4479-4722
Lei-Han Tang
1. South Bay Interdisciplinary Science Center, Songshan Lake Materials Laboratory, Dongguan, China
2. Department of Physics, Hong Kong Baptist University, Hong Kong, China
Present address
Center for Interdisciplinary Studies, Westlake University, Hangzhou, China

Contribution
Resources, Supervision, Funding acquisition, Validation, Investigation, Methodology, Project administration, Writing – review and editing

For correspondence
tangleihan@westlake.edu.cn

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-2925-4691
Chao Tang
1. School of Physics, Peking University, Beijing, China
2. Center for Quantitative Biology, Peking University, Beijing, China
3. Peking-Tsinghua Center for Life Sciences, Peking University, Beijing, China
Present address
Center for Interdisciplinary Studies, Westlake University, Hangzhou, China

Contribution
Resources, Supervision, Funding acquisition, Validation, Investigation, Methodology, Project administration, Writing – review and editing

For correspondence
tangc@pku.edu.cn

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0003-1474-3705

Funding

National Natural Science Foundation of China (12090053)

Chao Tang

National Natural Science Foundation of China (32088101)

Chao Tang

Research Grants Council of the Hong Kong SAR (12303219)

Lei-Han Tang

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We thank Prof. Wei Wang, Dr. Daimin Li, Peng Xie, Shiyu Shen, and Qi Ding for their assistance in the numerical simulation, and Prof. Hongli Wang for his constructive advice on the project. We appreciate Zhengyang Han for his assistance in improving the paper materials. We are grateful to Dr. Siyu Li, Dr. Baoshuang Shang, Dr. Shipu Xu, and Yuqing Zhong for their help during Yixuan Chen and Guoye Guan’s academic visit to Songshan Lake Materials Laboratory. Gratitude is extended to Prof. Zhongying Zhao, Dr. Dongying Xie, and Yiming Ma for their help during Guoye Guan’s academic visit to Hong Kong Baptist University. This work was supported by funding from the National Natural Science Foundation of China (12090053 and 32088101) and the Research Grants Council of the Hong Kong SAR (12303219). Computation was performed partly on the High-Performance Computing Platform at Peking University.

Version history

Preprint posted: January 14, 2024
Sent for peer review: March 1, 2024
Reviewed Preprint version 1: June 21, 2024
Reviewed Preprint version 2: April 16, 2025
Version of Record published: July 22, 2025

Cite all versions

You can cite all versions using the DOI https://doi.org/10.7554/eLife.96421. This DOI represents all versions, and will always resolve to the latest one.

Copyright

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

932

views
61

downloads
0

citations

Views, downloads and citations are aggregated across all versions of this paper published by eLife.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Mendeley

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

Yixuan Chen
Guoye Guan
Lei-Han Tang
Chao Tang

(2025)

Balancing reaction-diffusion network for cell polarization pattern with stability and asymmetry

eLife 13:RP96421.

https://doi.org/10.7554/eLife.96421.3

Categories and tags

Research organism

C. elegans

Share this article

Cite this article

The single modification on the antagonistic 2-node network causes the collapse of cell polarization pattern.

The combination of two opposite modifications recovers the cell polarization pattern stability.

The balance between system parameters is needed for maintaining pattern stability.

Adopting parameter sets corresponding to opposite interface velocities on two sides of the interface, the zero-velocity interface position of the polarity pattern is shiftable and regulatable.

The molecular interaction network in C. elegans zygote and its natural advantages in terms of pattern stability, balanced network configuration, viable parameter sets, and parameter robustness.

The control of the interface velocity and position by adjusting parameters in a multi-dimensional system.

The parameter description of the reaction-diffusion model for simulating a cell polarization network.

Viable parameter set searched by Monte Carlo method for the simple 2-node network.

Viable parameter set searched by Monte Carlo method for the C. elegans 5-node network.

The instructions to open the PolarSim GUI.

The instructions for the format of Network Parameter in an Excel table.

The examples of network parameters.

The progress bar, files, and subfolders of PolarSim.

The results of Example 1 (‘Simple Antagonistic 2-Node Network.xlsx’) and Example 2 (‘C. elegans Wild-Type 5-Node Network.xlsx’) of PolarSim.

The results of Example 3 (‘LGL-1 Mutant 4-Node Network.xlsx’) of PolarSim.

The effects of cell size (length) on the cell polarization pattern.

Author details

Yixuan Chen

Contribution

Contributed equally with

Competing interests

Guoye Guan

Present address

Contribution

Contributed equally with

Competing interests

Lei-Han Tang

Present address

Contribution

For correspondence

Competing interests

Chao Tang

Present address

Contribution

For correspondence

Competing interests

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

Categories and tags

Research organism