A Definition of the Truth Table Draws and Constraints on \{p^{(k,s)}_{i}\}

The Stability of Boolean Networks with Generalized Canalizing Rules


Boolean networks are discrete dynamical systems in which the state (zero or one) of each node is updated at each time to a state determined by the states at time of those nodes that have links to it. When these systems are used to model genetic control, the case of ‘canalizing’ update rules is of particular interest. A canalizing rule is one for which a node state at time is determined by the state at time of a single one of its inputs when that inputting node is in its canalizing state. Previous work on the order/disorder transition in Boolean networks considered complex, non-random network topology. In the current paper we extend this previous work to account for canalizing behavior.

complex networks, genetic networks, Boolean networks

I Introduction

Boolean networks have been extensively studied as models for genetic control of cells (1); (2). In this framework, the genetic regulatory network is modeled as a directed graph, where links correspond to the influence of one gene on the expression of another. Individual genes are either off or on, represented as 0 or 1, respectively, and the state of a gene at time is given by a Boolean update function of the states of its inputs at time . In early analyses, both the network topology and the update functions were randomly assigned. In particular, Kauffman’s network model (1); (3) has received significant study. According to this model, there are nodes (genes) in the network, each having the same number of input links, , and the nodes from which these input links originate are chosen randomly with uniform probability. Additionally, the update function determining the time evolution at each node is defined by a random, time-independent, -entry truth table that can be characterized (4) by the ‘bias’ , which, as discussed subsequently, is the probability that a one appears as the output of the update function. Using the Hamming distance between two states of the system (i.e., the number of nodes for which the states disagree) as the distance measure, these systems, when large, exhibit both disordered (unstable) behavior, where the distance between typical initially close states on average grows in time, as well as ordered (stable) behavior, where the distance decreases. Separating the two is a ‘critical’ regime.

So-called ‘canalizing functions’ are a significant modification of the random truth table model of previous work. Canalizing functions, believed to be relevant to genetic networks (5); (6), are those functions where an argument of the function (the ‘canalizing input’), having a certain value (the ‘canalizing value’), determines the value of the function independently of the values of the other arguments of the function (inputs) (2). If the canalizing input does not have the canalizing value, the function is determined by the other inputs. (Further refinements can include a hierarchy of canalization (7).) Previous work (7); (10); (8); (9) has considered random network topology and random choices of canalizing inputs and, with those assumptions, has demonstrated that canalizing functions often stabilize networks that would be disordered in their absence. A key quantity, defined in previous work by Shmulevich and Kauffman, is the ‘activity’ (11); (12) of a Boolean variable on a Boolean function. This quantity can be used to characterize the increased importance of canalizing inputs and plays a crucial role in the theory presented in this paper.

An important hypothesis in genetic networks is that biological systems exist at the transition between order and disorder (1) (i.e., the ‘life at the edge of chaos’ hypothesis). It has also been suggested that this transition may be relevant to the onset of cancer (13). This paper is concerned with deriving a condition determining the ‘edge of chaos’ in Boolean networks taking into account the specific topology of the network (as was treated in Ref. (13)) and canalizing update rules. In Sec. II, we present definitions that will be used in the remainder of the paper and summarize previous results on Boolean network stability. In Sec. III, we present a generalized probabilistic model of canalizing behavior, and we use this model with the Shmulevich-Kauffman activity to extend the results of Ref. (13) to the case of networks with canalizing functions. We derive a stability criterion that offers an advancement over previous work on Boolean network models with canalizing rules in two ways: (1) it accounts for the specific topological properties of a considered Boolean networks, rather than purely random topology; and (2) it accounts for possible correlation between network topology and canalizing behavior, e.g., correlation between the choice of canalizing input and local topology. In Sec. IV, we numerically test our derived stability criterion on a variety of complex network topologies, and in Sec. V we explore the effect of correlation between strict canalization and network topology.

Ii Preliminary Definitions and Previous Results

A Boolean network is defined by a state vector , where each , and a set of update functions , such that


where denote the indices of the nodes that input to node , and we denote this set of nodes by . [In the following discussion, , which is between 1 and , is used to label an input to node , or, equivalently an argument of ; , which is between 1 and , refers to the network index of the node corresponding to input ; and are used to switch between them. Similarly, is the state of node , and is the -th input to .] The number of input links to node is called its in-degree, and the number of output links from node is called its out-degree. The update function at each node is usually defined by a truth table, where the table for node has rows, one for each possible set of the states of the nodes that input to node , and each input state row is followed by its resulting update output state for node , thus forming a entry output column. Thus the input entries of the rows of the truth table represent all possible arguments of in Eq. (1), and the output column gives the value that assumes for that set of its arguments.

An important property of a given Boolean network is whether it is ‘ordered’ or ‘disordered.’ This property is defined in large networks by considering the trajectories resulting from two close initial states, and . To quantify their divergence, the Hamming distance of coding theory is used: , and the two initial conditions are ‘close’ if . If the network is ordered, on average as . In disordered networks, quickly increases to , while a ‘critical’ network is at the border separating the two regimes. (In this paper, we use the word ‘stable’ to refer to ordered networks.)

To obtain an analytically tractable model for the stability of Boolean networks, Derrida and Pomeau (4) considered an ‘annealing’ procedure applied to the original problem and calculated the probability that, after steps, a node state is the same on two trajectories that originated from initially close conditions. Annealing, in this context, refers to assuming that both the at each node and the network of connections are randomly reassigned at each time step, in contrast to the ‘quenched’ problem, in which both the update functions and topology are constant for all time. Later authors generalized the Derrida-Pomeau analysis to include variable in-degree (14); (15); (16) and joint in-degree/out-degree distributions (17). In the annealed situation, at each time step both the truth table outputs and the network of connections are randomly chosen. The actual situation of interest, however, is the case of ‘frozen-in’ networks, where the truth table and network of connections are fixed in time. It was hypothesized and subsequently numerically confirmed that, for large networks, results obtained in the analytically tractable annealed situation are the same as those in the analytically intractable frozen situation; this is referred to as the ‘annealing approximation.’

In Ref. (13), we developed a more general approximate technique for determining the stability of large Boolean networks. This technique used what we called a ‘semiannealing’ procedure, which differed from the situation considered by Derrida in that only the truth table outputs are randomly reassigned at each time step; the network of connections was kept frozen. By using this procedure, one can investigate the effect of given specific network topology on the stability of a Boolean network. Numerical experiments confirmed the predictions of the semiannealed theory (13) in cases of complex network topology, exploring such effects as correlation between the number of inputs and outputs at each node (18), assortativity (19), community structure (20), and the presence of small motifs similar to those found in biological systems (21). These cases cannot be accounted for using the full annealing approximation.

When considering annealed Boolean truth tables, the update function of each node, , is solely parameterized by the ‘bias’ of node , denoted . The bias is the probability that any given output entry in the truth table is a one. This can be recast in terms of the ‘sensitivity’ of node , denoted , which is the probability that differing inputs to yield a differing output of . In the annealed truth table case, .

However, this parameterization is not complete for canalizing functions. In Ref. (8), Kaufman et al. consider the case of networks with and non-annealed truth tables. For networks, the possible Boolean functions fall into four classes: frozen functions (i.e., those functions that output the same value, one or zero, independent of their inputs), non-canalizing functions, and two types of canalizing functions. Reference (8) shows that, when there are fewer nodes with non-canalizing functions than with frozen functions (by the addition of either frozen functions or canalizing functions), the network is in the ordered regime. Another approach to analyzing stability of Boolean networks with canalizing functions is to consider ‘damage spreading’ (i.e., the propagation of disagreement between and ). This approach was used in Ref. (9) for Boolean networks with canalizing rules, with results similar to (8).

In Sec. III, we refine the semiannealed procedure of Ref. (13), which relies heavily on the sensitivity, to account for canalizing behavior. However, as discussed above, the sensitivity does not include sufficient information to account for canalizing behavior. This quantity treats all input nodes equally, which is clearly not appropriate for canalizing functions. Instead, we make use of the ‘activity’ of node on node , (or, equivalently, the activity of input argument on update function , ). First described by Shmulevich and Kaufmann in Ref. (11), the ‘activity’ of input on a Boolean function is the probability that a bit flip of only input leads to a flip of the output of if all inputs other than are independent and have equal probabilities of being one or zero. This definition clearly emphasizes canalizing inputs over non-canalizing inputs (11). In our presentation of generalized canalization below, we calculate the activities for truth tables with this property and use the probabilistic interpretation of our definition to define a new semiannealing procedure, which has the activities as its defining parameters.

Iii Generalized Canalization and the Stability Criterion

The semiannealing procedure used in Ref. (13) independently and randomly reassigned the output elements of the truth table governing node to be one or zero with probability or , respectively. This is not true of canalizing functions: if the canalizing input takes its canalizing value in a given row of the table, the probability of a one appearing in the corresponding entry of the output column is zero (or one). We call this behavior ‘strictly canalizing.’ We now introduce a generalization of strictly canalizing behavior to the case of ‘quasicanalizing’ inputs, which we define as the case where the probability that a one appears in the output of node ’s truth table if input takes value , averaged over all other inputs, depends on . We illustrate this concept in Tables 1(a)-(c). (The average of , weighting the states and equally for all , is , which we call the ‘truth table bias.’ This quantity will be discussed in detail below.) Strict canalization with respect to input , illustrated in Table 1(a), is the case when or when is the canalizing value. If , is a non-canalizing input to , and Table 1(b) consists entirely of non-canalizing inputs. Quasicanalizing inputs are neither strictly canalizing nor non-canalizing: if and for , then input is said to be quasicanalizing. Table 1(c) illustrates this property. We call truth tables where all inputs are non-canalizing ‘unstructured’ (e.g., Table 1(b)), and those with any canalizing inputs, strict or quasi-, ‘structured’ (e.g., Tables 1(a) and (c)). We note that a complete, realizable set of is highly constrained, and the constraints are given in the Appendix.

10 0
11 0
Table 1: Tables illustrating three types of canalizing behavior in a function with two inputs. In (a), the first input is strictly canalizing: when it is one, the output of the function is guaranteed to be zero; if it is zero, the probability that the function takes the value one is . In (b), neither input is canalizing and, for each input, the output of the function is independently one with probability . In (c), the first input is quasicanalizing: if it one, the function is one with probability , otherwise it is one with probability . In all three cases, the second input is noncanalizing.

We now write down the procedure by which we can make random draws of the truth tables conforming to the set of that defines the update function . The details of the calculation are given in the Appendix, but we state here that, for a given draw the probability that is one, , is given by


This defines an ensemble of update functions for each node , where the sampling probability of each member can be derived from (2). Loosely speaking, we would like to be able to average the dynamics of all networks over this ensemble to derive the stability criterion (i.e., the boundary between the stable and unstable phases in the space of all ). However, this is analytically intractable since the dynamics of a given network depend crucially on the exact nodal update functions (e.g., frozen functions which are sampled with non-zero probability), and the vast number of sets of makes member-by-member calculation infeasible. Instead, we use a semiannealing sampling procedure to approximate this average by considering a system where, at each time step, a different randomly chosen member of the ensemble of is used to advance the system. Thus, at each time step, we randomly assign the truth table output of , when presented with inputs to be one with probability and zero with probability . The motivation for this semiannealed approach, as well as for other previous annealing approaches, is (i) that they are analytically tractable, and (ii) that it is supposed that, for large networks, the stability in the annealed situation approximates that of a typical (with respect to the ensemble) frozen-in case. Here by the frozen-in case we refer to a system with a given temporally constant network topology and set of truth tables, which are randomly assigned initially using the probabilities in the annealing procedure. We numerically confirm that the stability of frozen-in networks corresponds well with that of our semiannealed networks (Sec. IV).

We now follow Refs. (11); (12) and introduce a quantity that we call the ‘activity’, which we define as the probability that, at any time on an orbit, the output of changes if input is flipped (either from 0 to 1 or vice versa), keeping the other inputs unchanged. The activity will play an important role in the derivation of the stability criterion. In order to obtain the activities, we first make the assumption that orbits on the semiannealed network are ergodic, and, given such an orbit, we can define as the fraction of time that the state ; we call this the ‘dynamic bias’ since it is determined by the dynamics of the network, in contrast to the truth table bias above. Assuming independence of the probability of nodal input states to node (appropriate to our supposition of locally-tree like topology, discussed below), we have


where the sum is over all possible inputs to node , which is denoted ; is the probability that the set of states of the input nodes to node , , takes the values in the set (where again denotes the argument number to the update function , and refers to a node index in the network); and is a shorthand notation for the quantity in Eq. (2). Note that is a purely local quantity; i.e., it is determined from the random truth table annealing process for node and is independent of the truth table assignments at other nodes. The probability of the inputs to node taking the values can in turn be written in terms of the dynamic biases of the input nodes,


where and is used to convert from node index to argument number . In principle, one could insert Eq. (4) into Eq. (3) and iteratively solve for all (provided the iteration converged) or, alternatively, one could numerically generate a long orbit and approximate as the fraction of time that on the orbit. In this paper, as an example, we will restrict our consideration to a case where a very simple exact solution to Eqs. (3) and (4) is available, and we will base our numerical experiments in Sec. IV on that case.

We now define to be a input function that denotes the probability that the truth table output is one if input is given some element set of other inputs; again, we emphasize that the are local quantities. (The input to this function is a reduced set of inputs , which is related to by removing the -th element; we suppress the explicit dependence below.) With this, we calculate the activity as


where is the probability of node having the reduced input set . In the numerical tests in Sec. IV, we consider this probability to be uniform over all possibilities; i.e., we consider the trivial solution to Eqs. (3) and (4) given by and . This implies that both canalized and canalizing values are 0 or 1 with probability one-half and, further, that truth table biases are symmetrically distributed around one-half. However, in cases where there is a bias in the canalization, e.g., a node is more likely to be canalized to one, or truth table biases are not symmetrically distributed, is no longer uniform, and a procedure to calculate the for each node as described above may be employed.

Our semiannealing procedure defines an ensemble of nodal truth table time courses, and we can define probabilities of the dynamically evolving states for an arbitrary member of this ensemble. In particular, as in Ref. (13), we define the -dimensional vector , where each element tracks the probability that node differs between two initally close states after time-steps: . Our goal is to find an update equation for and perform linear stability analysis on the solution to obtain the stability criterion. The update equation will be derived under the assumption that the inputs are statistically independent of one another. This assumption holds in the case of locally tree-like topology (13); (22).

Since we are performing linear stability analysis, we can make several simplifying approximations. The probability of inputs to node being different between the trajectories of two initially close states is of order . Linear stability applies for small, so the case of multiple inputs to are flipped can be ignored, and we can approximate the probability that only input node to node is flipped as . The probability that the single bit flip of node causes a flip in the output of is the activity . Since the input flipping and the the probability that this leads to a flip in the output of are independent, the probability that a flip in node occurs and leads to a flip in the node in the next time step is thus . In the linear stability limit, we can sum these probabilities inputs to get the following approximate evolution equation for small perturbations from the solution :


This can be written in matrix form after discarding the higher-order terms as


where is the ‘activity matrix’ with elements if there is a link from to (), and zero otherwise. From this equation, we see that stability is determined by the largest eigenvalue of this matrix:


Note that, since the elements of the matrix are all non-negative, the Perron-Frobenius theorem guarantees that is real and non-negative.

Before concluding this section, we note the relationship of (8) to the case of unstructured truth tables (13). When there are no canalizing inputs present, the right hand side of Eq. (2) reduces to for all possible inputs, and is constant across each row of the matrix. This is the central result of Ref. (13).

Iv Numerical Results

In this section, we present numerical results testing our derived criterion for the stability of Boolean networks with canalizing truth tables. We test the theory by measuring the long-time Hamming distance between trajectories that differ in only a few initial states as a function of . We vary in two ways: (a) a varying proportion of nodes have a single, strictly canalizing input; and (b) each node has a single quasicanalizing input, and that input has varying activity on each node’s function. In Sec. IV.A, we treat case (a) with several network degree distributions. We consider networks, which is most comparable to previous work on canalizing inputs (8); (9), as well as networks with power-law in-degree distribution. We treated the cases of assortativity/disassortativity and community structure in Ref. (13), and we believe these cases can be treated in an analogous way using the activity matrix. In Sec. IV.B, we consider case (b) with networks of the same general topologies. The networks under consideration in Sec. IV.B, however, have more average inputs as compared to Sec. IV.A.

Our general approach to test the transition from order to disorder is the following. We begin with a given network topology and, using the techniques of Ref. (13), calculate a uniform sensitivity for all nodes that would yield a Boolean network slightly in the disordered regime. We then choose a uniform truth table bias for all nodes that satisfies Thus, for all networks under consideration in this section, Ref. (13) would predict that the network is slightly disordered. We then vary by methods (a) or (b) above. As discussed in Sec. V, the addition of canalizing behavior tends to stabilize a network, and thus the transition is approached from the disordered regime. We choose the canalizing input of each update function randomly.

Figures 1 and 2 demonstrate the results of our numerical tests. Each data point in Figs. 1 and 2 is the average steady-state Hamming distance measurement of 1000 different frozen realizations of the truth tables. Each network has nodes, and the steady-state Hamming distance is calculated as the average Hamming distance from time to between trajectories that have an initial Hamming distance of 10 (0.01% of the nodes are flipped). The important result of these figures is that the critical stability condition derived using our truth table annealing procedure agrees well with the numerical results from our simulations of frozen-in (quenched) systems.

iv.1 Networks with Strictly Canalizing Inputs

Figure 1: Steady state Hamming distance vs. for (a) networks with (squares), 4 (diamonds), or 5 (triangles) inputs; and (b) networks with power-law in-degree distibution with average (squares), 4 (diamonds), or 5 (triangles). Each datapoint consists of 1000 realizations of the truth table update functions, where the initial Hamming distance between close states is 10 nodes (0.01%). is varied by giving an increasing proportion of nodes a single, strictly canalizing input. The predicted transition is at .

Figure 1 is a numerical test of the theory of Sec. III on networks with an increasing proportion of nodes with a single strictly canalizing input, where the proportion increases from zero to one. For each node in the network, we choose whether the node will have a canalizing input with probability . When , takes its maximum value; when , takes its minimum value.

In order to calculate the activity with a uniform ensemble of inputs in Eq. (5), some care must be exercised when constructing the truth tables for each node, and the procedure is as follows. The canalized output for each truth table is chosen to be zero or one with equal probability, and the remaining values in the table are assigned one with probability or , depending on whether the canalized output is zero or one, in order to maintain a constant truth table bias over all nodes. That the canalized output is zero or one with equal probability is crucial, since we are assuming that any given set of input values to a node are equally likely in our evaluation of Eq. (5). A preponderance of canalized zeros or ones violates this assumption (although the presented calculation of the activity can be refined to take this into account as discussed in Sec. III). Similar care must be taken in the construction of unstructured truth tables. Since there are two solutions to for a given , the bias for each unstructured truth table is chosen to be either of the solutions with equal probability.

Figure 1(a) tests the theory on networks, where and 5. The networks are constructed by randomly drawing inputs to each node with uniform probability from the other nodes in the network, subject to the constraint that each node has exactly outputs. The sensitivities for these networks are and 0.21 for and 5 networks, respectively. According to the theory in Ref. (13), these parameters are in the disordered regime. We see that in all cases, the networks appear to undergo a transition from ordered to disordered near as we add more canalizing inputs. Additionally, we see that the scaling of the steady-state Hamming distance with is a strong function of , the number of inputs. This is to be contrasted with the case of power-law degree distribution discussed below.

Figure 1(b) considers networks where the in- and out-degrees are independently drawn from truncated power-law degree distributions: if , and 0 otherwise. Networks are then constructed by randomly making connections between nodes in accord with their degrees using the configuration model. The figure depicts cases where the average number of inputs are the same as the case considered in Fig. 1(a): for , , and ; for , , and ; and for , , and . The effective sensitivities are the same as the cases which again place them slightly in the disordered regime. We see again that there appears to be a transition to ordered dynamics as canalizing inputs are added. However we note that the scaling of steady-state Hamming distance has much weaker dependence on than in the case. We also note that for the same the steady-state Hamming distance is smaller than the case.

iv.2 Quasicanalizing Inputs

Figure 2: Steady state Hamming distance vs. for (a) networks with 3 (squares), 4 (diamonds), or 5 (triangles) inputs; and (b) networks with power-law in-degree distibution with average (squares), 4 (diamonds), or 5 (triangles). Each datapoint consists of 1000 realizations of the truth table update functions, where the initial Hamming distance between close states is 10 nodes (0.01%). is varied by giving each node a single quasicanalizing input and increasing the difference , where is the index of the canalizing input. The predicted transition is at .

Figure 2 demonstrates the results of the second method of varying , where each node has a single quasicanalizing input. The tests were done on the same network topologies as in the previous section. When assigning a generalized canalized truth table, we randomly choose the canalizing value to be zero or one with uniform probability. To vary , we vary from zero to : when , all nodes have a single strictly canalizing input (i.e., the same as the case where above); when , is no longer a canalizing input and the network is identical to the case where above.

Figure 2(a) considers networks where all parameters are the same as in Fig. 1(a). We see the same behavior in steady-state Hamming distance as a function of as in Fig. 1(a). In Fig. 2(b), we use slightly different parameters than those in Fig. 1(b). The parameters used in the truncated power-law distribution are as follows: for , , and ; for , , and ; and for , , and . The initial sensitivities are and 0.146, respectively. These parameters were chosen in order to get larger truth tables at each node. We note that the steady-state Hamming distance is not much larger than in the cases of Fig. 1(a), and that the scaling is an even weaker function of the average number of inputs.

V Effect of Correlations Between Canalization and Topology

A useful feature of the theory of Sec. III is that it allows us to analyze the interplay between choice of canalizing input and local topology. Here we use a perturbation approach to derive an expression for the marginal change of as canalization is added to the network. We will see that this result depends on the local topology at the node under consideration and the ‘amount’ of canalization added, viz., the difference between biases when the canalizing input takes its canalizing value and when it does not. In order to test the results in this section, we will consider networks where each node has a single strictly canalizing input, and we adjust which input is the canalizing input to vary .

In the case of no canalizing inputs, the right hand side of Eq. (2) reduces to and for all inputs, which is the same as the unstructured case. In the case of exactly one canalizing input and the rest non-canalizing, we get from Eq. (2)


where is the index of the canalizing input of node .

We now investigate the effect of adding canalizing inputs to a network. We assume that we start with a network that has some nodes with no canalizing inputs and is characterized by an activity matrix with largest eigenvalue . We choose a node that is not canalized and change its truth table so that node has canalizing input , but the same , by perturbation analysis. Let and denote the activity matrix and the maximum eigenvalue corresponding to the altered system. Given left and right eigenvectors and of the matrix , perturbation theory gives the change in the eigenvalue as (23)


Since , by Eq. (9) and substituting , we see that


Using (11) in (10), we get


The first term in the product depends on the amount of canalization; the second, however, captures the effect of the network topology on . We note that, since the first term is always positive, the sign of the effect on the largest eigenvalue is determined by the second term. Usually the second term is negative, but of varying magnitude based on the choice of . In some cases, when the canalizing input has a large compared to the other inputs, canalization may increase instability.

Figure 3: The effect of correlation between canalization and local network topology. Steady state Hamming distance vs. for networks with power-law in-degree distibution with average (squares), 4 (diamonds), or 5 (triangles). Each datapoint consists of 1000 realizations of the truth table update functions, where the initial Hamming distance between close states is 10 nodes (0.01%). Each node in the network has a single, strictly canalizing input. In this figure, is varied by the choice of canalizing input: a varying proportion of nodes have their canalizing input chosen to maximize Eq. (12), the remainder have their canalizing input chosen to minimize Eq. (12). The predicted transition is at .

We test whether effectively predicts the order/disorder transition in Fig. 3 using networks constructed as in Sec. IV. All networks under consideration have nodes and truncated power-law degree distributions corresponding to the parameters used in Fig. 1(b): for , , and ; for , , and ; and for , , and . Each node in the network has a single canalizing input, and, to vary , a varying proportion of nodes have their canalizing input chosen to maximize Eq. (12), with the remaining nodes having their canalizing input chosen to minimize Eq. (12). That is, given a list of input nodes to a node , the canalizing input is chosen such that or , where is the eigenvector centrality of node . Since every node has a canalizing input, we use slightly different values for the initial sensitivities: and 0.222 for and 5, respectively. Once again, as in Fig. 1(b), we see that the transition is approximately at .

Vi Conclusion

In this paper we have generalized previous work on the stability of large Boolean networks to account for canalization. Our generalization allows a continuum in the degree of canalization and a probabilistic interpretation, as opposed to the previous canalization model (5), where an input could only be strictly canalizing or not canalizing at all. We define a semiannealing procedure which we used to derive the condition, Eq. (8), under which Boolean networks that have canalizing functions are stable. The stability criterion derived in this paper offers two advantages to the study of genetic networks in particular, because it successfully handles complex network topologies that may be found in real genetic control networks, and because it can account for possible correlations between canalizing behavior and network topology. Given the likely prominance of canalizing behavior in gene networks, these results may offer insight into the understanding of these systems. As an example of some of the insights than can be gleaned, Ref. (24) reports that a mutant strain of macrophage that has a gene silenced actually exhibits mildly chaotic behavior, in constrast to the unmutated case which exhibits criticality. This may indicate that the silenced gene is a canalizing input to a large number of genes; when the canalized gene is silenced, the states of the genes it is connected to are free to evolve according to the dynamics of their other inputs, thus reducing stability and yielding a slightly chaotic network. Furthermore, since our technique allows analysis of any specified network (e.g., an experimentally determined network), our stability criterion may, with advances in gene network measurement technology, allow one to assess the criticality of real genetic networks directly.

We thank Wolfgang Losert and the anonymous reviewers for their comments. This work was partially supported by ONR Grant N00014-07-1-0734.

Appendix A Definition of the Truth Table Draws and Constraints on

In order to define the appropriate annealing procedure used on the truth tables, we must specify the probability that a given output value of is one. We derive this probability in two ways: first using Bayes’ Theorem and again using counting arguments. The Bayesian approach has the advantage of simplicity and clarity, however the counting approach yields the entire set of constraints on realizable sets of . In both presentations, a useful quantity is the ‘truth table bias’ , which is the probability that any truth table output is one, similar to the unstructured case. Letting be the number of rows in the truth table, the expected number of ones in the output column of the truth table of node is . For any given arbitrary input to , rows in the truth table have and have . By definition of , the expected number of ones in the output of the truth table with entries that have is . The total expected number of ones is the sum of the expected number of ones when and when , which leads to


Note that, since the expected number of ones does not depend on our choice of above, must be independent of . This provides a constraint on the set of possible values that describe a realizable truth table. Non-canalizing inputs have both and equal to the effective bias by definition, and unstructured truth tables have .

We now derive the probability that a given set of input values to node , , yields an output of one, and we denote this probability . Using Bayes’ Theorem, we have


where is the event that the -th input takes the value (i.e., is the event that , where denotes a specifc value, 0 or 1, of the node ’s state variable ). By definition, . Since we are considering an ensemble where every possible input string to has equal probability, . We note that since each of the events are independent, and we calculate again using Bayes’ Theorem:


Using these results in Eq. (14), we obtain 2.

Another method to derive 2 involves a counting argument. While this is a little less clear, it yields the full set of constraints for a realizable set of . This argument proceds from Eq. (13). We now choose a second arbitrary input, and divide the truth table into four sections, each of size : those rows with ; those with ; those with ; and . By definition of and , we have the total expected number of ones in the truth table as


This yields a new set of constraints, namely that for each pair of inputs and , . This disallows conflicting rules, such as two strictly canalizing inputs that canalize to different values. We can continue in a similar manner by considering a third input and deriving an equation like Eq. (16) and obtain constraints on the of triplets of inputs and so on until we have divided the truth table up into L sections of one row each.


  1. S.A. Kauffman, J. Theor. Biol. 22, 437 (1969).
  2. S.A. Kauffman, The Origins of Order (Oxford University Press, New York, 1993).
  3. M. Aldana, S. Coppersmith, and L. P. Kadanoff in Perspectives and Problems in Nonlinear Science, edited by E. Kaplan, J.E. Marsden, and K.R. Sreenivasan (Springer, Berlin, 2002).
  4. B. Derrida and Y. Pomeau, Europhys. Lett. 1, 45 (1986).
  5. S. A. Kauffman, C. Peterson, B. Samuelsson, and C. Troein, Proc. Natl. Acad. Sci. 100, 14796 (2003).
  6. S. E. Harris, B. K. Sawhill, A. Wuensche, and S. A. Kauffman, Complexity 7, 23 (2002).
  7. S. A. Kauffman, C. Peterson, B. Samuelsson, and C. Troein, Proc. Natl. Acad. Sci. 101, 17102 (2004).
  8. V. Kaufman, T. Mihaljev, and B. Drossel, Phys. Rev. E 72, 046124 (2005).
  9. B. Samuelsson, and J. E. S. Socolar, Phys. Rev. E 74, 036113 (2006).
  10. P. Ramo, J. Kesseli, and O. Yli-Harja, Chaos 15, 034101 (2005).
  11. I. Shmulevich and S. A. Kauffman, Phys. Rev. Lett. 93, 048701 (2004).
  12. I. Shmulevich, E. R. Dougherty, S. Kim, and W. Zhang, Bioinf. 18, 261-274 (2002).
  13. A. Pomerance, E. Ott, M. Girvan, and W. Losert, Proc. Natl. Acad. Sci. 106, 8209 (2009).
  14. R. V.  Solé and B. Luque, Phys. Lett. A 196, 331 (1995).
  15. B. Luque and R. V. Solé, Phys. Rev. E 55, 257 (1997).
  16. J. J. Fox and C. C. Hill, Chaos 11, 809 (2001).
  17. D. S. Lee and H. Rieger, J. Phys. A 41, 415001 (2008).
  18. J. G. Restrepo, E. Ott and B. R. Hunt, Phys. Rev. E 76, 056119 (2007).
  19. M. E. J. Newman, Phys. Rev. Lett. 89, 208701 (2002).
  20. M. E. J. Newman and M. Girvan, Phys. Rev. E 69, 026113 (2004).
  21. R. Milo, et al., Science 298, 5594 (2002).
  22. J. G. Restrepo, E. Ott, and B. R. Hunt, Phys. Rev. Lett. 100, 058701 (2008).
  23. J. G. Restrepo, E. Ott and B. R. Hunt, Phys. Rev. Lett. 97, 094102 (2006).
  24. M. Nykter, et al., Proc. Natl. Acad. Sci. 105, 1897, (2008).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description