New Research In
Physical Sciences
Social Sciences
Featured Portals
Articles by Topic
Biological Sciences
Featured Portals
Articles by Topic
 Agricultural Sciences
 Anthropology
 Applied Biological Sciences
 Biochemistry
 Biophysics and Computational Biology
 Cell Biology
 Developmental Biology
 Ecology
 Environmental Sciences
 Evolution
 Genetics
 Immunology and Inflammation
 Medical Sciences
 Microbiology
 Neuroscience
 Pharmacology
 Physiology
 Plant Biology
 Population Biology
 Psychological and Cognitive Sciences
 Sustainability Science
 Systems Biology
The effect of network topology on the stability of discrete state models of genetic control

Edited by Peter J. Bickel, University of California, Berkeley, CA, and approved February 26, 2009 (received for review January 7, 2009)
Abstract
Boolean networks have been proposed as potentially useful models for genetic control. An important aspect of these networks is the stability of their dynamics in response to small perturbations. Previous approaches to stability have assumed uncorrelated random network structure. Real gene networks typically have nontrivial topology significantly different from the random network paradigm. To address such situations, we present a general method for determining the stability of large Boolean networks of any specified network topology and predicting their steadystate behavior in response to small perturbations. Additionally, we generalize to the case where individual genes have a distribution of “expression biases,” and we consider a nonsynchronous update, as well as extension of our method to nonBoolean models in which there are more than two possible gene states. We find that stability is governed by the maximum eigenvalue of a modified adjacency matrix, and we test this result by comparison with numerical simulations. We also discuss the possible application of our work to experimentally inferred gene networks.
Boolean networks have been extensively investigated as a model for genetic control of cells (1, 2). In this model, each gene is represented by a node of a network, and each node has one of two states: on—i.e., producing (“expressing”) its target protein—or off. Directed links between genes indicate that one gene influences the expression of another, either through the expressed protein directly binding to DNA or to other signaling pathways that modulate transcription of a gene. In the standard Boolean network model, the system evolves in discrete time steps (t = 0,1,2,…), and at each step the state of every node is simultaneously updated according to some function of its inputs. This function approximates the action of activators (proteins that act to increase expression of a given gene) or inhibitors (proteins that act to reduce expression). Although this model might seem to be an oversimplification considering the complex kinetics involved in all steps of a transcription pathway, experimental evidence suggests that real biological systems are, in some cases, reasonably wellapproximated by Boolean networks (3).
In 1969, S. A. Kauffman (1) introduced a type of Boolean network known as an N − K network. In this model, there are N nodes each having exactly K input links, and the nodes from which these input links originate are chosen randomly with uniform probability. We refer to the number of input (output) links to (from) a node as the indegree (outdegree) of that node. At any given time t, the system state can be represented as an N vector whose ith component σ_{i}(t) is either zero or one, where zero or one corresponds to off or on and i = 1,2,…,N. There are 2^{N} possible states. The function determining the time evolution at each node is defined by a random, timeindependent, 2^{K}entry truth table. Because this is a finite, deterministic system, there is always an attractor: eventually, the system must return to a previously visited state (finiteness), after which the subsequent dynamics will be the same as for the previous visit (determinism). These attractors can be fixed points or periodic orbits. By using the Hamming distance between two states [i.e., the number of nodes for which the σ_{i}(t) disagree] as the distance measure, the system exhibits both what is termed a “chaotic” (or unstable) regime, where the distance between typical initially close states on average grows exponentially in time, as well as a stable regime, where the distance decreases exponentially. Between the two there is a “critical” regime. (Here, by “close” we mean that the Hamming distance is small compared with N.)
As a model of genetic control, these attractors have been postulated to represent a specific pattern of protein expression that defines the cell's character (1). In singlecelled organisms, these attractors might be taken to correspond to different cell states (growing, dividing, starving, heat or pHshocked). In multicellular organisms, different cell types (muscle, nerve, liver, etc.) have different expression patterns, and, within each type, a cell could be in a variety of states (resting, “activated,” dividing, etc.) that each correspond to different expression patterns. Boolean network approximations have been successful in predicting the gene expression time sequence of the segment polarity gene network in Drosophilia, a model for embryonic development where individual cells turn specific proteins on and off in patterns that guide the growth of certain organs and structures (3). Because the protein expression pattern of the cell is modeled from the state of the corresponding Boolean network, the question of the stability of the network then becomes important: do small perturbations in the expression pattern, due perhaps to chemical fluctuations, die out quickly, returning the cell to its original state, or do they quickly grow, pushing the cell into another state? The purpose of this article is to examine the stability of network dynamics in the context of discrete state models of gene networks.
One motivation for the consideration of dynamical stability is its possible relevance to cancer. Specifically, we hypothesize that dynamical instability of a gene network might be a causal mechanism contributing to the occurence of some cancers. We emphasize that this hypothesis is distinct from the previous hypothesis of “genomic instability” as a cause of cancer (4). In particular, genomic instability has been defined (as in the glossary of Nature Genetics) as “the failure to transmit an accurate copy of the entire genome from one cell to its two daughter cells.” In contrast, the instability we refer to is that of the dynamics of a given gene network, and we use the term “dynamical network instability” (DNI) to distinguish this condition. We speculate that DNI might arise from mutations and that, once established, as cells divide, DNI could lead to widely varying gene expression patterns from cell to cell. We emphasize that DNI implies that this variation would arise even in the absence of further mutation. That is, similar to the concept of chaos in continuousstate dynamical systems (e.g., ref. 5), DNI causes exponential sensitivity of typical system trajectories to small changes, which we speculate may lead to many different outcomes in the course of cell division. Recent microdissection results indicate wide variations in gene expression patterns, even for nearby cells within the same cancerous tissue (6). This variability provides a basis for understanding why cancer can adapt and evade treatment (7).
Another motivation for our study is the argument, put forward by Kauffman (2), that evolution favors gene networks that are on the border between stability and instability (8–11). Whether or not our cancer hypothesis or Kauffman's stabilityborder hypothesis holds, the question of dynamical stability of such networks is crucial to their understanding and use as models.
Although previous works have addressed the question of dynamical network stability in simple, specific types of random networks (e.g., N − K nets), in this article we address the question of dynamical network stability for general network topology and node attributes. We also consider nonsynchronous update and extend the considerations to nonBoolean models allowing for the possibility of nodes having more than two states. Thus, our work provides a potentially enhanced framework for modeling and using the discrete state network paradigm. In particular, we consider how our network stability considerations can be employed on experimentally derived gene networks.
In the original N − K nets as proposed by Kauffman, the truth table output governing node dynamics was randomly chosen with on and off having equal probability. Subsequently, it was shown that if the truth table output was biased such that p denotes the probability of randomly assigning an off output, the transition between the stable and chaotic regimes depends on p (12). We term p the “expression bias.” Additionally, networks with a distribution of indegrees, but no in/outdegree correlation, have been considered in refs. 13–16, and it has been shown that the nodal indegree average, 〈K^{in}〉, suffices to determine the stability. (Here, 〈·〉 indicates average of a nodal quantity over all nodes.) Specifically, the critical average number of connections, K_{c}, governing this transition is where the network is stable for 〈K^{in}〉 < K_{c}, unstable for 〈K^{in}〉 > K_{c}, and critical for 〈K^{in}〉 = K_{c}. Aldana and Cluzel (16) considered the consequences of Eq. 1 in the case of networks with scalefree topology (17), i.e., the indegree (outdegree) distribution P(K^{in}) [P(K^{out})] is a powerlaw: P(K) ∝ K^{−γ}. (Since every outlink from a node is an inlink to another, 〈K^{in}〉 = 〈K^{out}〉; thus, the result is unchanged whether it is the in or outdegree that has powerlaw scaling.)
Recently, some authors have noted, but not numerically tested, a generalization of Eq. 1 that takes into account nodal correlations between the indegree and outdegree characterized by the joint K^{in} – K^{out} degree distribution function P̃ (K^{in}, K^{out}). In this case, the critical transition occurs at (18)
We emphasize that Eq. 1 was derived in the annealed approximation (see later discussion) for networks with a given in or outdegree distribution P(K) and with the complimentary links completely random, and that Eq. 2 uses only the additional information contained in the nodal in/outdegree correlation. Furthermore, all nodes (“genes”) were taken to have the same p value. However, gene networks, in common with real networks occurring across a broad range of applications, can be expected to deviate substantially from the above simple network model. Examples of network properties that could make previous analyses of network stability inapplicable are assortativity/disassortativity (19) (the tendency for highly connected nodes to prefer or avoid linking to other highly connected nodes) and community structure (20) (the existence of highly connected, sparsely interconnected subgraphs), two properties that are not captured in the degree distributions. Additionally, these properties may have biological implications. For example, a recent article (21) examined gene interaction networks from cancerous tissue and found significant community structure, as well as positive correlation between the indegree and outdegree of nodes; additionally, protein interaction networks have been shown to exhibit significant disassortativity (19, 22). Furthermore, for modeling purposes, it might be important to allow the expression bias p to vary from node to node [as an extreme example, socalled housekeeping genes (23) have a predominant tendency to be on, corresponding to low p, unlike other genes]. In this article, we derive and test the stability criterion for large networks with arbitrary network topology and heterogeneous expression biases. In particular, our theory evaluates the stability of any given network with its specific topology (i.e., its adjacency matrix A defined subsequently) and its nodespecific expression biases. We show that stability is determined by the largest eigenvalue of a modified adjacency matrix, and we numerically test this criterion.
With respect to real gene networks, the synchronous update at integer times (t = 0,1,2,…) used in the above models represents an additional deviation from the real situation, where chemical kinetics and transport processes can be expected to introduce nontrivial dynamics. As a partial step toward remedying this (and to make Boolean approximations suitable for atmospheric and geophysical processes), Ghil and Mullhaupt (24) consider a generalization in which t is a continuous variable and σ_{i}(t) depends on σ_{j}(t − τ_{ij}), where τ_{ij} is a delay time that can be different for each link from j to i. The original formulation (e.g., in refs. 1, 12–16) corresponds to τ_{ij} = 1 for all i,j. We will argue and numerically confirm that the criterion determining the stability/instability border of this generalization of the Boolean network model is the same as that for the synchronous update models.
In addition to nonsynchronous update, another generalization of Boolean networks that we will examine is models in which each node i is allowed to have one of S_{i} possible discrete states (e.g., for S_{i} = 3, we label the states σ_{i} ∈{0,1,2}, and for Boolean networks S_{i} = 2 for all i). This model may be closer to the behavior of actual cells, and models with multiple states can be related to certain piecewise ODE models of transcription (25, 26). The general model using arbitrary, multivalued truth tables has been previously treated in the special case of N − K networks with all nodes having the same number of possible states S. In the case where each possible state is equally likely, the critical number of inputs has been shown by Sole et al. (see http://arxiv.org/abs/adaporg/9907011) to be
The applicability of our work to any specific network and set of nodewise expression biases may be of particular interest in situations where experimental data provide the possibility of estimating a gene network and expression biases. Such information could be used as input to our method which could give an indication of the stability of a given experimentally derived network. The possibility that such analyses may be feasible becomes more and more likely with the rapid technological advances in obtaining new types of highquality, quantitative data useful for deducing gene networks. For example, such analyses could be used to address the hypothesis that dynamical instability of gene networks is connected with the occurence of cancer.
Model
Deterministic Boolean networks are formally defined by a state vector Σ(t) = [σ_{1}(t)σ_{2}(t)…σ_{N}(t)]^{T}, where σ_{i} ∈ {0,1}, and a set of update functions f_{i}, such that where j(i,1), j(i,2),…,j(i,K_{i}^{in}) denote the indices of the K_{i}^{in} nodes that input to node i; we denote this set of nodes by K_{i} = {j(i,k)k = 1,2,…K_{i}^{in}}. The update function f_{i} is defined at each node i by specifying a truth table whose outputs are randomly populated. Previous analytic results assumed a constant expression bias for all nodes; however, we allow that, in the truth table for node i, output entries are randomly assigned zero with probability p_{i} or one with probability 1 − p_{i}. In the case of uniform expression bias, we drop the subscript and use the notation p ≡ p_{i}.
We consider the interaction structure of this system as a graph where the nodes represent individual elements of the state vector, and a directed edge is drawn from node j to node i if j ∈ K_{i}. An adjacency matrix A is defined in the usual way: A_{ij} is one if there is a directed edge from node j to node i and zero otherwise.
The stability of a large Boolean network is defined by considering the trajectories resulting from two close initial states, Σ(t) and Σ̃(t). To quantify their divergence, the Hamming distance of coding theory is used:
To study the stability of an N − K Boolean network, Derrida and Pomeau (12) considered an annealed situation and calculated the probability that, after t steps, a node state is the same on two trajectories that originated from initially close conditions. [This calculation was later generalized to variable indegree (13–15) and joint degree distribution (18).] In Derrida and Pomeau's “annealed” situation, at each time step t the truth table outputs and the network of connections are randomly chosen. The actual situation of interest, however, is the case of “frozenin” networks, where the truth table and network of connections are fixed in time. It has been commonly assumed that analytical results obtained for the annealed case are a good approximation to the frozenin case (e.g., refs. 12–15). We also adopt this view in a modified form, and we will test its predictions with numerical simulations.
The randomization of the network of connections at each time step while keeping the degree distribution fixed carries the implicit assumption that there is no additional dynamically relevant structure in the frozen network other than that contained in the joint degree distribution P̃(k^{in}k^{out}). To avoid this assumption, we obtain theoretical results for a different annealing protocol, which we term “semiannealed.” In this semiannealed procedure, we keep the network fixed (i.e., the adjacency matrix A does not change with time), and we envision randomly assigning the output entries of the truth table of each node i at every time t according to the timeindependent expression bias p_{i} assigned to node i. We then imagine tracking the probability that individual node states σ_{i}(t) and σ̃(t)differ over time with an N dimensional difference vector, whose components are
The case where both network states are exactly the same corresponds to y_{i}(t) = 0, which is a fixed point of Eq. 6. To determine the stability of this fixed point, we linearize Eq. 6 around y(t) = 0 for small perturbations: where A_{ij} are the elements of the adjacency matrix A. Eq. 7 can be written in matrix form as y(t + 1) = Qy(t) where Thus the largest eigenvalue of this matrix λ_{Q} governs stability: Since Q_{ij} ≥ 0, the Perron–Frobenius theorem (31) guarantees that λ_{Q} is real and positive. We also note that, for any given adjacency matrix A and assignment of q_{i}'s to nodes, Eq. 6 can be iterated numerically to predict the expected timeasymptotic saturation value of the difference in two initially nearby states when evolved to steady state. We numerically test this prediction, as well as the stability criterion in Eq. 9 in the next section.*
As a special case of interest, if the q_{i} are uniform, q_{i} ≡ q, then λ_{Q} = qλ, where λ is the maximum eigenvalue of the adjacency matrix. This yields the critical condition, Furthermore, for the case of a large network whose links are randomly assigned subject to a joint probability distribution P̃(K^{in}K^{out}) at each node (with no assortativity), the meanfield approximation for the largest eigenvalue is (32) where, since 〈K^{in}〉 = 〈K^{out}〉, we use the notation 〈K〉≡〈K^{in}〉 = 〈K^{out}〉. Eqs. 10 and 11 yield the same criterion as in Eq. 2. In the case where K^{in} and K^{out} are uncorrelated, P̃(K^{in}K^{out}=P_{in}(K^{in})P^{out} (K^{out})and 〈K^{in}K^{out}〉 = 〈K〉^{2}, yielding Eq. 1.
The eigenvalue of random network adjacency matrices with assortativity has been considered in ref. 32, which defines an assortativity measure ρ as where 〈K_{i}^{in}K_{j}^{out}〉_{e} denotes an average over all links (i, j) from node i to node j. The network is assortative (disassortative) if ρ > 1 (ρ < 1). For ρ near one, the largest eigenvalue λ is approximately given by (32) Thus, it is predicted that, for uniform q, assortativity (disassortativity) decreases (increases) the critical q value.
In the case of nonuniform q_{i}, we have recently generalized Eq. 11 to obtain an analogous meanfield approximation to λ_{Q} without assortativity or community structure, Our derivation of 14 will be published elsewhere. From Eq. 14, we see that correlation (anticorrelation) between q and K^{in}K^{out} decreases (increases) network stability and that, in the absence of correlation, the result is similar to that for a uniform q, λ_{Q} ≈ 〈q〉〈K^{in}K^{out}〉/〈K〉, with 〈q〉 replacing the uniform q (Eqs. 9 and 10).
We now consider the generalization to allow any number of discrete node states. We denote the number of possible states of node i by S_{i}, and we label the possible states 0,1,2,…,S_{i} − 1. The number of possible inputs to i from the set K_{i} of nodes that influence it is Π_{j ∈ Ki}S_{j}. For each of these possible inputs, the truthtable function f_{i} in Eq. 4 assigns one of the S_{i} possible states to node i. Similar to the Boolean case, we take the assignment to be random and to have an “expression bias” p_{i,s} for each of the s = 0,1,2,…S_{i} − 1 node states, where p_{i,s} denotes the probability that f_{i}, for a given set of inputs, assigns the state s to node i, and ∑_{s}p_{i,s} = 1 for all nodes i. As in the Boolean case, we can then introduce the sensitivity q_{i} giving the probability that two different sets of inputs result in a different updated state of node i, which, in the random truth table case, With this definition, we see that all our previous reasoning still applies, and Eqs. 6–9 hold with this generalized expression for the node sensitivities and with y_{i}(t) interpreted as the probability of disagreement between σ_{i}(t) and σ̃_{i}(t). In the case of uniform number of node states, S_{i} ≡ S, and equal expression biases, p_{i,s} ≡ p_{s} = 1/S, among these states, Eq. 15 becomes q = 1 − 1/S, which, when combined with Eq. 10 yields the previous result in Eq. 3.
Finally, we note that our criticality criterion, λ_{Q} = 1, is unchanged by the presence of delays, as in the models of refs. 24, and only a slight modification is required of Eq. 6 [i.e., y_{j}(t − τ_{ij}) replaces y_{j}(t − 1)]. The condition λ_{Q} = 1 implies that the components of y in Eq. 6 are timeindependent. Thus we predict that the delays τ_{ij} do not influence the result, and the criticality condition in Eq. 9 is independent of the synchronous update structure of the most commonly used random Boolean network models. Similarly, the timeasymptotic steady state obtained by repeated iteration of Eq. 6 is, by definition, timeindependent and thus also does not depend on the τ_{ij} (although the τ_{ij} will influence the timedependent approach to the asymptotic steady state; see supporting information (SI) Appendix).
Statistical Methods
We numerically test the above predictions on several classes of Boolean networks with uniform sensitivity (i.e., q_{i} = q is the same for all nodes):

Random networks with K^{in} = K^{out};

Random networks with imperfect correlation between K^{in} and K^{out};

Networks with assortativity or disassortativity; and

Networks constructed as in (a) but with a substantial number of feedforward loops.
We additionally test our predictions on two classes of networks with nonuniform sensitivities:

Networks constructed as in (a) but where nodes have different sensitivities correlated with the degrees of the nodes; and

Networks with significant community structure, where the two communities have different, uniform sensitivities.
Finally, we test our generalization to more than two node states on networks of type (a) but with S_{i} = 4 for all nodes. For types (a)–(c) and (e), we use networks with truncated powerlaw degree distributions. (Evidence for the presence of this type of distribution in gene networks has been seen in ref. 33.)
The algorithms for constructing the networks of types (a)–(c) are as follows. (i) Establish the in and outdegrees for each node, which are drawn from a distribution, where γ = 2.1 and K^{max} = 15 (Boolean case) or K^{max} = 8 (S_{i} = 4 case). The outdegree is initially set to the indegree. (ii) Randomly swap the outdegrees between pairs of nodes. If maximal correlation between in and outdegrees is desired, as in (a), this step is skipped so that K^{in} = K^{out} and 〈K^{in}K^{out}〉 is maximal. A completely uncorrelated network has every nodal outdegree swapped exactly once, yielding 〈K^{in}K^{out}〉 = 〈K〉^{2}. The quantity 〈K^{in}K^{out}〉, which approximately determines λ by Eq. 11, can thus be tuned by the number of nodes that have their outdegrees swapped. (iii) Place links randomly between nodes subject to the constraints of the specified in and outdegrees assigned at each node by the “configuration model” (34). (iv) If networks with assortativity (disassortativity) are desired, as in (c), perform a given number of link swaps, as in ref. 32, that increase (decrease) the assortativity ρ in Eq. 12. In all cases we employ networks with N = 10^{4} and two initial conditions separated by a Hamming distance of 100. In the SI Appendix we discuss finite size effects that can occur for smaller N.
We emphasize that, although we determine our networks randomly, in our numerical experiments we do not average over this randomness. Rather, we generate one random network for each experiment and examine the resulting behavior of that specific network.
Results
We test the steadystate predictions of Eq. 6 and the criticality condition of Eq. 9 in Fig. 1, and compare the calculated critical parameters to the meanfieldtype approximations of Eqs. 11, 12, and 14 in SI Appendix. To compare the theory (solid curves in Fig. 1) to experimental measurements of the Hamming distance from numerical evolution of true frozen Boolean dynamical systems (markers in Fig. 1), we simulate Eq. 6 to steadystate and calulate the expected fractional Hamming distance from the distance vector with When simulating test Boolean networks, we calculate the steadystate fractional Hamming distance as the average from t = 90 to t = 100 when all delays are the same ([τ]_{ij} = 1) or from t = 490 to t = 500 and normalize by the network size N. These times are well after the steadystate value is reached (see SI Appendix). Each experimental data point in Fig. 1 corresponds to a single realization of interconnections averaged over 100 realizations of the timeindependent truth table with specified sensitivity as before.
In/Outdegree Correlations and Heterogeneous Time Delay.
Fig. 1A shows the steadystate Hamming distance as a function of the sensitivity for one network of type (a) (λ = 4.4) and two of type (b) (λ = 2.9,2.3). The closed markers in the figure represent experiments with uniform delay τ_{ij} = 1 on all links, while the open markers correspond to experiments where half the links, randomly chosen, have τ_{ij} = 10 and the remainder have τ_{ij} = 1. Importantly, the degree distributions are the same for all three networks, and we attain different λ values by varying the correlation between the indegree and the outdegrees. We see from Fig. 1A that there is close agreement between the theoretical prediction and the experimental results and our prediction that the presence of delays does not change the stability is confirmed. Additionally, the measured steadystate Hamming distance is essentially zero below the critical value of the sensitivity, q_{crit} = 1/λ (this point is indicated by vertical downward arrows in Fig. 1A). We emphasize that the degree distributions (and hence 〈K〉) are the same for the networks in Fig. 1A, and thus, if the in/outdegree correlation were ignored, the observed difference between the stability conditions for these networks would not be predicted.
Assortativity/Disassortativity.
Fig. 1B shows results obtained when significant assortativity or disassortativity is present [type (c) networks]. In this experiment, as well as all those reported below, the delays are all uniform. The networks under consideration have the same joint degree distribution with K^{in} = K^{out}. However, each of the networks have very different assortativities (ρ = 0.52,1.0,1.7, defined in Eq. 12), which yield different largest eigenvalues (λ = 3.0,4.4,9.9). Since the joint degree distributions are the same, Eq. 2 would predict that the three networks have the same stability characteristics. However, since their eigenvalues are very different, we predict that, as observed, the transitions of the three networks occur at different values of q. Again the theoretical predictions of q_{crit} are indicated by vertical arrows.
Motifs.
Random construction of networks, as used in the networks above, is expected to yield networks that are locally treelike (32). However, we note that biological and other types of networks often have motifs (small subgraphs) that occur with higher frequency than in randomly constructed networks (30). For gene networks of Escherichia coli and Saccharomyces cerevisiae, it was found that the number of feedforward loop motifs (see Fig. 1C Inset) is significantly enhanced compared with the expected number in a randomly constructed network. In these networks, the number of feedforward loops per node c is roughly 0.1. Thus, we consider a network of type (a) (λ ≈ 2.9 and N = 10^{4}) after adding 1,000 (c = 0.1) and, in an extreme case, 2,000 (c = 0.2) feedfoward loops. To add a feedforward loop, we randomly choose a node A, follow a random output to node B, and follow a random output of B to node C. We then add a link from node A to node C. We do this a given number of times, avoiding nodes that already participate in an added feedforward loop. In Fig. 1C, we see that the semiannealed theory of Eq. 6 (solid curve) again agrees well with our numerical experiments (solid markers). Based on such results, we believe that the locally treelike network requirement does not invalidate application of our method to real gene networks. We also note that the critical point is essentially unchanged by the addition of loops (adding links only slightly increases the largest eigenvalue), however, more feedforward loops tend to increase the steadystate Hamming distance for q > q_{crit}.
Application to S. cerevisiae.
As a real biological example, we include in the SI Appendix a graph similar to those in Figs. 1A – C using a published network for the yeast S. cerevisiae (35).
Heterogeneous Correlated Sensitivities.
Fig. 1D demonstrates the effect of a distribution of q_{i}'s on the stability of a network with K_{i}^{in} = K_{i}^{out} = K_{i} and with correlation between the nodal values of q_{i} and K_{i}, i.e., type (d) networks. We consider two situations, one where 〈qK^{2}〉/(〈q〉〈K^{2}〉) is maximal, and one where it is minimal. The q_{i} 's are drawn from a uniform distribution centered at q_{0} (the abscissa in the figure), with width Δq = 0.1. Maximal (minimal) 〈qK^{2}〉 is attained by assigning the largest q_{i} to the node with the largest (smallest) K_{i}, the second largest q_{i} to the node with the second largest (second smallest) K_{i}, and so on. As can be seen from the figure, there is good agreement between the semiannealed theory and the numerical experiments, and the two networks become unstable at different values of q_{0}. (Vertical arrows again indicate the points where λ_{Q} = 1.)
Community Structure.
Fig. 1E shows our results for a case where there is community structure and communitydependent sensitivity. To construct the networks in Fig. 1E, consider the case where there are two communities, and we assign a link from node i in community a to node j in community b with probability θ_{ab}. We impose the additional constraints that θ_{aa} = θ_{bb} ≡ θ_{∪} and that θ_{ab} = θ_{ba} ≡ θ_{∩}, and the sizes of the two communities are the same, N/2. We take 〈K^{in}〉 = 〈K^{out}〉 = 〈K〉 = (θ_{∪} + θ_{∩})N to be the same for both communities, and we also assume that communities a and b have different sensitivities q_{a} and q_{b}, respectively. As θ_{∩} is increased from zero to θ_{∩} = θ_{∪}, λ_{Q} changes from the case of two completely separated communities to one of a single random network. Communities a and b have equal sizes of 5,000 nodes, community a has q_{a} = 0.5, and community b has q_{b} = 0.1. To vary λ_{Q}, we vary θ_{∪} and θ_{∩}, keeping their sum constant to maintain constant 〈K〉. As with the curves in Fig. 1A – C, the transition to chaos is governed by λ_{Q} (λ_{Q} = 1 at the vertical arrow), and Eq. 6 (solid curve) accurately predicts the numerically observed (solid circles) steadystate Hamming distance.
NonBoolean Models.
Fig. 1F illustrates an application to a case in which there are more than two possible states at each node. In particular, we consider S ≡ S_{i} = 4 possible states at each node. (Since the number of possible inputs to the truth table for node i in this case is 4^{Kiin}, we take K^{max} = 8 due to memory constraints.) Labeling the possible node states σ_{i} ∈ {0,1,2,3}, we take nodes to have uniform expression biases for occurrence of statelabel 0, p_{0} ≡ p_{i,0}, from 0 to 1. The three remaining labels (S = 1,2,3) also have uniform biases for all nodes i, p_{s} ≡ p_{i,s} = (1 − p_{i,0})/3. From Eq. 15, q ≡ q_{i} = 1 − [p_{0}^{2} + (1 − p_{0})^{2}/3], which has a maximum q_{max} = 0.75 at p_{0} = 0.25. As can be seen in the figure, the predicted fraction of nodes with differing states ȳ (solid curve) also has a maximum there. It is again seen that the measurements (markers) are well predicted by the theory.
Discussion
In this article, we have presented theoretical results (Eqs. 6 and 9) which predict the steadystate Hamming distance between states evolved from two nearby initial conditions and the stability of a given network. These results are derived using the hypothesis that a theory derived in the semiannealed case approximates the true situation, where by semiannealed we mean that the network of connections is frozen, but the truth table at each node is randomly reassigned at each time step. For large networks, this approximation was found to give excellent agreement with the true case of frozen connections and frozen truth tables. Our semiannealed hypothesis does not rely on gross statistical properties of the network, but instead uses the specific network topology, as characterized by the network adjacency matrix, and the individual node sensitivities, to make predictions.
We tested our theoretical predictions with numerical experiments. Previously unaddressed issues that we considered include the effects of assortativity, nonuniform time delay, nonuniform sensitivity, motifs, and community structure. In all cases tested we found good agreement with our theory.
The theory that we have presented and tested above may represent a step forward in facilitating the application of discrete state dynamical network models to biological systems. Given a specific genetic interaction network and an estimate of the node sensitivities, Eq. 9 predicts the stability of that particular network directly from the adjacency matrix. Curated networks already exist in the literature for model singlecellular systems, and new algorithms continue to be developed for inferring interaction networks from a wide range of data sources (microarray experiments, genome sequencing, etc.). We note that such a procedure has the advantage that, because the actual experimentally determined network is employed, topological aspects such as nodal in/out degree correlation, assortativity, community structure, etc., do not first have to be determined and then statistically modeled. Thus, by use of our stability criterion (Eq. 9), there is the potential that future analysis may be able to evaluate a supposed relationship between the stability characteristics of various networks and their functioning. For example, one might test whether cancer gene networks are less stable than those in healthy tissue. This could lead to the strong variations in gene expression observed in cancerous tissue (9), even when the underlying gene network is unchanged. We are currently pursuing research along this line.
Acknowledgments
We thank L. Staudt for discussions. The work was supported by the National Science Foundation (Physics) and by Office of Naval Research Contract N00001410734. A. P. was supported in part by the National Cancer Institute intramural program.
Footnotes
 ^{1}To whom correspondence should be addressed. Email: pomeranc{at}umd.edu

Author contributions: A.P. and E.O. designed research; A.P. and E.O. performed research; A.P. and E.O. analyzed data; and A.P., E.O., M.G., and W.L. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

This article contains supporting information online at www.pnas.org/cgi/content/full/0900142106/DCSupplemental.

↵* We note that Eq. 8 and the condition λ_{Q} = 1 also occurs in the treatment (28) of site percolation on directed networks where different sites have different removal properties. A similar condition involving 〈K^{2}〉/〈K〉 also arises in percolation on undirected networks (29).
References
 ↵
 ↵
 Kauffman S
 ↵
 ↵
 Coleman W,
 Tsongalis G
 ↵
 Ott E
 ↵
 GonzalezGarcia I,
 Sole RV,
 Costa J
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 Aldana M,
 Cluzel P
 ↵
 Barabasi AL,
 Albert R
 ↵
 ↵
 ↵
 ↵
 Cui Q,
 et al.
 ↵
 Maslov S,
 Sneppen K
 ↵
 ↵
 Ghil M,
 Mullhaupt A
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 Milo R,
 et al.
 ↵
 ↵
 ↵
 ↵
 ↵
Citation Manager Formats
More Articles of This Classification
Biological Sciences
Biophysics And Computational Biology
Related Content
 No related articles found.