Stability as a natural selection mechanism on interacting networks

 

Juan I. Perotti,¹' ²* Orlando V. Billoni,¹' Francisco A. Tamarit,¹' ², Sergio A. Cannas¹' ²

Email: juanpool@gmail.com
†Email: billoni@famaf.unc.edu.ar
‡Email: tamarit@famaf.unc.edu.ar
§Email: cannas@famaf.unc.edu.ar

Facultad de Matem´atica, Astronomia y Fisica, Universidad Nacional de C´ordoba, Argentina.

Instituto de Fisica Enrique Gaviola (IFEGCONICET), Ciudad Universitaria, 5000 Cordoba, Argentina.

 

Biological networks oí interacting agents exhibit similar topological properties for a wide range of scales, from cellular to ecological levéis, suggesting the existence of a common evolutionary origin. A general evolutionary mechanism based on global stability has been proposed recently [J I Perotti, et al., Phys. Rev. Lett. 103, 108701 (2009)]. This mechanism was incorporated into a model of a growing network of interacting agents in which each new agent's membership in the network is determined by the agent's effect on the network's global stability. In this work, we analyze different quantities that characterize the topology of the emerging networks, such as global connectivity, clustering and average nearest neighbors degree, showing that they reproduce scaling behaviors frequently observed in several biological systems. The influence of the stability selection mechanism on the dynamics associated to the resulting network, as well as the interplay between some topological and functional features are also analyzed.

 

These properties are observed for a wide range of scales, from the microscopic level of genetic, metabolic and proteins networks to the macroscopic level of communities of living beings (ecological networks). Such ubiquity suggests the existence of some natural selection process that promotes the development of those particular structures [3]. One possible constraint general enough to act across such a range of scales is the proper stability of the underlying dynamics.

Growing biological networks involve the coupling of at least two dynamical processes. The ñrst one concerns the addition of new nodes, attached during a slow evolutionary (i.e., species lifetime) process. A second one is the node dynamics which affects and in turn is affected by the growing processes. It is reasonable to expect that the network topologies we finally witness could have emerged out of these coupled processes. Consider, for example, the case of an ecological network like a food web, where nodes are species within an ecosystem and edges are consumerresource relationships between them. New nodes are added during evolutionary time scales, through speciation or migration of new species. Then, the network grows through community assembly rules, strongly influenced by the underlying dynamics of species and speciñc interactions among them [15, 16]. The consequence of adding a new member with a given connectivity affecting a global in/stability, is represented in this case by the aboundance/lack of food [17]. Notice that each new member may not only result in its own addition/rejection to the system, but it can also promote avalanches of extinctions amongst existing members.

The above ingredients were recently incorporated into a simple model of growing networks under stability constraints [19]. Numerical simulations on this model showed that, indeed, complex topology can emerge out of a stability selection pressure. In the present work, we further explore different topological and dynamical properties predicted by the model, whose deñnition is reviewed in section II. The results are presented in sections III. and IV. In section III., we analyze the topological features that emerge in growing networks under stability constraint. In section IV., we show that this constraint not only induces topological features of the resulting networks but also influences the associated dynamics. A discussion of the results is presented in section V.

II. The Model

Let us consider a system of n interactive agents, whose dynamics is given by a set of differential equations dx/dt = F(x), where x is an ncomponent vector describing the relevant state variables of each agent and F is an arbitrary nonlinear function. One could imagine that x in different systems may represent concentrations of some hormones, the average density populations in a food web, the concentration of chemicals in a biochemical network, or the activity of genes in a gene regulation net, etc. We assume that a given agent i interacts only with a limited set of ki < n other agents; thus, Fi depends only on the variables belonging to that set. This defines the interaction network.

We assume that there are two time scales in the dynamics. Let f_m be the average frequency of the incoming flux of new agents (migration, mutation, etc.). This defines a characteristic time r_m = J"¹. On the long time scale í »r_m (much larger than the observation time) new agents arrive to the system and start to interact with some of the previous ones. Some of them can be incorporated into the system or not, so n (and the whole set of differential equations) can change. Once a new agent starts to interact with the system, we will assume that the enlarged system evolves towards some stationary state with characteristic relaxation time T_re¡ <C r_m. Then, in the short time scale T_re¡ < í « r_m we can assume that n is constant and the dynamics already led the system to a particular stable stationary state x* deñned by F(x*) = 0. Following May's ideas [18], we assume that the only attractors of the dynamics are ñxed points. Nevertheless, the proposed mechanism is expected to work, as well, for more complex attractors (e.g, limit cycles).

The stability of the solution x* is determined by the eigenvalue with máximum real part of the Jacobian matrix

A new agent will be incorporated to the network if its inclusión results in a new stable ñxed point, that is, if the valúes of the interaction matrix ajeare such that the eigenvalue with máximum real part A_m of the enlarged Jacobian matrix is negative (A_m < 0). Assuming that isolated agents will reach stable states by themselves after certain characteristic relaxation time, the diagonal elements of the matrix a¿ ¿ are negative and given unity valué to further simplify the treatment [18]. The interaction valúes, (i.e., the nondiagonal matrix elements aij) will take random valúes (both positive and negative) taken from some statistical distribution. In this way, we have an unbounded enseñable of systems [18] characterized by a "growing through stability" history. Randomness would be selfgenerated through the addition of new agents processes. Each speciñc set of matrix elements, after addition, defines a particular dynamical system and the subsequent analysis for time scales between successive migrations is purely deterministic.

The model is then defined by the following algoritmia [19]. At every step, the network can either grow or shrink. In each step, an attempt is made to add a new node to the existing network, starting from a single agent (n = 1). Based on the stability criteria already discussed, the attempt can be successful or not. If successful, the agent is accepted, so the existing n x n matrix grows its size by one column and one row. Otherwise, the nóvate agent will have a probability to be deleted together with some other nodes, as further explained below.

More specifically, suppose that we have an already created network with n nodes, such that the n x n associated interaction matrix a¿j is stable. Then, for the attachment of the (n + í)th node, we first choose its degree k_n+i randomly between 1 and n with equal probability. Then, the new agent interaction with the existing network member i is chosen, such that nondiagonal matrix elements (C'i,n+i, 0"n+i,i) (* = 1, . . ., «) are zero with probability 1 k_n+i/n and different from zero with probability k_n+i/n] to each nonzero matrix element we assign a different real random valué uniformly distributed in [6,6]. 6 determines the interaction range variability and it is one of the two parameters of the model [20].

Then, we calcúlate numerically A_m for the resulting (n + 1) x (n + 1) matrix. If A_m < 0, the new node is accepted. If A_m > 0, it means that the introduction of the new node destabilized the entire system and we will impose that, the new agent is either eliminated or it remains but produces the extinction of a certain number of previous existing agents. In order to further simplify the numerical treatment, we allow up to q < k_n+i extinctions, taken from the set of k_n+i nodes connected to the new one; q is the other parameter of the model. To choose which nodes are to be eliminated, we first select one with equal probability in the set of k_n+1 and remove it. If the resulting n x n matrix is stable, we start a new trial; otherwise, another node among the remaining k_n+i 1 is chosen and removed, repeating the previous procedure. If after q remováis the matrix remains unstable, the new node is removed (we return to the original n x n matrix and start a new trial). The process is repeated until the network reaches a máximum size n = n_max (typically n_max = 200) and restarted M times from n = 1 to obtain statistics of the networks (typically M = 10⁵).

III. Topological properties

i. Connectivity

First, we analyzed the average connectivity C(n), defined as the fraction of nondiagonal matrix elements different from zero, averaged over different runs. In Fig. 1, we show the typical behavior of C{n) for different valúes of b (we found that C{n) is completely independent of q). The connectivity presents a power law tail for large valúes of n. Prom a ñtting of the tail with a power law (see insets in Fig. 1) we obtain the scaling behavior

for large valúes of n, where a is the variance of the nondiagonal elements of the stability matrix (a = 5²/3 for the uniform distribution) and u> = 0.7±0.1. Prom the inset of Fig. 1, we see that the exponent e shows a weak dependency on 5, taking valúes in the range (0.1,0.3) . It is interesting to compare Eq. (2) with May's stability line for random networks [18] C(n) = (an)¹. It is easy to see that Eq. (2) lies above May's stability line for network sizes up to ~ 10⁶ [21]. This shows that networks growing under stability constraint develop particular structures whose probability in a completely random ensemble is almost zero. In other words, the associated matrices belong to a subset of the random ensemble with zero measure and therefore they are only attainable through a constrained development process. In the next sections, we explore the characteristics of those networks.

In Fig. 2, we plotted the connectivity for different biological networks across three orders of magnitude of network size scales, using data collected from the literature. We see that the data are very well ñtted by a single power law C(n) ~ n¹'², in a nice agreement with the average valué e = 0.2 predicted by the present model. It is worth mentioning that the behavior C(n) ~ n(^1+e> has also been obtained in a self organized criticality model of Food Webs [26].

ii. Degree distribution

The degree distribution P(k) of the network was analyzed in detail in Ref. [19]. We briefly summarize the main results here. In Fig. 3, we illustrate the typical behavior oí P(k). It presents a power law tail P(k) ~ k '' for valúes of k > 20, with a ñnite size drop at k = n_max. The degree exponent 7 takes valúes between 2 and 3 for valúes of b in the interval 5 G (1.5, 3.5), which become almost independent of q as it increases. The exponent 7 can also fall below 2 when the global stability constraint

 

Figure 1: Connectivity as a function of the network size for q = 3, n_max = 200 and different valúes of b. The symbols correspond to numerical simulations and the dashed lines to power law ñttings of the tails C(n) = Bn'^1+£l The insets show the ñtting valúes B and e as a function of b

is replaced by a local one. The qualitative structure of P(k) remains when the stability criterium A_m < 0 is relaxed by the condition A_m < A, with A some small positive number. In other words, the power law tail emerges also when the addition of new nodes destabilizes the dynamics, provided that the characteristic time to leave the ñxed point t = A¹ is large enough to become comparable to the migration time scale r_m [19].

iii. Network growth and clustering properties

Networks grown under stability constraint also display small world properties. The average clustering coefñcient decays with the network size as Cc(n) ~ no.75 (wjjjcjj js glower than the í/n decay in a random net), while the average path length L between two nodes increases as L(n) ~ A ln (n + C) [19]. A similar behavior is observed in the BarabasiAlbert model [1], where the clustering can be approximated by a power law with the same exponent, although the exact scaling is [27] Cc(n) ~ (lnn)²/n (therefore that behavior cannot be excluded in the present model). While this suggests the presence of an underlying preferential attachment rule mechanism, a detailed analysis has shown that this is not the dominant mechanism [19]. The behavior of Ce and L is linked with the selection dynamics ruling which node is accepted or rejected. The stability constraint favors the nodes with few links, since they modify the matrix ai¿ stability much less than new nodes with many links (of course this is reflected in the P(k) density). Thus, most frequently the network grows at the expense of adding nodes with one or few links, producing an increase of L and a decrease of Ce, but sporadically, a highly connected node is accepted, decreasing L and increasing Cc(n) [19]. Those fluctuations lead to a slow diffusivelike growth of the network size n(t) ~ t¹'² (See Fig. 4).

Another quantity of interest is the average clustering Cc(k) as a function of the degree k. A typical example is shown in Fig. 5. We see that Cc(k) decreases monotonously with k and displays a power law tail Cc(k) ~ k@ with an exponent ¡3 « 0.9, cióse to one. The exponent appears to be completely independent of b and q. This behavior is indicative of a modular structure with hierarchical organization [13]. Notice that this power law decay appears for degrees k > 20, precisely the same range of valúes for which the degree distribution P(k) displays a power law tail (see subsection ii.).

iv. Mixing by degree patterns

To analyze the mixing by degree properties of the networks selected by the stability constraint, we calculated the average degree k_nn among the nearest neighbors of a node with degree k. In Fig. 6, we see that k_nn decays with a power law k_nn ~ kfor k > 20, with an exponent S cióse to 0.25, in a clear disassortative behavior. This result is also consistent with previous works showing that assortative mixing by degree decreases the stability of a network, Le., the máximum real part A_m of the eigenvalues of random matrices of the type here considered increases faster on assortative networks than on disassortative ones [29].

IV. Dynamical properties

In the previous section, we analyzed different topological properties that are selected by the stability constraint, Le., properties associated to the underlying adjacency matrix, regardless of the valúes of the interaction strengths. We now analyze the characteristics of the dynamics associated to the networks emerging from such constraint. In other words, we investígate the statistics of valúes of the non nuil elements a¿j = 0.

First of all, we calcúlate the probability distribution of valúes for a single non nuil matrix element a,ij of the ñnal network with size n = n_max. The typical behavior is shown in Fig. 7. We see that Piflij) is an even function, almost uniform in the interval [6,6], with a small cusp around a¿¿ = 0. This shows that stability is not enhanced by a particular sign or absolute valué of the individual interaction coefñcients. It has been shown recently that the presence of anticorrelated links between pairs of nodes (Le., links between pairs of nodes (i,j) such that signiciij) = sign(aji)) signiñcatively enhances the stability of random matrices [31]. In an ecological network, this typically corresponds to a predatorprey or parasitehost interaction. To check for the presence of such type of interactions, we calculated the correlation (a¿jaj¿), where the average is taken over pairs of nodes with a double link (a¿j = 0 and a^ = 0).

In Fig. 8, we show (a,ijaji) as a function of the network size n. We see that this correlation is negative for any valué of n and saturates into a valué ¿) ra 0.65 for large valúes of n. In the inset of Fig. 8, we compare the average fraction of double links (rj) with the corresponding quantity for a completely random network with the same connectivity C(n), that is, a network where all edges are independently distributed with probabil

Then for large valúes of n, we have (rj)_ran ~ C(n) ~ n^(1+e). From the inset of Fig. 8, we see that (rj) ~ n⁰'⁶⁸ when íi > 1 in the present case. The fraction of double links is considerable larger than in a random network. The two results of Fig. 8 together show that the present networks have indeed a signiñcantly large number of anticorrelated pair interactions.

Next, we calculated the correlation (o>íjo>jí) / (\o>ij\)² between the matrix elements, linking a node i and its neighbors j, as a function of its degree k_iy where the average is taken only on the double links. From Fig. 9, we see that the absolute valué of the correlation presents a máximum around ki = 25 and tends to zero as the degree increases. The inset of Fig. 9 shows that the average fraction of anticorrelated links (k) (i.e., # anticorrelated links/total =$= double links) tends to 1/2 as the degree increases. We can conclude from these results that the interactions strengths between the hubs and their neighbors are almost uncorrelated. This suggests that the influence of hubs in the stabilization of the dynamics is mainly associated to their topological role (e.g., reduction of the average length L) rather than to the nature of their associated interactions.

V. Discussion

The recent advances in the research on networks theory in biological systems have called for a deeper understanding about the relationship between network structure and function, based on evolutionary grounds [3]. In this work, we have shown that a key factor to explain the emergeney of many of the complex topological features commonly observed in biological networks could be just the stability of the underlying dynamics. Stability can then be considered as an effective ñtness acting in all biological situations. The results presented in Fig. 2 for the connectivity of real biological networks at different network size scales support this conclusión. In addition, the present approach (although based on a very simple model) allows to draw some conclusions about the interplay between network structure and function that could be of general applicability. The present results suggest that hubs play mainly a topological role of linking modules (disassortativity low clustering, uncorrelated links), while low connected nodes inside modules enhance stability through the presence of many anticorrelated interactions.

The stabilizing effects of some of the topological and functional network features here analyzed have been previously addressed separately (small world [32], dissasortative mixing [29, 30], anticorrelated interactions [31]). However, the present analysis suggests that the simultaneous observance of all of them is highly unlikely to be a result of a purely random process. Such delicate balance of speciñc topological and functional features would only be attainable through a slow, evolutionary stability selection process.

In particular, the above scenario agrees very well with the observed structures in cellular networks. For instance, the scaling behavior of Cc(k), displayed in Fig. 5, has been observed in metabolic [28] and protein [6, 7] networks. Disassortative

mixing by degree is another ubiquous property of those systems and indeed a very similar behavior to that shown in Fig. 6 has been observed in certain proteinprotein interaction networks [7]. Also, the available data for the degree distribution in all those cases are consistent with a power law behavior with an exponent 7 between 2 and 2.5 [1]. The agreement with the whole set of properties predicted by the model suggests that stability could be a key evolutionary factor in the development of cellular networks.

The situation is a bit different in the case of ecological networks, where the predictions of the model do not completely agree with the observations, specially those related to food webs. On the one hand, food webs usually display also disassortative mixing by degree, modularity and relatively low small worldness [33] (rather low valúes of clustering, compared with other biological networks), in agreement with the present predictions. Regarding the scaling behavior of C{n) [24], this is the topic of an oíd debate in ecology (see Dunne's review in Ref. [23] for a summary of the debate). While in general it is expected a power law behavior, the valué of the exponent (and the associated interpretations) is controversial, due to the large dispersión of the available data, the rather small range of network sizes available and, in some cases, the low resolution of the data [25]. The consistency of the scaling shown in Fig. 2 for a broad range of size scales suggests that the ecological debate should be reconsidered in a broader context of evolutionary growth under stability constraints.

On the other hand, the degree distribution of food webs is not always a power law, but it frequently exhibits an exponential cutoff at some máximum characteristic degree k_max [9, 10]. Such varianee between food webs and other biological networks is probably related to the way ecosystems assemble and evolve compared with other systems. While the hypothesis of the present model are general enough to apply in principie to any biological system, that difference suggests that stability is not enough to explain the observed structures in food webs, but further constraints should be included to account for them. For instance, at least two different (although closely related) constraints are known that can genérate a cutoff in the degree distribution: aging and limited capacity of the nodes [34]. In the former, nodes can become inactive with some probability through time (in the sense that they stop interacting with new agents), while in the latter they systematically pay a "cost" every time a new link is established with them, so that they become inactive when some maximum degree is reached. One can easily imagine different situations in which mechanisms of that type become important in the evolution of ecological webs, ther by limitations in the available resources or by dynamical changes in the diet of species due to external perturbations (for instance, there are many factors that constrain a predator's diet; see Ref. [9] and references therein for a related discussion). Mechanisms of these kind can be easily incorporated into the model, serving as a base for the description of more complex behaviors in particular systems like food webs.

Finally, it would be interesting to analyze the relationship between dynamical stability in evolving complex networks and synchronization, a topic about which closely related results have been recently published [35].

 

Acknowledgements This work was supported by CONICET, Universidad Nacional de C´ordoba, and FONCyT grant PICT2005 33305 (Argentina). We thank useful discussions with P. Gleiser and D. R. Chialvo. We acknowledge useful comments and criticisms of the referee.

[1] R Albert, AL Barab´asi, Statistical mechanics of complex networks, Rev. Mod. Phys. 74, 47 (2002).

[2] M E J Newman, The Structure and Function of Complex Networks, SIAM Review 45, 167 (2003).

[3] S R Proulx, D E L Promislow, P C Phillips, Network thinking in ecology and evolution, Trends Ecol. Evol. 20, 345 (2005).

[4] S Boccaletti,   V Latora, Y Moreno, M Chavez, D U  Hwang, Complex networks: Structure and  dynamics, Phys. Rep. 424, 175 (2006).

[5] N A Alves, A S Martinez, Inferring topological features of proteins from amino acid residue networks, Physica A 375, 336 (2007).

[6] S H Yook, Z N Oltvai, AL Barab´asi, Functional and topological characterization of protein interaction networks, Proteomics 4, 928 (2004).

[7] V Colizza, A Flammini, A Maritan, A Vespignani, Characterization and modeling of proteinprotein interaction networks, Physica A 352, 1 (2005).

[8] H Jeong, B Tombor, R Albert, Z N Oltvai, AL Barab´asi, The largescale organization of metabolic networks, Nature 407, 651 (2000).

[9] J M Montoya, S L Pimm, R V Sol´e, Ecological networks and their fragility, Nature 442, 259 (2006).

[10] J A Dunne, R J Williams, N D Martinez, Food web structure and network theory: The role of connectance and size, P. Natl. Acad. Sci. USA 99, 12917 (2002).

[11] O Sporns, D R Chialvo, M Kaiser, C C Hilgetag, Organization, development and function of complex brain networks, Trends Cogn. Sci. 8, 418 (2004).

[12] V M Eguiluz, D R Chialvo, G A Cecchi, M Baliki, A V Apkarian, ScaleFree Brain Functional Networks, Phys. Rev. Lett. 94, 018102 (2005).

[13] E Ravasz, AL Barab´asi, Hierarchical organization in complex networks, Phys. Rev. E 67, 026112 (2003).

[14] M E J Newman, Mixing patterns in networks, Phys. Rev. E 67, 026126 (2003).

[15] E Weiher, P Keddy, Ecological Assembly Rules, Cambridge U. Press, Cambridge (1999).

[16] S L Pimm, The balance of Nature, The University of Chicago Press, Chicago (1991).

[17] In particular cases, like foodwebs, there could be other coupled process almost probably as important as the two considered here, such as the dynamical change in previously existing links weights. For instance, if the abundace of one species is reduced, this could trigger the adaptation of predators to feed on other existing species, thus creating new links not directly related to the incoming species. The present scheme should be considered as a minimal model about stability as an evolutionary constraint.

[18] R M May, Will a large complex system be stable?, Nature 238, 413 (1972).

[19] J I Perotti, O V Billoni, F A Tamarit, D R Chialvo, S A Cannas, Emergent SelfOrganized Complex Network Topology out of Stability Constraints, Phys. Rev. Lett. 103, 108701 (2009).

[20] Some of the numerical calculations performed in section III. were repeated replacing the uniform distribution for the non diagonal elements by a Gaussian one with the same variance. The results were indistinguishable.

[21] An extrapolation of the present results suggest that the probability of growing goes to zero above some maximum network size. Preliminary estimations give a value 10⁵ for such maximum size and therefore all the obatined networks would be above May's stability line. This topic is the subject of present investigations and the corresponding results will be published in the near future.

[22] C Melia, J Bascompte, Food web cohesion, Ecology 85, 352 (2004).

[23] J A Dunne, Food webs, In: Encyclopedia of Complexity and Systems Science, Ed. R. A. Meyers, pag. 3661, Springer, New York (2009).

[24] The quantity usually considered in ecology is the connectance, defined as the number of links in the web divided by n². It is proportional to the connectivity C(n), with a proportionality factor ≈ 1/2.

[25] For instance, a power law fitting of the food web data of Fig. 2 alone gives a smaller exponent with a relative error of about 50% and correlation coefficient R² = 0.5.

[26] R V Sol´e, D Alonso, A McKane, Physica A 286, 337 (2000).

[27] K Klemm, V M Eguiluz, Growing scalefree networks with smallworld behavior, Phys. Rev. E 65, 057102 (2002).

[28] E Ravasz, A L Somera, D A Mongru, Z N Oltvai, AL Barab´asi, Hierarchical Organization of Modularity in Metabolic Networks, Science 297, 1551 (2002).

[29] M Brede, S Sinha, Assortative mixing by degree makes a network more unstable, arXiv:condmat/0507710 (2005).

[30] S Sinha, S Sinha, Robust emergent activity in dynamical networks, Phys. Rev. E 74, 066117 (2006).

[31] S Allesina, M Pascual, Network structure, predator prey modules, and stability in large food webs, Theor. Ecol. 1, 55 (2008).

[32] S Sinha, Complexity vs. stability in small world networks, Physica A 346, 147 (2005).

[33] J A Dunne, The Network Structure of Food Webs, In: Ecological networks: Linking structure to dynamics in food webs, Eds. M. Pascual, J A Dunne, pag. 27, Oxford University Press, Oxford (2006).

[34] L A N Amaral, A Scala, M Barth´el´emy, H E Stanley, P. Natl. Acad. Sci. USA 97, 11149 (2000).

[35] T Nishikawa, A Motter, Network synchronization landscape reveals compensatory structures, quantization, and the positive effect of negative interactions, P. Natl. Acad. Sci. USA 107, 10342 (2010).