Papers in Physics, vol. 2, art. 020005 (2010)

Received: 17 July 2010, Accepted: 27 September 2010
Edited by: D. H. Zanette
Reviewed by: V. M. Eguiluz, Inst. F́ısica Interdisciplinar y Sist. Complejos,

Palma de Mallorca, Spain
Licence: Creative Commons Attribution 3.0
DOI: 10.4279/PIP.020005

www.papersinphysics.org

ISSN 1852-4249

Stability as a natural selection mechanism on interacting networks

Juan I. Perotti,1, 2∗ Orlando V. Billoni,1, 2† Francisco A. Tamarit,1, 2‡

Sergio A. Cannas1, 2§

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, et al., Phys. Rev. Lett. 103, 108701 (2009)]. This mech-
anism 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 ob-
served 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.

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-4]). Although not new, with the ad-
vent of high performance computing, this theoreti-
cal construction opened a new door for the statisti-
cal 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 in-
sights about the dynamical behavior of systems as

∗E-mail: juanpool@gmail.com
†E-mail: billoni@famaf.unc.edu.ar
‡E-mail: tamarit@famaf.unc.edu.ar
§E-mail: cannas@famaf.unc.edu.ar

1 Facultad de Matemática, Astronomı́a y F́ısica, Universi-
dad Nacional de Córdoba, Argentina.

2 Instituto de F́ısica Enrique Gaviola (IFEG-CONICET),
Ciudad Universitaria, 5000 Córdoba, Argentina.

complex as biological and social systems. In ad-
dition, 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 interac-
tion network for a given system is frequently not
unique (see for example the case of protein-protein
interaction networks [5-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 consid-
ered 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

020005-1



Papers in Physics, vol. 2, art. 020005 (2010) / J. I. Perotti et al.

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 kmax; in
the latter, the degree distribution is generally well
described by P (k) ∼ k−γ e−k/kmax [1, 7-12]. In
any case, the networks present a nonhomogeneous
structure, very different from that expected in a
random network. (c) Scaling of the clustering co-
efficient: 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 β tak-
ing 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 con-
nected 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 macro-
scopic level of communities of living beings (eco-
logical networks). Such ubiquity suggests the ex-
istence of some natural selection process that pro-
motes the development of those particular struc-
tures [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 dur-
ing a slow evolutionary (i.e., species lifetime) pro-
cess. A second one is the node dynamics which
affects and in turn is affected by the growing pro-
cesses. It is reasonable to expect that the network
topologies we finally witness could have emerged
out of these coupled processes. Consider, for ex-
ample, the case of an ecological network like a food
web, where nodes are species within an ecosystem
and edges are consumer-resource relationships be-
tween them. New nodes are added during evo-
lutionary time scales, through speciation or mi-
gration of new species. Then, the network grows
through community assembly rules, strongly influ-
enced 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 repre-
sented 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 sta-
bility 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 fea-
tures that emerge in growing networks under sta-
bility constraint. In section IV., we show that this
constraint not only induces topological features of
the resulting networks but also influences the as-
sociated 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 differen-
tial equations d~x/dt = ~F (~x), where ~x is an n-
component vector describing the relevant state
variables of each agent and ~F is an arbitrary non-
linear function. One could imagine that ~x in differ-
ent systems may represent concentrations of some
hormones, the average density populations in a
food web, the concentration of chemicals in a bio-
chemical 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 be-
longing to that set. This defines the interaction
network.

We assume that there are two time scales in the
dynamics. Let fm be the average frequency of the
incoming flux of new agents (migration, mutation,
etc.). This defines a characteristic time τm = f−1m .
On the long time scale t � τm (much larger than
the observation time) new agents arrive to the sys-
tem and start to interact with some of the previous
ones. Some of them can be incorporated into the

020005-2



Papers in Physics, vol. 2, art. 020005 (2010) / J. I. Perotti et al.

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 sta-
tionary 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 Ja-
cobian matrix

ai,j ≡
(

∂Fi
∂xj

)
~x∗

. (1)

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 ai,j
are such that the eigenvalue with maximum real
part λm of the enlarged Jacobian matrix is nega-
tive (λm < 0). Assuming that isolated agents will
reach stable states by themselves after certain char-
acteristic relaxation time, the diagonal elements of
the matrix ai,i are negative and given unity value
to further simplify the treatment [18]. The inter-
action values, (i.e., the non-diagonal matrix ele-
ments ai,j ) will take random values (both positive
and negative) taken from some statistical distri-
bution. In this way, we have an unbounded en-
semble 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, af-
ter 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 algo-
rithm [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 suc-
cessful 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 al-
ready created network with n nodes, such that the
n × n associated interaction matrix ai,j is stable.
Then, for the attachment of the (n + 1)th node,
we first choose its degree kn+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
(ai,n+1, an+1,i) (i = 1, . . . , n) are zero with proba-
bility 1−kn+1/n and different from zero with prob-
ability kn+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 result-
ing (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 en-
tire 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 ≤ kn+1 extinctions,
taken from the set of kn+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 kn+1
and remove it. If the resulting n × n matrix is sta-
ble, we start a new trial; otherwise, another node
among the remaining kn+1 − 1 is chosen and re-
moved, 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 re-
peated until the network reaches a maximum size
n = nmax (typically nmax = 200) and restarted M
times from n = 1 to obtain statistics of the net-
works (typically M = 105).

III. Topological properties

i. Connectivity

First, we analyzed the average connectivity C(n),
defined as the fraction of non-diagonal matrix el-
ements different from zero, averaged over different
runs. In Fig. 1, we show the typical behavior of

020005-3



Papers in Physics, vol. 2, art. 020005 (2010) / J. I. Perotti et al.

C(n) for different values of b (we found that C(n)
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

C(n) ∼ α−ω n−(1+�) (2)

for large values of n, where α is the variance of the
non-diagonal elements of the stability matrix (α =
b2/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 net-
works [18] C(n) = (αn)−1. It is easy to see that
Eq. (2) lies above May’s stability line for network
sizes up to ∼ 106 [21]. This shows that networks
growing under stability constraint develop partic-
ular 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 de-
velopment process. In the next sections, we explore
the characteristics of those networks.

In Fig. 2, we plotted the connectivity for differ-
ent 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 men-
tioning 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 summa-
rize 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 = nmax. The degree exponent
γ takes values between 2 and 3 for values of b in
the interval b ∈ (1.5, 3.5), which become almost in-
dependent of q as it increases. The exponent γ can
also fall below 2 when the global stability constraint

Figure 1: Connectivity as a function of the network
size for q = 3, nmax = 200 and different values of b.
The symbols correspond to numerical simulations
and the dashed lines to power law fittings of the
tails C(n) = B n−(1+�). The insets show the fitting
values B and � as a function of b

Figure 2: Connectivity as a function of the network
size for different biological networks. The straight
line is a power law fitting C(n) = A n−(1+�), giv-
ing an exponent � = 0.2 ± 0.1 (R2 = 0.92). Data
extracted from: [6, 7] (protein-protein interaction
networks); [8] (metabolic networks); [22, 23] (food
webs).

020005-4



Papers in Physics, vol. 2, art. 020005 (2010) / J. I. Perotti et al.

Figure 3: Degree distribution P (k) for q = 3,
nmax = 200 and different values of b; the dashed
lines correspond to power law fittings of the tail
P (k) ∼ k−γ . Logarithmic binning has been used
to smooth the curves.

is replaced by a local one. The qualitative struc-
ture 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
that the characteristic time to leave the fixed point
τ = λ−1m is large enough to become comparable to
the migration time scale τm [19].

iii. Network growth and clustering proper-
ties

Networks grown under stability constraint also dis-
play 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 be-
tween 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 ap-
proximated 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 mecha-
nism, a detailed analysis has shown that this is not

Figure 4: Average network size as a function of the
time measured in number of trials for b = 2 and
q = 3. The continuous red line corresponds to the
power law t1/2,

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 ai,j stability much less than
new nodes with many links (of course this is re-
flected 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 con-
nected node is accepted, decreasing L and increas-
ing Cc(n) [19]. Those fluctuations lead to a slow
diffusive-like growth of the network size n(t) ∼ t1/2
(See Fig. 4).

Another quantity of interest is the average clus-
tering Cc(k) as a function of the degree k. A typical
example is shown in Fig. 5. We see that Cc(k) de-
creases 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 com-
pletely independent of b and q. This behavior is
indicative of a modular structure with hierarchical
organization [13]. Notice that this power law de-
cay 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.).

020005-5



Papers in Physics, vol. 2, art. 020005 (2010) / J. I. Perotti et al.

Figure 5: Average clustering coefficient Cc(k) as a
function of the degree k for b = 2, q = 3 and dif-
ferent values of nmax. The dashed black line corre-
sponds to the power law k−0.9.

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 knn among the near-
est neighbors of a node with degree k. In Fig. 6, we
see that knn decays with a power law knn ∼ 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 topo-
logical properties that are selected by the stabil-
ity constraint, i.e., properties associated to the un-
derlying 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 aij 6= 0.

First of all, we calculate the probability distribu-

Figure 6: Average nearest neighbors degree knn(k)
of a node with degree k for q = 3 and different
values of b and nmax. The dashed line corresponds
to the power law k−0.25.

tion of values for a single non null matrix element
aij of the final network with size n = nmax. The
typical behavior is shown in Fig. 7. We see that
P (aij ) is an even function, almost uniform in the
interval [−b, b], with a small cusp around aij = 0.
This shows that stability is not enhanced by a par-
ticular sign or absolute value of the individual in-
teraction 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(aij ) = −sign(aji)) significa-
tively enhances the stability of random matrices
[31]. In an ecological network, this typically cor-
responds to a predator-prey or parasite-host inter-
action. To check for the presence of such type of
interactions, we calculated the correlation 〈aij aji〉,
where the average is taken over pairs of nodes with
a double link (aij 6= 0 and aji 6= 0).

In Fig. 8, we show 〈aij aji〉 as a function of the
network size n. We see that this correlation is neg-
ative for any value of n and saturates into a value
〈aij aji〉 ≈ −0.65 for large values of n. In the in-
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 probabil-

020005-6



Papers in Physics, vol. 2, art. 020005 (2010) / J. I. Perotti et al.

Figure 7: Probability density of the matrix ele-
ments aij for b = 2, q = 3 and nmax = 100.

ity P (aij 6= 0) = C(n). Then, the probability of
having a link between an arbitrary pair of sites is
pd = C(n)(2 − C(n)) and the average degree per
node 〈k〉 = pdn. Hence

〈η〉ran =
C2(n) n
〈k〉

=
C(n)

2 − C(n)
(3)

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 to-
gether show that the present networks have indeed
a significantly large number of anticorrelated pair
interactions.

Next, we calculated the correlation
〈aij aji〉/〈|aij|〉2 between the matrix elements,
linking a node i and its neighbors j, as a function
of its degree ki, 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 ki = 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

Figure 8: Correlation function 〈aij aji〉 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 aver-
age was calculated over all pair of sites with dou-
ble 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). Yel-
low triangles correspond to a numerical calculation
of 〈η〉ran. The dashed line corresponds to a power
law n−0.68.

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 net-
work 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 consid-
ered 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 ad-

020005-7



Papers in Physics, vol. 2, art. 020005 (2010) / J. I. Perotti et al.

Figure 9: Correlation function 〈aij aji〉/〈|aij|〉2 be-
tween the matrix elements linking a node i and its
neighbors j as a function of its degree ki, for b = 2,
q = 3 and different values of nmax. The inset shows
the average fraction of anticorrelated links 〈κ〉 as a
function of ki.

dition, 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 (disassortativ-
ity, low clustering, uncorrelated links), while low
connected nodes inside modules enhance stability
through the presence of many anticorrelated inter-
actions.

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 se-
lection process.

In particular, the above scenario agrees very well
with the observed structures in cellular networks.
For instance, the scaling behavior of Cc(k), dis-
played 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 cer-
tain protein-protein interaction networks [7]. Also,
the available data for the degree distribution in all
those cases are consistent with a power law behav-
ior with an exponent γ between 2 and 2.5 [1]. The
agreement with the whole set of properties pre-
dicted 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 eco-
logical networks, where the predictions of the model
do not completely agree with the observations, spe-
cially those related to food webs. On the one hand,
food webs usually display also disassortative mix-
ing 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 scal-
ing 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 resolu-
tion of the data [25]. The consistency of the scal-
ing shown in Fig. 2 for a broad range of size scales
suggests that the ecological debate should be recon-
sidered 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 fre-
quently exhibits an exponential cutoff at some max-
imum characteristic degree kmax [9, 10]. Such vari-
ance between food webs and other biological net-
works is probably related to the way ecosystems
assemble and evolve compared with other systems.
While the hypothesis of the present model are gen-
eral 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 dis-
tribution: aging and limited capacity of the nodes
[34]. In the former, nodes can become inactive with

020005-8



Papers in Physics, vol. 2, art. 020005 (2010) / J. I. Perotti et al.

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 situ-
ations in which mechanisms of that type become
important in the evolution of ecological webs, ei-
ther by limitations in the available resources or by
dynamical changes in the diet of species due to ex-
ternal 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 incorpo-
rated into the model, serving as a base for the de-
scription of more complex behaviors in particular
systems like food webs.

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

Acknowledgements - This work was supported
by CONICET, Universidad Nacional de Córdoba,
and FONCyT grant PICT-2005 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, A-L Barabási, Statistical mechan-
ics 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 topo-
logical features of proteins from amino acid

residue networks, Physica A 375, 336
(2007).

[6] S H Yook, Z N Oltvai, A-L Barabási, Func-
tional and topological characterization of
protein interaction networks, Proteomics 4,
928 (2004).

[7] V Colizza, A Flammini, A Maritan, A
Vespignani, Characterization and model-
ing of protein–protein interaction networks,
Physica A 352, 1 (2005).

[8] H Jeong, B Tombor, R Albert, Z N Olt-
vai, A-L Barabási, The large-scale organiza-
tion of metabolic networks, Nature 407, 651
(2000).

[9] J M Montoya, S L Pimm, R V Solé, Ecologi-
cal 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, Scale-Free Brain
Functional Networks, Phys. Rev. Lett. 94,
018102 (2005).

[13] E Ravasz, A-L Barabási, Hierarchical orga-
nization in complex networks, Phys. Rev. E
67, 026112 (2003).

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

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

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

[17] In particular cases, like food-webs, there
could be other coupled process almost prob-
ably as important as the two considered

020005-9



Papers in Physics, vol. 2, art. 020005 (2010) / J. I. Perotti et al.

here, such as the dynamical change in previ-
ously 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 incom-
ing species. The present scheme should be
considered as a minimal model about sta-
bility 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 Self-
Organized Complex Network Topology out of
Stability Constraints, Phys. Rev. Lett. 103,
108701 (2009).

[20] Some of the numerical calculations per-
formed in section III. were repeated replac-
ing the uniform distribution for the non di-
agonal elements by a Gaussian one with the
same variance. The results were indistin-
guishable.

[21] An extrapolation of the present results sug-
gest that the probability of growing goes
to zero above some maximum network size.
Preliminary estimations give a value ≈ 105
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 correspond-
ing results will be published in the near fu-
ture.

[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 n2. It is pro-
portional 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 ex-
ponent with a relative error of about 50%
and correlation coefficient R2 = 0.5.

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

[27] K Klemm, V M Eguiluz, Growing scale-free
networks with small-world behavior, Phys.
Rev. E 65, 057102 (2002).

[28] E Ravasz, A L Somera, D A Mongru, Z N
Oltvai, A-L Barabási, Hierarchical Organi-
zation of Modularity in Metabolic Networks,
Science 297, 1551 (2002).

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

[30] S Sinha, S Sinha, Robust emergent activ-
ity 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 Uni-
versity Press, Oxford (2006).

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

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

020005-10