1. 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 et al., Nat. Rev. Mol. Cell Biol., 2001; Goodrich et al., Development, 2011]. Polarity establishment is an essential step in a wide range of biological phenomena, including embryonic development, wound healing, immune activity, and so forth [Etienne-Manneville et al., Nature, 2002; Bilder et al. Science, 2000]. During cytokinesis, a polarized cell can allocate its molecular contents unequally to its daughter cells, leading to asymmetric cell size and cell fate [Macara, Nat. Rev. Mol. Cell. Biol., 2004]. At the multicellular level, a group of polarized cells can undergo collective movements and constitute stereotypic architectures [Etienne-Manneville et al., Nature, 2002; Bilder et al. Science, 2000]. Thus, loss or disorder of cell polarization could severely violate normal biological processes, for example, resulting in embryonic lethality and cancerous tumor [Bilder et al., Science, 2000; Kim et al., J. Cell. Sci., 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 paid to uncover both the regulatory pathways and design principles involved [Doe et al., Curr. Opin. Cell Biol., 2001; Koorman, Nat. Cell Biol., 2016; Tostevin et al., Biophys J, 2008; Chau et al., Cell, 2012].

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 et al., WormBook, 2014; Lang et al., Development, 2017]. After fertilization, the C. elegans zygote is polarized where the entry side of the sperm (one pole of the ellipsoidal egg) turns into the posterior of the embryo [Goldstein et al., Development, 1996; Motegi et al., Nat. Cell Biol., 2011]. Driven by cell polarization, the embryo proceeds tough four successive rounds of asymmetric cell divisions and produces four somatic founder cells sequentially [Sulston et al., Dev. Biol., 1983; Hubatsch et al., Nat. Phys., 2019].

The asymmetric division in the C. elegans zygote is governed by a protein family termed partitioning-defective protein (PAR) [Kemphues et al., Cell, 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., Development, 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 through 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 et al., WormBook, 2014; Lang et al., Development, 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., J. Cell Biol., 2011; Rose et al., WormBook, 2014; Beatty et al., Development, 2013]. The cell polarization patterns, defined by the 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., J. Cell Biol., 2011; Rose et al., WormBook, 2014; Jankele et al., eLife, 2021].

A series of experiments (including 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., Development, 2003; Motegi et al., Nat. Cell Biol., 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 et al., Biophys J, 2008; Chau et al., Cell, 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., Mol. Cell, 2000; Wang et al., Adv. Exp. Med. Biol., 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 the CDC-42, LGL-1, and, CHIN-1 [Kumfer et al., Mol. Biol. Cell, 2010; Sailer et al., Dev. Cell, 2015; Beatty et al., Development, 2010]. The addition of these players increases the complexity of the cell polarity 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 et al., Mol. Hum. Reprod., 2016; Lim et al., Cell Rep., 2021]. Therefore, a more general understanding of the stability and asymmetry of cell polarization patterns will help to streamline experimental data analysis and interpretation. In a separate development, 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 since over a decade ago [Lin et al., Nat. Commun., 2021; Watson et al., Cell, 2023; Tostevin et al., Biophys J, 2008; Chau et al., Cell, 2012].

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 when translational symmetry is assumed. Unbalanced modification of kinetic parameters triggers movement of the interface in favor of one of the coexisting domains. There are three types of such unbalanced modification: 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.

2. Results

2.1 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., Science, 2011; Lang et al., Development, 2017], we propose a computational pipeline for numerical exploration of the dynamics of a given reaction-diffusion network's dynamics, specifically targeting the maintenance phase of stable cell polarization after its initial establishment [Motegi et al., Nat. Cell Biol., 2011; Goehring et al., Science, 2011; Seirin-Lee et al., Cells, 2020]. 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 the cell cytosol, inhibits the association of the other group in the same region, thereby creating two types of distinct domains [Knoblich et al., Nat. Rev. Mol. Cell Biol., 2001; Doe et al., Curr. Opin. Cell Biol., 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:

1. 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 (Fig. 1a) [Gross et al., Nat. Phys., 2019; Seirin-Lee, Dev. Growth Differ., 2020].

2. A molecular species of type [X] can associate with the cell membrane from cell cytosol at a rate and dissociate from the cell membrane into cell cytosol at a rate , where t represents time. Both rates include leaky term and regulatory pathways affected by other molecular species [Seirin-Lee et al., J. Theor. Biol., 2015; Seirin-Lee et al., Cells, 2020].

3. 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., J. Cell Biol., 2011; Lim et al., Cell Rep., 2021]. Consequently, we adopted 2.748×10^-5 and 1.472×10^-5 from previous research (nondimensionalized based on experimental measurements [Goehring et al., J. Cell Biol., 2011; Goehring et al., Science, 2011; Seirin-Lee et al., Cells, 2020], as diffusion coefficients of anterior and posterior molecular species, while the cytoplasmic molecules are regarded as well-mixed with infinite diffusion [Kravtsova et al., Bull. Math. Biol., 2014; Gross et al., Nat. Phys., 2019; Morita et al., J. Math. Biol., 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 et al., Bull. Math. Biol., 2014; Goehring et al., Science, 2011].

For the concentration of each molecular species [X] on the cell membrane [X_m], the evolution equation is described by:

where 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 the cytosol; while the third term represents the dissociation from the cell membrane into the cytosol. These regulated association and dissociation are given by the Hill equation [Seirin-Lee et al., J. Theor. Biol., 2015; Seirin-Lee, Bull. Math. Biol., 2021]:

where consists of a basal on-rate (i.e., association rate) ϒ_X and recruitment by other membrane-bound molecular species (Y ≠ X) and itself (Y = X, i.e., self-activation), quantified by positive parameters consists of a basal off-rate (i.e., dissociation rate) α_X and exclusion by other molecular species (Y ≠ X) and itself (Y = X, i.e., self-inhibition), quantified by positive parameters ; the Hill coefficients and are set to 2 as used before [Seirin-Lee et al., J. Theor. Biol., 2015; Seirin-Lee, Bull. Math. Biol., 2021]. To simplify the model for numerical exploration with an affordable computational cost, we set the same responsive concentration for both activation and inhibition pathways; the inhibition intensities are set to 1 as a numerically normalized value. We also set self and mutual activation parameters 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₁, and q₂. A detailed description of the parameters is listed in Table S1. This dimensionally reduced parameter configuration is sufficient to describe the temporal evolution of polarization distribution 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₁, q₂) on a three-dimensional (3D) grid ϒ∈ [0,0.05] in steps Δϒ = 0.001, k₁ ∈ [0,5] in steps Δk₁ = 0.05, and q₂ ∈ [0,0.05] in steps Δq₂ = 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 (Supplemental Text). 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 ϒ ∈ [0,0.05], α ∈ [0,0.05], k₁ ∈ [0,5], q₂ ∈ [0,5], q₁ ∈ [0,5], q₂ ∈ [0,0.05], and [X_c] ∈ [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, Springer New York, 1995].

Focusing on the stability and asymmetry of cell polarization 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., Dev. Cell, 2004; Kravtsova et al., Bull. Math. Biol., 2014], using the sigmoid function:

Then the set of partial differential equations (1) evolves for 500 steps with a time step of 1 (Fig. S1). The solution at t = 500 is saved for further analysis if:

1. All the molecular species still have a polarized pattern, as defined by [X_m](x = −L/2, t = 500) > 0.5 and [X_m](x = L/2, t = 500) < 0.5 for the anterior molecular species and [X_m](x = −L/2, t = 500) < 0.5 and [X_m](x = L/2, t = 500) > 0.5 for the posterior molecular species.

2. The cell polarization pattern is stable or nearly intact over time for each molecular species, as defined by

   

Implemented on Matlab 2022b [The MathWorks Inc., 2022b], such a computational pipeline can describe the spatiotemporal dynamics of any reaction-diffusion network design, for evaluating its capability of maintaining cell polarization (Fig. S1). Meanwhile, we construct an automatic software, PolarSim, for users to explore networks with arbitrarily set node numbers, parameter values, and regulatory pathways (Supplemental Text).

2.2 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., Cell, 1988; Cuenca et al., Development, 2003; Tostevin et al., Biophys. J., 2008; Chau et al., Cell, 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₂ is not considered, the computational pipeline above establishes a total of 122 viable parameter sets (ϒ, k₁) (~2.37% among 5,151 sets) that can achieve stable cell polarization. We use (ϒ = 0.05, k₁ = 0.05) as a representative to show how the corresponding cell polarization pattern behaves under alternative initial condition and elementary modification (Fig. 1, Fig. S2).

Regarding the initial conditions, while uniform initial distributions maintain a homogeneous pattern (Fig. S2a, b), a polarized initial distribution reliably maintains the cell polarization regardless of variations in the initial distribution, such as shifts in transition planes (Fig. S2c, d), weak polarization (Fig. S2e, f), or asymmetric curve shapes between the two molecular species (Fig. S2f, g). In contrast, an initial distribution with only Gaussian noise produces a multipolar pattern, with the two molecular species still locally excluding each other (Fig. S2h). 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., Development, 2003; Gross et al., Nat. Phys., 2019].

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

1. Single-sided self-regulation (Fig. 1b): a self-activation () or self-inhibition () is added on the node [A].

2. Single-sided additional regulation (Fig. 1c and Fig. S3): a new node [L] is added in the anterior or posterior, with mutual activation () or inhibition () 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.

3. Unequal system parameters (Fig. 1d): on one hand, the inhibition intensity k₂ 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 protein distribution dynamics of eight modified conditions are simulated as shown in Fig. 1 and Movie S1, 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₁, k₂, and [X_c]) assigned with various values concerning all nodes and regulations (Fig. S4, Fig. S5 - 1^st row, and Table S1). 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.

Fig. 1.

2.3 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 [Hubatsch et al., Nat. Phys., 2019; Wang et al., Adv. Exp. Med. Biol., 2013]. The simplest idea for stability recovery is to combine two kinds of modifications with opposite trends, for example, adding self-activation and posteriorly inhibition on [A] simultaneously. For the three types of modification, we arbitrarily select one subtype within each of them from Fig. 1b-d so three combinatorial networks are constituted, all of whose interfaces turn out to be stabilized finally (Fig. 2 and Movie S2).

Fig. 2.

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 Fig. 2), where the two corresponding intensities and compose a phase diagram that distinguishes the final state of the reaction-diffusion pattern. The moving velocity of the pattern is defined as follows:

where N is the total number of the molecular species. Here, we classify the final state by calculating the moving velocity of the pattern at the final time t = 500; the region with v < 10⁻⁴ 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 (Fig. 3a). Therefore, the triphase diagram and the exemplary patterns generated by the parameter assignments ①②③ shown in Fig. 3a suggest the necessity of the balance between two opposite modifications for setting up a stable cell polarization pattern.

Apart from the unequal inhibition intensity and cytoplasmic concentration mentioned in Fig. 2d, we further ask if all the kinetic parameters with biophysical significance (including inhibition intensity k₂, responsive concentration k₁, 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 . Fascinatingly, monotonic correlations exist between all those system parameters, suggesting that they can be tuned to maintain pattern stability (Fig. 3b). In the grey region, the symmetric-broken parameters can be weighed against each other to realize the interface velocity close to zero.

Fig. 3.

2.4 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, there exists influence from cellular actomyosin flow and extracellular/intercellular signals, leading to changeable molecular concentration and parameter values in different parts of space [Arata et al., Development, 2010; Beatty et al., Development, 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 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₁ = 0.05, Fig. 4a, the same simulation as the one in Fig. 1a), we alter its initial condition by shifting the interface position or exerting a homogeneous or Gaussian noise distribution without any positional information (Fig. S2); it is found that the initial polarized pattern is essential for constricting only one single interface with predetermined position, while the homogeneous or Gaussian noise distribution results in no or multiple interfaces (Fig. S2a, 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., Development, 2003; Motegi et al., Nat. Cell Biol., 2011].

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

1. Increasing the inhibitory intensity of [A] on [P], , from 1 to 1.5, leads to the interface continuously shifting to the posterior, and [A] finally dominates.

2. Reducing the basal on-rate of [A], ϒ_A, from 0.05 to 0.01, results in the interface moving to the anteriority, 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 (Fig. 4b). By moving the step position of the parameter function to x = −0.1 (Fig. 4c) and x = 0.1 (Fig. 4d), the interface stabilizes around the step, which means the interface position is tunable (Movie S3). This indicates that by using parameter sets with values corresponding to opposite interface velocities on two sides of the interface, the interface velocity and position of the polarity pattern can be precisely regulated. The nonuniform 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.

Fig. 4.

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

In the 2-node network, the cell polarization pattern stability could be broken by a single modification (Fig. 1) and recovered by two combinatorial modifications (Fig. 2, Fig. 3), or with the zero-velocity interface position regulated by spatially inhomogeneous parameters (Fig. 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 3 decades [Kemphues et al., Cell, 1988; Lang et al., Development, 2017; Goehring et al., Science, 2011; Lim et al., Cell Rep., 2021]. By conducting an exhaustive literature search (a total of 19 references), we summarize the mutual 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 (Fig. 5a). The detailed description of the biochemical mechanism of each regulatory pathway as well as the corresponding supporting references are listed in Table S2. Based on the network obtained experimentally, the computational pipeline described in Section 2.1 establishes a total of 602 viable parameter sets (ϒ, k₁, q₂) (~0.229% among 262,701) that can achieve stable cell polarization (Fig. 5d, e).

To verify whether the computational pipeline is reliable enough to simulate the dynamics of protein distribution 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 (Fig. 5b) with full connections, and the [A]~[C] mutual activation (Fig. 5c) and [A]~[L] mutual inhibition (Fig. 5d). As the antagonistic 2-node network has been revealed to be a fundamental structure capable of stable cell polarization [Tostevin et al., Biophys J, 2008; Goehring et al., Science, 2011; Chau et al., Cell, 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 (Fig. S6).

Taking [L] as the first example for its current absence of theoretical model analysis [Seirin-Lee et al., Cells, 2020, Lim et al., Cell Rep., 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 and proven to be a consequence of polarity loss in the zygote, which is essential for asymmetric cell division and cell differentiation) [Beatty et al., Development, 2010; Beatty et al., Development, 2013; Rodriguez et al., Dev. Cell, 2017; Jankele et al., eLife, 2021]. Here, we utilize the pattern error in a mutant (abbr., MT) embryo compared to the wild-type (abbr., WT) one, , to represent level of polarity loss in simulation and lethality in experiment (measured from previous literature), where M = 602 is the number of viable parameter sets (~0.229% among 262,701). The first group of experiments 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 Fig. S7a, 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., Curr. Biol., 2010]. The second group of experiments 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 Fig. S7c, 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., Development, 2010; Beatty et al., Development, 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 pattern polarity in embryo individuals. The first group of experiments is that the single depletion on [C] results in polarity defects and the mislocalization of PAR-6 and PAR-2, which is recapitulated by Fig. S8a, b (cell polarization pattern characteristics in simulation: less polarized distributions of [A] and [P], compared to WT) [Gotta et al., Curr. Biol., 2001; Aceto et al., Dev. Biol., 2006]. The second group of experiments is that the overexpression of [C] increases the size of the anterior domain, which is recapitulated by Fig. S8a, c (cell polarization pattern characteristics in simulation: posteriorly shifted interfaces of [A] and [P], compared to WT) [Aceto et al., Dev. Biol., 2006]. The third group of experiments 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 Fig. S8a, d (cell polarization pattern characteristics in simulation: less polarized distribution of [A] in both conditions, compared to WT) [Aceto et al., Dev. Biol., 2006].

To sum up, apart from the cell polarization in a wild-type embryo, our modeling framework is further validated by reproducing two groups of perturbation experiments that haven’t been theoretically explained before, allowing further computational investigation of the network dynamics.

Fig. 5.

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 the posterior (i.e., [A]~[L]) subsequently, turning into an asymmetric 5-node network (Fig. 5b-d). 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 Fig. 2 and Fig. 3a.

Here, we seek to compare the concentration distribution and moving velocity of the patterns generated by the three successive network structures termed “4-Node” (the symmetric structure) (Fig. 5b), “LGL-1^-” (the symmetric structure with [A]~[C] added, i.e., the mutant network with [A]~[L] depleted from the wild-type network) (Fig. 5c), and “WT” (the symmetric structure with both [A]~[C] and [A]~[L] added, i.e., the wild-type network) (Fig. 5d). As expected, both simulations on the 4-Node and WT networks successfully pass the computation pipeline described in Section 2.1 with 62 (~1.20% among 5,151) and 602 (~0.229% among 262,701) viable parameter sets respectively (Fig. 5e), all of which stabilize into a cell polarization pattern with the concentration distribution highly intact and the moving velocity of each parameter set continuously declining below 10⁻⁴ (Fig. 5b, d and Movie S4). Nonetheless, the intermediate one (generated by depleting [A]~[L] from the WT network) with [A]~[C] but without [A]~[L] fails (Fig. 5c and Movie S4), exhibiting a dispersed concentration distribution and stubbornly high moving 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₁, k₂, q₁, q₂ and [X_c]) assigned with various values concerning all nodes and regulations (Fig. S5 - 2^nd row, Fig. S9, Table S1), just like the simple 2-node network (Fig. S4, Fig. S5 – 1^st row, Table S1). 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 (Fig. 5f). To evaluate how the interface moves, we quantify the interface velocity by the rate of the transition plane 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.

Scanning the intensities of [A]~[C] mutual activation (i.e., ) and [A]~[L] mutual inhibition (i.e., ), 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 (Fig. 5f). The redder the contour color is, the faster the interface moves posteriorly, with an increasing aPARs domain; 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., Cells, 2020; Lim et al., Cell Rep., 2021], then the existence, as well as the form of the feed loops between [L] and a preexisting molecular species, is a supporting modification that consolidates stable 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 (Fig. S10). 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₁, k₂, q₁, and, q₂, and initiate 1000 independent simulations, where the noise amplification is represented by the standard deviation σ of the Gaussian noise [Guan et al., Phys. Rev. E, 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 the viable parameter sets, and all independent simulations, is defined as follows:

where M, Q, N represents the number of the viable parameter sets (i.e. 602 (~0.229% among 262,701) for WT, 62 (~1.20% among 5,151) for 4-Node and 62 (~0.024% among 262,701) for LGL-1^-), independent simulations (i.e. 1000) and molecular species (i.e. five for WT, four for 4-Node and LGL-1^-). 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 (Fig. 5g).

2.6 A protocol to identify responsive parameters for interface positioning

For a polarization pattern to be stable, fine-tuning of the kinetic parameters is required to reach the vanishing velocity of the pattern interface (numerically identified by Eq. (5)). Using the 5-node C. elegans network as an example, we outline a method to delineate iso-velocity surfaces v_I(Ω) = constant 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 Sec. 2.5, 602 (~0.229% among 262,701) groups of parameters capable of reaching stable polarity patterns are found by scanning through parameter space. When plotted in the three-dimensional parameter space shown in Fig. 5e, they span an iso-velocity surface (blue points in Fig. 5e) at vanishing velocity, in other words, their pattern interfaces are stabilized with little or no moving. To illustrate the local orientation of the surface, we selected 9,261 points (enumerated on a 3D grid ϒ∈ [0.034,0.044] in steps Δϒ = 0.0005, k₁ ∈ [1.05,2.05] in steps Δk₁ = 0.05, and q₂ ∈ [0.045,0.055] in steps Δq₂ = 0.0005) in the neighborhood of a benchmark point , as shown by the brown box in Fig. 6a, top. The least-square fit of computed interface velocity to a linear function yields:

with the coefficient of determination:

The coefficients in Eq. (9) give the ascending gradient of the interface velocity (Fig. 6a, bottom).

Within the assumed relationships adopted for parameter reduction (see Sec. 2.1), Eq. (9) shows a strong dependence of the interface velocity on basal on/off-rates ϒ and activation intensity q₂, but weak sensitivity to responsive concentration k₁ in the inspected parameter region. Interestingly, increase in the basal on/off-rates ϒ and the self-activation rate q₂ have opposite effects on the interface velocity, which can be attributed to the asymmetric 5-node network topology. In a more realistic scenario, one may consider 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 Eq. (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 polarity circuit is incomplete. In particular, the parameter-dependent interface velocity picture potentially enables biologists to manipulate and even synthesize the cell polarization 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., Nat. Phys., 2019; Guan et al., Phys. Rev. E, 2021].

Fig. 6.

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 et al., Annu. Rev. Genet., 2003; Ajduk et al., Mol. Hum. Reprod., 2016], while the acquired cell volume has been reported to have a chain effect in its cell cycle, cell position, and other behaviors [Fickentscher et al., New J. Phys., 2018; Guan et al., Phys. Rev. E, 2021; Fickentscher et al., Phys. Rev. Lett., 2016; Tian et al., Phys. Biol., 2020]. As in the 2-node network, the precise control of the stable position of the interface 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 most stable pattern (Fig. 6b and Movie S5), two modifications are adopted to generate patterns with opposite interface velocity:

1. Increasing the activation intensity between [A] and [C] () to 0.12 leads to the interface traveling backward, and [A] and [C] finally dominate the whole domain of the membrane.

2. 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 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 interface settles around x = 0 (Fig. 6c and Movie S5). The interface position turns tunable as the step switch moves to x = −0.1 (Fig. 6d and Movie S5) and x = 0.1 (Fig. 6e and Movie S5). 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), inducing polarity reversal by affecting the on-rate of PAR-2 [Arata et al., Development, 2010; Seirin-Lee, J. Theor. Biol., 2016]. Such a scheme depicts the introduction of extracellular or intracellular cues that break the spatial uniformity of parameters and tune the stable localization of interface, serving as a theoretical basis for controlling oriented cell division with designated volume partition and fate differentiation [Arata et al., Development, 2010; Hubatsch et al., Nat. Phys., 2019; Schubert et al., Mol. Cell, 2000].

3. Discussion and conclusion

Cell polarization is a fundamental issue in both prokaryotes and eukaryotes, playing crucial roles in diverse biological phenomena ranging from chemotaxis to embryogenesis [Nance, J. Cell Biol., 2014; Kondo et al., J. Leukoc. Biol., 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 [Koorman et al., Nat. Cell Biol., 2016; Lang et al., Development, 2017; Tostevin et al., Biophys J, 2008; Chau et al., Cell, 2012; Lin et al., Nat. Commun., 2021; Watson et al., Cell, 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 3 unbalanced modifications (i.e., single-sided self-regulation, single-sided additional regulation, and unequal system parameters) were introduced. To recover stable polarization, two strategies are proposed: (i) the combination of two unbalanced modifications with opposite effects; (ii) the 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. New J. Phys., 2018; Guan et al., Phys. Rev. E, 2021; Fickentscher et al., Phys. Rev. Lett., 2016; Tian et al., Phys. Biol., 2020]. Analogous conclusions are further identified in the C. elegans network summarized experimentally which obtains larger parameter space and is higher 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 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 system with desired functions. 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.

1. 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 et al., Biophys. J., 2008; Gross et al., Nat. Phys., 2019; Seirin-Lee et al., Cells, 2020; Lim et al., Cell Rep., 2021], might not be ideally identical concerning their slightly different polarized patterns [Wang et al., Nat. Cell Biol., 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., Cell Rep., 2021]. Future experimental validation is needed to detail whether and how the molecular species described in Table S2 (and beyond) interact with each other.

2. 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 et al., Nat. Commun., 2020].

3. Quantitative biological details, such as the concentration distributions, diffusions, productions, and degradations of molecular species on the cell membrane versus in the cell cytosol [Tostevin et al., Biophys. J., 2008; Goehring et al., Science, 2011; Kravtsova et al., Bull. Math. Biol., 2014; Sailer et al., Dev. Cell, 2015; Seirin-Lee, et al., J. Theor. Biol., 2015; Seirin-Lee et al., J. Theor. Biol., 2016; Gross et al., Nat. Phys., 2019; Seirin-Lee et al., J. Math. Biol., 2020; Lim et al., Cell Rep., 2021] and the precise mathematical forms of regulatory interactions (e.g., Hill equations or product equations) [Goehring et al., Science, 2011; Seirin-Lee et al., Cells, 2020; Lim et al., Cell Rep., 2021], challenge model assumptions and hinder the unification of parameter value with unit between theory and experiment [Tostevin et al., Biophys. J., 2008; Kravtsova et al., Bull. Math. Biol., 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 et al., Nature, 2000; Gardner et al., Nature, 2000; Ma et al., Cell, 2009; Chau et al., Cell, 2012].

4. 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 sperm entry and the establishment phase, which involve active actomyosin contractility and flow drive 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., Development, 2003; Motegi et al., Nat. Cell Biol., 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., J. Cell Biol., 2011; Rose et al., WormBook, 2014; Geβele et al., Nat. Commun., 2020]) and cellular (asymmetric segregation of cell volume and molecular content, mediated by molecular species like MES-1 and SRC-1 [Arata et al., Development, 2010; Seirin-Lee, J. Theor. Biol., 2016]) effect 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 [Rose et al., WormBook, 2014; Knoblich et al., Nat. Rev. Mol. Cell Biol., 2001]. Governed by the cell polarization, the C. elegans embryo actually proceeds 4 rounds of asymmetric division to produce 4 somatic founder cells, then such reaction-diffusion pattern will lose its sharp transition plane as it is unscalable over cell size, leading to symmetric cell division at last [Hubatsch et al., Nat. Phys., 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., Cell, 2014; Guan et al., Sci. Adv., 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., Science, 2023; Zhu et al., Sci. Adv., 2023]. Like toggle switch, oscillation, and adaptation [Gardner et al., Nature, 2000; Zhou et al., Science, 2023; Sun et al., Nucleic Acids Res., 2022], the polarized distribution of molecular species (accounting for cell polarization) is also another elementary behavior in a cell, and the corresponding circuits have been synthesized successfully around a decade ago [Chau et al., Cell, 2012]. However, it remains unclear whether the distribution pattern (e.g., the transition plane of a molecular species) 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 a real reaction-diffusion network with pattern stability and asymmetry, and has been packed as a user-friendly software, PolarSim (Supplemental Text); thus we believe it can be used for not only understanding the natural molecular networks reported before but also designing new molecular networks in the polarized cells with more physiological functions, for example, with tunable transition plane of molecular species to guide the allocation of downstream fate determinants and unequal cleavage during cytokinesis [Fickentscher et al., Sci. Rep., 2017]. In the future, our computational framework for cell polarization network could be linked to the one about cytoskeleton activity, achieving a more comprehensive modeling description for cell division [Pavin et al., New J. Phys., 2012; Ma et al., New J. Phys., 2014].

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 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. We appreciate Zhengyang Han for his assistance in improving the paper materials. Gratitude is extended to Prof. Zhongying Zhao for his 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.

Additional information

Author contributions

Y.C. and G.G. designed and conceived the project, and wrote the paper; Y.C., G.G., L.- H.T., and C.T. performed research and conducted analysis; L.-H.T. and C.T. supervised the project and revised the paper.

Additional files

Supplementary material.

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

Table S2. The literature summary on the cell polarization network in the C. elegans zygote.

Movie S1. The concentration distribution of each molecular species on the cell membrane over time generated by the antagonistic 2-node network and the ones with a single modification added, corresponding to Fig. 1

Movie S2. The concentration distribution of each molecular species on the cell membrane over time generated by the antagonistic 2-node network and the ones with a single modification or two combinatorial modifications added, corresponding to Fig. 2.

Movie S3. The concentration distribution of each molecular species on the cell membrane over time generated by the antagonistic 2-node network with spatially inhomogeneous parameters, corresponding to Fig. 4.

Movie S4. The concentration distribution of each molecular species on the cell membrane over time generated by the 4-Node, LGL-1-, and WT networks, corresponding to Fig. 5b-d.

Movie S5. The concentration distribution of each molecular species on the cell membrane over time generated by the C. elegans 5-node network with spatially inhomogeneous parameters, corresponding to Fig. 6b-e.