Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Phys Rev Lett. Author manuscript; available in PMC 2010 August 23.
Published in final edited form as:
PMCID: PMC2925242

Microdynamics and Criticality of Adaptive Regulatory Networks


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 [set membership] 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 S is a digraph in which each edge υi ~ υj [set membership] E additionally has a unique sign σij [set membership] [−1, +1] depending on whether it is “activating” (σij = +1) or “inhibiting” (σij = −1). The adjacency matrix A = aij of a signed digraph has the form aij = σij if υi ~ υj [set membership] E and aij = 0 otherwise. When considering structural features of S without regard for signs we shall also use the absolute adjacency matrix à = |aij|. The in-degree (out-degree) of a vertex is the number of incoming (outgoing) edges it has, without regard for sign. The net degree dneti) = |dini) − douti)| of a vertex υi is the absolute difference of its incoming and outgoing degree. Intuitively, net degree measures how “sourcelike” or “sinklike” a vertex is. By extension, we define the imbalance of a vertex pair as the absolute difference of their net degrees, Ii, υυ) = |dneti) − dnetj)|. It has recently been observed that many empirical networks contain significantly more source and sink vertices than expected by chance, and that this degree imbalance naturally leads to depletion of feedback loops which, in turn, confers enhanced stability properties [8]. Thus, degree imbalance and dynamic stability are intrinsically related, a fact that our model exploits.


We begin at t = 0 with a random signed digraph 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 S(t) at successive time steps according to the following rules. (1) Randomly and uniformly choose an edge eold = υa ~ υb connecting two vertices in S(t) such that S(t) − eold is not disconnected and an ordered pair of nonadjacent vertices υc, υd ≠ υc. (2) Calculate the pair-wise imbalances Ia, υb) and Ic, υd). (3) Delete eold and create a new edge enew = υ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 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 − hmax(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].

FIG. 1
(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 × 106 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 Dnet(t) = Σidneti(t)], a measure of overall degree imbalance in 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) = ncyc(t)/n where ncyc(t) is the number of vertices which participate in a cycle in S(t), measures overall cyclic structure without regard for details such as cycle numbers or distribution of cycle lengths. The second,


where [lambda with tilde]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 S(t). In particular, Ψ(t) is a sum of all closed walks in S(t) weighted in decreasing order by length [note that Ψ(t) + n may be thought of as the partition function of S(t)] [14].

FIG. 2
(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 S(t) is acyclic, this indicates that periods of stability occur primarily when 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 Hmax(t)] = 0.100. The mutual information between these series is Mmax(t), Φ(t)] = Mmax(t), Ψ(t)] = 0.088, indicating that changes in stability are strongly related to changes in cyclic structure [15].

FIG. 3
(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 S are disjoint (that is, each vertex υ [set membership] 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 S is cycle-disjoint, then its spectrum has a particularly simple form. Specifically, if S contains ck+ positive cycles and ck negative cycles of length k (for k = 3,…, n) [17] and all cycles are disjoint, then its spectrum is the zero eigenvalue with multiplicity (nncyc), along with the eigenvalues of each of the cycles considered separately as induced subgraphs (that is, the union of ck+ copies of the kth roots of +1, and ck copies of the kth roots of −1, for k = 3,…, n). An immediate consequence of this result is that if S is cycle-disjoint and possesses at least one positive cycle then λmax = 1, while if all cycles are negative then λmax = Reeπi/l = cos(π/l), where l is the length of the longest cycle in S. Examination of the time-series data shows that λmax(t) = 0, 1 and cos(π/l) for some 3 ≤ ln [set membership] Z+ 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 S is cycle-disjoint. In particular, if S contains ck(=ck++ck) disjoint cycles of length k for k = 3,…, n then, using Eq. (1),


where Hk,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 S(t) is acyclic, λmax(t) = 0 and μmax(t) = − minεi < 0, and the system is stable. Second, if 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).

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

For completeness it should be noted that if 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.

Supplementary Material

Supporting materials


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 × 106 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?
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.