Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2925242

Formats

Article sections

Authors

Related links

Phys Rev Lett. Author manuscript; available in PMC 2010 August 23.

Published in final edited form as:

Published online 2010 April 19. doi: 10.1103/PhysRevLett.104.168701

PMCID: PMC2925242

NIHMSID: NIHMS227652

The publisher's final edited version of this article is available at Phys Rev Lett

See other articles in PMC that cite the published article.

We present a model of adaptive regulatory networks consisting of a simple biologically motivated rewiring procedure coupled to an elementary stability criterion. The resulting networks exhibit a characteristic stationary heavy-tailed degree distribution, show complex structural microdynamics, and self-organize to a dynamically critical state. We show analytically that the observed criticality results from the formation and breaking of transient feedback loops during the adaptive process.

Much recent research attention has focused on understanding the structure of naturally occurring empirical networks and associated random graph models [1]. An overarching aim of many of these studies is to determine the relationships between network structure and dynamics. For instance, modularity and sparsity have long been known to contribute to global stability, while feedback is a well-studied prerequisite for the support of complex dynamics such as oscillations, multistability, and chaos [2]. However, although much work has, so far, focused on networks which are static in their topology, many real-world complex systems evolve both structurally and dynamically over time [3]. For instance, neural networks change in structure depending on synaptic activity while genetic regulatory networks change structurally on the evolutionary time scale in a fitness-dependent manner. Consequently, adaptive networks—in which changes in network topology and dynamics continually feedback on each other—are now attracting increasing research interest [4]. Since many biological regulatory systems, such as neural and genetic regulatory networks, are also thought to optimally balance stability and adaptability by operating at, or near to, criticality [5], a number of prior studies have sought to elucidate mechanisms of self-organized criticality (SOC) in adaptive networks [6]. For example, important early results were obtained by Christensen *et al.* and Bornholdt and Rohlf, who showed that adaptive networks may self-organize to a critical state by a simple mechanism in which “quiet nodes grow links [and] active nodes lose links” [7]. However, despite the apparent ubiquity of critical adaptive networks in nature, the mechanisms of adaptive SOC remain to be fully determined.

In this Letter we outline a simple new adaptive network model which reproduces characteristic features of biological systems, including a heavy-tailed degree distribution and self-organization to a dynamically critical state. To fix ideas our model may be thought of as describing adaptive changes in a genetic regulatory network, although it may also be applied to other systems which undergo adaptive rewiring. In genetic regulatory networks, genetic mutations cause changes in protein structure which, in turn, not only alter local network connectivity but also global system stability. Consequently, our model is a simple scheme intended to describe, albeit in an idealized way, mutation-driven local rewiring in the face of a global stability (fitness) constraint: mutations are allowed to accumulate during times of stability, but harmful mutations are suppressed during times of instability.

Mathematically, a network is a graph consisting of a set of vertices (of size *n*) and a set of edges *E*. A directed graph (digraph) is a graph in which each edge υ_{i} ~ υ_{j} *E* has a unique orientation (υ_{i} → υ_{j}). Although digraphs describe well structural relationships in complex systems, in many cases relationships also have an intrinsic sign. To cope with such systems, a natural framework is that of signed digraphs. A signed digraph *$\stackrel{\u20d7}{S}$* is a digraph in which each edge υ_{i} ~ υ_{j} *E* additionally has a unique sign σ_{ij} [−1, +1] depending on whether it is “activating” (σ_{ij} = +1) or “inhibiting” (σ_{ij} = −1). The adjacency matrix **A** = *a _{ij}* of a signed digraph has the form

We begin at *t* = 0 with a random signed digraph *$\stackrel{\u20d7}{S}$*(*t* = 0) with Erdős-Rényi connectivity, in which edge orientations and signs have been assigned independently in an equiprobable random manner [9]. We then rewire *$\stackrel{\u20d7}{S}$*(*t*) at successive time steps according to the following rules. (1) Randomly and uniformly choose an edge *e*_{old} = υ_{a} ~ υ_{b} connecting two vertices in *$\stackrel{\u20d7}{S}$*(*t*) such that *$\stackrel{\u20d7}{S}$*(*t*) − *e*_{old} is not disconnected and an ordered pair of nonadjacent vertices υ_{c}, υ_{d} ≠ υ_{c}. (2) Calculate the pair-wise imbalances *I*(υ_{a}, υ_{b}) and *I*(υ_{c}, υ_{d}). (3) Delete *e*_{old} and create a new edge *e*_{new} = υ_{c} → υ_{d}, choosing its sign randomly and uniformly, and recalculate the imbalances. (4) If the sum of the two imbalances after the switch is greater than that before then accept the switch unconditionally, otherwise accept with probability ρ(*t*).

In order to couple structural rearrangement to dynamics we allow ρ(*t*) to vary in a manner which takes into account the changing stability of the system. To do so we assume that, in addition to regulatory links defined by *$\stackrel{\u20d7}{S}$*(*t*), each species (vertex υ_{i}) also decays at a constant characteristic rate ε_{i} which we fix at *t* = 0 independently, randomly, and uniformly on the unit interval. Thus, at each evolutionary time point we obtain a modified adjacency matrix **B**(*t*) = **A**(*t*) − diag(ε_{i}). Global stability is then given by the magnitude of μ_{max}(*t*) = max Re μ_{i}(*t*), where μ_{i}(*t*) for *i* = 1,…, *n* are the eigenvalues of **B**(*t*). In particular, the system is stable when μ_{max}(*t*) < 0 and unstable when μ_{max}(*t*) > 0 [10]. Therefore, we set ρ(*t*) = 1 − *h*[μ_{max}(*t*)], where *h*[*x*] is the Heaviside step function, allowing defective switches when the system is stable and suppressing defective switches when the system is unstable. The key property of this coupling is that it makes global information available to the local structural reorganizing process, providing continual feedback between structure and dynamics.

The networks produced by this simple model are characterized by a stationary heavy-tailed degree distribution (see Fig. 1, left) indicating the presence of hub source and sink vertices, a well-known feature of real-world networks [8,11]. However, since our model allows for periods of random structural rearrangement, this macroscopic stationarity masks complex structural *microdynamics* in which individual vertices continually accumulate and lose edges (see Fig. 1, right). This kind of “mixing” microdynamics is not produced by classical rich-get-richer models of hub formation [11], but has recently been highlighted as an important characteristic of real-world evolving networks [3].

(color online). Graphs resulting from the evolutionary process exhibit a stationary heavy-tailed degree distribution and complex microdynamics. Left: The net degree distribution at 1000 time-step intervals for a period of 1 × 10^{6} time steps at **...**

Figure 2 gives a plot of μ_{max}(*t*) at equilibrium [12] for a representative system showing that dynamics on the evolutionary time scale are characterized by periods of stability [μ_{max}(*t*) < 0] and instability [μ_{max}(*t*) > 0] punctuating back and forth. To help interpret these dynamics, also shown is λ_{max}(*t*) = max Re λ_{i}(*t*), where λ_{i}(*t*) are the eigenvalues of the graph adjacency matrix **A**(*t*) and three measures of network structure. The first structural measure shown is total net degree *D*_{net}(*t*) = Σ_{i}*d*_{net}[υ_{i}(*t*)], a measure of overall degree imbalance in *$\stackrel{\u20d7}{S}$*(*t*). It is apparent that changes in total net degree correlate poorly with changes in stability, suggesting that although fluctuations in net degree are observed during the evolutionary process, it is not degree imbalance *per se* that drives the characteristic dynamics of μ_{max}(*t*). In order to identify more precisely the structural origin of the observed bursting dynamics, and based upon the observation that degree imbalance naturally leads to feedback loop depletion [8], also shown are two measures of network cyclic structure [13]. The first, Φ(*t*) = *n*_{cyc}(*t*)/*n* where *n*_{cyc}(*t*) is the number of vertices which participate in a cycle in *$\stackrel{\u20d7}{S}$*(*t*), measures overall cyclic structure without regard for details such as cycle numbers or distribution of cycle lengths. The second,

$$\mathrm{\Psi}(t)=\text{Tr}{e}^{\mathbf{\xc3}(t)}-n={\displaystyle \sum _{i=1}^{n}{e}^{{\tilde{\lambda}}_{i}(t)}-n,}$$

(1)

where _{i} are the eigenvalues of the absolute adjacency matrix **Ã**(*t*), is an indirect measure of “returnability,” which takes into account details of closed walks in *$\stackrel{\u20d7}{S}$*(*t*). In particular, Ψ(*t*) is a sum of all closed walks in *$\stackrel{\u20d7}{S}$*(*t*) weighted in decreasing order by length [note that Ψ(*t*) + *n* may be thought of as the partition function of *$\stackrel{\u20d7}{S}$*(*t*)] [14].

(color online). Evolutionary dynamics are characterized by punctuated equilibrium. Dynamics of a representative 100 vertex network are given. Panel (a) shows the time series for μ_{max}(*t*). Note that μ_{max}(*t*) = − minε_{i} < **...**

Examining the time series of Ψ(*t*) and Φ(*t*) it is apparent that, unlike total net degree, both Ψ(*t*) and Φ(*t*) exhibit similar bursting behavior to that of μ_{max}(*t*). In particular, periods of stability [μ_{max}(*t*) < 0] generally correspond to periods when both Ψ(*t*) = 0 and Φ(*t*) = 0 (in Fig. 2 this occurs >90% of the time). Since Ψ(*t*) = 0 and Φ(*t*) = 0 if and only if *$\stackrel{\u20d7}{S}$*(*t*) is acyclic, this indicates that periods of stability occur primarily when *$\stackrel{\u20d7}{S}$*(*t*) is acyclic. Furthermore, changes in stability predominantly occur concordantly with changes in Ψ(*t*) and Φ(*t*) (for instance, in Figs. 2 and and33 this occurs >99% of the time). Considering the time series as binary variables (“stable or unstable” and “cyclic or acyclic”) and calculating entropies gives *H*[Φ(*t*)] = *H*[Ψ(*t*)] = 0.093 and *H*[μ_{max}(*t*)] = 0.100. The mutual information between these series is *M*[μ_{max}(*t*), Φ(*t*)] = *M*[μ_{max}(*t*), Ψ(*t*)] = 0.088, indicating that changes in stability are strongly related to changes in cyclic structure [15].

(color online). Transient cycles trigger bursts of instability. Left: A plot of ΔΨ(*t*) = Ψ(*t*) − Ψ(*t* − 1) against Δμ_{max}(*t*) = μ_{max}(*t*) − μ_{max}(*t* − 1) using the **...**

In order to better understand this relationship, we now derive some analytical results relating cycles and spectra of signed digraphs which will help interpret these numerics. To obtain exact results we shall focus on deriving analytic formulas for Ψ and λ_{max} in the particular case that all cycles in *$\stackrel{\u20d7}{S}$* are disjoint (that is, each vertex υ *V*; belongs to at most one cycle). Although this is a strong condition to impose, and most real-world networks are not expected to be cycle-disjoint, this case is analytically tractable and, since our evolutionary scheme favors the minimization of cycles, yields results which shed light on the observed dynamics. Full proofs of all analytic results are provided in the supplementary materials [16].

First we observe that if a signed digraph *$\stackrel{\u20d7}{S}$* is cycle-disjoint, then its spectrum has a particularly simple form. Specifically, if *$\stackrel{\u20d7}{S}$* contains
${c}_{k}^{+}$ positive cycles and
${c}_{k}^{-}$ negative cycles of length *k* (for *k* = 3,…, *n*) [17] and all cycles are disjoint, then its spectrum is the zero eigenvalue with multiplicity (*n* − *n*_{cyc}), along with the eigenvalues of each of the cycles considered separately as induced subgraphs (that is, the union of
${c}_{k}^{+}$ copies of the *k*th roots of +1, and
${c}_{k}^{-}$ copies of the *k*th roots of −1, for *k* = 3,…, *n*). An immediate consequence of this result is that if *$\stackrel{\u20d7}{S}$* is cycle-disjoint and possesses at least one positive cycle then λ_{max} = 1, while if all cycles are negative then λ_{max} = Re*e*^{πi/l} = cos(π/*l*), where *l* is the length of the longest cycle in *$\stackrel{\u20d7}{S}$*. Examination of the time-series data shows that λ_{max}(*t*) = 0, 1 and cos(π/*l*) for some 3 ≤ *l* ≤ *n* ^{+} do indeed occur commonly during evolution (for instance, in Fig. 2 this occurs ≈ 53% of the time), indicating the continual formation and breaking of isolated cycles by the evolutionary scheme.

This result is also useful since it allows us to calculate Ψ analytically in the case that *$\stackrel{\u20d7}{S}$* is cycle-disjoint. In particular, if *$\stackrel{\u20d7}{S}$* contains
${c}_{k}\phantom{\rule{thinmathspace}{0ex}}(={c}_{k}^{+}+{c}_{k}^{-})$ disjoint cycles of length *k* for *k* = 3,…, *n* then, using Eq. (1),

$$\mathrm{\Psi}={\displaystyle \sum _{k=3}^{n}{c}_{k}{H}_{k,0}(1)-{n}_{\text{cyc}},}$$

(2)

where *H*_{k,0}(*z*) is the generalized hyperbolic function of order *k* and kind 0 [18]. Figure 3 shows that values of Ψ(*t*) calculated using Eq. (2) often occur during evolution, again indicating that isolated cycles are continually formed and broken by the evolutionary scheme.

These analytical results may be used to interpret numerics by making use of two further results which relate λ_{max} to μ_{max} in the cycle-disjoint case. First, note that in the special case that *$\stackrel{\u20d7}{S}$*(*t*) is acyclic, λ_{max}(*t*) = 0 and μ_{max}(*t*) = − minε_{i} < 0, and the system is stable. Second, if *$\stackrel{\u20d7}{S}$*(*t*) is cycle-disjoint then μ_{max}(*t*) < λ_{max}(*t*) and this bound is tight [μ_{max}(*t*) → λ_{max}(*t*) as ε_{i} → 0 for all *i*]. Consequently, if vertex decay rates are all small then the completion of a single cycle in an otherwise acyclic network is sufficient to trigger a burst of instability, as seen in Fig. 2. When this occurs the evolutionary process responds by suppressing any further defective switches and rearranging local network structure to remove the cause of the instability. Typically, this is quickly achieved and the burst of instability is relatively short. However, occasionally cycles may accumulate more rapidly than they are removed, giving rise to extended bursts of instability and heavy-tailed statistics characteristic of a critical state (see Fig. 4).

(color online). Bursts of instability have heavy-tailed statistics. The distribution of burst durations is shown over an interval of 4 × 10^{6} evolutionary time steps for a representative system with 100 vertices at equilibrium. A power law with **...**

For completeness it should be noted that if *$\stackrel{\u20d7}{S}$*(*t*) is not cycle-disjoint then the relationship between cycles and stability can be complex: it is not necessarily true that μ_{max}(*t*) < λ_{max}(*t*) and, in rare cases, changes in cyclic structure and stability may occur discordantly. Further details are included in the supplementary materials [16].

Many biological regulatory systems are thought to balance stability and adaptability by self-organizing to a dynamically critical state [5]. In this Letter we have presented a simple adaptive network model which reproduces characteristic features of biological systems, including a heavy-tailed connectivity distribution, microdynamics, and robust self-organization to criticality. Previous models have shown that adaptive networks may self-organize to a critical state due to rewiring based upon local activity [6,7]. Here, the mechanism of self-organization is somewhat different and relies on the fact that feedback and stability are generally inversely related: by employing a flexible rewiring scheme which allows feedback loops to be formed during periods of stability and eliminated during periods of instability, criticality naturally arises in our model. It seems plausible that these (and other, as yet unknown) adaptive processes may be responsible for the criticality observed in nature.

PACS numbers: 89.75.Fb, 02.10.Ox, 05.65.+b, 89.75.Hc

1. Albert R, Barabási A-L. Rev. Mod. Phys. 2002;74:47.Newman MEJ. SIAM Rev. 2003;45:167.

2. May RM. Nature (London) 1972;238:413. [PubMed]Thomas R, D’Ari R. Biological Feedback. Boca Raton, FL: CRC Press; 1990.

3. Batty M. Nature (London) 2006;444:592. [PubMed]Luscombe NM, Madan Babu M, Yu H, Snyder M, Teichmann SA, Gerstein M. ibid. 2004;431:308. [PubMed]Gautreau A, Barrat A, Barthélemy M. Proc. Natl. Acad. Sci. U.S.A. 2009;106:8847. [PubMed]

4. Gross T, Blasius B. J. R. Soc. Interface. 2008;5:259. [PubMed]

5. Beggs JM, Plenz D. J. Neurosci. 2003;23:11 167. [PubMed]Shmulevich I, Kauffman SA, Aldana M. Proc. Natl. Acad. Sci. U.S.A. 2005;102:13 439, 11 167. [PubMed]Nykter M, Price ND, Aldana M, Ramsey SA, Kauffman SA, Hood LE, Yli-Harja O, Shmulevich I. ibid. 2008;105:1897. [PubMed]Balleza E, Alvarez-Buylla ER, Chaos A, Kauffman S, Shmulevich I, Aldana M. PLoS ONE. 2008;3:e2456. [PubMed]

6. Bornholdt S, Sneppen K. Phys. Rev. Lett. 1998;81:236.Luque B, Ballesteros FJ, Muro EM. Phys. Rev. E. 2001;63:051913. [PubMed]Bornholdt S, Röhl T. ibid. 2003;67:066118. [PubMed]Liu M, Bassler KE. ibid. 2006;74:041910. [PubMed]Garlaschelli D, Capocci A, Caldarelli G. Nature Phys. 2007;3:813.Rohlf T. Europhys. Lett. 2008;84:10 004.Magnasco MO, Piro O, Cecchi GA. Phys. Rev. Lett. 2009;102:258102. [PubMed]Meisel C, Gross T. Phys. Rev. E. 2009;80:061917. [PubMed]

7. Christensen K, Donangelo R, Koiller B, Sneppen K. Phys. Rev. Lett. 1998;81:2380.Bornholdt S, Rohlf T. ibid. 2000;84:6114. [PubMed]

8. Ma’ayan A, Cecchi GA, Wagner J, Rao AR, Iyengar R, Stolovitzky G. Proc. Natl. Acad. Sci. U.S.A. 2008;105:19235. [PubMed]

9.
For efficiency in the simulations shown we set the edge-inclusion probability ≈ ln(*n*)/*n* and consider a maximally sparse connected random graph. Qualitatively the same results may be achieved for more dense graphs.

10.
This stability criterion assumes that **B**(*t*) is the Jacobian matrix of a dynamical system evaluated at a fixed point.

11. Barabási A-L, Albert R. Science. 1999;286:509. [PubMed]

12.
That is, after an initial transient “settling-down” period (4 × 10^{6} time steps prior to data shown).

13.
A cycle of length *k* is a nonintersecting path of length *k* from a vertex back to itself respecting edge directions.

14. Estrada E, Hatano N. Chem. Phys. Lett. 2007;439:247.Linear Algebra Appl. 2009;430:1886.

15.
By discretizing the data we are asking, How much does knowing whether the network is cyclic or not tell us about whether the system is stable or not?

16.
See supplementary material at http://link.aps.org/supplemental/10.1103/PhysRevLett.104.168701 for proofs and examples.

17.
A cycle *c* is positive (negative) if the product of the edge signs in *c* equals +1 (−1).

18. Ungar A. Am. Math. Mon. 1982;89:688.

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |