Long-range prisoner's dilemma game on a cycle
Jiwon Bahk Seung Ki Baek Hyeong-Chai Jeong







truemm

1/2
#1
#1

 2 
 4 

    =   =   =̱ =  
    P^0_s_k→_CP_s_k→_Cmember i member j p__Lp__Rp__Cp__Dx__Ccthcth,Ncthidcthnnccth,Nnncc̃thkeffr__ijr__kjζ(α)ζ(2α)ζ_N/2(α)ζ(β)ζ(2β)ζ_N/2(2β)ζ_N/2(β)ζ(α+β)ζ(α)ζ(β)f(α,β)g(α,β)#1#2
        < g r a p h i c s >
    #1#2
        < g r a p h i c s >
    Π̃]Long-range prisoner's dilemma game on a cycleDepartment of Physics and Astronomy, Sejong University, Seoul
    05006, Korea[]seungki@pknu.ac.krDepartment of Physics, Pukyong National University, Busan 48513,
    Korea[]hcj@sejong.ac.krDepartment of Physics and Astronomy, Sejong University, Seoul
    05006, Korea
    We investigate evolutionary dynamics of altruism with long-range
    interaction on a cycle. The interaction between individuals
    is described by a simplified version of the prisoner's dilemma (PD) game
    in which the payoffs are parameterized by c, the cost of a cooperative
    action.
    In our model, the probabilities of the game interaction and competition
    decay algebraically with r_AB, the distance between two players A and B,
    but with different exponents: That is, the probability to play the PD game is
    proportional to r_AB^-α. If player A is chosen for death, on the
    other hand, the probability for B to occupy the empty site is
    proportional to r_AB^-β.
    In a limiting case of β→∞, where the competition for an empty site
    occurs between its nearest neighbors only,
    we analytically find the condition for the proliferation
    of altruism in terms of , a threshold of c below which altruism
    prevails. For finite β, we conjecture a formula for  as a function
    of α and β. We also propose a numerical method to locate ,
    according to which we observe excellent agreement with the conjecture
    even when the selection strength is of considerable magnitude.
    [
        Hyeong-Chai Jeong
        December 30, 2023
    =====================
    
    
    
    § INTRODUCTION
    
    
    As Charles Darwin remarked in  On the Origin of Species,
    “Natural selection will never produce in a being anything injurious
    to itself, for natural selection acts solely by and for the good of each.”
    The existence of altruism is seemingly at odds with this statement, as far as
    altruists are defined as those who benefit others at the cost to themselves.
    For this reason,
    there have been extensive studies to explain the evolution of altruism.
    The explanations can be categorized (arguably) into five rules, i.e., direct
    reciprocity, indirect reciprocity, spatial reciprocity, kin selection,
    and group
    selection <cit.>.
    An essential common factor in these mechanisms is the
    assortment of individuals carrying the cooperative
    phenotype <cit.>.
    Let us consider a population of defectors.
    If there emerges a group within which the members cooperate
    themselves,
    its size will increase because the members earn
    higher fitness than the population average. Furthermore,
    if the assortment can be maintained while the size of the group grows,
    the whole population will become cooperative.
    On the other hand, if the assortment breaks, making
    every individual interact with every other, cooperators
    will disappear because they earn less than the population average by
    construction.
    
    Considering that the defectors always win in the mean-field regime,
    one may expect that spatial reciprocity becomes weaker as the
    dimensionality of the space d increases. The logic goes as
    follows: In a population of size N, the average chemical distance
    between a pair of individuals behaves as ∼ N^1/d and thus
    decreases dramatically as d grows. Hence, it is reasonable to expect
    a relatively homogeneous state in a high-dimensional structure.
    In many statistical-physical models, the mean-field theory becomes
    exact as d exceeds a certain value d_c<cit.>.
    Another way to look at the crossover between a low-dimensional spin
    system and the mean-field limit is to consider long-range interactions
    on a one-dimensional (1D) lattice 
    structure <cit.>. 
    For example, when the interaction strength decays as r^-1-σ
    between a pair of 1D Ising spins with distance r,
    the spins show the mean-field-like critical behaviors for
    0 < σ≤ 1/2 as those with short range 
    interactions in high dimensions d>d_c=4 do.   
    For 1/2 ≤σ≤ 1, 1D Ising spin system shows the
    order-disorder transitions at finite temperatures but
    its critical behaviors are different from mean-field transition.
    For σ>1, the interactions become short-ranged effectively
    and they behave as the 1D spins with finite-range interactions.  
    
    The story becomes more complicated when it comes to
    evolutionary games: Let us consider
    the Prisoner's Dilemma (PD) game with the death-birth (DB)
    reproduction process in the weak selection limit. It has a certain population
    structure, and the game is played between players at one step away.
    On regular structures, the condition for cooperation turns out to be
    determined by the number of nearest neighbors regardless
    of d <cit.> (see also
    Ref. allen2017evolutionary for more general results).
    Furthermore, spatial reciprocity breaks down for any population structures if
    we employ the birth-death (BD) process instead.
    What is the difference between the DB and BD processes? The point is that
    we actually have another kind of interaction, that is,
    competition for reproduction. For social behavior to evolve, these two should
    have different length scales <cit.>.
    Even if the game interaction shares the same underlying structure with
    competition for reproduction, the DB process naturally separates the scales
    because the competition for reproduction effectively occurs between players at
    two steps away.
    
    In this paper, we wish to investigate the roles of these two length scales
    explicitly by considering the PD game on a cycle with long-range
    interactions, whose dynamics is governed by the DB process.
    The game is parameterized by the cost of cooperation, denoted by
    c. We are interested in its threshold, ,
    above which cooperators are outnumbered by defectors in the
    infinite population with zero mutation limit. 
    The length scale of the game interaction is determined by an exponent
    α, so that the game result between two individuals at distance r
    has a weighting factor of r^-α. On the other hand, the length scale of
    competition for reproduction is determined by another exponent β,
    so the probability for an individual to occupy an empty site at distance r
    decays algebraically as r^-β.
    Our first main result is to derive  for β→∞, where
    only the nearest neighbors compete for reproduction. Then, we
    consider the case of α = β, where an analytic result for 
    is available in the weak-selection limit <cit.>.
    Based on these two special cases, we conjecture a
    general form of (α, β). It agrees well with
    numerical results.
    
    This work is organized as follows: In the next section, we give a detailed
    description of our model. In Sec. <ref>, we present  as a
    function of α in case of β→∞ and review the exact 
    when α=β in the weak-selection limit. Based on these two results,
    we conjecture  for general α and β. Section <ref>
    compares this conjecture with Monte Carlo data. The conjectured threshold nicely
    matches with our numerical calculation, and this holds even if the selection
    strength increases. We discuss this finding and conclude this work in
    Sec. <ref>.
    
    
    
    § MODEL
    
    The population structure is a cycle, i.e., a 1D lattice with the
    periodic boundary condition. For the sake of convenience, the total
    number of sites is assumed to be an add number in this section and we
    set N ≡ 2n+1. Of course, it will be irrelevant
    whether N is even or odd if we deal with a sufficiently large system.
    Every site i is occupied by a player, who is also denoted by
    i=-n,-n+1,…, -1, 0, 1, …, n-1, n, and each player i can
    choose an action s_i between cooperation C and defection D. The PD game is
    formulated as the donation game, where a cooperator benefits the co-player by an
    amount of b > 0 with paying cost c ∈ (0,b). The payoff matrix π is
    thus defined as
    [                                           [ C D ];                 [ C; D ] ( [ b-c  -c;   b   0 ]). ]
    Without loss of generality, we will set b=1 throughout this work.
    The expected payoff of player i per game is given by the following
    weighted sum, 
    Π_i
     = ∑_j≠ i^-α/2ζ_n (α)π_s_i s_j=̱1/2ζ_n (α)∑_j≠ iπ_s_i s_j/^α
    where ζ_n(α) ≡∑_k=1^n k^-α is the
    incomplete Riemann zeta function and
    ^-α/2ζ_n(α)
    is the probability that player i play a game with player j
    where  is the distance between them. On a cycle,
     is given by the minimum of
    |i-j|, |i-j+N|, and |i-j-N|. 
    We consider the expected payoff per game (rather than the total
    payoff, ∑_j≠ i^-απ_s_i s_j) 
    to make the average payoff independent of the population size. 
    This normalized payoff is also better to compare two systems with
    different α fairly. The systems with smaller α would
    show effectively more strong selections without the normalization.
    
    
    
    Fig. <ref>(a) illustrates the way to calculate payoff
    Π_0 of player 0 for α=2 and n=3. 
    Player 0 at the top plays the PD game with
    every other in the population, k=±3, ±2, and k=±1. 
    From Eq. (<ref>), we see that π_0k is
    -c for k=-2, -1 and 1 and b-c for k=-3, 3, and 2.
    Π_0 is the weighted sum of these values with the wighting factors,  
    |k|^-2/2ζ_3(2) where
    2 ζ_3(2) = 21/1^2+ 1/2^2+1/3^2= 98/36.  
    Since |k|^-2/2ζ_3(2) is
    36/98 for n=±1,
    9/98 for n=±2, 
    and 4/98 for n=±3,
    Π_0 is given by
    π_0= 17/98b - c.
    We can calculate the payoffs of the other players similarly.
    They are 
    Π_-3 = 49/98b - c,
    Π_-2 = 58/98b,
    Π_-1 = 53/98b,
    Π_1 = 85/98b,
    Π_2 = 54/98b - c,
    Π_3 = 49/98b - c.
    
    Player i's fitness is given by the exponential form of the payoff,
    
    f_i = e^ωΠ_i,
    
    where the selection strength ω represents how strongly the
    fitness depends on the payoff. Exponential fitness of Eq. (<ref>)
    becomes the usual linear fitness of 1+ωΠ_i in the
    weak selection limit. It has an advantage over linear fitness in the
    strong selection limit. We can interpret the exponential fitness
    as fecundity directly even when the payoff is negative. Furthermore,
    DB dynamics is invariant under the addition of a constant to all 
    elements in the pay off matrix which certainly preserve the Nash
    equilibrium.   
    
    We employ a long-range version of the DB process as a model of
    reproduction. The first step is to choose a death site k randomly,
    irrespective of fitness, and make it empty. The probability to choose
    a particular site k is therefore inversely proportional to N. 
    The difference from the usual DB process for the nearest-neighbor
    competition is that any player i in the population can produce an
    offspring at the empty site. For given k, the conditional
    probability p_ik that player i's offspring takes over the site
    k is proportional to the player's fitness f_i but decreases as 
    r_ik^-β where r_ik is the distance between them.
    In other words, p_ik is given by
    
     p_ik= 1/Z_ke^ωΠ_i/r_ik^β,
     
    where Z_k ≡∑_i≠k  e^ωΠ_i× r_ik^-β is a
    normalization factor. Note that the player at the death site k cannot leave an
    offspring. In particular, s_k will become C with probability
    = 1/Z_k ∑_i≠ke^ωΠ_i/r_ik^β δ_s_i,C,
    
    where δ_x,y is the Kronecker delta.
    To remove trivial absorbing states where every player takes the same
    action, we introduce a small mutation probability μ. 
    When player i reproduces its offspring at the site k,
    s_k becomes s_i with probability 1-μ and by the opposite of s_i
    with probability μ. If we combine the DB process and mutation
    together, the actual probability that s_k becomes C is given 
    by
    =  1-μ + μ1-.
    
    It is straightforward to see that the complementary event s_k → D
    occurs with probability 1-. 
    
    We demonstrate the reproduction process in Fig. <ref>(b) for
    β=2 and n=3. 
    When a player, say, i=0, is chosen for death, the
    other players k≠0 compete to produce their offspring at this
    empty site. Player k's winning probability is given by
    1/|k|^2 Z_0  e^ωΠ_k with
    Π_k given in the caption of (a). Here
    Z_0 = e^ωΠ_1+e^ωΠ_-1
    +1/4e^ωΠ_2+e^ωΠ_-2
    +1/9e^ωΠ_3+e^ωΠ_-3
    so that the normalization factor 1/Z_0 makes
    the total sum of probabilities be 1. The winner's offspring inherits
    the winner's strategy with probability (1-μ) or the opposite one
    with μ≪ 1. 
    
    
    
    § ANALYSIS
    
    
    We begin by considering how the threshold cost (α, β)
    behaves in special cases. The first one is the case with the
    nearest-neighbor competition, represented by β→∞.
    The second one is when
    α = β, and (α, α) in the weak-selection limit
    is readily obtained by applying the result of
    Ref. allen2017evolutionary.
    
    
    
     §.§ Nearest-neighbor competition (β→∞)
    
    
    If β→∞, only the nearest neighbors of the death site will compete
    for reproduction. Let us define a domain as a set of consecutive sites with the
    same action, bounded by sites with the opposite action. When a death site k is
    picked up inside a domain, therefore, the action of its neighbors must be the
    same as s_k. If μ is negligibly small, we can say that the configuration
    { s_i } will not change. It implies that the system shows
    dynamics only at domain boundaries. Let us assume that the system is divided
    into two equally large domains so that
    
        s_i =
        {[                  C                for i=0, -1, -2, …, -n;                  D                for      i=1, 2, …, n, ].
    
    as shown in Fig. <ref>(a). Assume that k=1 is chosen as for death.
    Only its neighboring players k=0 and k=2 compete to take over the empty
    site, and the domain wall will move right if the offspring of player i=0
    wins. The payoffs Π_0 and Π_2 are calculated as
    Π_0
    = 1/2ζ_n(α)∑_j≠0π_s_0 s_j/r_0j^α= 1/2ζ_n(α)∑_r=1^n1-c/r^α+-c/r^α=̱1/2 - c, 
    Π_2
    = 1/2ζ_n(α)∑_j≠2π_s_2
    s_j/r_2j^α= 1/2ζ_n(α)∑_r=2^n1/r^α=̱1/2 - 1/2ζ_n(α),
    
    and the probability that the wall moves right is given by
    
    P^0_s_1 → C= e^ωΠ_0/e^ωΠ_0+e^ωΠ_2.
    
    Similarly, we can calculate the probability that the wall moves left,
    P^0_s_0 → D, from Π_-1 and Π_1.
    The population will become cooperative when
    P^0_s_1 → C > P^0_s_0 → D.
    If N ≫ 1, this inequality reduces to
    
    c < 1/2ζ_n(α)
    and we have
    (α, ∞)=1/2ζ_n(α).
    
    
    
     §.§ Analytic solution of α = β
    
    
    If the game interaction and competition share the same population
    structure, strategies selected by the network reciprocity can be
    obtained on any population structure in the weak-selection limit of
    ω→ 0 using
    coalescent theory <cit.>.
    When we denote the expected coalescence
    time <cit.>
    from the two ends of an n-step random walk
    by τ_n, cooperation is favored when c is smaller than
    
        = τ_3-τ_1/τ_2
    
    in the weak-selection limit.
    In a weighted undirected graph, if site i is connected to j by an
    edge, the edge has a weight w_ij = w_ji. The weighted degree of
    site i is the sum of its edge weights:
    w_i ≡∑_j w_ij. The probability for a random 
    walker to transit from i to j is denoted by p_ij≡ w_ij / w_i.
    A useful quantity is the Simpson degree
    κ_i = ( ∑_j p_ij^2 )^-1,
    a generalization of the topological degree.
    For a regular graph, the site index i is unnecessary since
    every vertex has the same Simpson degree k_i=k. 
    Equation (<ref>) is then rewritten as
    
        = N/κ - 2/N - 2.
    
    In our case, noting that w_ij∝ r_ij^-α for i≠ j and
    w_ii = 0, we find the threshold as a combination of the
    zeta functions <cit.>,
    
        (α, α) = ζ_n(2α)/2ζ_nα^2,
    
    where we have taken into account N ≫ 1 in the sense of the
    ω N-limit <cit.>.
    
    
    
     §.§ Conjecture for (α, β)
    
    To recap, we have two predictions:
    First, if β→∞ and N ≫ 1,
    the threshold cost is predicted to be
    
        (α, ∞) =1/2ζ_n(α),
    
    so that ζ_n(α) can be interpreted as the effective
    number of neighbors to play the game with on one side of the focal player.
    Second, when α = β,
    the threshold cost  has been derived as
    
        (α, α) =
        ζ_n(α + α)/2ζ_n(α)ζ_n(α)
    
    for ω→ 0. To incorporate Eqs. (<ref>) and (<ref>),
    we conjecture a connection formula (α, β) as follows:
     (α, β) =
    ζ_n(α + β)/2 ζ_n(α)ζ_n(β).
    
    In particular, for α→∞, the conjecture predicts
    
        (∞, β) = 1/2ζ_n(β).
    
    
    § NUMERICAL CALCULATION
    
    
    For Monte Carlo simulation, we construct a cycle of size N with a player at
    each site. Initial s_i's are randomly chosen between C and D equally
    probably. Then, according to the DB process, we randomly select a player (say,
    k) for death and choose another player i with probability p_ik in
    Eq. (<ref>).
    One Monte Carlo step consists of N trials to update the configuration
    {s_i }. In addition, we introduce mutation with small probability
    μ≪ 1 so that the system has a unique stationary
    distribution over configurations { s_i }. The quantity of
    interest is cooperator abundance , the average frequency of cooperators in
    the steady state. For each data point in a sample, we run the Monte Carlo
    simulation for 6× 10^3 Monte Carlo steps, and the time average of 
    is taken for the last 10^2 Monte Carlo steps, at which the simulation has
    already reached a steady state. The total number of samples amounts to 4×
    10^3.
    
    
    
     §.§ Weak selection (w=0.02)
    
    
    Noting that most of the analytic results have been obtained in the
    weak-selection limit, we have to check this limit first and proceed to stronger
    selection strengths thereafter.
    Let us begin by setting β→∞. If α also goes to infinity,
    the threshold at which () = 1/2 is analytically given as
    (∞, ∞) = 1/2 in the ω N-limit <cit.>.
    However, this limit cannot be taken exactly in the Monte Carlo simulation,
    and μ cannot be infinitesimally small either. For c=1/2,
    we actually observe = 0.495 ± 0.001. We regard this value as a
    reference point to estimate (α,∞), assuming that it takes into
    account all the corrections due to the difference from the ω N-limit as
    well as the finiteness of μ.
    The procedure to estimate (α,∞) now goes as follows: For given
    α, we calculate  as a function of c, together with the standard
    error from statistically independent samples. As a lower bound of (α,
    ∞), we pick up c_min, defined as
    c below which all the error bars of (c) lie above 0.495 ±
    0.001. Likewise, the upper bound of (α, ∞) is given by
    c_max, defined as c above which all the error bars of (c) lie
    below 0.495 ± 0.001 [Fig. <ref>(a)]. The numerical estimate of (α, ∞) is thus represented by an interval [c_min, c_max] as
    shown in Fig. <ref>(b). The agreement with Eq. (<ref>) is
    compelling.
    The above procedure can be applied to the case of α = β as well
    [Fig. <ref>(a)]. Here, the estimation is more difficult
    because the slope of  becomes steeper as α decreases. Nevertheless,
    Fig. <ref>(b) shows that our estimates are
    consistent with the analytic solution [Eq. (<ref>)].
    
    
    
    Let us move on to check whether the conjecture of Eq. (<ref>) correctly
    predicts the threshold behavior even when no analytic proof is available.
    One specific example that we can test is an extreme case of α→∞,
    where the conjecture reduces to Eq. (<ref>).
    Comparing Fig. <ref>(a) with Fig. <ref>(a),
    we see that α and β do not play the same role in determining the
    cooperator abundance. Nevertheless, Fig. <ref>(b) substantiates the
    conjectured symmetry between α and β in the threshold behavior.
    To make both α and β finite yet unequal, we consider another case
    of α = 2β [Fig. <ref>(a)]. Again, our numerical results
    support the conjecture in Eq. (<ref>). The overall behavior might look
    similar to the previous case of α=∞, but a closer look shows that
    our numerical estimates are sharp enough to distinguish the difference between
    them [Fig. <ref>(b)].
    
    
    
     §.§ Stronger selection (w=1 and w=5)
    
    
    The conjectured formula [Eq. (<ref>)] works even beyond the
    weak-selection limit. Our numerical estimation procedure hinges on the fact that
    we can still take c=1/2 as the threshold in a double limit of
    α→∞ and β→∞ for arbitrary selection
    strength <cit.>. In this
    limit, our numerical calculation gives x_c = 0.2923 ± 0.0045 at c=1/2,
    provided that we have set w=1 [Fig. <ref>(a)]. As
    we have done above, we will take this value as a reference point to locate the
    threshold in other combinations of α and β for which one can
    suitably define an extrapolation to α→∞ and β→∞:
    Let us take a look at the results shown in Fig. <ref>(b). Here, we
    demonstrate three cases, α = ∞, α=2β, and α=β.
    Note that the limit can be approached by taking β→∞ for all these
    cases, so that the abundance levels may be compared with the reference value in
    the double limit. Although none of the cases have been solved analytically, the
    threshold behavior is indeed captured by Eq. (<ref>) with high
    precision.
    
    As w increases, the quality of fit deteriorates, however.
    For w=5, for example, the reference value is lowered to x_c = 0.0280 ±
    0.0007 [Fig. <ref>(a)]. It is a tiny value, and the slope of (c) is
    steep, which makes the estimate more challenging than the previous
    weaker-selection cases. It is not surprising that deviations from the conjecture
    become visible [Fig. <ref>(b)], but the overall agreement is still
    remarkable.
    
    
    
    § CONCLUDING REMARKS
    
    
    To summarize, we have studied the PD game on a cycle, where the
    probabilities for two players at distance r to play the game and to
    compete for reproduction decay algebraically as r^-α and
    r^-β, respectively. When the benefit of cooperation is fixed
    as unity, there exists a threshold cost (α, β) below
    which the population becomes cooperative. Based on analytically
    tractable cases [Eqs. (<ref>) and (<ref>)], we have
    conjectured its functional form as given in Eq. (<ref>).
    The apparent symmetry between α and β in Eq. (<ref>) is
    nontrivial as one sees that Fig. <ref>(a) clearly
    differs from <ref>(a).
    Although we are unable to prove our conjecture at the moment,
    the striking agreement with numerical results suggests that the exact
    result of c(α,α) can be generalized in a simple form, at least
    in this particular long-ranged interaction model on a cycle.
    
    The threshold costs in our model do not seem to depend on the
    selection strength ω strongly. For the nearest neighbor competitions
    (β→∞), this independence can be understood easily since the
    threshold cost is obtained from the domain boundary dynamics
    as discussed in III A. For finite β, our conjecture relies on the
    analytic solution of weak selection limit. However, the threshold
    costs are almost independent of selection strength up to
    ω≈ 5.
    Further studies on the robustness of the weak selection
    limit as well as the threshold costs for general 
    population structures are needed.
    
    
    
    S.K.B. was supported by Basic Science Research Program through the
    National Research Foundation of Korea (NRF) funded by the Ministry of Science,
    ICT and Future Planning (NRF-2017R1A1A1A05001482). H.C.J. was
    supported by Basic Science Research Program through the National
    Research Foundation of Korea (NRF) funded by the Ministry of Education
    (NRF-2018R1D1A1A02086101).
    
    
    
    
    
    
    
    
    
    
    
    16
    fxundefined [1]
     ifx#1
    fnum [1]
      #1firstoftwo
     secondoftwo
     
    fx [1]
      #1firstoftwo
     secondoftwo
     
    
    
    
    
    
    noop [0]secondoftwo
    
    ref[1]@startlink#1@href
    href[1]#1@endlink
    anitize@url [0] `
    12 `$12
      `&12 `#12 `1̂2 `_12 `%12
    startlink[1]
    endlink[0]
    
    rl [1]href #1
    
    
    
    
    
    
    
    
    
    
    
    
    @bib@innerbibempty
    
     [Nowak(2006)]nowak2006five
      authorauthorM. A. Nowak, @noop journaljournalScience volume314, pages1560
      (year2006)NoStop
     [Nowak et al.(2010)Nowak,
      Tarnita, and Wilson]nowak2010evolution
      authorauthorM. A. Nowak, authorC. E. Tarnita,
       and authorE. O. Wilson,@noop journaljournalNature volume466, pages1057 (year2010)NoStop
     [Abbot et al.(2011)Abbot,
      Abe, Alcock, Alizon,
      Alpedrinha, Andersson, Andre,
      Van Baalen, Balloux, Balshineet al.]abbot2011inclusive
      authorauthorP. Abbot, authorJ. Abe, authorJ. Alcock, authorS. Alizon, authorJ. A. Alpedrinha, authorM. Andersson, authorJ.-B.Andre, authorM. Van Baalen, authorF. Balloux,
      authorS. Balshine,  et al., @noop journaljournalNature volume471, pagesE1
      (year2011)NoStop
     [Marshall(2011)]marshall2011group
      authorauthorJ. A. Marshall, @noop journaljournalTrends Ecol. Evol. volume26, pages325 (year2011)NoStop
     [Fletcher and Doebeli(2009)]fletcher2009simple
      authorauthorJ. A. Fletcher and authorM. Doebeli, @noop journaljournalProc.
      Royal Soc. B volume276, pages13
      (year2009)NoStop
     [Park and Jeong(2012)]park2012emergence
      authorauthorS. Park and authorH.-C. Jeong, @noop journaljournalJ.
      Korean Phys. Soc. volume60, pages311 (year2012)NoStop
     [Débarre et al.(2014)Débarre, Hauert, and Doebeli]debarre2014social
      authorauthorF. Débarre, authorC. Hauert,
       and authorM. Doebeli,@noop journaljournalNat. Commun.volume5, pages3409 (year2014)NoStop
     [Landau and Lifshitz(1980)]landau1980statistical
      authorauthorL. D. Landau and authorE. M. Lifshitz, @noop titleStatistical Physics,
      Part 1 (publisherPergamon Press, addressOxford, year1980)NoStop
     [Thouless(1969)]thouless1969long
      authorauthorD. Thouless, @noop journaljournalPhys. Rev. volume187, pages732
      (year1969)NoStop
     [Fisher et al.(1972)Fisher,
      Ma, and Nickel]fisher1972critical
      authorauthorM. E. Fisher, authorS.-k. Ma,   and authorB. Nickel,@noop journaljournalPhys. Rev. Lett.volume29, pages917 (year1972)NoStop
     [Luijten(2000)]luijten2000monte
      authorauthorE. Luijten, in @noop booktitleComputer
      Simulation Studies in Condensed-Matter Physics XII (publisherSpringer, year2000) pp. pages86–99NoStop
     [Ohtsuki et al.(2006)Ohtsuki, Hauert, Lieberman, andNowak]ohtsuki2006simple
      authorauthorH. Ohtsuki, authorC. Hauert,
      authorE. Lieberman,  andauthorM. A. Nowak, @noop
      journaljournalNature volume441, pages502 (year2006)NoStop
     [Allen et al.(2017)Allen,
      Lippner, Chen, Fotouhi,
      Momeni, Yau, and Nowak]allen2017evolutionary
      authorauthorB. Allen, authorG. Lippner,
      authorY.-T. Chen, authorB. Fotouhi, authorN. Momeni, authorS.-T. Yau,  and authorM. A. Nowak, @noop journaljournalNature volume544, pages227 (year2017)NoStop
     [Wakeley(2008)]wakeley2009coalescent
      authorauthorJ. Wakeley, @noop titleCoalescent Theory: An
      Introduction (publisherRoberts & Co., addressGreen Village, Colorado, year2008)NoStop
     [Hardy and Wright(1979)]hardy1979introduction
      authorauthorG. H. Hardy and authorE. M. Wright, @noop titleAn Introduction to the
      Theory of Numbers, edition5th ed. (publisherOxford Univ. Press, year1979) p. pages254NoStop
     [Jeong et al.(2014)Jeong,
      Oh, Allen, and Nowak]jeong2014optional
      authorauthorH.-C. Jeong, authorS.-Y. Oh,
      authorB. Allen,  and authorM. A. Nowak, @noop journaljournalJ. Theor. Biol. volume356, pages98 (year2014)NoStop
