A mathematical framework for understanding the spontaneous emergence of complexity applicable to growing multicellular systems

We began our investigation by studying the simplest division-decision system, in which a one-dimensional array is formed by consecutive divisions starting from a single cell. The binary expression states of a single gene A, or 1, characterize the state of the i-th cell as . The rules of gene expression updating were described in Eq 2 in the Method section.

Gene expression at time t + 1 is determined by a step function of the regulation-weighted summation of gene expressions at time t, with the threshold of the step function set to 0 [38][39]:
(Eq 5)

For one dimensional system, there are nine possible genetic networks (Figs 1B and S2A). As the gene expression is binary, the maximum positional information in this single-gene system is 1.0 bit. We simulated all nine genetic networks starting from a single cell to assess the maximum positional information achievable by the stable states. In each generation, the state of each cell was updated twenty times before simultaneous division of all cells in that generation.

Fig 1. One-dimensional division-decision growing model and its maximum information pattern in the single-gene system.

(a). A diagram illustrating our one-dimensional division-decision growing model. In our model, cells with the identical genome divide and differentiate from a single cell state to produce a multicellular system with several cell states. Distinct colors are used to indicate different cell states. (b). Possible intracellular and intercellular interactions in the single-gene system. There are nine possible types of genetic networks. The black arrow represents the possible interaction relationship, like an activation, an inhibition or an absence of interaction. Solid and dotted lines represent intracellular and intercellular interaction relationships, respectively. (c). A distribution diagram of the maximum positional information in the single-gene system. (d). The exclusive genetic network which controls maximum positional information pattern formation in the single-gene system. (e). The pattern of gene expression varies over time in a single-gene system regulated by the genetic network O1. At step 8, the system obtains 1.0 bit positional information. The time-course of positional information is calculated and shown in the right of the pattern. The binarized gene A expression level determines the cell state. Green circles represent gene A expressed while yellow circles represent gene A unexpressed. (f). Explanation of increased positional information in the single-gene system. The maximum positional information 1.0 bit is created by the new space network resulted from cell division. Cell 4, cell 5, cell 6, and cell 7 inherit cell states from their mothers when cell 2 and cell 3 divide. The expression state matrix (S) and the cell-contacting state matrix (CS) comprised the concatenated matrix ([S,C∙S]), which changed after cell division. The new concatenated matrix contains 1.0 bit potential information which could potentially distinguish different two cell states. With the time moving forward, updated cell states according to the genetic network O1 in (d), cell 5 and cell 6 acquires new cell states. At same time, positional information in 4-cell stage increases to 1.0 bit. (g). A diagram illustrating the mechanics of information growth in division-decision systems. The cell contacting matrix (C), together with the current cell expression states (S), form the potential information source manifested in the concatenated matrix ([S,C∙S]). Appropriate genetic networks convert the potential information contained in the concatenated matrix ([S,C∙S]) into the updated expression states (S), resulting in an increase in the positional information. Cell division may change the cell contacting matrix (C), which therefore increases the potential information.

https://doi.org/10.1371/journal.pcbi.1011882.g001

We discovered that only one genetic network, O1, could achieve maximum positional information of 1.0 bit, for the initial cell state ). In contrast, the other eight genetic networks remained at 0.0 bit throughout division and decision (Fig 1C). In addition, all genetic networks maintained 0.0 bit positional information throughout development if the initial cell state is . The O1 genetic network comprises one intracellular self-activation and one intercellular inhibition (Fig 1D), and an example of the gene expression pattern during system growth governed by the O1 was shown in Fig 1E. The development pattern of the system regulated by network O1 consists of three updates in each generation due to the fact that the system can reach stable states via three updates of gene expression in each generation. This self-activation/cross-cell-inhibition connection resembles the lateral inhibition regularly reported in multiple biological patterns [27,40,41], which prompted us to examine the mechanism of information genesis in greater depth.

We observed that the increase in information content in the O1-directed development is due to the update of the contact matrix C induced by cell division. The positional information increases from 0.0 bit to 1.0 bit at step 8, in the 4-cell stage. As shown in Fig 1E, from generation 2 to generation 3, cell 2 and cell 3 divides into cell 4, 5 and cell 6, 7, respectively. These four daughter cells inherit cell states from cell 2 and cell 3, with gene A in the ON state for all four of them. Meanwhile, the action of a division changes the contact relationship between cells: cells 4 and 7 only have one neighbor with gene A ON, whereas cells 5 and 6 have two neighbors with gene A ON (Fig 1F). The contacting state (CS) differentiates four cells at step 7 into two groups (i.e., those contacting one g^A = 1 cell, vs. those contacting two g^A = 1 cells), providing the system with potential information of 1.0 bit. At step 8, regulated by the rules of genetic network O1, the cells contacting two neighbors with high A expressions are strongly inhibited and turn off their gene A expressions. Therefore, cells 5 and 6 update their gene A expression states from ON to OFF. This transition results in four cells with two distinct cell states, so that the system’s positional information reaches 1.0 bit. In other words, the genetic network O1 transforms 1.0 bit of potential information at step 7 into 1.0 bit positional information at step 8 (Fig 1F). In contrast, the other eight genetic networks could not efficiently transfer potential information into the system’s actual positional information (S2A Fig). Usually, the process by which genetic networks convert potential information into positional information has some degeneracy. Multiple inputs in concatenated matrix [S,C∙S] might result in the same output in gene expression states S (Fig 1G). For example, as shown in Fig 1E, from generation 3 to generation 4, the neighbors of cells 10 and 11 are in distinct states. Yet, the degeneracy prevented the genetic network from updating these two cells into two distinct cell states. Genetic networks exhibiting a smaller degree of degeneracy are more capable of converting information.

In addition to Eq 5, we also attempted to control gene expression using additional gene regulatory logics and updating functions [42–44]. We found another genetic network, O2 (Eq S1 and S3 Fig), of the similar topology with O1, has the same performance as O1 (details in S1 Appendix).

Taken together, in the single-gene system, the results above demonstrate that cell division can change the spatial contacting state of a system, offering potential information that can be transformed into positional information with the appropriate genetic network (Fig 1G).

Next, we investigated the principles of information amplification in two-gene systems. With two genes A and B with binary ON and OFF states, a cell’s state can be expressed as , and there are four possible states for each cell: none expressed , both expressed , and only one gene expressed . We used Eq 5 to update the gene expression of the cell. As mentioned in Method, the number of possible intracellular and intercellular regulatory networks (i.e., genetic networks) composed by two genes is 6561 (Fig 2A). Starting from four distinct initial states in a single cell, we simulated all 6561 genetic networks in the division-decision system, and recorded the maximal positional information each network could reach (Fig 2B).

Fig 2. Maximum positional information pattern formation in the two-gene system.

(a). Possible intracellular and intercellular interactions in the two-gene system. The number of possible genetic networks is 6561. The black arrow represents the possible interaction relationship, like an activation, an inhibition or an absence of interaction. Solid and dotted lines represent intracellular and intercellular interaction relationships, respectively. (b). The distribution diagram of the number of genetic networks corresponding to different maximum positional information under four different initial cell states. (c). The exclusive genetic network (T1) in the two-gene system that regulates maximum positional information pattern formation in a stable state. (d). The expression pattern of two genes throughout time and space. This pattern formation regulated by the genetic network T1 in (c). At step 14, system reaches 2.0 bit positional information. (e). Mechanisms for increasing positional information from 1.0 bit to 2.0 bit. When cells 4 to 7 divide, cells 8 to 15 inherit cell states from their mothers. At this time, the potential information in the concatenated matrix ([S,C∙S]) is 2.0 bit. With the time moving forward, updated cell states according to the genetic network T1 in (c), cell 8, cell 11, cell 12 and cell 15 alter their states. At step 14, positional information increased from 1.0 bit to 2.0 bit. (f). Transition linkages across four cell states under the regulation of genetic network T1. The number and color of the circles next to the line indicate the number and state of the cell’s neighbors, and an arrow points from the initial state of the cell to the final state of the cell. The similar cell state transition that took place in (d) is shown by circles with a gray background. (g). The relationship between maximum positional information and the proportion of times a regulatory motif occurs. The proportion of self-inhibition and direct-lateral-activation motif occurrence are negatively correlated with the maximum positional information. While the proportion of direct-lateral-inhibition motif occurrence is positively connected with the maximum positional information. (h). According to the maximum positional information achievable in the development system directed by genetic networks, networks with the same maximum positional information are grouped together. The first column represents the maximum positional information value for each genetic network group. The second column indicates the number of genetic networks in each group. The third column represents the ratio of direct lateral inhibition motif present in the group. The fourth column represents the ratio of indirect lateral inhibition motif present in the group. The fifth column represents the proportion of networks in the group that do not exhibit lateral inhibition motif.

https://doi.org/10.1371/journal.pcbi.1011882.g002

With binary gene activities, the maximum positional information in the two-gene system is 2.0 bit. Our simulations revealed that only two genetic networks, named T1 and T1’, are capable of driving systems to reach this maximal positional information at stable states (Figs 2C, 2D and S4) (S2 Appendix). Since systems can reach stable states via four updates of gene expression in each generation, development patterns of systems regulated by network T1 and T1’ consist of four updates in each generation. Network T1 comprises four regulatory edges: intracellular self-activation of gene A, lateral inhibition between A genes across contacting cells, activation of gene B by gene A across contacting cells, and lateral inhibition between B genes across contacting cells (Fig 2C). The network T1’ can be generated by swapping the A and B genes in the network T1 (S4A Fig). When the initial cell states are or , systems regulated by network T1 reach multicellular stable states with 2.0 bit of positional information (Figs 2C, 2D and S4B). While when the initial cell states are or , systems regulated by network T1’ reach multicellular stable states with 2.0 bit positional information (S4A and S4C Fig). Due to the symmetry between networks T1 and T1’, the ensuing analysis will only focus on network T1. Under the regulation of network T1, the positional information undergoes two leaps during the system’s growth: the first leap occurs at step 10, as the system divides from the 2-cell stage into the 4-cell stage, increasing from 0.0 bit to 1.0 bit; the second leap occurs at step 14, as the system divides from the 4-cell stage into the 8-cell stage, increasing from 1.0 bit to 2.0 bit (Fig 2D).

The second leap of information is influenced not only by the contacting state CS but also by the current gene expression state S. From generation 3 to 4, cells 4 to 7 give birth to cells 8 to 15. Cells 8, 9, 14, and 15 inherit their mother cells’ GREEN state (GREEN state, [A = 1,B = 0]), while cells 10 to 13 inherit the RED state (RED state, [A = 0,B = 1]). Immediately after division at step 13, the concatenated matrix ([S,CS]) is shown in Fig 2: Cells 8 and 9 share the same cell state (GREEN) but have different neighboring states (GREEN vs. GREEN and RED). On the other hand, cells 9 and 10 share identical neighboring states (GREEN and RED), but their cell states differ (GREEN vs. RED). As a result of these differences in neighboring and expression states, the eight rows of the concatenated matrix ([S,CS]) can be divided into four distinct classes with a potential information of 2.0 bit. Under the regulation of T1, at step 14, the states of cells 8, 15, 11, and 12 transition into BLUE and YELLOW. As a result, the system exhibits four distinct cell states with a positional information of 2.0 bit. This transition from two to four states demonstrates how the system converts the combined cell and neighboring states ([S,CS]) into updated positional information. As the differences between cells 9 and 10 are inherited from mother cells, this source of information come from cell lineage.

From the perspective of transition linkages across cell states [45], we observe that physical clustering of cells in the same state plays a crucial role in the emergence of new cell states during the development process in Fig 2D. As shown in Fig 2F, a BLUE cell ([A = 1,B = 1]) transitions to a RED state ([A = 0,B = 1]) when it encounters two neighboring cells with the same BLUE state, which occurs between steps 9 and 10 in Fig 2D. Similarly, a RED ([A = 0,B = 1]) cell transitions to a YELLOW state ([A = 0,B = 0]) when it contacts two cells in RED, which happens between steps 13 and 14. This cascade of new state emergence is facilitated by cell division, in which cells with identical states proliferate to promote the formation of the next new state.

In the two-gene system, we also simulated the growth of systems with different gene regulatory logics and updating functions. Using Eq S1 (S3A Fig), we discovered three genetic networks (T2, T3, and T4) that enabled controlled systems to attain multicellular stable states with 2.0 bit positional information (S5 Fig, See S1 Appendix for details).

Among the 6561 genetic networks we examined, their capacity for increasing positional information during the division-decision process varied. To explore this, we investigated relationships with specific gene regulatory motifs, including direct positive feedback, mutual inhibition, negative feedback, and coherent and incoherent feedforward loops [39,46]. Focusing on six intracellular regulatory motifs and four intercellular regulatory motifs (S6A Fig), we considered only two initial cell states for simplicity. We grouped the genetic networks with the same maximum positional information after running for 5 generations. For each group of genetic networks, we calculated the “motif ratio” by dividing the number of each motif’s occurrences by the total number of genetic networks in the group. We found two types of strong correlations between motif ratio and the maximum positional information (Figs 2G, S6B and S6C). Intercellular direct lateral inhibition positively correlated with maximum positional information, while intracellular self-inhibition and intercellular direct lateral activation negatively correlated (Figs 2G, S6B and S6C). Analyzing genetic networks categorized by their maximum positional information, we observed that direct or indirect lateral inhibition motifs were predominant in networks achieving positional information larger than 1.0 bit (Figs 2H and S6D). These results highlight the crucial role of intercellular lateral inhibition in generating high positional information patterns.

In actual biological processes, various embryonic development stages are initially guided by specific positional cues. In mimicking such initial asymmetry, we also examined the relationship between positional information and regulatory motifs in a one-dimensional multicellular model with one cell fixed at the left end of the system. The fixed cell is not allowed to divide nor change cell state, imposing an asymmetry on the system analogous to the P cell linages in C.elegans development (S7A Fig). After evaluating all possible two-gene networks, we identified two capable of producing a four-striped pattern at the 9-cell stage, and fourteen capable of three distinct stripes (S7D Fig). Under a screened genetic network T5 (S7B Fig), the division-decision system produced a unique pattern of four stripes at the 9-cell stage (S7C Fig). While motif ratios showed no significant correlation with stripe number (S7E Fig), positional information negatively correlated with intercellular direct lateral activation and positively correlated with intercellular direct lateral inhibition motif ratios (S7F Fig). Moreover, all networks achieving positional information higher than 1.0 bit at the 9-cell stable state featured direct or indirect lateral inhibition motifs (S7G Fig). These findings underscore the pivotal role of lateral inhibition in amplifying positional information during development.

For accurate dynamic models in reality, a reasonable basin size is crucial for stable states at a biological level. The basin size represents the total initial cell states evolving into a stable state under specific genetic network regulation. For instance, in network T1, with two cells, there are 16 initial states, and only 10 represent stable multicellular states, suggesting a mean basin size of 1.6. Increasing the cell number to four results in 256 initial states, with a mean basin size of 3.9 (S8A and S8B Fig). Calculations for basin size and positional information at different cell numbers (regulated by T1) revealed a rapid increase in the largest basin size, reaching 30 for 4-cell states and 1275 for 8-cell states (Fig 3A and 3B). The high-information state, evolving through division-decision, shows a decreasing basin size ratio relative to the largest basin: 11 at the 4-cell stage and 71 at the 8-cell stage (Fig 3A and 3B). Basin size and positional information recordings for 4-cell and 8-cell stable states demonstrated that higher positional information correlated with smaller basin sizes (Fig 3C and 3D). Therefore, if the number of cells in the system remains constant, randomized initial states tend to evolve into a state with a large basin size, keeping positional information low. However, as the system grows from a single cell, early stages provide a larger basin ratio for the high-information state, “escorting” the system into its small basin through division. The pedigree procedure’s cell division prevents the system from evolving into a stable state with limited positional information by eliminating numerous unstable starting states (Fig 3E).

Fig 3. Basin size and stability of multicellular stable states.

(a). The basin size of 4-cell stable states and the distribution diagram of the corresponding number of stable states under the regulation of genetic network T1. The red dotted box and red arrow denote the position of the 4-cell stable state formed by the growing system regulated by the genetic network T1. (b). In logarithmic coordinates, the basin size of 8-cell stable states and the distribution diagram of the corresponding number of stable states under the regulation of genetic network T1. The red arrow denotes the position of the 8-cell stable state formed by the growing system regulated by the genetic network T1. (c). Basin size and positional information corresponding to distinct 4-cell stable states under the regulation of genetic network T1. (d). In logarithmic coordinates, basin size and positional information corresponding to distinct 8-cell stable states under the regulation of genetic network T1. (e). Under the regulation of genetic network T1, multicellular stable states corresponding to various basin sizes were displayed in each generation in the schematic diagram. The areas of the smaller, colored circles represent the basin sizes of the distinct multicellular stable states, while the large gray circle represents the cumulative basin size of all stable states for each generation. And the circle’s color corresponds to the positional information of the multicellular stable state. Red lines indicate stable states that appeared in the one-dimensional division-decision growing system under the regulation of genetic network T1. (f). To make the cell state update with a certain probability of mistake, noise is introduced to the growing system. The noise that occurred during cell state update is depicted on the left panel. And the relationship between the probability of cell state update mistakes and pattern accuracy at various cell cycle lengths was illustrated in the right panel. (g). The left panel illustrates a representative instance of asynchronous cell division, while the right panel delineates the relationship between the standard deviation (σ) of cell cycle lengths and the accuracy of patterning under conditions of asynchronous division.

https://doi.org/10.1371/journal.pcbi.1011882.g003

Noisy cell state decision involves erroneous gene expression updates with a certain probability (Fig 3F), while differences in cell cycle length contribute to a divergence in the timing of cell divisions, leading to an asynchronous division pattern among cells within the same generation (Fig 3G). We simulate the division-decision system 10,000 times under these two distinct conditions for various genetic networks. Identifying the standard multicellular stable state corresponding to the maximum positional information (2.0 bit) in a noise-free and synchronous environment, we assess how likely systems can reach this state under conditions characterized by stochastic perturbations or asynchrony in cell division. The “accuracy” of the system, determined by the fraction of simulations reaching the standard state given gene expression variability and cell cycle length diversity, decreases for both forms. However, the impact differs slightly: in the case of the first condition, accuracy decreases as the cell cycle length increases, suggesting more mistakes with longer cell cycles. Conversely, for the second condition, where cell division is asynchronous, longer cell cycles enhance accuracy as they provide cells more time to rectify aberrant states. This pattern holds for other screened genetic networks using Eq S1, and T1-regulated systems exhibit better resistance to both conditions than other genetic networks (S8C and S8D Fig).

In our exploration of increasing positional information, we initially employed a basic Boolean model for gene expression updates. However, concerns about potential artificial outcomes prompted us to replace the Boolean model with an ordinary differential equation in this section [47,48]. This substitution aims to assess the continued applicability of the principle of increased positional information within a continuous model.

In the single-gene system, based on the network O1 (Fig 1D) and O2 (S3B Fig), we constructed different ordinary differential equations, by adopting AND or OR logics in combinatory regulations (Fig 4A and 4B, See Eqs S2-S3 in S1 Appendix). When the system reached the 4-cell stage (generation 3), we observed that the stable expression levels of gene A were different across cells, indicating the increase of the positional information in the system and the embellishment of differential cell states (Fig 4A and 4B). Gene A expression, influenced by neighboring cells, led to varied stable states in cells with differing degrees of inhibition, consistent with Boolean model findings (S9A and S9C Fig). Moreover, transforming gene A’s intracellular self-activation to constitutive expression aligned with the Boolean model’s network O2 (Figs 4C, S3B and Eq S4). Despite lacking self-activation, the system still exhibited spontaneous pattern formation, highlighting the significance of the lateral inhibition motif (Fig 4B and 4C, see S1 Appendix for detail).

Fig 4. The spatiotemporal expression patterns of genes in the continuous model and the noise resistance of particular patterns.

(a). The AND-gate regulatory function of the genetic network O1 in the continuous model and the spatiotemporal expression pattern of gene A simulated by this function. (b). The OR-gate regulatory function of the genetic network O1 in the continuous model and the spatiotemporal expression pattern of gene A simulated by this function. (c). The OR-gate regulatory function of the genetic network O2 in the continuous model and the spatiotemporal expression pattern of gene A simulated by this function. (d). Resistance to noise (λ) of patterns generated in generation 3 under various regulatory settings in single-gene systems. (e). The spatiotemporal expression patterns regulated by the genetic network T1 in the continuous model, and corresponding equations which regulate this pattern formation. (f). Resistance to noise (λ) of patterns generated in generation 4 under various regulatory settings in two-gene systems.

https://doi.org/10.1371/journal.pcbi.1011882.g004

Analyzing three genetic networks under continuous expressions revealed that noise resistance varied (S1 Appendix). In generation 3, the system regulated by constitutive expression showed the highest noise resistance. In network O1, when intracellular self-activation and intercellular lateral inhibition are integrated through an AND-gate instead of an OR-gate, the pattern generated is more noise-resistant (Fig 4D). All networks experienced decreased noise resistance in generation 4 (S9F Fig).

We also applied the continuous model to simulate the two-gene growing system. By scanning the parameter space for network T1, we found that the stable expression pattern of genes A and B in generation 4 would match that of the Boolean model in the continuous model if gene A was AND-gate and gene B was AND-gate or OR-gate (S1 Appendix). The simulation results in pattern formation are consistent with the Boolean model (Figs 4E and S10). However, the continuous adaptation of the T1 model did not mirror the bi-stability pattern exhibited by gene B in the Boolean model, unless an intracellular self-activation was introduced (see S1 Appendix). This modified network, referred to as T6, showcased even more distinctive cell states in generation 4 (see S11D Fig). Nevertheless, when it comes to noise resistance, T1 demonstrated superior performance compared to T6 (Fig 4F).

We sought to investigate the applicability of our concise mathematical framework to biological phenomena, particularly in relation to the development of C. elegans. This organism has long been recognized as a model for studying developmental processes due to its highly predictable lineage specifications that occur independently of external cues. To assess the accuracy of our framework, we utilized published experimental data from studies investigating cell-cell contact relationships and gene expression patterns during the early stages of C. elegans embryo development [18,49–53]. Using this available data, we digitalized the cellular states by binarizing their gene expressions, then performed calculations to determine the positional information and potential information at different time points during C. elegans embryo development, starting from the fertilized egg up until the 200-cell stage (Fig 5A, See S1 Appendix for details). Remarkably, our calculations revealed that the embryo possesses more than 5.0 bits of positional information by the time it reaches the 200-cell stage (See S1 Appendix for details). Notably, the positional information exhibits a clear and incremental trend, closely following the increases in potential information as development progresses (Fig 5A). This evident and consistent rise in positional information throughout development, without external sources, demands mechanistic explanations.

Fig 5. Examples of positional information increase during C.elegans embryonic development.

(a). Positional information and potential information are estimated at each step of the development of the C.elegans embryo, from the fertilized egg to the 200-cell stage. The cell state is defined by the binarized gene expression profiles. Different states are marked by different colors, and named by the lower case of the first cell entering this state. (b). ABa cell transitions from the old cell state “ab” to the new cell state “aba” because it is not in contact with P2 cell during the 4-cell stage. (c). At the 12-cell stage, cell ABarp inherited the cell state “aba” and changed to “abarp” upon contact with C cell. Although the cell ABpra is also in contact with the C cell, its inherited cell state is “ab”, so its cell state does not change.

https://doi.org/10.1371/journal.pcbi.1011882.g005

During the incipient stages of C. elegans embryo development, evidence suggests that neighboring cells provide cells with potential information (Fig 5B). Also, the significance of asymmetric cell divisions cannot be ignored, particularly in how cell polarity established prior to division. Our model integrates these elements by accounting for the inherent cell polarity and examining the influence of intercellular interactions on fate determination post-division. Specifically, at the 4-cell stage, cells ABa and ABp, progenies of progenitor cell AB, inherit the same “ab” cell state and equally express the Notch receptor. The anteroposterior disposition of ABa and ABp is contingent upon the polarity that progenitor cell AB establishes before division. Concurrently, the caudal cell P2 expresses the Delta ligand. Prior to the signaling event, the four-cell system exhibits three cell states, representing 1.5 bit of positional information. Meanwhile, cell P2 interacts with cell ABp, but not with cell ABa. Therefore, in a concatenated matrix that includes the states of each cell and their neighboring cells, there is 2.0 bit of potential information. During the signaling event, the Notch-Delta signal pathway prevents cell ABp from activating a specific set of genes associated with the “aba” state. As a result of lacking contact with cell P2, cell ABa transitions from the “ab” cell state to the “aba” cell state. Subsequently, all four cells demonstrate different gene expression states, resulting in the embryo possessing 2.0 bit of positional information. Hence, the spatial contact between cells conveys potential information regarding this event. This also underscores the interplay between intrinsic asymmetric divisions and extrinsic contact-mediated signaling in the developmental process.

In the 12-cell stage of C. elegans development, a noteworthy example involves the cells ABala, ABalp, ABara, and ABarp, all inheriting the common cell state “aba” from the ABal cell. The subsequent divergence in cell fate arises from the distinct interactions each cell has: ABalp interacts with MS, while ABarp interacts with the cell C. Consequently, ABalp adopts the “abalp” cell state, while ABarp transitions into the “abarp” cell state (depicted in Fig 5C). This process exemplifies the acquisition of positional information by the difference in neighboring cell states.

Simultaneously, another fate specification demonstrates how lineage inheritance contributes to the system’s acquisition of positional information. Descendants from ABar, including cells ABpla, ABplp, ABpra, and ABprp, all share the “ab” cell state due to their lineage connection with ABar. While cells inheriting the “aba” cell state transform into the “abarp” cell state upon contacting cell C, as exemplified by ABarp, those inheriting the “ab” state remain unchanged even in contact with cell C, such as cell ABpra (Fig 5C). In this scenario, the positional information distinguishing ABarp and ABpra arises not only from the state of neighboring cells but also from the inherent state of the cells themselves, inherited through their maternal lineages.

We developed a one-dimensional multicellular model in which all cells were derived from a single ancestor cell and shared identical genetic codes for cell state decisions (Fig 1A). The gene expression profiles of each cell were influenced by the states of its neighboring cells. Additionally, cell division occurred over time, and daughter cells inherited the expression states of their mothers. The system was designed to simulate the emergence of spatial patterns in one-dimensional space as the system grows.

We enumerated genetic networks consisting of one and two genes. In the single-gene network, the state of a cell was determined by only one gene within the cell: . We considered all possible one-gene regulations, including all possible intracellular and intercellular interactions, as described by each of its two edges (Fig 1B). There is one edge in the intracellular interaction and one edge between genes in neighboring cells. In this binary model, each of these edges can be an activation interaction (1), an inhibition interaction (-1), or an absence of interaction (0) between the two genes. There are therefore 3² = 9 possible types of genetic networks (S2A Fig). Based on the aforementioned division-decision model, we simulated this dynamic system starting from a single cell under all possible one-gene networks.

In the two-gene system, the state of one cell was determined by the expression profile of two genes in the cell: . We also considered all possible two-gene networks (Fig 2A). There are eight edges among four genes in two adjacent cells, and each edge can indicate an activation (1), inhibition (-1), or absence (0) of interaction between the two genes it is associated with. These forms 3⁸ = 6561 possible genetic networks. Similarly, we simulated this two-gene division-decision system under all possible genetic networks.

At each step of the system’s growth, we calculated the positional information corresponding to the current cell states (i = 1, 2, …, m; m represents the current number of cells in the system). Positional information (PI(state,position)), defined as the reduced uncertainty of gene expression states given the cellular position, quantifies both diversity and reproducibility [16].

(Eq 4)

In Eq 4, “state” refers to the gene expression profile of the cell, denoted by , while “position” refers to the spatial location of the cell. H(state) term represents the diversity of cell states in the system. Shannon entropy is used to measure the uncertainty or randomness of a probability distribution. In this case, H(state) is calculated based on the states of m cells, where each cell has a gene expression profile denoted by . H(state|position) term represents the conditional entropy of the gene expression states given the cellular position. Conditional entropy measures the remaining uncertainty of a random variable (cell state) given the knowledge of another random variable (position). In this context, H(state|position) captures how much uncertainty in cell states can be reduced by knowing the spatial location of the cells. By subtracting H(state|position) from H(state), we obtain the positional information (PI) for the current system. PI(state, position) quantifies the reduced uncertainty or increased information about cell states when considering the cellular position. It represents both the diversity and reproducibility aspects of the system.

S1A Fig depicts two scenarios of calculating positional information for a simple one-dimensional system with 24 cells, each having 4 potential states. To measure the uncertainty for the two systems we assess the frequency of each state and calculate the total entropy (H(state) = 2.00) bit for the left system and 1.98 bit for the right systems). In the left system, since the gene expression levels within cells are rigidly governed by regulatory rules without any variance from internal or external noise, cell states are entirely determined by their positions, resulting in a conditional entropy of 0 bit. Conversely, the right system incorporates variability through both regulatory rules and stochastic influences, manifesting in reduced consistency of cell states at identical positions and a non-zero conditional entropy of 0.68 bit. The positional information for each system is then determined by subtracting the conditional entropy from the total entropy. Consequently, the first system has 2.0 bit of positional information, while the second system has 1.32 bit.

In addition to positional information, we also defined the concept of “potential information”. At each step of the system update, the concatenated input matrix [S,C∙S] in Eq 3 provides each cell with the cell-specific input for cell state decisions. Therefore, the mutual information between [S,C∙S] and cellular position is defined as “potential information”, as it offers potential differences between cells. This potential information needs to be processed by the identical genetic network A to finally update the expression profiles.