How Criticality of Gene Regulatory Networks Affects the Resulting Morphogenesis under Genetic Perturbations

Many complex biological structures, including tissues and organs, are formed through a developmental process from a single fertilized cell to multicellular embryos [5]. The morphogenesis of those complex patterns is driven by a gene regulatory network (GRN) that exists within each cell and responds to cell-cell interactions [35, 36]. Random Boolean networks (RBNs) were proposed by Kauffman as a theoretical model of such GRNs [28]. In RBNs, genes (nodes) have binary states (either on or off) whose dynamics are determined by a set of Boolean functions over the states of other genes. Although RBNs are a highly simplified model, they have been extensively utilized in artificial life and complex systems research [1, 3, 21, 22, 38, 45, 55, 60].

In the context of GRNs, the concept of criticality of RBNs has been discussed as a phase transition point between ordered and chaotic regimes for the dynamics of those networks [29, 30]. The criticality of GRNs has been recognized as a property that makes robustness and adaptability coexist in living organisms [2]. When perturbations are added to GRNs, ordered GRNs are so robust that they just sustain existing cellular functions. On the contrary, chaotic GRNs are so adaptable that they vigorously create new functions rather than conserving existing ones. In between, critical GRNs stably sustain their functions against perturbations, and at the same time flexibly generate new phenotypes, which help organisms to adapt to new environments because they have optimal balance between robustness and adaptability.

Whereas many studies have been performed to elucidate the relationship between criticality of GRNs and dynamics of GRNs at the single-cell level based on RBNs [6, 43, 47, 52, 53, 56], the relationship between their criticality and the properties of organisms at a higher level has not been fully explored. Only a few studies have determined how the properties of intracellular GRNs influence properties at a multicellular level, using RBNs as GRNs [12, 13, 18, 40, 58]. They all used multicellular models where all the cells have the same RBNs in discrete spaces like cellular automata. Such a model is not realistic as a morphogenetic model representing the development from a single cell to a multicellular embryo. Moreover, research on the relationship between the criticality of GRNs and properties of multicellular systems under genetic perturbations (e.g., mutations) has not been conducted yet, although mutations do occur in cells of living organisms by stochasticity or through environmental factors [4, 7–9, 11, 49, 54].

To investigate a potential role of criticality of GRNs in morphogenesis, which is hard to uncover through the single-cell-level studies, we have recently proposed GRN-based morphogenetic systems operating in a continuous space and found that the criticality of GRNs facilitated the formation of nontrivial morphologies [31, 32]. To make our finding more solid from an evolutionary viewpoint, in this study we include an evolutionary process in our morphogenetic model by introducing genetic perturbations into intracellular GRNs. We add perturbations to the GRNs by adding, deleting, or switching a regulatory link between a pair of genes. We look into whether the role of the criticality of GRNs can be maintained even in the presence of the evolutionary perturbations. Also, we examine what the resulting morphologies are like and what kind of biological implications this model can offer from the epigenetic viewpoint in morphology.

The rest of the article is structured as follows. In Section 2, we describe our GRN-based morphogenetic model. In Section 3, parameter setting for simulations and various measures for the analysis of results are explained. In Section 4, morphologies obtained from the simulations of the morphogenetic model are analyzed. Section 5 summarizes and concludes the article.

We have developed a computational model of morphogenetic processes of cell aggregation, in which all the cells have an identical intracellular GRN. Figure 1 shows the algorithm for our model. The simulation starts with one seed cell. It imitates the process in which a single zygote divides and grows into multicellular form during embryonic development. In the cell, a RBN is produced as an intracellular GRN. Neighboring cells are detected within a fixed neighborhood radius. Through the interaction with neighbors, cells' fates are determined by the GRN. We assume that there are four fundamental cell fates in our model: proliferation, differentiation, apoptosis, and quiescence. Cells expressing proliferation, differentiation, or quiescence can switch their fates through cell-cell interaction. The cells are positioned in a two-dimensional continuous space by spring-mass-damper (SMD) kinetics. Until the termination condition of the simulation is satisfied, the initial seed cell grows into an aggregation iterating the processes of finding neighboring cells and positioning cells in the space in each time step. The simulator of our model was implemented in Java.

A RBN (a.k.a. NK Boolean network) was suggested as a GRN model by Kauffman [28–30]. Here N is the number of nodes, and K is the number of input links per node. A node represents a gene. The state of a node can be either on (1, activated) or off (0, inhibited). The node state is determined by the states of input nodes and a Boolean function assigned to each node. A state space that is constructed from the topology of a RBN and assigned Boolean functions refers to the set of all the possible configurations and all the transitions among them. Figure 2 shows schematic diagrams for an example RBN and its state space. In the state space, stationary configurations are called attractors, and the others are called basins of attraction of the attractors. The dynamics of RBNs are divided into three regimes, depending on the structure of their state space: ordered, critical, and chaotic. Using the node indegree (K), internal homogeneity (p), or canalizing functions, the dynamics of RBNs can be systematically varied. The dynamics of a RBN is known to be determined by K, namely, K = 1 is ordered, K = 2 is critical, and K > 2 is chaotic, on average [29, 30]. For our morphogenetic model, we use a RBN that consists of 16 nodes (N = 16) as an intracellular GRN. Adjusting K from 1 to 4, we vary the dynamics of RBNs: ordered (K = 1), critical (K = 2), and chaotic (K = 3, 4).

In view of in vitro experimental data showing that attractors of GRNs represent cell types or cell fates, Huang suggested a conceptual framework to explain stochastic and reversible transitions between cell fates, using NK Boolean networks [10, 25–27]. Our morphogenetic systems are based on his framework. We randomly assign the four cell fates to attractors of GRNs. Specifically, if there is only one attractor, proliferation is assigned to the attractor. If there are two attractors, proliferation and differentiation are randomly assigned to the two attractors. Likewise, if there are three attractors, proliferation, differentiation, and apoptosis are randomly assigned to those attractors. If there are four or more attractors, proliferation, differentiation, and apoptosis are randomly assigned to three attractors and quiescence is assigned to the rest of the attractors (Figure 3).

With regard to cellular behaviors, cells in proliferation are divided into two, and the daughter cells are placed within a fixed neighborhood radius (R) centering on the mother cells. Cells in differentiation are labeled as differentiated. Cells in apoptosis die and disappear. Cells in quiescence do not show any behaviors.

Switching between cell fates occurs by perturbations of internal gene expression values of an intracellular GRN through cell-cell interaction. Our mechanism for cell-cell interaction is based on cell signaling of Damiani et al.'s multiple random Boolean networks on 2D cellular automata [12, 13]. In our model, an intracellular GRN has n genes, which consist of normal genes (g) and special genes (r) as shown in Figure 4a. The special genes exist in pairs where the genes producing signaling molecules (r₁) and receptors (r₂) are matched one to one. This one-to-one correspondence indicates signal transduction specificity by which certain signaling molecules respond to particular receptors. The special genes r₁ synthesize signaling molecules and release them. Then, those molecules bind to the corresponding receptors r₂ on cells within the neighborhood radius R.

The signal transduction has two mechanisms: autocrine and paracrine. Autocrine means a cell produces signaling molecules that bind to receptors on the same cell. Paracrine means a cell produces signaling molecules that bind to receptors on its neighboring cells. Figure 4b explains the two mechanisms. In our model, when there are no neighboring cells, autocrine is used. When there are neighbors, paracrine is used.

The gene expression values of an intracellular GRN are updated as follows:

• • 

  Normal genes g: The states of the normal genes are updated by the states of input nodes and randomly assigned Boolean functions.

• • 

  Genes producing signaling moleculesr₁: Like the normal genes g, the states of the genes producing signaling molecules r₁ are updated by the states of input nodes and randomly assigned Boolean functions. If the states of r₁ are 1, the genes produce signaling molecules. If the states are 0, they do not.

• • 

  Receptors r₂: The states of the receptors r₂ are determined by the average concentration of the signaling molecules within the neighborhood radius R. Figure 5 shows an example of calculating the average concentration of the signaling molecules from the neighboring cells to determine the state of the receptor gene 2 of cell i. Based on a certain threshold (τ_th), the state of the receptor gene 2 is updated. If the average concentration value is bigger than τ_th, the state becomes 1. Otherwise, it becomes 0.

The following steps are taken to update the gene expression values and to determine a cell fate from the change of the gene values:

• (1) 

  Look into whether there are neighboring cells within the neighborhood radius R. If there are neighbors, paracrine signaling is used. Otherwise, autocrine signaling is used.

• (2) 

  Determine the states of the receptors r₂ through the average concentrations of signaling molecules within the neighborhood radius R and the threshold value τ_th.

• (3) 

  Activate or inhibit genes that have the receptors as input nodes in an intracellular GRN based on the states of the receptors r₂. If the states of the receptors are 1, the states of genes become on (1, activated). Otherwise, the states become off (0, inhibited).

• (4) 

  Check the attractor that the updated gene states finally converge to.

• (5) 

  Express the cell fate that is assigned to the attractor for cellular behaviors.

• (6) 

  Assign the states of the attractor as gene expression values for the next time step.

To acquire much more diverse spatial patterns, we introduce the dependence of parameters k, l, and c on cell fates and perturbations to the cell position (x, y). In the case of k, l, and c, we determine the values depending on six types of cell fates possible between two cells: [proli-proli], [proli-diff], [proli-qui], [diff-qui], [diff-diff], and [qui-qui], where proli is proliferation, diff is differentiation, and qui is quiescence (Figure 6). Here apoptosis is excluded because cells due to apoptosis disappear from the space. In each simulation run, the values of k, l, and c are randomly chosen in certain ranges given in Table 1. For the perturbations, small perturbation values are added to the updated cell positions.

We performed 10,000 independent simulation runs for each value of K, that is, K = 1, 2, 3, and 4. Specifications of parameters for the simulations were the following:

• • 

  Space: The cells were positioned in a 2D continuous 700 × 700 (in arbitrary units) square area.

• • 

  Limitation of cell population: The population growth was limited to 200 cells to keep computational loads reasonable in each run.

• • 

  Simulation termination condition: The simulations were terminated when the time step t reached 1,000 or there was no cell remaining in the space due to apoptosis.

• • 

  Parameter values: The values of parameters concerning GRNs, cell-cell interaction, and cellular motion are given in Table 1. The number of special genes (r) was determined according to Damiani et al.'s model [13] and biological evidence [44]. Damiani et al. showed that the diversity of cellular behaviors was maximized when the coupling strength (fraction of genes that are affected by neighboring cells) was around 0.1 [13]. This value is also similar to the ratio of the number of genes related to cell signaling to the number of human genes [44]. In our model, the coupling strength is set to 0.125, because the number of special genes (r) is 2 while the number of nodes of a GRN (N) is 16.

As mentioned in the Introduction (Section 1), we previously presented a GRN-based morphogenetic model and revealed the role of the criticality of GRNs in morphogenesis [31, 32]. To make our finding more solid from an evolutionary viewpoint, we have included genetic perturbations in our morphogenetic model in this study. Specifically, we assumed that the genetic perturbation was due to a germinal mutation occurring in a pre-zygotic cell, which is a small-scale mutation at a genetic level. The germinal mutation is passed on to offspring, and it is present in all resulting cells during embryo development [19, 23]. We perturbed the intracellular GRN of a seed cell at the initial time step by adding, deleting, or switching one regulatory link between a pair of genes [15, 24, 37]. Because cells were duplicated through the process of cell division from the perturbed GRN in the seed cell, all the cells composing a morphogenetic pattern had the same perturbed GRNs.

Such a small regulatory link perturbation did not significantly change the average number of input links per node (K). For example, if K = 2, the total number of links of a GRN is 32, because the node size is 16. If one regulatory link is deleted as a genetic perturbation, the GRN consists of 31 links, making the value of K 1.94.

Based on the studies showing that biological systems have robustness and evolvability [14, 33, 42, 46, 57, 59], we focused on analyzing robust and evolvable GRNs among the perturbed GRNs. In our model, if the GRN conserved its existing attractors and created new attractors simultaneously after perturbation, we considered the GRN a robust and evolvable GRN (Figure 7) [2].

To compare how morphogenetic patterns are different between groups K = 1, 2, 3, and 4, we used the 12 measures (i–xii) described below [50, 51]. Among those measures, vi–xii are measures regarding network topology. To apply them to our morphologies, we constructed a network from each morphogenetic pattern by connecting each cell to other cells within the neighborhood radius R. Figure 8 shows an example morphology and a network constructed from it using our network construction method. Such network construction allowed for detection of topological differences more effectively. All the morphogenetic measures were obtained from the final configuration of each simulation.

• Number of cells (numOfCells). This is the total number of cells in a morphogenetic pattern.

• Average distance of cells from center of mass (massDistance). This is the mean of Euclidean distances between each cell position and the center of mass (⁠$x¯$⁠, $y¯$⁠), that is, the point with the average coordinates of all the cells.

• Average pairwise distance (pairDistance). This is the mean of Euclidean distances between two randomly sampled cells' positions. The mean was calculated based on 10,000 pairs, which were sampled with replacement.

• Kullback-Leibler divergence between pairwise particle distance distributions of a morphogenetic pattern and a random pattern (kld). This quantifies how nontrivial morphogenetic patterns are, compared to randomly distributed patterns. It was calculated as the Kullback-Leibler (KL) divergence between pairwise particle distance distributions of a morphogenetic pattern (Figure 9a) and a randomly distributed pattern (Figure 9b) made of the same number of cells within the same spatial dimensions. Each pairwise particle distance distribution was obtained through 10,000 random samplings with replacement of a pair of coordinates of cells (Figure 9c). Thus, the larger kld is, the more structured (nonrandom) the morphogenetic pattern is.

  Both pairDistance and kld used pairwise particle distances. pairDistance measures a rough size of a morphogenetic pattern, while kld quantifies the nontriviality of its morphology. Two morphogenetic patterns may have similar pairDistance values but very different kld values at the same time.

• Mutual information of cell fates between neighboring cells (MI). This examines nonlinear correlation of cell fates between neighboring cells in a morphogenetic pattern. It was calculated using the frequencies of three neighboring cell fates (except for apoptosis, because cells expressing apoptosis die and disappear from the space). Figure 10 shows an example calculating MI in a morphogenetic pattern. Counting the frequencies of combinations of fates of neighboring cells, MI captures how much informational correlation would exist between the fate of a cell and that of its neighbors. If there was only one cell remaining, the value of MI was set to 0.

• Average clustering coefficient (avgCluster). This explains how densely connected the nodes (cells) are to each other in a network. The clustering coefficient C(i) of node ν_i in a network is defined as follows:

  $Ci=ejk:νj,νk∈Ni,ejk∈Ekiki−1/2.$

  Here, e_jk is a link that connects node ν_j and node ν_k within the set of neighboring cells N_i around a node ν_i, E is a set of links in the network, and k_i is the degree of ν_i (i.e., the size of N_i). The denominator is the total number of possible node pairs within node ν_i's neighborhood. The numerator is the number of actually connected node pairs among them. The average clustering coefficient is given by

  $C¯=∑icin,$

  where n is the total number of nodes.
• Link density (linkDensity). This describes the density of connections in a network. For a network G composed of nodes n and m links, the link density ρ(G) is given by

  $ρG=mnn−12=2mnn−1$

  for an undirected network, where m is the number of links.

• Number of connected components (numConnComp). This is the number of connected components in a network. A connected component is a subgraph in which there is a path between every pair of nodes. A single isolated cell was also considered one connected component by itself.

• Average size of connected components (meanSizeConnComp). This is the mean of the numbers of nodes in each connected component in a network. When there was no connected component, the value was set to 0.

• Homogeneity of sizes of connected components (homoSizeConnComp). This quantifies how similar the sizes of connected components are in a network. This measure was calculated as one minus the normalized entropy in the distribution of sizes of connected components. Figure 11 shows an example calculating homoSizeConnComp in a morphogenetic pattern. When there was only one connected component, the value was set to 1.

• Size of the largest connected components (sizeLarConnComp). This is the maximum size of the connected components in a network.

• Average size of connected components smaller than the largest one (meanSizeSmaller). This is the mean of the sizes of connected components except for the largest connected component in a network. If there was only one connected component in the network, the value was set to 0.

To investigate the relationship between intracellular GRNs and expressed cell fates, we calculated the basin entropy and cell fate entropy from the sizes of basins of attractions and cell fates distributed in a morphogenetic pattern. We thought that the numbers of actually expressed cell fates in a morphology might be proportional to the basin sizes of attractors where each cell fate was assigned. We calculated basin and cell fate entropies to look into whether or not our expectation would be correct. As in the computation of MI, only three cell fates (proliferation, differentiation, and quiescence) were considered (i.e., apoptosis was ignored).

• Basin entropy:

  $Hbasin=−∑ρPρ⋅log2Pρ,$

  where P_ρ is the size of the basin of the attractor ρ (to which proliferation, differentiation, or quiescence was assigned) divided by the sum of sizes of all the basins except for the basin size of the attractor for apoptosis. Thus,

  $∑ρPρ=1.$

  Originally, basin entropy was suggested by Krawitz as a measure of the complexity of information that a system can store in NK Boolean networks [34]. We used it as a measure to examine the versatility of the three cell fates (proliferation, differentiation, quiescence).
• Cell Fate Entropy:

  $Hcellfate=−∑fPf⋅log2Pf,$

  where P_f is the number of cells expressing a cell fate f (proliferation, differentiation, quiescence), divided by the numbers of cells (except for those in apoptosis) in a morphogenetic pattern at the final time step. Hence,

  $∑fPf=1.$

  In the case that there were no cells expressing a fate (proliferation, differentiation, quiescence), its log value was set to 0. Also, when there was no cell in the space, the value was set to 0.

Figure 12a shows probabilities of producing robust and evolvable GRNs against perturbations for K = 1, 2, 3, and 4 in 10,000 simulation runs. We found that robust and evolvable GRNs were generated with the highest probability at K = 2. Figure 12b shows samples of visualized morphogenetic patterns produced by robust and evolvable GRNs for each K.

We calculated the 12 morphological measures of morphogenetic patterns generated by robust and evolvable GRNs. Because the robust and evolvable GRNs were generated with different probabilities for different values of K, we applied bootstrap sampling 1,000 times to the values of the 12 measures for comparison between groups (K = 1–4) with unequal sample sizes. Figure 13 indicates the comparison of means between groups for the measures, where Kruskal-Wallis and Nemenyi (as post hoc analysis) tests were performed to show statistically significant differences among the groups.

Furthermore, we produced a correlation matrix to investigate correlations between the 12 measures (Figure 14). We found that the following six measures were highly correlated with numOfCells:avgCluster, massDistance, meanSizeConnComp, meanSizeSmaller, pairDistance, and sizeLarConnComp. These correlations were found in Figure 13 as well. The values of numOfCells at K = 1, 2, 3 were similar but fell sharply at K = 4. This trend was also shown in the six measures highly correlated with numOfCells. Meanwhile, kld, MI, homoSizeConnComp, linkDensity, and numConnComp showed different trends. In Figure 13c, kld was highest at K = 2, which means that nontrivial morphogenetic patterns were generated most frequently when the GRNs were critical under the genetic perturbations. This result demonstrates that the role of criticality of GRNs is maintained even in the presence of evolutionary perturbations.

In addition, from MI, homoSizeConnComp, linkDensity, and numConnComp, we found two interesting properties of the nontrivial morphologies at criticality. Firstly, certain combinations of cell fates between neighboring cells occurred most frequently. In Figure 13l, MI was highest at K = 2. It indicates that the fate of a cell is strongly correlated with the fate of its neighboring cells in a morphogenetic pattern generated at the criticality. To examine the relationship between intracellular GRNs and expressed cell fates, we measured basin entropy and cell fate entropy (Figure 15). Our original expectation was that if the basins of attraction for the three cell fates were most evenly distributed at K = 2, the expressions of different cell fates would be maximally balanced in a morphogenetic pattern. However, the cell fate entropy was highest at K = 1, although the basin entropy was highest at K = 2. This means that the distribution of cell fates in a morphology was not a simple reflection of the basin sizes of a GRN at a single-cell level, but more like an emergent property at a multicellular level obtained through the developmental process involving cell-cell interactions.

Secondly, the nontrivial morphologies emerged typically in topologically homogeneous cell clusters. In Figure 13b, d, h, K = 1 showed relatively high homoSizeConnComp, low linkDensity, and low numConnComp values, on average, compared to the corresponding measures at K = 2, 3, and 4. Here we simply express the observations qualitatively as (high, low, low). Similarly, in the same order, K = 2 showed (high, high, high), K = 3 showed (low, high, high), and K = 4 showed (low, low, low). These can be interpreted as follows: The morphologies at K = 1 consisted of homogeneous large-size connected components, which had a shape like a long chain. The morphologies at K = 2 were composed of homogeneous-size connected components where cells were interconnected. The morphologies at K = 3 were composed of heterogeneous-size connected components where cells were interconnected. The morphologies at K = 4 were composed of heterogeneous-size small connected components that had a shape like a short chain. Figure 16 summarizes the typical topological properties for K = 1–4 schematically.

In our morphogenetic model, there is a feedback relationship (Figure 17). Interactions with neighboring cells determine a cell fate. Depending on the cell fates, the parameters of SMD kinetics are applied. The cells are positioned by SMD kinetics. The positions of cells influence the number of neighboring cells. In the feedback relationship, we found that the nontrivial morphologies were produced most frequently when the GRNs were critical under the genetic perturbations. Besides, the nontrivial morphologies at criticality had the most frequent occurrence of certain combinations of cell fates between neighbors, and were composed of topologically homogeneous cell clusters. Because the parameter values of SMD kinetics determining cells' positions depended on the cell fates, the more frequent those combinations of cell fates between neighbors were, the more likely the same SMD parameter values were to be applied among the cells. Thus, the most frequent combinations of cell fates between neighbors would naturally produce more homogeneous-size connected components where cells were interconnected. Such spatial arrangements of multicellular patterns due to the criticality of GRNs are not easily predictable from the single-cell level. This justifies the value of the present study on the role of criticality of GRNs in the context of multicellular settings.

In this study, we introduced genetic perturbations into our morphogenetic model, which uses Kauffman's RBNs as intracellular GRNs, and SMD kinetics for cellular movements. We looked into whether the previously reported role of the criticality of GRNs could be maintained even in the presence of evolutionary perturbations, and what the resulting morphogenesis would look like. We found that nontrivial morphologies were generated most frequently when the GRNs were critical under the genetic perturbations, which is consistent with the previous result obtained for morphogenetic systems without evolutionary perturbations. Moreover, we found that the nontrivial morphologies at criticality tended to be made of topologically homogeneous cell clusters due to the spatial arrangements in which certain combinations of cell fates between neighboring cells occurred most frequently. Based on these findings, we conclude that the criticality of GRNs facilitates the formation of nontrivial morphologies by adjusting the spatial arrangements of cells in GRN-based morphogenetic systems, even under the genetic perturbations.

Our findings have implications from an epigenetic viewpoint. Researchers in epigenesis have suggested that heterogeneous and complex features emerge from homogeneous and less complex components through the interactions among them [41, 48]. In our model, we showed that the nontrivial morphologies were produced most frequently at criticality, typically with topologically homogeneous cell clusters. Thus, the result not only supports the theory of epigenesis in developmental biology, but also implies that highly structured tissues or organs in morphogenesis of multicellular organisms might stem from cell aggregations with critical GRNs.

The present study has limitations. The properties of morphogenetic patterns were explored only using our artificial model based on RBNs as GRNs. To make our findings more relevant to real biological systems, we need to develop more biologically plausible models of pattern formation in various tissues, such as the embryonic epithelium of Xenopus and Drosophila brain tumors, using empirically obtained biological Boolean networks. For example, one could investigate the epigenetic mechanisms of tumorigenesis depending on micro-environmental factors such as nutrient or oxygen gradients. One could build a genome-scale metabolic Boolean network model in a cancer cell, measure which regime the dynamics of that network would be in, and simulate cancer growth processes in a given micro environment. Such studies would allow for the analysis of a potential relationship between the properties of metabolic Boolean networks at the single-cell level and morphologies at a collective level, using biologically more realistic model assumptions.

For future work, we plan to track the growing processes from the seed cell to the cell aggregation to fully account for the spatial and temporal distribution of cells in terms of the criticality of GRNs in morphogenesis. Also, in order to solidify our findings on the nontrivial patterns, we will measure the complexity of our morphogenetic patterns using other measures such as information complexity. Furthermore, we will introduce somatic mosaic mutations in post-zygotic cells that are not inherited genetic alterations in the course of cell division, which are known to be present in actual embryo development [17, 20, 39]. Adding this kind of mutation to our model may reveal different roles of the criticality of GRNs under diverse forms of mutations.

This material is based upon work supported by the US National Science Foundation under Grant No. 1319152.