Evolutionary games on minimally structured populations
Evolutionary games on minimally
structured populations
Gergely J. Szöllősi∗ and Imre Derényi†
Biological Physics Department Eötvös University, Budapest
Abstract
Population structure induced by both spatial embedding and more general networks of interaction, such
as model social networks, have been shown to have a fundamental effect on the dynamics and outcome of
evolutionary games. These effects have, however, proved to be sensitive to the details of the underlying
topology and dynamics. Here we introduce a minimal population structure that is described by two distinct
hierarchical levels of interaction, similar to the structured metapopulation concept of ecology and island
models in population genetics. We believe this model is able to identify effects of spatial structure that
do not depend on the details of the topology. While effects depending on such details clearly lie outside
the scope of our approach, we expect that those we are able to reproduce should be generally applicable
to a wide range of models. We derive the dynamics governing the evolution of a system starting from
fundamental individual level stochastic processes through two successive meanfield approximations. In
our model of population structure the topology of interactions is described by only two parameters: the
effective population size at the local scale and the relative strength of local dynamics to global mixing. We
demonstrate, for example, the existence of a continuous transition leading to the dominance of cooperation
in populations with hierarchical levels of unstructured mixing as the benefit to cost ratio becomes smaller
then the local population size. Applying our model of spatial structure to the repeated prisoner’s dilemma
we uncover a novel and counterintuitive mechanism by which the constant influx of defectors sustains
cooperation. Further exploring the phase space of the repeated prisoner’s dilemma and also of the “rock-
paper-scissor” game we find indications of rich structure and are able to reproduce several effects observed
in other models with explicit spatial embedding, such as the maintenance of biodiversity and the emergence
of global oscillations.
PACS numbers: 87.10.+e 87.23.-n
∗ssolo@angel.elte.hu; angel.elte.hu/˜ssolo
†derenyi@angel.elte.hu; angel.elte.hu/˜derenyi
http://arxiv.org/abs/0704.0357v3
mailto:ssolo@angel.elte.hu
angel.elte.hu/~ssolo
mailto:derenyi@angel.elte.hu
angel.elte.hu/~derenyi
I. INTRODUCTION
The dynamics of Darwinian evolution is intrinsically frequency dependent, the fitness of in-
dividuals is tightly coupled to the type and number of competitors. Evolutionary dynamics acts,
however, on populations, not individuals and as a consequence depends on not only population
composition, but also population size and structure. Evolutionary game theory came about as the
result of the realization that frequency dependent fitness introduces strategic aspects to evolution
[1, 2, 3]. More recently the investigation of the evolutionary dynamics of structured populations,
where individuals only compete with some subset of the population, e.g. their neighbors in space
or more generally in some graph [4, 5], has lead to the recognition that the success of different
strategies can be greatly influenced by the topology of interactions within the population. Funda-
mental differences were found – compared to well-mixed populations, where individuals interact
with randomly chosen partners – in models that describe the evolution of cooperation (variants of
the prisoner’s dilemma game [4, 6, 7, 8, 9]) or deal with the maintenance of biodiversity in the
context of competitive cycles (variants of the rock-paper-scissors game [3, 10, 11, 12, 13, 14]).
In order to investigate the coevolutionary dynamics of games on structured populations the
full set of connections between a potentially very large number of individuals must be specified.
This is only possible by reducing the number of degrees of freedom considered, either through
postulating a highly symmetrical (such as lattices [4, 8, 16, 17, 18, 19, 20]) or fundamentally
random connection structure (such as some random graph ensemble [21, 22]). The question of how
one goes about the task of reducing the number of degrees of freedom – of choosing the relevant
parameters to describe the population structure constrained to which individuals undergo evolution
– is not trivial. Both the explicit spatial as well as the random graph ensemble approach have clear
precedents in condensed matter physics and network theory, respectively. It is not, however, clear
which – if either – approach best describes natural populations of cyclically competing species or
societies composed of individuals playing the prisoner’s dilemma game.
As an example let us consider colicin producing bacteria, that play the so called ”rock-paper-
scissors” (RPS) game (for details see below). This system has recently been the subject of two
experimental studies aimed at demonstrating the role of structured populations in the maintenance
of diversity. In the first study [10, 11] bacteria were cultured in vitro in Petri dishes, effectively re-
stricting competition between bacteria to neighbors on the (2D) Petri dish surface (Fig.1 top left),
while in the second experiment [12] in vivo bacterial colonies were established in co-caged mice
and their development was subsequently followed. In the case of the first experiment the analogy
with explicit 2D spatial embedding (present by construction) is clear (Fig.1 bottom left). The pop-
ulation structure of the second experiment is, however, clearly different. The bacteria in individual
mice can be readily considered as locally well-mixed populations, the coevolutionary dynamics
of which reduces in the standard meanfield limit to a system of non-linear differential equations
(the adjusted replicator equations [24]). As the experiments show, however, migration of bacteria
between mice may also occur – resulting in the observed cyclic presence of the three strains in
individuals. There are two distinct scales of mixing present in the system. Bacteria within each
mice compete with each other forming local populations – an unstructured neighborhood com-
posed of individual bacteria, while also being exposed to migrants from mice with whom they
share the cage, together forming a global population – an unstructured neighborhood composed of
individual local populations (Fig.1 top and bottom right). This setup is referred to in the ecology
literature – albeit in significantly different contexts – as a ”structured metapopulation” [26, 27]
where structured here refers to the detailed consideration of the population dynamics of the indi-
vidual populations (often called ”patches”) comprising the metapopulation and is also related to
the finite island models of population genetics [28].
The above example of co-caged mice is not unique, we may readily think of other ecological or
sociological examples where an approximation with hierarchical scales of mixing with no internal
structure can be relevant (such as human societies with two distinct scales of mixing present, the
first within individual nations the between them at an international level). We have, also, recently
used a similar approach to construct a model of genetic exchange among bacteria of the same
species (the bacterial equivalent of sex) with which we were able to take into account the effects of
spatial and temporal fluctuations in a manner that can explain the benefit of such genetic exchange
at the level of the individual [31].
In this paper we construct a hierarchical meanfield theory where the two distinct (i.e. local and
global) scales of mixing are each taken into account in terms of two separate meanfield approx-
imations and fluctuations resulting from finite population size on the local scale of mixing are
also considered. We subsequently explore the similarities and differences between this and other
models of structured populations in the case of the ”rock-paper-scissors” and prisoner’s dilemma
games. Through these examples we suggest that our approach allows the separation of the effects
of structured populations on coevolutionary dynamics into effects which are highly sensitive to
and dependent on the details of the topology and those which only require the minimal structure
Explicit spatial embeding Two distinct scales of mixing
Petri dish experiment Co-caged mice experiment
FIG. 1: (Color online) In the colicin version of the RPS game, strains that produce colicins (red/dark
grey) kill sensitive (green/light grey) strains, that outcompete resistant (blue/black) strains, that outcom-
pete colicin producing strains (toxin production involves bacterial suicide). Experiments [10] show that
colicin-producing strains cannot coexist with sensitive or resistant strains in a well-mixed culture, yet all
three phenotypes are recovered in natural populations. Two recent experiments have examined the role of
population structure in the maintenance of diversity among colicin-producing bacteria. In the first [10] in
vitro colonies were established on an agar substrate in Petri dishes, a setup which effectively limits compe-
tition to neighbors on the petri dish in analogy with explicit spatial embedding in 2D. In the second [12] in
vivo colonies were established in the intestines of co-caged mice, a setup which has two distinct scales of
mixing, with no explicit structure on either scale.
present in our approximation and can consequently (in terms of sensitivity to the details of the
topology) be considered more robust.
II. HIERARCHICAL MEANFIELD THEORY FOR TWO DISTINCT SCALES
Let us consider an evolutionary game between d types (strategies) described by the d×d payoff
matrix A with elements αkj . Assuming finite and constant population size, natural selection can
be described at the level of the individual by the so called the Moran process [30], during which at
each time step an individual is selected randomly from the population to be replaced (death) by the
offspring of an individual that is chosen proportional to its fitness to reproduce (birth). This models
a population in equilibrium, where the time scale of the population dynamics is set by the rate at
which ”vacancies” become available in the population. The fitness of each individual depends on
the payoff received from playing the game described by A with competitors (an individual of type
k receiving a payoff αkj when playing with an individual of type j). In well-mixed populations,
individuals can be considered to come into contact (compete) with equal probability with any
member of the population excluding themselves – this allows one to calculate the fitness of an
individual of type k in a meanfield manner, yielding
πk = πbase +
αkj(nj − δkj)
N − 1
, (1)
where nk is the number of individuals of type k in the population,
k=1 nk = N is the size of
the population, πbase is some baseline fitness and the Kronecker delta symbol δkj is equal to unity
if k = j and is zero otherwise. From this we may calculate the transition probabilities of our
stochastic process, i.e., the probability of an individual of type i being replaced by an offspring of
an individual of type k is given by
Tik =
, (2)
where π̄ =
k=1 πknk/N . The state of any population is completely described by the frequency
of the different strategies xk = nk/N . Due to the normalization
k=1 xk = 1, the values of xk
are restricted to the unit simplex Sd [3]. For d = 2 this is the interval [0, 1], S3 is the triangle with
vertices {(1, 0, 0), (0, 1, 0), (0, 0, 1)} while S4 is a tetrahedron etc.
As Traulsen et al. have recently shown [24, 25] for sufficiently large, but finite populations
the above stochastic process can be well approximated by a set of stochastic differential equa-
tions combining deterministic dynamics and diffusion (population drift) referred to as Langevin
dynamics:
ẋk = ak(x) +
ckj(x)ξj(t), (3)
where the effective deterministic terms ak(x) are given by
ak(x) =
(Tjk − Tkj) = xk
πk(x)− π̄(x)
π̄(x)
, (4)
ckj(x) are effective diffusion terms, that can also be expressed in terms of the transition probabil-
ities as described in [25], and ξj are delta correlated 〈ξk(t)ξj(t′)〉 = δkjδ(t − t′) Gaussian white
noise terms. As N → ∞ the diffusion term tends to zero as 1/
N and we are left with the
modified replicator equation.
In the context of our hierarchical mixing model the topology of connections can be described
by two parameters, the populations size at the local scale of mixing N , and a second parameter
µ, which tunes the strength of global mixing relative to the local dynamics. We take into account
the second (global) scale of mixing – mixing among local populations – by introducing a modified
version of the Moran process. In the modified process a random individual is replaced at each
time step either with the offspring of an individual from the same population (local reproduction)
or with an individual from the global population (global mixing). This is equivalent to considering
the global population to be well-mixed at the scale of local populations.
Let us consider a global population that is composed of M local populations of size N . In each
local population vacancies become available that local reproduction and global mixing compete to
fill. In any local population l the probability of an individual of some type k filling a new vacancy
due to local reproduction must be proportional to the number of individuals of type k multiplied by
their fitness i.e. πlkn
k, where we consider π
k to be determined only by interactions with individuals
in the same local population according to equation (1). To describe the tendency of individuals of
some type k in local population l to contribute to global mixing we introduce the parameters σlk.
The choice of appropriate σlk depends on the details of the global mixing mechanism, for systems
where only the offspring of individuals mix globally it is proportional to the fitness of a given type,
while for mechanisms such as physical mixing, by e.g. wind or ocean currents, it may be identical
for each type. Irrespective of the details, however, the probability of an individual of some type
k filling in a new vacancy due to global mixing should be proportional to the global average of
the number of individuals of type k multiplied by their mixing tendency, which we denoted as
〈σknk〉 =
l=1 σ
k/M , and the strength of global mixing µ. These consideration lead to the
new transition probabilities:
T̂ lik =
k + µ〈σknk〉
k=1(π
+ µ〈σknk〉)
k + µ〈σknk〉
N(π̄l + µ ¯〈σ〉)
, (5)
where π̄l =
k/N and
¯〈σ〉 =
〈σknk〉/N .
We have found that the results presented below are qualitatively the same for both the fitness
dependent choice of σlk = π
k and the fitness independent choice of σ
k = 1. Therefore, in the fol-
lowing we restrict ourselves to the somewhat simpler fitness independent choice of σlk = 1, which
can be considered to correspond to some form of physical mixing mechanism. The transition
probabilities (5) then reduce to:
T̄ lik =
π̄l + µ
π̄l + µ
. (6)
We can see that after a vacancy appears either local reproduction occurs, with probability π̄l/(π̄l+
µ), or global mixing, with probability µ/(π̄l + µ). From (6) we may derive the Langevin equation
describing the coevolutionary dynamics of population l from the
ẋlk = âk(x
l, 〈x〉) +
ĉkj(x
l, 〈x〉)ξj(t), (7)
with the modified deterministic terms given by
âk(x
l, 〈x〉) = x
k(πk(x
l)− π̄(xl)) + µ(〈xk〉 − xlk)
π̄(xl) + µ
, (8)
where the vector 〈x〉 =
l=1 x
l/M with components 〈xk〉 =
l=1 x
k/M describes the fre-
quencies of the individual types in the global population and the diffusion terms ĉ(xl, 〈x〉) can be
expressed in terms of the modified transition probabilities T̂ lik as above.
Equations (7) describe the coevolutionary dynamics of the global population through the cou-
pled evolution of the {x1, . . . ,xM} local populations. In the limit of a large number of local popu-
lations (M → ∞) the distribution of the local populations over the space of population states (the
simplex Sd) is described by a density function ρ(x) that is normalized over Sd, i.e.,
ρ(x) = 1.
The time evolution of ρ(x) follows a d−1 dimensional advection-diffusion equation – the Fokker-
Planck equation corresponding to eq. (7):
ρ̇(x) = −∇{â(x, 〈x〉)ρ(x)}+ 1
b̂(x, 〈x〉)ρ(x)
, (9)
with the global averages 〈xk〉 =
xkρ(x) coupled back in a self-consistent manner into the
deterministic terms âk(x, 〈x〉) and the diffusion matrix b̂kj(x, 〈x〉) =
i=1 ĉki(x, 〈x〉)ĉij(x, 〈x〉).
For large local populations (N → ∞) the diffusion term vanishes as 1/N .
The above advection-diffusion equation (9) presents an intuitive picture of the coevolutionary
dynamics of the population at a global scale. We can see that local populations each attempt to
follow the trajectories corresponding to the deterministic replicator dynamics, while under the
influence of two additional opposing forces: (i) global mixing, which attempts to synchronize
local dynamics and (ii) diffusion resulting from finite population size effects, which attempts to
smear them out over the simplex. The strength of these forces are tuned by two parameters µ and
N , respectively.
If, further, the effects of synchronization are irrelevant, as for example in the case of populations
where selection is externally driven by independent environmental fluctuations, we may replace
the global population average with the time average of any single population. This is the approach
we used in our study of genetic mixing in bacteria [31].
During our numerical investigations we found solving the advection-diffusion equation (9) nu-
merically challenging, particularly in the N → ∞ limit. We resorted instead to solving the coupled
Langevin equations (7) for large M = 104 − 105 to simulate the time evolution of ρ(x).
III. COOPERATION IN POPULATIONS WITH HIERARCHICAL LEVELS OF MIXING
The evolution of cooperation is a fundamental problem in biology, as natural selection under
most conditions favors individuals who defect. Despite of this, cooperation is widespread in na-
ture. A cooperator is an individual who pays a cost c to provide another individual with some
benefit b. A defector pays no cost and does not distribute any benefits. This implies the payoff
matrix
b− c −c
 , (10)
where b is the benefit derived from playing with a cooperator while c is the cost for cooperation.
From the perspective of evolutionary game theory, which equates payoff with fitness, the apparent
dominance of defection is simply the expression of the fact that natural selection a priori selects
for fitness of individuals and not the fitness of groups.
Defection dominates cooperation in any well-mixed population [3]. Population structure in-
duced by spatial structure [4, 18] and more general networks of interactions [21, 22, 23]) has,
however, been found to facilitate the emergence and maintenance of cooperation. The mecha-
nism responsible, termed spatial, or more generally, network reciprocity[32] depends strongly on
biased infux
movement of global average
density of 
local populations
drift
deterministic dynamics
x = 0 x = 1
ρ(x, t)
ρ(x, t + ∆t)biased influx
 0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9  1
b/c = 95, 97, 100, 111, 125
 1000
 10  100  1000
d b/c
 60  70  80  90  100  110
FIG. 2: a In an infinitely large well-mixed population evolutionary dynamics is deterministic and leads
to the extinction of cooperators as average fitness monotonically declines. The only stable fixed point
corresponds to the point where the fraction of cooperators is zero (x = 0). To understand qualitatively the
mechanism favoring cooperation in hierarchically mixed populations let us consider some density of local
populations (ρ(x, t)) that is symmetric around its mean at time t. Due to global mixing all local populations
are being driven toward the global average. Due to the influx bias, however, populations with a lower than
average number of cooperators will be driven stronger (faster) than those on the other side of the average.
Examining the density of local populations at some time t+∆t, this results in a net movement of the global
average toward a larger fraction of cooperators. This is, of course, opposed by local reproduction that favors
an increase in the number of defectors. For the global average to keep moving toward a higher number of
cooperators and eventually to keep balance with local reproduction bias a density of local population with
finite width is needed over which the effect of the influx bias can exert itself. It is drift caused by local
population size that maintains this finite width, and this is the reason that the b/c threshold above which
cooperation dominates depends on local population size. b Stationary density of local populations ρ(x) for
different values of b/c with N = 100, µ = 0.1. c Transition toward a global dominance of cooperation for
µ = 10. (triangles), µ = 1 (crosses), µ = 0.1 (squares), µ = 0.01 (circles) with N = 100. The critical value
of b/c depends only weakly on µ changing by 20% over four orders of magnitude d Critical values of b/c
as a function of N for different values of µ (notation as before). The dashed line corresponds to b/c = N .
The critical b/c values were determined by numerically finding the inflection point of the transition curves.
M = 103 was used throughout.
the details of local topology. In particular, it seems that lattice like connectivity structures where
three-site clique percolation occurs [17] and more general interaction graphs where the degree of
nodes k does not exceed the ratio of benefit to cost (i.e. k < b/c) [22] are required for cooperation
to be favored.
Examining the effects of hierarchical mixing on the evolutionary dynamics of cooperation we
found that a sharp, but continuous transition leads to the dominance of cooperation as the benefit
to cost ratio becomes smaller then the local population size, i.e. b/c < N . If the benefit to cost
ratio is larger then the local population size the global population is dominated by defectors. The
mechanism leading to the dominance of cooperation arises due to the competition between local
reproduction and global mixing. In local populations with lower average fitness – larger number of
defectors – the influx of individuals from the global scale will be larger than in local populations
with higher average fitness (cf. eq. (6) where the relative strength of the two terms on the left hand
side depends on the sum of the average fitness of population l and µ). The crucial ingredient for
cooperation to be successful is population drift introduced by finite local population size. It is bi-
ased influx coupled with drift that can result in cooperation being favored in the global population
(Fig 2.).
IV. THE RPS GAME
To explore the effects of hierarchical mixing in the context of games with three strategies we
first turn to the case of the so called ”rock-paper-scissors” (RPS) game. In the original popular
version of the game two players are afforded the chance to simultaneously display either rock
(fist), paper (flat hand) or scissors (two fingers). If player one displays a flat hand while player
two displays a fist, player one wins as paper wraps rock. Similarly scissors cut paper, and rocks
smashes scissors. Several examples of this game have been found in nature (e.g. among lizards
[33] ), but it is bacteria that have received the most experimental and theoretical attention.
In ecology the often high diversity among microbial organisms in seemingly uniform envi-
ronments, referred to as the ”paradox of the plankton”, has been difficult to understand. Several
models based on spatially explicit game theoretical models have been proposed to explain this di-
versity [10, 11, 13, 14]. These models are all variants of the RPS game played by colicin producing
bacteria. Colicins are antibiotics produced by some strains of Echerichia coli. In experiments (see
Fig.1) typically three strains are used: colicin producing (C), sensitive (S) and resistant (R). The
(a) (b) ALLD
µ = 0
µ = 0.01
µ = 0.2
µ = 0
µ = 0.1
FIG. 3: a Deterministic replicator dynamics (the N → ∞ limit) of the symmetric RPS game consists of
neutrally stable orbits along which the product of the strategy frequencies xRxPxS is conserved. If global
mixing is present (µ > 0) local populations deviate from these neutral orbits toward the global average
〈x〉. Considering the simplest system with global mixing, that consisting of M = 2 local populations we
see that in the presence of global mixing population x1 and population x2 move toward each other, respec-
tively moving closer and further from the barycentre of the triangle until they become synchronized and
subsequently pursue a common orbit. For deterministic local dynamics (N → ∞) such synchronization
invariably occurs for any M if µ > 0 and typically converges to the barycentre of the simplex for suffi-
ciently homogeneous initial conditions. b The deterministic replicator dynamics of the repeated PD game is
markedly different from that of the RPS game in that the internal fixed point is unstable and in the absence
of global mixing only ALLD survives. Again turning to the simplest scenario with M = 2 we see that
if µ = 0 any pair of populations x1 and x2 (gray and black lines) converge to the to the ALLD corner.
As µ is increased above a critical value a second, stable configuration emerges: for a large subset of the
possible initial conditions (all, but the left most x2) we see that one of the populations (x1) converges to
ALLD , while the second (x1) approaches a limit cycle. If µ is increased further, the above configuration
ceases to be stable, the population which initially converges to ALLD (x1) is subsequently ”pulled out”
by global mixing, following which the two populations synchronize and are finally absorbed together in
ALLD. Simulations, however, show that synchronization may be avoided for M > 2 if µ is not too large.
coevolutionary dynamics of the three strains can be cast in terms of an RPS game, C strains kill S
strains, but are outcompeted, by R strains, because toxin production involves the suicide of bacte-
ria. The cycle is closed by S strains that outcompete R strains, because resistance requires mutant
versions of certain membrane protein, which are less efficient than the wild type [10]. Despite
the cyclic dynamics colicin-producing strains cannot coexist with sensitive or resistant strains in a
well-mixed culture, yet all three phenotypes are recovered in natural populations. Local dispersal
(modeled as explicit spatial embedding) has widely been credited with promoting the maintenance
of diversity in this system [10, 11, 13, 14].
In its most symmetric form the RPS game is described by the payoff matrix
0 −ǫ ǫ
ǫ 0 −ǫ
−ǫ ǫ 0
, (11)
and some πbase > ǫ. The dynamics of this game in an infinitely large well mixed population
consists of neutral orbits along which the product xRxPxS is conserved. For any finite N , however,
fluctuations lead to the inevitable extinction of all but one of the strategies [15]. Spatial population
structure can avert this reduction in diversity [10, 13] through the emergence of a stable fixed
point at the barycentre of the simplex . The effect of the gradual randomization of different lattice
topologies (where a small number of edges are randomly rewired) on the dynamics of the game
has also been investigated. A Hopf bifurcation leading to global oscillations was observed [34, 35]
as the fraction of rewired links was increased above some critical value.
Examining the dynamics of the symmetric RPS game in terms of our hierarchical meanfield
approximation we observed that an internal fixed point emerged for N → ∞ (Fig.3a). More
importantly, diversity was also maintained for finite local population sizes if global mixing was
present. Simulations of the time evolution of ρ(x) also revealed a Hopf bifurcation leading to
the oscillation of the global average as µ was increased above a critical value µc depending on
N (Fig.4a). These results show that previous results obtained from simulations of populations
constrained to different lattice topologies can be considered universal in the sense that not only
lattices, but any population structure that can be approximated by two distinct internally unstruc-
tured scales of mixing are sufficient for their existence. In the context of the ”paradox of the
plankton” these results imply that aside of local dispersal (modeled as explicit spatial embedding)
a minimal metapopulation structure (with local competition and global migration) can also facili-
tate the maintenance of diversity in cyclic competition systems.
 0.045  0.055 0.05  0.06  0.065  0.07
 0.05 0.07 0.09
b d f
µ = 0.05
 1600 1400 1200 1000 800 600 400
µ = 0.05
 0.16 0.12 0.08 0.04 0
 0.08
 0.06
 0.08 0.06 0.04 0.02
N = 1000
N = 2000
N = 4000
N = 1200
small N 
cyclic
TFTALLC
N = 1000
µ = 0.05µ = 0.15
N = 4000
FIG. 4: (Color online) a In the case of the rock-paper-scissors game a Hopf bifurcation similar to that
observed for populations evolving on gradually randomized lattices [34, 35] leads to the emergence of
global oscillations (the red line indicates the trajectory of 〈x〉) if µ is larger than a critical value µc(N) (see
video S1 [36]). The density ρ(x) is indicated with a blue color scale. b The ratio A of the area of the global
limit cycle and the area of the simplex is plotted as a function of µ for three different values of N . For the
repeated prisoner’s dilemma game the combination of finite local population size and global mixing µ > 0
can lead to a stationary solution (c) qualitatively similar to that observed for explicit spatial embedding e.
This state is characterized by a stable global average (large dot), just as the lattice system (data not shown)
and sustained local cycles of cooperation, defection and reciprocity, also similar to the lattice case where
groups of ALLD (red, dark grey) individuals are chased by those playing TFT (blue, black), which are
gradually outcompeted by ALLC (green, light grey). d As µ is decreased a discontinuous transition can be
observed to the ALLD phase. The ratio I of populations on the internal cycle is plotted as a function of
µ. The inset shows the transition for different values of N . f The same critical line in the µ-N plane can
be approached by increasing N with µ fixed. A large hysteresis can be observed as N is decreased below
the critical value indicating the discontinuous nature of the transition. We numerically simulated the time
evolution of ρ(x) by integrating the stochastic differential equation system defined by eq. (7) for large M
(104−105) throughout. For the RPS game we used πbase = 1 and ǫ = 0.5, while in the case of the repeated
PD game we followed ref. [38], setting T = 5, R = 3, P = 1, S = 0.1,m = 10 and c = 0.8. Lattice
simulations (e) where performed on 1000 × 1000 square lattice with an asynchronous local Moran process
between neighbors and periodic boundary conditions.
V. THE REPEATED PRISONER’S DILEMMA GAME
In the general formulation of the prisoner’s dilemma (PD) game, two players have the choice
to cooperate or to defect. Both obtain some payoff R for mutual cooperation and some lower
payoff P for mutual defection. If only one of the players defects, while the other cooperates,
the defector receives the highest payoff T and the cooperator receives the lowest payoff S. That
is T > R > P > S and defection dominates cooperation in any well-mixed population. New
strategies become possible, however if the game is repeated, and players are allowed to chose
whether to defect or cooperate based on the previous actions of the opponent. In the following
we consider, similar to refs. [37] and [38] that recently examined the role of finite population
size and mutation and finite population size, respectively in terms of the repeated PD game with
three strategies: always defect (ALLD), always cooperate (ALLC), and tit-for-tat (TFT). TFT
cooperates in the first move and then does whatever the opponent did in the previous move. TFT
has been a world champion in the repeated prisoner’s dilemma ever since Axelrod conducted his
celebrated computer tournaments [7], although it does have weaknesses and may be defeated by
other more complex strategies [39].
Previous results indicate that if only the two pure strategies are present (players who either
always defect or ones who always cooperate) explicit spatial embedding [4] and some sufficiently
sparse interaction graphs [22, 40] allow cooperation to survive and the behavior of populations is
highly sensitive to the underlying topology of the embedding [17]. We have found that introducing
global mixing into the PD game with only the two pure strategies present also allows cooperation
to survive. The mechanism responsible for favoring cooperation in this case, however, depends
on the details of the competition between local reproduction and global mixing. For more than
two strategies these details are much less relevant and do not qualitatively influence the dynamics.
We will, therefore, consider the delicate issues concerning the PD game with only the two pure
strategies in a separate publication, and concentrate here on the repeated PD game with three
strategies.
To investigate the effect of global mixing on the repeated PD game with three possible strate-
gies: ALLD, ALLC and TFT following Imhof et al. [38] we considered the payoff matrix:
large N 
cyclic
small N 
cyclic
TFTALLC
TFTALLC
TFTALLC
TFTALLC
N = 10
N = 10
µ = 0.05
µ = 0.05
(iii)
(iii)
deterministic local populations
(i) (iii) (i)(ii)
effective single population
N = 10
N = 10
µ = 0.15
µ = 0.15
(iii)
FIG. 5: (Color online) Phase space for the repeated prisoner’s dilemma game on a population structure
with two distinct scales (see video S2 [36]). Three different phases are possible depending on the values of
µ and N : (i) only ALLD survives (ii) an internal limit cycles is maintained by global mixing due to a large
density of local populations around the ALLD corner (iii) a globally oscillating self maintaining limit cycle
is formed. For extreme values of µ the global dynamics reduces to that of some well-mixed population
where only ALLD survives: As µ becomes negligible (µ ≪ πk for all k) we approach the limit of isolated
local populations, while for µ ≫ πk we are left with a single synchronized population. Similarly for N = 2
– the smallest system with competition – the system can be described as a single well mixed population for
any µ and ALLD again prevails. In the limit of deterministic local populations (N → ∞) all three phases
can be found depending on the value of µ. The density ρ(x) is indicated with the color scale. A figure
illustrating the phase space of the repeated prisoner’s dilemma game with fitness dependent global mixing
is included in the supplementary material [36].
Rm Sm Rm
Tm Pm T + P (m− 1)
Rm− c S + P (m− 1)− c Rm− c
, (12)
where the strategies are considered in the order ALLC, ALLD, TFT, m corresponds to the number
of rounds played and c to the complexity cost associated with conditional strategies (TFT). The
dynamics of this game has a single unstable internal fixed point and the state where each member
of the population plays ALLD is the only nontrivial stable equilibrium (Fig.3b).
Introducing global mixing, between local well-mixed populations, however, causes new sta-
tionary states to emerge . Three phases can be identified: (i) ALLD wins (ii) large fraction of local
populations in the ALLD corner maintains local cycles of cooperation defection and reciprocity
through providing an influx of defectors that prevent TFT players from being outcompeted by
ALLC playing individuals (iii) a self maintaining internal globally oscillating cycle emerges. The
simplest scenario of two (M = 2) deterministic (N → ∞) local populations coupled by global
mixing (µ > 0) already leads to the emergence of phase (ii) as demonstrated in Fig.3b while phase
(iii) only emerges for larger M . For larger M simulations show that in the limit of large local pop-
ulations all global configurations with less than some maximum ratio of the populations I on the
internal cycle are stable in phase (ii). A transition from phase (ii) to (i) happens as µ is decreased
below a critical value µii→ic and I approaches zero as I = (1−µii→ic /µ) (data not shown). This can
be understood if we considered that near the transition point a critical proportion C = µ(1− I) of
ALLD individuals needs to arrive to stabilize local cycles of cooperation defection and reciprocity.
At the critical point I = 0 and µ = µii→ic which implies C = µ
c giving I = (1− µii→ic /µ)
Exploring the N − µ phase space (Fig.5) we see that the transition from phase (i) to (ii) be-
comes discontinuous for finite N (Fig.4d,e). Further, for any given value of N and µ the global
configuration is described by a unique I due to the presence of diffusion. For appropriate values
of the parameters the global average converges to a stationary value in phase (ii) similarly to case
of explicit spatial embedding (Fig.4c).
For very small (µ ≪ πk for all k) and very large (µ ≪ πk) values of µ the global dynamics can
be reduced to that of some well-mixed population where only ALLD persists (Fig 5.). For small N
we again have an effective well-mixed population – the only limit were defectors do not dominate
is N → ∞. In comparison with previous results of Imhof et al. we can see that evolutionary
cycles of cooperation defection and reciprocity can be maintained not only by mutation, but also
by population structures with hierarchical levels of mixing.
VI. DISCUSSION
While it is, of course, clear that the reduction of any realistic population structure to a manage-
able construction is always an approximation, it has not been clearly established what the relevant
degrees of freedom are in terms of evolutionary dynamics. Meanfield approximations are a classic
method of statistical and condensed matter physics and are routinely used to circumvent intractable
combinatorial problems which arise in many-body systems. Cluster-meanfield approximations of
sufficient precision [18, 19] have been developed that adequately describe the evolutionary dynam-
ics of explicitly structured populations through systematically approximating the combinatorial
complexity of the entire topology with that of small motif of appropriate symmetry. The effects of
more minimal effective topologies have, however, not been investigated previously. In the above
we have shown that straightforward hierarchical application of the meanfield approximation (the
assumption of a well-mixed system) surprisingly unveils a new level of complexity.
In the broader context of ecological and population genetics research on structured populations
our model can be described as a metapopulation model. The term ’metapopulation’ is, however,
often used for any spatially structured population [27], and models thereof. More restrictive defi-
nitions of the term are often implied in the context of ecology and population genetics literature.
The foundations of the classic metapopulation concept where laid down by Levin’s vision of
a ”metapopulation” as a population of ephemeral local populations prone to extinction. A classic
metapopulation persists, like an ordinary population of mortal individuals, in a balance between
’deaths’ (local extinctions) and ’births’ (establishment of new populations at unoccupied sites)
[27]. This classic framework is most wide spread in the ecology literature, a less often employed
extension is the concept of a structured metapopulation where the state of the individual popula-
tions is considered in more detail, this is more similar to our concept of hierarchical mixing, but
differs in considering the possibility of local extinctions.
The effects of finite population size and migration, which our model considers, has been of
more central concern in the population genetics literature. The analog of Levin’s classic metapop-
ulation concept is often referred to as the ’finite-island’ model [28] the effective population genetic
parameters describing which, have been explored in detail[29]. The study of the population ge-
netics of spatially subdivided populations in fact predates Levin, Wright having emphasised the
capacity of drift in small populations to bring about genetic differentiation in the face of selection
and/or migration several decades prior[28].
Our hierarchical mixing model treats the coevolutionary dynamics of evolutionary games on
structured populations in a manner similar to the most simple population genetic models of spa-
tially subdivided populations, focusing on the parallel effects of selection, drift and migration. It
goes beyond these models both in considering the effects of frequency dependent selection (and
the strategic aspects of the evolutionary dynamics this implies) and in using a self-consistent ap-
proach to describe the global state of the subdivided population. Also, in order to maintain a
connection with previous work on the effects of spatial structure on evolutionary games, which
rely on Nowak’s concept of spatial games [4], with individuals restricted to interact, and hence
compete, only with neighbours as defined by some topology of interaction, we develop our model
from the level of the individual by introducing a modified version of the Moran process – and not
by extending the Wright-Fisher process (which considers discrete generations and binomial sam-
pling to account for finite population size). The effective population structure described by our
hierarchical mixing model can be thought of as a population of individuals, interactions among
which are specified by the edges of a hierarchically organized random graph. The fundamental
difference in our picture is that the edges of this graph of interactions are not considered to be
fixed, but are instead in a constant state of change, being present with a different probability be-
tween pairs of individuals who share the same local population and between pairs of individuals
who do not (Fig.1.). We consider annealed randomness, which in contrast to the usual quenched
picture of fixed edges is insensitive to the details of topology. Our approach we believe best facil-
itates the exploration of the effects of changing the relative strengths of drift and migration in the
context of evolutionary games on structured populations.
Examining the effects of hierarchical mixing in the context of the evolution of robustness we
demonstrated that biased influx coupled with drift can result in cooperation being favored, pro-
vided the ratio of benefit to cost exceeds the local population size. This result bears striking
resemblance to that of Ohtsuki et al. [22], who were able to calculate the fixation probability of a
randomly placed mutant for any two-person, two-strategy game on a regular graph and found that
cooperation is favored provided the ratio of benefit to cost exceeds the degree of the graph. Our
results demonstrate that this rule extends to the minimal spatial structure induced by hierarchical
levels of mixing.
Applying our model of spatial structure to the repeated prisoners dilemma revealed that a con-
stant influx of defectors can help to stabilize cycles of cooperation, defection, and reciprocity
through preventing the emergence of an intermittent period of ALLC domination in the popula-
tion, which would present a situation that ”leaves the door wide open” to domination by defectors.
While previous work has been done on the effects of ”forcing” cooperation [41] the idea that an
influx of defectors can in fact stabilize the role of reciprocity in promoting cooperation has not
been proposed previously. It seems highly unlikely that this mechanism can be explained in terms
of kin or multilevel (group) selection, the similarities between which in structured populations
have recently been the subject of intensive debate (see e.g. [42] and [43] or [44] and [45]). Kin
selection can operate whenever interactions occurring among individuals who share a more recent
common ancestor than individuals sampled randomly from the whole population [45] are relevant.
In our case it is the interaction between defectors, arriving from the global scale, and TFT players
present at the local scale that is important, and not the interaction between individuals in the local
population, who may be thought of as sharing a recent common ancestor due to local dispersal.
Also, while the concept of multilevel selection presents a promising framework for the study of
evolution of cooperation, it must nonetheless be possible to derive it from ”first principles” – just
as kin selection can be cast as an emergent effect of local dispersal.
While there has been considerable work on studying the evolutionary games on graphs and
highly symmetric spatial structures very little attention has been paid to the effects of more min-
imal effective population structures, despite their widespread application in ecology and popula-
tion genetics, fields from which evolutionary game theory was born and must ultimately reconnect
with. We believe that the minimal population structure that such a hierarchical meanfield theory
describes is potentially more relevant in a wide range of natural systems, than more subtle setups
with a delicate dependence on the details and symmetries of the topology. We showed through
two examples that such structure is sufficient for the emergence of some phenomena previously
only observed for explicit spatial embedding, demonstrating the potential of our model to identify
robust effects of population structure on the dynamics of evolutionary games that do not depend
on the details of the underlying topology. The practical advantage of our approach, lies in its
ability to readily determine whether or not some feature of a structured population depends on the
topological details of local interactions.
Recent simulation result concerning the dynamics of public goods games on different popula-
tion structures [9, 46] and experiments where global mixing in an RPS like bacteria-phage system
lead to the emergence of a ”Tragedy of the commons” scenario [47] should all be amicable to
analysis in terms of our method.
VII. ACKNOWLEDGMENTS
This work was partially supported by the Hungarian Scientific Research Fund under grant No:
OTKA 60665.
VIII. APPENDIX
Our approach readily generalizes for an arbitrary number of hierarchical mixing levels. For
three levels of mixing we may consider the global population to be comprised of M subpopula-
tions each of which is in turn subdivided into M local populations. With m ∈ {1, · · · ,M} running
over subpopulations and l ∈ {1, · · · ,M} over local populations the transition probabilities can be
written as:
T̂mlik =
πmlk n
k + µ
(1)〈(1)σml′k 〉l′ + µ(2)〈〈(2)σm
k 〉l′〉m′
k=1(π
+ µ(1)〈(1)σml′
〉l′ + µ(2)〈〈(2)σm
〉l′〉m′)
, (13)
where primed indices indicate the scale of mixing over which the average is taken, µ(1) describes
the strength of mixing, and the (1)σmlk the tendencies of mixing among local populations within a
subpopulation, while µ(2) describes the strength of mixing, and the (2)σmlk the tendencies of mixing
among subpopulations in the global population.
[1] Fisher R.A. The Genetical Theory of Natural Selection (Oxford University Press, 1930).
[2] Maynard Smith, J. & Price, G. Nature (London) 246, 15-18 (1930).
(doi:10.1038/246015a0)
[3] Hofbauer, J. & Sigmund, K. Evolutionary Game Theory and Population Dynamics (Cambridge Uni-
versity Press, 1998).
[4] Nowak, M.A. & May, R.M. Nature (London) 359, 826-829 (1992).
(doi:10.1038/359826a0)
[5] Lieberman, E., Hauert, C. & Nowak, M.A. Nature (London) 433, 312-316 (2005).
(doi:10.1038/nature03204)
[6] Axelrod, R. & Hamilton, W.D. Science 211, 1390-1396 (1981).
(doi:10.1126/science.7466396)
[7] Axelrod, R. The Evolution of Cooperation (Basic Books, New York, 1984).
[8] Nowak, M.A., Bonhoeffer S. & May, R.M. Proc. Natl. Acad. Sci. USA 91, 4877-4881 (1994).
(doi:10.1073/pnas.91.11.4877)
[9] Szabó, Gy. & Hauert, C. Phys. Rev. Lett. 89, 118101 (2002).
(doi:10.1103/PhysRevLett.89.118101)
[10] Kerr, B., Riley, M.A., Feldman, M.W. & Bohannan, B.J. Nature (London) 418, 171-174 (2002).
(doi:10.1038/nature00823)
[11] Nowak, M.A. & Sigmund, K. Nature (London) 418, 138-139 (2002).
(doi:10.1038/418138a)
[12] Kirkup, B.C. & Riley, M.A. Nature (London) 428, 412-414 (2004).
(doi:10.1038/nature02429)
[13] Czárán, T.L., Hoekstra, R.F. & Pagie, L. Proc. Natl. Acad. Sci. USA 99, 786-790 (2002).
(doi:10.1073/pnas.012399899)
[14] Lenski, R.E. & Riley, M.A. 99, 556-558 (2002).
(doi:10.1073/pnas.02264199)
[15] Reichenbach, T., Mobilia, M. & Frey, E. Phys. Rev. E 74, 051907 (2006).
(doi:10.1103/PhysRevE.74.051907)
[16] Szabó, Gy. & Tőke, Cs. Phys. Rev. E 58, 69 (1998).
(doi:10.1103/PhysRevE.58.69)
[17] Szabó, Gy., Vukov, J. & Szolnoki, A. Phys. Rev. E 72, 047107 (2005).
(doi:10.1103/PhysRevE.72.047107)
[18] Hauert, C. & Szabo, Gy. Am. J. Phys. 73, 405-414 (2005).
(doi:10.1119/1.1848514)
[19] Hui, C. & McGeochb, M.A. Bull. Math. Biol. 69, 659-676 (2007).
(doi:10.1007/s11538-006-9145-1)
[20] Hui, C., Zhang, F., Han, X. & Lid, Z. Ecol. Mod. 184, 397-412 (2005).
(doi:10.1016/j.ecolmodel.2004.11.004)
[21] Santos, F.C. & Pacheco, J.M. Phys. Rev. Lett. 95, 098104 (2005).
(doi:10.1103/PhysRevLett.95.098104)
[22] Ohtsuki, H., Hauert, C., Lieberman, E. & Nowak, M.A.
Nature (London) 441, 502-505 (2006).
(doi:10.1038/nature04605)
[23] Wang, S., Szalay, M.S., Zhang, C. & Csermely, P.
PLoS ONE 3 e1917 (2008).
(doi:10.1371/journal.pone.0001917)
[24] Traulsen, A., Claussen, J.C. & Hauert, C. Phys. Rev. Lett. 95, 238701 (2005).
(doi:10.1103/PhysRevLett.95.238701)
[25] Traulsen, A., Claussen, J.C. & Hauert, C. Phys. Rev. E 74, 011901 (2006).
(doi:10.1103/PhysRevE.74.011901)
[26] Hanski, I.A. & Gilpin, M.E. in Metapopultaion Biology (Academic Press, 1997), pp. 93-108.
[27] Hanski, I.A. Nature 396, 41 (1998).
(doi:10.1038/23876)
[28] Pannell, R.P. & Charlesworth B. Philosophical Transactions: Biological Sciences 355, 1851-1864
(2000).
(doi:10.1098/rstb.2000.0740)
[29] Whitlock M.C. & Barton N.H. Genetics 146, 427-41 (1997).
[30] Moran, P.A.P. The Statistical Processes of Evolutionary Theory (Clarendon, 1962)
[31] Szöllősi, G.J., Derényi, I. & Vellai, T. Genetics 174, 2173-2180 (2006).
(doi:10.1534/genetics.106.063412)
[32] Nowak, M.
Science 314, 1560-1563 (2006).
(doi:10.1126/science.1133755)
[33] Sinervo, B. & Lively, C.M. Nature (London) 380, 240-243 (1996).
(doi:10.1038/380240a0)
[34] Szolnoki, A. & Szabó, Gy. Phys. Rev. E 70, 037102 (2004).
(doi:10.1103/PhysRevE.70.037102)
[35] Szabó, Gy., Szolnoki, A. & Izsák, R. J. Phys. A 37, 2599-2606 (2004).
(doi:10.1088/0305-4470/37/7/006)
[36] See EPAPS Document No. [number will be inserted by publisher] for videos S1,S2 and S3.
[37] Nowak M.A., Sasaki, A., Taylor, C. & Fudenberg, D. Nature (London) 428, 646-650 (2004).
(doi:10.1038/nature02414)
[38] Imhof, L.A., Fudenberg, D. & Nowak, M.A. Proc. Natl. Acad. Sci. USA 102, 10797-10800 (2005).
(doi:10.1073/pnas.0502589102)
[39] Molander, P. J. Conflict Resolut. 29, 611-618 (1985).
(doi:10.1177/0022002785029004004)
[40] Taylor, P.D., Day T. & Wild G. Nature (London) 447, 469-472 (2007).
(doi:10.1038/nature05784)
[41] Szabó, Gy., Antal, T., Szabó, P. and Droz, M. Phys. Rev. E 62, 1095-1103 (2000).
(doi:10.1103/PhysRevE.62.1095)
[42] Killingback, T., Bieri, J., & Flatt, T. Proc. R. Soc. B. 273 1477-1481 (2006).
(doi:10.1098/rspb.2006.3476)
[43] Grafen, A. Proc. R. Soc. B 274, 713-719 (2007).
(doi:10.1098/rspb.2006.0140)
[44] Traulsen, A. & Nowak, M.A. Proc. Natl. Acad. Sci. USA 103, 10952-10955 (2007).
(doi:10.1073/pnas.0602530103)
[45] Lehmann L., Keller, L., West, S. & Roze, D. Proc. Natl. Acad. Sci. USA 104, 6736-6739 (2007).
(doi:10.1073/pnas.0700662104)
[46] Hauert, C., Silvia De Monte, S. , Hofbauer, J. & Sigmund K. Science 296, 1129-1132 (2002).
(doi:10.1126/science.1070582)
[47] Kerr, B., Neuhauser, C., Bohannan, B.J.M. & Dean A.M. Nature (London) 442, 75-78 (2006).
(doi:10.1038/nature04864)
	Introduction
	 Hierarchical Meanfield Theory for Two Distinct Scales
	Cooperation in populations with hierarchical levels of mixing
	The RPS Game
	The Repeated Prisoner's Dilemma Game
	Discussion
	Acknowledgments
	Appendix
	References
