Unifying Evolutionary and Network Dynamics
APS/123-QED
Unifying Evolutionary and Network Dynamics
Samarth Swarup∗
Department of Computer Science,
University of Illinois at Urbana-Champaign
Les Gasser
Graduate School of Library and Information Science, and
Department of Computer Science,
University of Illinois at Urbana-Champaign.
(Dated: November 17, 2018)
Many important real-world networks manifest “small-world” properties such as scale-free degree
distributions, small diameters, and clustering. The most common model of growth for these networks
is “preferential attachment”, where nodes acquire new links with probability proportional to the
number of links they already have. We show that preferential attachment is a special case of the
process of molecular evolution. We present a new single-parameter model of network growth that
unifies varieties of preferential attachment with the quasispecies equation (which models molecular
evolution), and also with the Erdõs-Rényi random graph model. We suggest some properties of
evolutionary models that might be applied to the study of networks. We also derive the form of the
degree distribution resulting from our algorithm, and we show through simulations that the process
also models aspects of network growth. The unification allows mathematical machinery developed
for evolutionary dynamics to be applied in the study of network dynamics, and vice versa.
PACS numbers: 89.75.Hc, 89.75.Da, 87.23.Kg
Keywords: Evolutionary dynamics, Small-world networks, Scale-free networks, Preferential attachment,
Quasi-species, Urn models.
I. INTRODUCTION
The study of networks has become a very active area
of research since the discovery of “small-world” networks
[1, 2]. Small-world networks are characterized by scale-
free degree distributions, small diameters, and high clus-
tering coefficients. Many real networks, such as neuronal
networks [2], power grids [3], the world wide web [4] and
human language [5], have been shown to be small-world.
Small-worldness has important consequences. For exam-
ple, such networks are found to be resistant to random
attacks, but susceptible to targeted attacks, because of
the power-law nature of the degree distribution.
The process most commonly invoked for the genera-
tion of such networks is called “preferential attachment”
[6, 7]. Briefly, new links attach preferentially to nodes
with more existing links. Simon analyzed this stochas-
tic process, and derived the resulting distribution [8].
This simple process has been shown to generate networks
with many of the characteristics of small-world networks,
and has largely replaced the Erdõs-Rényi random graph
model [9] in modeling and simulation work.
Another major area of research in recent years has been
the consolidation of evolutionary dynamics [10], and its
application to alternate areas of research, such as lan-
guage [11]. This work rests on the foundation of quasi-
species theory [12, 13], which forms the basis of much
subsequent mathematical modeling in theoretical biology.
∗Electronic address: swarup@uiuc.edu
In this paper we bring together network generation
models and evolutionary dynamics models (and partic-
ularly quasi-species theory) by showing that they have
a common underlying probabilistic model. This unified
model relates both processes through a single parameter,
called a transfer matrix. The unification allows mathe-
matical machinery developed for evolutionary dynamics
to be applied in the study of network dynamics, and vice
versa. The rest of this paper is organized as follows: first
we describe the preferential attachment algorithm and
the quasispecies model of evolutionary dynamics. Then
we show that we can describe both of these with a single
probabilistic model. This is followed by a brief analy-
sis, and some simulations, which show that power-law
degree distributions can be generated by the model, and
that the process can also be used to model some aspects
of network growth, such as densification power laws and
shrinking diameters.
II. PREFERENTIAL ATTACHMENT
The Preferential Attachment algorithm specifies a pro-
cess of network growth in which the addition of new (in-
)links to nodes is random, but biased according to the
number of (in-)links the node already has. We identify
each node by a unique type i, and let xi indicate the
proportion of the total number of links in the graph that
is already assigned to node i. Then equation 1 gives the
probablity P (i) of adding a new link to node i [6].
P (i) = αx
i . (1)
http://arxiv.org/abs/0704.1811v1
mailto:swarup@uiuc.edu
where α is a normalizing term, and γ is a constant. As γ
approaches 0 the preference bias disappears; γ > 1 causes
exponentially greater bias from the existing in-degree of
the node.
III. EVOLUTIONARY DYNAMICS AND
QUASISPECIES
Evolutionary dynamics describes a population of types
(species, for example) undergoing change through repli-
cation, mutation, and selection[28]. Suppose there are N
possible types, and let si,t denote the number of individ-
uals of type i in the population at time t. Each type has
a fitness, fi which determines its probability of repro-
duction. At each time step, we select, with probability
proportional to fitness, one individual for reproduction.
Reproduction is noisy, however, and there is a probability
qij that an individual of type j will generate an individ-
ual of type i. The expected value of the change in the
number of individuals of type i at time t is given by,
∆si,t =
j fjsjqij
j fjsj
This is known as the quasispecies equation [13]. The fit-
ness, fi, is a constant for each i. Fitness can also be
frequency-dependent, i.e. it can depend on which other
types are present in the population. In this case the
above equation is known as the replicator-mutator equa-
tion (RME) [10],[14].
IV. A GENERALIZED POLYA’S URN MODEL
THAT DESCRIBES BOTH PROCESSES
Urn models have been used to describe both prefer-
ential attachment [15], and evolutionary processes [16].
Here we describe an urn process derived from the quasis-
pecies equation that also gives a model of network gener-
ation. In addition, this model of network generation will
be seen to unify the Erdõs-Rényi random graph model
[9] with the preferential attachment model.
Our urn process is as follows:
• We have a set of n urns, which are all initially
empty except for one, which has one ball in it.
• We add balls one by one, and a ball goes into urn i
with probability proportional to fimi, where fi is
the “fitness” of urn i, and mi is the number of balls
already in urn i.
• If the ball is put into urn j, then a ball is taken out
of urn j, and moved to urn k with probability qkj .
The matrix Q = [qij ], which we call the transfer matrix,
is the same as the mutation matrix in the quasispecies
equation.
This process describes the preferential attachment
model if we set the fitness, fi, to be proportional tom
where γ is a constant (as in equation 1). Now we get a
network generation algorithm in much the same way as
Chung et al. did [15], where each ball corresponds to a
half-edge, and each urn corresponds to a node. Placing
a ball in an urn corresponds to linking to a node, and
moving a ball from one urn to another corresponds to
rewiring. We call this algorithm Noisy Preferential At-
tachment (NPA). If the transfer matrix is set to be the
identity matrix, Noisy Preferential Attachment reduces
to pure preferential attachment.
In the NPA algorithm, just like in the preferential at-
tachment algorithm, the probability of linking to a node
depends only on the number of in-links to that node. The
“from” node for a new edge is chosen uniformly randomly.
In keeping with standard practice, the graphs in the next
section show only the in-degree distribution. However,
since the “from” nodes are chosen uniformly randomly,
the total degree distribution has the same form. Consider
the case where the transfer matrix is almost diagonal, i.e.
qii is close to 1, and the same ∀i, and all the qij are small
and equal, ∀i 6= j. Let qii = p and
qij =
= q, ∀i 6= j. (3)
Then, the probability of the new ball being placed in bin
P (i) = αm
i p+ (1− αm
i )q, (4)
where α is a normalizing constant. That is, the ball could
be placed in bin i with probability αm
i and then replaced
in bin i with probability p, or it could be placed in some
other bin with probability (1 − αm
i ), and then trans-
ferred to bin i with probability q. Rearranging, we get,
P (i) = αm
i (p− q) + q. (5)
In this case, NPA reduces to preferential attachment with
initial attractiveness [17], where the initial attractiveness
(q, here) is the same for each node. We can get differ-
ent values of initial attractiveness by setting the transfer
matrix to be non-uniform. We can get the Erdõs-Rényi
model by setting the transfer matrix to be entirely uni-
form, i.e. qij = 1/n, ∀i, j. Thus the Erdõs-Rényi model
and the preferential attachment model are seen as two ex-
tremes of the same process, which differ with the transfer
matrix, Q.
This process also obviously describes the evolutionary
process when γ = 1. In this case, we can assume that at
each step we first select a ball from among all the balls in
all the urns with probability proportional to the fitness
of the ball (assuming that the fitness of a ball is the same
as the fitness of the urn in which it is). The probability
that we will choose a ball from urn i is proportional to
fimi. We then replace this ball and add another ball
to the same urn. This is the replication step. This is
followed by a mutation step as before, where we choose
a ball from the urn and either replace it in the urn with
with probability p or move it to any one of the remaining
urns. If we assume that all urns (i.e. all types or species)
have the same intrinsic fitness, then this process reduces
to the preferential attachment process.
Having developed the unified NPA model, we can now
point towards several concepts in quasi-species theory
that are missing from the study of networks, that NPA
makes it possible to investigate:
• Quasi-species theory assumes a genome, a bit string
for example. This allows the use of a distance mea-
sure on the space of types.
• Mutations are often assumed to be point mutations,
i.e. they can flip one bit. This means that a mu-
tation cannot result in just any type being intro-
duced into the population, only a neighbor of the
type that gets mutated.
• This leads to the notion of a quasi-species, which
is a cloud of mutants that are close to the most-fit
type in genome space.
• Quasi-species theory also assumes a fitness land-
scape. This may in fact be flat, leading to neutral
evolution [18]. Another (toy) fitness landscape is
the Sharply Peaked Landscape (SPL), which has
only one peak and therefore does not suffer from
problems of local optima. In general, though, fit-
ness landscapes have many peaks, and the rugged-
ness of the landscape (and how to evaluate it) is
an important concept in evolutionary theory. The
notion of (node) fitness is largely missing from net-
work theory (with a couple of exceptions: [19],
[20]), though the study of networks might benefit
greatly from it.
• The event of a new type entering the population
and “taking over” is known as fixation. This means
that the entire population eventually consists of
this new type. Typically we speak of gene fixa-
tion, i.e. the probability that a single new gene
gets incorporated into all genomes present in the
population. Fixation can occur due to drift (neu-
tral evolution) as well as due to selection.
V. ANALYSIS AND SIMULATIONS
We next derive the degree distribution of the network.
Since there is no “link death” in the NPA algorithm and
the number of nodes is finite, the limiting behavior in our
model is not the same as that of the preferential attach-
ment model (which allows introduction of new nodes).
This means that we cannot re-use Simon’s result [8] di-
rectly to derive the degree distribution of the network
that results from NPA.
A. Derivation of the degree distribution
Suppose there are N urns and n balls at time t. Let
xi,t denote the fraction of urns with i balls at time t.
We choose a ball uniformly at random and “replicate” it,
i.e. we add a new ball (and replace the chosen ball) into
the same urn. Uniformly random choice corresponds to a
model where all the urns have equal intrinsic fitness. We
follow this up by drawing another ball from this urn and
moving it to a uniformly randomly chosen urn (from the
N − 1 other urns) with probability q = (1− p)/(N − 1),
where p is the probability of putting it back in the same
urn. Let P1(i) be the probability that the ball to be
replicated is chosen from an urn with i balls. Let P2(i)
be the probability that the new ball is placed in an urn
with i balls. The net probability that the new ball ends
up in an urn with i balls,
P (i) = P1(i) and P2(i) or P̄1(i) and P2(i). (6)
The probability of selecting a ball from an urn with i
balls,
P1(i) =
Nxi,ti
n0 + t
where n0 is the number of balls in the urns initially. P2(i)
depends on the outcome of the first step.
P2(i) =
p+ (Nxi,t − 1)q when step 1 is “successful”,
Nxi,tq when step 1 is a “failure”.
Putting these together, we get,
P (i) =
Nxi,ti
n0 + t
(p+ (Nxi,t − 1)q) +
Nxi,ti
n0 + t
Nxi,tq
Nxi,ti
n0 + t
(p− q) +Nxi,tq.
Now we calculate the expected value of xi,t+1. xi,t will
increase if the ball goes into an urn with i−1 balls. Sim-
ilarly it will decrease if the ball ends up in an urn with
i balls. Otherwise it will remain unchanged. Remember-
ing that xi,t is the fraction of urns with i balls at time t,
we write,
Nxi,t+1 =
Nxi,t + 1 w. p.
Nxi−1,t(i−1)
(p− q) +Nxi−1,tq,
Nxi,t − 1 w. p.
Nxi,ti
(p− q) +Nxi,tq,
Nxi,t otherwise.
From this, the expected value of xi,t+1 works out to be,
xi,t+1 =
i(p− q)
n0 + t
xi,t+
[(i− 1)(p− q)
n0 + t
xi−1,t.
We can show the approximate solution for xi,t to be,
xi,t =
ri−1Γ(i)
k=1(kr + 1)
(t+ 1)(1− q)t−1, (8)
 0  200  400  600  800  1000
0.5t(0.99)t
t(0.99)t
FIG. 1: Example xi,t curves.
 1  10  100
r=0.33
FIG. 2: The form of the degree distribution.
where r = (p − q)/(1 − q). This approximation is valid
while t << N . See Appendix A for details. For any
particular i, the shape of this curve is given by t(1− q)t.
An example curve is shown in fig 1. This matches our
intuition. Initially, xi,t = 0 for i > 1. As t increases, xi,t
increases through mutations. However, since N is finite
and we keep adding balls, eventually the number of bins
with i balls must go to zero for any particular i. Thus
xi,t must eventually start decreasing, which is what we
see in figure 1. The middle term can be simplified further
k=1(kr + 1)
∏i+1/r
k=1+1/r
∏i+1/r
k=1+1/r
Γ(1/r)
r2Γ(i + 1+ 1/r)
Therefore, in terms of i, equation 8 can be written as
(for fixed t),
xi = C
Γ(i+ 1 + 1
, (9)
where C is a constant. This is the form of the degree
distribution. This is a power law, because as i → ∞,
equation 9 tends to i−(1+1/r) (see discussion of eq. 1.4
in [8, pg 426]). This is also demonstrated in the sample
plots in figure 2.
These results are confirmed through simulation. We
did an experiment where the number of possible nodes
was set to 100000, and 10000 links were added. The
experiment was repeated for values of p ranging from
0.01 to 0.99, in steps of 0.01. Figure 3 shows a plot of
 0  10  20  30  40  50  60  70  80  90  100
p (= 1 - mutation probability)
FIG. 3: N = 100000, number of edges = 10000.
 1000
 10000
 1  10  100  1000  10000
Indegree
FIG. 4: p = 0.8, N = 100000, number of edges = 10000.
coherence, φ, which is defined as,
x2i . (10)
Coherence is a measure of the non-uniformity of the de-
gree distribution. It is 1 when a single node has all the
links. When all nodes have one link each, coherence has
its lowest value, 1/N . We see that as p increases (i.e.
mutation rate decreases), coherence also increases. This
is borne out by the degree distribution plots (figures 4
through 6). The degree distribution is steeper for lower
values of p.
B. Stability
We can rewrite equation 2 as
∆si =
j fjsj
(fisiqii +
j 6=i
fjsjqij) (11)
The first term in the parentheses represents the change in
si due to selection. Some of the copies of type i are lost
due to mutation. The fraction that are retained are given
 1000
 10000
 1  10  100  1000  10000
Indegree
FIG. 5: p = 0.6, N = 100000, number of edges = 10000.
 1000
 10000
 1  10  100  1000  10000
Indegree
FIG. 6: p = 0.4, N = 100000, number of edges = 10000.
by the product fiqii. If this product is greater than 1, the
proportion of type i will increase due to selection, oth-
erwise it will decrease. The second term represents the
contribution to type i due to mutation from all the other
types in the population. Thus, if si decreases towards
zero due to a selective disadvantage, it will be maintained
in the population at “noise” level due to mutations.
This leads to the notion of an error threshold. Sup-
pose that the fitness landscape has only one peak. This
is known as the Sharply Peaked Landscape, or SPL. Sup-
pose further that mutations only alter one position on the
genome at a time. Then it can be shown that if the mu-
tation rate is small enough the population will be closely
clustered about the fittest type. The fittest type keeps
getting regenerated due to selection, and mutations gen-
erate a cloud of individuals with genomes very close to
the genome of the fittest type. This cloud is known as a
quasi-species [21].
If, on the other hand, the mutation rate is above a
certain threshold (essentially 1/fi, where i is the fittest
type) then all types will persist in the population in equal
proportions. This threshold is known as the error thresh-
VI. FITNESS LANDSCAPES AND NEUTRAL
EVOLUTION
We have seen above that noisy preferential attachment
is equivalent to molecular evolution where all intrinsic
fitnesses are equal. If node fitnesses are allowed to be
different, we get standard quasi-species behavior. If the
mutation rate is low enough, the fittest node dominates
the network and acquires nearly all the links. If the mu-
tation rate is high enough to be over the error threshold,
no single node dominates.
Figures 7 and 8 show simulations where nodes are as-
signed intrinsic fitness values uniformly randomly in the
range (0, 1), for different values of p. We see that when p
is high (0.9), i.e. mutation rate is low, the degree distri-
bution stretches out along the bottom, and one or a few
nodes acquire nearly all the links. When p = 0.4, though,
we don’t get this behavior, because the mutation rate is
over the error threshold.
Since we generally don’t see a single node dominating
in real-world networks, we are led to one of two conclu-
sions: either mutation rates in real-world networks are
 1000
 10000
 1  10  100  1000  10000
Indegree
FIG. 7: p = 0.4, N = 100000, number of edges = 10000, node
fitnesses are uniformly randomly distributed between 0 and
 1000
 1  10  100  1000  10000
Indegree
dominant
 node
FIG. 8: p = 0.9, N = 100000, number of edges = 10000, node
fitnesses are uniformly randomly distributed between 0 and
rather high, or the intrinsic fitnesses of the nodes are
all equal. The former seems somewhat untenable. The
latter suggests that most networks undergo neutral evo-
lution [18].
Fitness landscapes can also be dynamic. Golder and
Huberman give examples of short term dynamics in col-
laborative tagging systems (in particular Del.icio.us) [22].
Figures 9 and 10, which are taken from their paper, show
two instances of the rate at which two different web sites
acquired bookmarks. The first one shows a peak right
after it appears, before the rate of bookmarking drops to
a baseline level. The second instance shows a web site
existing for a while before it suddenly shows a peak in
the rate of bookmarking. Both are examples of dynamic,
i.e. changing, fitness. Wilke et al. have shown that in
the case of molecular evolution a rapidly changing fit-
FIG. 9: This is figure 6a from [22]. It shows number of book-
marks received against time (day number). This particular
site acquires a lot of bookmarks almost immediately after it
appears, but thereafter receives few bookmarks.
FIG. 10: This is figure 6b from [22]. It shows number of
bookmarks received against time (day number). This par-
ticular site suddenly acquires a lot of bookmarks in a short
period of time, though it has existed for a long time.
ness landscape is equivalent to the time-averaged fitness
landscape [23]. Thus while short term dynamics show
peaks in link (or bookmark) acquisition, the long-term
dynamics could still be neutral or nearly neutral.
VII. DYNAMICAL PROPERTIES OF
REAL-WORLD NETWORKS
Leskovec et al. point out that though models like pref-
erential attachment are good at generating networks that
match static “snapshots” of real-world networks, they do
not appropriately model how real-world networks change
over time [24]. They point out two main properties which
are observed for several real-world networks over time:
densification power laws, and shrinking diameters. The
term densification power law refers to the fact that the
number of edges grows super-linearly with respect to the
number of nodes in the network. In particular, it grows as
a power law. This means that these networks are getting
more densely connected over time. The second surprising
property of the dynamics of growing real-world networks
is that the diameter (or 90th percentile distance, which
is called the effective diameter) decreases over time. In
most existing models of scale-free network generation, it
has been shown that the diameter increases very slowly
over time [25]. Leskovec et al. stress the importance
of modeling these dynamical aspects of network growth,
and they present an alternate algorithm that displays
both the above properties.
Noisy preferential attachment can also show these
properties if we slowly decrease the mutation rate over
time. Figures 11 and 12 show the effective diameter of the
network and the rate of change of the number of nodes
with respect to the number of edges for a simulation in
which the mutation rate was changed from 0.3 to 0.01
over the course of the simulation run.
 0  20  40  60  80  100  120  140  160  180  200
Number of edges (x100)
FIG. 11: The effective diameter of the network when the mu-
tation rate decreases over time from 0.3 to 0.01. It increases
quickly at first and then decreases slowly over time.
 1000
 10000
 100  1000  10000  100000
Number of edges
FIG. 12: The number of nodes grows as a power law with re-
spect to the number of edges (or time, since one edge is added
at each time step). The slope of the line is approximately 0.86.
VIII. CONCLUSIONS
We have shown that, when modeled appropriately, the
preferential attachment model of network generation can
be seen as a special case of the process of molecular evo-
lution because they share a common underlying proba-
bilistic model. We have presented a new, more general,
model of network generation, based on this underlying
probabilistic model. Further, this new model of network
generation, which we call Noisy Preferential Attachment,
unifies the Erdõs-Rényi random graph model with the
preferential attachment model.
The preferential attachment algorithm assumes that
the fitness of a node depends only on the number of links
it has. This is not true of most real networks. On the
world wide web, for instance, the likelihood of linking to
an existing webpage depends also on the content of that
webpage. Some websites also experience sudden spurts
of popularity, after which they may cease to acquire new
links. Thus the probability of acquiring new links de-
pends on more than the existing degree. This kind of
behavior can be modeled by the Noisy Preferential At-
tachment algorithm by including intrinsic fitness values
for nodes.
The Noisy Preferential Attachment algorithm can also
be used to model some dynamical aspects of network
growth such as densification power laws and shrinking di-
ameters by gradually decreasing mutation rate over time.
If true, this brings up the intriguing question of why mu-
tation rate would decrease over time in real-world net-
works. On the world wide web, for example, this may
have to do with better quality information being avail-
able through the emergence of improved search engines
etc. However, the fact that many different kinds of net-
works exhibit densification and shrinking diameters sug-
gests that there may be some deeper explanation to be
found.
From a design point of view, intentional modulation
of the mutation rate can provide a useful means of trad-
ing off between exploration and exploitation of network
structure. We have been exploring this in the context of
convergence in a population of artificial language learners
[26].
The larger contribution of this work, however, is to
bring together the fields of study of networks and evo-
lutionary dynamics, and we believe that many further
connections can be made.
IX. ACKNOWLEDGEMENTS
We appreciate the helpful comments of Roberto Al-
dunate and Jun Wang. Work supported under NSF
Grant IIS-0340996.
APPENDIX A
Here we solve the difference equation,
xi,t+1 =
i(p− q)
n0 + t
xi,t+
[(i− 1)(p− q)
n0 + t
xi−1,t.
x0,t is a special case.
Nx0,t+1 =
Nx0,t − 1 w. p. Nx0,tq,
Nx0,t otherwise.
Expanding and simplifying as above, we get,
x0,t+1 = (1− q)x0,t.
The solution to this difference equation is simply,
x0,t = (1− q)
tx0,0, (A2)
where x0,0 = (N−1)/N is the initial value of the number
of empty urns. Note that here, and henceforth, we are
assuming that initially all the urns are empty except for
one, which has one ball in it. Therefore x1,0 = 1, and
xi,0 = 0 ∀i > 1. This also means that n0 = 1. These
conditions together specify the entire initial state of the
system.
Equation A1 is difficult to solve directly, so we shall
take the approach of finding the solution to x1,t and x2,t
and then simply guessing the solution to xi,t.
Substituting i = 1 in equation 7 gives us,
x1,t+1 =
(p− q)
n0 + t
x1,t + qx0,t.
Substituting the solution for x0,t from equation A2 gives
x1,t+1 =
(p− q)
n0 + t
x1,t + q(1 − q)
tx0,0. (A3)
The complete solution for x1,t is (see Appendix B),
x1,t = (1 − q)
A(t+ 1) +
, (A4)
where A =
qx0,0
1+p−2q
and B =
2(p−q)
(1+p−2q)NΓ(1−r)
are con-
stants. Let us now use this result to derive the solution
for x2,t. Substituting i = 2 in equation A1, we get,
x2,t+1 =
2(p− q)
n0 + t
x2,t +
[ p− q
n0 + t
x1,t.
Substituting the solution for x1,t from equation A4 and
replacing n0 by 1 for convenience gives us,
x2,t+1 =
2(p− q)
1 + t
x2,t+
(1− q)t
A(t+ 1) +
][p− q
1 + t
. (A5)
The solution to this (after some work) turns out to be
(see Appendix B),
x2,t = (1− q)
A(t+ 1)
2r + 1
q(1− q)t
1 + p− 2q
A(t+ 1)
2rt+ t+ 2r
2(2r + 1)
(t+ 2)
In the above expression, compared to the first term, the
remaining terms are negligible. To see this, consider that
B/tr can be at most B (as r → 0), and at least B/t (as
r → 1). B itself is less than 1/N . Therefore the con-
tribution of the second term is upper-bounded by 1/N .
A similar observation will hold for D/t2r. This is far
less than the contribution due to the first term, since A
(which is also close to 1/N) is multiplied by (t+1). The
remaining terms are approximately of the form t2/N2
(and higher i will contain higher powers). We can ignore
these as long as t << N . Thus, we can write the solution
for x2,t approximately as,
x2,t =
2r + 1
(t+ 1)(1− q)t
1 + p− 2q
N − 1
(t+ 1)(1− q)t
(r + 1)(2r + 1)
(t+ 1)(1− q)t−1.
We can continue on with x3,t:
x3,t+1 =
3(p− q)
1 + t
x3,t +
[2(p− q)
1 + t
x2,t.
If we follow through with this as for x2,t, we will see the
2 from the constant in the second term (
) appear as
a factor in the first term of the solution for x3,t. In the
general expression for the solution, this appears as Γ(i).
Therefore, we can guess the approximate expression for
xi,t to be,
xi,t =
ri−1Γ(i)
k=1(kr + 1)
(t+ 1)(1− q)t−1, (A7)
which is the same as equation 8
APPENDIX B
Equation A3 is,
x1,t+1 =
(p− q)
n0 + t
x1,t + q(1− q)
tx0,0.
This equation is of the form y(t + 1) = p(t)y(t) + r(t).
The general form of the solution is,
y(t) = u(t)
∑ r(t)
Eu(t)
, (B1)
where u(t) is the solution of the homogeneous part of
the above equation, i.e. u(t + 1) = p(t)u(t), and E is
the time-shift operator, i.e. Eu(t) = u(t+ 1). Now, the
homogeneous part of equation A3 is,
u(t+ 1) =
1− q −
n0 + t
( (1− q)t+ (1 − q)n0 − (p− q)
n0 + t
= (1− q)
t+ n0 −
t+ n0
u(t).
The solution to this difference equation is,
u(t) = C(1− q)t
Γ(t+ n0 − r)
Γ(t+ n0)
, (B2)
where r = (p − q)/(1 − q), C is a constant, and Γ(·) is
the gamma-function, which is a “generalization” of the
factorial to the complex plane. It is defined recursively
as Γ(n + 1) = nΓ(n). The derivation of equation B2 is
given in Appendix C. From equations A3, B1, and B2,
we get,
x1,t =
C(1− q)t
Γ(t+ n0 − r)
Γ(t+ n0)
∑ qx0,0(1− q)
tΓ(t+ 1 + n0)
C(1 − q)t+1Γ(t+ 1 + n0 − r)
C(1 − q)t
(t+ n0 − 1)r
[ qx0,0
C(1 − q)
(t+ n0)
r +D1
(tr is read as “t to the r falling”)
q(1− q)t−1x0,0
(t+ n0 − 1)r
(t+ n0)
r + 1
D(1− q)t
(t+ n0 − 1)r
(where D = CD1 is another constant)
q(1− q)tx0,0
1 + p− 2q
Γ(t+ n0 − r)
Γ(t+ n0)
Γ(t+ n0 + 1)
Γ(t+ n0 − r)
D(1− q)t
(t+ n0 − 1)r
q(1− q)tx0,0(t+ n0)
1 + p− 2q
D(1 − q)t
(t+ n0 − 1)r
Let us evaluate the constant by applying the initial con-
ditions t = 0, x0,0 = (N−1)/N , x1,0 = 1/N , and n0 = 1.
We get,
1 + p− 2q
+DΓ(1− r)
q(N − 1)
1 + p− 2q
+NDΓ(1− r).
Therefore, D =
2(p− q)
(1 + p− 2q)NΓ(1− r)
. (B3)
This gives us the complete solution for x1,t as,
x1,t = (1 − q)
A(t+ 1) +
where A =
qx0,0
1+p−2q
and B = D =
2(p−q)
(1+p−2q)NΓ(1−r)
constants. This is the same as equation A4.
1. Solution to equation A5
Equation A5 is,
x2,t+1 =
2(p− q)
1 + t
+(1− q)t
A(t+ 1) +
][p− q
1 + t
Again, this equation is of the form of equation B1. The
solution to the homogeneous part in this case is,
u(t) = C(1− q)t
Γ(t+ 1−
2(p−q)
Γ(t+ 1)
. (B4)
This is found in exactly the same way as equation B2
(see Appendix B). Now, from equations B1, A5, and B4,
we get,
x2,t =
C(1 − q)t
∑ (1− q)t(A(t+ 1) + B
)(p−q
C(1− q)t+1 1
(t+1)2r
C(1 − q)t
C(1− q)
A(p− q)
(t+ 1)2r
(t+ 1)(t+ 1)2r +B(p− q)
∑ (t+ 1)2r
tr(t+ 1)
∑ (t+ 1)2r
Solving the summations (see Appendix C), we get,
x2,t =
C(1 − q)t
C(1− q)
[A(p− q)(t+ 1)2r+1
2r + 1
( t(t+ 1)2r+1
2r + 1
(t+ 1)2r+2
(2r + 1)(2r + 2)
+B(p− q)
(t+ 2)t2r
(1 + r)tr
Simplifying,
x2,t = (1− q)
[Ar(t + 1)
2r + 1
Aq(t+ 1)(2rt+ t+ 2r)
(1 − q)(2r + 1)(2r + 2)
Bq(t+ 2)
(1− q)(1 + r)tr
D(1 − q)t
= (1− q)t
A(t+ 1)
2r + 1
q(1− q)t
1 + p− 2q
A(t+ 1)
2rt+ t+ 2r
2(2r + 1)
(t+ 2)
This is the same as equation B5.
APPENDIX C
1. Derivation of equation B2
Equation B2 is the solution to the following difference
equation:
u(t+ 1) = (1 − q)
t+ n0 −
t+ n0
u(t).
Note that all the factors in this equation are positive.
Taking log, we get,
log u(t+ 1) = log
(1− q)
( t+ n0 − r
t+ n0
+ log u(t),
∆log u(t) = log
(1− q)
( t+ n0 − r
t+ n0
log u(t) =
log(1− q) + log(t+ n0 − r)
−log(t+ n0)
Remembering that
a = ta, and
log(t+a) = logΓ(t+
a), we get,
log u(t) = tlog(1− q) + logΓ(t+ n0 − r)
−logΓ(t+ n0) +D,
Therefore, u(t) = C(1− q)t
Γ(t+ n0 − r)
Γ(t+ n0)
This is the same as equation B2.
2. Derivation of equation B5
Equation B5 is the solution to the following difference
equation:
x2,t =
C(1 − q)t
C(1− q)
A(p− q)
(t+ 1)2r
(t+ 1)(t+ 1)2r +B(p− q)
∑ (t+ 1)2r
tr(t+ 1)
∑ (t+ 1)2r
We shall solve each of the summations individually. At
several points, we will use the summation by parts for-
mula,
Ey(t)∆z(t)
= y(t)z(t)−
z(t)∆y(t)
. (C1)
The first summation term can be obtained directly:
(t+ 1)2r =
(t+ 1)2r+1
2r + 1
+ C1. (C2)
The second summation term can be obtained using the
summation by parts formula. Let Ey(t) = t + 1. Then
y(t) = t, and ∆y(t) = 1. Let ∆z(t) = (t + 1)2r. Then
z(t) =
(t+1)2r+1
. We get,
(t+1)(t+1)2r =
(t+ 1)(t+ 1)2r+1
2r + 1
∑ (t+ 1)2r+1
2r + 1
(t+1)(t+1)2r =
(t+ 1)(t+ 1)2r+1
2r + 1
(t+ 1)2r+2
(2r + 1)(2r + 2)
Before proceeding, we pause to calculate
(1/tr). Note
that,
(t+ 1)r
t+ 1− r
(t+ 1)tr
(t+ 1)tr
Taking summation, we get,
Using the summation by parts formula, we get,
(1− r)tr
We now proceed to the third summation term in the dif-
ference equation for x2,t.
∑ (t+ 1)2r
tr(t+ 1)
∑ t2r−1
We shall again use the summation by parts formula. Let
Ey(t) = t2r−1. Therefore y(t) = (t−1)2r−1, and ∆y(t) =
(2r − 1)(t − 1)2r−2. Let ∆z(t) = 1/tr. Therefore z(t) =
t/(1− r)tr (from equation C4). We get,
∑ t2r−1
t(t− 1)2r−1
(1− r)tr
∑ 2r − 1
t(t− 1)2r−2
t(t− 1)2r−1
(1− r)tr
2r − 1
∑ t2r−1
2r − 1
∑ t2r−1
(t− 1)2r−1
∑ t2r−1
Therefore,
∑ (t+ 1)2r
tr(t+ 1)
The fourth summation term in the difference equation
for x2,t is similar to the third one.
∑ (t+ 1)2r
∑ (t+ 1)2r
tr(t+ 1)
(t+ 1)
Let Ey(t) = (t + 1). Then y(t) = t, and ∆y(t) = 1. Let
∆z(t) =
∑ (t+1)2r
tr(t+1)
. Then z(t) = t
(from equation C5).
Therefore, using the summation by parts rule, we get,
∑ (t+ 1)2r
∑ t2r
∑ t2r
∑ (t+ 1− 2r)t2r−1
(t− 2r)t2r
∑ t2r
t− 2r
1 + r
Substituting back in equation C6, we get,
∑ (t+ 1)2r
( t− 2r
1 + r
t− 2r
1 + r
Therefore, we have,
∑ (t+ 1)2r
(t+ 2)t2r
(1 + r)tr
Combining equations C2, C3, C5, and C7, we get the
solution for x2,t, i.e. equation B5.
[1] S. Milgram, Psychology Today 2, 60 (1967).
[2] D. J. Watts and S. H. Strogatz, Nature 393, 440 (1998).
[3] C. Asavathiratham, S. Roy, B. Lesieutre, and G. Vergh-
ese, IEEE Control Systems (2001).
[4] R. Albert, H. Jeong, and A.-L. Barabási, Nature 401,
130 (1999).
[5] R. Ferrer i Cancho and R. V. Solé, Proceedings of the
Royal Society of London B 268, 2261 (2001).
[6] A.-L. Barabási and R. Albert, Science 286, 509 (1999).
[7] R. Albert and A.-L. Barabási, Physical Review Letters
85, 5234 (2000).
[8] H. A. Simon, Biometrika 42, 425 (1955).
[9] P. Erdõs and A. Rényi, Publicationes Mathematicae De-
brecen 6, 290 (1959).
[10] K. M. Page and M. A. Nowak, Journal of theoretical
biology 219, 93 (2002).
[11] M. A. Nowak, Z. Phys. Chem. 16, 5 (2002).
[12] M. Eigen and P. Schuster, Naturwissenschaften 64, 541
(1977).
[13] M. Eigen, J. McCaskill, and P. Schuster, Journal of Phys-
ical Chemistry 92, 6881 (1988).
[14] N. L. Komarova, Journal of Theoretical Biology 230, 227
(2004).
[15] F. Chung, S. Handjani, and D. Jungreis, Annals of Com-
binatorics 7, 141 (2003).
[16] M. Benäım, S. Schreiber, and P. Tarrès, Annals of Ap-
plied Probability 14, 1455 (2004).
[17] S. N. Dorogovtsev, J. F. F. Mendes, and A. N. Samukhin,
Physical Review Letters 85, 4633 (2000).
[18] M. Kimura, The Neutral Theory of Molecular Evolution
(Cambridge University Press, Cambridge, 1983).
[19] G. Caldarelli, A. Capocci, P. De Los Rios, and M. A.
Muñoz, Physical Review Letters 89 (2002).
[20] G. Bianconi and A.-L. Barabási, Europhysics Letters 54,
436 (2001).
[21] M. Eigen, J. McCaskill, and P. Schuster, Adv. Chem.
Phys. 75, 149 (1989).
[22] S. A. Golder and B. A. Huberman, Journal of Informa-
tion Science 32, 198 (2006).
[23] C. O. Wilke, C. Ronnewinkel, and T. Martinetz, Phys.
Rep. 349, 395 (2001).
[24] J. Leskovec, J. Kleinberg, and C. Faloutsos, in Proceed-
ings of KDD’05 (Chicago, Illinois, USA, 2005).
[25] B. Bollobás and O. Riordan, Combinatorica 24, 5 (2004).
[26] S. Swarup and L. Gasser, in From Animals to Animats 9:
Proceedings of the Ninth International Conference on the
Simulation of Adaptive Behavior (Rome, Italy, 2006).
[27] G. U. Yule, Philosophical Transactions of the Royal Soci-
ety of London. Series B, Containing Papers of a Biological
Character 213, 21 (1925).
[28] Simon (and Yule [27] before him) applied their stochas-
tic model to the estimation of numbers of species within
genera, but the notion of quasi-species was unknown at
the time, and it addresses a much wider range of issues
than species frequency.
