# Evolution of activity-dependent adaptive Boolean networks towards criticality: an analytic approach

Evolution of activity-dependent adaptive Boolean networks towards criticality: an analytic approach Abstract We propose new activity-dependent adaptive Boolean networks (BNs) inspired by the cis-regulatory mechanism in gene regulatory networks (GRNs). 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 BN dynamics and that evolved BNs 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 BNs 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 GRNs. 1. Introduction Boolean networks (BNs) [1] were originally proposed as a model of gene regulatory networks (GRNs) by Kauffman in 1969 [2]. Since then, it also has been used as a useful representation for modelling 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 behaviour towards the boundary between the two phases, criticality [6–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 GRNs 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 BNs 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 BN 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 BNs 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 article, we propose a new activity-dependent adaptive BN 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 Section 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 the normalization condition for the mean-field solution holds. For example, the condition is satisfied when the parameter is smaller than the sensitivity of a node in the limit of large in-degree and an extra technical assumption holds. 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 BNs consist of a directed network with |$N$| nodes that can take two states |$0$| 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$|⁠: $$x_i(t+1)=f_i({\boldsymbol x}_i(t)), \label{eq:01}$$ (1) where |${\boldsymbol 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 article, 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 |${\boldsymbol 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({\boldsymbol x}_i) \neq f_i(\tilde{{\boldsymbol x}}_i)$| where |$\tilde{{\boldsymbol 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 behaviour 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 article is as follows: (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. (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'$| time steps, the last |$T$| steps are stored and a state is chosen randomly from the |$T$| states. In this article, we set |$T=1000$| and |$T'=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 Section 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. (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. (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. (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 every 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 modelling 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 analysed 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 $$b_{t+1}=\sum_k \beta_k(b_t) P_{\rm in}(k), \label{eq:02}$$ (2) 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 $$\tilde{d}_{t+1}= \lambda \tilde{d}_t \label{eq:03}$$ (3) for the time evolution of |$\tilde{d}_t$| by a similar reasoning with previous work [35, 37], where |$\lambda=\sum_{k,l} \frac{l P(k,l)}{z} \lambda_k$| which we call average sensitivity, |$z=\sum_k k P_{\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 $$\lambda=\bar{\lambda}=\sum_k P_{\rm in}(k)\lambda_k. \label{eq:04}$$ (4) 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 $$P_{\rm in}(e+1,k)= \left( 1 - \pi_n -\pi_a \lambda_k \right) P_{\rm in}(e,k) + \pi_n P_{\rm in}(e,k-1) + \pi_a \lambda_{k+1} P_{\rm in}(e,k+1) \label{eq:05}$$ (5) for |$k \geq 1$| and $$P_{\rm in}(e+1,0)=\left( 1 - \pi_n \right) P_{\rm in}(e,0) + \pi_a \lambda_1 P_{\rm in}(e,1). \label{eq:06}$$ (6) 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 $$\pi_n P_{\rm in}^s(k)= \pi_a \lambda_{k+1} P_{\rm in}^s(k+1) \label{eq:07}$$ (7) for |$k \geq 0$| if it exists. When the stationary solution exists, we obtain $$\lambda=\pi_n/\pi_a \label{eq:08}$$ (8) 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 $$P_{\rm in}^s(k)=P_{\rm in}^s(0) \sigma^k \left( \prod_{l=1}^k \lambda_l \right)^{-1}. \label{eq:09}$$ (9) 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 Section 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 |$0$| with probability |$1-p$| for each input |${\boldsymbol x}_i$|⁠. (b) Threshold functions (TF) [38]: Only threshold functions are considered. |$f_i({\boldsymbol x}_i)=1$| if |$\sum_{l=1}^{k_i} w_{j_l i}(2x_{j_l}-1)+h_i \geq 0$| or |$0$| otherwise, where |${\boldsymbol 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 NCF is given by $$f({\boldsymbol x}_i)= \begin{cases} s_1 & \text{if x_{j_1}=c_1} \\ s_2 & \text{if x_{j_1} \neq c_1 and x_{j_2}=c_2} \\ s_3 & \text{if x_{j_1} \neq c_1 and x_{j_2} \neq c_2 and x_{j_3}=c_3} \\ \vdots & \\ s_{k_i} & \text{if x_{j_1} \neq c_1 and}\,\,\,\dots\,\,\,\text{and x_{j_{k_i}}=c_{k_i}} \\ s_d & \text{otherwise} \end{cases} \label{eq:10}$$ (10) for |${\boldsymbol 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 behaviour 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. Figure 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$|⁠. Fig. 1. View largeDownload slide (Color online) Time evolution of the average sensitivities for (a) BFs, (b) TFs, (c) HBFs and (d) NCFs. Insets are enlarged views from epoch |$20000$| to |$30000$| for the first three update rules and that from |$50000$| to |$60000$| for NCFs. BNs with NCFs were simulated for a longer period because their convergence is slower than the others. Fig. 1. View largeDownload slide (Color online) Time evolution of the average sensitivities for (a) BFs, (b) TFs, (c) HBFs and (d) NCFs. Insets are enlarged views from epoch |$20000$| to |$30000$| for the first three update rules and that from |$50000$| to |$60000$| for NCFs. BNs with NCFs were simulated for a longer period because their convergence is slower than the others. 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$|⁠. Fig. 2. View largeDownload slide (Color online) Comparison between numerical stationary in-degree distributions (symbols) and theoretical stationary in-degree distributions (lines) for (a) BFs, (b) TFs, (c) HBFs and (d) NCFs. Fig. 2. View largeDownload slide (Color online) Comparison between numerical stationary in-degree distributions (symbols) and theoretical stationary in-degree distributions (lines) for (a) BFs, (b) TFs, (c) HBFs and (d) NCFs. 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. Fig. 3. View largeDownload slide (Color online) 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$| for (a) |$\sigma=0.95$|⁠, (b) |$\sigma=1.00$| and (c) |$\sigma=1.05$|⁠. Dotted lines have a slope equal to |$\sigma$|⁠. Insets are enlarged views of the region |$0 \leq d_t, d_{t+1} \leq 0.1$|⁠. Fig. 3. View largeDownload slide (Color online) 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$| for (a) |$\sigma=0.95$|⁠, (b) |$\sigma=1.00$| and (c) |$\sigma=1.05$|⁠. Dotted lines have a slope equal to |$\sigma$|⁠. Insets are enlarged views of the region |$0 \leq d_t, d_{t+1} \leq 0.1$|⁠. 5. Discussion In this article, we proposed a new activity-dependent adaptive BN 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 the normalization condition for the mean-field solution holds. In particular, the condition is satisfied when TAS is smaller than |$\lambda_k$| in the limit of large |$k$| and |$b^*$| is independent of |$P_{\text{in}}^s$|⁠. We 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 BNs, 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 modelled several real-world GRNs including the bacterium GRN by BFs 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 four inputs do not affect the near critical dynamical behaviour of the bacterium GRN. On the other hand, most of genes in the bacterium GRN have at most three 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 article. 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. Acknowledgements The author thanks the anonymous reviewers for their helpful comments to improve the manuscript. Funding JSPS KAKENHI Grant Number 25280091, in part. 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 Section 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 Section 4 except that we only show results for |$k_0=3$| here. In Fig. A1, 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. A2) and the Derrida plot (Fig. A3). On the other hand, theoretical predictions and results of numerical simulation for the other cases agree well in both Figs A2 and A3. These results provide further support for the claim at the end of Section 3.1 on the applicability and the limitation of our mean-field theory. Fig. A1. View largeDownload slide (Color online) Time evolution of the average sensitivity. (a) BFs. The values of TAS |$\sigma$| are |$0.5$|⁠, |$1.0$|⁠, |$1.5$|⁠, |$2.0$| and |$2.5$| from below. (b) HBFs. The values of TAS |$\sigma$| are |$0.5$|⁠, |$1.0$|⁠, |$1.5$| and |$2.0$| from below. Fig. A1. View largeDownload slide (Color online) Time evolution of the average sensitivity. (a) BFs. The values of TAS |$\sigma$| are |$0.5$|⁠, |$1.0$|⁠, |$1.5$|⁠, |$2.0$| and |$2.5$| from below. (b) HBFs. The values of TAS |$\sigma$| are |$0.5$|⁠, |$1.0$|⁠, |$1.5$| and |$2.0$| from below. Fig. A2. View largeDownload slide (Color online) Comparison between numerical stationary in-degree distributions (symbols) and theoretical stationary in-degree distributions (lines) for (a) BFs and (b) HBFs. The theoretical stationary in-degree distribution does not exist for HBFs with |$\sigma=2.0$|⁠. The shown theoretical line is calculated by Eq. (9) and is truncated by |$N=200$|⁠. Fig. A2. View largeDownload slide (Color online) Comparison between numerical stationary in-degree distributions (symbols) and theoretical stationary in-degree distributions (lines) for (a) BFs and (b) HBFs. The theoretical stationary in-degree distribution does not exist for HBFs with |$\sigma=2.0$|⁠. The shown theoretical line is calculated by Eq. (9) and is truncated by |$N=200$|⁠. Fig. A3. View largeDownload slide (Color online) 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$| for (a) BFs and (b) HBFs. Dotted lines have the same slope as corresponding value of TAS |$\sigma$|⁠. Insets are enlarged views of the region |$0 \leq d_t, d_{t+1} \leq 0.1$|⁠. Fig. A3. View largeDownload slide (Color online) 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$| for (a) BFs and (b) HBFs. Dotted lines have the same slope as corresponding value of TAS |$\sigma$|⁠. Insets are enlarged views of the region |$0 \leq d_t, d_{t+1} \leq 0.1$|⁠. References 1. Drossel B. ( 2008 ) Random Boolean networks. Reviews of Nonlinear Dynamics and Complexity ( Schuster H. G. ed.). Weinheim, Germany : Wiley-VCH Verlag GmbH & Co. KGaA . 2. Kauffman S. A. ( 1969 ) Metablic stability and epigenesis in randomly constructed genetic nets. J. Theor. Biol. , 22 , 437 – 467 . Google Scholar Crossref Search ADS PubMed 3. Kürten K. E. ( 1988 ) Critical phenomena in model neural networks. Phys. Lett. A , 129 , 157 – 160 . Google Scholar Crossref Search ADS 4. Paczuski M. , Bassler K. E. & Corral Á. ( 2000 ) Self-organized networks of competing Boolean agents. Phys. Rev. Lett. , 84 , 3185 – 3188 . Google Scholar Crossref Search ADS PubMed 5. Derrida B. & Pomeau Y. ( 1986 ) Random networks of automata: a simple annealed approximation. Europhys. Lett. , 1 , 45 – 49 . Google Scholar Crossref Search ADS 6. Beggs J. M. & Plenz D. ( 2003 ) Neuronal avalanches in neocortical circuits. J. Neurosci. , 23 , 11167 – 11177 . Google Scholar Crossref Search ADS PubMed 7. Petermann T. , Thiagarajan T. C. , Lebedev M. A. , Nicolelis M. A. L. , Chialvo D. R. & Plenz D. ( 2009 ) Spontaneous cortical activity in awake monkeys composed of neuronal avalanches. Proc. Natl. Acad. Sci. USA , 106 , 15921 – 15926 . Google Scholar Crossref Search ADS 8. Balleza E. , Alvarez-Buylla E. R. , Chaos A. , Kauffman S. , Shmulevich I. & Aldana M. ( 2008 ) Critical dynamics in genetic regulatory networks: examples from four kingdoms. PLoS One , 3 , e2456 . Google Scholar Crossref Search ADS PubMed 9. Nykter M. , Price N. D. , Aldana M. , Ramsey S. A. , Kauffman S. A. , Hood L. E. , Yli-Harja O. & Shmulevich I. ( 2008 ) Gene expression dynamics in the macrophage exhibit criticality. Proc. Natl. Acad. Sci. USA , 105 , 1897 – 1900 . Google Scholar Crossref Search ADS 10. Valverde S. , Ohse S. , Turalska M. , Garcia-Ojalvo J. & West B. J. ( 2015 ) Structural determinants of criticality in biological networks. Front. Physiol. , 6 , 127 . Google Scholar Crossref Search ADS PubMed 11. Bertschinger N. & Natschläger T. ( 2004 ) Real-time computation at the edge of chaos in recurrent neural networks. Neural Comput. , 16 , 1413 – 1436 . Google Scholar Crossref Search ADS PubMed 12. Goudarzi A. , Teuscher C. , Gulbahce N. & Rohlf T. ( 2012 ) Emergent criticality through adaptive information processing in Boolean networks. Phys. Rev. Lett. , 108 , 128702 . Google Scholar Crossref Search ADS PubMed 13. Kinouchi O. & Copelli A. M. ( 2006 ) Optimal dynamical range of excitable networks at criticality. Nat. Phys. , 2 , 348 – 352 . Google Scholar Crossref Search ADS 14. Haldeman C. & Beggs J. M. ( 2005 ) Critical branching captures activity in living neural networks and maximizes the number of metastable states. Phys. Rev. Lett. , 94 , 058101 . Google Scholar Crossref Search ADS PubMed 15. Bornholdt S. & Röhl T. ( 2003 ) Self-organized critical neural networks. Phys. Rev. E , 67 , 066118 . Google Scholar Crossref Search ADS 16. Rybarsch M. & Bornholdt S. ( 2014 ) Avalanches in self-organized critical neural networks: a minimal model for the neural SOC universality class. PLoS One , 9 , e93090 . Google Scholar Crossref Search ADS PubMed 17. Meisel C. & Gross T. ( 2009 ) Adaptive self-organization in a realistic neural network model. Phys. Rev. E , 80 , 061917 . Google Scholar Crossref Search ADS 18. Rubinov M. , Sporns O. , Thivierge J.-P. & Breakspear M. ( 2011 ) Neurobiologically realistic determinants of self-organized criticality in networks of spiking neurons. PLoS Comput. Biol. , 7 , e1002038 . Google Scholar Crossref Search ADS PubMed 19. Levina A. , Herrmann J. M. & Geisel T. ( 2007 ) Dynamical synapses causing self-organized criticality in neural networks. Nat. Phys. , 3 , 857 – 860 . Google Scholar Crossref Search ADS 20. Levina A. , Herrmann J. M. & Geisel T. ( 2009 ) Phase transitions towards criticality in a neural system with adaptive interactions. Phys. Rev. Lett. , 102 , 118110 . Google Scholar Crossref Search ADS PubMed 21. Droste F. , Do A.-L. & Gross T. ( 2013 ) Analytical investigation of self-organized criticality in neural networks. J. R. Soc. Interface , 10 , 20120558 . Google Scholar Crossref Search ADS PubMed 22. MacArthur B. D. , Sánchez-García R. J. & Ma’ayan A. ( 2010 ) Microdynamics and criticality of adaptive regulatory networks. Phys. Rev. Lett. , 104 , 168701 . Google Scholar Crossref Search ADS PubMed 23. Lee D.-S. ( 2014 ) Evolution of regulatory networks towards adaptability and stability in a changing environment. Phys. Rev. E , 90 , 052822 . Google Scholar Crossref Search ADS 24. Gross T. & Sayama H. (eds) ( 2009 ) Adaptive Networks: Theory, Models and Applications . Heidelberg : Springer . 25. Sayama H. , Pestov I. , Schmidt J. , Bush B. J. , Wong C. , Yamanoi J. & Gross T. ( 2013 ) Modeling complex systems with adaptive networks. Comput. Math. Appl. , 65 , 1645 – 1664 . Google Scholar Crossref Search ADS 26. Bornholdt S. & Rohlf T. ( 2000 ) Topological evolution of dynamical networks: global criticality from local dynamics. Phys. Rev. Lett. , 84 , 6114 – 6117 . Google Scholar Crossref Search ADS PubMed 27. Christensen K. , Donangelo R. , Koiller B. & Sneppen K. ( 1998 ) Evolution of random networks. Phys. Rev. Lett. , 81 , 2380 – 2383 . Google Scholar Crossref Search ADS 28. Liu M. & Bassler K. E. ( 2006 ) Emergent criticality from coevolution in random Boolean networks. Phys. Rev. E , 74 , 041910 . Google Scholar Crossref Search ADS 29. Górski P. J. , Czaplicka A. & Hołyst A. ( 2016 ) Coevolution of information processing and topology in hierarchical adaptive random Boolean networks. Eur. Phys. J. B , 89 , 33 . Google Scholar Crossref Search ADS 30. Rohlf T. ( 2008 ) Self-organization of heterogeneous topology and symmetry breaking in networks with adaptive thresholds and rewiring. Europhys. Lett. , 84 , 10004 . Google Scholar Crossref Search ADS 31. Peter I. S. & Davidson E. H. ( 2011 ) Evolution of gene regulatory networks controlling body plan development. Cell , 144 , 970 – 985 . Google Scholar Crossref Search ADS PubMed 32. Wittkopp P. J. & Kalay G. ( 2012 ) Cis-regulatory elements: molecular mechanisms and evolutionary processes underlying divergence. Nat. Rev. Genet. , 13 , 59 – 69 . Google Scholar Crossref Search ADS 33. Rohlf T. & Bornholdt S. ( 2009 ) Self-organized criticality and adaptation in discrete dynamical networks. In Adaptive Networks ( Gross T. & Sayama H. eds). Heidelberg, Germany : Springer , pp. 73 – 106 . 34. Newman M. E. J. , Strogatz S. H. & Watts D. J. ( 2001 ) Random graphs with arbitrary degree distributions and their applications. Phys. Rev. E , 64 , 026118 . Google Scholar Crossref Search ADS 35. Lee D.-S. & Rieger H. ( 2007 ) Comparative study of the transcriptional regulatory networks of E. coli and yeast: structural characteristics leading to marginal dynamic stability. J. Theor. Biol. , 248 , 618 – 626 . Google Scholar Crossref Search ADS PubMed 36. Shmulevich I. & Kauffman S. A. ( 2004 ) Activities and sensitivities in Boolean network models. Phys. Rev. Lett. , 93 , 048701 . Google Scholar Crossref Search ADS PubMed 37. Squires S. , Ott E. & Girvan M. ( 2012 ) Dynamical instability in Boolean networks as a percolation problem. Phys. Rev. Lett. , 109 , 085701 . Google Scholar Crossref Search ADS PubMed 38. Rohlf T. & Bornholdt S. ( 2002 ) Criticality in random threshold networks: annealed approximation and beyond. Physica A , 310 , 245 – 259 . Google Scholar Crossref Search ADS 39. Pomerance A. , Ott E. , Girvan M. & Losert W. ( 2009 ) The effect of network topology on the stability of discrete state models of genetic control. Proc. Natl. Acad. Sci. USA , 106 , 8209 – 8214 . Google Scholar Crossref Search ADS 40. Kauffman S. , Peterson C. , Samuelsson B. & Troein C. ( 2003 ) Random Boolean network models and the yeast transcriptional network. Proc. Natl. Acad. Sci. USA , 100 , 14796 – 14799 . Google Scholar Crossref Search ADS 41. Peixoto T. P. ( 2010 ) The phase diagram of random Boolean networks with nested canalizing functions. Eur. Phys. J. B , 78 , 187 – 192 . Google Scholar Crossref Search ADS 42. Derrida B. & Weisbuch G. ( 1986 ) Evolution of overlaps between configurations in random Boolean networks. J. Phys. , 47 , 1297 – 1303 . Google Scholar Crossref Search ADS 43. Kesseli J. , Rämö P. & Yli-Harja O. ( 2006 ) Iterated maps for annealed Boolean networks. Phys. Rev. E , 74 , 046104 . Google Scholar Crossref Search ADS 44. Hesse J. & Gross T. ( 2014 ) Self-organized criticality as a fundamental property of neural systems. Front. Syst. Neurosci. , 8 , 1 – 14 . Google Scholar Crossref Search ADS PubMed 45. Marković D. & Gros C. ( 2014 ) Power laws and self-organized criticality in theory and nature. Phys. Rep. , 536 , 41 – 74 . Google Scholar Crossref Search ADS 46. Bak P. , Tang C. & Wiesenfeld K. ( 1987 ) Self-organized criticality: an explanation of |$1/f$| noise. Phys. Rev. Lett. , 59 , 381 – 384 . Google Scholar Crossref Search ADS PubMed 47. Jensen H. J. ( 1998 ) Self-organized Criticality: Emergent Complex Behavior in Physical and Biological Systems . Cambridge : Cambridge University Press . 48. Aldana M. , Balleza E. , Kauffman S. & Resendiz O. ( 2007 ) Robustness and evolvability in genetic regulatory networks. J. Theor. Biol. , 245 , 433 – 448 . Google Scholar Crossref Search ADS PubMed 49. Harris S. E. , Sawhill B. K. , Wuensche A. & Kauffman S. ( 2002 ) A model of transcriptional regulatory networks based on biases in the observed regulation rules. Complexity , 7 , 23 – 40 . Google Scholar Crossref Search ADS 50. Haruna T. & Tanaka S. ( 2014 ) On the relationship between local rewiring rules and stationary out-degree distributions in adaptive random Boolean network models. In Artificial Life 14: Proceedings of the Fourteenth International Conference on the Synthesis and Simulation of Living Systems ( Sayama H. Rieffel J. Risi S. Doursat R. & Lipson H. eds). Cambridge : MIT Press , pp. 419 – 426 . 51. Squires S. , Pomerance A. , Girvan M. & Ott E. ( 2014 ) Stability of Boolean networks: the joint effects of topology and update rules. Phys. Rev. E , 90 , 022814 . Google Scholar Crossref Search ADS © The author 2018. Published by Oxford University Press. All rights reserved. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model) http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Journal of Complex Networks Oxford University Press

# Evolution of activity-dependent adaptive Boolean networks towards criticality: an analytic approach

Journal of Complex Networks, Volume 6 (6) – Dec 1, 2018
13 pages

Publisher
Oxford University Press
ISSN
2051-1310
eISSN
2051-1329
D.O.I.
10.1093/comnet/cny006
Publisher site
See Article on Publisher Site

### Abstract

Abstract We propose new activity-dependent adaptive Boolean networks (BNs) inspired by the cis-regulatory mechanism in gene regulatory networks (GRNs). 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 BN dynamics and that evolved BNs 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 BNs 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 GRNs. 1. Introduction Boolean networks (BNs) [1] were originally proposed as a model of gene regulatory networks (GRNs) by Kauffman in 1969 [2]. Since then, it also has been used as a useful representation for modelling 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 behaviour towards the boundary between the two phases, criticality [6–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 GRNs 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 BNs 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 BN 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 BNs 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 article, we propose a new activity-dependent adaptive BN 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 Section 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 the normalization condition for the mean-field solution holds. For example, the condition is satisfied when the parameter is smaller than the sensitivity of a node in the limit of large in-degree and an extra technical assumption holds. 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 BNs consist of a directed network with |$N$| nodes that can take two states |$0$| 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$|⁠: $$x_i(t+1)=f_i({\boldsymbol x}_i(t)), \label{eq:01}$$ (1) where |${\boldsymbol 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 article, 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 |${\boldsymbol 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({\boldsymbol x}_i) \neq f_i(\tilde{{\boldsymbol x}}_i)$| where |$\tilde{{\boldsymbol 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 behaviour 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 article is as follows: (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. (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'$| time steps, the last |$T$| steps are stored and a state is chosen randomly from the |$T$| states. In this article, we set |$T=1000$| and |$T'=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 Section 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. (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. (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. (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 every 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 modelling 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 analysed 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 $$b_{t+1}=\sum_k \beta_k(b_t) P_{\rm in}(k), \label{eq:02}$$ (2) 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 $$\tilde{d}_{t+1}= \lambda \tilde{d}_t \label{eq:03}$$ (3) for the time evolution of |$\tilde{d}_t$| by a similar reasoning with previous work [35, 37], where |$\lambda=\sum_{k,l} \frac{l P(k,l)}{z} \lambda_k$| which we call average sensitivity, |$z=\sum_k k P_{\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 $$\lambda=\bar{\lambda}=\sum_k P_{\rm in}(k)\lambda_k. \label{eq:04}$$ (4) 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 $$P_{\rm in}(e+1,k)= \left( 1 - \pi_n -\pi_a \lambda_k \right) P_{\rm in}(e,k) + \pi_n P_{\rm in}(e,k-1) + \pi_a \lambda_{k+1} P_{\rm in}(e,k+1) \label{eq:05}$$ (5) for |$k \geq 1$| and $$P_{\rm in}(e+1,0)=\left( 1 - \pi_n \right) P_{\rm in}(e,0) + \pi_a \lambda_1 P_{\rm in}(e,1). \label{eq:06}$$ (6) 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 $$\pi_n P_{\rm in}^s(k)= \pi_a \lambda_{k+1} P_{\rm in}^s(k+1) \label{eq:07}$$ (7) for |$k \geq 0$| if it exists. When the stationary solution exists, we obtain $$\lambda=\pi_n/\pi_a \label{eq:08}$$ (8) 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 $$P_{\rm in}^s(k)=P_{\rm in}^s(0) \sigma^k \left( \prod_{l=1}^k \lambda_l \right)^{-1}. \label{eq:09}$$ (9) 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 Section 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 |$0$| with probability |$1-p$| for each input |${\boldsymbol x}_i$|⁠. (b) Threshold functions (TF) [38]: Only threshold functions are considered. |$f_i({\boldsymbol x}_i)=1$| if |$\sum_{l=1}^{k_i} w_{j_l i}(2x_{j_l}-1)+h_i \geq 0$| or |$0$| otherwise, where |${\boldsymbol 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 NCF is given by $$f({\boldsymbol x}_i)= \begin{cases} s_1 & \text{if x_{j_1}=c_1} \\ s_2 & \text{if x_{j_1} \neq c_1 and x_{j_2}=c_2} \\ s_3 & \text{if x_{j_1} \neq c_1 and x_{j_2} \neq c_2 and x_{j_3}=c_3} \\ \vdots & \\ s_{k_i} & \text{if x_{j_1} \neq c_1 and}\,\,\,\dots\,\,\,\text{and x_{j_{k_i}}=c_{k_i}} \\ s_d & \text{otherwise} \end{cases} \label{eq:10}$$ (10) for |${\boldsymbol 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 behaviour 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. Figure 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$|⁠. Fig. 1. View largeDownload slide (Color online) Time evolution of the average sensitivities for (a) BFs, (b) TFs, (c) HBFs and (d) NCFs. Insets are enlarged views from epoch |$20000$| to |$30000$| for the first three update rules and that from |$50000$| to |$60000$| for NCFs. BNs with NCFs were simulated for a longer period because their convergence is slower than the others. Fig. 1. View largeDownload slide (Color online) Time evolution of the average sensitivities for (a) BFs, (b) TFs, (c) HBFs and (d) NCFs. Insets are enlarged views from epoch |$20000$| to |$30000$| for the first three update rules and that from |$50000$| to |$60000$| for NCFs. BNs with NCFs were simulated for a longer period because their convergence is slower than the others. 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$|⁠. Fig. 2. View largeDownload slide (Color online) Comparison between numerical stationary in-degree distributions (symbols) and theoretical stationary in-degree distributions (lines) for (a) BFs, (b) TFs, (c) HBFs and (d) NCFs. Fig. 2. View largeDownload slide (Color online) Comparison between numerical stationary in-degree distributions (symbols) and theoretical stationary in-degree distributions (lines) for (a) BFs, (b) TFs, (c) HBFs and (d) NCFs. 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. Fig. 3. View largeDownload slide (Color online) 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$| for (a) |$\sigma=0.95$|⁠, (b) |$\sigma=1.00$| and (c) |$\sigma=1.05$|⁠. Dotted lines have a slope equal to |$\sigma$|⁠. Insets are enlarged views of the region |$0 \leq d_t, d_{t+1} \leq 0.1$|⁠. Fig. 3. View largeDownload slide (Color online) 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$| for (a) |$\sigma=0.95$|⁠, (b) |$\sigma=1.00$| and (c) |$\sigma=1.05$|⁠. Dotted lines have a slope equal to |$\sigma$|⁠. Insets are enlarged views of the region |$0 \leq d_t, d_{t+1} \leq 0.1$|⁠. 5. Discussion In this article, we proposed a new activity-dependent adaptive BN 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 the normalization condition for the mean-field solution holds. In particular, the condition is satisfied when TAS is smaller than |$\lambda_k$| in the limit of large |$k$| and |$b^*$| is independent of |$P_{\text{in}}^s$|⁠. We 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 BNs, 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 modelled several real-world GRNs including the bacterium GRN by BFs 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 four inputs do not affect the near critical dynamical behaviour of the bacterium GRN. On the other hand, most of genes in the bacterium GRN have at most three 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 article. 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. Acknowledgements The author thanks the anonymous reviewers for their helpful comments to improve the manuscript. Funding JSPS KAKENHI Grant Number 25280091, in part. 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 Section 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 Section 4 except that we only show results for |$k_0=3$| here. In Fig. A1, 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. A2) and the Derrida plot (Fig. A3). On the other hand, theoretical predictions and results of numerical simulation for the other cases agree well in both Figs A2 and A3. These results provide further support for the claim at the end of Section 3.1 on the applicability and the limitation of our mean-field theory. Fig. A1. View largeDownload slide (Color online) Time evolution of the average sensitivity. (a) BFs. The values of TAS |$\sigma$| are |$0.5$|⁠, |$1.0$|⁠, |$1.5$|⁠, |$2.0$| and |$2.5$| from below. (b) HBFs. The values of TAS |$\sigma$| are |$0.5$|⁠, |$1.0$|⁠, |$1.5$| and |$2.0$| from below. Fig. A1. View largeDownload slide (Color online) Time evolution of the average sensitivity. (a) BFs. The values of TAS |$\sigma$| are |$0.5$|⁠, |$1.0$|⁠, |$1.5$|⁠, |$2.0$| and |$2.5$| from below. (b) HBFs. The values of TAS |$\sigma$| are |$0.5$|⁠, |$1.0$|⁠, |$1.5$| and |$2.0$| from below. Fig. A2. View largeDownload slide (Color online) Comparison between numerical stationary in-degree distributions (symbols) and theoretical stationary in-degree distributions (lines) for (a) BFs and (b) HBFs. The theoretical stationary in-degree distribution does not exist for HBFs with |$\sigma=2.0$|⁠. The shown theoretical line is calculated by Eq. (9) and is truncated by |$N=200$|⁠. Fig. A2. View largeDownload slide (Color online) Comparison between numerical stationary in-degree distributions (symbols) and theoretical stationary in-degree distributions (lines) for (a) BFs and (b) HBFs. The theoretical stationary in-degree distribution does not exist for HBFs with |$\sigma=2.0$|⁠. The shown theoretical line is calculated by Eq. (9) and is truncated by |$N=200$|⁠. Fig. A3. View largeDownload slide (Color online) 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$| for (a) BFs and (b) HBFs. Dotted lines have the same slope as corresponding value of TAS |$\sigma$|⁠. Insets are enlarged views of the region |$0 \leq d_t, d_{t+1} \leq 0.1$|⁠. Fig. A3. View largeDownload slide (Color online) 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$| for (a) BFs and (b) HBFs. Dotted lines have the same slope as corresponding value of TAS |$\sigma$|⁠. Insets are enlarged views of the region |$0 \leq d_t, d_{t+1} \leq 0.1$|⁠. References 1. Drossel B. ( 2008 ) Random Boolean networks. Reviews of Nonlinear Dynamics and Complexity ( Schuster H. G. ed.). Weinheim, Germany : Wiley-VCH Verlag GmbH & Co. KGaA . 2. Kauffman S. A. ( 1969 ) Metablic stability and epigenesis in randomly constructed genetic nets. J. Theor. Biol. , 22 , 437 – 467 . Google Scholar Crossref Search ADS PubMed 3. Kürten K. E. ( 1988 ) Critical phenomena in model neural networks. Phys. Lett. A , 129 , 157 – 160 . Google Scholar Crossref Search ADS 4. Paczuski M. , Bassler K. E. & Corral Á. ( 2000 ) Self-organized networks of competing Boolean agents. Phys. Rev. Lett. , 84 , 3185 – 3188 . Google Scholar Crossref Search ADS PubMed 5. Derrida B. & Pomeau Y. ( 1986 ) Random networks of automata: a simple annealed approximation. Europhys. Lett. , 1 , 45 – 49 . Google Scholar Crossref Search ADS 6. Beggs J. M. & Plenz D. ( 2003 ) Neuronal avalanches in neocortical circuits. J. Neurosci. , 23 , 11167 – 11177 . Google Scholar Crossref Search ADS PubMed 7. Petermann T. , Thiagarajan T. C. , Lebedev M. A. , Nicolelis M. A. L. , Chialvo D. R. & Plenz D. ( 2009 ) Spontaneous cortical activity in awake monkeys composed of neuronal avalanches. Proc. Natl. Acad. Sci. USA , 106 , 15921 – 15926 . Google Scholar Crossref Search ADS 8. Balleza E. , Alvarez-Buylla E. R. , Chaos A. , Kauffman S. , Shmulevich I. & Aldana M. ( 2008 ) Critical dynamics in genetic regulatory networks: examples from four kingdoms. PLoS One , 3 , e2456 . Google Scholar Crossref Search ADS PubMed 9. Nykter M. , Price N. D. , Aldana M. , Ramsey S. A. , Kauffman S. A. , Hood L. E. , Yli-Harja O. & Shmulevich I. ( 2008 ) Gene expression dynamics in the macrophage exhibit criticality. Proc. Natl. Acad. Sci. USA , 105 , 1897 – 1900 . Google Scholar Crossref Search ADS 10. Valverde S. , Ohse S. , Turalska M. , Garcia-Ojalvo J. & West B. J. ( 2015 ) Structural determinants of criticality in biological networks. Front. Physiol. , 6 , 127 . Google Scholar Crossref Search ADS PubMed 11. Bertschinger N. & Natschläger T. ( 2004 ) Real-time computation at the edge of chaos in recurrent neural networks. Neural Comput. , 16 , 1413 – 1436 . Google Scholar Crossref Search ADS PubMed 12. Goudarzi A. , Teuscher C. , Gulbahce N. & Rohlf T. ( 2012 ) Emergent criticality through adaptive information processing in Boolean networks. Phys. Rev. Lett. , 108 , 128702 . Google Scholar Crossref Search ADS PubMed 13. Kinouchi O. & Copelli A. M. ( 2006 ) Optimal dynamical range of excitable networks at criticality. Nat. Phys. , 2 , 348 – 352 . Google Scholar Crossref Search ADS 14. Haldeman C. & Beggs J. M. ( 2005 ) Critical branching captures activity in living neural networks and maximizes the number of metastable states. Phys. Rev. Lett. , 94 , 058101 . Google Scholar Crossref Search ADS PubMed 15. Bornholdt S. & Röhl T. ( 2003 ) Self-organized critical neural networks. Phys. Rev. E , 67 , 066118 . Google Scholar Crossref Search ADS 16. Rybarsch M. & Bornholdt S. ( 2014 ) Avalanches in self-organized critical neural networks: a minimal model for the neural SOC universality class. PLoS One , 9 , e93090 . Google Scholar Crossref Search ADS PubMed 17. Meisel C. & Gross T. ( 2009 ) Adaptive self-organization in a realistic neural network model. Phys. Rev. E , 80 , 061917 . Google Scholar Crossref Search ADS 18. Rubinov M. , Sporns O. , Thivierge J.-P. & Breakspear M. ( 2011 ) Neurobiologically realistic determinants of self-organized criticality in networks of spiking neurons. PLoS Comput. Biol. , 7 , e1002038 . Google Scholar Crossref Search ADS PubMed 19. Levina A. , Herrmann J. M. & Geisel T. ( 2007 ) Dynamical synapses causing self-organized criticality in neural networks. Nat. Phys. , 3 , 857 – 860 . Google Scholar Crossref Search ADS 20. Levina A. , Herrmann J. M. & Geisel T. ( 2009 ) Phase transitions towards criticality in a neural system with adaptive interactions. Phys. Rev. Lett. , 102 , 118110 . Google Scholar Crossref Search ADS PubMed 21. Droste F. , Do A.-L. & Gross T. ( 2013 ) Analytical investigation of self-organized criticality in neural networks. J. R. Soc. Interface , 10 , 20120558 . Google Scholar Crossref Search ADS PubMed 22. MacArthur B. D. , Sánchez-García R. J. & Ma’ayan A. ( 2010 ) Microdynamics and criticality of adaptive regulatory networks. Phys. Rev. Lett. , 104 , 168701 . Google Scholar Crossref Search ADS PubMed 23. Lee D.-S. ( 2014 ) Evolution of regulatory networks towards adaptability and stability in a changing environment. Phys. Rev. E , 90 , 052822 . Google Scholar Crossref Search ADS 24. Gross T. & Sayama H. (eds) ( 2009 ) Adaptive Networks: Theory, Models and Applications . Heidelberg : Springer . 25. Sayama H. , Pestov I. , Schmidt J. , Bush B. J. , Wong C. , Yamanoi J. & Gross T. ( 2013 ) Modeling complex systems with adaptive networks. Comput. Math. Appl. , 65 , 1645 – 1664 . Google Scholar Crossref Search ADS 26. Bornholdt S. & Rohlf T. ( 2000 ) Topological evolution of dynamical networks: global criticality from local dynamics. Phys. Rev. Lett. , 84 , 6114 – 6117 . Google Scholar Crossref Search ADS PubMed 27. Christensen K. , Donangelo R. , Koiller B. & Sneppen K. ( 1998 ) Evolution of random networks. Phys. Rev. Lett. , 81 , 2380 – 2383 . Google Scholar Crossref Search ADS 28. Liu M. & Bassler K. E. ( 2006 ) Emergent criticality from coevolution in random Boolean networks. Phys. Rev. E , 74 , 041910 . Google Scholar Crossref Search ADS 29. Górski P. J. , Czaplicka A. & Hołyst A. ( 2016 ) Coevolution of information processing and topology in hierarchical adaptive random Boolean networks. Eur. Phys. J. B , 89 , 33 . Google Scholar Crossref Search ADS 30. Rohlf T. ( 2008 ) Self-organization of heterogeneous topology and symmetry breaking in networks with adaptive thresholds and rewiring. Europhys. Lett. , 84 , 10004 . Google Scholar Crossref Search ADS 31. Peter I. S. & Davidson E. H. ( 2011 ) Evolution of gene regulatory networks controlling body plan development. Cell , 144 , 970 – 985 . Google Scholar Crossref Search ADS PubMed 32. Wittkopp P. J. & Kalay G. ( 2012 ) Cis-regulatory elements: molecular mechanisms and evolutionary processes underlying divergence. Nat. Rev. Genet. , 13 , 59 – 69 . Google Scholar Crossref Search ADS 33. Rohlf T. & Bornholdt S. ( 2009 ) Self-organized criticality and adaptation in discrete dynamical networks. In Adaptive Networks ( Gross T. & Sayama H. eds). Heidelberg, Germany : Springer , pp. 73 – 106 . 34. Newman M. E. J. , Strogatz S. H. & Watts D. J. ( 2001 ) Random graphs with arbitrary degree distributions and their applications. Phys. Rev. E , 64 , 026118 . Google Scholar Crossref Search ADS 35. Lee D.-S. & Rieger H. ( 2007 ) Comparative study of the transcriptional regulatory networks of E. coli and yeast: structural characteristics leading to marginal dynamic stability. J. Theor. Biol. , 248 , 618 – 626 . Google Scholar Crossref Search ADS PubMed 36. Shmulevich I. & Kauffman S. A. ( 2004 ) Activities and sensitivities in Boolean network models. Phys. Rev. Lett. , 93 , 048701 . Google Scholar Crossref Search ADS PubMed 37. Squires S. , Ott E. & Girvan M. ( 2012 ) Dynamical instability in Boolean networks as a percolation problem. Phys. Rev. Lett. , 109 , 085701 . Google Scholar Crossref Search ADS PubMed 38. Rohlf T. & Bornholdt S. ( 2002 ) Criticality in random threshold networks: annealed approximation and beyond. Physica A , 310 , 245 – 259 . Google Scholar Crossref Search ADS 39. Pomerance A. , Ott E. , Girvan M. & Losert W. ( 2009 ) The effect of network topology on the stability of discrete state models of genetic control. Proc. Natl. Acad. Sci. USA , 106 , 8209 – 8214 . Google Scholar Crossref Search ADS 40. Kauffman S. , Peterson C. , Samuelsson B. & Troein C. ( 2003 ) Random Boolean network models and the yeast transcriptional network. Proc. Natl. Acad. Sci. USA , 100 , 14796 – 14799 . Google Scholar Crossref Search ADS 41. Peixoto T. P. ( 2010 ) The phase diagram of random Boolean networks with nested canalizing functions. Eur. Phys. J. B , 78 , 187 – 192 . Google Scholar Crossref Search ADS 42. Derrida B. & Weisbuch G. ( 1986 ) Evolution of overlaps between configurations in random Boolean networks. J. Phys. , 47 , 1297 – 1303 . Google Scholar Crossref Search ADS 43. Kesseli J. , Rämö P. & Yli-Harja O. ( 2006 ) Iterated maps for annealed Boolean networks. Phys. Rev. E , 74 , 046104 . Google Scholar Crossref Search ADS 44. Hesse J. & Gross T. ( 2014 ) Self-organized criticality as a fundamental property of neural systems. Front. Syst. Neurosci. , 8 , 1 – 14 . Google Scholar Crossref Search ADS PubMed 45. Marković D. & Gros C. ( 2014 ) Power laws and self-organized criticality in theory and nature. Phys. Rep. , 536 , 41 – 74 . Google Scholar Crossref Search ADS 46. Bak P. , Tang C. & Wiesenfeld K. ( 1987 ) Self-organized criticality: an explanation of |$1/f$| noise. Phys. Rev. Lett. , 59 , 381 – 384 . Google Scholar Crossref Search ADS PubMed 47. Jensen H. J. ( 1998 ) Self-organized Criticality: Emergent Complex Behavior in Physical and Biological Systems . Cambridge : Cambridge University Press . 48. Aldana M. , Balleza E. , Kauffman S. & Resendiz O. ( 2007 ) Robustness and evolvability in genetic regulatory networks. J. Theor. Biol. , 245 , 433 – 448 . Google Scholar Crossref Search ADS PubMed 49. Harris S. E. , Sawhill B. K. , Wuensche A. & Kauffman S. ( 2002 ) A model of transcriptional regulatory networks based on biases in the observed regulation rules. Complexity , 7 , 23 – 40 . Google Scholar Crossref Search ADS 50. Haruna T. & Tanaka S. ( 2014 ) On the relationship between local rewiring rules and stationary out-degree distributions in adaptive random Boolean network models. In Artificial Life 14: Proceedings of the Fourteenth International Conference on the Synthesis and Simulation of Living Systems ( Sayama H. Rieffel J. Risi S. Doursat R. & Lipson H. eds). Cambridge : MIT Press , pp. 419 – 426 . 51. Squires S. , Pomerance A. , Girvan M. & Ott E. ( 2014 ) Stability of Boolean networks: the joint effects of topology and update rules. Phys. Rev. E , 90 , 022814 . Google Scholar Crossref Search ADS © The author 2018. Published by Oxford University Press. All rights reserved. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)

### Journal

Journal of Complex NetworksOxford University Press

Published: Dec 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations