Stability as a natural selection mechanism on interacting networks

Biological networks of interacting agents exhibit similar topological properties for a wide range of scales, from cellular to ecological levels, suggesting the existence of a common evolutionary origin. A general evolutionary mechanism based on global stability has been proposed recently [J I Perotti, O V Billoni, F A Tamarit, D R Chialvo, S A Cannas, Phys. Rev. Lett. 103, 108701 (2009)]. This mechanism is 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. We show that, out of this stability constraint, several topological properties observed in biological networks emerge in a self organized manner. The influence of the stability selection mechanism on the dynamics associated to the resulting network is analyzed as well.


I. Introduction
The concept of networks of interacting agents has proven, in the last decade, to be a powerful tool in the analysis of complex systems (for reviews, see Refs. [1][2][3][4]). Although not new, with the advent of high performance computing, this theoretical construction opened a new door for the statistical physics methodology in the analysis of systems composed by a large number of units that interact in a complicated way. This allowed to get new insights about the dynamical behavior of systems as complex as biological and social systems. In addition, it constitutes a basic backbone upon which relatively simple models can be constructed in a bottom-up strategy.
As a modeling tool, the definition of an interaction network for a given system is frequently not unique (see for example the case of protein-protein interaction networks [5][6][7]), depending on the coarse grain level of the approach. Nevertheless, many topological properties appear to be independent of the definition of the network. Moreover, some of those properties have emerged in the last years as universal features among systems otherwise considered very different from each other. In particular, the following properties are characteristic of most biological networks. (a) Small worldness: all of them exhibit high clustering Cc and relatively short path length L, compared with random networks. L is defined as the minimum number of links needed to connect any pair of nodes in the network and Cc is defined as the fraction of connections between topological neighbors of any site [1]. (b) Scale free degree distribution: the degree distribution P (k) (the probability of a node to be connected to k other ones) presents a broad tail for large values of k. In some cases, the tail can be approached by a power law P (k) ∼ k −γ with degree exponents γ < 3 for a wide range of scales, while in others, a cutoff appears for some maximum degree k max ; in the latter, the degree distribution is generally well described by P (k) ∼ k −γ e −k/kmax [1,[7][8][9][10][11][12]. In any case, the networks present a nonhomogeneous structure, very different from that expected in a random network. (c) Scaling of the clustering coefficient : in many natural networks, it is observed that the clustering coefficient of a node with degree k follows the scaling law Cc(k) ∼ k −β , with β taking values close to one. This has been interpreted as an evidence for a modular structure organized in a hierarchical way [13]. (d) Disassortative mixing by degree: in most biological networks, highly connected nodes tend to be preferentially connected to nodes with low degree and vice versa [14].
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 first 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 consumer-resource 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 specific interactions among them [15,16]. The con-sequence 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 definition 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 d x/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 k i < n other agents; thus, F i 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 τ m = f −1 m . On the long time scale t ≫ τ 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 τ rel ≪ τ m . Then, in the short time scale τ rel ≪ t ≪ τ m we can assume that n is constant and the dynamics already led the system to a particular stable stationary state x * defined by F ( x * ) = 0. Following May's ideas [18], we assume that the only attractors of the dynamics are fixed 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 maximum real part of the Jacobian matrix A new agent will be incorporated to the network if its inclusion results in a new stable fixed point, that is, if the values of the interaction matrix a i,j are such that the eigenvalue with maximum real part λ m of the enlarged Jacobian matrix is negative (λ m < 0). Assuming that isolated agents will reach stable states by themselves after certain characteristic relaxation time, the diagonal elements of the matrix a i,i are negative and given unity value to further simplify the treatment [18]. The interaction values, (i.e., the non-diagonal matrix elements a i,j ) will take random values (both positive and negative) taken from some statistical distribution. In this way, we have an unbounded ensemble of systems [18] characterized by a "growing through stability" history. Randomness would be self-generated through the addition of new agents processes. Each specific 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 algorithm [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 × n matrix grows its size by one column and one row. Otherwise, the novate 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 × n associated interaction matrix a i,j is stable. Then, for the attachment of the (n + 1) th node, we first choose its degree k n+1 randomly between 1 and n with equal probability. Then, the new agent interaction with the existing network member i is chosen, such that non-diagonal matrix elements (a i,n+1 , a n+1,i ) (i = 1, . . . , n) are zero with probability 1−k n+1 /n and different from zero with probability k n+1 /n; to each non-zero matrix element we assign a different real random value uniformly distributed in [−b, b]. b determines the interaction range variability and it is one of the two parameters of the model [20].
Then, we calculate numerically λ m for the resulting (n + 1) × (n + 1) matrix. If λ m < 0, the new node is accepted. If λ 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+1 extinctions, taken from the set of k n+1 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 × n matrix is stable, we start a new trial; otherwise, another node among the remaining k n+1 − 1 is chosen and removed, repeating the previous procedure. If after q removals the matrix remains unstable, the new node is removed (we return to the original n × n matrix and start a new trial). The process is repeated until the network reaches a maximum size n = n max (typically n max = 200) and restarted M times from n = 1 to obtain statistics of the networks (typically M = 10 5 ).

III. Topological properties i. Connectivity
First, we analyzed the average connectivity C(n), defined as the fraction of non-diagonal matrix elements different from zero, averaged over different runs. In Fig. 1, we show the typical behavior of is completely independent of q). The connectivity presents a power law tail for large values of n. From a fitting of the tail with a power law (see insets in Fig. 1) we obtain the scaling behavior for large values of n, where α is the variance of the non-diagonal elements of the stability matrix (α = b 2 /3 for the uniform distribution) and ω = 0.7±0.1. From the inset of Fig. 1, we see that the exponent ǫ shows a weak dependency on b, taking values 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) = (αn) −1 . It is easy to see that Eq. (2) lies above May's stability line for network sizes up to ∼ 10 6 [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 mag- nitude of network size scales, using data collected from the literature. We see that the data are very well fitted by a single power law C(n) ∼ n −1.2 , in a nice agreement with the average value ǫ = 0.2 predicted by the present model. It is worth mentioning that the behavior C(n) ∼ n −(1+ǫ) 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 of P (k). It presents a power law tail P (k) ∼ k −γ for values of k > 20, with a finite size drop at k = n max . The degree exponent γ takes values between 2 and 3 for values of b in the interval b ∈ (1.5, 3.5), which become almost independent of q as it increases. The exponent γ can also fall below 2 when the global stability constraint is replaced by a local one. The qualitative structure of P (k) remains when the stability criterium λ m < 0 is relaxed by the condition λ m < ∆, with ∆ some small positive number. In other words, the power law tail emerges also when the addition of new nodes destabilizes the dynamics, provided

iii. Network growth and clustering properties
Networks grown under stability constraint also display small world properties. The average clustering coefficient decays with the network size as Cc(n) ∼ n −0.75 (which is slower than the 1/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 Barabasi-Albert model [1], where the clustering can be approximated by a power law with the same exponent, although the exact scaling is [27] Cc(n) ∼ (ln n) 2 /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 Cc 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 a i,j 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 Cc, but sporadically, a highly connected node is accepted, decreasing L and increasing Cc(n) [19]. Those fluctuations lead to a slow diffusive-like growth of the network size n(t) ∼ t 1/2 (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 β ≈ 0.9, close 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 values 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 ∼ k −δ for k > 20, with an exponent δ close to −0.25, in a clear disassortative behavior. This result is also consistent with previous works showing that assor- tative mixing by degree decreases the stability of a network, i.e., the maximum real part λ 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, i.e., properties associated to the underlying adjacency matrix, regardless of the values of the interaction strengths. We now analyze the characteristics of the dynamics associated to the networks emerging from such constraint. In other words, we investigate the statistics of values of the non null elements a ij = 0.
First of all, we calculate the probability distribution of values for a single non null matrix element a ij of the final network with size n = n max . The typical behavior is shown in Fig. 7. We see that P (a ij ) is an even function, almost uniform in the interval [−b, b], with a small cusp around a ij = 0. This shows that stability is not enhanced by a particular sign or absolute value of the individual interaction coefficients. It has been shown recently that the presence of anticorrelated links between pairs of nodes (i.e., links between pairs of nodes (i, j) such that sign(a ij ) = −sign(a ji )) significa-  tively enhances the stability of random matrices [31]. In an ecological network, this typically corresponds to a predator-prey or parasite-host interaction. To check for the presence of such type of interactions, we calculated the correlation a ij a ji , where the average is taken over pairs of nodes with a double link (a ij = 0 and a ji = 0).
In Fig. 8, we show a ij a ji as a function of the network size n. We see that this correlation is negative for any value of n and saturates into a value a ij a ji ≈ −0.65 for large values of n. In the in- Figure 8: Correlation function a ij a ji between a pair of nodes (i, j) with a double link as a function of the network size, for b = 2 and q = 3. The average was calculated over all pair of sites with double link in networks with the same size. The inset shows a comparison between the fraction of double links in the present network η (green circles) and a random network of the same size and connectivity C(n): η ran = C/(2 − C) (continuous line). Yellow triangles correspond to a numerical calculation of η ran . The dashed line corresponds to a power law n −0.68 . set of Fig. 8, we compare the average fraction of double links η 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 probability P (a ij = 0) = C(n). Then, the probability of having a link between an arbitrary pair of sites is p d = C(n)(2 − C(n)) and the average degree per node k = p d n. Hence Then for large values of n, we have η ran ∼ C(n) ∼ n −(1+ǫ) . From the inset of Fig. 8, we see that η ∼ n −0.68 when n ≫ 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 significantly large number of anticorrelated pair interactions. Next, we calculated the correlation Figure 9: Correlation function a ij a ji / |a ij | 2 between the matrix elements linking a node i and its neighbors j as a function of its degree k i , for b = 2, q = 3 and different values of n max . The inset shows the average fraction of anticorrelated links κ as a function of k i .
a ij a ji / |a ij | 2 between the matrix elements, linking a node i and its neighbors j, as a function of its degree k i , where the average is taken only on the double links. From Fig. 9, we see that the absolute value of the correlation presents a maximum around k i = 25 and tends to zero as the degree increases. The inset of Fig. 9 shows that the average fraction of anticorrelated links κ (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 emergency 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 fitness 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 conclusion. 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 specific 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 protein-protein 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 γ 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 values of clustering, com-pared with other biological networks), in agreement with the present predictions. Regarding the scaling behavior of C(n) [24], this is the topic of an old 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 value of the exponent (and the associated interpretations) is controversial, due to the large dispersion 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 maximum characteristic degree k max [9,10]. Such variance 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 principle 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 generate 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, either 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].