A natural class of robust networks
See allHide authors and affiliations

Edited by Peter G. Wolynes, University of California at San Diego, La Jolla, CA, and approved May 1, 2003 (received for review November 6, 2002)
Abstract
As biological studies shift from molecular description to system analysis we need to identify the design principles of large intracellular networks. In particular, without knowing the molecular details, we want to determine how cells reliably perform essential intracellular tasks. Recent analyses of signaling pathways and regulatory transcription networks have revealed a common network architecture, termed scalefree topology. Although the structural properties of such networks have been thoroughly studied, their dynamical properties remain largely unexplored. We present a prototype for the study of dynamical systems to predict the functional robustness of intracellular networks against variations of their internal parameters. We demonstrate that the dynamical robustness of these complex networks is a direct consequence of their scalefree topology. By contrast, networks with homogeneous random topologies require finetuning of their internal parameters to sustain stable dynamical activity. Considering the ubiquity of scalefree networks in nature, we hypothesize that this topology is not only the result of aggregation processes such as preferential attachment; it may also be the result of evolutionary selective processes.
Complex protein and genetic networks are large systems of interacting molecular elements that control the propagation and regulation of various biological signals (1–3). They constitute essential classes of biological computations reflecting vital cellular processes such as the regulation of the cell cycle and gene expression. Because some intracellular processes are crucial for the survival of the cell they need to be achieved with reliability. For example, the variation of the concentration of network components in a metabolic network may affect numerous processes. The intricate architecture of such networks raises the question of the stability of their functioning. How can such complex dynamical systems achieve important cellular tasks and remain stable against the variations of their internal parameters?
Three decades ago, Savageau (4, 5) hypothesized that robustness was an essential property of some genetic networks whose functioning would be preserved even if some of their components were produced in various quantities. More recently, several compelling theoretical (6, 7) and experimental (8, 9) studies demonstrated that key processes of specific intracellular networks exhibited a robust behavior to variations of biochemical parameters. Robustness has emerged as a fundamental concept for the characterization of the dynamical stability of biological systems.
Are there universal design principles that would determine whether a given class of networks achieves important intracellular processing with reliability (10)? It would then be possible to predict dynamical properties underlying cell functions without a full knowledge of the molecular details. In particular, it would be important to characterize dynamical robustness as the ability for the network to perform a sequence of biological tasks in the presence of perturbations (3). Driven by such considerations, we analyzed when dynamical robustness may be a direct consequence of the network's architecture.
As a testbed for the study of the dynamical properties of complex biological networks we chose to model molecular activities by a simple twostate model closely related to the one proposed by Kauffman (11, 12). This model was originally created for the study of the dynamics of genetic networks, but was later applied to evolution and social models and became a prototype for the study of dynamical systems. In this model, the network's architecture follows a random topology in which every element interacts, on average, with K other elements. Each network element has two functional states, active and inactive. The state of a given element is determined by its interactions with other elements of the network (see Appendix). A generic biochemical parameter for the whole system is the probability ρ that an arbitrary element of the network becomes active after interacting with other elements (13). This probability can be inferred in principle from the reaction rates and the relative concentration of each element. The nature of the dynamical process performed by a given network is determined by its topological parameter K and its biochemical parameter ρ. Under the simple assumption of a random network topology, a very rich and unexpected dynamical behavior of the network was found (11, 12, 14). During the time evolution of the system, the network elements pass through different states until they reach a cyclic behavior. Different cycles are possible. Each cycle represents a variety of intracellular tasks. Two regimes of network activity exist: chaotic and robust. In the chaotic regime a perturbation in the state of a single element can make the system jump from one cycle to another. In the robust regime all such perturbations die out over time. The transition from one regime to the other is controlled by the parameters K and ρ and is determined by the equation [1] Fig. 1a defines the two regimes of network activity with drastically different dynamical properties. In the chaotic regime, which entirely dominates the parameter space (Kρ), small variations of the initial concentration of only one component would ruin the network function. Conversely, the robust regime would give rise to a reliable dynamics for which the network activity would be insensitive to perturbations in the initial state of the network elements. With a random network topology this regime is attained only for a narrow range of the parameters K and ρ.
Unfortunately, this model is inadequate to account for the functioning of biological networks that have heterogeneous architecture. Such topological heterogeneity is present in the human genome. For example, the expression of the βglobin gene is controlled by >20 regulatory proteins. The fibroblast or the plateletderived growth factor activation induces the activation of >60 other genes (15). Also, the zincfinger protein Sp1 controls the expression of >300 genes (16). Under an assumption of a random homogeneous architecture with a mean connectivity K = 20, robust behavior would be virtually attained only for ≈5% of the possible values of ρ. Fig. 1b shows the fraction of the interval (0,1) for which the parameter ρ gives rise to networks with a robust behavior. It is apparent that for a random topology only the finetuning of the parameters K and ρ in a narrow interval would allow networks to operate with a robust behavior. An alternative approach is to substitute the random network topology for a more realistic one (17).
In recent years, the analyses of the topology of large intracellular networks have revealed a common architecture (18–20). In these natural networks not all elements are equivalent. A small, but significant, fraction of the elements are highly connected, whereas the majority of the elements are poorly connected. This architecture is called scalefree topology and was found to be ubiquitous in networks as diverse as social (21), ecological (22), and genetic, protein networks (18, 23), and the World Wide Web (24, 25). Although the effect of deleterious perturbations on the topology of these networks has been thoroughly studied (23, 26), their dynamical properties need to be explored.
In view of the ubiquity of scalefree networks in nature, we chose to implement the scalefree topology into the twostate model. This topology is defined such that the probability P(k) that an arbitrary element of the networks interacts with exactly k other elements follows the power law P(k) = [Z(γ)k^{γ}]^{–}^{1}. The parameter γ is called the scalefree exponent, and Z(γ) is the normalization factor. When γ increases both the fraction of highly connected elements and the mean connectivity of the network decrease. In the limit γ → ∞ the network becomes a homogeneous random network with K = 1.
Our approach aims to characterize how the network architecture shapes its dynamical properties and how one particular topology favors a class of networks with robust dynamical behavior (10). Although the average connectivity K is a relevant parameter to characterize the homogeneous topology of random networks, it becomes irrelevant to describe the highly heterogeneous scalefree topology. Thus, the only relevant parameter that characterizes the architecture of scalefree networks is the scalefree exponent γ.
For networks with scalefree topology, the transition from chaotic to robust dynamics is determined by the equation (see Appendix) [2] The values of γ and ρ that satisfy this transcendental equation are plotted on Fig. 2a. The results of our model reported in the diagram reveal a large class of robust networks. This class is defined for networks with topological parameters γ > 2. These networks can operate with a robust dynamical behavior for which the mean biochemical parameter ρ spans over a large range of values. In particular, networks with a parameter γ > 2.5 would display a robust dynamics for any value of ρ. The transition from the chaotic to the robust regimes occurs for values of γ in the interval [2,2.5] and different values of ρ. Fig. 2a reveals another class of networks with a topological parameter γ < 2 that exhibit chaotic behavior. This dynamical behavior is characterized by an extreme sensitivity to small variations of the initial states of the network elements. Our model shows that the emergence of a large class of robust networks is a direct consequence of the network architecture.
It is interesting to note that experimental values of γ, extracted from real intracellular networks, range in the interval [2, 2.5]. To illustrate this general remark, we indiscriminately report on a histogram 46 published values of scalefree exponents (Fig. 2c). These values were obtained not only from real biological systems but also from many other scalefree networks (20, 23, 25). Remarkably, the majority of these exponents fall in the predicted interval where the transition from chaotic to robust behavior occurs.
This observation leads us to explore the dynamical response of the network to perturbations in the three different regimes: chaotic, critical, and robust. In a scalefree network not all elements are equivalent. We do not expect that perturbing highly connected elements have the same effect on the network dynamics than perturbing poorly connected elements. To test this idea we analyzed numerically how a perturbation produced on only one element affects the dynamics of the network. We computed the overlap (see Appendix) between two trajectories. One trajectory is computed with one element σ_{i} with k_{i} connections being perturbed at every time step of the temporal evolution of the network. The second trajectory is computed from the same initial condition but with no perturbation. Such overlap is a measure of the dynamical stability of the system under sustained perturbations. Fig. 3 shows that the stability of the network decreases with the connectivity k_{i} of the perturbed element σ_{i}. Remarkably, we found that networks operating in the robust regime are sensitive to perturbations applied to the highly connected elements. In contrast, perturbations produced on arbitrary elements, which in general are poorly connected, will have little effect on the network functioning. The heterogeneity of scalefree networks operating in the robust regime allows different dynamical responses to perturbations depending on the particular element being perturbed. This behavior cannot be observed in homogeneous random networks in which all of the elements are equivalent.
For a homogeneous network topology, stable behavior was achieved only for very restricted values of the parameters ρ and K. Conversely, the scalefree topology does not require a finetuning of the system parameters to attain robustness. This topology allows the existence of robust dynamics in heterogeneous networks. Furthermore, the dynamics of scalefree networks can change if the highly connected elements are perturbed. Thus, scalefree networks operating in the robust regime exhibit both, robustness and sensitivity to perturbations. The above results suggest that the dynamical robustness of a wide variety of biological processes may be a direct consequence of the network architecture. The ubiquity of scalefree networks in nature with a topological parameter 2 < γ < 3 could be the result of evolutionary processes (27). Such a large class of networks, whose key dynamical processes are robust, may have an essential evolutionary advantage in possessing a larger parameter space to adapt to new environmental conditions (27, 28). In our opinion, it would be surprising if nature did not exploit this convenient dynamical property of the scalefree topology.
This work describes the dynamics of networks where complex molecular activities are idealized with a twostate model (29). There exists a fair number of intracellular regulatory systems with components obeying the rules of a twostate model. For example, in signal transduction networks, the signal is processed through phosphorylation cascades where signaling molecules can be active or inactive upon phosphorylation. It is also common to describe regulatory transcription networks in this simple language: A transcription factor regulates positively or negatively the transcription of a given gene. Moreover, available data on the activity of molecular components from large intracellular networks are produced by highthroughput experiments such as twohybrid assays and CDNA microarrays. Because of technical constraints, the readout of these experiments is fashioned accordingly to a twostate format that facilitates the applicability of our approach.
The identification of a relationship between topology and dynamical properties of large biological networks primarily motivated our work. But it seems appropriate to think that in view of the significance of scalefree networks in other fields such as sociology and economics our results could be of some interest beyond the scope of biology.
Acknowledgments
We thank Leo P. Kadanoff, Leo Silbert, Tao Pan, Tamara Griggs, and Calin Guet for useful comments. This work was supported by the Materials Research Science and Engineering Center program of the National Science Foundation under Award 9808595 and National Science Foundation Grant DMR 0094569. M.A. also acknowledges the Santa Fe Institute of Complex Systems for partial support through the David and Lucile Packard Foundation Program in the Study of Robustness.
Appendix
A Boolean network is represented by a set of N elements, {σ_{1}, σ_{2},..., σ_{N}}. Each element has two possible states: active (1) or inactive (0). The value of each σ_{i} is controlled by k_{i} other elements of the network (see Fig. 4), where k_{i} is a random variable chosen from a probability distribution P(k). This probability distribution determines the topology of the network. We define K as the mean value of P(k). Let {σ_{i}_{1}, σ_{i}_{2},..., σ_{i}_{ki}} be the set of controlling elements of σ_{i}. We then assign to each σ_{i} a Boolean function f_{i}(σ_{i}_{1}, σ_{i}_{2},..., σ_{i}_{ki}). For each configuration of the controlling elements, f_{i}(σ_{i}_{1}, σ_{i}_{2},..., σ_{i}_{ki}) = 1 with probability ρ and f_{i}(σ_{i}_{1}, σ_{i}_{2},..., σ_{i}_{ki}) = 0 with probability 1 – ρ. Once the controlling elements and the Boolean functions have been assigned to every element in the network, the dynamics of the system are given by [A.1] We will denote as Σ(t) the configuration of the entire system at time t: Σ(t) = {σ_{1}(t), σ_{2}(t),..., σ_{N}(t)}.
The dynamical robustness of the network is analyzed by considering the trajectories Σ(0) → Σ(1) → · · · → Σ(t) and → → · · · → produced by two slightly different initial configurations, Σ(0) and , respectively. These trajectories can always be different under the temporal evolution of the system, or they can eventually converge to the same trajectory. An important quantity that reveals the robustness of the dynamics against perturbations in the initial configuration is the overlap between Σ(t) = {σ_{1}(t), σ_{2}(t),..., σ_{N}(t)} and , defined as where 〈·〉 represents the average over all possible initial configurations and network realizations. If lim_{t}_{→}_{∞} x(t) = 1, all perturbations in a given initial configuration die out over time (robust behavior). In contrast, if lim_{t}_{→}_{∞} x(t) ≠ 1 even small perturbations in the initial configuration propagate across the entire system and never disappear (chaotic behavior).
It has been shown in ref. 17 that for a large system (N → ∞), the temporal evolution of the overlap is given by [A.2] In the limit t → ∞ the above equation becomes a fixedpoint equation for the stationary value of the overlap x = lim_{t}_{→}_{∞}x(t). It follows from Eq. A.2 that if 2ρ(1 – ρ)K ≤ 1, the only stable fixed point is x = 1 [K is mean value of P(k)]. In this case the system is in the robust regime. On the other hand, if 2ρ(1 – ρ)K > 1 there is a stable fixed point x ≠ 1 and the dynamics of the network are chaotic. The critical value K_{c} of the mean connectivity at which the transition from robust to chaotic dynamics occurs is given as a function of ρ by the equation [A.3] The values of K_{c} and ρ that satisfy this equation are plotted on Fig. 1. It is interesting to note that the transition from robust to chaotic dynamics is governed by K, the mean value of P(k). However, this parameter is meaningful only when the distribution P(k) has a welldefined variance. In this case the fluctuations in the number of connections of the individual elements around the mean connectivity K are bounded. This is actually the case when P(k) is a Poissonian or a Gaussian distribution. For these distributions the mean connectivity K is the relevant parameter that characterizes the network topology. Originally Kauffman (11, 12) proposed this model under the assumption that K is a relevant parameter. However, the variance of the scalefree distribution P(k) = [Z(γ)k^{γ}]^{–}^{1} is infinite for γ ≤ 3. This means that the actual number of connections of the individual elements can vary from 1 up to N. Consequently, the mean value of the connectivity is no longer a meaningful parameter to characterize the network topology. For a highly heterogeneous network with scalefree topology, the only relevant parameter that characterizes the network topology is the scalefree exponent γ itself. For these kinds of networks it is better to represent the mean connectivity K as a function of the scalefree exponent γ, which gives K = Z(γ – 1)/Z(γ). In the above expression is the Riemann Zeta function. With this representation, the transition from robust to chaotic behavior is then given by [A.4] The values of γ and ρ for which this transcendental equation is satisfied are plotted on Fig. 2.
The network that we have considered is a directed graph. The value of each element σ_{i} is controlled by a set of k_{i} elements. But σ_{i} can in turn control the values of a number of other elements. Therefore, each element has a set of input connections and a set of output connections. The number of inputs and outputs of each element are not necessarily distributed with the same probability distribution P(k). However, the total number of inputs in the network equals the total number of outputs, so that on average, every element has the same number of inputs than outputs. Because the phase transition is governed only by the first moment of P(k), it follows that Eq. A.4 is valid if either the distribution of input connections or that of the output connections (or both) is scale free.
Footnotes

↵† To whom correspondence should be addressed. Email: maximino{at}control.uchicago.edu.

This paper was submitted directly (Track II) to the PNAS office.
 Received November 6, 2002.
 Copyright © 2003, The National Academy of Sciences
References
 ↵

Kitano, H. (2002) Science 295, 1662–1664.pmid:11872829
 ↵
Ideker, T., Thorsson, V., Ranish, J. A., Christmas, R., Buhler, J., Eng, J. K., Bumgarner, R., Goodlett, D. R., Aebersold, R. & Hood, L. (2001) Science 292, 929–934.pmid:11340206
 ↵
 ↵
 ↵
Yi, T. M., Huang, Y., Simon, M. I. & Doyle, J. (2000) Proc. Natl. Acad. Sci. USA 97, 4649–4653.pmid:10781070
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
Luque, B. & Sole, R. V. (1997) Phys. Rev. E Stat. Phys. Plasmas Fluids Relat. Interdiscip. Top. 55, 257–260.
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
Ravasz, E., Somera, A. L., Mongru, D. A., Oltvai, Z. N. & Barabasi, A. L. (2002) Science 297, 1551–1555.pmid:12202830
 ↵
Newman, M. E. J. (2001) Proc. Natl. Acad. Sci. USA 98, 404–409.pmid:11149952
 ↵
Sole, R. V. & Montoya, J. M. (2001) Proc. R. Soc. London Ser. B 268, 2039–2045.
 ↵
Maslov, S. & Sneppen, K. (2002) Science 296, 910–913.pmid:11988575
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵

Rzhetsky, A. & Gomez, S. M. (2001) Bioinformatics 17, 988–996.pmid:11673244

Cancho, R. F. I. & Sole, R. V. (2001) Proc. R. Soc. London Ser. B 268, 2261–2265.
 ↵