1. Introduction

Cell polarization is a biophysical and biochemical process where a cell acquires spatial anisotropy by establishing directional gradients of molecules across its membrane or in the cytosol [1, 2]. Polarity establishment is an essential step in a wide range of biological phenomena, including embryonic development, wound healing, immune activity, and so forth [3, 4]. During cytokinesis, a polarized cell can allocate its molecular contents unequally to its daughter cells, leading to asymmetric cell size and cell fate [5]. At the multicellular level, a group of polarized cells can undergo collective movements and constitute stereotypic architectures [3, 4]. Thus, loss or disorder of cell polarization could severely violate normal biological processes, for example, resulting in embryonic lethality and cancerous tumor [4, 6]. 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 [7–10].

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 [11, 12]. 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 [13, 14]. Driven by cell polarization, the embryo proceeds through four successive rounds of asymmetric cell divisions and produces four somatic founder cells sequentially [15, 16].

The asymmetric division in the C. elegans zygote is governed by a protein family termed partitioning-defective protein (PAR) [17]. 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 [18]. A series of experiments (including RNA interference, immunoprecipitations, fluorescence recovery after photobleaching, and time-lapse single-molecule imaging) further uncovered the mutual inhibition of those two groups of proteins upon their association with the membrane and demonstrated its essential role in the robust spatial separation of PAR proteins [14, 18]. 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 [9, 10]. 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 [19, 20].

In recent years, more proteins have been identified to significantly interact with the antagonistic 2-node network, such as the CDC-42, LGL-1, and, CHIN-1 in the cell polarization of C. elegans zygote [21–23]. 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 [24, 25]. Therefore a more general understanding of the stability and asymmetry of cell polarity 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 [9, 10, 26, 27].

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 simplified network and the realistic network consisting of various node numbers and regulatory pathways [12, 28], we propose a computational pipeline for numerical exploration of the dynamics of a given reaction-diffusion network capable of maintaining stable cell polarization. Numerous biological experiments on different organisms have demonstrated that such a reaction-diffusion network typically comprises two groups of molecules. 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 [1, 7]. 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 molecules accumulated on the cell membrane is denoted by [A_m]) and 0.25 (posterior, where the molecule accumulated on the cell membrane is denoted by [P_m]) are its two poles (Fig. 1) [29, 30].

2. A molecule 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 molecules [31, 32].

3. Based on previous experimental measurements, it has been reported that the diffusion rate of molecules 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 [25, 33]. Consequently, the cytoplasmic molecules can be regarded as well-mixed. We further assume that these molecules 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 [28, 34].

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 molecule [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 [25, 35]:

where consists of a basal on-rate (i.e., association rate) γ_X and recruitment by other membrane-bound molecules (Y ≠ X) and itself (Y = X, i.e., self-activation), quantified by positive parameters (i = 1,2); consists of a basal off-rate (i.e., dissociation rate) α_X and exclusion by other molecules (Y ≠ X) and itself (Y = X, i.e., self-inhibition), quantified by positive parameters (i = 1,2); the Hill coefficients and are set to 2 as used before [31, 35]. To further simplify the model for numerical exploration, we set the same responsive concentration for both activation and inhibition pathways; the inhibition intensities are set to 1, which also defines our unit of time, with a nominal value of 1 sec in the observation experiment [28, 36]. We also set self and mutual activation parameters = q₂ 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 molecules 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 initial state of the molecule distribution [X_m](x, 0) is set as a polarized pattern with the sigmoid function as follows:

Then the set of partial differential equations (1) evolves for 500 steps with a time step of 1 (or equivalent to 500 steps × 1 sec per step = 8 min 20 sec in reality), approaching the period of maintenance phase in cell polarization, ~10 min, measured in vivo [37]. The solution at t = 500 is saved for further analysis if:

1. All the molecules 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 molecule(s) and [X_m](x = –L/2, t = 500) < 0.5 and [X_m](x = L/2, t = 500) > 0.5 for the posterior molecule(s).

2. The cell polarization pattern is stable or nearly intact over time, as defined by < 10^–4 (x₀ ∈ [–0.25,0.25]) for each molecule.

Fig. 1.

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 [9, 10, 17, 18]. 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., molecules) 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₁) 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 elementary modification (Fig. 1a). In total, three types of modification are exerted on the node [A]:

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

2. Single-sided additional regulation (Fig. 1c): a new node [L] is added in the anterior or posterior, with mutual activation or inhibition ( = 0.025) with [A].

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. 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.

2.2. 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 [16, 20]. 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 molecules. 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^–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 (Fig. 3a). Therefore, the triphase diagram and the exemplary patterns generated by the parameter assignments ①②③ shown in Fig. 3a suggest the necessity of the detailed 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 have a detailed balance under the requirement of 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.3. The speed and position of the interface can be adjusted by setting up spatial cues

The detailed balance between parameters provides insight into controlling the pattern stability. However, in biological organisms, there exists interference from extracellular/intercellular signals, leading to changeable parameter values in different parts of space. How do cells maintain pattern stability in the case of non-uniform parameters, namely finding zero-velocity solutions at interfaces? Based on the symmetric 2-node network (γ = 0.05, k₁ = 0.05, Fig. 4a), we employ the following adjustment on the parameter set:

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 recovers 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 localization is tunable. This indicates that by using parameter sets with values corresponding to opposite interface velocities on two sides of the interface, the interface speed 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 localization in response to variable cues.

Fig. 4.

2.4. 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. 2) and recovered by two combinatorial modifications (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 [12, 17, 25, 28]. 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 molecules or molecular complexes: 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₂) 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 the perturbation experiments on [L] (the protein accumulated in the posterior and with a mutual inhibition with [A] in the anterior). A group of experiments were conducted by knocking down or overexpressing [L] in different conditions, followed by a measurement of the lethality (defined by death rate) in embryo individuals. Here, we utilize the pattern error in a mutant (abbr., MT) embryo compared to the wild-type (abbr., WT) one, Error = , to represent lethality qualitatively, where M = 602 is the number of parameter sets. 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 faithfully recapitulated by Fig. S2a, b [38]. 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 faithfully recapitulated by Fig. S2c, d [23, 39]. 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 and 602 solution numbers respectively (Fig. 5e), all of which stabilize into a cell polarization pattern the concentration distribution highly intact and the moving velocity of each parameter set continuously declining below 10^-4 (Fig. 5b, d). Nonetheless, the intermediate one (generated by depleting [A]~[L] from the WT network) with [A]~[C] but without [A]~[L] fails (Fig. 5c), exhibiting a dispersed concentration distribution and stubbornly high moving velocity for its pattern. Therefore, WT shows the largest variable parameter sets in the limited 3D grid and the strongest stability among the three network structures.

Remarkably, the detailed balance between these two modifications is required to achieve stable cell polarization (Fig. 5f). To evaluate how the interface moves, we quantify the position of the transition plane x_T(t) of the pattern, defined by the mean position of all the molecules where the absolute value of the curve derivative reaches its maximum. 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., ): the system can be stably polarized when they are in detailed 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 [25, 32], then the existence, as well as the form of the feed loops between [L] and a preexisting molecule, 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 34 possible network structures with an additional regulation (Fig. S3). 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 [40]. 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 molecules, 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 for WT, 62 for 4-Node and LGL-1⁻), independent simulations (i.e. = 1000) and molecules (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.5. 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 vanishing interface velocity. Using the 5-node C. elegans network as an example, we outline a method to delineate iso-velocity surfaces v_I(P) = constant in the high-dimensional space of the parameter set P. 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.

In Sec. 2.4, 602 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 at vanishing speed. To illustrate the local orientation of the surface, we selected 9261 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 P*(γ* = 0.039, = 1.55, = 0.05), 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. (8) give the ascending gradient of the interface velocity (Fig. 6a, bottom).

Within the assumed relationships adopted for parameter reduction (see Sec. 2.1), Eq. (8) 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 molecules and their pairwise interactions separately and explore the relationships in a much larger parameter space. An expression similar to Eq. (8) 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 [16, 40].

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 [24, 41], while the acquired cell volume has been reported to have a chain effect in its cell cycle, cell position, and other behaviors [40, 42–44]. 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 P* which gives out the most stable pattern (Fig. 6b), two modifications are adopted to generate patterns with opposite interface velocity: (i) increasing the activation intensity between [A] and [C] to 0.12 leads to the interface traveling backward, (ii) increasing the basal on-rate of posterior proteins (γ_P, γ_L, γ_H) to 0.06 results in the interface traveling forward. 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. 6b). The interface position turns tunable as the step switch moves to x = –0.1 (Fig. 6c) and x = 0.1 (Fig. 6d). 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 [45, 46]. 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 [16, 19, 45].

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 [47, 48]. 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 [8–10, 12, 26, 27]. 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 (molecule) 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 [40, 42–44]. 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 molecules 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, there might be more experimental details remaining to be uncovered, which are indispensable for achieving a more comprehensive understanding of this problem. For instance, in this paper, we presumed that the parameter values were the same for each molecule and each pathway at first and drew a phase diagram to visualize the trade-off between any two parameters, regardless of the others. Such simplification is unavoidable because the in vivo measurement of those massive parameters (>50 in total) is not only timeconsuming but also technically difficult, as the molecules and pathways interact with each other in a coupled manner. Besides, some molecules and pathways may have not been quantitively identified yet, especially for the later stages where the cell polarization requests more dynamic functions. Furthermore, previous works have shown that the polarized pattern is self-organized with its interface stabilizing automatically under a mass-conserving framework [9, 28, 29, 49]. However, here we provide a more robust strategy to control the location of the interface artificially by a step-like parameter shift, which can be realized by any parameter with sufficient biological dynamic range. Although nonuniform parameters imply potential signaling, additional revision on the model is needed to describe the concrete effects of external cues like MES-1 and SRC-1, which is necessary for the orientation of the cleavage spindle and polarity reversal [45, 46]. With more information on the molecules involved, regulatory pathways, and parameter values, a more comprehensive and reliable model can be built to explain more experimental observations and unveil more design principles for networks in the future.

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 [1, 11]. 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 [16]. 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 [50, 51].

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 [52, 53]. Like toggle switch, oscillation, and adaptation [52, 54, 55], the polarized distribution of molecules (accounting for cell polarization) is also another elementary behavior in a cell, and the corresponding circuits have been synthesized successfully around a decade ago [10]. However, it remains unclear whether the distribution pattern (e.g., the transition plane of a molecule) 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; 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 a molecule to guide the allocation of downstream fate determinants and unequal cleavage during cytokinesis [56]. In the future, our computational framework for cell polarity network could be linked to the one about cytoskeleton activity, achieving a more comprehensive modeling description for cell division [57, 58].

Conflict of interest

The authors declare that they have no conflict of interest.

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 their 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.

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.

Supplementary Materials

Supplemental Figures

Fig. S1.

Fig. S2.

Fig. S3.

Instructions for PolarSim

1. Introduction

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

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 (Fig. G1).

  Fig. G1.

• With the following parameters inputted, the GUI gives out a group of “Pattern_*.mat” files containing pattern information.

  1. 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 (Fig. G2). The Excel table for parameter value assignments should follow the format below:

     Table G1

     Fig. G2.

  2. 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 Fig. G6.

  3. Import: Output Folder

     Give a folder name for storing the output results (e.g. “Output 2-Node”)

  4. Setting: Time for Simulation

     Simulation duration “500” is used in this paper. Any positive number is allowed in this box.

  5. Setting: Time Step Length

     We take “1” as an example. Larger values are not recommended, in consideration of possible time failure, and one may try a smaller step length while the error tolerance is exceeded.

  6. 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 (Fig. G3a). A folder named “Import: Output Folder” is generated in the current path to store the output “Pattern_*.mat” data containing the name, location, and concentration distribution on the cell membrane of the node (molecule), where “*” denotes the in silico time corresponding to each file (Fig. G4a, right).

• Import the pathway of the file outputted by PolarSim into the box “Import: Pattern”.

• Click “Plot” and then a figure comes out 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”.

Fig. G3.

Fig. G4.

3. Examples

• Example 1: Simple Antagonistic 2-Node Network simulation as Fig. 1a.

  1. 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”. (Fig. G4a, ①)

  2. Click “Run”, and then 6 “Pattern_*.mat” files at different time points are saved in the subfolder “Output 2-Node” (Fig. G4a, ②-③).

  3. Input “Output 2-Node\Pattern_500.mat” into the box “Import: Pattern”, and a figure of Simple 2-Node Antagonistic Network at t = 500 appears. The position of the transition plane X_T = 0.00252525 is given in the box “Output: Transition Plane” (Fig. G4a, ④-⑥).

  4. 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, and then the mean interface velocity between t = 300 and t = 500 is given in the box “Output: Interface Velocity” (Fig. G4a, ⑦-⑨). The 2-node network approaches a stable polarized pattern with its interface velocity being 0.

• Example 2: C. elegans Wild-Type 5-Node Network simulation as Fig. 6b.

  1. 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” (Fig. G4b, ①).

  2. Click “Run”, and then 6 “Pattern_*.mat” files at different time points are saved in the subfolder “Output 5-Node” (Fig. G4b, ②-③).

  3. Input “Output 5-Node\Pattern_500.mat” into the box “Import: Pattern”, and a figure of C. elegans Wild-Type 5-Node Network at t = 500 appears. The position of the transition plane XT = -0.0176768 is given in the box “Output: Transition Plane” (Fig. G4b, ④-⑥).

  4. 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, and then the mean interface velocity between t = 300 and t = 500 is given in the box “Output: Interface Velocity” (Fig. G4b, ⑦-⑨). The 5-node network approaches a stable polarized pattern with its interface velocity being 0.

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

  1. 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” (Fig. G5a).

  2. Click “Run”, and then 11 “Pattern_*.mat” files at different time points are saved in the subfolder “Output 4-Node”.

  3. Input “Output 4-Node\Pattern_500.mat” into the box “Import: Pattern”, and then a figure of LGL-1 Mutant 4-Node Network at t = 500 appears with the transition plane close to the posterior pole. The position of the transition plane X_T = 0.207071 is given in the box “Output: Transition Plane” (Fig. G5a).

  4. 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, and then the mean interface velocity between t = 300 and t = 500 is given in the box “Output: Interface Velocity” (Fig. G5a). The interface of the 4-node network keeps moving toward the posterior with V_I = 0.000429293.

  5. Input “Output 4-Node\Pattern_1000.mat” in the box “Import: Pattern”, and then a figure of LGL-1 Mutant 4-Node Network at t = 1000 appears. The string “Transition plane doesn’t exist” appears in the box “Output: Transition Plane” as the pattern reaches a homogeneous state dominated by [A] and [C] (Fig. G5b).

  6. 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. The mean interface velocity between t = 800 and t = 1000 doesn’t exist with the string “It isn’t a polarized pattern” appearing in the box “Output: Interface Velocity” to hint (Fig. G5b). Note that the interface velocity can’t be calculated when either pattern at start time or end time is homogeneous.

Fig. G5.

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 two examples to see whether there is a cell size threshold limiting cell polarization as discovered before [16]. Patterns at t = 500 are plotted with the cell length ranging from 0.1 to 0.5 in steps 0.1 (Fig. G6).

Fig. G6.

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. It’s worth noting that the collapse of cell polarization pattern over cell size decrease may be attributed to two different mechanisms - while the curve of concentration distribution changes with the interface becoming indistinguishable as proposed in [16], our new results further suggest that the whole curve of concentration distribution will even turn fully homogeneous when the cell size is too small. In all, PolarSim provides a user-friendly tool for more applications on the studies in cell polarization.

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.