Interplay between degree and Boolean rules in the stability of Boolean networks

The random Boolean network proposed by Kauffman in $1969$ kauffman1969 has been widely used in physics, biology, and computer science for modeling biological regulatory systems in an abstract manner karlebach2008 ; albert2008 ; li2004 ; thomas2001 ; shmulevich2002 . Many functions in living systems can be modeled by Boolean networks, including genetic regulation kauffman1969 ; li2004 ; thomas2001 , neural firing rosin2013 , and social activity sayama2013 . The dynamical patterns of Boolean networks fall into two phases, namely stable and unstable (chaotic) phases. In a stable phase, most nodes rapidly reach a steady state and remain unchanged. In contrast, in an unstable phase most nodes change their states in a chaotic manner. It has been suggested that many biological regulatory systems ranging from genetic systems to neural systems tend to hover near the borderline between these two phases, achieving both stability and evolvability kauffmanbook ; munoz2018 ; krotov2014 . Empirical evidence supports the hypothesis that biological networks remain near criticality kauffmanbook ; daniels2018 , especially for knockout experiments for single genes in Saccharomyces cerevisiae serra2004 and gene expression dynamics in macrophages nykter2008 .

Theoretical predictions regarding the dynamics of Boolean networks are based on stability against perturbations derrida1986 ; luque1997 . While damage caused by perturbations dies out quickly in a stable phase, it can spread through an entire system in an unstable phase. Pioneering research on the theoretical prediction of the stability of random Boolean networks has revealed that the mean degree $\langle k\rangle$ of a network and the mean bias $\langle p\rangle$ of Boolean functions typically determine the location of criticality derrida1986 ; luque1997 . Since then, many studies have attempted to assess the effects of structural and dynamical features, including scale-free structures aldana2002 ; lee2008 , noise peixoto2009 ; villegas2016 , multi-level interactions cozzo2012 , asynchronous updates klemm2005 , continuous dynamics ghanbarnejad2011 , veto functions ebadi2014 , bipartite interactions lee2012 , knockout boldhaus2010 ; wang2018 , adaptive dynamics goudarzi2012 , and canalizing functions moreira2005 .

Recently, an exhaustive search of real-world networks revealed that biological regulatory circuits achieve stable and adaptive functionality based on the interplay between logical variables and causal structures daniels2018 . To be specific, anti-correlated sensitivity to the local degree of nodes and abundant canalizing functions were observed in the real-world biological networks across diverse regulatory systems daniels2018 ; moreira2005 . Therefore, it is expected that real-world biological networks achieve near-critical behavior in a distinguished way than random Boolean networks without any correlation between Boolean rules and topology. There have been several attempts to study the role of the correlations between structural and dynamic properties in Boolean networks shmulevich2004 ; squires2014 ; pomerance2009 , but this role is still not fully understood, particularly for the correlations between node degree and Boolean functions.

In this study, we analyze the stability of a random Boolean network incorporating interplay between network topology and Boolean variables. Specifically, we aim to assess the effects of the correlation between local node degree and Boolean functions in terms of sensitivity shmulevich2004 ; squires2014 and canalizing inputs moreira2005 . In this paper, we elucidate the role of the coupling between local topology and Boolean rules in promoting the stability of Boolean networks. Specifically, we demonstrate that negatively correlated sensitivity to degree enhances the stability of Boolean networks. Additionally we find that coordination between high-degree nodes and canalizing inputs can enhance stability. Numerical simulations are conducted to verify our analytical predictions, revealing excellent agreement.

We analyze the effects of the correlation between node degree and sensitivity using the general framework described above. First, we constructed an Erdös-Rényi (ER) graph and assign the sensitivity $q_{i}=2p_{i}(1-p_{i})$, where $p_{i}$ is the bias. Here we used a directed version of ER graphs meaning that we assign directionality randomly at each link. In addition, we assume that the degree distribution $P(k^{in},k^{out})$ of ER graphs follows the Poisson distribution as the expected distribution for $N\rightarrow\infty$. We consider three representative cases of the coupling between sensitivity and node in-degree: uncorrelated (UC), positively correlated (PC), and negatively correlated (NC). For UC case, we assign the same $\langle p\rangle$ to each node to obtain uncorrelated coupling. For the PC case, we assign the linearly correlated bias $p_{i}$ of node $i$ to its in-degree $k_{i}^{in}$ as $p_{i}=C_{P}k_{i}^{in}$. Here, $C_{P}$ determines the average bias for a given mean in-degree as $\langle p\rangle=C_{P}\langle k^{in}\rangle$. In contrast, for the NC case, $p_{i}$ is assigned as $p_{i}=-C_{N}k_{i}^{in}+1/2$, where $C_{N}$ determines the average bias as $\langle p\rangle=-C_{N}\langle k^{in}\rangle+1/2$. The factor of $1/2$ ensures that the maximum value of bias is $1/2$. The exact linear relationships in the PC and NC cases do not sustain all possible ranges of $\langle p\rangle$ because $0\leq p_{i}\leq 1/2$. However, the range of linear dependency is still sufficiently broad to examine the impact of the correlation. By substituting all of these parameters into Eq. 7, we can derive the self-consistency equations for $H_{a}$ for the three coupling cases as follows:

By solving these self-consistency equations, we can obtain the average Hamming distances and identify the critical points.

We implemented numerical simulations on ER networks with $\langle k^{in}\rangle=\langle k^{out}\rangle=4$ without any degree-degree correlation. We assigned the biases and corresponding Boolean variables according to the process described above. From initial states selected randomly, the state of each node are updated synchronously according to the Boolean variables. After a transient period, the trajectory will fall into a point or cycle attractor since the dynamics are deterministic. To simulate a perturbation, we flipped a fraction of $0.01$ of the nodes by force. When the system reached a point or cycle attractor again following the perturbation, we measured the Hamming distance over all nodes.

In Fig. 1, we compare the analytical predictions from Eq. 9 to the numerical simulations. The agreement between the theory and the simulations is excellent. We find that a negatively correlated sensitivity to degree enhances stability when comparing the UC and PC cases. For the NC case, the transition point of mean bias $\langle p\rangle_{c}$ is delayed and $H$ is lower than in the other cases. In contrast, the PC case exhibits an enlarged chaotic region compared to the other cases, making it more vulnerable to perturbations. These results demonstrate that the correlation between sensitivity and degree can significantly affect the stability in terms of the location of the critical point and the size of Hamming distance.

In order to examine the effect of the correlation in more realistic circuits, we consider two examples beyond ER networks: i) an empirical regulatory network conroy and ii) a network with a heavy-tailed out-degree distribution. First, we consider a real-world regulatory network of CD4⁺ T-cell with $N=188$ and $\langle k\rangle=3.73$ conroy . Each node and link in the network respectively represents proteins and protein-protein interactions. We assume every link in the network is bi-directional. We choose the CD4⁺ T-cell network as a representative example because the network is the second largest example in the all tested networks in the survey of empirical Boolean networks daniels2018 and also has a sufficient link density to examine the stability of the network. Second, we study networks with scale-free out-degrees and Poisson distributed in-degrees. The rationale behind these networks is that real transcriptional regulatory networks commonly possess a heavy-tailed out-degree distribution and a short-tailed in-degree distribution guelzim ; dobrin . We constructed $10^{3}$ different network realizations with size $N=10^{4}$ having scale-free out-degrees with the degree exponent $\gamma=2.5$ and Poisson distributed in-degrees with $\langle k_{in}\rangle=4.72\ldots$, according to the configuration model. We implemented numerical simulations on the two networks with three couplings as described above.

In Fig. 2, we show the analytical predictions and the numerical simulations (symbols) for the CD4⁺ T-cell network and scale-free out-degree distribution network. The theory (lines) based on Eqs. 2 and 3 shows well agreement with the simulations in a qualitative manner. Deviation from the theory in the CD4⁺ network shown in Fig. 2(a) is due to the short loops and small size of the real network. We find again a negatively correlated sensitivity to degree increases the stability of Boolean dynamics when comparing the UC and PC cases. We confirm the delayed transition $\langle p\rangle_{c}$ to chaotic phase and lower $H$ for the NC coupling than the other cases. On the other hand, the PC coupling shows less stable behaviors to perturbations. There is an inversion between $H$ for UC and PC cases for $\langle p\rangle\gtrsim 0.25$, far above $\langle p\rangle_{c}$ in Fig. 2(b). In this region, low degree nodes with low bias in PC coupling suppress $H$, so that $H$ for PC is above that for UC. The similar crossing behaviors are often observed when one considers degree correlations such as assortative mixing assortative and interlayer degree correlations interlayer .

To evaluate the impact of the interplay between node degree and sensitivity more clearly, we consider a transparent example with a bimodal in-degree distribution $P(k^{in})=(1/2)\delta_{k^{in},K_{1}}+(1/2)\delta_{k^{in},K_{2}}$, where $\delta_{i,j}$ represents the Kronecker delta. We construct networks with the bimodal in-degree and Poisson distributed out-degree distribution according to the configuration model. We assign a bias drawn from a bimodal distribution $Q(p)=(1/2)\delta_{p,\phi_{1}}+(1/2)\delta_{p,\phi_{2}}$. Similar to the analysis above, we consider three types of correlated coupling: UC, PC, and NC. For the UC case, we assign a bias to each node at random, independent of the node degree. Positively (negatively) correlated coupling can be achieved by ensuring that higher (lower) degree nodes have greater bias values. In our examples, we use $K_{1}=10$, $K_{2}=6$, $\phi_{1}=2\langle p\rangle-1/2$, and $\phi_{2}=0.1$, where $0.1\leq\phi_{1}\leq 0.5$. Note that the range of $\langle p\rangle$ is $0.1\leq\langle p\rangle\leq 0.3$. By annealing the probability $H_{i}$, we can derive the self-consistency equation for $H_{a}$ as follows:

where $q_{1}=2\phi_{1}(1-\phi_{1})$ and $q_{2}=2\phi_{2}(1-\phi_{2})$.

As shown in Fig. 3, negatively correlated coupling is more resilient to external perturbations compared to the other types of coupling. The average Hamming distance clearly highlights the effect of the correlation between sensitivity and local node in-degree. Specifically, negative correlation between sensitivity and node degree enhances the global stability of Boolean networks [Fig. 3(a)]. For the NC case, the majority of incoming links are connected to unbiased nodes ($p=1/2$), leading to more stable Boolean dynamics. In contrast, for the PC case, damage can easily spread through an entire network because an adequate fraction of high-degree nodes have high sensitivity. Fig. 3(b) presents a phase diagram as functions of the mean node degree $\langle k^{in}\rangle$ and $\langle p\rangle$, where $K_{1}=\langle k^{in}\rangle+2$ and $K_{2}=\langle k^{in}\rangle-2$. An increasing value of $p_{c}$ for the NC coupling can be observed consistently over a wide range of parameter sets.

In addition to the global stability of Boolean networks, we can also assess the stability of each node in a given network topology using Eq. 2. Fig. 4 reveals perfect agreement between the numerical results $H_{i}^{sim}$ and theoretical predictions $H_{i}^{th}$ for the probability that a node $i$ changes its state when a perturbation occurs. We computed the average Hamming distances $\langle H(k^{in}=10)\rangle$ and $\langle H(k^{in}=6)\rangle$ for nodes with high in-degree ($k^{in}=10$) and low in-degree ($k^{in}=6$), respectively, as indicated by the red and blue arrows in Fig. 4, respectively. The average Hamming distance over all nodes is denoted by a filled circle. One can see two clearly separated groups of nodes with different stability values and Hamming distances $H_{i}$ for the PC coupling. We can confirm that in the PC coupling, damage can spread through high in-degree nodes with high sensitivity, which are prone to instability. However, for the NC coupling, these two groups merge and perturbations terminate quickly. From the perspective of a percolation problem, NC coupling corresponds to the case where high in-degree nodes have low occupation probabilities, leading to stable dynamics, which is analogous to degree-based removal in network percolation attack .

Finally, we consider the role of the interplay between local node in-degree and canalizing inputs moreira2005 ; daniels2018 . Canalizing functions have a single input that forces the corresponding output to a specific value, regardless of the values of other inputs. The Hamming distance $H_{i}$ for each node $i$ with a fraction $c_{i}$ of canalizing inputs is calculated as squires2012

where $H_{ic}$ is the Hamming distance of the neighbor of node $i$ connected by a canalizing input, $q_{j}=2p_{j}(1-p_{j})$, $\partial j/c$ define a set of inputs excluding a canalizing input squires2012 . By definition, canalizing functions lead to stable dynamics moreira2005 because they effectively reduce the sensitivity of non-canalizing inputs connected to nodes shared by canalizing inputs. However, the effects of canalizing inputs are not solely determined by the fraction of canalizing inputs. The topological locations of canalizing links also affect stability, which can be predicted using Eq. 11.

We consider three different correlations between node in-degree and the locations of canalizing inputs, which are again denoted as UC, PC, and NC. For the UC case, the canalizing inputs are distributed randomly. For the PC (NC) case nodes with high (low) in-degree tend to have a canalizing input. For the sake of simplicity, we consider a random network with a bimodal in-degree distribution defined as $P(k^{in})=(1/2)\delta_{k^{in},K_{1}}+(1/2)\delta_{k^{in},K_{2}}$, where $K_{1}=6$ and $K_{2}=2$. In this example, we assume that $1/4$ of the nodes have a canalizing input and each node can have at most one canalizing input. For the UC case, a quarter of nodes chosen at random have a canalizing input. For the PC case, the low in-degree nodes ($k^{in}=2$) do not have a canalizing input and a half of high in-degree nodes ($k^{in}=6$) have a canalizing input. On the other hand, for the NC case only low in-degree nodes have a canalizing input with the probability $1/2$ [see Fig. 5(a)].

As shown in Fig. 5(b), correlation between canalizing inputs and local node in-degree can increase global stability. Specifically, the PC coupling enhances the stability of Boolean networks. In contrast, the NC coupling decreases stability, leading to smaller $p_{c}$ and larger $H$ values. In this example, one can see that correlation between local node degree and canalizing inputs alters Boolean dynamics significantly. When a canalizing input becomes active, all other connections are ineffective. Therefore, for the PC case, a larger fraction of non-canalizing inputs lose their influence on Boolean dynamics to the canalizing inputs. In the NC case, the effects of canalizing inputs are minimized because they only affect low in-degree nodes.

We analyzed the stability of random Boolean networks incorporating dynamic rules, as well as the topological properties of each node. We find that correlation between node degree and Boolean functions plays an important role in determining Hamming distances and critical point. Specifically, negatively correlated sensitivity to the in-degree of each node increases stability. We also find that a correlation between high-degree nodes and canalizing inputs can increase global stability. Our results reveal that analysis based on the naive mean-field approach may fail to predict dynamical consequences in Boolean networks with intertwined structural and functional properties, which are often observed in real-world biological systems. Further study is required to examine the effects of more complex features in Boolean networks, such as loops and feedback in network topologies kinoshita ; cho , and hierarchical dynamics.

This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korean government (MSIT) (NRF-2018R1C1B5044202 and NRF-2020R1I1A3068803).

The data that supports the findings of this study are available within the article.