Functional Resilience of Mutually Repressing Motifs Embedded in Larger Networks

Elucidating the design principles of regulatory networks driving cellular decision-making has important implications for understanding cell differentiation and guiding the design of synthetic circuits. Mutually repressing feedback loops between ‘master regulators’ of cell fates can exhibit multistable dynamics enabling “single-positive” phenotypes: (high A, low B) and (low A, high B) for a toggle switch, and (high A, low B, low C), (low A, high B, low C) and (low A, low B, high C) for a toggle triad. However, the dynamics of these two motifs have been interrogated in isolation in silico, but in vitro and in vivo, they often operate while embedded in larger regulatory networks. Here, we embed these motifs in complex larger networks of varying sizes and connectivity to identify hallmarks under which these motifs maintain their canonical dynamical behavior. We show that an increased number of incoming edges onto a motif leads to a decay in their canonical stand-alone behaviors. We also show that this decay can be exacerbated by adding self-inhibition but not self-activation loops on the ‘master regulators’. These observations offer insights into the design principles of biological networks containing these motifs and can help devise optimal strategies for the integration of these motifs into larger synthetic networks.


Introduction
Gene Regulatory Networks (GRNs) are an integral part of the control structure involved in various cellular processes, such as cell-fate decisions made during embryonic development, cellular reprogramming, and phenotypic switching among two or more cell types. A pluripotent cell is capable of differentiating into more than one cell type in response to varying stimuli. This property of coexistence of more than one stable steady state (phenotypes) is referred to as multistability, and it underlies the dynamics of many GRNs involved in decision-making during differentiation [1]. Such multistability has been seen during cellular reprogramming as well as phenotypic switching under many circumstances. Thus, elucidating the dynamical principles of multistable GRNs and network motifs holds promise for understanding many biological processes and control applications in synthetic biology [2][3][4][5].
One of the most frequently observed and extensively investigated network motifs is the 'Toggle Switch' (TS), i.e., two mutually repressing regulators, A and B, each driving a different cell fate [6][7][8]. The TS enables two mutually exclusive "single-positive" outcomes-(high A, low B) and (low A, high B), thus showing bistable dynamics and allowing a pluripotent cell to choose from two cell fates [2,8,9]. For instance, PU.1 and GATA1 form a TS that drives hematopoietic stem cells to either a common myeloid progenitor (PU.1 high, GATA1 low) or an erythroid one (PU.1 low, GATA1 high) [2,10]. Furthermore,

Random Circuit Perturbation (RACIPE)
Random Circuit Perturbation (RACIPE) formalism generates a system of ordinary differential equations (ODEs) for a given network topology and simulates the ODEs by pooling parameters from a randomized predetermined range to identify dynamical properties of a network topology [23].
For a node T in the network having P i activating and N i inhibiting nodes with incoming edges, the ODE generated by RACIPE to represent the node T will be given by: where the terms T, P i , and N j are the concentrations of the nodes at time t, n is the Hill coefficient showing the influence of P i or N j on T, λ is the fold change in expression caused by node P i or N j upon acting on node T, P i 0 or N j 0 are the threshold values of Hill function, G T is the production rate, and k T is the degradation rate of the node T.
H S represents the shifted Hill equation and is defined by: For a particular topology file, RACIPE generates multiple randomized parameter sets and simulates them over multiple initial conditions to identify the steady-state levels of the nodes. The parameters are randomized by sampling from their respective pre-defined ranges (Table 1) given below: Simulations for all the networks were performed in triplicates, with 10,000 parameter sets per replicate and 100 initial conditions for each parameter set.

Normalizing the Steady State Values
The following formula first normalized the steady state values obtained from RACIPE simulations for each network: Here, S i N is the normalized steady-state value of the ith node; S i is the steady state value of ith node given in RACIPE output; G i is the production rate, and k i is the degradation rate parameter value of the ith node. The normalized values were then converted to Zscores using the zscore function of the SciPy library of python3. If a node showed a Z-score above zero, it was considered to be showing higher expression and was considered to be in "ON" state represented by 1, and if it showed a Z-score below zero, it was inferred to be in "OFF" state represented by 0. These criteria were then used to convert the steady state values into a string of zeros and ones representing the binarized steady state shown for a particular parameter set.

Sarle's Bimodality Coefficient (BiC)
Sarle's bimodality coefficient (BiC) [24] was used to identify the nature of the distribution of nodes' steady state values as either bimodal or not. The formula used to calculate bimodality is given by: where m 3 is the skew of the distribution, m 4 is the excess kurtosis, and the sample size is denoted by n. BiC varies from 0 to 1, with values above 0.55 (i.e., 5/9) representing bimodality in the distribution and those below this threshold indicating a unimodal distribution [24]. Skew and excess kurtosis were calculated using functions of the SciPy library of python3. The values were then plugged into the above equation to obtain the bimodality coefficients.

Correlation Coefficient (CC)
Correlation coefficients (CC) were calculated using the function spearmanr of the SciPy library of python3.

Stand-Alone Properties of Toggle Switch and Toggle Triad
We first investigated the stand-alone properties of a toggle TS and TT before embedding them in larger networks. RACIPE formalism [23] (see Section 2) was used to simulate these networks for 10,000 randomized parameter sets sampled from a predetermined parameter space. Three such replicates (each with 10,000 parameter sets × 100 initial conditions per parameter set) were performed for each motif. The resultant steady-state values for each node were normalized and converted to Z-scores. Pairwise correlation coefficients (CC) were then calculated between the steady state values of nodes in a TS or TT. Furthermore, for the TS, we calculated Sarle's bimodality coefficient (BiC) [24] values for each node. BiC values range from 0 to 1, with values closer to 1 representing higher bimodality, and any value above 0.55 is considered to represent a bimodal distribution [24].
The TS motif ( Figure 1A) mainly showed two "single-positive" steady states: ((low A, high B); (A, B) = (0, 1)) and ((high A, low B); (A, B) = (1, 0)), as observed in the bivariate plot ( Figure 1B). The steady-state values of both nodes in a TS showed a bimodal distribution (BiC A = BiC B = 0.78), with the peaks representing the corresponding high and low steady-state values. Because the two nodes of a TS repress each other, the correlation coefficient between the steady state values of the nodes was strongly negative (CC AB = −0.83) ( Figure 1C). The steady-state values for a TT motif ( Figure 1D), for any given pair of TT nodes, had three distinct clusters, which represent the three "single-positive" stable, steady states, namely ((low A, low B, high C); (A, B, C) = (0, 0, 1)) state, ((low A, high B, low C); (A, B, C) = (0, 1, 0)) state and ((high A, low B, low C); (A, B, C) = (1, 0, 0)) state. One node shows higher expression in these states, while the other two nodes have repressed expression ( Figures 1E and S4A,B). In contrast to TS, the average BiC of a node in TT is 0.43 (standard deviation = 0.004), implying a more unimodal-like distribution with little difference in the high and low steady-state values of a node in a TT. Although negative, the magnitude of the pairwise correlation coefficient between the steady state values of any two nodes of a TT was less than that of a TS (CC AB = CC AC = −0.39, CC BC = −0.40) ( Figures 1F and S4C,D). This decrease could be because although any pair of two nodes mutually repress, due to the dominance of "single-positive" steady states, two nodes of  Figures 1E and S4A,B). This trend leads to a relatively lower magnitude of the correlation coefficient between any two nodes of a TT motif. On the other hand, two nodes of a TS are strictly confined to having opposing expression profiles, leading to a strongly negative correlation. expression ( Figures 1E and S4A,B). In contrast to TS, the average BiC of a node in TT is 0.43 (standard deviation = 0.004), implying a more unimodal-like distribution with little difference in the high and low steady-state values of a node in a TT. Although negative, the magnitude of the pairwise correlation coefficient between the steady state values of any two nodes of a TT was less than that of a TS (CC AB = CC AC = −0.39, CC BC = −0.40) ( Figures 1F and S4C,D). This decrease could be because although any pair of two nodes mutually repress, due to the dominance of "single-positive" steady states, two nodes of TT can still show similar low-expression steady-state values ( Figures 1E and S4A,B). This trend leads to a relatively lower magnitude of the correlation coefficient between any two nodes of a TT motif. On the other hand, two nodes of a TS are strictly confined to having opposing expression profiles, leading to a strongly negative correlation.

Functional Traits of Toggle Switch Depend on Density Rather Than the Size of the Larger Networks It Is Embedded in
Next, we embedded TS and TT motifs in different larger networks having combinations of four different network orders and three distinct densities (mean connectivity) to understand how the abovementioned stand-alone dynamic traits of TS and TT change. The four network orders are 5N, 10N, 15N, and 20N, where N is the number of nodes in a network in which these motifs were embedded. Mean connectivity, defined as the ratio of the number of edges to the number of nodes in a network, was used to represent the density of the network. The three mean connectivity values are E:2N, E:4N, and E:6N, where E:xN signifies that the number of edges (E) is x times the number of nodes (N). The combinations of the four network orders (5N, 10N, 15N, and 20N) and three mean connectivity (E:2N, E:4N, and E:6N) resulted in twelve different types of networks ( Figure  2A). For each type of network, n = 100 random network topologies were generated. TS

Functional Traits of Toggle Switch Depend on Density Rather Than the Size of the Larger Networks It Is Embedded in
Next, we embedded TS and TT motifs in different larger networks having combinations of four different network orders and three distinct densities (mean connectivity) to understand how the abovementioned stand-alone dynamic traits of TS and TT change. The four network orders are 5N, 10N, 15N, and 20N, where N is the number of nodes in a network in which these motifs were embedded. Mean connectivity, defined as the ratio of the number of edges to the number of nodes in a network, was used to represent the density of the network. The three mean connectivity values are E:2N, E:4N, and E:6N, where E:xN signifies that the number of edges (E) is x times the number of nodes (N). The combinations of the four network orders (5N, 10N, 15N, and 20N) and three mean connectivity (E:2N, E:4N, and E:6N) resulted in twelve different types of networks ( Figure 2A). For each type of network, n = 100 random network topologies were generated. TS motifs were then embedded into these 1200 (12 types × 100 random networks) randomly generated networks to study the motif's behavior. For instance, a TS embedded in a 5N, E:4N network will have seven nodes (5 nodes + 2 nodes of TS) and 22 edges (5 × 4 edges + 2 edges of TS) ( Figure 2B). The same process was repeated for the TT motif to generate 1200 random networks. RACIPE was used to simulate the dynamics of these larger networks of varying sizes and mean connectivity, with each network being simulated three times ( Figure 2C).
The generated outputs were normalized by Z-scoring and then analyzed to characterize the behavior of TS and TT motifs upon embedding. Three metrics were assessed to quantify the dynamic resilience of TS and TT motifs when embedded in larger networks: bimodality coefficient (BiC), correlation coefficient (CC), and frequency of canonical 'single-positive' states as a fraction of all steady states observed (F1). motifs were then embedded into these 1200 (12 types × 100 random networks) randomly generated networks to study the motif's behavior. For instance, a TS embedded in a 5N, E:4N network will have seven nodes (5 nodes + 2 nodes of TS) and 22 edges (5 × 4 edges + 2 edges of TS) ( Figure 2B). The same process was repeated for the TT motif to generate 1200 random networks. RACIPE was used to simulate the dynamics of these larger networks of varying sizes and mean connectivity, with each network being simulated three times ( Figure 2C). The generated outputs were normalized by Z-scoring and then analyzed to characterize the behavior of TS and TT motifs upon embedding. Three metrics were assessed to quantify the dynamic resilience of TS and TT motifs when embedded in larger networks: bimodality coefficient (BiC), correlation coefficient (CC), and frequency of canonical 'single-positive' states as a fraction of all steady states observed (F1). To evaluate how network size and density can independently influence the behavior of a TS when embedded in larger networks of varying orders and mean connectivity, we compared the behavior for networks sharing the same mean connectivity but having different network orders or vice versa. Interestingly, for networks with the same mean connectivity, CC between the two nodes in a TS (CC AB) did not show any significant variation for varied network orders ( Figure 3A(i)). However, when controlling for network order, CC AB reduced as the mean connectivity increased ( Figure 3A(ii)). A similar trend, i.e., the dependence on network density rather than on network size and a decrease in magnitude with increasing network density, was also observed in distributions of BiC values: BiC A and BiC B (Figures 3B(i,ii) and S5A(i,ii)) and for F1 ( Figure 3C(i,ii)).
For some cases, the decrease observed in metrics between E:4N and E:6N mean connectivity values was not significantly different, potentially because the TS dynamics were compromised enough in the E:4N case but not in the E:2N case. As a control case, we plotted the distributions of the CC and F1 value for any two randomly chosen nodes which did not belong to the TS motif to compare them with those of TS nodes for networks of order 5N and mean connectivity E:2N. The distributions for the randomly chosen nodes To evaluate how network size and density can independently influence the behavior of a TS when embedded in larger networks of varying orders and mean connectivity, we compared the behavior for networks sharing the same mean connectivity but having different network orders or vice versa. Interestingly, for networks with the same mean connectivity, CC between the two nodes in a TS (CC AB) did not show any significant variation for varied network orders ( Figure 3A(i)). However, when controlling for network order, CC AB reduced as the mean connectivity increased ( Figure 3A(ii)). A similar trend, i.e., the dependence on network density rather than on network size and a decrease in magnitude with increasing network density, was also observed in distributions of BiC values: BiC A and BiC B (Figures 3B(i,ii) and S5A(i,ii)) and for F1 ( Figure 3C(i,ii)).
For some cases, the decrease observed in metrics between E:4N and E:6N mean connectivity values was not significantly different, potentially because the TS dynamics were compromised enough in the E:4N case but not in the E:2N case. As a control case, we plotted the distributions of the CC and F1 value for any two randomly chosen nodes which did not belong to the TS motif to compare them with those of TS nodes for networks of order 5N and mean connectivity E:2N. The distributions for the randomly chosen nodes differed significantly when compared to those of the TS nodes ( Figure S5B(i,ii)), showing that the distributions of metrics obtained for TS nodes when embedded were indeed unique to the TS. Put together, all the three metrics considered here to capture the behavior of a TS-CC (how strongly the two nodes in a TS are anti-correlated), BiC (how clearly the high and low levels of a node are segregated), and F1 (how strong is the dominance of mutually exclusivity of the two nodes)-tend to show trends indicating a weakening of the dynamical behavior of a TS, as it is embedded in increasingly denser large networks. In other words, the gradual decrease in a switch-like behavior, as noted through these metrics, is seen in the case of TS being embedded in increasing network density (i.e., mean connectivity) but not in the scenario of increasing network size (order).
Biomolecules 2022, 12, x FOR PEER REVIEW 7 of 21 differed significantly when compared to those of the TS nodes ( Figure S5B(i,ii)), showing that the distributions of metrics obtained for TS nodes when embedded were indeed unique to the TS. Put together, all the three metrics considered here to capture the behavior of a TS-CC (how strongly the two nodes in a TS are anti-correlated), BiC (how clearly the high and low levels of a node are segregated), and F1 (how strong is the dominance of mutually exclusivity of the two nodes)-tend to show trends indicating a weakening of the dynamical behavior of a TS, as it is embedded in increasingly denser large networks. In other words, the gradual decrease in a switch-like behavior, as noted through these metrics, is seen in the case of TS being embedded in increasing network density (i.e., mean connectivity) but not in the scenario of increasing network size (order).

Local Density around a Toggle Switch Impacts its Dynamic Behavior
The mean connectivity of a network is the ratio of the total number of edges to the number of nodes in the network, i.e., a measure of global network density. Thus, assuming that the network, on average, is equally sparse or dense, with an increase in mean connectivity of the network, the average in-degree of the nodes of a TS embedded in the network also increases. To ascertain whether this increase in the in-degree (in A and in B) for the TS nodes (as a consequence of the increased mean connectivity of the network) contributed to divergence from stand-alone TS dynamics, we analyzed the variation in the three metrics (CC, BiC, and F1) with a change in the in-degrees of the TS nodes. We observed that as the in-degree of both nodes of a TS increased, the mean CC AB values decreased

Local Density around a Toggle Switch Impacts Its Dynamic Behavior
The mean connectivity of a network is the ratio of the total number of edges to the number of nodes in the network, i.e., a measure of global network density. Thus, assuming that the network, on average, is equally sparse or dense, with an increase in mean connectivity of the network, the average in-degree of the nodes of a TS embedded in the network also increases. To ascertain whether this increase in the in-degree (in A and in B) for the TS nodes (as a consequence of the increased mean connectivity of the network) contributed to divergence from stand-alone TS dynamics, we analyzed the variation in the three metrics (CC, BiC, and F1) with a change in the in-degrees of the TS nodes. We observed that as the in-degree of both nodes of a TS increased, the mean CC AB values decreased in magnitude, i.e., the TS nodes A and B were not as strongly negatively correlated with one another as in a stand-alone case ( Figure 4A(i)). For an in-degree of one for both the nodes in a TS, i.e., the case when the two nodes only had outgoing edges apart from their mutual inhibitions, we noticed a mean CC AB value of −0.83 ( Figure 4A(i)), the same as that of an isolated TS motif ( Figure 1B). This trend is observed because the TS in this scenario acts like an independent motif as it is not being affected by the network due lack of incoming edges. Furthermore, the magnitude of CC AB showed the fastest decline when both the nodes had equally increased in degrees (along the diagonal of the heatmap shown in Figure 4A(i)). Similarly, F1 decreased steadily with in-degrees increasing equally for the two nodes ( Figure 4A(ii)). On the other hand, the BiC of a given node in TS changed only with the in-degree for that node and not with the in-degree for the other node or with the overall in-degree of a TS ( Figure S6A(i,ii)). Additionally, the mean BiC values for nodes with an in-degree more than two were lower than the typical cut-offs considered for bimodality (~0.55), indicating the compromised canonical bimodal distributions observed in the nodes of an isolated TS.
Further, we quantified the impact of asymmetry in terms of incoming edges on a TS by considering the impact of the ratio of in-degree of A to that of B (log 2 (in A/in B)) on CC AB and relative BiC values simultaneously. We noticed that the higher the asymmetry in terms of in-degree (log 2 (in A/in B) > 1 or log 2 (in A/in B) < −1), the stronger the negative correlation between the two TS nodes (CC AB < −0.5) ( Figure 4A(iii)). Importantly, as the magnitude of log 2 (in A/in B) increased, the range of CC AB values (initially even spanning positive values; above the red dotted horizontal line) narrowed to highly negative values close to those observed in isolated TS nodes (−0.83) ( Figure 4A(iii)). Similarly, the value of F1 approached closer to that observed in an isolated TS, as the magnitude of (log 2 (in A/in B)) increased ( Figure S6B(i)). Furthermore, we noted that the higher the in-degree of a node, the more likely it becomes for that node to lose its bimodality seen in a stand-alone TS (log 2 (BiC A/BiC B) < −1 for log 2 (in A /in B) > 1 and log 2 (BiC A/BiC B) > 1 for log 2 (in A/in B) < −1) ( Figure 4A(iii)). Together, this analysis suggests that while F1 (fraction of single positive states) and CC (AB) (correlation coefficient) depend on the in-degree of a TS, the BiC of individual nodes depends on the in-degree of the respective node.
To substantiate this trend further, we investigated representative cases of varied indegrees of A and B. When the TS has only outgoing connections, and there is no asymmetry between the in-degree of A and B (in A = in B= 1; the only incoming links on A and B are from each other), the bivariate plot of A and B is very similar to that of an isolated TS (compare Figure 1B with Figure 4B(i)). However, upon asymmetry in the in-degrees of nodes in a TS (in A = 2, in B = 6, in TS = 2 + 6 = 8), the node with higher in-degree (B) starts to lose its switch-like behavior and shows a more unimodal distribution of its steady state values ( Figure 4B(ii)). However, the strong negative correlation between the TS nodes and, concomitantly, the fraction of "single-positive" states does not decrease as sharply compared to those for the stand-alone case (compare F1 and CC (AB) in Figure 4B(ii) with those in Figure 4B(i)). To deconvolute the impact of a higher in-degree of TS vs. asymmetry in the in-degree of both the nodes, we considered a case with the same in-degree for both nodes without changing the net in-degree for a TS (in A = 4, in B = 4, in TS = 4 + 4 =8).
Here, the switch-like behavior of both nodes is largely lost, and they show a unimodal distribution of their respective steady-state values ( Figure 4B(iii)). However, CC (AB) and F1 are comparable to the case of asymmetric in-degrees (compare CC (AB) and F1 values in Figure 4B(ii) with those in Figure 4B(iii), further supporting that while the higher the in-degree of a toggle switch, the weaker the negative correlation between nodes and the smaller the fraction of "single-positive" states, the bimodality patterns for each node depend on in-degree of that individual node and not on in-degree for TS.
After looking at these representative trends showcasing an increasing in-degree of a node leading to loss of bimodality in steady-state distributions of the corresponding node, we investigated how generic these trends were for embedded TS motifs. We hypothesized that with an increasing in-degree of TS, the frequency of co-occurring 'single positive' states (01 and 10 states in a bistable setting) should decrease, with a concomitant increase in having either of these two states, i.e., 01 or 10 states in a monostable setup. This feature can be quantified by the ratio of the fraction of bistable parameter sets showing 01 and 10 steady states to the fraction of monostable parameter sets showing 01 or 10 steady states (B:M). When the TS has only outgoing connections to the network (in A = in B = 1; the only incoming links on A and B are from each other), the B:M ratio is greater than one ( Figure 4C). However, as the in-degree for the TS increases, B:M decreases to values below one, with the sharpest decline when both nodes have equal in-degrees (along the diagonal of the heatmap in Figure 4C). These results show that as the in-degree of TS increases, the canonical bistable behavior (co-existing 'single positive' states) starts to decrease, and simultaneously, the fraction of 'single-positive' steady states increases, implying a loss of bistable traits of the TS motif.
one another as in a stand-alone case ( Figure 4A(i)). For an in-degree of one for both the nodes in a TS, i.e., the case when the two nodes only had outgoing edges apart from their mutual inhibitions, we noticed a mean CC AB value of −0.83 ( Figure 4A(i)), the same as that of an isolated TS motif ( Figure 1B). This trend is observed because the TS in this scenario acts like an independent motif as it is not being affected by the network due lack of incoming edges. Furthermore, the magnitude of CC AB showed the fastest decline when both the nodes had equally increased in degrees (along the diagonal of the heatmap shown in Figure 4A(i)). Similarly, F1 decreased steadily with in-degrees increasing equally for the two nodes ( Figure 4A(ii)). On the other hand, the BiC of a given node in TS changed only with the in-degree for that node and not with the in-degree for the other node or with the overall in-degree of a TS ( Figure S6A(i,ii)). Additionally, the mean BiC values for nodes with an in-degree more than two were lower than the typical cut-offs considered for bimodality (~0.55), indicating the compromised canonical bimodal distributions observed in the nodes of an isolated TS.  Finally, across larger networks of varying sizes and mean connectivity values, we interrogated how in-degrees for individual nodes as well as for a TS (In A, In B, and In TS) correlates with various metrics-F1, CC AB, BiC A, BiC B, and B:M ratio. Net in-degree of the TS (in TS) was found to best explain the decline in the magnitude of CC AB, F1, and the B:M ratio ( Figure 4D). The bimodality coefficients (BiC A and BiC B), on the other hand, were more influenced by the in-degree of their respective nodes and not in TS ( Figure 4D). Therefore, it is the local density on the TS motif (given by in TS) that drives the divergence from TS-like behavior rather than the properties of the whole network in which a TS is embedded.

In a Toggle Triad, the Fraction of Single-Positive States Captures Its Functional Resilience
After investigating the patterns seen in a TS embedded in large networks, we focused our attempt on understanding the functional resilience of the TT motif. We embedded it into the previously described 12 types of large random networks. Similar to observations in TS, the distributions of pairwise correlation coefficients-CC AB, CC BC, and CC AC-did not show any significant, consistent variation when they were grouped by mean connectivity and compared across the different network orders (Figures 5A(i) and S7A(i),B(i)). Intriguingly, unlike the observations in TS, we did not observe any significant differences when the CC values between the TT nodes were grouped by their network orders and compared across the three mean connectivity either (Figures 5A(ii) and S7A(ii),7B(ii)), despite a visible increase in the range of values. Thus, we investigated how the maximum of the three pairwise correlation values (MaxCC) between the TT nodes-CC AB, CC BC, and CC AC-varied as a function of network order and/or mean connectivity. We observed that when grouped by order, the higher the mean connectivity, the higher the average MaxCC value; however, no such trend was seen when grouped by mean connectivity (Figure 5B(i,ii)). Reminiscent of observations in TS, the fraction of "single-positive" (010, 100, and 001) steady states (F1) also decreased overall when comparisons were made across their mean connectivity ( Figure 5C(i)), but not across their network orders ( Figure 5C(ii)). Consistently, the fraction of "double-positive" (011, 110, 101) states (F2) and the fraction of "all-positive" or "all-negative" (111, 000) states (F3) increased across different mean connectivity values when grouped by network orders (Figure 5D(i,ii)), but not when grouped by mean connectivity and compared across network orders. (Figure S8A(i,ii)). The ratio of the fraction of "single-high" to that of "double-high" states (F1/F2) also showed the same trend, asymptotically reaching the value of one (Figures 5D(iii) and S8(iii)). These observations help understand the trend seen for MaxCC tending towards positive values with increasing mean connectivity. With a decreasing frequency of "single-positive" states, the negative pairwise correlation between the TT nodes starts to weaken, and in some cases, one or more of the correlation coefficient values can be positive, suggesting a decay of stand-alone TT dynamics.
To understand if the distributions of MaxCC and F1/F2 values were unique to the TT motif, the distributions of these metrics for TT embedded in 5N and E:2N networks were compared with the distributions of metric values of three randomly chosen nodes which were not a part of the TT motif. We observed that the distributions of both MaxCC ( Figure S9A(i)) and F1/F2 ( Figure S9A(ii)) were significantly different for the randomized nodes when compared to those of TT nodes showing that the distributions of these metrics are indeed unique to TT motif nodes.
We next investigated how the in-degree of the embedded TT motif affected its behavior. MaxCC values correlated positively (ρ = 0.41, p < 0.05) with the in-degree of TT, indicating that the higher the in-degree of TT, the stronger the decay of TT dynamics ( Figure 6A(i)). Consistently, F1/F2 values decreased as the in-degree of TT increased (ρ = −0.51, p < 0.05), tending towards a value of one for high in-degrees of TT ( Figure 6A(ii)), driven by the decrease in F1 and increase in F2 and F3 ( Figure S10A). Moreover, Max CC correlated negatively with F1 (ρ = −0.48, p < 0.05) and F1/F2 (ρ = −0.45, p < 0.05) but positively with F2 (ρ = 0.37, p < 0.05) and F3 (ρ = 0.35, p < 0.05) (Figures 6A(iii) and S10B). Therefore, with an increasing in-degree of a TT, the fraction of "single-positive" states diminishes as they are replaced by "double-positive" (and, to some extent, by "all-positive" or "all-negative") states. Subsequently, this change in frequencies of different states can weaken the canonical mutual inhibition among nodes in a stand-alone TT, driving one or more pairwise correlation values to be positive, thus validating our choice of Max CC as a metric to assess the decay of TT dynamics. To understand if the distributions of MaxCC and F1/F2 values were unique to the TT motif, the distributions of these metrics for TT embedded in 5N and E:2N networks were compared with the distributions of metric values of three randomly chosen nodes which were not a part of the TT motif. We observed that the distributions of both MaxCC ( Figure  S9A(i)) and F1/F2 ( Figure S9A(ii)) were significantly different for the randomized nodes when compared to those of TT nodes showing that the distributions of these metrics are with F2 (ρ = 0.37, p < 0.05) and F3 (ρ = 0.35, p < 0.05) (Figures 6A(iii) and S10B). Therefore, with an increasing in-degree of a TT, the fraction of "single-positive" states diminishes as they are replaced by "double-positive" (and, to some extent, by "all-positive" or "all-negative") states. Subsequently, this change in frequencies of different states can weaken the canonical mutual inhibition among nodes in a stand-alone TT, driving one or more pairwise correlation values to be positive, thus validating our choice of Max CC as a metric to assess the decay of TT dynamics. After characterizing the effect of the net in-degree of TT on its behavior, we investigated how the three pairwise correlations between the TT nodes (CC AB, CC BC, and CC AC) varied with varying in-degrees for the nodes. Unlike the observations for embedded TS, where the CC AB decreased with network mean connectivity and in-degree for motif nodes (Figures 4A(i) and S11), the pairwise CCs did not show any discernible trend in an embedded TT; therefore, we excluded them from any further analysis. On the other hand, F1, F2, and F3 changed significantly as the in-degrees of any two nodes in the TT increased (F1 decreased while F2 and F3 increased); at high in-degrees for any two nodes, F2 is approximately equal to F1 and six to seven times the corresponding F3 values ( Figure S12). Next, we quantified changes in these metrics brought about by simultaneously varying the in-degree of the third node as well. An increase in C, while maintaining the values of in A and in B, led to a lower F1, higher F2, and F3 values, as expected (compare corresponding cells in Figure S13B with Figure S13A), thereby showcasing the impact of increasing net in-degree on the breakdown of stand-alone TT dynamics.
Because pairwise correlation coefficients between TT nodes were unable to gauge the changes in TT dynamics, we performed multiple linear regression (MLR) on steady-state values of the TT nodes to understand the changes in inter-node dependence brought about by embedding the TT. MLR was performed by taking steady-state values of A and B as independent variables and steady states of C as the dependent variable. The distributions of coefficients of A and B (ACoeff and BCoeff, respectively, i.e., the slopes of the regression plane) and the model's scores (R 2 values) were then compared across network orders by grouping them according to common mean connectivity and vice versa. As the mean connectivity of the network increased, the magnitude of the mean values of ACoeff and BCoeff decreased, indicating that the repressing effect of the two nodes on C decreased as the network density increased ( Figure S14A,B(i,ii)). Furthermore, the scores of MLR models declined as the mean connectivity of the network increased, pointing towards a decline in the predictive power of the MLR model as the TT dynamics deviated from their canonical behavior with the increasing network density ( Figure S14A(iii),B(iii)). Further, the higher the in-degree of TT, the lower the model's score and the magnitudes of ACoeff and BCoeff values ( Figure S14C), reiterating the disruptive impact of the in-degree of TT on its canonical behavior. Intriguingly, for E:6N, we noticed some ACoeff and BCoeff values to be positive ( Figure S14C), reminiscent of MaxCC values also becoming greater than zero for high mean connectivity ( Figure 6A(iii)).
Next, we performed a meta-analysis across network orders and mean connectivity in terms of correlation of various metrics with the in-degree of individual nodes as well as the net in-degree of TT. With all metrics considered here (F1, F2, F3, Max CC, F1/F2), In TT showed a stronger correlation (ρ > 0 for F2, F3; ρ < 0 for F1, F2/F2, MaxCC) than the in-degree of any of the individual nodes (in A, in B, in C), indicating that net in-degree was the best predictor of embedded TT dynamics. (Figure 6B), and reinforcing the trends seen in embedded TS ( Figure 4D). This trend was consistent even when the individual in-degree or combined in-degree of two nodes, normalized by net in-degree (in A/in TT, in B/in TT, in C/in TT, in AB/in TT, in BC/in TT, in AC/in TT) was considered ( Figure 6C). Intriguingly, while TT correlated strongly with MaxCC (i.e., the pair of nodes whose "mutual exclusion" is the weakest), it did not correlate as strongly with individual pairwise correlation coefficients (CC AB, CC BC, CC AC) ( Figure 6C). Furthermore, while MaxCC correlated strongly with individual in-degree as well as with net in-degree, Min CC (the minimum value among three pairwise correlation coefficients, i.e., the pair of nodes showing the strongest mutual inhibition properties) did not show any such association ( Figure 6B,C). This difference in terms of MaxCC vs. MinCC is consistent with observations that all three individual pairwise correlation coefficients show only a weak correlation with both in-degree of TT and with F1 ( Figure 6D), thus justifying our choice of using MaxCC as a metric to track the dynamics of embedded TT.
Put together, these results highlight that while the total in-degree of a TT motif is a good indicator of the decay of TT dynamics (decrease in "single-positive" states and simultaneous increase in "double-positive"), it often does not contain precise information on which out of the three possible pairs of nodes in a TT have their mutual repression compromised and thus drive a decay in the stand-alone dynamical behavior of a TT. Despite this limitation, similar to the results seen for TS, the dynamics of embedded TT were determined not by network order, network density, or individual in-degrees of the motif nodes but by the total in-degree of the motif (In TT) ( Figure 6E). Moreover, both for the TS and TT motifs, the nature of these in-degrees-being either activating or repressingdid not affect the nature of divergence from canonical behavior ( Figure S15), highlighting that the motif is sensitive to the total number of incoming edges, not their distributions in terms of their sign/effect.

Effect of Self-Activation and Self-Inhibition of Nodes on the Modularity of Motifs
Many "master regulators" involved in driving cell fate decisions form TS or TT motifs, such as PU.1 and GATA1 in hematopoietic stem cells, and CDX2 and OCT4 in embryonic stem cells can also self-activate to stabilize an "undecided" multipotent state characterized by promiscuous expression of motif nodes [2,13]. In order to evaluate the effects of such self-activations and self-inhibitions on the motif nodes, the TS and TT motifs containing self-activation (TS-SA and TT-SA) and self-inhibition edges (TS-SI and TT-SI) on all nodes were embedded into combinations of two network orders (5N and 20N) and two mean connectivity (E :2N and E:6N), resulting in four types of networks. The same pipeline as before was followed to generate the networks, with n = 100 for each type of network (i.e., 400 networks per motif). The motifs TS-SA, TS-SI, TT-SA, and TT-SI were then embedded into these networks and simulated with RACIPE, and corresponding metrics were compared. First, we compared the TS, TS-SA, and TS-SI embedded in networks having the same mean connectivity. For all three metrics (CC AB, F1, and BiC), TS, TS-SA, and TS-SI showed significantly different distributions ( Figure 7A(i-iii)). While including self-activation moderately strengthened the correlation between A and B, self-inhibition significantly weakened it ( Figure 7A(i)). Similarly, while including self-activation led to slightly increased F1, and BiC values, conversely, including self-inhibition on TS nodes reduced the F1 and BiC values ( Figure 7A(ii,iii)). These results suggest that while adding self-activation on nodes can preserve the dynamical features of a TS embedded in large networks, self-inhibition can accelerate the decay of TS dynamics, possibly offering a reason for an observed higher frequency of self-activating 'master regulators' rather than the self-inhibiting ones [25].
Next, we probed the impact of self-activation and self-inhibition in the case of an embedded TT. Similar to observations in TS, self-inhibition had an opposite and a stronger impact in influencing the metrics as compared to self-activation. Median F1/F2 values noted for TT-SA were higher than those for both TT and TT-SI, irrespective of the mean connectivity ( Figure 7B(i)). Consistently, median F3 values are higher for TT-SI than for TT and TT-SA ( Figure 7B(ii)). Together, these results indicate that while adding self-activation on nodes of a TT can enrich for canonical "single-positive states", adding self-inhibition can enrich for "all-positive" or "all-negative" states instead. Similarly, for Max CC, TT-SI had higher median values than both TT and TT-SA cases ( Figure 7B(iii)), resulting in a faster decay of canonical properties of TT. In other words, the self-inhibiting "master regulators" can exhibit compromised phenotypic decision-making as compared to non-selfinhibiting ones.
The effect of the addition of self-activation and self-inhibition edges on the nodes of the motifs was also gauged by embedding the motif, its self-activation, and self-inhibition version in the same ensemble of random networks having network orders 5N and mean connectivity of E:2N and E:6N ( Figure S16A). These networks were then simulated using RACIPE [23]. The metrics of both the motifs showed significantly different distributions when grouped by common network mean connectivity ( Figure S16B(i-iii),C(i,ii)). A similar trend of the addition of self-activation leading to preserving and the addition of self-inhibition leading to a decay in the dynamical properties of motifs was observed. This trend was also reflected in the ratios of the metric values of TT-SA or TT-SI motifs to the corresponding metric values of TT motifs embedded in the same larger network ( Figure S16C(iii)), indicating that the addition of self-activation led to a decay of the dynamical properties even when the larger networks in which the motifs were embedded were maintained to be the same. Next, we probed the impact of self-activation and self-inhibition in the case of an embedded TT. Similar to observations in TS, self-inhibition had an opposite and a stronger impact in influencing the metrics as compared to self-activation. Median F1/F2 values noted for TT-SA were higher than those for both TT and TT-SI, irrespective of the mean connectivity ( Figure 7B(i)). Consistently, median F3 values are higher for TT-SI than for TT and TT-SA ( Figure 7B(ii)). Together, these results indicate that while adding self-activation on nodes of a TT can enrich for canonical "single-positive states", adding self-inhibition can enrich for "all-positive" or "all-negative" states instead. Similarly, for Max CC, TT-SI had higher median values than both TT and TT-SA cases ( Figure 7B(iii)), resulting in a faster decay of canonical properties of TT. In other words, the self-inhibiting "master regulators" can exhibit compromised phenotypic decision-making as compared to nonself-inhibiting ones.
The effect of the addition of self-activation and self-inhibition edges on the nodes of the motifs was also gauged by embedding the motif, its self-activation, and self-inhibition Finally, we assessed whether adding self-regulation preserves the correlation of the in-degree of the motif (TS or TT) with various metrics (BiC A, BiC B, B:M, F1 and CC AB for TS; F1, F2, F1/F2 and MaxCC for TT). We consistently observed that motifs with selfinhibition had a lower magnitude of correlation with all these motif properties as compared to motifs with no self-regulation or with self-activation ( Figure 7C(i,ii)), indicating a faster loss and subsequent saturation of the motif's dynamical properties. These trends reaffirm that the addition of self-activating loops preserves the canonical behavior of the motif while adding self-inhibition loops diminishes it remarkably ( Figure 7D).

Discussion
Investigating the modularity of biological networks has been an active area of research. Modularity has been loosely defined as corresponding to a highly interconnected set of nodes such that the density of connections within that module is significantly higher than that of the density of connections of this module with other modules [26]. Thus, modularity has been mostly studied from a network topology perspective rather than a functional one. Similarly, the concept of network motifs-recurring sets of regulatory interactions that appear more frequently than expected in a given network-also highlights network substructures based on their topology [27]. While the dynamics of such motifs have received extensive attention [27], how insular or intact the dynamics of these network motifs are when embedded in a large network remains largely underexplored.
Here, we investigated how 'modular' the behavior of a TS or TT is when embedded in large random networks of varying sizes and densities. For both these motifs, we observed that an increase in local density around them (i.e., the number of incoming edges on TS or TT) was capable of changing their dynamical behavior rather than any global topological properties associated with the large network. Although we witnessed how increasing mean connectivity of the network changed the distributions of the metrics we have used to characterize the dynamical behavior of TS or TT, further analysis revealed that it was the increasing in-degree of the motif (reflected, in part, by mean connectivity as well) that was driving this change in motif behavior. For the TS, it was found that all three metrics, BiC, CC, and the fraction of 'single-positive' (01, 10) states were suitable to gauge its change in dynamics. On the other hand, for a TT motif, CC metrics were not suitable to gauge its behavior, as they did not show extensive variation upon being embedded in large networks. This trend could be because CC, being a pairwise metric, is not optimal for capturing the variations in the steady state values of all three nodes. On the other hand, the maximum value of all three pairwise correlation coefficients between the TT nodes, MaxCC, was found to be suitable for gauging the change in the dynamics of TT. MaxCC values, being correlated positively with double-positive states (F2), were able to capture the enrichment of these states as the in-degree of TT increased. Similar to TS, the fraction of 'single-positive' states was also found to be a good metric for observing the variations in the dynamics of an embedded TT. We found that as the in-degree of an embedded TT increased, the single-positive states and double-positive states were almost comparable in frequency. This observation suggests that in the case of CD4+ T-helper cell differentiation-a case study of TT dynamics-the "double positive" cell states (hybrid Th1/Th2, Th1/Th17, and Th2/ Th17 phenotypes) seen experimentally [28,29] could exist due to the TT between GATA3, RORγT, and T-bet being driven by various other stimuli that impinge on these nodes via activation or repression. Therefore, besides self-activation [2,30], embedding in large networks can be an additional way to enrich such 'hybrid' states, as being increasingly reported in various biological systems [31][32][33], both in cases of TS and TT.
The bistable dynamics of a TS motif have been extensively investigated. The nodes of a TS motif show a bimodal distribution [34][35][36] of their steady states, which are dependent upon a balance between the kinetic parameters of the two nodes [14,37]. An asymmetry in these parameters bringing an imbalance in inhibitory strength can lead to a one-way cause-effect relation between the two TS nodes, with only one node showing bimodality while the other node shows unimodal behavior [14,20,34,37]. Moreover, an increase in the number of downstream interacting elements can also induce competition between produced proteins to bind to the promoter sites of the TS genes, changing the dynamics of motifs and potentially leading to a loss of bimodality [38]. Here, we show that the relative in-degrees of the nodes can also contribute to this behavior. A skew in the in-degrees led to a divergence from the bistable behavior of a TS, much like the skew in the kinetic parameters. Thus, our analysis uncovers an important design principle of gene regulatory networks (GRNs)-in order to maintain bistable features for a TS, the in-degree, which represents the number of regulators acting at a given point of time, should be minimal. This observation is reminiscent of previous studies demonstrating that in an E. coli transcriptional network, no transcription factor had an in-degree or out-degree greater than 2, and this feature played a key role in enabling the robustness of the network [39]. Thus, the in-degrees we report here for robustness are in good agreement with those seen in networks for GRNs of various organisms [39]. Similarly, another recent study associated the dynamical robustness of networks with their low in-degree [40]. Therefore, our results explain why network motifs such as TF, despite being embedded in much larger networks, can often exhibit their stand-alone behavior of multistability and "all-or-none" behavior, owing to their resilient dynamical features. Together, these dynamical traits can help us design optimal strategies to design and integrate synthetic circuits into the GRN of a cell rather than developing stand-alone/isolated modules.
Another aspect of many 'master regulator' TFs engaged in TS or TT is self-activation (positive auto-regulation). Our results here highlight that self-activation, but not selfinhibition, can help preserve the functional resilience of the network TFs. Thus, besides self-sustaining their levels to enable cell-state stability, self-activation in a TS or TT can aid in enabling robustness to fluctuations in a dynamic intra-and extra-cellular environment.
Notably, many 'master regulator' TFs are known to contain intrinsically disordered protein regions (IDPRs) or are intrinsically disordered proteins (IDPs) [41][42][43]. This disorder can govern the binding specificity of these TFs [44]; they can have time-varying interactions with other biomolecules and promiscuously interact with many of them, leading to possible rewiring of the protein interaction network [45], driven through the layer of critical nodes connecting the IDPs to the network at large [46]. Therefore, due to such 'conformational noise' scenarios, not each incoming edge to a TS or TT is necessarily active in each time step. This aspect further adds another layer of robustness to the dynamical hallmarks of the motifs, thus putatively conferring more robustness to these motifs.
Often, various cellular processes are viewed as generally being robust to noise and perturbations. However, in diseased states such as cancer, owing to a repertoire of changes (network topology, protein production rates, or stability) can potentially perturb these steady states and encourage the progression of disease into states which can be hard to reverse [37,[47][48][49]. Our results show that in addition to the abovementioned factors, the local density of a motif can also deviate from its canonical behavior. This insight can help us develop better algorithms to identify potential drug targets which are more susceptible to alterations in their local neighborhood, to potentially escape these diseased states [36,50]. Elucidating the dynamics of such multistable network motifs and the ability to modulate their behavior in larger networks can help us better understand and control the diverse trajectories of cellular differentiation to reprogram cells into the desired cell fate(s) [4,36,51].