Evolution of Activity-Dependent Adaptive Boolean Networks towards Criticality: An Analytic Approach

TL;DR

This paper introduces activity-dependent adaptive Boolean networks inspired by gene regulation, providing an analytical framework to understand their evolution towards criticality and matching real-world gene network properties.

Contribution

It presents a novel analytical approach to model adaptive Boolean networks, demonstrating how they can evolve towards criticality with tunable sensitivity, supported by numerical validation.

Findings

Analytical solution for stationary in-degree distribution.

Network rewiring drives networks towards criticality when sensitivity is set to 1.

Good agreement between theory and numerical simulations.

Abstract

We propose new activity-dependent adaptive Boolean networks inspired by the cis-regulatory mechanism in gene regulatory networks. We analytically show that our model can be solved for stationary in-degree distribution for a wide class of update rules by employing the annealed approximation of Boolean network dynamics and that evolved Boolean networks have a preassigned average sensitivity that can be set independently of update rules if certain conditions are satisfied. In particular, when it is set to 1, our theory predicts that the proposed network rewiring algorithm drives Boolean networks towards criticality. We verify that these analytic results agree well with numerical simulations for four representative update rules. We also discuss the relationship between sensitivity of update rules and stationary in-degree distributions and compare it with that in real-world gene regulatory networks.

Equations

$$x_i(t+1)=f_i({\mathbf{x}}_i(t)),$$

$$b_{t+1}=\sum _k{\beta }_k(b_t)P_{\mathrm{in}}(k),$$

$${\tilde{d}}_{t+1}=\lambda {\tilde{d}}_t$$

$$\lambda =\bar{\lambda }=\sum _k{P}_{\mathrm{in}}(k){\lambda }_k.$$

$$P_{\mathrm{in}}(e+1,k)=\left(1-{\pi }_n-{\pi }_a{\lambda }_k\right)P_{\mathrm{in}}(e,k)+{\pi }_n{P}_{\mathrm{in}}(e,k-1)+{\pi }_a{\lambda }_{k+1}{P}_{\mathrm{in}}(e,k+1)$$

$$P_{\mathrm{in}}(e+1,0)=\left(1-{\pi }_n\right)P_{\mathrm{in}}(e,0)+{\pi }_a{\lambda }_1{P}_{\mathrm{in}}(e,1).$$

$${\pi }_n{P}_{\mathrm{in}}^s(k)={\pi }_a{\lambda }_{k+1}{P}_{\mathrm{in}}^s(k+1)$$

$$\lambda ={\pi }_n/{\pi }_a$$

$$P_{\mathrm{in}}^s(k)=P_{\mathrm{in}}^s(0){\sigma }^k{\left(\prod _{l=1}^k{\lambda }_l\right)}^{-1}.$$

$$f({\mathbf{x}}_i)=\{\begin{array}{ll}{s}_1& \text{if xj1=c1}\\ s_2& \text{if xj1=c1 and xj2=c2}\\ s_3& \text{if xj1=c1 and xj2=c2 and xj3=c3}\\ ⋮& \\ s_{k_i}& \text{if xj1=c1 and …and xjki=cki}\\ s_d& \text{otherwise}\end{array}$$

Full text

1 Introduction

Boolean networks (BNs) [1] were originally proposed as a model of gene regulatory networks (GRNs) by S. Kauffman in 1969 [2]. Since then, it also has been used as a useful representation for modeling other complex systems such as neuronal networks [3] and social networks [4]. Although Boolean abstraction of real-world complex systems ignores fine details of them, it enables us to study some important aspects of their generic features. One such feature of BNs is the phase transition between ordered phase and disordered phase [5]. It has often been argued, but still is controversial, that real-world living systems such as GRNs and neuronal networks adjust their dynamical behavior towards the boundary between the two phases, criticality [6, 7, 8, 9, 10]. The advantages of criticality also have been studied: Optimal computational ability [11, 12], maximal sensitivity to external stimuli [13], maximal memory capacity [14] and so on.

So far, many plausible theoretical models of network self-organization towards criticality have been proposed although the exact mechanisms in gene regulatory or neuronal networks have not yet been known. For example, the following literatures discuss biologically inspired mechanisms: Hebbian learning [15, 16], spike-timing dependent plasticity [17, 18], dynamical synapses [19, 20] and homeostatic plasticity [21] for neuronal networks and local control of feedback loops [22] and adaptation towards both adaptability and stability [23] for gene regulatory networks. Such models have been collectively called adaptive networks, in which network structure and network state coevolve, and have been paid much attention recently [24, 25]. The study of adaptive networks originates from the work by Bornholdt and Rohlf [26], which is also motivated by a preliminary work on the relationship between network structure and network state [27]. They showed that a simple activity-dependent rewiring rule based on measurement of local dynamics drives random threshold networks towards criticality by numerical simulation. The model has been extended to different situations: Liu and Bassler [28] reported that the activity-dependent rewiring rule drives random Boolean networks towards criticality by numerical simulation. Recently, this model was extended to networks with modular structure [29]. In spite of the structural constraint, self-organization towards criticality was shown to be preserved. Rohlf introduced an activity-dependent threshold change into the original Bornholdt-Rohlf model [30]. In his model, the threshold change and rewiring are switched stochastically. It was shown that the adaptive thresholds yield a new class of self-organized networks. However, it was confirmed numerically that networks still evolve towards criticality in the large size limit. In summary, these previous works based on numerical simulation suggest that the activity dependent rewiring rule robustly drives networks towards criticality in different conditions.

However, in these adaptive Boolean network models, no analytic approach has been reported so far to the best of the author’s knowledge. One reason for this would be the fact that the definition of activity is dependent on attractors which are usually avoided to discuss the phase transition of Boolean networks in the limit of large system size [5]. In the activity-dependent rewiring rule of Bornholdt and Rohlf [26], a node is defined to be active if it does not change its state on the attractor reached from a random initial condition. Otherwise, the node is said to be static. The rewiring rule is as follows: The active node loses one of its incoming link randomly and the static node acquires a new incoming link randomly. Indeed, it seems that the activity on attractors is crucial for self-organization towards criticality. Bornholdt and Rohlf [26] numerically identified a first-order-like transition of the frozen component defined as the fraction of static nodes and argued that this transition is the main mechanism of robust self-organization of networks towards criticality.

In this paper, we propose a new activity-dependent adaptive Boolean network model inspired by the cis-regulatory mechanism of real-world GRNs in which activity does not dependent on attractors but is defined by typical states that will be defined in Sec. 2. By this change of the definition of activity, we expect that our model admits analysis based on a mean-field theory called the annealed approximation in the limit of large system size. In the following, we show that our model can be solved for stationary in-degree distribution for a wide class of update rules to which the annealed approximation of the Boolean dynamics is applicable. At first sight, one would suspect that our network rewiring rule is designed towards a desired result, namely, criticality. However, it turns out that whether our model can self-organize towards criticality depends on a parameter of our network rewiring rule independent of update rules. We analytically show that the average sensitivity of stationary BN dynamics is equal to the parameter if certain conditions are satisfied. Thus, only when the value of the parameter is set to $1$, we expect that BNs evolve towards criticality. The analytic result is verified by numerical simulation in four representative update rules. We also discuss the relationship between sensitivity of update rules and the tail of stationary in-degree distributions and compare it with that in real-world GRNs.

2 Model

Boolean networks (BNs) consist of a directed network with $N$ nodes that can take two states [math] and $1$. The state of node $i$ at time step $t$ is denoted by $x_{i}(t)$ and is updated by a rule $f_{i}$ selected from a given ensemble of Boolean functions $\mathcal{E}_{i}$:

[TABLE]

where ${\bm{x}}_{i}(t)=(x_{j_{1}}(t),\dots,x_{j_{k_{i}}}(t))$ and $j_{1},\dots,j_{k_{i}}$ are nodes from which node $i$ receives inputs. The number of inputs $k_{i}$ is called in-degree of $i$. In this paper, all nodes are updated simultaneously. We also assume that the ensemble of Boolean functions $\mathcal{E}_{i}$ associated with node $i$ only depends on its in-degree $k_{i}$.

Our activity-dependent rewiring rule for network evolution is different from those proposed in previous work [26, 28] in the following two respects. First, both nodes and arcs can be selected at each time step of network evolution, in contrast to the previous models where only nodes are assumed to be selected. Second, we consider activity of arcs rather than that of nodes. In the previous models, activity of a selected node is measured by time-averaging its state value along a reached attractor and the decision whether the selected node gets a new incoming arc or loses an existing arc is made depending on the value of activity. In our model, when a node is selected, the node gets a new incoming arc. On the other hand, when an arc is selected, it is deleted when it is active. Here, activity of the arc is evaluated by the response of the target node $i$ to perturbations on the arc given a typical state. That is, given an input ${\bm{x}}_{i}=(x_{j_{1}},\dots,x_{j_{k_{i}}})$ sampled randomly from a collection of states after sufficiently long time steps starting from a random initial condition, the arc is said to be active if $f_{i}({\bm{x}}_{i})\neq f_{i}(\tilde{{\bm{x}}}_{i})$ where $\tilde{{\bm{x}}}_{i}=(\tilde{x}_{j_{1}},\dots,\tilde{x}_{j_{k_{i}}})$ is given by $\tilde{x}_{j_{l}}=1-x_{j_{l}}$ if $j_{l}$ is the source of the selected arc and $\tilde{x}_{j_{l}}=x_{j_{l}}$ otherwise. These modifications are motivated by the following biological consideration: Deletion of an arc in a GRN of an organism can be caused by mutations in cis-regulatory elements (CREs) [31, 32] of a gene that are nearby non-coding regions of DNA where a number of proteins called transcription factors (TFs) that are themselves products of other genes can bind. TFs regulate expression of the gene by increasing or decreasing the frequency of transcription initiation. If mutations in existing CREs of a gene change the binding pattern of TFs and the expression level of the gene, it could result in undesirable behavior of the organism and the corresponding arcs in its gene regulatory network are deleted in an evolutionary time scale [33]. On the other hand, mutations in a non-coding region of DNA within functional interaction range that is not involved in existing CREs could give rise to binding of a new TF. This means addition of a new incoming arc to the node representing the gene. Thus, nodes in a GRN can be conceived as carrying capacity to accept new incoming arcs incarnated by non-coding regions of DNA rather than coding DNA. In summary, when considering rewiring of a GRN, it is natural to treat nodes and arcs on the same footing because the physical basis of them is the same.

In detail, our algorithm for network evolution in this paper is as follows:

1. (i)

An initial BN with a given ensemble of Boolean functions is generated. The in-degree of each node is sampled from a Poisson distribution with mean $k_{0}$ and the source of each arc is chosen randomly. 2. (ii)

The state of the BN is evolved from a random initial state for sufficiently long time steps to find a typical state. For any BN of finite size $N$, its state trajectory eventually falls onto an attractor. Hence, it is ideal to choose a state randomly from the attractor. However, when numerically simulating the model, it is difficult to find an attractor in a reasonable time if the BN is in the disordered phase. For efficient numerical simulation, we limit the maximum length of attractors to be detected as $T$. If no attractor is found within $2T+T^{\prime}$ time steps, the last $T$ steps are stored and a state is chosen randomly from the $T$ states. In this paper, we set $T=1000$ and $T^{\prime}=100$. We expect that this way of sampling a state approximates that of sampling from true typical states in the limit of large $N$ because correlations between nodes are negligible for $N\gg 1$ if the underlying network is locally tree-like and thus whether a state is on an attractor or not does not matter if it is reached after many time steps from a random initial state [1]. Indeed, this expectation accommodates to the assumptions of the mean-field theory used in Sec. 3 and we will see that the numerically obtained in-degree distributions by this network rewiring algorithm agree well with the theoretical predictions based on the mean-field theory. 3. (iii)

A particular node or a particular arc is chosen with probability $\pi_{n}$ or $\pi_{a}$, respectively. Here, we fix the ratio $\sigma:=\pi_{n}/\pi_{a}$ throughout the network evolution. If a node is chosen, then a new incoming arc is added to the node. The source of the new arc is chosen uniformly at random. If an arc is chosen, then its activity in the state chosen in step (ii) is assessed. If the arc is active, then it is deleted. Otherwise, do nothing. 4. (iv)

The Boolean function on the chosen node or the target of the chosen arc in step (iii) is re-assigned following the given ensemble of Boolean functions. 5. (v)

Go back to step (ii).

The steps (ii)-(v) constitute time unit of network evolution. We call it epoch after [28]. Note that $\pi_{n}N+\pi_{a}z(e)N=1$ should hold for all epoch $e$ where $z(e)$ is the average in-degree of the underlying directed network of BN at epoch $e$. Thus, $\pi_{n}=\sigma/[(\sigma+z(e))N]$ and $\pi_{a}=1/[(\sigma+z(e))N]$.

In each epoch, the network topology and Boolean functions assigned are fixed as in typical applications of BNs for modeling real-world complex systems. Thus, in the above model, the time scale separation between BN dynamics and network evolution is taken for granted.

3 Analytic results

In this section, first we develop a general mean-field theory of network evolution that can be applied to any update rule which satisfies certain conditions mentioned below. Second, we apply the analytic result derived from the mean-field theory to four update rules that have been paid attention in the literature.

3.1 Mean-field theory

If the large system size limit $N\to\infty$ is taken and the underlying directed network is random networks with a specified degree distribution $P(k,l)$ [34], where $P(k,l)$ is the probability that a randomly chosen node has in-degree $k$ and out-degree $l$, the stability of BN dynamics can be analyzed by a mean-field theory so-called annealed approximation [5, 35]. In the annealed approximation, correlations between nodes are neglected. This is manifested as the following ansatz taken in the mean-field calculation of BN dynamics [1]: The sources of incoming arcs to a node are chosen randomly at each time step and the Boolean functions are also re-assigned randomly at each time step.

We apply the annealed approximation to BN dynamics in each epoch and assess its stability. For this purpose, we need to calculate sensitivity of Boolean functions selected from a given ensemble for each input [36]. Let $\lambda_{k,j}$ be the probability that the output of an assigned Boolean function with $k$ inputs changes when $j-$th input is flipped for $1\leq j\leq k$. We put $\lambda_{k}:=\sum_{j=1}^{k}\lambda_{k,j}$. In general, $\lambda_{k,j}$ depends on the fraction $b_{t}$ of nodes with state $1$ at time step $t$. $b_{t}$ evolves by the following equation

[TABLE]

where $\beta_{k}(b_{t})$ is the probability that the output of a node with $k$ inputs is $1$ at time step $t+1$ and $P_{\rm in}(k)=\sum_{l}P(k,l)$ is the in-degree distribution. Although Eq. (2) can have periodic or chaotic solutions depending on update rules [1], we only consider the case that Eq. (2) has a unique stable stationary solution $b^{*}$ in the following.

Now let us suppose that the dynamics of a BN settle down to the stationary regime and apply a small perturbation. Let $\tilde{d}_{t}$ be the fraction of damaged inputs at time step $t$. That is, $\tilde{d}_{t}$ is the probability that the source node of a randomly chosen arc is flipped. Neglecting the higher order terms of $\tilde{d}_{t}$, we obtain

[TABLE]

for the time evolution of $\tilde{d}_{t}$ by a similar reasoning with previous work [35, 37], where $\lambda=\sum_{k,l}\frac{lP(k,l)}{z}\lambda_{k}$ which we call average sensitivity, $z=\sum_{k}kP_{\rm in}(k)$ is the average in-degree and $\lambda_{k}$ is evaluated at $b^{*}$. Let $d_{t}$ be the fraction of damaged nodes at time step $t$. Since $d_{t+1}=\bar{\lambda}\tilde{d}_{t}$ where $\bar{\lambda}=\sum_{k}P_{\rm in}(k)\lambda_{k}$, $d_{t}$ also follows Eq. (3). When in-degree and out-degree are independent as we expect for networks evolved by the proposed network rewiring algorithm, we have

[TABLE]

When $\lambda<1$, $d_{t}$ dies out eventually and the dynamics are said to be ordered or stable. If $\lambda>1$, $d_{t}$ grows exponentially at first and the dynamics are said to be disordered or unstable. $\lambda=1$ is the boundary between the two cases and the dynamics are said to be critical.

Now let us write down the equation for the time evolution of in-degree distribution by assuming the annealed approximation for the dynamics of BN at each epoch. Let $P_{\rm in}(e,k)$ be the in-degree distribution at epoch $e$. According to the proposed network rewiring algorithm, we have

[TABLE]

for $k\geq 1$ and

[TABLE]

In order to iteratively solve Eqs. (5) and (6), in each iteration one must calculate $\lambda_{k}$ which is in general a function of $b^{*}$, which in turn depends on the entire in-degree distribution at epoch $e$ through Eq. (2). In addition, $\pi_{n}$ and $\pi_{a}$ are functions of average in-degree $z(e)$. A stationary solution $P_{\rm in}^{s}(k)$ of Eqs. (5) and (6) should satisfy

[TABLE]

for $k\geq 0$ if it exists. When the stationary solution exists, we obtain

[TABLE]

by substituting Eq. (7) into Eq. (4). Thus, we predict that we can control the stability of evolved BNs by adjusting the ratio $\sigma=\pi_{n}/\pi_{a}$ which we call target average sensitivity (TAS) hereafter. Note that $\sigma$ can be given independently of update rules. In particular, when $\sigma=1$, that is, when a node or an arc is selected uniformly at random, the proposed network rewiring algorithm is expected to drive BNs towards criticality.

The limitation of our mean-field theory arises from the normalization condition for the stationary in-degree distribution. If Eq. (7) has a solution, it is solved by

[TABLE]

Hence the infinite series $\sum_{k=0}^{\infty}r_{k}$ must be convergent, where $r_{k}=\sigma^{k}\left(\prod_{l=1}^{k}\lambda_{l}\right)^{-1}$. Since $r_{k+1}/r_{k}=\sigma/\lambda_{k+1}$, this is always the case when $\lambda_{k}$ diverges as $k\to\infty$ by d’Alembert’s ratio test. However, when $\lambda_{k}$ converges to a number $\alpha$ as $k\to\infty$, it must hold that $\sigma\leq\alpha$. When $b^{*}$ is independent of $P_{\rm in}^{s}$, we can give the condition for the existence of $P_{\rm in}^{s}$ as follows: (i) If $\lambda_{k}\to\infty$ as $k\to\infty$, then $P_{\rm in}^{s}$ exists. (ii) If $\lambda_{k}\to\alpha<\infty$ as $k\to\infty$, then $P_{\rm in}^{s}$ exists if $\sigma<\alpha$. If $\sigma>\alpha$, then $P_{\rm in}^{s}$ does not exist. If $\sigma=\alpha$, then the existence of $P_{\rm in}^{s}$ depends on the precise form of $\lambda_{k}$. Even when $P_{\rm in}^{s}$ does not exist in the mean-field theory, we can formally obtain $P_{\rm in}^{s}$ by truncating Eq. (9) at $k=N$ for BNs of finite size $N$. However, it is not guaranteed that the truncated $P_{\rm in}^{s}$ can reproduce the stationary in-degree distribution of the evolved finite size BNs. This is because the assumption of the absence of correlations between nodes in the annealed approximation of BN dynamics will be violated in such case due to the existence of non-negligible amount of nodes with in-degree proportional to system size $N$.

3.2 Examples

In this subsection, we apply the analytic result presented in Sec. 3.1 to four ensembles of Boolean functions: (a) Biased functions (BF) [5]: All Boolean functions with $k_{i}$ inputs are weighted with bias $p$. The value of output of $f_{i}$ is assigned to be $1$ with probability $p$ or [math] with probability $1-p$ for each input ${\bm{x}}_{i}$. (b) Threshold functions (TF) [38]: Only threshold functions are considered. $f_{i}({\bm{x}}_{i})=1$ if $\sum_{l=1}^{k_{i}}w_{j_{l}i}(2x_{j_{l}}-1)+h_{i}\geq 0$ or [math] otherwise, where ${\bm{x}}_{i}=(x_{j_{1}},\dots,x_{j_{k_{i}}})\in\{0,1\}^{k_{i}}$ and $w_{j_{l}i}=\pm 1$ with equal probability. In the following, we only consider the case $h_{i}=0$ for all $i$. (c) Heterogeneous biased functions (HBF) [39]: In this update rule, we allow the bias of BFs to depend on in-degree. That is, a BF with bias $p_{k_{i}}$ is selected for node $i$ with in-degree $k_{i}$. (d) Nested Canalizing functions (NCF) [40]: A nested canalizing function is given by

[TABLE]

for ${\bm{x}}_{i}=(x_{j_{1}},\dots,x_{j_{k_{i}}})\in\{0,1\}^{k_{i}}$, where $c_{l}\in\{0,1\}$ is the canalizing value for input from node $j_{l}$ and $s_{l}\in\{0,1\}$ is the corresponding output value for $l=1,\dots,k_{i}$. Here, we consider a weight on NCFs defined by the following parameters [41]: $s_{l}=1$ with probability $a$ and $c_{l}=1$ with probability $c$ for $l=1,\dots,k_{i}$, and $s_{d}=1$ with probability $d$.

The formula of $\lambda_{k}$ for BFs, TFs and HBFs are given by $\lambda_{k}=2p(1-p)k$, $\lambda_{k}=k2^{-(k-1)}\binom{k-1}{\lfloor k/2\rfloor}\sim\sqrt{2/\pi}\sqrt{k}$ [38] and $\lambda_{k}=2p_{k}(1-p_{k})k$, respectively. For these three rules, $\lambda_{k}$ is independent of $b^{*}$. However, $\lambda_{k}$ of NCFs depends on $b^{*}$. We have $\beta_{k}(b_{t})=a+(d-a)(1-\gamma(b_{t}))^{k}$ in Eq. (2) where $\gamma(b_{t})=b_{t}c+(1-b_{t})(1-c)$ is the probability that a randomly chosen input is at its canalizing value [41]. $\lambda_{k}$ of NCFs is shown to be $\lambda_{k}=(1-\eta)(1-(1-\gamma(b^{*}))^{k})/\gamma(b^{*})+k(1-\gamma(b^{*}))^{k-1}(\eta-\eta_{0})\sim(1-\eta)/\gamma(b^{*})$ when $0<\gamma(b^{*})<1$, where $\eta=a^{2}+(1-a)^{2}$ and $\eta_{0}=ad+(1-a)(1-d)$ at stationarity [41].

By substituting $\lambda_{k}$ into the right-hand side of Eq. (9), we obtain stationary in-degree distributions. For BFs, we get a Poisson stationary in-degree distribution $P_{\rm in}^{s}(k)=e^{-z_{s}}z_{s}^{k}/k!$ with the stationary average in-degree $z_{s}=\sigma/[2p(1-p)]$. The tail of the stationary in-degree distribution for TFs decays slower than that of any Poisson distribution but does faster than that of any exponential distribution. HBFs have different stationary in-degree distributions depending on the functional form of $p_{k}$ if it exists. For NCFs, the stationary in-degree distribution exists and is asymptotically equal to an exponential distribution provided that $0<\gamma(b^{*})<1$ and $\sigma<(1-\eta)/\gamma(b^{*})$ where $b^{*}$ satisfies $b^{*}=\sum_{k}\beta_{k}(b^{*})P_{\rm in}^{s}(k)$.

In next section, we test these analytic predictions for TAS $\sigma$ close to $1$ since our primary interest is evolution towards criticality. The behavior of our model for a wider range of $\sigma$ is investigated in Appendix where we also present an example in which our mean-field theory fails.

4 Numerical results

We compared analytic results with numerical simulations for the above four ensembles of Boolean functions. We simulated evolution of BNs with $N=200$ for three different values of TAS: $\sigma=0.95,1.00$ and $1.05$. Parameters used are $p=0.7$ for BFs, $p_{k}=(1+\sqrt{1-2q_{k}})/2$ with $q_{k}=1/2$ if $1\leq k\leq 3$ and $q_{k}=2/k$ if $k\geq 4$ for HBFs (thus, we have $\lambda_{k}=2$ for $k\geq 4$) and $a=1/3$, $c=0.95$ and $d=0$ for NCFs. The condition for the existence of the stationary in-degree distribution for HBFs is $\sigma<2$ and is satisfied in the numerical simulation here. For NCFs, we numerically checked that $0<\gamma(b^{*})<1$ and $\sigma<(1-\eta)/\gamma(b^{*})$ hold for the above parameter values.

Fig. 1 shows time evolution of the average sensitivities for each update rule from five different initial average in-degree $1\leq k_{0}\leq 5$. For each pair of values of $\sigma$ and $k_{0}$, $100$ realizations were averaged. In Fig. 1, the average sensitivity of a BN at epoch $e$ was calculated by Eq. (4) with a numerical in-degree distribution at epoch $e$ and analytic values of $\lambda_{k}$. We can clearly see that the average sensitivities approach to given values of $\sigma$ independent of $k_{0}$.

The numerical stationary in-degree distributions agree well with the theoretical predictions (Eq. (9)) for all three values of TAS $\sigma$ (Fig. 2). Here, they were obtained by averaging numerical in-degree distributions over last $10000$ epochs in Fig. 1 of $100$ realizations for each $k_{0}$.

Finally, we verified numerically that Eq. (3) (with replacing $\tilde{d}_{t}$ and $\tilde{d}_{t+1}$ by $d_{t}$ and $d_{t+1}$, respectively) holds in evolved BNs for all three values of TAS $\sigma$ by constructing so-called Derrida plots (Fig. 3) [42]. Derrida plots show the fraction of damaged nodes $d_{t+1}$ at time step $t+1$ as a function of the fraction of damaged nodes $d_{t}$ at time step $t$. In Fig. 3, the value of $d_{t+1}$ was averaged over $200$ states of $500$ realizations of evolved BNs (those at the last step in Fig. 1) for each value of $d_{t}$. We can see that for all three values of TAS, the slope at the origin agrees well between numerical calculations and theoretical predictions. In constructing Derrida plots numerically, a subtlety arises when $\lambda_{k}$ depends on $b^{*}$ as in case of NCFs. For BFs, TFs and HBFs, we can choose a random state and randomly flip its fraction of $d_{t}$ nodes to compute $d_{t+1}$ because $\lambda_{k}$ is independent of $b^{*}$ in these update rules. On the other hand, for NCFs, we must choose a typical state and then randomly flip its fraction of $d_{t}$ nodes. It was predicted that this procedure produces the correct slope at the origin of Derrida plots [43]. However, in order for a Derrida plot to be correct for larger values of $d_{t}$, the perturbed state must also be a random sample of typical states (This does not guarantee that the Derrida plot is correct over all the range of $d_{t}$ as shown in [43]). Here, we are interested in only the slope of the Derrida plots at the origin. Hence, it suffices for our purpose to adopt the above procedure.

5 Discussion

In this paper, we proposed a new activity-dependent adaptive Boolean network model and presented its analytic solutions for stationary in-degree distribution by employing the annealed approximation of Boolean dynamics. We showed analytically that stationary BNs evolved by the proposed network rewiring algorithm have in-degree distributions whose average sensitivity is equal to TAS if certain conditions are satisfied and verified the analytic solutions agree well with numerical simulations for four representative update rules. We emphasize that TAS can be given independently of update rules. In particular, if it is set to $1$, our mean-field theory predicts that BNs evolve towards criticality.

In previous work [44, 45], network self-organization towards criticality has been explained by the self-organized criticality picture [46, 47]. That is, criticality is achieved by slowly adding links in the subcritical phase and rapidly deleting links in the supercritical phase of an absorbing transition of network activity. In particular, Droste et al. [21] analytically demonstrated this mechanism based on the pair-approximation of the network activity dynamics. They showed that two different time-scale separations are necessary to realize self-organization towards criticality: one is that between state dynamics on networks and topological changes of networks and the other is that between deletion of links and addition of links. In our model, the former time-scale separation is incorporated. However, the latter does not hold because the ratio of the probability of link addition to that of link deletion is finite. Thus, the self-organized criticality picture seems not to hold. In our model, the criticality is realized by stochastically balancing the mutually opposed processes, addition and deletion of links.

In previous work on activity-dependent adaptive Boolean networks, influence of the update rule on the structure of evolved networks is assessed by only numerical simulations [28, 33]. In our model, we have a simple relationship between the sensitivity of update rules represented by $\lambda_{k}$ and the stationary in-degree distribution as shown above. Although our model is parsimonious, it is worth to compare our result with real-world GRNs. The in-degree distribution of the prokaryote Escherichia coli is best fitted by a Poisson distribution, whereas that of the eukaryote Saccharomyces cerevisiae is best fitted by an exponential distribution [48]. As for update rules, NCFs were introduced to model the yeast GRN [40] because NCFs are found abundantly in eukaryotic GRNs by an extensive literature study [49]. On the other hand, the analysis by Balleza et al. [8] suggested that BFs are enough to model the GRN of E. coli. They modeled several real-world GRNs including the bacterium GRN by biased functions to reveal whether they operate close to criticality or not and showed that changes in the fraction of canalizing functions for genes with at least 4 inputs do not affect the near critical dynamical behavior of the bacterium GRN. On the other hand, most of genes in the bacterium GRN have at most 3 inputs and canalizing functions are abundant just by chance for such genes [8]. Thus, there is no need for the bacterium to bias the sampling strategy of update rules towards canalizing functions even if they have an evolutionary advantage. Our model predicts Poisson and exponential stationary in-degree distributions for BFs and NCFs, respectively, and thus is consistent with the real-world GRNs.

We are almost ignorant of out-degree distributions in this paper. Under the proposed network rewiring algorithm, the stationary out-degree distribution becomes a Poisson distribution independent of update rules. This disagrees with real-world GRNs because they have heavy-tailed out-degree distributions [48]. However, we can control the shape of stationary out-degree distribution by modifying step (iii) of the algorithm without changing the value of average sensitivity: selecting the source of a new arc following an appropriate weight depending on the out-degree of each node [50].

Finally, we note that it is an interesting open question whether our model can be extended to the network ensembles to which the semi-annealed approximation of Boolean dynamics [39, 51] is applicable.

Acknowledgments

This work was partially supported by JSPS KAKENHI Grant Number 25280091. The author thanks the anonymous reviewers for their helpful comments to improve the manuscript.

Appendix

In this appendix, we compare our theoretical results with numerical simulation for BFs and HBFs for TAS $\sigma$ apart from criticality. The parameters of the update rules are the same as those in Sec. 4. Our theory predicts that $P_{\rm in}^{s}$ exists for any $\sigma$ for BFs, while exists only for $\sigma<2$ for HBFs since $\lambda_{k}=2$ for large $k$. The condition of numerical simulation is the same as that in Sec. 4 except that we only show results for $k_{0}=3$ here.

In Fig. 4, time evolution of the average sensitivity for BFs and HBFs is shown. We can see that the average sensitivity approaches to each specified value of $\sigma$ except $\sigma=2.0$ for HBFs. The failure of evolution of BNs with HBFs towards TAS $\sigma=2.0$ can also be seen from its in-degree distribution (Fig. 5) and the Derrida plot (Fig. 6). On the other hand, theoretical predictions and results of numerical simulation for the other cases agree well in both Fig. 5 and Fig. 6. These results provide further support for the claim at the end of Sec. 3.1 on the applicability and the limitation of our mean-field theory.