#  Evolutionary game with stochastic payoffs in a finite island model

Dhaker Kroumi1111Author for correspondence, and e-mail:
<EMAIL_ADDRESS>and Sabin Lessard2
1Department of Mathematics and Statistics
King Fahd University of Petroleum and Minerals
Dhahran 31261, Saudi Arabia
2Department of Mathematics and Statistics
University of Montreal
Montreal H3C 3J7, Canada

## Abstract

In this paper, we consider a two-player two-strategy game with random payoffs
in a population subdivided into $d$ demes, each containing $N$ individuals at
the beginning of any given generation and experiencing local extinction and
recolonization with some fixed probability $m$ after reproduction and
selection among offspring. Within each deme, offspring engage in random
pairwise interactions, and the payoffs are assumed to have means and variances
proportional to the inverse of the population size. By verifying the
conditions given in Ethier and Nagylaki (1980) to approximate Markov chains
with two time scales, we establish that the discrete-time evolutionary
dynamics with $Nd$ generations as unit of time converges to a continuous-time
diffusion as $d\rightarrow\infty$. The infinitesimal mean and variance of this
diffusion are expressed in terms of the population-scaled means and variances
of the payoffs besides identity-by-descent measures between offspring in the
same deme in a neutral population. We show that the probability for a strategy
to fix in the population starting from an initial frequency $(Nd)^{-1}$
generally increases as the payoffs to that strategy exhibit less variability
or the payoffs to the other strategy more variability. As a result,
differences in variability can make this fixation probability for cooperation
larger than the corresponding one for defection. As the deme-scaled extinction
rate $\nu=mN$ decreases for $N$ large enough and $m$ small enough, creating a
higher level of identity among offspring within demes, the differences between
the population-scaled variances of the payoffs for interacting offspring of
different types increases this effect to a greater extent than the differences
for interacting offspring of the same type. Furthermore, assuming no
variability in the payoffs to cooperation, we explore in more details the
relationships between the levels of variability in the payoffs to defection as
functions of $\nu$ that are necessary for the fixation probability of
cooperation to be above $(Nd)^{-1}$, the corresponding fixation probability of
defection to be below $(Nd)^{-1}$, and the former fixation probability to
exceed the latter.

Keywords and phrases: Prisoner’s Dilemma; Random payoffs; Fixation
probability; Identity by descent; Island Model; Diffusion approximation.

Mathematics Subject Classification (2010): Primary 92D25; Secondary 60J70

##  1 Introduction

Evolutionary game dynamics is a process of frequency-dependent selection, in
which the reproductive success of an individual, also called its fitness, is
determined by the payoff that the individual receives according to the
strategy used by the individual in interactions with others. Individuals with
higher fitnesses are more likely to reproduce and pass on their strategies
(Maynard Smith [36], Maynard Smith and Price [37], Hofbauer and Sigmund [15],
Nowak [44]). The reproductive success of an individual is not determined by
its absolute fitness but rather by its relative fitness compared to the other
individuals in the population. This framework can be used to shed light on the
emergence of cooperation or competition in natural populations (Axelrod and
Hamilton [3], Sugden [52]).

Early studies of evolutionary games have primarily focused on an infinite
well-mixed population, using a straightforward differential equation known as
the replicator equation which assumes that the growth rate of a strategy is
given its average payoff (Taylor and Jonker [54], Zeeman [61]). However, real
biological species are finite in number, and their dynamics cannot be
adequately captured by a deterministic differential equation. In order to
address this fact, investigations have turned towards finite well-mixed
populations as modeled by finite Markov chains (see, e.g., Nowak et al. [43],
Taylor and Wild [59], Imhof and Nowak [17], Antal et al. [1]). These
stochastic processes play a crucial role in comprehending the emergence and
maintenance of cooperative behaviors, the influence of random drift, and the
impact of population size on evolutionary dynamics. The Moran model (Moran
[39]) assumes that one individual at a time is replaced by another following a
birth-death event, while the Wright-Fisher model (Fisher [13], Wright [60])
assumes non-overlapping generations with a binomial scheme from one generation
to the next. More general exchangeable selection models from one time step to
the next extending the neutral Cannings model (Cannings [4]) have also been
considered (Lessard and Ladret [33]).

Well-mixed populations oversimplify the reality of biological species, which
are more complex. Evidence of geographic structure in genetic data (Slatkin
[51]) highlights the importance of studying populations subdivided into
different groups subject to specific environmental and ecological effects.
Evolutionary game dynamics in structured populations explore the dynamics of
strategy frequencies among interacting individuals in populations that exhibit
spatial or social structure. In such populations, individuals are not randomly
mixed but instead interact primarily with their neighbors or with individuals
in close geographic proximity. These structural contraints on interactions can
significantly influence the evolutionary outcomes such as the emergence of
cooperation.

The island model, which serves as a basic representation of a structured
population, is commonly used to investigate the impact of limited migration,
local adaptation, and population subdivision on the process of evolution.
Ethier and Nagylaki [11] have established conditions for a discrete-time
Markov chain with two time scales to converge weakly to a diffusion process in
the limit of a large population. The strong migration limit as the deme size
tends to infinity and the number of demes is kept fixed was studied in
Nagylaki [40, 41, 42], while Slatkin [50] and Takahata [53] considered the
island model under conditions of low migration.

Numerical simulations have been conducted by Cherry [5, 6] and Cherry and
Wakeley [7] to evaluate the effectiveness of diffusion methods in scenarios
involving a large number of small demes. A similar approach using diffusion
methods was employed by Whitlock [58] to investigate the stepping-stone model
of population structure. This model assumes a population subdivided into demes
of fixed size connected by migration. Roze and Rousset [48] have proposed a
general method for constructing diffusion approximations in structured
population genetics models. They achieved this by using approximations for the
expectation and variance in allele frequency change over one generation, based
on partial derivatives of fitness functions with respect to the trait value
and probabilities of genetic identity under neutrality. Wakeley [58], Wakeley
and Takahashi [59] and Lessard [32] have studied an island model consisting of
a finite number of demes, all of the same finite size in the limit of a large
number of demes, using a diffusion approximation based on the conditions in
Ethier and Nagylaki [11]. This is an ideal framework to study the effects of
kin selection and group selection.

All the aforementioned studies assume a constant environment over time.
However, this assumption is unrealistic as environmental conditions can
randomly fluctuate, impacting competition abilities and birth-death rates
(Kaplan et al. [18], Lande et al. [30], May [38]). These fluctuations can
influence population size and composition, leading researchers to investigate
their effects. For instance, Lambert [29], Parsons and Quince [46, 47], and
Otto and Whitlock [45] have studied the fixation probability of a mutant type
in an unstructured population with dynamic size fluctuations due to
demographic scenarios of growth or decline. Conversely, Uecker and Hermisson
[55] have focused on a single beneficial allele in a population that
experiences temporal variation in both size and selection pressure.

A random payoff matrix in pairwise interactions provides an alternative
approach for examining fluctuations in the environment, reflecting the
stochastic nature of real-world interactions. It encompasses variations in the
environment, the presence of noise or decision-making errors, and the inherent
unpredictability of outcomes.

Recently, the effect of stochastic changes in payoffs in discrete time were
studied with particular attention to stochastic local stability of fixation
states and constant polymorphic equilibria in an infinite population (Zheng
_et al._ [63, 64]), the fixation probability in a large population that
reproduces according to a Wright-Fisher model (Li and Lessard [34]) or in a
finite well-mixed population that reproduces according to a Moran model
(Kroumi et al. [27, 26]). Moreover, Kroumi and Lessard [24, 25] explored the
effects of random payoffs on the average abundance of strategies such as
cooperation and defection under recurrent mutation in finite well-mixed
populations using the Moran model for reproduction. The general findings
reveal that the inclusion of uncertainty in the payoffs received by defectors
or the reduction of uncertainty in the payoffs received by cooperators can
facilitate the evolution of cooperation.

The objective of this paper is to examine the influence of stochastic
variability in a population structured into a large number of isolated demes
of the same finite size. Our analysis will focus on exploring how the
interplay between the second moments in payoffs and the migration rate affects
the evolutionary dynamics of cooperation, specifically its probability of
becoming fixed under different scenarios.

The structure of this paper is as follows: Section 2 presents the island model
with random pairwise interactions within demes and random payoffs for a two-
player game, while Section 3 establishes the diffusion approximation in the
limit of a large number of demes based on the two time scales for the changes
in strategy frequencies and deme type frequencies. In Sections 4 and 5, we
develop the expressions of the infinitesimal mean and variance in the
diffusion approximation and we examine conditions for the evolution of a
strategy to be favoured by selection based the fixation probability following
its introduction as a single mutant. Sections 6-9 investigate these conditions
under various scenarios on the variances of the payoffs. Finally, in Section
10, we discuss our results and their relationship with the existing
literature. Technical details are relegated to Appendices A-J.

##  2 Model

Consider a haploid population that is structured into $d$ isolated demes.
Generations are discrete, non-overlapping, and every generation starts with
$N$ individuals in each deme. These produce a large number of offspring which
is assumed to be the same for each individual in the population. This is
followed by random pairwise interactions between the offspring produced in the
same deme.

Each individual in the population can be of type $S_{1}$ or type $S_{2}$, and
an offspring inherits the type of its parent. In the case where an offspring
of type $S_{1}$ interacts with another offspring of type $S_{1}$, both
offspring receive a payoff $a_{1,1}$. If an offspring of type $S_{1}$
interacts with an offspring of type $S_{2}$, the former receives a payoff
$a_{1,2}$ while the latter receives $a_{2,1}$. Finally, when an offspring of
type $S_{2}$ interacts with another offspring of type $S_{2}$, both offspring
receive a payoff $a_{2,2}$.

The payoffs are given by the entries of the game matrix

$\bordermatrix{&S_{1}&S_{2}\cr S_{1}&a_{1,1}&a_{1,2}\cr
S_{2}&a_{2,1}&a_{2,2}\cr}=(a_{k,l}).$ (1)

In this paper, we suppose that the payoffs randomly fluctuate over successive
generations in such a way that $|a_{k,l}|<1$ and

$\displaystyle
E\left[a_{k,l}\right]=\frac{\mu_{k,l}}{Nd}+o\left(\frac{1}{d}\right),$ (2a)
$\displaystyle
E\left[a_{k,l}^{2}\right]=\frac{\sigma^{2}_{k,l}}{Nd}+o\left(\frac{1}{d}\right),$
(2b) $\displaystyle E\left[|a_{k,l}|^{n}\right]=o\left(\frac{1}{d}\right),$
(2c)

for $k,l=1,2$ and $n\geq 3$. The parameters $\mu_{k,l}$ and $\sigma^{2}_{k,l}$
are the population-scaled mean and variance of the payoffs, respectively, for
$k,l=1,2$. In addition, we suppose that $a_{k_{1},l_{1}}$ and
$a_{k_{2},l_{2}}$ are independent random variables as long as
$(k_{1},l_{1})\not=(k_{2},l_{2})$.

Note that the payoff received by an offspring $I$ can be written as

$a_{I}=\sum_{k,l=1}^{2}a_{k,l}q_{k,I}q_{l,I^{\prime}},$ (3)

where $I^{\prime}$ represents a randomly chosen offspring produced in the same
deme as $I$, while $q_{k,l}$ and $q_{l,I^{\prime}}$ denote the frequencies
($0$ or $1$) of $S_{k}$ and $S_{l}$, respectively, in $I$ and $I^{\prime}$,
respectively, for $k,l=1,2$. This payoff is added to the reference value $1$
to obtain the viability of $I$, given by

$\omega_{I}=1+a_{I}.$ (4)

This quantity, always positive, is proportional to the survival probability of
$I$ from conception to maturity.

Note that the average viability of all offspring produced by an individual can
be interpreted as the reproductive success, or fitness, of this individual.
Moreover, the average viability of all offspring produced in a deme is
proportional to the number of offspring in the deme after viability selection.
It is assumed that a fixed proportion of these offspring go to a migrant pool,
while the other proportion stay in the deme where they were produced.

On the other hand, each deme goes to extinction with a fixed probability
$m\in(0,1)$ and, when this occurs, the deme is repopulated by a large number
of offspring coming from the migrant pool. Note that the contributions of the
demes to the migrant pool being proportional to the average viabilities within
demes, this is a case of hard selection (Christiansen [8]). Moreover, the
limiting case $m=1$ corresponds to a well-mixed population as in Hilbe [14].

Finally, the individuals to start the next generation are obtained by sampling
at random $N$ offspring within each deme. Note that the number of individuals
of type $S_{1}$ to start a deme can vary from $0$ to $N$, leading to a total
of $N+1$ possible deme states. A deme is considered to be of type $i$ if it
starts with $i$ individuals of type $S_{1}$ and $N-i$ individuals of type
$S_{2}$, for $i=0,1,\ldots,N$. Note that the initial frequency of $S_{1}$ in a
deme of type $i$ is given by

$\displaystyle x_{i}=\frac{i}{N},$ (5)

and this is also the frequency of $S_{1}$ among all offspring produced in the
deme, for $i=0,1,\ldots,N$.

Let $Z_{i}(t)$ be the fraction of demes of type $i$ at the beginning of
generation $t\geq 0$. Then, the population state is described by the frequency
vector $\mathbf{Z}(t)=(Z_{0}(t),Z_{1}(t),\ldots,Z_{N}(t))$, where
$\sum_{i=0}^{N}Z_{i}(t)=1$. Note that

$X(t)=\sum_{i=0}^{N}Z_{i}(t)x_{i}$ (6)

is the frequency of type $S_{1}$ in the whole population at the beginning of
generation $t\geq 0$.

Suppose that the initial population state is
$\mathbf{Z}(0)=\mathbf{z}=(z_{0},z_{1},\ldots,z_{N})$ with corresponding
frequency of type $S_{1}$ in the population given by

$X(0)=x=\sum_{i=0}^{N}z_{i}x_{i}.$ (7)

Following reproduction and random pairwise interactions, the average payoffs
received by an offspring of type $S_{k}$ and by offspring of all types in a
deme of type $i$ are

$\displaystyle\tilde{a}_{k,i}$
$\displaystyle=x_{i}a_{k,1}+\left(1-x_{i}\right)a_{k,2},$ (8a)
$\displaystyle\tilde{a}_{i}$
$\displaystyle=x_{i}\tilde{a}_{1,i}+\left(1-x_{i}\right)\tilde{a}_{2,i},$ (8b)

respectively, and the corresponding average viabilities are

$\displaystyle\tilde{\omega}_{k,i}$ $\displaystyle=1+\tilde{a}_{k,i},$ (9a)
$\displaystyle\tilde{\omega}_{i}$
$\displaystyle=x_{i}\tilde{\omega}_{1,i}+\left(1-x_{i}\right)\tilde{\omega}_{2,i}=1+\tilde{a}_{i},$
(9b)

respectively, for $k=1,2$ and $i=0,1,\ldots,N$. Therefore, the frequency of
type $S_{1}$ among offspring after viability selection in a deme of type $i$
can be expressed as

$\tilde{x}_{i}=\frac{x_{i}\tilde{\omega}_{1,i}}{\tilde{\omega}_{i}},$ (10)

for $i=0,1,\ldots,N$.

Note that the weighted frequency of $S_{1}$ among all offspring after
viability selection is

$\tilde{x}=\sum_{i=0}^{N}z_{i}\tilde{x}_{i}=\sum_{i=0}^{N}z_{i}x_{i}\frac{\tilde{\omega}_{1,i}}{\tilde{\omega}_{i}}$
(11)

if all demes have the same weight, but

$\tilde{\tilde{x}}=\sum_{i=0}^{N}z_{i}\tilde{x}_{i}\frac{\tilde{\omega}_{i}}{\tilde{\omega}}=\sum_{i=0}^{N}z_{i}x_{i}\frac{\tilde{\omega}_{1,i}}{\tilde{\omega}}$
(12)

where

$\tilde{\omega}=\sum_{i=0}^{N}z_{i}\tilde{\omega}_{i}=1+\sum_{i=0}^{N}z_{i}\tilde{a}_{i}=1+\tilde{a},$
(13)

if the weights of the demes are given by the average viabilities of the
offspring within demes. The latter is the frequency of $S_{1}$ in the migrant
pool. Therefore, this is the frequency of $S_{1}$ among offspring after
recolonization in a deme of type $i$ that goes extinct, which occurs with
probability $m$. In such a case, the number of offspring of type $S_{1}$ drawn
at random in that deme to start the next generation follows a binomial
distribution with parameters $N$ and $\tilde{\tilde{x}}$. Otherwise, with
probability $1-m$, the parameters are $N$ and $\tilde{x}_{i}$, for
$i=0,1,\ldots,N$.

We conclude that a deme that was of type $i$ at the beginning of generation
$0$ has a probability

$P_{ij}(\mathbf{z})=(1-m)\binom{N}{j}\left(\tilde{x}_{i}\right)^{j}\left(1-\tilde{x}_{i}\right)^{N-j}+m\binom{N}{j}\tilde{\tilde{x}}^{j}\left(1-\tilde{\tilde{x}}\right)^{N-j}$
(14)

to become a deme of type $j$ at the beginning of generation $1$, for
$i,j=0,1,\ldots,N$.

##  3 Diffusion approximation

In this section, we will derive a diffusion approximation that relies on the
existence of two timescales in the stochastic process. The first timescale is
long and pertains to the changes in the individual type frequencies. The
second timescale, that will be described later on, is short and relates to the
changes in the deviations of the deme type frequencies from their equilibrium
values in a neutral population subdivided into an infinite number of demes.

Let $Z_{ij}(1)$ be the proportion of demes of type $j$ at the beginning of
generation $1$ that were of type $i$ at the beginning of generation $0$. Note
that the variables $Z_{i_{1}j_{1}}(1)$ and $Z_{i_{2}j_{2}}(1)$ are independent
for all $i_{2},j_{2}=0,1,\ldots,N$ as long as $i_{1}\neq i_{2}$. Moreover,
given $\mathbf{Z}(0)=\mathbf{z}$ and the payoff matrix $(a_{k,l})$, the vector
$dz_{i}\Big{(}Z_{i0}(1),Z_{i1}(1),\ldots,Z_{iN}(1)\Big{)}$ follows a
conditional multinomial distribution with parameters $dz_{i}$ and
$\Big{(}P_{i0}(\mathbf{z}),P_{i1}(\mathbf{z}),\ldots,P_{iN}(\mathbf{z})\Big{)}$,
for $i=0,1,\ldots,N$. Therefore, we have

$\displaystyle
E_{\mathbf{z}}\left[Z_{ij}(1))\Big{|}(a_{k,l})\right]=P_{ij}(\mathbf{z}),$
(15a) $\displaystyle
Var_{\mathbf{z}}\left[Z_{ij}(1))\Big{|}(a_{k,l})\right]=\frac{P_{ij}(\mathbf{z})\left(1-P_{ij}(\mathbf{z})\right)}{dz_{i}},$
(15b) $\displaystyle
Cov_{\mathbf{z}}\left(Z_{ij_{1}}(1)),Z_{ij_{2}}(\mathbf{z}))\Big{|}(a_{k,l})\right)=-\frac{P_{ij_{1}}(\mathbf{z})P_{ij_{2}}(\mathbf{z})}{dz_{i}},$
(15c)

for $j_{1}\not=j_{2}$. Here, we use $E_{\mathbf{z}}$, $Var_{\mathbf{z}}$, and
$Cov_{\mathbf{z}}$ to represent the conditional expectation, variance, and
covariance, respectively, given that $\mathbf{Z}(0)=\mathbf{z}$.

Note that, for the fraction of demes of type $j$ at the beginning of
generation $1$ given by

$Z_{j}(1)=\sum_{i=0}^{N}z_{i}Z_{ij}(1),$ (16)

we have

$\displaystyle
E_{\mathbf{z}}\left[Z_{j}(1)\Big{|}(a_{k,l})\right]=\sum_{i=0}^{N}z_{i}P_{ij}(\mathbf{z}),$
(17a) $\displaystyle
Var_{\mathbf{z}}\left(Z_{j}(1)\Big{|}(a_{k,l})\right)=\frac{1}{d}\sum_{i=0}^{N}z_{i}P_{ij}(\mathbf{z})\left(1-P_{ij}(\mathbf{z})\right),$
(17b) $\displaystyle
Cov_{\mathbf{z}}\left(Z_{j_{1}}(1),Z_{j_{2}}(1)\Big{|}(a_{k,l})\right)=-\frac{1}{d}\sum_{i=0}^{N}z_{i}P_{ij_{1}}(\mathbf{z})P_{ij_{2}}(\mathbf{z}),$
(17c)

for $j,j_{1},j_{2}=0,1,\ldots,N$ with $j_{1}\neq j_{2}$.

The first three propositions below address changes in the individual type
frequencies. Throughout the paper, $I$ and $J$ will be used to designate
randomly chosen offspring produced in demes chosen at random and independently
before viability selection and recolonization. For offspring produced in the
same deme chosen at random, the notation $I_{1}$, $I_{2}$, etc. will be used,
and similarly $J_{1}$, $J_{2}$, etc. for offspring produced in the same deme
chosen at random and independently from the former deme. Moreover, we
introduce the population-scaled parameters

$\displaystyle\mu_{I_{1}}$
$\displaystyle=\sum_{k,l=1}^{2}q_{k,I_{1}}q_{l,I_{2}}\mu_{k,l},$ (18a)
$\displaystyle\sigma_{I_{1},I_{2}}$
$\displaystyle=\sum_{k,l=1}^{2}q_{k,I_{1}}q_{k,I_{2}}q_{l,I_{3}}q_{l,I_{4}}\sigma^{2}_{k,l},$
(18b) $\displaystyle\sigma_{I_{1},J_{1}}$
$\displaystyle=\sum_{k,l=1}^{2}q_{k,I_{1}}q_{k,J_{1}}q_{l,I_{2}}q_{l,J_{2}}\sigma^{2}_{k,l}.$
(18c)

On the other hand, note that

$\displaystyle E_{\mathbf{z}}\left(q_{1,I_{1}}\right)=x,$ (19a) $\displaystyle
Var_{\mathbf{z}}(q_{1,I_{1}})=x-x^{2},$ (19b) $\displaystyle
Cov_{\mathbf{z}}\left(q_{1,I_{1}},q_{1,I_{2}}\right)=\sum_{i=0}^{N}z_{i}x^{2}_{i}-x^{2},$
(19c)

where $I_{1}$ and $I_{2}$ are two randomly chosen offspring produced in the
same deme at the beginning of generation $0$.

Our first result concern the expected change in the frequency of of $S_{1}$.

######  Proposition 1.

The first conditional moment of the change in the frequency of $S_{1}$ in the
whole population from the beginning of generation $0$ to the beginning of
generation $1$ is

$E_{\mathbf{z}}\left[X(1)-X(0)\right]=\frac{1}{Nd}M(\mathbf{z})+o\left(d^{-1}\right),$
(20)

where

$\displaystyle M(\mathbf{z})$
$\displaystyle=Cov_{\mathbf{z}}(q_{1,I_{1}},\mu_{I_{1}})+(1-m)\Big{[}Cov_{\mathbf{z}}(q_{1,I_{1}},\sigma_{I_{2},I_{3}})-Cov_{\mathbf{z}}(q_{1,I_{1}},\mu_{I_{2}})-Cov_{\mathbf{z}}(q_{1,I_{1}},\sigma_{I_{1},I_{2}})\Big{]}$
$\displaystyle\quad-mCov_{\mathbf{z}}(q_{1,I_{1}},\sigma_{I_{1},J_{1}}).$ (21)

In order to prove this proposition, we will use two technical lemmas whose
proofs are relegated to Appendices A and B.

######  Lemma 1.

For the uniformly weighted frequency of $S_{1}$ after viability selection
within demes given in Eq. (11), we have

$\displaystyle E_{\mathbf{z}}\left[\tilde{x}\right]$
$\displaystyle=x+\frac{1}{Nd}\Big{[}Cov_{\mathbf{z}}(q_{1,I_{1}},\mu_{I_{1}})+Cov_{\mathbf{z}}(q_{1,I_{1}},\sigma_{I_{2},I_{3}})-Cov_{\mathbf{z}}(q_{1,I_{1}},\mu_{I_{2}})$
$\displaystyle\quad\quad\quad\quad\quad-
Cov_{\mathbf{z}}(q_{1,I_{1}},\sigma_{I_{1},I_{2}})\Big{]}+o(d^{-1}).$ (22)

######  Lemma 2.

For the frequency of $S_{1}$ in the migrant pool given in Eq. (12), we have

$E_{\mathbf{z}}\left[\tilde{\tilde{x}}\right]=x+\frac{1}{Nd}\Big{[}Cov_{\mathbf{z}}(q_{1,I_{1}},\mu_{I_{1}})-Cov_{\mathbf{z}}(q_{1,I_{1}},\sigma_{I_{1},J_{1}})\Big{]}+o\left(d^{-1}\right).$
(23)

###### Proof of Proposition 1.

Note first that the probability distribution given by
$(P_{ij}(\mathbf{z}))_{j=0}^{N}$ is a mixture of a binomial distribution with
parameters $N$ and $\tilde{x}_{i}$ and a binomial distribution with parameters
$N$ and $\tilde{\tilde{x}}$ whose weights are $1-m$ and $m$, respectively.
Therefore, we have

$\sum_{j=0}^{N}jP_{ij}(\mathbf{z})=(1-m)N\tilde{x}_{i}+mN\tilde{\tilde{x}}.$
(24)

Using this identity and Eq. (17), we find that

$\displaystyle E_{\mathbf{z}}\left[X(1)-X(0)\Big{|}(a_{k,l})\right]$
$\displaystyle=\sum_{j=0}^{N}x_{j}E_{\mathbf{z}}\left[Z_{j}(1)\Big{|}(a_{k,l})\right]-x$
$\displaystyle=\sum_{i=0}^{N}z_{i}\sum_{j=0}^{N}x_{j}P_{ij}(\mathbf{z})-x$
$\displaystyle=(1-m)\tilde{x}+m\tilde{\tilde{x}}-x.$ (25)

Taking the expected value and using Lemmas 1 and 2 lead to Proposition 1. ∎

Our second result concerns the variance of the change in the frequency of
$S_{1}$.

######  Proposition 2.

The conditional variance of the change in the frequency of $S_{1}$ in the
whole population from the beginning of generation $0$ to the beginning of
generation $1$ is

$Var_{\mathbf{z}}\left[X(1)-X(0)\right]=\frac{1}{Nd}Q(\mathbf{z})+o\left(d^{-1}\right),$
(26)

where

$\displaystyle Q(\mathbf{z})$
$\displaystyle=Var_{\mathbf{z}}(q_{1,I_{1}})+(1-m)(Nm-1)Cov_{\mathbf{z}}(q_{1,I_{1}},q_{1,I_{2}})+m^{2}Cov_{\mathbf{z}}(q_{1,I_{1}}q_{1,J_{1}},\sigma_{I_{1},J_{1}})$
$\displaystyle\quad-2mx\,Cov_{\mathbf{z}}(q_{1,I_{1}},\sigma_{I_{1},J_{1}})+(1-m^{2})Cov_{\mathbf{z}}(q_{1,I_{1}},q_{1,J_{1}}\sigma_{I_{1},J_{1}})$
$\displaystyle\quad-(1-m^{2})Cov_{\mathbf{z}}(q_{1,I_{1}},q_{1,J_{1}}\sigma_{I_{2},J_{1}})+2m(1-m)x\,Cov_{\mathbf{z}}(q_{1,I_{1}},\sigma_{I_{2},J_{1}})$
$\displaystyle\quad+(1-m)^{2}Cov_{\mathbf{z}}(q_{1,I_{1}},q_{1,J_{1}}\sigma_{I_{2},J_{2}})-(1-m)^{2}Cov_{\mathbf{z}}(q_{1,I_{1}},q_{1,J_{1}}\sigma_{I_{1},J_{2}}).$
(27)

The proof of this proposition relies on three technical lemmas proved in
Appendices C, D and E.

######  Lemma 3.

The conditional variance of the uniformly weighted frequency of $S_{1}$ after
viability selection within demes is given by

$\displaystyle Var_{\mathbf{z}}\left[\tilde{x}\right]$
$\displaystyle=\frac{1}{Nd}\Big{[}Cov_{\mathbf{z}}(q_{1,I_{1}},q_{1,J_{1}}\sigma_{I_{1},J_{1}})+Cov_{\mathbf{z}}(q_{1,I_{1}},q_{1,J_{1}}\sigma_{I_{2},J_{2}})-Cov_{\mathbf{z}}(q_{1,I_{1}},q_{1,J_{1}}\sigma_{I_{2},J_{1}})$
$\displaystyle\quad-
Cov_{\mathbf{z}}(q_{1,I_{1}},q_{1,J_{1}}\sigma_{I_{1},J_{2}})\Big{]}+o(d^{-1}).$
(28)

######  Lemma 4.

The conditional variance of the frequency of $S_{1}$ in the migrant pool is
given by

$Var_{\mathbf{z}}(\tilde{\tilde{x}})=\frac{1}{Nd}\Big{[}Cov_{\mathbf{z}}(q_{1,I_{1}}q_{1,J_{1}},\sigma_{I_{1},J_{1}})-2x\,Cov_{\mathbf{z}}(q_{1,I_{1}},\sigma_{I_{1},J_{1}})\Big{]}+o(d^{-1}).$
(29)

######  Lemma 5.

The conditional covariance between the uniformly weighted frequency of $S_{1}$
and the frequency of $S_{1}$ in the migrant pool is given by

$\displaystyle Cov_{\mathbf{z}}\left[\tilde{x},\tilde{\tilde{x}}\right]$
$\displaystyle=\frac{1}{Nd}\Big{[}Cov_{\mathbf{z}}(q_{1,I_{1}},q_{1,J_{1}}\sigma_{I_{1},J_{1}})+x\,Cov_{\mathbf{z}}(q_{1,I_{1}},\sigma_{I_{2},J_{1}})-x\,Cov_{\mathbf{z}}(q_{1,I_{1}},\sigma_{I_{1},J_{1}})$
$\displaystyle\quad-
Cov_{\mathbf{z}}(q_{1,I_{1}},q_{1,J_{1}}\sigma_{I_{2},J_{1}})\Big{]}+o(d^{-1}).$
(30)

###### Proof of Proposition 2.

First, we note that

$\displaystyle\sum_{j=0}^{N}j^{2}P_{ij}(\mathbf{z})$
$\displaystyle=(1-m)\sum_{j=0}^{N}j^{2}\binom{N}{j}\left(\tilde{x}_{i}\right)^{j}\left(1-\tilde{x}_{i}\right)^{N-j}+m\sum_{j=0}^{N}j^{2}\binom{N}{j}\tilde{\tilde{x}}^{j}\left(1-\tilde{\tilde{x}}\right)^{N-j}$
$\displaystyle=(1-m)N\tilde{x}_{i}\left((N-1)\tilde{x}_{i}+1\right)+mN\tilde{\tilde{x}}\left((N-1)\tilde{\tilde{x}}+1\right).$
(31)

Combining the above identity with Eqs. (17) and (24), we get

$\displaystyle Var_{\mathbf{z}}\left[X(1)-X(0)\Big{|}(a_{k,l})\right]$
$\displaystyle=Var_{\mathbf{z}}\left[\sum_{j=0}^{N}\frac{j}{N}Z_{j}(1)\Big{|}(a_{k,l})\right]$
$\displaystyle=\sum_{j=0}^{N}\frac{j^{2}}{N^{2}}Var_{\mathbf{z}}\left[Z_{j}(1)\Big{|}(a_{k,l})\right]+\sum_{j_{1}\not=j_{2}}\frac{j_{1}j_{2}}{N^{2}}Cov_{\mathbf{z}}\left(Z_{j_{1}}(1),Z_{j_{2}}(1)\Big{|}(a_{k,l})\right)$
$\displaystyle=\sum_{j=0}^{N}\frac{j^{2}}{dN^{2}}\sum_{i=0}^{N}z_{i}P_{ij}(\mathbf{z})\left(1-P_{ij}(\mathbf{z})\right)-\sum_{j_{1}\not=j_{2}}\frac{j_{1}j_{2}}{dN^{2}}\sum_{i=0}^{N}z_{i}P_{ij_{1}}(\mathbf{z})P_{ij_{2}}(\mathbf{z})$
$\displaystyle=\frac{1}{d}\sum_{i=0}^{N}z_{i}\left[\sum_{j=0}^{N}\frac{j^{2}}{N^{2}}P_{ij}(\mathbf{z})-\left(\sum_{j=0}^{N}\frac{j}{N}P_{ij}(\mathbf{z})\right)^{2}\right]$
$\displaystyle=\frac{1}{Nd}\sum_{i=0}^{N}z_{i}\big{[}(1-m)(Nm-1)\tilde{x}^{2}_{i}+m(N(1-m)-1)\tilde{\tilde{x}}^{2}$
$\displaystyle\quad\quad\quad\quad\quad\quad+(1-m)\tilde{x}_{i}+m\tilde{\tilde{x}}-2Nm(1-m)\tilde{x}_{i}\tilde{\tilde{x}}\big{]}$
$\displaystyle=\frac{1}{Nd}\Big{[}(1-m)(Nm-1)\sum_{i=0}^{N}z_{i}\tilde{x}^{2}_{i}+m(N(1-m)-1)\tilde{\tilde{x}}^{2}$
$\displaystyle\quad\quad\quad+(1-m)\tilde{x}+m\tilde{\tilde{x}}-2Nm(1-m)\tilde{x}\tilde{\tilde{x}}\Big{]}.$
(32)

Using this equation with Eq. (3) and the facts that

$\displaystyle
E_{\mathbf{z}}\left(\tilde{x}_{i}^{n}\right)=x_{i}^{n}+o\left(1\right),$ (33a)
$\displaystyle
E_{\mathbf{z}}\left(\tilde{x}\tilde{\tilde{x}}\right)=x^{2}+o\left(1\right),$
(33b)

for $n\geq 1$, we obtain

$\displaystyle Var_{\mathbf{z}}\left[X(1)-X(0)\right]$
$\displaystyle=Var_{\mathbf{z}}\left(E_{\mathbf{z}}\left[X(1)-X(0)\Big{|}(a_{k,l})\right]\right)+E_{\mathbf{z}}\left(Var_{\mathbf{z}}\left[X(1)-X(0)\Big{|}(a_{k,l})\right]\right)$
$\displaystyle=Var_{\mathbf{z}}\Big{(}(1-m)\tilde{x}+m\tilde{\tilde{x}}-x\Big{)}+\frac{1}{Nd}\Big{[}(1-m)(Nm-1)\sum_{i=0}^{N}z_{i}x^{2}_{i}$
$\displaystyle\quad+m(N(1-m)-1)x^{2}+(1-m)x+mx-2Nm(1-m)x^{2}\Big{]}+o\left(d^{-1}\right)$
$\displaystyle=(1-m)^{2}Var_{\mathbf{z}}\left[\tilde{x}\right]+m^{2}Var_{\mathbf{z}}\left[\tilde{\tilde{x}}\right]+2m(1-m)Cov_{\mathbf{z}}\left[\tilde{x},\tilde{\tilde{x}}\right]+\frac{1}{Nd}$
(34)
$\displaystyle\quad\times\Big{[}(1-m)(Nm-1)\sum_{i=0}^{N}z_{i}x^{2}_{i}-m(N(1-m)+1)x^{2}+x\Big{]}+o\left(d^{-1}\right).$

Using the identities in Eq. (19) and Lemmas 3-5 leads to the statement of
Proposition 2. ∎

The next result, proved in Appendix F, implies that the fourth conditional
moment of the change in the frequency of $S_{1}$ is negligible compared to the
conditional moments of order one and two when the number of demes is large.

######  Proposition 3.

For the fourth conditional moment of the change in the frequency of $S_{1}$ in
the whole population from the beginning of generation $0$ to the beginning of
generation $1$, we have

$E_{\mathbf{z}}\left[\Big{(}X(1)-X(0)\Big{)}^{4}\right]=o(d^{-1}).$ (35)

Let us now focus on the second timescale for the deviations of the deme type
frequencies from their equilibrium values in a neutral population subdivided
into an infinite number of demes. We first note that, as a consequence of Eq.
(14), we have

$E_{\mathbf{z}}\Big{[}P_{ij}(\mathbf{z})\Big{]}=P^{*}_{ij}(x)+o(1),$ (36)

as $S_{2}$ approaches infinity, where

$P^{*}_{ij}(x)=(1-m)\binom{N}{j}x_{i}^{j}\left(1-x_{i}\right)^{N-j}+m\binom{N}{j}x^{j}\left(1-x\right)^{N-j},$
(37)

for $i,j=0,1,\ldots,N$, with
$x=\sum_{j=0}^{N}x_{j}z_{j}=\sum_{j=0}^{N}jz_{j}/N$. Here, $P^{*}_{ij}(x)$
represents the transition probability for a deme of type $i$ at a the
beginning of generation $0$ to become a deme of type $j$ at the beginning of
generation $1$ in an infinite neutral population subdivided into an infinite
number of demes each of size $N$, where the frequency of $S_{1}$ is constant
and given by $x$.

Let $\mathbf{v}(x)=(v_{0}(x),v_{1}(x),\ldots,v_{N}(x))$ be the frequency
vector that is the solution to the linear system of equations

$v_{j}(x)=\sum_{i=0}^{N}v_{i}(x)P_{ij}^{*}(x)$ (38)

for $j=0,1,\ldots,N$. In other words, $\mathbf{v}(x)$ is the stationary
probability distribution of the process $\left\\{\mathbf{Z}(t)\right\\}_{t\geq
0}$ for the deme type frequencies in the absence of selection, when the
population is subdivided into an infinite number of demes of the same finite
size $N$, and the total frequency of type $S_{1}$ is constant and equal to
$x$.

The focus of our analysis will now shift to the process
$\left\\{\mathbf{Y}(t)=\left(Y_{0}(t),Y_{1}(t),\ldots,Y_{N}(t)\right)\right\\}_{t\geq
0}$, where

$Y_{j}(t)=Z_{j}(t)-v_{j}(X(t))$ (39)

for $j=0,1,\ldots,N$ define the deviations of the deme type frequencies from
their equilibrium values as the number of demes $S_{2}$ goes to infinity.

######  Proposition 4.

The first conditional moment of the changes in the deviations of the deme type
frequencies from their equilibrium values as $d\rightarrow\infty$ from the
beginning of generation $0$ to the beginning of generation $1$ is given by

$E_{\mathbf{z}}\left[Y_{i}(1)-Y_{i}(0)\right]=c_{i}(x,\mathbf{y})+o\left(1\right),$
(40)

where

$c_{i}(x,\mathbf{y})=\sum_{j=0}^{N}y_{j}P_{ji}^{*}(x)-y_{i}$ (41)

with $x=\sum_{j=0}^{N}x_{j}z_{j}$ and $y_{i}=z_{i}-v_{i}(x)$, for
$i=0,1,\ldots,N$.

###### Proof.

It has been shown in Lessard [32] that

$v_{i}(X(1))=v_{i}(x)+\sum_{j=1}^{N}r_{j}\left(X(1)-X(0)\right)^{j},$ (42)

for $i=0,1,\ldots,N$, where $r_{j}$ is a constant that depends only on $N$,
$j$, $m$ and $x$, for $j=1,\ldots,N$. Furthermore, we have

$\left|E_{\mathbf{z}}\left[\left(X(1)-X(0)\right)^{k}\right]\right|\leq
E_{\mathbf{z}}\left[\left|X(1)-X(0)\right|^{k}\right]\leq
E_{\mathbf{z}}\left[\left(X(1)-X(0)\right)^{2}\right]$ (43)

for any integer $k\geq 2$. Hence, owing to Propositions 1 and 2, we have

$E_{\mathbf{z}}\left[v_{i}(X(1))\right]=v_{i}(x)+\sum_{j=1}^{N}r_{j}E_{\mathbf{z}}\left[\left(X(1)-X(0)\right)^{j}\right]=v_{i}(x)+o(1).$
(44)

On the other hand, we have

$\displaystyle E_{\mathbf{z}}\left[Z_{i}(1)\right]$
$\displaystyle=\sum_{j=0}^{N}z_{j}E_{\mathbf{z}}\left[P_{ji}(\mathbf{z})\right]$
$\displaystyle=\sum_{j=0}^{N}z_{j}P_{ji}^{*}(x)+o(1)$
$\displaystyle=\sum_{j=0}^{N}(y_{j}+v_{j}(x))P_{ji}^{*}(x)+o(1)$
$\displaystyle=\sum_{j=0}^{N}y_{j}P_{ji}^{*}(x)+v_{i}(x)+o(1),$ (45)

where $y_{j}=z_{j}-v_{j}(x)$ for $j=0,1,\ldots,N$. By combining (44) and (3),
the first conditional moment of the change in $Y_{i}$ from the beginning of
generation 0 to the beginning of generation 1 can be expressed as

$E_{\mathbf{z}}\left[Y_{i}(1)-Y_{i}(0)\right]=E_{\mathbf{z}}\left[Z_{i}(1)\right]-E_{\mathbf{z}}\left[v_{i}(X(1))\right]-y_{i}=c_{i}(x,\mathbf{y})+o\left(1\right),$
(46)

with $c_{i}(x,\mathbf{y})$ as given in the statement of Proposition 4, for
$i=0,1,\ldots,N$. ∎

On the other hand, it can be shown (see Appendix G for a proof) that the
conditional variance of the change in $Y_{i}$ is negligible compared to the
first conditional moment when the number of demes is large.

######  Proposition 5.

For the conditional variance of the change in $Y_{i}$ from the beginning of
generation 0 to the beginning of generation 1, we have

$Var_{\mathbf{z}}\left[Y_{i}(1)-Y_{i}(0)\right]=o\left(1\right)$ (47)

for $i=0,1,\ldots,N$.

Finally, we state and prove a convergence result that holds for an infinite
neutral population.

######  Proposition 6.

For the function
$\mathbf{c}(x,\mathbf{y})=(c_{0}(x,\mathbf{y}),c_{1}(x,\mathbf{y}),\ldots,c_{N}(x,\mathbf{y}))$
whose components are defined in Eq. (41), we have
$\mathbf{c}(x,\mathbf{0})=\mathbf{0}$ and the zero solution of the recurrence
equation

$\mathbf{Y}(k+1,x,\mathbf{y})-\mathbf{Y}(k,x,\mathbf{y})=\mathbf{c}(x,\mathbf{Y}(k,x,\mathbf{y})),$
(48)

for $k\geq 0$ with the initial condition
$\mathbf{Y}(0,x,\mathbf{y})=\mathbf{y}$ is globally asymptotically stable.

###### Proof.

Note that the recurrence equation Eq. (48) can be rewritten as

$\mathbf{Z}(k+1,x,\mathbf{y})=\mathbf{Z}(k,x,\mathbf{y})P^{*}(x),$ (49)

where $\mathbf{Z}(k,x,\mathbf{y})=\mathbf{Y}(k,x,\mathbf{y})+\mathbf{v}(x)$
and $P^{*}(x)=(P_{i,j}^{*}(x))$ is as defined in Eq. (37). Moreover, owing to
Eq. (38), the frequency vector $\mathbf{v}(x)$ represents the stationary
distribution associated with $P^{*}(x)$ which is the transition matrix of an
irreducible Markov chain on a finite state space. Then, the ergodic theorem
ensures that

$\lim_{k\rightarrow\infty}\mathbf{Z}(k,x,\mathbf{y})=\mathbf{v}(x)$ (50)

for any $\mathbf{y}$, which corresponds to the statement of Proposition 6. ∎

Owing to Propositions (1)-(6), which are actually valid uniformly with respect
to the population state, and assuming the existence and uniqueness of a
diffusion process with infinitesimal mean and variance $M(\mathbf{v}(x))$ and
$Q(\mathbf{v}(x))$, respectively, whose closure determines a strongly
continuous semigroup (see Ethier [10]), Theorem 3.3 in Ethier and Nagylaki
[11] ensures the following convergence theorem.

######  Theorem 1.

Let $X(\lfloor Nd\tau\rfloor)$ be the frequency of type $S_{1}$ at the
beginning of generation $\lfloor Nd\tau\rfloor$ in the model of Section 2.
Here, $\lfloor Nd\tau\rfloor$ denotes the integer part of $Nd\tau$ for any
real $\tau\geq 0$. Then, the process $\\{X(\lfloor Nd\tau\rfloor)\\}_{\tau\geq
0}$ converges in distribution to a diffusion process in $[0,1]$ whose
infinitesimal generator is

$\mathcal{L}=\frac{1}{2}Q(\mathbf{v}(x))\frac{d^{2}}{dx^{2}}+M(\mathbf{v}(x))\frac{d}{dx},$
(51)

where the functions $M$ and $Q$ are given in Propositions 1 and 2,
respectively, while the frequency vector $\mathbf{v}(x)$ is defined by Eqs.
(37) and (38).

##  4 Infinitesimal mean and variance

The infinitesimal mean and variance $M(\mathbf{v}(x))$ and $Q(\mathbf{v}(x))$
in the statement of Theorem 1 can be expressed in terms of identity-by-descent
measures in a neutral population.

Let
$I_{1}^{(1)},\ldots,I_{n_{1}}^{(1)},\ldots,I_{1}^{(k)},\ldots,I_{n_{k}}^{(k)}$
be an ordered sample of $n=n_{1}+\cdots+n_{k}$ offspring chosen at random in
the same deme at the beginning of a given generation in a neutral population
of an infinite number of demes that is in the stationary state. Denote by

$\displaystyle
f_{I_{1}^{(1)}\ldots\,I_{n_{1}}^{(1)}|\,\ldots\,|I_{1}^{(k)}\ldots\,I_{n_{k}}^{(k)}}=f_{n_{1}n_{2}\ldots
n_{k}}$ (52)

the probability that $I_{1}^{(i)},\ldots,I_{n_{i}}^{(i)}$ are identical by
descent for $i=1,\ldots,k$, but $I_{1}^{(i)},\ldots,I_{n_{i}}^{(i)}$ and
$I_{1}^{(j)},\ldots,I_{n_{j}}^{(j)}$ are not identical by descent for
$i,j=1,\ldots,k$ with $i\neq j$. The exact expressions of these identity-by-
descent measures for up to five offspring are given in Appendix H.

The above identity-by-descent measures lead to exact expressions for
$M(\mathbf{v}(x))$ and $Q(\mathbf{v}(x))$ deduced in Appendix I.
Simplifications are possible in the case where the deme size is large and the
probability of extinction is small.

######  Proposition 7.

If $N$ is large enough and $m$ small enough, the infinitesimal mean and
variance of the diffusion process in Theorem 1 can be approximated as

$\displaystyle M(\mathbf{v}(x))\approx$
$\displaystyle\;x(1-x)\Big{[}(\mu_{1,1}-\mu_{2,1})\psi(x)+(\mu_{1,2}-\mu_{2,2})\psi(1-x)-\sigma_{1,1}^{2}g(x)-\sigma_{1,2}^{2}h(x)$
$\displaystyle\quad\quad\quad\quad+\sigma_{2,1}^{2}h(1-x)+\sigma_{2,2}^{2}g(1-x)\Big{]},$
(53a) $\displaystyle Q(\mathbf{v}(x))\approx$
$\displaystyle\;x(1-x)\Big{[}f_{11}+\nu
f_{2}+x(1-x)\psi^{2}(x)\left(\sigma^{2}_{1,1}+\sigma^{2}_{2,1}\right)$
$\displaystyle\quad\quad\quad\quad+x(1-x)\psi^{2}(1-x)\left(\sigma^{2}_{1,2}+\sigma^{2}_{2,2}\right)\Big{]},$
(53b)

where $\nu=Nm$ is a deme-scaled extinction rate and

$\displaystyle\psi(x)$ $\displaystyle=f_{21}+xf_{111},$ (54a) $\displaystyle
g(x)$
$\displaystyle=f_{41}+(4f_{311}+3f_{221})x+6f_{2111}x^{2}+f_{11111}x^{3},$
(54b) $\displaystyle h(x)$
$\displaystyle=f_{32}+\left(3f_{221}+f_{2111}\right)(1-x)+xf_{311}+2x(1-x)f_{2111}+x(1-x)^{2}f_{11111}.$
(54c)

Moreover, we have the following ancillary result for the identity-by-descent
measures, proved in Appendix J.

######  Proposition 8.

For a deme size $N$ large enough and a probability of extinction $m$ small
enough, we have the approximation

$f_{n_{1}n_{2}\ldots
n_{k}}\approx\binom{n_{1}+\cdots+n_{k}-k}{n_{1}-1,\ldots,n_{k}-1}\frac{\prod_{i=1}^{k}\prod_{j=2}^{n_{i}}\binom{j}{2}}{\prod_{j=k+1}^{n_{1}+\cdots+n_{k}}\left(\binom{j}{2}+\nu\right)}\times\frac{\nu}{\binom{k}{2}+\nu}$
(55)

for $n_{1},\ldots,n_{k}\geq 1$ for $k\geq 1$, where $\nu=Nm$.

From now on, the conditions for Propositions 7 and 8 will be assumed to hold.

##  5 Fixation probability

In the absence of mutation, the frequency of $S_{1}$ in the model of Section 2
is a Markov chain that fixes into either $x=1$, all individuals are of type
$S_{1}$, or $x=0$, all individuals are of type $S_{2}$. On the other hand, the
probability of absorption in state $x=1$ at or before time $t$ starting in
state $x=p$ at time $0$ in the diffusion process of Theorem 1, denoted by
$P_{1}(p,t)$, is known to satisfy the backward Kolmogorv equation

$-\frac{\partial P_{1}(p,t)}{\partial t}=M(\mathbf{v}(x))\frac{\partial
P_{1}(p,t)}{\partial
p}+\frac{Q(\mathbf{v}(x))}{2}\frac{\partial^{2}P_{1}(p,t)}{\partial p^{2}},$
(56)

with the boundary conditions $P_{1}(0,t)=0$ and $P_{1}(1,t)=1$. See Karlin and
Taylor [19] or Ewens [12] for more details. Moreover,

$\displaystyle u_{1}(p)=\lim_{t\rightarrow\infty}P_{1}(1,t),$ (57)

which is the fixation probability of $S_{1}$ given that $p$ is the initial
frequency of $S_{1}$, satisfies

$0=M(\mathbf{v}(x))\frac{du_{1}(p)}{dp}+\frac{Q(\mathbf{v}(x))}{2}\frac{d^{2}u_{1}(p)}{dp^{2}},$
(58)

with the boundary conditions $u_{1}(0)=0$ and $u_{1}(1)=1$. This equation has
the exact solution

$u_{1}(p)=\frac{\int_{0}^{p}\Psi(y)dy}{\int_{0}^{1}\Psi(y)dy},$ (59)

where

$\displaystyle\Psi(y)=\exp\left\\{-2\int_{0}^{y}\frac{M(\mathbf{v}(x))}{Q(\mathbf{v}(x))}dx\right\\}.$
(60)

Similarly,

$u_{2}(p)=1-u_{1}(1-p)=\frac{\int_{1-p}^{1}\Psi(y)dy}{\int_{0}^{1}\Psi(y)dy}$
(61)

is the fixation probability of $S_{2}$ given that the initial frequency of
$S_{2}$ is $p$.

Let $F_{1}$ (respectively, $F_{2}$) denote the fixation probability of $S_{1}$
(respectively, $S_{2}$) introduced as a single mutant in a population that was
initially fixed for $S_{2}$ (respectively, $S_{1}$). Assuming a sufficiently
large number of demes and using the diffusion approximation established in the
previous section, we obtain

$\displaystyle F_{1}\approx
u_{1}\left(\frac{1}{Nd}\right)=\frac{\int_{0}^{(Nd)^{-1}}\Psi(y)dy}{\int_{0}^{1}\Psi(y)dy}\approx\frac{\Psi(0)}{Nd\int_{0}^{1}\Psi(y)dy},$
(62a) $\displaystyle F_{2}\approx
u_{2}\left(\frac{1}{Nd}\right)=\frac{\int_{1-(Nd)^{-1}}^{1}\Psi(y)dy}{\int_{0}^{1}\Psi(y)dy}\approx\frac{\Psi(1)}{Nd\int_{0}^{1}\Psi(y)dy},$
(62b)

since

$\displaystyle\int_{0}^{(Nd)^{-1}}\Psi(y)dy\approx\frac{\Psi(0)}{Nd},$ (63a)
$\displaystyle\int_{1-(Nd)^{-1}}^{1}\Psi(y)dy\approx\frac{\Psi(1)}{Nd}.$ (63b)

In the absence of selection, it can be noted that $\Psi(y)=1$ for any
$y\in[0,1]$. Then, the above fixation probabilities are equal and given by the
initial frequency, that is,

$F_{1}=F_{2}=\frac{1}{Nd}.$ (64)

In the presence of selection, the evolution of $S_{i}$ is said to be favoured
if its fixation probability when introduced as a single mutant is greater than
what it would be in the absence of selection, that is, $F_{i}>(Nd)^{-1}$, for
$i=1,2$ (Nowak et al. [43]). If the inequality is reversed, that is,
$F_{i}<(Nd)^{-1}$, then the evolution of $S_{i}$ is said to be disfavoured,
for $i=1,2$. It should be noted that the evolution of both types can be
favoured or disfavoured by selection. Thus, a third criterion is introduced to
determine which type has a higher probability of invading a population
consisting entirely of the other type. Selection is said to favour the
evolution of $S_{1}$ more than the evolution of $S_{2}$ if $F_{1}>F_{2}$.

Using Eq. (62) with the fact that $\Psi(0)=1$ and introducing the function

$\Phi(y)=\frac{\Psi(y)}{\Psi(1)}=\exp\left(2\int_{y}^{1}\frac{M(\mathbf{v}(x))}{Q(\mathbf{v}(x))}dx\right),$
(65)

where $\Psi$ is defined in Eq. (60), we can conclude the following.

######  Proposition 9.

Selection favours the evolution of $S_{1}$ in the sense that $F_{1}>(Nd)^{-1}$
as long as

$\int_{0}^{1}\Psi(y)dy<1.$ (66)

On the other hand, the condition for selection to disfavour the evolution of
$S_{2}$ in the sense that $F_{2}<(Nd)^{-1}$ is

$\int_{0}^{1}\Phi(y)dy>1.$ (67)

Finally, the condition for selection to favour the evolution of $S_{1}$ more
than the evolution of $S_{2}$ in the sense that $F_{1}>F_{2}$ is

$\Psi(1)<1.$ (68)

If the above three conditions are satisfied, then $F_{1}>(Nd)^{-1}>F_{2}$ and
we say that selection fully favours the evolution of $S_{1}$.

##  6 Small population-scaled means and variances

Let us assume that the population-scaled means $\mu_{1,1}$, $\mu_{1,2}$,
$\mu_{2,1}$ and $\mu_{2,2}$ as well as population-scaled variances
$\sigma_{1,1}^{2}$, $\sigma_{1,2}^{2}$, $\sigma_{2,1}^{2}$ and
$\sigma_{2,2}^{2}$ are all of sufficiently small order. Then, using
Proposition 7 and the fact that $\nu f_{2}=\nu/(1+\nu)=f_{11}$, we get the
approximation

$\displaystyle\frac{M(\mathbf{v}(x))}{Q(\mathbf{v}(x))}\approx\frac{M(\mathbf{v}(x))}{2f_{11}x(1-x)}$
$\displaystyle=\frac{\psi(x)}{2f_{11}}(\mu_{1,1}-\mu_{2,1})+\frac{\psi(1-x)}{2f_{11}}(\mu_{1,2}-\mu_{2,2})$
$\displaystyle\quad-\frac{g(x)}{2f_{11}}\sigma_{1,1}^{2}-\frac{h(x)}{2f_{11}}\sigma_{1,2}^{2}+\frac{h(1-x)}{2f_{11}}\sigma_{2,1}^{2}+\frac{g(1-x)}{2f_{11}}\sigma_{2,2}^{2}.$
(69)

By using this approximation, along with $e^{u}\approx 1+u$ when $u$ is
sufficiently small, we obtain

$\displaystyle\Psi(y)$
$\displaystyle=\exp\left(-2\int_{0}^{y}\frac{M(\mathbf{v}(x))}{Q(\mathbf{v}(x))}\right)\approx
1-\int_{0}^{y}\frac{M(\mathbf{v}(x))}{f_{11}x(1-x)}dx,$ (70a)
$\displaystyle\Phi(y)$
$\displaystyle=\exp\left(2\int_{y}^{1}\frac{M(\mathbf{v}(x))}{Q(\mathbf{v}(x))}\right)\approx
1+\int_{y}^{1}\frac{M(\mathbf{v}(x))}{f_{11}x(1-x)}dx.$ (70b)

Then, applying the conditions in Proposition 9 and using the approximations in
Proposition 8 for the identity-by-descent measures yield our first result.

##### Result 1.

Assuming that the population-scaled means and variances of the payoffs are of
the same small order, selection will favour the evolution of $S_{1}$ if

$\displaystyle\frac{4+\nu}{3(3+\nu)}(\mu_{1,1}-\mu_{2,1})+\frac{5+2\nu}{3(3+\nu)}(\mu_{1,2}-\mu_{2,2})-\frac{288+127\nu+20\nu^{2}+\nu^{3}}{10(3+\nu)(6+\nu)(10+\nu)}\sigma^{2}_{1,1}$
$\displaystyle-\frac{188+117\nu+20\nu^{2}+\nu^{3}}{10(3+\nu)(6+\nu)(10+\nu)}\sigma^{2}_{1,2}+\frac{243+127\nu+20\nu^{2}+\nu^{3}}{15(3+\nu)(6+\nu)(10+\nu)}\sigma^{2}_{2,1}$
$\displaystyle+\frac{231+174\nu+35\nu^{2}+2\nu^{3}}{5(3+\nu)(6+\nu)(10+\nu)}\sigma^{2}_{2,2}>0,$
(71)

disfavour the evolution of $S_{2}$ if

$\displaystyle\frac{5+2\nu}{3(3+\nu)}(\mu_{1,1}-\mu_{2,1})+\frac{4+\nu}{3(3+\nu)}(\mu_{1,2}-\mu_{2,2})-\frac{231+174\nu+35\nu^{2}+2\nu^{3}}{5(3+\nu)(6+\nu)(10+\nu)}\sigma^{2}_{1,1}$
$\displaystyle-\frac{243+127\nu+20\nu^{2}+\nu^{3}}{15(3+\nu)(6+\nu)(10+\nu)}\sigma^{2}_{1,2}+\frac{188+117\nu+20\nu^{2}+\nu^{3}}{10(3+\nu)(6+\nu)(10+\nu)}\sigma^{2}_{2,1}$
$\displaystyle+\frac{288+127\nu+20\nu^{2}+\nu^{3}}{10(3+\nu)(6+\nu)(10+\nu)}\sigma^{2}_{2,2}>0,$
(72)

and favour the evolution of $S_{1}$ more than the evolution of $S_{2}$ if

$\mu_{1,1}+\mu_{1,2}-\mu_{2,1}-\mu_{2,2}+\frac{\nu+5}{2\left(\nu+6\right)}(\sigma^{2}_{2,2}-\sigma^{2}_{1,1})+\frac{\nu+7}{6\left(\nu+6\right)}(\sigma^{2}_{2,1}-\sigma^{2}_{1,2})>0.$
(73)

It should be noted that the coefficients of $\sigma^{2}_{1,1}$ and
$\sigma^{2}_{1,2}$ in the above conditions are all negative, while the
coefficients of $\sigma^{2}_{2,1}$ and $\sigma^{2}_{2,2}$ are all positive. As
a consequence, reducing the population-scaled variance of any payoff to
$S_{1}$ or increasing the population-scaled variance of any payoff to $S_{2}$
makes all the above conditions for the evolution of $S_{1}$ easier to be
satisfied. This is even the only way to make the conditions possible to hold
when the population-scaled means of the payoffs verify $\mu_{1,1}<\mu_{2,1}$
and $\mu_{1,2}<\mu_{2,2}$ as in a deterministic Prisoner’s dilemma.

As for the effect of the deme-scaled extinction rate, we note that the
coefficient of $\sigma^{2}_{2,2}-\sigma^{2}_{1,1}$ in the condition for
selection to favour the evolution of $S_{1}$ more than the evolution of
$S_{2}$ is an increasing positive function with respect to $\nu$, while the
coefficient of $\sigma^{2}_{2,1}-\sigma^{2}_{1,2}$ is a decreasing positive
function. This means that the effect of variability in payoffs to interacting
offspring of the same type on this condition decreases as $\nu$ increases,
while the effect of variability in payoffs to interacting offspring of
different types decreases. This is surprising since, as $\nu$ increases,
interacting offspring are less likely to be of the same type.

Finally, as $\nu\rightarrow\infty$, the above conditions for
$F_{1}>(Nd)^{-1}$, $F_{2}<(Nd)^{-1}$ and $F_{1}>F_{2}$ become

$\displaystyle\frac{1}{3}(\mu_{1,1}-\mu_{2,1})+\frac{2}{3}(\mu_{1,2}-\mu_{2,2})-\frac{1}{10}\sigma^{2}_{1,1}-\frac{1}{10}\sigma^{2}_{1,2}+\frac{1}{15}\sigma^{2}_{2,1}+\frac{2}{5}\sigma^{2}_{2,2}>0,$
(74a)
$\displaystyle\frac{2}{3}(\mu_{1,1}-\mu_{2,1})+\frac{1}{3}(\mu_{1,2}-\mu_{2,2})-\frac{2}{5}\sigma^{2}_{1,1}-\frac{1}{15}\sigma^{2}_{1,2}+\frac{1}{10}\sigma^{2}_{2,1}+\frac{1}{10}\sigma^{2}_{2,2}>0,$
(74b)
$\displaystyle\mu_{1,1}+\mu_{1,2}-\mu_{2,1}-\mu_{2,2}+\frac{1}{2}(\sigma^{2}_{2,2}-\sigma^{2}_{1,1})+\frac{1}{6}(\sigma^{2}_{2,1}-\sigma^{2}_{1,2})>0,$
(74c)

respectively. Note that similar conditions were derived for well-mixed
populations (see Eqs. ($38$, $39$, $40$) in Kroumi et al. [27] for a Moran
model, and Eqs. ($95$, $100$, $104$) in Li and Lessard [34] for a Wright-
Fisher model). Note also that, if the population-scaled variances are all
equal to $\sigma^{2}$, then the condition for selection to favour the
evolution of $S_{1}$ more than the evolution of $S_{2}$ ($F_{1}>F_{2}$) does
not depend on $\sigma^{2}$, which is not surprising by symmetry, but it
becomes easier for selection to favour of the evolution of $S_{1}$
($F_{1}>(Nd)^{-1}$) as well as the evolution of $S_{2}$ ($F_{2}>(Nd)^{-1}$) as
$\sigma^{2}$ increases.

##  7 Case 1: $\sigma^{2}_{1,1}=\sigma^{2}_{1,2}=\sigma^{2}_{2,1}=0$ and
$\sigma^{2}_{2,2}=\sigma^{2}>0$

Here, we suppose that the population-scaled variances of the payoffs different
from $\sigma^{2}_{2,2}=\sigma^{2}>0$ are insignificant. In this case, we have

$\frac{M(\mathbf{v}(x))}{Q(\mathbf{v}(x))}=\frac{(\mu_{1,1}-\mu_{2,1})\psi(x)+(\mu_{1,2}-\mu_{2,2})\psi(1-x)+\sigma^{2}g(1-x)}{2f_{11}+x(1-x)\psi^{2}(1-x)\sigma^{2}},$
(75)

whose first derivative with respect to $\sigma^{2}$ is given by

$\frac{2f_{11}g(1-x)+x(1-x)\psi^{2}(1-x)\Big{[}(\mu_{2,1}-\mu_{1,1})\psi(x)+(\mu_{2,2}-\mu_{1,2})\psi(1-x)\Big{]}}{\left(2f_{11}+x(1-x)\sigma^{2}\psi^{2}(1-x)\right)^{2}}>0,$
(76)

for $x\in[0,1]$. Using Leibniz integral rule, it can be shown that for every
$y\in(0,1]$, the functions

$\displaystyle\sigma^{2}\longmapsto\int_{0}^{y}\frac{M(\mathbf{v}(x))}{Q(\mathbf{v}(x))}dx=-\frac{\ln(\Psi(y))}{2}$
(77)

and

$\displaystyle\sigma^{2}\longmapsto\int_{y}^{1}\frac{M(\mathbf{v}(x))}{Q(\mathbf{v}(x))}dx=\frac{\ln(\Phi(y))}{2}$
(78)

are strictly increasing on $[0,\infty)$, so that the functions
$\sigma^{2}\longmapsto\Psi(y)$ and $\sigma^{2}\longmapsto-\Phi(y)$ are
strictly decreasing on $[0,\infty)$. Hence, an increase in $\sigma^{2}$ will
increase the fixation probability $F_{1}$ and decrease the fixation
probability $F_{2}$. In addition, increasing $\sigma^{2}$ will make it
possible for selection to favour the evolution of $S_{1}$.

Using Fatou’s lemma, we get

$\liminf_{\sigma^{2}\rightarrow\infty}\int_{0}^{y}\frac{M(\mathbf{v}(x))}{Q(\mathbf{v}(x))}dx\geq\int_{0}^{y}\lim_{\sigma^{2}\rightarrow\infty}\frac{M(\mathbf{v}(x))}{Q(\mathbf{v}(x))}dx=\int_{0}^{y}\frac{g(1-x)}{x(1-x)\psi^{2}(1-x)}dx=\infty$
(79)

for every $y\in(0,1]$, since

$\frac{g(1-x)}{x(1-x)\psi^{2}(1-x)}\sim\frac{g(1)}{x\psi^{2}(1)}$ (80)

as $x\rightarrow 0^{+}$. Therefore, we have

$\lim_{\sigma^{2}\rightarrow\infty}\int_{0}^{y}\frac{M(\mathbf{v}(x))}{Q(\mathbf{v}(x))}dx=\infty,$
(81)

from which

$\lim_{\sigma^{2}\rightarrow\infty}\Psi(y)=\lim_{\sigma^{2}\rightarrow\infty}\exp\left(-2\int_{0}^{y}\frac{M(\mathbf{v}(x))}{Q(\mathbf{v}(x))}dx\right)=0,$
(82)

for every $y\in(0,1]$. On the other hand, when $\sigma^{2}=0$, we have
$M(\mathbf{v}(x))/Q(\mathbf{v}(x))<0$ for $x\in[0,1]$, which implies that

$\Psi(y)|_{\sigma^{2}=0}>1\;\mbox{ for every }y\in(0,1].$ (83)

Combining (82) and (83) for $y=1$, and using the fact that the function
$\sigma^{2}\longmapsto\Psi(1)$ is continuous and strictly decreasing on
$[0,\infty)$, we deduce that the equation

$\Psi(1)=1,$ (84)

with respect to $\sigma^{2}>0$, has a unique solution
$\sigma^{2}_{**}(\nu)>0$, so that $F_{1}>F_{2}$ for
$\sigma^{2}>\sigma^{2}_{**}(\nu)$ and $F_{1}<F_{2}$ for
$\sigma^{2}<\sigma^{2}_{**}(\nu)$.

Similarly, Fatou’s lemma yields

$\liminf_{\sigma^{2}\rightarrow\infty}\int_{y}^{1}\frac{M(\mathbf{v}(x))}{Q(\mathbf{v}(x))}dx\geq\int_{y}^{1}\lim_{\sigma^{2}\rightarrow\infty}\frac{M(\mathbf{v}(x))}{Q(\mathbf{v}(x))}dx=\int_{y}^{1}\frac{g(1-x)}{x(1-x)\psi^{2}(1-x)}dx=\infty$
(85)

for every $y\in[0,1)$, since

$\frac{g(1-x)}{x(1-x)\psi^{2}(1-x)}\sim\frac{g(0)}{(1-x)\psi^{2}(0)}$ (86)

as $x\rightarrow 1^{-}$, which implies

$\lim_{\sigma^{2}\rightarrow\infty}\int_{y}^{1}\frac{M(\mathbf{v}(x))}{Q(\mathbf{v}(x))}dx=\infty$
(87)

and

$\lim_{\sigma^{2}\rightarrow\infty}\Phi(y)=\lim_{\sigma^{2}\rightarrow\infty}\exp\left(2\int_{y}^{1}\frac{M(\mathbf{v}(x))}{Q(\mathbf{v}(x))}dx\right)=\infty.$
(88)

Using again Fatou’s lemma, we have

$\liminf_{\sigma^{2}\rightarrow\infty}\int_{0}^{1}\Phi(y)dy\geq\int_{0}^{1}\lim_{\sigma^{2}\rightarrow\infty}\Phi(x)dx=\infty,$
(89)

from which

$\lim_{\sigma^{2}\rightarrow\infty}\int_{0}^{1}\Phi(y)dy=\infty.$ (90)

On the other hand, the fact that $M(\mathbf{v}(x))/Q(\mathbf{v}(x))<0$ for
$x\in[0,1]$ when $\sigma^{2}=0$ implies that

$\int_{0}^{1}\Phi(y)dy|_{\sigma^{2}=0}<1.$ (91)

Combining (90) and (91), and using the fact that the function
$\sigma^{2}\longmapsto\int_{0}^{1}\Phi(y)dy$ is continuous and strictly
increasing on $[0,\infty)$, the equation

$\int_{0}^{1}\Phi(y)dy=1$ (92)

with respect to $\sigma^{2}>0$, has a unique solution
$\sigma^{2}_{***}(\nu)>0$, so that $F_{2}>(Nd)^{-1}$ for
$\sigma^{2}<\sigma^{2}_{***}(\nu)$ and $F_{2}<(Nd)^{-1}$ for
$\sigma^{2}>\sigma^{2}_{***}(\nu)$.

Note that, owing to Eq. (75), there exists a $\sigma^{2}_{0}>0$ such that if
$\sigma^{2}>\sigma^{2}_{0}$, then $M(\mathbf{v}(x))>0$ for $x\in[0,1]$, from
which $|\Psi(y)|<1$ for every $y\in(0,1]$. Then, adding this conclusion to Eq.
(82), the Lebesgue dominated convergence theorem ensures that

$\lim_{\sigma^{2}\rightarrow\infty}\int_{0}^{1}\Psi(y)dy=0.$ (93)

On the other hand, a consequence of Eq. (83) is that

$\int_{0}^{1}\Psi(y)dy|_{\sigma^{2}=0}>1.$ (94)

Combining (93) and (94), and using the fact that the function
$\sigma^{2}\longmapsto\int_{0}^{1}\Psi(y)dy$ is continuous and strictly
decreasing on $[0,\infty)$, the equation

$\int_{0}^{1}\Psi(y)dy=1$ (95)

with respect to $\sigma^{2}>0$, has a unique solution $\sigma^{2}_{*}(\nu)>0$,
so that $F_{1}<(Nd)^{-1}$ for $\sigma^{2}<\sigma^{2}_{*}(\nu)$ and
$F_{1}>(Nd)^{-1}$ for $\sigma^{2}>\sigma^{2}_{*}(\nu)$.

Next, we study the relationships between the threshold values
$\sigma^{2}_{*}(\nu)$, $\sigma^{2}_{**}(\nu)$ and $\sigma^{2}_{***}(\nu)$ for
$\sigma^{2}=\sigma^{2}_{2,2}$ and the effect of the deme-scaled extinction
rate $\nu$ on these values when $\sigma^{2}$ and the population-scaled means
of all the payoffs are of the same small order, while
$\sigma^{2}_{1,1}=\sigma^{2}_{1,2}=\sigma^{2}_{2,1}=0$. By using Eqs.
(6)-(73), the threshold values are found to be

$\displaystyle\sigma^{2}_{*}(\nu)$
$\displaystyle=\frac{5(6+\nu)(10+\nu)}{3(231+174\nu+35\nu^{2}+2\nu^{3})}\Big{[}(4+\nu)(\mu_{2,1}-\mu_{1,1})+(5+2\nu)(\mu_{2,2}-\mu_{1,2})\Big{]},$
(96a) $\displaystyle\sigma^{2}_{**}(\nu)$
$\displaystyle=\frac{2(6+\nu)}{5+\nu}\Big{(}\mu_{2,1}+\mu_{2,2}-\mu_{1,1}-\mu_{1,2}\Big{)},$
(96b) $\displaystyle\sigma^{2}_{***}(\nu)$
$\displaystyle=\frac{10(6+\nu)(10+\nu)}{3(288+127\nu+20\nu^{2}+\nu^{3})}\Big{[}(5+2\nu)(\mu_{2,1}-\mu_{1,1})+(4+\nu)(\mu_{2,2}-\mu_{1,2})\Big{]}.$
(96c)

It a worth noting, which can be easily checked, that the coefficients of
$\mu_{2,1}-\mu_{1,1}$ and $\mu_{2,2}-\mu_{1,2}$ in the first two expressions
above are decreasing functions with respect to $\nu$ in $(0,\infty)$, while
those on the third are increasing functions. Therefore, in the case where
$\mu_{2,1}-\mu_{1,1}>0$ and $\mu_{2,2}-\mu_{1,2}>0$ as in a Prisoner’s
dilemma, the conditions for $F_{1}>(Nd)^{-1}$ and $F_{1}>F_{2}$, given by
$\sigma^{2}>\sigma^{2}_{*}(\nu)>0$ and $\sigma^{2}>\sigma^{2}_{**}(\nu)>0$,
respectively, become less stringent as $\nu$ increases, while the condition
for $F_{2}<(Nd)^{-1}$, given by $\sigma^{2}>\sigma^{2}_{***}(\nu)>0$, becomes
more stringent. Furthermore, in this case, we have that

$\displaystyle\sigma^{2}_{***}(\nu)-\sigma^{2}_{**}(\nu)$
$\displaystyle=\frac{2(1+\nu)(6+\nu)\Big{[}(136+33\nu+2\nu^{2})(\mu_{2,1}-\mu_{1,1})+(386+108\nu+7\nu^{2})(\mu_{2,2}-\mu_{1,2})\Big{]}}{3(5+\nu)(288+127\nu+20\nu^{2}+\nu^{3})}$
(97)

and

$\displaystyle\sigma^{2}_{**}(\nu)-\sigma^{2}_{*}(\nu)$
$\displaystyle=\frac{(1+\nu)(6+\nu)\Big{[}(386+108\nu+7\nu^{2})(\mu_{2,1}-\mu_{1,1})+(136+33\nu+2\nu^{2})(\mu_{2,2}-\mu_{1,2})\Big{]}}{3(5+\nu)(231+174\nu+35\nu^{2}+2\nu^{3})}$
(98)

are positive, from which

$\sigma^{2}_{*}(\nu)<\sigma^{2}_{**}(\nu)<\sigma^{2}_{***}(\nu),$ (99)

for $\nu>0$. These inequalities imply that the conditions for
$F_{1}<(Nd)^{-1}$ and $F_{2}<(Nd)^{-1}$ cannot occur simultaneously for any
$\sigma_{2}^{2}>0$ and any $\nu>0$.

The above properties are illustrated in Figure 1 in the special case where the
population-scaled means of the payoffs are given by
$\mu_{1,1}=\mu_{b}-\mu_{c}$, $\mu_{1,2}=-\mu_{c}$, $\mu_{2,1}=\mu_{b}$ and
$\mu_{2,2}=0$, where $\mu_{b}>\mu_{c}>0$, which corresponds to an additive
Prisoner’s dilemma.

Figure 1: The curves represent the threshold values
$\sigma^{2}_{*}(\nu)/\mu_{c}$, $\sigma^{2}_{**}(\nu)/\mu_{c}$ and
$\sigma^{2}_{***}(\nu)/\mu_{c}$ with respect to the deme-scaled extinction
rate $\nu$ in Case 1 where
$\sigma^{2}_{1,1}=\sigma^{2}_{1,2}=\sigma^{2}_{2,1}=0$,
$\sigma^{2}_{2,2}=\sigma^{2}>0$ with $\mu_{1,1}=\mu_{b}-\mu_{c}$,
$\mu_{1,2}=-\mu_{c}$, $\mu_{2,1}=\mu_{b}$, $\mu_{2,2}=0$ and
$\mu_{b}>\mu_{c}>0$. These curves indicate the minimum values that
$\sigma^{2}/\mu_{c}$ must exceed in order to have $F_{1}>(Nd)^{-1}$,
$F_{1}>F_{2}$ and $F_{2}<(Nd)^{-1}$, respectively. Since the inequalities
$\sigma^{2}_{*}(\nu)<\sigma^{2}_{**}(\nu)<\sigma^{2}_{***}(\nu)$ hold for
every $\nu>0$, the scenario where $F_{1}<(Nd)^{-1}$ and $F_{2}<(Nd)^{-1}$
cannot occur.

Let us summarize our findings in Case 1.

##### Result 2.

If the population-scaled variance of the payoff to $S_{2}$ against $S_{2}$,
represented by $\sigma^{2}$, is large enough, while the population-scaled
variances of the other payoffs are negligible, then the inequalities
$F_{1}>(Nd)^{-1}$, $F_{1}>F_{2}$ and $F_{2}<(Nd)^{-1}$ will hold. Moreover,
when the population-scaled means of the payoffs to $S_{1}$ against $S_{1}$ and
$S_{2}$ are less than the corresponding ones to $S_{2}$ and all the
population-scaled parameters are of the same small order, these inequalities
will occur in the listed order as $\sigma^{2}$ increases, earlier for the
first two and later in the third one as the deme-scaled extinction rate $\nu$
increases.

##  8 Case 2: $\sigma^{2}_{1,1}=\sigma^{2}_{1,2}=\sigma^{2}_{2,2}=0$ and
$\sigma^{2}_{2,1}=\sigma^{2}>0$

Proceeding as previously in the case where the population-scaled variances of
the payoffs different from $\sigma^{2}_{2,1}=\sigma^{2}>0$ are insignificant,
it can be shown that an increase in $\sigma^{2}$ increases $F_{1}$ and
decreases $F_{2}$. Moreover, there exist three threshold values such that

$\displaystyle F_{1}>(Nd)^{-1}\mbox{ as soon as
}\sigma^{2}>\sigma^{2}_{*}(\nu),$ (100a) $\displaystyle F_{1}>F_{2}\mbox{ as
soon as }\sigma^{2}>\sigma^{2}_{**}(\nu),$ (100b) $\displaystyle
F_{2}<(Nd)^{-1}\mbox{ as soon as }\sigma^{2}>\sigma^{2}_{***}(\nu).$ (100c)

If $\sigma^{2}$ and the population-scaled means of all the payoffs are of the
same small order, then the threshold values are given by

$\displaystyle\sigma^{2}_{*}(\nu)$
$\displaystyle=\frac{5(6+\nu)(10+\nu)}{243+127\nu+20\nu^{2}+\nu^{3}}\Big{[}(4+\nu)(\mu_{2,1}-\mu_{1,1})+(5+2\nu)(\mu_{2,2}-\mu_{1,2})\Big{]},$
(101a) $\displaystyle\sigma^{2}_{**}(\nu)$
$\displaystyle=\frac{6(6+\nu)}{7+\nu}\Big{(}\mu_{2,1}+\mu_{2,2}-\mu_{1,1}-\mu_{1,2}\Big{)},$
(101b) $\displaystyle\sigma^{2}_{***}(\nu)$
$\displaystyle=\frac{10(6+\nu)(10+\nu)}{3(188+117\nu+20\nu^{2}+\nu^{3})}\Big{[}(5+2\nu)(\mu_{2,1}-\mu_{1,1})+(4+\nu)(\mu_{2,2}-\mu_{1,2})\Big{]}.$
(101c)

Even in the case where $\mu_{2,1}-\mu_{1,1}>0$ and $\mu_{2,2}-\mu_{1,2}>0$,
the order of these threshold values is not clear. In the particular case where
$\mu_{1,1}=\mu_{b}-\mu_{c}$, $\mu_{1,2}=-\mu_{c}$, $\mu_{2,1}=\mu_{b}$ and
$\mu_{2,2}=0$ with $\mu_{b}>\mu_{c}>0$, however, we have

$\displaystyle\sigma^{2}_{*}(\nu)-\sigma^{2}_{**}(\nu)$
$\displaystyle=\frac{3(1+\nu)(6+\nu)^{2}(13+\nu)}{(7+\nu)(243+127\nu+20\nu^{2}+\nu^{3})}\mu_{c}>0,$
(102a) $\displaystyle\sigma^{2}_{**}(\nu)-\sigma^{2}_{***}(\nu)$
$\displaystyle=\frac{2(1+\nu)(6+\nu)^{2}(13+\nu)}{(7+\nu)(188+117\nu+20\nu^{2}+\nu^{3})}\mu_{c}>0,$
(102b)

so that

$\sigma^{2}_{***}(\nu)<\sigma^{2}_{**}(\nu)<\sigma^{2}_{*}(\nu),$ (103)

for $\nu>0$. These inequalities imply that the conditions for
$F_{1}>(Nd)^{-1}$ and $F_{2}>(Nd)^{-1}$ cannot simultaneously occur. Moreover,
while $\sigma^{2}_{*}(\nu)$ and $\sigma^{2}_{**}(\nu)$ are monotonic
functions, actually monotonically increasing functions, with respect to $\nu$
in the whole interval $(0,\infty)$, this is not necessarily the case for
$\sigma^{2}_{***}(\nu)$ which can decrease first and then increase to a global
maximum as $\nu$ increases. This is illustrated in Figure 2.

We are now ready to state our next result.

##### Result 3.

If the population-scaled variance of the payoff to $S_{2}$ against $S_{1}$,
represented by $\sigma^{2}$, is large enough, while the population-scaled
variances of the other payoffs are negligible, then the inequalities
$F_{1}>(Nd)^{-1}$, $F_{1}>F_{2}$ and $F_{2}<(Nd)^{-1}$ will hold. Moreover,
when all the population-scaled parameters are of the same small order, these
inequalities may occur in the reversed order as $\sigma^{2}$ increases, and
all later as the deme-scaled extinction rate $\nu$ tends to infinity.

Figure 2: The curves represent the minimum threshold values
$\sigma^{2}_{*}(\nu)/\mu_{c}$, $\sigma^{2}_{**}(\nu)/\mu_{c}$ and
$\sigma^{2}_{***}(\nu)/\mu_{c}$ for $\sigma^{2}/\mu_{c}$ to have
$F_{1}>(Nd)^{-1}$, $F_{1}>F_{2}$ and $F_{2}<(Nd)^{-1}$, respectively, in Case
2 where $\sigma^{2}_{1,1}=\sigma^{2}_{1,2}=\sigma^{2}_{2,2}=0$,
$\sigma^{2}_{2,1}=\sigma^{2}>0$ with $\mu_{1,1}=\mu_{b}-\mu_{c}$,
$\mu_{1,2}=-\mu_{c}$, $\mu_{2,1}=\mu_{b}$, $\mu_{2,2}=0$ and
$\mu_{b}>\mu_{c}>0$. Note that the worst scenario for the three conditions to
hold occurs when the deme-scaled extinction rate $\nu\rightarrow\infty$.

##  9 Case 3: $\sigma^{2}_{1,1}=\sigma^{2}_{1,2}=0$ and
$\sigma^{2}_{2,2}=\sigma^{2}_{2,1}=\sigma^{2}>0$

Finally, in the case where the population-scaled variances of the payoffs to
$S_{1}$ against $S_{1}$ and $S_{2}$ are negligible compared to those to
$S_{2}$ given by $\sigma^{2}>0$, it can be shown as previously that $F_{1}$
increases and $F_{2}$ decreases as $\sigma^{2}$ increases such that
$F_{1}>(Nd)^{-1}$, $F_{1}>F_{2}$ and $F_{2}<(Nd)^{-1}$ occur as soon as
$\sigma^{2}>\sigma^{2}_{*}(\nu),$ $\sigma^{2}_{**}(\nu)$ and
$\sigma^{2}_{***}(\nu)$, respectively. Moreover, if the population-scaled
parameters are of the same small order, then the threshold values are given by

$\displaystyle\sigma^{2}_{*}(\nu)$
$\displaystyle=\frac{5(6+\nu)(10+\nu)}{936+649\nu+125\nu^{2}+7\nu^{3}}\Big{[}(4+\nu)(\mu_{2,1}-\mu_{1,1})+(5+2\nu)(\mu_{2,2}-\mu_{1,2})\Big{]},$
(104a) $\displaystyle\sigma^{2}_{**}(\nu)$
$\displaystyle=\frac{3(6+\nu)}{11+2\nu}\Big{(}\mu_{2,1}+\mu_{2,2}-\mu_{1,1}-\mu_{1,2}\Big{)},$
(104b) $\displaystyle\sigma^{2}_{***}(\nu)$
$\displaystyle=\frac{5(6+\nu)(10+\nu)}{3(238+122\nu+20\nu^{2}+\nu^{3})}\Big{[}(5+2\nu)(\mu_{2,1}-\mu_{1,1})+(4+\nu)(\mu_{2,2}-\mu_{1,2})\Big{]}.$
(104c)

Here, the coefficients of $\mu_{2,1}-\mu_{1,1}$ and $\mu_{2,2}-\mu_{1,2}$ in
the expressions of $\sigma^{2}_{*}(\nu)$ and $\sigma^{2}_{**}(\nu)$ are
decreasing functions on $(0,\infty)$, while the first coefficient is
increasing and the second decreasing in the expression of
$\sigma^{2}_{***}(\nu)$. Moreover, when $\mu_{2,1}>\mu_{1,1}$ and
$\mu_{2,2}>\mu_{1,2}$, the differences

$\displaystyle\sigma^{2}_{**}(\nu)-\sigma^{2}_{*}(\nu)$
$\displaystyle\quad=\frac{(1+\nu)(6+\nu)\Big{[}(608+169\nu+11\nu^{2})(\mu_{2,1}-\mu_{1,1})+(58+14\nu+\nu^{2})(\mu_{2,2}-\mu_{1,2})\Big{]}}{(11+2\nu)(936+649\nu+125\nu^{2}+7\nu^{3})}$
(105)

and

$\displaystyle\sigma^{2}_{***}(\nu)-\sigma^{2}_{**}(\nu)$
$\displaystyle\quad=\frac{(1+\nu)(6+\nu)\Big{[}(608+169\nu+11\nu^{2})(\mu_{2,1}-\mu_{1,1})+(58+14\nu+\nu^{2})(\mu_{2,2}-\mu_{1,2})\Big{]}}{3(11+2\nu)(238+122\nu+20\nu^{2}+\nu^{3})}$
(106)

are both positive, from which we have

$\sigma^{2}_{*}(\nu)<\sigma^{2}_{**}(\nu)<\sigma^{2}_{***}(\nu),$ (107)

for $\nu>0$. This is similar to Case 1.

Figure 3 shows the curves of $\sigma^{2}_{*}(\nu)/\mu_{c}$,
$\sigma^{2}_{**}(\nu)/\mu_{c}$ and $\sigma^{2}_{***}(\nu)\mu_{c}$ in the
particular case where $\mu_{1,1}=\mu_{b}-\mu_{c}$, $\mu_{1,2}=-\mu_{c}$,
$\mu_{2,1}=\mu_{b}$ and $\mu_{2,2}=0$ with $\mu_{b}>\mu_{c}>0$. There is
little difference with Figure 1 in Case 1, which suggests that the variability
in the payoff to $S_{2}$ against $S_{2}$ is more important than the
variability in the payoff to $S_{2}$ against $S_{1}$ in conditions for the
evolution of $S_{1}$.

##### Result 4.

The conditions for $F_{1}>(Nd)^{-1}$, $F_{1}>F_{2}$ and $F_{2}<(Nd)^{-1}$ in
Case 3, where the population-scaled variances of the payoffs to $S_{2}$
against $S_{1}$ and $S_{2}$ are equal to $\sigma^{2}>0$, while the population-
scaled variances of the payoffs to $S_{1}$ against $S_{1}$ and $S_{2}$ are
negligible in comparison, are similar to the conditions obtained in Case 1.

Figure 3: The plotted curves represent the minimum threshold values
$\sigma^{2}_{*}(\nu)/\mu_{c}$, $\sigma^{2}_{**}(\nu)/\mu_{c}$ and
$\sigma^{2}_{***}(\nu)/\mu_{c}$ for $\sigma^{2}/\mu_{c}$ to have
$F_{1}>(Nd)^{-1}$, $F_{1}>F_{2}$ and $F_{2}<(Nd)^{-1}$, respectively, in Case
3 where $\sigma^{2}_{1,1}=\sigma^{2}_{1,2}=0$,
$\sigma^{2}_{2,2}=\sigma^{2}_{2,1}=\sigma^{2}>0$ with
$\mu_{1,1}=\mu_{b}-\mu_{c}$, $\mu_{1,2}=-\mu_{c}$, $\mu_{2,1}=\mu_{b}$,
$\mu_{2,2}=0$ and $\mu_{b}>\mu_{c}>0$. Like in Case 1, the best scenario for
the first two conditions to hold occurs when $\nu\rightarrow\infty$, while the
best scenario for the third condition occurs when $\nu\rightarrow 0$.

##  10 Discussion

Previous studies examining the island model primarily focused on constant
payoffs. This paper examines a haploid population that is subdivided into $d$
demes, each containing $N\geq 2$ individuals. The payoffs received by
offspring using the strategies $S_{1}$ and $S_{2}$ in random pairwise
interactions within demes are subject to stochastic fluctuations from one
generation to the next such that their means and variances are proportional to
the inverse of the population size $Nd$. Following viability selection within
demes and allowing for local extinction and uniform recolonization with
probability $m\in(0,1)$, the next generation starts with $N$ offspring chosen
at random in each deme according to a Wright-Fisher model.

Having examined the conditions outlined in Ethier and Nagylaki ([11]), we have
shown that the exact discrete-time Markov chain, which is determined by the
deme type frequencies, involves two time scales so that it converges to a
diffusion process for the strategy frequencies as $d$ goes to infinity. This
convergence allows us to compute the probability for strategy $S_{1}$ to
become fixed in the population following its introduction as a single mutant,
represented by $F_{1}$, and the corresponding fixation probability for
$S_{2}$, represented by $F_{2}$. As a result, we have examined how stochastic
fluctuations in the payoffs impact the evolutionary dynamics. It is important
to highlight that the infinitesimal mean and variance of the limiting
diffusion process are expressed in terms of the population-scaled means and
variances of the payoffs, represented by $\mu_{k,l}$ and $\sigma^{2}_{k,l}$,
respectively, for $k,l=1,2$, with coefficients that depend on identity-by-
descent measures between offspring under neutrality besides. These have been
deduced exactly and shown to depend only on the deme-scaled extinction rate
$\nu=Nm$ for $N$ large enough and $m$ small enough (see Appendices H and J).
It is worth noting that if $\nu\rightarrow\infty$, the limiting diffusion
process is exactly the one obtained in a large well-mixed population (Li and
Lessard [34]).

In the particular context of the deterministic Prisoner’s dilemma with $S_{1}$
for cooperation and $S_{2}$ for defection, which means that
$\sigma^{2}_{1,1}=\sigma^{2}_{1,2}=\sigma^{2}_{2,1}=\sigma^{2}_{2,2}=0$ with
$\mu_{1,1}<\mu_{2,1}$ and $\mu_{1,2}<\mu_{2,2}$, our findings indicate that
selection fully disfavours the evolution of cooperation in the sense that
$F_{1}<(Nd)^{-1}<F_{2}$ for any $\nu>0$. Here $(Nd)^{-1}$ is the probability
of ultimate fixation of a new mutant under neutrality.

We have investigated the effects of fluctuations in the payoffs for defection,
represented by $S_{2}$, across various scenarios. It is reasonable to expect
that there would be greater uncertainty in the payoffs associated with
defection compared to cooperation. Several factors contribute to this
difference, one being that defectors often find themselves more isolated than
cooperators. Consequently, they may not benefit or suffer as much from
fluctuations in the surrounding resources caused by environmental variations.
Furthermore, defectors may face punishments or a lack of rewards from other
individuals in the population, adding to the uncertainty about their payoffs.
While defection might offer higher expected payoffs compared to cooperation,
this increased isolation and exposure to external factors can also result in
higher variances in their payoffs. As a consequence, this elevated variance in
the payoffs for defection can significantly influence the evolutionary
dynamics and the overall prevalence of cooperative strategies in the
population.

In the first scenario, we assume that all the population-scaled variances of
the payoffs are negligible compared to $\sigma_{2,2}^{2}$, which represents
the population-scaled variance of the payoff to $S_{2}$ against $S_{2}$. As we
increase $\sigma_{2,2}^{2}$, the fixation probability $F_{1}$ increases, while
the fixation probability $F_{2}$ decreases for any deme-scaled extinction rate
$\nu>0$. Additionally, we have identified three thresholds, namely
$\sigma^{2}_{*}(\nu)$, $\sigma^{2}_{**}(\nu)$ and $\sigma^{2}_{***}(\nu)$,
such that the conditions $F_{1}>(Nd)^{-1}$, $F_{1}>F_{2}$, and
$F_{2}<(Nd)^{-1}$ hold when $\sigma_{2,2}^{2}>\sigma^{2}_{*}(\nu)$,
$\sigma_{2,2}^{2}>\sigma^{2}_{**}(\nu)$ and
$\sigma_{2,2}^{2}>\sigma^{2}_{***}(\nu)$, respectively. Furthermore, when all
the population-scaled means of the payoffs, represented by
$\mu_{1,1},\mu_{1,2},\mu_{2,1},\mu_{2,2}$, are of the same small order, these
thresholds satisfy the inequalities
$\sigma^{2}_{*}(\nu)<\sigma^{2}_{**}(\nu)<\sigma^{2}_{***}(\nu)$. In the case
where the population-scaled means determine an additive prisoner’s dilemma
($\mu_{1,1}=\mu_{b}-\mu_{c}$, $\mu_{1,2}=-\mu_{c}$, $\mu_{2,1}=\mu_{b}$, and
$\mu_{2,2}=0$), we have shown that $\sigma^{2}_{*}(\nu)$ and
$\sigma^{2}_{**}(\nu)$ are decreasing with respect to $\nu>0$ while
$\sigma^{2}_{***}(\nu)$ is increasing. This means that increasing the deme-
scaled extinction rate will lessen the conditions for selection to favour the
evolution of $S_{1}$ and to favour the evolution of $S_{1}$ more than the
evolution of $S_{2}$. Conversely, this will reinforce the condition for
selection to disfavour the evolution of $S_{2}$ (see Figure 1). Under such a
scenario, it is impossible for selection to disfavour the evolution of both
cooperation and defection for any $\sigma_{2,2}^{2}>0$ and $\nu>0$.

In the second scenario, we consider that all the population-scaled variances
of the payoffs other than $\sigma_{2,1}^{2}$, for the payoff to $S_{2}$
against $S_{1}$, are insignificant. Similarly to the first scenario, we have
shown that an increase in $\sigma_{2,1}^{2}$ will increase the fixation
probability $F_{1}$ and decrease the fixation probability $F_{2}$. Moreover,
this scenario will make it possible for selection to favour the evolution of
$S_{1}$ in any sense. More precisely, there exist three thresholds
$\sigma^{2}_{*}(\nu)$, $\sigma^{2}_{**}(\nu)$ and $\sigma^{2}_{***}(\nu)$,
where $F_{1}>(Nd)^{-1}$, $F_{1}>F_{2}$ and $F_{2}<(Nd)^{-1}$ as long as
$\sigma_{2,1}^{2}>\sigma^{2}_{*}(\nu)$,
$\sigma_{2,1}^{2}>\sigma^{2}_{**}(\nu)$ and
$\sigma_{2,1}^{2}>\sigma^{2}_{***}(\nu)$, respectively. Here, there is no
clear relationship between these thresholds. However, when the population-
scaled means of the payoffs are of the same small order and determine an
additive Prisoner’s dilemma, we find that
$\sigma^{2}_{***}(\nu)<\sigma^{2}_{**}(\nu)<\sigma^{2}_{*}(\nu)$ for any
$\nu>0$. Furthermore, the functions $\sigma^{2}_{*}(\nu)$ and
$\sigma^{2}_{**}(\nu)$ are increasing, while $\sigma^{2}_{***}(\nu)$ is
decreasing. This means that increasing the deme-scaled extinction rate will
strengthen the conditions for $F_{1}>(Nd)^{-1}$ and $F_{1}>F_{2}$ and lessen
the condition for $F_{2}<(Nd)^{-1}$ (see Figure 2). Under such a scenario, it
is impossible for selection to favour the evolution of both cooperation and
defection for any $\sigma_{2,2}^{2}>0$ and $\nu>0$.

In the third scenario, we assume that the population-scaled variances of the
payoffs to $S_{1}$ are insignificant compared to the population-scaled
variances of the payoffs to $S_{2}$ against $S_{1}$ and to $S_{2}$ against
$S_{2}$, given by $\sigma^{2}$. It is important to highlight that increasing
$\sigma^{2}$ will lead to an increase in $F_{1}$ and a decrease in $F_{2}$. In
addition, it is crucial to note that all the conclusions from the first
scenario remain valid (refer to Figure 3). By comparing Figures 1, 2 and 3, it
becomes evident that for increasing $F_{1}$, it is more effective to increase
the population-scaled variances of the payoffs to both $S_{2}$ against $S_{1}$
and $S_{2}$ against $S_{1}$ simultaneously, rather than concentrating the
increase on just one of these variances.

In general, assuming that all population-scaled means and variances of the
payoffs are of the same small order, we have shown that decreasing the
population-scaled variance of any payoff to $S_{1}$ or increasing the
population-scaled variance of any payoff to $S_{2}$ increases the fixation
probability $F_{1}$ and decreases the fixation probability $F_{2}$. Therefore,
either increases the difference $F_{1}-F_{2}$. As the deme-scaled extinction
rate $\nu$ decreases, which rises the level of identity by descent among
offspring within demes, the weight of $\sigma_{2,1}^{2}-\sigma_{1,2}^{2}$ in
this difference increases, while the weight of
$\sigma_{2,2}^{2}-\sigma_{1,1}^{2}$ decreases. This may be counterintuitive,
since then offspring of the same type are more likely to interact than
offspring of different types. But we have to remind ourselves that then
offspring are also more likely to be in competition with offspring of the same
type in the same deme. Note also that $F_{1}-F_{2}$ does not depend on $\nu$
in the case of all deterministic payoffs.

Our results extend previous studies that have explored the stochastic
Prisoner’s dilemma in a large but finite well-mixed population which
corresponds to the case where $\nu\rightarrow\infty$ (see Li and Lessard [34]
for a Wright-Fisher model and by Kroumi et al. [27] for a Moran model). In
these studies, it was shown that more variability in the payoffs for defection
than in the payoffs for cooperation is required for selection to favour the
evolution of cooperation in any sense. Note that a similar conclusion was
drawn in the context of an infinite well-mixed population (Zheng et al. [64,
63]). In this setting, the same condition was found to be necessary for the
fixation of cooperation to be stochastically locally stable, and for the
fixation of defection to be stochastically locally unstable. In an island
model, the weights of the variances of the different payoffs may differ, but
the general conclusion still holds.

Finally, the island model with the possibility of extinction and
recolonization can be extended to study groups distributional dynamics
stemming from the possibility of opting out in repeated games (Cressman and
Krivan [9], Krivan and Cressman [21], Li and Lessard [35], Zhang _et al._
[62]). In the general case, however, the extinction probability or rate
depends on the group composition, which makes the analysis more challenging.

## Appendices

### Appendix A: Proof of Lemma 1

By using Taylor’s theorem with the Lagrange form of the remainder, there
exists a real random variable $\xi_{i}$ such that
$|\xi_{i}|<|\tilde{a}_{i}|<1$ and

$\frac{1}{1+\tilde{a}_{i}}=1-\tilde{a}_{i}+\tilde{a}_{i}^{2}-\frac{\tilde{a}_{i}^{3}}{(1+\xi_{i})^{3}},$
(108)

from which we deduce

$\displaystyle E_{\mathbf{z}}\left[\tilde{x}\right]$
$\displaystyle=\sum_{i=0}^{N}z_{i}x_{i}E_{\mathbf{z}}\left[\frac{1+\tilde{a}_{1,i}}{1+\tilde{a}_{i}}\right]$
$\displaystyle=\sum_{i=0}^{N}z_{i}x_{i}E_{\mathbf{z}}\left[\left(1+\tilde{a}_{1,i}\right)\left(1-\tilde{a}_{i}+\tilde{a}_{i}^{2}\right)\right]+o(d^{-1})$
$\displaystyle=x+\sum_{i=0}^{N}z_{i}x_{i}E_{\mathbf{z}}\left[\tilde{a}_{1,i}-\tilde{a}_{i}\right]-\sum_{i=0}^{N}z_{i}x_{i}E_{\mathbf{z}}\left[\left(\tilde{a}_{1,i}-\tilde{a}_{i}\right)\tilde{a}_{i}\right]+o(d^{-1}).$
(109)

On the other hand, we have

$\displaystyle\sum_{i=0}^{N}z_{i}x_{i}E_{\mathbf{z}}\left[\tilde{a}_{1,i}-\tilde{a}_{i}\right]$
$\displaystyle=\sum_{i=0}^{N}z_{i}x_{i}E_{\mathbf{z}}\left[\tilde{a}_{1,i}\right]-\sum_{i=0}^{N}z_{i}x_{i}\sum_{j=0}^{N}z_{j}E_{\mathbf{z}}\left[\tilde{a}_{j}\right]-\left(\sum_{i=0}^{N}z_{i}x_{i}E_{\mathbf{z}}\left[\tilde{a}_{i}\right]-\sum_{i=0}^{N}z_{i}x_{i}\sum_{j=0}^{N}z_{j}E_{\mathbf{z}}\left[\tilde{a}_{j}\right]\right)$
$\displaystyle=Cov_{\mathbf{z}}(q_{1,I_{1}},a_{I_{1}})-Cov_{\mathbf{z}}(q_{1,I_{1}},a_{I_{2}})$
(110)

and

$\displaystyle\sum_{i=0}^{N}z_{i}x_{i}E_{\mathbf{z}}\left[\left(\tilde{a}_{1,i}-\tilde{a}_{i}\right)\tilde{a}_{i}\right]$
$\displaystyle=\sum_{i=0}^{N}z_{i}x_{i}E_{\mathbf{z}}\left[\tilde{a}_{1,i}\tilde{a}_{i}\right]-\sum_{i=0}^{N}z_{i}x_{i}\sum_{j=0}^{N}z_{j}E_{\mathbf{z}}\left[\tilde{a}_{j}^{2}\right]-\left(\sum_{i=0}^{N}z_{i}x_{i}E_{\mathbf{z}}\left[\tilde{a}_{i}^{2}\right]-\sum_{i=0}^{N}z_{i}x_{i}\sum_{j=0}^{N}z_{j}E_{\mathbf{z}}\left[\tilde{a}_{j}^{2}\right]\right)$
$\displaystyle=Cov_{\mathbf{z}}(q_{1,I_{1}},a_{I_{1}}a_{I_{2}})-Cov_{\mathbf{z}}(q_{1,I_{1}},a_{I_{2}}a_{I_{3}}),$
(111)

where $I_{1}$, $I_{2}$, and $I_{3}$ are randomly chosen offspring produced in
the same deme chosen at random before viability selection and recolonization.
Inserting Eqs. (Appendix A: Proof of Lemma 1) and (Appendix A: Proof of Lemma
1) in Eq. (Appendix A: Proof of Lemma 1) and then using the expression for the
payoff of an offspring in Eq. (3) and the assumptions in Eq. (2) lead to the
stated result.

### Appendix B: Proof of Lemma 2

By a similar argument as in the previous proof, we have

$\begin{split}E_{\mathbf{z}}[\tilde{\tilde{x}}]&=\sum_{i=0}^{N}z_{i}x_{i}E_{\mathbf{z}}\left[\frac{1+\tilde{a}_{1,i}}{1+\tilde{a}}\right]\\\
&=\sum_{i=0}^{N}z_{i}x_{i}E_{\mathbf{z}}\left[\left(1+\tilde{a}_{1,i}\right)\left(1-\tilde{a}+\tilde{a}^{2}\right)\right]+o(d^{-1})\\\
&=x+\sum_{i=0}^{N}z_{i}x_{i}E_{\mathbf{z}}\left[\tilde{a}_{1,i}-\tilde{a}-(\tilde{a}_{1,i}-\tilde{a})\tilde{a}\right]+o(d^{-1}).\end{split}$
(112)

On the other hand, we have the identities

$\displaystyle Cov_{\mathbf{z}}(q_{1,I_{1}},a_{I_{1}})$
$\displaystyle=\sum_{i=0}^{N}z_{i}x_{i}E_{\mathbf{z}}\left[\tilde{a}_{1,i}-\tilde{a}\right],$
(113a) $\displaystyle Cov_{\mathbf{z}}(q_{1,I_{1}},a_{I_{1}}a_{J_{1}})$
$\displaystyle=\sum_{i=0}^{N}z_{i}x_{i}E_{\mathbf{z}}\left[(\tilde{a}_{1,i}-\tilde{a})\tilde{a}\right],$
(113b)

where $I_{1}$ and $J_{1}$ are randomly chosen offspring produced in demes
chosen at random and independently before viability selection and
recolonization. The stated result is obtained by substituting these identities
into Eq. (112) and then using Eqs. (3 ) and (2).

### Appendix C: Proof of Lemma 3

Proceeding as in the proof of Lemma 1, we have

$\displaystyle E_{\mathbf{z}}\left[\tilde{x}^{2}\right]$
$\displaystyle=\sum_{i,j=0}^{N}z_{i}z_{j}x_{i}x_{j}E_{\mathbf{z}}\left[\frac{1+\tilde{a}_{1,i}}{1+\tilde{a}_{i}}\frac{1+\tilde{a}_{1,j}}{1+\tilde{a}_{j}}\right]$
$\displaystyle=\sum_{i,j=0}^{N}z_{i}z_{j}x_{i}x_{j}E_{\mathbf{z}}\Big{[}\left(1+\tilde{a}_{1,i}\right)\left(1+\tilde{a}_{1,j}\right)\left(1-\tilde{a}_{i}+\tilde{a}_{i}^{2}\right)\left(1-\tilde{a}_{j}+\tilde{a}_{j}^{2}\right)\Big{]}+o(d^{-1})$
$\displaystyle=x^{2}+\sum_{i,j=0}^{N}z_{i}z_{j}x_{i}x_{j}E_{\mathbf{z}}\Big{[}\tilde{a}_{1,i}-\tilde{a}_{i}+\tilde{a}_{1,j}-\tilde{a}_{j}-\left(\tilde{a}_{1,i}-\tilde{a}_{i}\right)\tilde{a}_{i}-\left(\tilde{a}_{1,j}-\tilde{a}_{j}\right)\tilde{a}_{j}$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad+\left(\tilde{a}_{1,i}-\tilde{a}_{i}\right)\tilde{a}_{1,j}-\left(\tilde{a}_{1,i}-\tilde{a}_{i}\right)\tilde{a}_{j}\Big{]}+o(d^{-1}).$
(114)

Moreover, we have the identities

$\displaystyle\sum_{i,j=0}^{N}z_{i}z_{j}x_{i}x_{j}E_{\mathbf{z}}\left[\tilde{a}_{1,i}-\tilde{a}_{i}\right]$
$\displaystyle=x\Big{[}Cov_{\mathbf{z}}\left(q_{1,I_{1}},a_{I_{1}}\right)-Cov_{\mathbf{z}}\left(q_{1,I_{1}},a_{I_{2}}\right)\Big{]},$
(115a)
$\displaystyle\sum_{i,j=0}^{N}z_{i}z_{j}x_{i}x_{j}E_{\mathbf{z}}\left[(\tilde{a}_{1,i}-\tilde{a}_{i})\tilde{a}_{i}\right]$
$\displaystyle=x\Big{[}Cov_{\mathbf{z}}\left(q_{1,I_{1}},a_{I_{1}}a_{I_{2}}\right)-Cov_{\mathbf{z}}\left(q_{1,I_{1}},a_{I_{2}}a_{I_{3}}\right)\Big{]},$
(115b)
$\displaystyle\sum_{i,j=0}^{N}z_{i}z_{j}x_{i}x_{j}E_{\mathbf{z}}\left[(\tilde{a}_{1,i}-\tilde{a}_{i})\tilde{a}_{1,j}\right]$
$\displaystyle=Cov_{\mathbf{z}}\left(q_{1,I_{1}},q_{1,J_{1}}a_{I_{1}}a_{J_{1}}\right)-Cov_{\mathbf{z}}\left(q_{1,I_{1}},q_{1,J_{1}}a_{I_{2}}a_{J_{1}}\right),$
(115c)
$\displaystyle\sum_{i,j=0}^{N}z_{i}z_{j}x_{i}x_{j}E_{\mathbf{z}}\left[(\tilde{a}_{1,i}-\tilde{a}_{i})\tilde{a}_{j}\right]$
$\displaystyle=Cov_{\mathbf{z}}\left(q_{1,I_{1}},q_{1,J_{1}}a_{I_{1}}a_{J_{2}}\right)-Cov_{\mathbf{z}}\left(q_{1,I_{1}},q_{1,J_{1}}a_{I_{2}}a_{J_{2}}\right),$
(115d)

where $I_{1}$, $I_{2}$ are two randomly chosen offspring produced in the same
deme chosen at random, and similarly for $J_{1}$, $J_{2}$ in an independent
way. Inserting these identities into Eq. (Appendix C: Proof of Lemma 3) yields

$\displaystyle E_{\mathbf{z}}\left[\tilde{x}^{2}\right]$
$\displaystyle=x^{2}+2x\Big{[}Cov_{\mathbf{z}}\left(q_{1,I_{1}},a_{I_{1}}\right)-Cov_{\mathbf{z}}\left(q_{1,I_{1}},a_{I_{2}}\right)\Big{]}-2x\Big{[}Cov_{\mathbf{z}}\left(q_{1,I_{1}},a_{I_{1}}a_{I_{2}}\right)$
$\displaystyle\quad\quad\quad-
Cov_{\mathbf{z}}\left(q_{1,I_{1}},a_{I_{2}}a_{I_{3}}\right)\Big{]}+Cov_{\mathbf{z}}\left(q_{1,I_{1}},q_{1,J_{1}}a_{I_{1}}a_{J_{1}}\right)+Cov_{\mathbf{z}}\left(q_{1,I_{1}},q_{1,J_{1}}a_{I_{2}}a_{J_{2}}\right)$
$\displaystyle\quad\quad\quad-
Cov_{\mathbf{z}}\left(q_{1,I_{1}},q_{1,J_{1}}a_{I_{2}}a_{J_{1}}\right)-Cov_{\mathbf{z}}\left(q_{1,I_{1}},q_{1,J_{1}}a_{I_{1}}a_{J_{2}}\right)+o\left(d^{-1}\right)$
$\displaystyle=x^{2}+\frac{1}{Nd}\Big{[}2x\Big{(}Cov_{\mathbf{z}}\left(q_{1,I_{1}},\mu_{I_{1}}\right)+Cov_{\mathbf{z}}\left(q_{1,I_{1}},\sigma_{I_{2},I_{3}}\right)-Cov_{\mathbf{z}}\left(q_{1,I_{1}},\mu_{I_{2}}\right)$
$\displaystyle\quad\quad\quad-
Cov_{\mathbf{z}}\left(q_{1,I_{1}},\sigma_{I_{1},I_{2}}\right)\Big{)}+Cov_{\mathbf{z}}\left(q_{1,I_{1}},q_{1,J_{1}}\sigma_{I_{1},J_{1}}\right)+Cov_{\mathbf{z}}\left(q_{1,I_{1}},q_{1,J_{1}}\sigma_{I_{2},J_{2}}\right)$
$\displaystyle\quad\quad\quad-
Cov_{\mathbf{z}}\left(q_{1,I_{1}},q_{1,J_{1}}\sigma_{I_{2},J_{1}}\right)-Cov_{\mathbf{z}}\left(q_{1,I_{1}},q_{1,J_{1}}\sigma_{I_{1},J_{2}}\right)\Big{]}+o\left(d^{-1}\right).$
(116)

Combining the last equation and Lemma 1 completes the proof of Lemma 3.

### Appendix D: Proof of Lemma 4

Proceeding as in the proof of Lemma 2, we obtain

$\displaystyle E_{\mathbf{z}}\left[\tilde{\tilde{x}}^{2}\right]$
$\displaystyle=\sum_{i,j=0}^{N}z_{i}z_{j}x_{i}x_{j}E\Big{[}\left(1+\tilde{a}_{1,i}\right)\left(1+\tilde{a}_{1,j}\right)\left(1-\tilde{a}+\tilde{a}^{2}\right)^{2}\Big{]}+o(d^{-1})$
$\displaystyle=x^{2}+\sum_{i,j=0}^{N}z_{i}z_{j}x_{i}x_{j}E_{\mathbf{z}}\Big{[}\tilde{a}_{1,i}-\tilde{a}+\tilde{a}_{1,j}-\tilde{a}-2\left(\tilde{a}_{1,i}-\tilde{a}\right)\tilde{a}$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad-2\left(\tilde{a}_{1,j}-\tilde{a}\right)\tilde{a}+\tilde{a}_{1,i}\tilde{a}_{1,j}-\tilde{a}\tilde{a}\Big{]}+o(d^{-1})$
$\displaystyle=x^{2}+\sum_{i,j=0}^{N}z_{i}z_{j}x_{i}x_{j}E_{\mathbf{z}}\Big{[}2(\tilde{a}_{1,i}-\tilde{a})-4(\tilde{a}_{1,j}-\tilde{a})\tilde{a}+\tilde{a}_{1,i}\tilde{a}_{1,j}-\tilde{a}^{2}\Big{]}+o(d^{-1})$
$\displaystyle=x^{2}+2x\,Cov_{\mathbf{z}}\left(q_{1,I_{1}},a_{I_{1}}\right)+Cov_{\mathbf{z}}\left(q_{1,I_{1}}q_{1,J_{1}},a_{I_{1}}a_{J_{1}}\right)-4x\,Cov_{\mathbf{z}}\left(q_{1,I_{1}},a_{I_{1}}a_{J_{1}}\right)+o(d^{-1})$
$\displaystyle=x^{2}+\frac{1}{Nd}\Big{[}2x\,Cov_{\mathbf{z}}\left(q_{1,I_{1}},\mu_{I_{1}}\right)+Cov_{\mathbf{z}}\left(q_{1,I_{1}}q_{1,J_{1}},\sigma_{I_{1},J_{1}}\right)$
$\displaystyle\quad\quad\quad\quad\quad-4x\,Cov_{\mathbf{z}}\left(q_{1,I_{1}},\sigma_{I_{1},J_{1}}\right)\Big{]}+o\left(d^{-1}\right),$
(117)

where $I_{1}$ and $J_{1}$ are randomly chosen offspring produced in demes
independently chosen at random. The proof of Lemma 4 can be completed by using
Lemma 2 in combination with the above equation.

### Appendix E: Proof of Lemma 5

Proceeding as in the proofs of Lemmas 3 and 4, we get

$\displaystyle E_{\mathbf{z}}\left[\tilde{x}\tilde{\tilde{x}}\right]$
$\displaystyle=\sum_{i,j=0}^{N}z_{i}z_{j}x_{i}x_{j}E_{\mathbf{z}}\Big{[}\left(1+\tilde{a}_{1,i}\right)\left(1+\tilde{a}_{1,j}\right)\left(1-\tilde{a}_{i}+\tilde{a}_{i}^{2}\right)\left(1-\tilde{a}+\tilde{a}^{2}\right)\Big{]}+o(d^{-1})$
$\displaystyle=x^{2}+\sum_{i,j=0}^{N}z_{i}z_{j}x_{i}x_{j}E_{\mathbf{z}}\Big{[}\tilde{a}_{1,i}-\tilde{a}_{i}+\tilde{a}_{1,j}-\tilde{a}-\left(\tilde{a}_{1,i}-\tilde{a}_{i}\right)\tilde{a}_{i}-\left(\tilde{a}_{1,j}-\tilde{a}\right)\tilde{a}$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad+\left(\tilde{a}_{1,i}-\tilde{a}_{i}\right)\tilde{a}_{1,j}-\left(\tilde{a}_{1,i}-\tilde{a}_{i}\right)\tilde{a}\Big{]}+o(d^{-1}).$
$\displaystyle=x^{2}+x\,Cov_{\mathbf{z}}(q_{1,I_{1}},a_{I_{1}})-x\,Cov_{\mathbf{z}}(q_{1,I_{1}},a_{I_{2}})+x\,Cov_{\mathbf{z}}(q_{1,I_{1}},a_{I_{2}}a_{I_{3}})$
$\displaystyle\quad-x\,Cov_{\mathbf{z}}(q_{1,I_{1}},a_{I_{1}}a_{I_{2}})+Cov_{\mathbf{z}}(q_{1,I_{1}},q_{1,J_{1}}a_{I_{1}}a_{J_{1}})+x\,Cov_{\mathbf{z}}(q_{1,I_{1}},a_{I_{1}})$
$\displaystyle\quad+x\,Cov_{\mathbf{z}}(q_{1,I_{1}},a_{I_{2}}a_{J_{1}})-x\,Cov_{\mathbf{z}}(q_{1,I_{1}},a_{I_{1}}a_{J_{1}})-x\,Cov_{\mathbf{z}}(q_{1,I_{1}},a_{I_{1}}a_{J_{1}})$
$\displaystyle\quad-
Cov_{\mathbf{z}}(q_{1,I_{1}},q_{1,J_{1}}a_{I_{2}}a_{J_{1}})+o(d^{-1})$
$\displaystyle=x^{2}+\frac{1}{Nd}\Big{[}2x\,Cov_{\mathbf{z}}(q_{I_{1}},\mu_{I_{1}})-2x\,Cov_{\mathbf{z}}(q_{1,I_{1}},\sigma_{I_{1},J_{1}})+x\,Cov_{\mathbf{z}}(q_{1,I_{1}},\sigma_{I_{2},I_{3}})$
$\displaystyle\quad\quad\quad\quad\quad-x\,Cov_{\mathbf{z}}(q_{1,I_{1}},\sigma_{I_{1},I_{2}})-Cov_{\mathbf{z}}(q_{1,I_{1}},\mu_{I_{2}})+x\,Cov_{\mathbf{z}}(q_{1,I_{1}},\sigma_{I_{2},J_{1}})$
$\displaystyle\quad\quad\quad\quad\quad+Cov_{\mathbf{z}}(q_{1,I_{1}},q_{1,J_{1}}\sigma_{I_{1},J_{1}})-Cov_{\mathbf{z}}(q_{1,I_{1}},q_{1,J_{1}}\sigma_{I_{2},J_{1}})\Big{]}+o(d^{-1}).$
(118)

Here, $I_{1}$ and $I_{2}$ are randomly chosen offspring produced in the same
deme chosen at random, while $J_{2}$ is an independently chosen offspring
produced in the whole population. The proof of Lemma 5 is completed by
combining Lemmas 1 and 2 with the above equation.

### Appendix F: Proof of Proposition 3

Let each deme at the beginning of generation $0$ be labeled by $(i,d_{i})$
where $i$ is the type of the deme and $d_{i}$ its rank when the demes of the
same type are placed in an arbitrary order, for $i=0,1,\ldots,N$ and
$d_{i}=1,\ldots,dz_{i}$. Thus, the frequency of $S_{1}$ at the beginning of
generation $1$ can be expressed as

$X(1)=\frac{1}{d}\sum_{(i,d_{i})}X_{i,d_{i}}(1),$ (119)

where $X_{i,d_{i}}(1)$ represents the frequency of type $S_{1}$ in deme
$(i,d_{i})$ at the beginning of generation $1$. Note that this random variable
is independent of $X_{j,d_{j}}(1)$ for all $(j,d_{j})\neq(i,d_{i})$ and
satisfies

$E_{\mathbf{z}}\left[X_{i,d_{i}}(1)\Big{|}(a_{k,l})\right]=\tilde{\tilde{x}}_{i}.$
(120)

Defining

$U_{i,d_{i}}=X_{i,d_{i}}(1)-\tilde{\tilde{x}}_{i},$ (121)

we have $\left|U_{i,d_{i}}\right|\leq 1$ and
$E_{\mathbf{z}}\left[U_{i,d_{i}}\right]=0$. Moreover, for the third central
moment of the frequency of type $S_{1}$ at the beginning of generation $1$, we
have

$\displaystyle
E_{\mathbf{z}}\left[\Big{(}X(1)-E_{\mathbf{z}}\left[X(1)\right]\Big{)}^{3}\right]$
$\displaystyle=E_{\mathbf{z}}\Bigg{[}\Bigg{(}\frac{1}{d}\sum_{(i,d_{i})}U_{i,d_{i}}\Bigg{)}^{3}\Bigg{]}$
$\displaystyle=\frac{1}{d^{3}}\Bigg{(}\sum_{(i,d_{i})}E_{\mathbf{z}}\left[U_{i,d_{i}}^{3}\right]+3\sum_{(i,d_{i}),(j,d_{j})}E_{\mathbf{z}}\left[U_{i,d_{i}}^{2}\right]E_{\mathbf{z}}\left[U_{j,d_{j}}\right]$
$\displaystyle\quad\quad\quad+\sum_{\textrm{different
}(i,d_{i}),(j,d_{j}),(k,d_{k})}E_{\mathbf{z}}\left[U_{i,d_{i}}\right]E_{\mathbf{z}}\left[U_{j,d_{j}}\right]E_{\mathbf{z}}\left[U_{k,d_{k}}\right]\Bigg{)}$
$\displaystyle=\frac{1}{d^{3}}\sum_{(i,d_{i})}E_{\mathbf{z}}\left[U_{i,d_{i}}^{3}\right],$
(122)

which is bounded in absolute value by $1/d^{2}$. Similarly, for the fourth
central moment of the frequency of $S_{1}$ at the beginning of generation $1$,
we have

$\displaystyle
E_{\mathbf{z}}\left[\left(X(1)-E_{\mathbf{z}}\left[X(1)\right]\right)^{4}\right]=E_{\mathbf{z}}\Bigg{[}\Bigg{(}\frac{1}{d}\sum_{(i,d_{i})}U_{i,d_{i}}\Bigg{)}^{4}\Bigg{]}$
$\displaystyle=\frac{1}{d^{4}}\Bigg{[}\sum_{(i,d_{i})}E_{\mathbf{z}}\left[U_{i,d_{i}}^{4}\right]+4\sum_{(i,d_{i})\neq(j,d_{j})}E_{\mathbf{z}}\left[U_{i,d_{i}}^{3}\right]E_{\mathbf{z}}\left[U_{j,d_{j}}\right]+6\sum_{(i,d_{i})\neq(j,d_{j})}E_{\mathbf{z}}\left[U_{i,d_{i}}^{2}\right]E_{\mathbf{z}}\left[U_{j,d_{j}}^{2}\right]$
$\displaystyle\quad\quad+12\sum_{\textrm{different
}(i,d_{i}),(j,d_{j}),(k,d_{k})}E_{\mathbf{z}}\left[U_{i,d_{i}}^{2}\right]E_{\mathbf{z}}\left[U_{j,d_{j}}\right]E_{\mathbf{z}}\left[U_{k,d_{k}}\right]$
$\displaystyle\quad\quad+\sum_{\textrm{different
}(i,d_{i}),(j,d_{j}),(k,d_{k}),(l,d_{l})}E_{\mathbf{z}}\left[U_{i,d_{i}}\right]E_{\mathbf{z}}\left[U_{j,d_{j}}\right]E_{\mathbf{z}}\left[U_{k,d_{k}}\right]E_{\mathbf{z}}\left[U_{l,d_{l}}\right]\Bigg{]}$
$\displaystyle=\frac{1}{d^{4}}\Bigg{[}\sum_{(i,d_{i})}E_{\mathbf{z}}\left[U_{i,d_{i}}^{4}\right]+6\sum_{(i,d_{i})\neq(j,d_{j})}E_{\mathbf{z}}\left[U_{i,d_{i}}^{2}\right]E_{\mathbf{z}}\left[U_{j,d_{j}}^{2}\right]\Bigg{]}$
$\displaystyle\leq\frac{1}{d^{4}}\Big{[}d+6d(d-1)\Bigg{]}.$ (123)

On the other hand, we have

$\displaystyle E_{\mathbf{z}}\left[\left(X(1)-X(0)\right)^{4}\right]$
$\displaystyle=E_{\mathbf{z}}\left[\left(X(1)-E_{\mathbf{z}}[X(1)]+E_{\mathbf{z}}[X(1)-X(0)]\right)^{4}\right]$
$\displaystyle=E_{\mathbf{z}}\left[\left(X(1)-E_{\mathbf{z}}[X(1)]\right)^{4}\right]+4E_{\mathbf{z}}\left[X(1)-X(0)\right]E_{\mathbf{z}}\left[\left(X(1)-E_{\mathbf{z}}[X(1)]\right)^{3}\right]$
$\displaystyle\quad+6\left(E_{\mathbf{z}}\left[X(1)-X(0)\right]\right)^{2}E_{\mathbf{z}}\left[\left(X(1)-E_{\mathbf{z}}[X(1)]\right)^{2}\right]$
$\displaystyle\quad+\left(E_{\mathbf{z}}\left[X(1)-X(0)\right]\right)^{4}.$
(124)

This is a function $o(d^{-1})$ owing to Eqs. (20), (Appendix F: Proof of
Proposition 3), (Appendix F: Proof of Proposition 3) and the fact that $X(1)$
is bounded in absolute value by $1$.

### Appendix G: Proof of Proposition 5

Note that we have

$\begin{split}Var_{\mathbf{z}}\left[Y_{i}(1)-Y_{i}(0)\right]=Var_{\mathbf{z}}\left[Z_{i}(1)-v_{i}(X(1))\right]\leq
2Var_{\mathbf{z}}\left[Z_{i}(1)\right]+2Var_{\mathbf{z}}\left[v_{i}(X(1))\right].\end{split}$
(125)

Therefore, in order to establish (47), it suffices to show
$Var_{\mathbf{z}}\left[Z_{i}(1)\right]=o(1)$, which is a direct consequence of
Eq. (17), and that $Var_{\mathbf{z}}\left[v_{i}(X(1))\right]=o(1)$. For this
latter variance, we have

$\displaystyle Var_{\mathbf{z}}\left[v_{i}(X(1))\right]$
$\displaystyle=Var_{\mathbf{z}}\left[v_{i}(x)+\sum_{j=1}^{N}r_{j}\left(X(1)-X(0)\right)^{j}\right]$
$\displaystyle=\sum_{j_{1},j_{2}=1}^{N}r_{j_{1}}r_{j_{2}}Cov_{\mathbf{z}}\left[\left(X(1)-X(0)\right)^{j_{1}},\left(X(1)-X(0)\right)^{j_{2}}\right]$
$\displaystyle=\sum_{j_{1},j_{2}=1}^{N}r_{j_{1}}r_{j_{2}}\Bigg{(}E_{\mathbf{z}}\left[\left(X(1)-X(0)\right)^{j_{1}+j_{2}}\right]$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad-
E_{\mathbf{z}}\left[\left(X(1)-X(0)\right)^{j_{1}}\right]E_{\mathbf{z}}\left[\left(X(1)-X(0)\right)^{j_{2}}\right]\Bigg{)}.$
(126)

This is a function $o(1)$ owing to Propositions 1 and 2 besides the fact that
$X(1)-X(0)$ is bounded in absolute value by $1$.

### Appendix H: Exact identity-by-descent measures

In this Appendix, we deduce the exact identity-by-descent measures defined in
Section 4 for up to five offspring $I_{1},\ldots,I_{5}$ chosen at random in
the same deme at the beginning of a given generation in a neutral population
of an infinite number of demes of fixed size $N$ that is in the stationary
state.

#### H1: identity-by-descent measures for two offspring

Considering the parents of two offspring $I_{1}$ and $I_{2}$ one generation
back, three potential scenarios are possible:

  * •

The offspring have different parents that come from the same deme. This
scenario occurs with probability $\left(1-\frac{1}{N}\right)(1-m)$.

  * •

The offspring have different parents that come from different demes. This
situation has probability $\left(1-\frac{1}{N}\right)m$.

  * •

The offspring have the same parent, a situation that occurs with probability
$\frac{1}{N}$.

Then, the probabilities $f_{2}$ and $f_{11}$ that $I_{1}$ and $I_{2}$ are
identical by descent and not identical by descent, respectively, must satisfy

$\displaystyle f_{2}=\left(1-\frac{1}{N}\right)(1-m)\times
f_{2}+\left(1-\frac{1}{N}\right)m\times 0+\frac{1}{N}\times 1,$ (127a)
$\displaystyle f_{11}=\left(1-\frac{1}{N}\right)(1-m)\times
f_{11}+\left(1-\frac{1}{N}\right)m\times 1+\frac{1}{N}\times 0,$ (127b)

from which

$\displaystyle f_{2}=\frac{1}{mN+1-m},$ (128a) $\displaystyle
f_{11}=\frac{m(N-1)}{mN+1-m}.$ (128b)

#### H2: identity-by-descent measures for three offspring

##### Calculation of $f_{3}$:

For three offspring to be identical by descent, three situations are possible
for their parents one generation back:

  * •

The offspring have three different parents that come from the same deme. This
scenario has probability
$\left(1-\frac{1}{N}\right)\left(1-\frac{2}{N}\right)(1-m)$.

  * •

The offspring have exactly two different parents that come from the same deme,
a situation whose probability is
$\frac{\binom{3}{2}}{N}\left(1-\frac{1}{N}\right)(1-m)$.

  * •

The offspring have the same parent, which occurs with probability
$\frac{1}{N^{2}}$.

Note that all other scenarios are not compatible with $I_{1}$, $I_{2}$ and
$I_{3}$ being identical by descent. By conditioning on the parents one
generation back, we find

$f_{3}=\left(1-\frac{1}{N}\right)\left(1-\frac{2}{N}\right)(1-m)\times
f_{3}+\frac{3}{N}\left(1-\frac{1}{N}\right)(1-m)\times
f_{2}+\frac{1}{N^{2}}\times 1,$ (129)

from which we deduce

$\displaystyle f_{3}$
$\displaystyle=\frac{1+3(N-1)(1-m)f_{2}}{1+(N-1)(3-2m+mN)}$
$\displaystyle=\frac{1+(N-1)(3-2m)}{(mN+1-m)(mN^{2}+(3N-2)(1-m))}.$ (130)

##### Calculation of $f_{21}$:

For an offspring $I_{1}$ to be identical by descent to an offspring $I_{2}$
but not to an offspring $I_{3}$, the scenarios to consider for the parents of
the offspring one generation back and their probabilities are the following:

  * •

The offspring have different parents that come from the same deme, which has
probability $\left(1-\frac{1}{N}\right)\left(1-\frac{2}{N}\right)(1-m)$.

  * •

The offspring $I_{1}$ and $I_{2}$ have the same parent different from the
parent of $I_{3}$ but coming from the same deme, whose probability is
$\frac{1}{N}\left(1-\frac{1}{N}\right)(1-m)$.

  * •

The offspring $I_{1}$ and $I_{2}$ have the same parent different from the
parent of $I_{3}$ and coming from a different deme, which occurs with
probability $\frac{1}{N}\left(1-\frac{1}{N}\right)m$.

By conditioning on these possible scenarios one generation back in time, we
get the following relationship

$f_{21}=\left(1-\frac{1}{N}\right)\left(1-\frac{2}{N}\right)(1-m)\times
f_{21}+\frac{1}{N}\left(1-\frac{1}{N}\right)\Big{[}(1-m)\times f_{11}+m\times
1\Big{]},$ (131)

from which we obtain

$\displaystyle f_{21}$
$\displaystyle=\frac{(N-1)[(1-m)\phi_{I_{1}|I_{3}}+m]}{mN^{2}+(3N-2)(1-m)}$
$\displaystyle=\frac{mN(N-1)}{(mN+1-m)(mN^{2}+(3N-2)(1-m))}.$ (132)

##### Calculation of $f_{111}$:

For no offspring to be identical by descent to another among three, two
scenarios for the parents one generation back are possible:

  * •

The offspring have different that come from the same deme, with probability
$\left(1-\frac{1}{N}\right)\left(1-\frac{2}{N}\right)(1-m)$.

  * •

The offspring have different that come from different deme, with probability
$\left(1-\frac{1}{N}\right)\left(1-\frac{2}{N}\right)m$.

By conditioning on these scenarios, we get

$f_{111}=\left(1-\frac{1}{N}\right)\left(1-\frac{2}{N}\right)(1-m)\times
f_{111}+\left(1-\frac{1}{N}\right)\left(1-\frac{2}{N}\right)m\times 1,$ (133)

from which we obtain

$f_{111}=\frac{(N-1)(N-2)m}{mN^{2}+(3N-2)(1-m)}.$ (134)

#### H3: identity-by-descent measures for four offspring

Similarly, conditioning on the number of parents of four offspring and the
demes where they come from, we obtain the system of equations

$\displaystyle f_{4}=$
$\displaystyle\left(1-\frac{1}{N}\right)\left(1-\frac{2}{N}\right)\left(1-\frac{3}{N}\right)(1-m)f_{4}+\frac{6}{N}\left(1-\frac{1}{N}\right)\left(1-\frac{2}{N}\right)(1-m)f_{3}$
$\displaystyle+\frac{7}{N^{2}}\left(1-\frac{1}{N}\right)(1-m)f_{2}+\frac{1}{N^{3}},$
(135a) $\displaystyle f_{31}=$
$\displaystyle\left(1-\frac{1}{N}\right)\left(1-\frac{2}{N}\right)\left(1-\frac{3}{N}\right)(1-m)f_{31}+\frac{3}{N}\left(1-\frac{1}{N}\right)\left(1-\frac{2}{N}\right)(1-m)f_{21}$
$\displaystyle+\frac{1}{N^{2}}\left(1-\frac{1}{N}\right)\Big{[}(1-m)f_{11}+m\Big{]},$
(135b) $\displaystyle f_{22}=$
$\displaystyle\left(1-\frac{1}{N}\right)\left(1-\frac{2}{N}\right)\left(1-\frac{3}{N}\right)(1-m)f_{22}+\frac{2}{N}\left(1-\frac{1}{N}\right)\left(1-\frac{2}{N}\right)(1-m)f_{21}$
$\displaystyle+\frac{1}{N^{2}}\left(1-\frac{1}{N}\right)\Big{[}(1-m)f_{11}+m\Big{]},$
(135c) $\displaystyle f_{211}=$
$\displaystyle\left(1-\frac{1}{N}\right)\left(1-\frac{2}{N}\right)\left(1-\frac{3}{N}\right)(1-m)f_{211}$
$\displaystyle+\frac{1}{N}\left(1-\frac{1}{N}\right)\left(1-\frac{2}{N}\right)\Big{[}(1-m)f_{111}+m\Big{]},$
(135d) $\displaystyle f_{1111}=$
$\displaystyle\left(1-\frac{1}{N}\right)\left(1-\frac{2}{N}\right)\left(1-\frac{3}{N}\right)\Big{[}(1-m)f_{1111}+m\Big{]}.$
(135e)

The solution of this system is given by

$\displaystyle
f_{4}=\frac{1+7(N-1)(1-m)f_{2}+6(N-1)(N-2)(1-m)f_{3}}{mN^{3}+(6N^{2}-11N+6)(1-m)},$
(136a) $\displaystyle
f_{31}=\frac{(N-1)m+(N-1)(1-m)f_{11}+3(N-1)(N-2)(1-m)f_{21}}{mN^{3}+(6N^{2}-11N+6)(1-m)},$
(136b) $\displaystyle
f_{22}=\frac{(N-1)m+(N-1)(1-m)f_{11}+2(N-1)(N-2)(1-m)f_{21}}{mN^{3}+(6N^{2}-11N+6)(1-m)},$
(136c) $\displaystyle
f_{2111}=\frac{(N-1)(N-2)m+(N-1)(N-2)(1-m)f_{111}}{mN^{3}+(6N^{2}-11N+6)(1-m)},$
(136d) $\displaystyle
f_{1111}=\frac{(N-1)(N-2)(N-3)m}{mN^{3}+(6N^{2}-11N+6)(1-m)}.$ (136e)

#### H4: identity-by-descent measures for five offspring

Proceeding as previously, the identity-by-descent measures for five offspring
are found to be

$\displaystyle
f_{5}=\frac{1+(N-1)(1-m)[15f_{2}+25(N-2)f_{3}+10(N-2)(N-3)f_{4}]}{mN^{4}+(1-m)(10N^{3}-35N^{2}+50N-24)},$
(137a) $\displaystyle
f_{41}=\frac{m(N-1)+(N-1)(1-m)[f_{11}+7(N-2)f_{21}+6(N-2)(N-3)f_{31}]}{mN^{4}+(1-m)(10N^{3}-35N^{2}+50N-24)},$
(137b) $\displaystyle
f_{32}=\frac{m(N-1)+(N-1)(1-m)[f_{11}+4(N-2)f_{21}+(N-2)(N-3)(f_{31}+3f_{22})]}{mN^{4}+(1-m)(10N^{3}-35N^{2}+50N-24)},$
(137c) $\displaystyle
f_{311}=\frac{(N-1)(N-2)m+(N-1)(N-2)(1-m)[f_{111}+3(N-3)f_{211}]}{mN^{4}+(1-m)(10N^{3}-35N^{2}+50N-24)},$
(137d) $\displaystyle
f_{221}=\frac{(N-1)(N-2)m+(N-1)(N-2)(1-m)[f_{111}+2(N-3)f_{211}]}{mN^{4}+(1-m)(10N^{3}-35N^{2}+50N-24)},$
(137e) $\displaystyle
f_{2111}=\frac{(N-1)(N-2)(N-3)[m+(1-m)f_{1111}]}{mN^{4}+(1-m)(10N^{3}-35N^{2}+50N-24)},$
(137f) $\displaystyle
f_{11111}=\frac{(N-1)(N-2)(N-3)(N-4)m}{mN^{4}+(1-m)(10N^{3}-35N^{2}+50N-24)}.$
(137g)

### Appendix I: Proof of Proposition 7

In order to derive explicit expressions for $M(\mathbf{v}(x))$ and
$Q(\mathbf{v}(x))$ in Theorem 1, we will write the covariances in the
expressions of $M(\mathbf{z})$ and $Q(\mathbf{z})$ given in Propositions 1 and
2 in terms of the identity-by-descent measures defined in Eq. (52) for
offspring chosen at random in the same deme in an infinite neutral population
in the stationary state $\mathbf{z}=\mathbf{v}(x)$ defined by Eq. (38). Note
that an offspring chosen at random and all offspring identical by descent to
this offspring are of type $S_{1}$ with probability $x$ and of type $S_{2}$
with probability $1-x$ independently of all offspring that are not identical
by descent to this offspring, where $x$ is the constant frequency of $S_{1}$
in the neutral population at equilibrium. Moreover, offspring in different
demes are necessarily non identical by descent since they cannot have a common
ancestor in the same deme in the case of an infinite number of demes.

By conditioning on the identity-by-descent status of two offspring $I_{1}$ and
$I_{2}$ randomly chosen in the same deme, we have

$\displaystyle
Cov_{\mathbf{v}(x)}\left(q_{1,I_{1}},q_{k,I_{1}}q_{l,I_{2}}\right)$
$\displaystyle=xf_{I_{1}|I_{2}}\left(\delta_{1k}x_{(l)}-x_{(k)}x_{(l)}\right)+xf_{I_{1}I_{2}}\left(\delta_{1kl}-\delta_{kl}x_{(k)}\right)$
$\displaystyle=x\left(\delta_{1k}-x_{(k)}\right)\left(\delta_{kl}f_{2}+f_{11}x_{(l)}\right),$
(138)

where $x_{(k)}=x$ if $k=1$ and $1-x$ if $k=2$, while $\delta_{l_{1}l_{2}\ldots
l_{n}}=1$ if $l_{1}=\cdots=l_{n}$ and $0$ otherwise.

Similarly, by conditioning on the identity-by-descent relationship between
three offspring $I_{1}$, $I_{2}$ and $I_{3}$ randomly selected in the same
deme, as well as between two offspring $J_{1}$ and $J_{2}$ independently
chosen in a same deme, we get

$\displaystyle
Cov_{\mathbf{v}(x)}\left(q_{1,I_{1}},q_{k,I_{2}}q_{l,I_{3}}\right)$
$\displaystyle\quad\quad\quad\quad=xf_{I_{1}I_{2}I_{3}}\left(\delta_{1kl}-\delta_{kl}x_{(k)}\right)+xf_{I_{1}I_{2}|I_{3}}\left(\delta_{1k}x_{(l)}-x_{(k)}x_{(l)}\right)+xf_{I_{1}I_{3}|I_{2}}\left(\delta_{1l}x_{(k)}-x_{(k)}x_{(l)}\right)$
$\displaystyle\quad\quad\quad\quad=x\left(\delta_{1k}-x_{(k)}\right)\left(f_{3}\delta_{kl}+f_{21}x_{(l)}\right)+xx_{(k)}f_{21}\left(\delta_{1l}-x_{(l)}\right),$
(139) $\displaystyle
Cov_{\mathbf{v}(x)}\left(q_{1,I_{1}},q_{k,I_{1}}q_{l,I_{2}}q_{k,J_{1}}q_{l,J_{2}}\right)$
$\displaystyle\quad\quad\quad\quad=Cov_{\mathbf{v}(x)}\left(q_{1,I_{1}},q_{k,I_{1}}q_{l,I_{2}}\right)E_{\mathbf{v}(x)}\left(q_{k,J_{1}}q_{l,J_{2}}\right)$
$\displaystyle\quad\quad\quad\quad=x\left(\delta_{1k}-x_{(k)}\right)\Big{[}f_{I_{1}I_{2}}\delta_{kl}+f_{I_{1}|I_{2}}x_{(l)}\Big{]}\Big{[}f_{J_{1}J_{2}}\delta_{kl}x_{(k)}+f_{J_{1}|J_{2}}x_{(k)}x_{(l)}\Big{]}$
$\displaystyle\quad\quad\quad\quad=x\left(\delta_{1k}-x_{(k)}\right)x_{(k)}\Big{[}f_{2}\delta_{kl}+f_{11}x_{(l)}\Big{]}^{2},$
(140) $\displaystyle
Cov_{\mathbf{v}(x)}\left(q_{1,I_{1}},q_{k,I_{2}}q_{l,I_{3}}q_{k,J_{1}}q_{l,J_{2}}\right)$
$\displaystyle\quad\quad\quad\quad=Cov_{\mathbf{v}(x)}\left(q_{1,I_{1}},q_{k,I_{2}}q_{l,I_{3}}\right)E_{\mathbf{v}(x)}\left(q_{k,J_{1}}q_{l,J_{2}}\right)$
$\displaystyle\quad\quad\quad\quad=xx_{(k)}\Big{[}\left(\delta_{1k}-x_{(k)}\right)\left(f_{3}\delta_{kl}+f_{21}x_{(l)}\right)+x_{(k)}f_{21}\left(\delta_{1l}-x_{(l)}\right)\Big{]}\Big{[}f_{2}\delta_{kl}+f_{11}x_{(l)}\Big{]},$
(141) $\displaystyle
Cov_{\mathbf{v}(x)}\left(q_{1,I_{1}}q_{1,J_{1}},q_{k,I_{1}}q_{l,I_{2}}q_{k,J_{1}}q_{l,J_{2}}\right)$
$\displaystyle\quad\quad\quad\quad=E_{\mathbf{v}(x)}\left(q_{1,I_{1}}q_{k,I_{1}}q_{l,I_{2}}\right)^{2}-x^{2}E_{\mathbf{v}(x)}\left(q_{k,I_{1}}q_{l,I_{2}}\right)^{2}$
$\displaystyle\quad\quad\quad\quad=x^{2}\Big{[}f_{11}\delta_{1k}x_{(l)}+f_{2}\delta_{1kl}\Big{]}^{2}-x^{2}x_{(k)}^{2}\Big{[}f_{11}x_{(l)}+f_{2}\delta_{kl}\Big{]}^{2}$
$\displaystyle\quad\quad\quad\quad=x^{2}\left(\delta_{1k}-x_{(k)}^{2}\right)\Big{[}f_{2}\delta_{kl}+f_{11}x_{(l)}\Big{]}^{2},$
(142) $\displaystyle
Cov_{\mathbf{v}(x)}\left(q_{1,I_{1}},q_{1,J_{1}}q_{k,I_{1}}q_{l,I_{2}}q_{k,J_{1}}q_{l,J_{2}}\right)$
$\displaystyle\quad\quad\quad\quad=Cov_{\mathbf{v}(x)}\left(q_{1,I_{1}},q_{k,I_{1}}q_{l,I_{2}}\right)E_{\mathbf{v}(x)}\left(q_{1,J_{1}}q_{k,J_{1}}q_{l,J_{2}}\right)$
$\displaystyle\quad\quad\quad\quad=x\left(\delta_{1k}-x_{(k)}\right)\Big{[}f_{I_{1}I_{2}}\delta_{kl}+f_{I_{1}|I_{2}}x_{(l)}\Big{]}\Big{[}xx_{(l)}f_{I_{1}|I_{2}}\delta_{1k}+xf_{I_{1}I_{2}}\delta_{1kl}\Big{]},$
$\displaystyle\quad\quad\quad\quad=x^{2}\delta_{1k}\left(\delta_{1k}-x_{(k)}\right)\Big{[}f_{2}\delta_{kl}+f_{11}x_{(l)}\Big{]}^{2},$
(143) $\displaystyle
Cov_{\mathbf{v}(x)}\left(q_{1,I_{1}},q_{1,J_{1}}q_{k,I_{2}}q_{l,I_{3}}q_{k,J_{1}}q_{l,J_{2}}\right)$
$\displaystyle\quad\quad\quad\quad=Cov_{\mathbf{v}(x)}\left(q_{1,I_{1}},q_{k,I_{2}}q_{l,I_{3}}\right)E_{\mathbf{v}(x)}\left(q_{1,J_{1}}q_{k,J_{1}}q_{l,J_{2}}\right)$
$\displaystyle\quad\quad\quad\quad=x^{2}\delta_{1k}\Big{[}\left(\delta_{1k}-x_{(k)}\right)\left(f_{3}\delta_{kl}+f_{21}x_{(l)}\right)+x_{(k)}f_{21}\left(\delta_{1l}-x_{(l)}\right)\Big{]}\Big{[}f_{2}\delta_{kl}+f_{11}x_{(l)}\Big{]},$
(144) $\displaystyle
Cov_{\mathbf{v}(x)}\left(q_{1,I_{1}},q_{1,J_{1}}q_{k,I_{1}}q_{l,I_{2}}q_{k,J_{2}}q_{l,J_{3}}\right)$
$\displaystyle\quad\quad\quad\quad=Cov_{\mathbf{v}(x)}\left(q_{1,I_{1}},q_{k,I_{1}}q_{l,I_{2}}\right)E_{\mathbf{v}(x)}\left(q_{1,J_{1}}q_{k,J_{2}}q_{l,J_{3}}\right)$
$\displaystyle\quad\quad\quad\quad=x^{2}\left(\delta_{1k}-x_{(k)}\right)\left(f_{2}\delta_{kl}+f_{11}x_{(l)}\right)$
$\displaystyle\quad\quad\quad\quad\quad\times\Big{[}\delta_{1k}\left(f_{3}\delta_{kl}+f_{21}x_{(l)}\right)+\left(f_{21}\delta_{1l}+f_{21}\delta_{kl}+f_{111}x_{(l)}\right)x_{(k)}\Big{]},$
(145) $\displaystyle
Cov_{\mathbf{v}(x)}\left(q_{1,I_{1}},q_{1,J_{1}}q_{k,I_{2}}q_{l,I_{3}}q_{k,J_{2}}q_{l,J_{3}}\right)$
$\displaystyle\quad\quad\quad\quad=Cov_{\mathbf{v}(x)}\left(q_{1,I_{1}},q_{k,I_{2}}q_{l,I_{3}}\right)E_{\mathbf{v}(x)}\left(q_{1,J_{1}}q_{k,J_{2}}q_{l,J_{3}}\right)$
$\displaystyle\quad\quad\quad\quad=x^{2}\Big{[}\left(\delta_{1k}-x_{(k)}\right)\left(f_{3}\delta_{kl}+f_{21}x_{(l)}\right)+x_{(k)}f_{21}\left(\delta_{1l}-x_{(l)}\right)\Big{]}$
$\displaystyle\quad\quad\quad\quad\quad\times\Big{[}\delta_{1k}\left(f_{3}\delta_{kl}+f_{21}x_{(l)}\right)+\left(f_{21}\delta_{1l}+f_{21}\delta_{kl}+f_{111}x_{(l)}\right)x_{(k)}\Big{]}.$
(146)

Finally, by conditioning on the identity-by-descent relationship between four
or five individuals randomly selected from the same deme, we obtain

$\displaystyle
Cov_{\mathbf{v}(x)}\left(q_{1,I_{1}},q_{k,I_{1}}q_{l,I_{2}}q_{k,I_{3}}q_{l,I_{4}}\right)$
$\displaystyle\quad\quad\quad\quad=x\left(\delta_{1k}-x_{(k)}\right)\Big{[}\delta_{kl}f_{4}+x_{(l)}f_{22}+x^{2}_{(l)}f_{211}\Big{]}$
$\displaystyle\quad\quad\quad\quad\quad+x\left(\delta_{1k}x-x^{2}_{(k)}\right)\Big{[}\delta_{kl}\left(2f_{22}+4f_{31}\right)+x_{(l)}f_{211}+x^{2}_{(l)}f_{1111}\Big{]}$
$\displaystyle\quad\quad\quad\quad\quad+4x\delta_{kl}\left(\delta_{1k}x^{2}-x^{3}_{(k)}\right)f_{211}$
(147)

and

$\displaystyle
Cov_{\mathbf{v}(x)}\left(q_{1,I_{1}},q_{k,I_{2}}q_{k,I_{3}}q_{l,I_{4}}q_{l,I_{5}}\right)$
$\displaystyle\quad\quad\quad\quad=x\left(\delta_{1k}-x_{(k)}\right)\left[\delta_{kl}f_{5}+x_{(l)}f_{32}+x_{(l)}^{2}f_{311}\right]+4x\delta_{kl}\left(\delta_{1k}x^{2}-x_{(k)}^{3}\right)\Big{[}2f_{221}+f_{311}\Big{]}$
$\displaystyle\quad\quad\quad\quad\quad+2x\left(\delta_{1k}x-x_{(k)}^{2}\right)\Big{[}\delta_{kl}\left(2f_{41}+4f_{32}\right)+x_{(l)}f_{221}+x_{(l)}^{2}f_{2111}\Big{]}$
$\displaystyle\quad\quad\quad\quad\quad+2xx_{(k)}\left(\delta_{1l}x-x_{(l)}^{2}\right)\Big{[}f_{221}+x_{(k)}f_{2111}\Big{]}$
$\displaystyle\quad\quad\quad\quad\quad+xx_{(k)}\left(\delta_{1l}-x_{(l)}\right)\Big{[}f_{32}+x_{(k)}f_{311}\Big{]}.$
(148)

Using the above expressions for the covariances and the relationships

$\displaystyle f_{2}=f_{3}+f_{21},$ (149a) $\displaystyle
f_{11}=2f_{21}+f_{111},$ (149b) $\displaystyle f_{4}=f_{5}+f_{41},$ (149c)
$\displaystyle f_{31}=f_{41}+f_{32}+f_{311},$ (149d) $\displaystyle
f_{22}=2f_{32}+f_{221},$ (149e) $\displaystyle
f_{211}=f_{311}+2f_{221}+f_{2111},$ (149f)

the infinitesimal mean $M(\mathbf{v}(x))$ in Theorem 1, given by

$\displaystyle\sum_{k,l=1}^{2}\mu_{k,l}\Big{[}Cov_{\mathbf{v}(x)}\left(q_{1,I_{1}},q_{k,I_{1}}q_{l,I_{2}}\right)-(1-m)\,Cov_{\mathbf{v}(x)}\left(q_{1,I_{1}},q_{k,I_{2}}q_{l,I_{3}}\right)\Big{]}$
$\displaystyle+\sum_{k,l=1}^{2}\sigma^{2}_{k,l}\Big{[}(1-m)\Big{(}Cov_{\mathbf{v}(x)}\left(q_{1,I_{1}},q_{k,I_{1}}q_{k,I_{2}}q_{l,I_{3}}q_{l,I_{4}}\right)$
$\displaystyle\quad\quad\quad\quad\quad-
Cov_{\mathbf{v}(x)}\left(q_{1,I_{1}},q_{k,I_{2}}q_{k,I_{3}}q_{l,I_{4}}q_{l,I_{5}}\right)\Big{)}-m\,Cov_{\mathbf{v}(x)}\left(q_{1,I_{1}},q_{k,I_{1}}q_{l,I_{2}}q_{k,J_{1}}q_{l,J_{2}}\right)\Big{]},$
(150)

can be expressed as

$\displaystyle\sum_{k,l=1}^{2}\mu_{k,l}\Bigg{[}x\left(\delta_{1k}-x_{(k)}\right)\left(\delta_{kl}f_{2}+f_{11}x_{(l)}\right)-(1-m)\Big{(}x\left(\delta_{1k}-x_{(k)}\right)\left(f_{3}\delta_{kl}+f_{21}x_{(l)}\right)$
$\displaystyle+xx_{(k)}f_{21}\left(\delta_{1l}-x_{(l)}\right)\Big{)}\Bigg{]}+\sum_{k,l=1}^{2}\sigma^{2}_{k,l}\Bigg{[}(1-m)\Big{(}x\left(\delta_{1k}-x_{(k)}\right)\Big{[}\delta_{kl}f_{4}+x_{(l)}f_{22}+x^{2}_{(l)}f_{211}\Big{]}$
$\displaystyle+x(\delta_{1k}x-x_{(k)}^{2})\Big{[}\delta_{kl}\left(2f_{22}+4f_{31}\right)+x_{(l)}f_{211}+x^{2}_{(l)}f_{1111}\Big{]}+4x\delta_{kl}(\delta_{1k}x^{2}-x_{(k)}^{3})f_{211}\Big{)}$
$\displaystyle-(1-m)\Big{(}x\left(\delta_{1k}-x_{(k)}\right)\Big{[}\delta_{kl}f_{5}+x_{(l)}f_{32}+x_{(l)}^{2}f_{311}\Big{]}+4x\delta_{kl}(\delta_{1k}x^{2}-x_{(k)}^{3})[2f_{221}+f_{311}]$
$\displaystyle+2x(\delta_{1k}x-x_{(k)}^{2})\Big{[}\delta_{kl}\left(2f_{41}+4f_{32}\right)+x_{(l)}f_{221}+x_{(l)}^{2}f_{2111}\Big{]}+xx_{(k)}(\delta_{1l}-x_{(l)})[f_{32}+x_{(k)}f_{311}]$
$\displaystyle+2xx_{(k)}(\delta_{1l}x-x_{(l)}^{2})[f_{221}+x_{(k)}f_{2111}]\Big{)}-mx\left(\delta_{1k}-x_{(k)}\right)x_{(k)}\left[f_{2}\delta_{kl}+f_{11}x_{(l)}\right]^{2}\Bigg{]}.$
(151)

As for the infinitesimal variance $Q(\mathbf{v}(x))$ in Theorem 1, given by

$\displaystyle x(1-x)\left(1+(1-m)(Nm-1)f_{2}\right)$
$\displaystyle+\sum_{k,l=1}^{2}\sigma^{2}_{k,l}\Bigg{[}m^{2}Cov_{\mathbf{v}(x)}\left(q_{1,I_{1}}q_{1,J_{1}},q_{k,I_{1}}q_{l,I_{2}}q_{k,J_{1}}q_{l,J_{2}}\right)-2mx\,Cov_{\mathbf{v}(x)}\left(q_{1,I_{1}},q_{k,I_{1}}q_{l,I_{2}}q_{k,J_{1}}q_{l,J_{2}}\right)$
$\displaystyle+(1-m^{2})\Big{(}Cov_{\mathbf{v}(x)}\left(q_{1,I_{1}},q_{1,J_{1}}q_{k,I_{1}}q_{l,I_{2}}q_{k,J_{1}}q_{l,J_{2}}\right)-Cov_{\mathbf{v}(x)}\left(q_{1,I_{1}},q_{1,J_{1}}q_{k,I_{2}}q_{l,I_{3}}q_{k,J_{1}}q_{l,J_{2}}\right)\Big{)}$
$\displaystyle+2m(1-m)\,Cov_{\mathbf{v}(x)}\left(q_{1,I_{1}},q_{k,I_{2}}q_{l,I_{3}}q_{k,J_{1}}q_{l,J_{2}}\right)+(1-m)^{2}\Big{(}Cov_{\mathbf{v}(x)}\left(q_{1,I_{1}},q_{1,J_{1}}q_{k,I_{2}}q_{l,I_{3}}q_{k,J_{2}}q_{l,J_{3}}\right)$
$\displaystyle-
Cov_{\mathbf{v}(x)}\left(q_{1,I_{1}},q_{1,J_{1}}q_{k,I_{1}}q_{l,I_{2}}q_{k,J_{2}}q_{l,J_{3}}\right)\Big{)}\Bigg{]},$
(152)

it takes the form

$\displaystyle x(1-x)\left(1+(1-m)(Nm-1)f_{2}\right)$
$\displaystyle+\sum_{k,l=1}^{2}\sigma^{2}_{k,l}\Bigg{[}m^{2}x^{2}(\delta_{1k}-x_{(k)}^{2})\left(f_{2}\delta_{kl}+f_{11}x_{(l)}\right)^{2}-2mx^{2}\left(\delta_{1k}-x_{(k)}\right)x_{(k)}\left(f_{2}\delta_{kl}+f_{11}x_{(l)}\right)^{2}$
$\displaystyle+\delta_{1k}x^{2}(1-m^{2})(f_{2}\delta_{kl}+f_{11}x_{(l)})\Big{(}\left(\delta_{1k}-x_{(k)}\right)(f_{21}\delta_{kl}+(f_{21}+f_{111})x_{(l)})-x_{(k)}f_{21}\left(\delta_{1l}-x_{(l)}\right)\Big{)}$
$\displaystyle+2m(1-m)xx_{(k)}(f_{2}\delta_{kl}+f_{11}x_{(l)})\Big{(}\left(\delta_{1k}-x_{(k)}\right)\left(f_{3}\delta_{kl}+f_{21}x_{(l)}\right)+x_{(k)}f_{21}\left(\delta_{1l}-x_{(l)}\right)\Big{)}$
$\displaystyle+x^{2}(1-m)^{2}\Big{(}\delta_{1k}\left(f_{3}\delta_{kl}+f_{21}x_{(l)}\right)+\left(f_{21}\delta_{1l}+f_{21}\delta_{kl}+f_{111}x_{(l)}\right)x_{(k)}\Big{)}$
$\displaystyle\times\Big{(}x_{(k)}f_{21}(\delta_{1l}-x_{(l)})-(\delta_{1k}-x_{(k)})(f_{21}\delta_{kl}+(f_{21}+f_{111})x_{(l)})\Big{)}\Bigg{]}.$
(153)

Here, besides the previous expressions for the covariances, we have used the
identities

$\displaystyle Var_{\mathbf{v}(x)}(q_{I_{1}})=x(1-x),$ (154a) $\displaystyle
Cov_{\mathbf{v}(x)}(q_{1,I_{1}},q_{1,I_{2}})=x(1-x)f_{2}.$ (154b)

Finally, letting $N\rightarrow\infty$ and $m\rightarrow 0$ in a way such that
$Nm\rightarrow\nu$, we get the approximations

$\displaystyle M(\mathbf{v}(x))$ $\displaystyle\approx
x(1-x)\Bigg{[}\Big{(}f_{21}+xf_{111}\Big{)}\left(\mu_{1,1}-\mu_{2,1}\right)+\Big{(}f_{21}+(1-x)f_{111}\Big{)}\left(\mu_{1,2}-\mu_{2,2}\right)$
$\displaystyle\quad-\Big{(}f_{41}+(4f_{311}+3f_{221})x+6f_{2111}x^{2}+f_{11111}x^{3}\Big{)}\sigma^{2}_{1,1}$
$\displaystyle\quad-\Big{(}f_{32}+\left(3f_{221}+f_{2111}\right)(1-x)+xf_{311}+2x(1-x)f_{2111}+x(1-x)^{2}f_{11111}\Big{)}\sigma^{2}_{1,2}$
$\displaystyle\quad+\Big{(}f_{32}+\left(3f_{221}+f_{2111}\right)x+f_{311}(1-x)+2x(1-x)f_{2111}+x^{2}(1-x)f_{11111}\Big{)}\sigma^{2}_{2,1}$
$\displaystyle\quad+\Big{(}f_{41}+(4f_{311}+3f_{221})(1-x)+6f_{2111}(1-x)^{2}+f_{11111}(1-x)^{3}\Big{)}\sigma^{2}_{2,2}\Bigg{]}$
(155)

and

$\displaystyle Q(\mathbf{v}(x))$ $\displaystyle\approx x(1-x)(f_{11}+\nu
f_{2})+x^{2}(1-x)^{2}\left(f_{21}+xf_{111}\right)^{2}\left(\sigma^{2}_{1,1}+\sigma^{2}_{2,1}\right)$
$\displaystyle\quad+x^{2}(1-x)^{2}\left(f_{21}+(1-x)f_{111}\right)^{2}\left(\sigma^{2}_{1,2}+\sigma^{2}_{2,2}\right).$
(156)

### Appendix J: Proof of Proposition 8

In a neutral population of an infinite number of demes of size $N$ with $N$
generations as unit of time as $N$ approaches infinity and $m$ approaches zero
such that $Nm\rightarrow\nu$, each pair of lineages of offspring within the
same deme coalesces backward in time at the rate $1$ independently of all the
other pairs of lineages like in the Kingmann coalescent (Kingmann [20]). On
the other hand, each deme goes extinct backward in time at the rate $\nu$.
When this occurs, the lineages within the deme migrate to all different demes
and stay in different demes afterward so that they cannot coalesce anymore.
Therefore, offspring are identical by descent only if they are in the same
deme and their lineages coalesce backward in time before the deme goes extinct
(see, e.g., Hudson [16] and Rousset [49]).

We are interested in the probability for an ordered sample of offspring
selected at random from the same deme, represented by
$I^{(1)}_{1},\ldots,I^{(1)}_{n_{1}},\ldots,I^{(k)}_{1},\ldots,I^{(k)}_{n_{k}}$,
to be such that :

  * •

$I^{(i)}_{1},\ldots,I^{(i)}_{n_{i}}$ are identical by descent for
$i=1,\ldots,k$, and

  * •

$I^{(i)}_{1}$ and $I^{(j)}_{1}$ are not identical by descent if $i\neq j$ for
$i,j=1,\ldots,k$.

This event is represented by
$I^{(1)}_{1}\ldots\,I^{(1)}_{n_{1}}|\ldots\,|I^{(k)}_{1}\ldots\,I^{(k)}_{n_{k}}$
and its probability by
$f_{I^{(1)}_{1}\ldots\,I^{(1)}_{n_{1}}|\ldots\,|I^{(k)}_{1}\ldots\,I^{(k)}_{n_{k}}}$.

The offspring $I_{1},\ldots,I_{n}$ will be identical by descent if and only if
all their lineages coalesce before the deme goes extinct. Note that the
probability of a coalescence event involving two lineages among $j$ in the
same deme before an extinction event involving the deme occurs is given by the
rate of coalescence, $\binom{j}{2}$, divided by the rate of coalescence or
extinction, $\binom{j}{2}+\nu$, for $j=2,\ldots,n$. Therefore, we have

$f_{I_{1}\ldots\,I_{n}}=\prod_{j=2}^{n}\frac{\binom{j}{2}}{\binom{j}{2}+\nu}.$
(157)

On the other hand, none of the offspring $I_{1},\ldots,I_{n}$ will be
identical by descent to another if and only if the first event to occur is the
extinction of the deme, which gives

$f_{I_{1}|\ldots\,|I_{n}}=\frac{\nu}{\binom{n}{2}+\nu}.$ (158)

More generally, the occurrence of the event
$I^{(1)}_{1}\ldots\,I^{(1)}_{n_{1}}|\ldots\,|I^{(k)}_{1}\cdots
I^{(k)}_{n_{k}}$ can be explained by a series of $n_{1}+\cdots+n_{k}-k$
coalescent events followed by an extinction event. Specifically, $n_{i}-1$
coalescence events must occur among the ancestral lines of
$I^{(i)}_{1},\ldots,I^{(i)}_{n_{i}}$ until they reach their most recent common
ancestor, denoted by $J^{(i)}$, for $i=1,\ldots,k$, before the deme goes
extinct. Moreover, the ancestral lines of $J^{(i_{1})}$ and $J^{(i_{2})}$ must
never experience a coalescence event, for $i_{1}\not=i_{2}$, before the deme
goes extinct. As a result, we obtain

$f_{I^{(1)}_{1}\ldots\,I^{(1)}_{n_{1}}|\ldots\,|I^{(k)}_{1}\ldots\,I^{(k)}_{n_{k}}}=\binom{n_{1}+\cdots+n_{k}-k}{n_{1}-1,\ldots,n_{k}-1}\frac{\prod_{i=1}^{k}\prod_{j=2}^{n_{i}}\binom{j}{2}}{\prod_{j=k+1}^{n_{1}+\cdots+n_{k}}\left(\binom{j}{2}+\nu\right)}\times\frac{\nu}{\binom{k}{2}+\nu},$
(159)

where

$\binom{n_{1}+\cdots+n_{k}-k}{n_{1}-1,\ldots,n_{k}-1}=\dfrac{(n_{1}+\cdots+n_{k}-k)!}{(n_{1}-1)!\cdots(n_{k}-1)!},$
(160)

for $n_{1},\ldots,n_{k}\geq 1$ for $k\geq 2$.

## References

  * [1] Antal, T., Nowak, M.A., Traulsen, A. Strategy abundance in games for arbitrary mutation rates. Journal of Theoretical Biology 257, 340-344, 2009.
  * [2] Avery, P.J. The effect of random selection coefficients on populations of finite size - some particular models. Genetics Research 29, 97-112, 1977.
  * [3] Axelrod, R., Hamilton, W.D. The evolution of cooperation. Science 211, 1390-1396, 1981.
  * [4] Cannings, C. The latent roots of certain Markov chains arising in genetics: a new approach. Advances in Applied Probability 6, 260-290, 1974.
  * [5] Cherry, J.L. Selection in a subdivided population with dominance or local frequency dependence. Genetics 163, 1511-1518, 2003a.
  * [6] Cherry, J.L. Selection in a subdivided population with local extinction and recolonization. Genetics 164, 789-795,, 2003b.
  * [7] Cherry, J.L., Wakeley, J. A diffusion approximation for selection and drift in a subdivided population.Genetics 163, 421-428, 2003.
  * [8] Christiansen, F.B. Hard and soft selection in a subdivided population. American Naturalist 109, 11-16, 1975.
  * [9] Cressman, R., Kr̆ivan, V. Using chemical reaction network theory to show stability of distributional dynamics in game theory. Journal of Dynamics and Games 9, 351-371, 2022.
  * [10] Ethier, S. A class of degenerate diffusion processes occurring in population genetics. Communications on Pure and Applied Mathematics 29, 483-493, 1976.
  * [11] Ethier, S.N., Nagylaki, T. Diffusion approximations of Markov chains with two timescales and applications to population genetics. Advances in Applied Probability 12, 14-49, 1980.
  * [12] Ewens, W.J. Mathematical Population Genetics: I Theoretical Introduction. Springer, New York, 2004.
  * [13] Fisher, R.A. The Genetical Theory of Natural Selection. Clarendon Press, Oxford, 1930.
  * [14] Hilbe, C. Local replicator dynamics: A simple link between deterministic and stochastic models of evolutionary game theory. Bulletin of Mathematical Biology 73, 2068-2087, 2011.
  * [15] Hofbauer, J., Sigmund, K. The Theory of Evolution and Dynamical Systems; Cambridge University Press: Cambridge, 1988.
  * [16] Hudson, R.R. Island models and the coalescent process. Molecular Ecology 7, 413-418, 1998.
  * [17] Imhof, L.A., Nowak, M.A. Evolutionary game dynamics in a Wright-Fisher process. Journal of Mathematical Biology 52, 667-681, 2006.
  * [18] Kaplan, H., Hill, K., Hurtado, A.M. Risk, foraging and food sharing among the Ache. In: Cashdan E (ed), _Risk and uncertainty in tribal and peasant economies_. Westview Press, Boulder, pp. 107-144, 1990.
  * [19] Karlin, S., Taylor, P. A Second Course in Stochastic Processes, Academic Press, New York, 1981.
  * [20] Kingman, J.F.C. The coalescent. Stochastic Processes and their Applications 13, 235-248, 1982.
  * [21] Kr̆ivan, V., Cressman, R. Defectors’ intolerance of others promotes cooperation in the repeated public goods game with opting out. Scientific Reports 10, 19511, 2020.
  * [22] Kroumi, D., Lessard, S. Conditions for cooperation to be more abundant than defection in a hierarchically structured population. Dynamic Games and its Applications  5, 239-262, 2015.
  * [23] Kroumi, D., Lessard, S. Strong migration limit for games in structured populations: Applications to dominance hierarchy and set structure. Games 3, 318-346, 2015.
  * [24] Kroumi, D., Lessard S. The effect of variability in payoffs on average abundance in two-player linear games under symmetric mutation. Journal of Theoretical Biology 513,110569, 2021.
  * [25] Kroumi, D., Lessard, S. Average abundancy of cooperation in multi-player games with random payoffs. Journal of Mathematical Biology 85, 27, 2022.
  * [26] Kroumi, D., Martin, É., Lessard, S. Evolution of cooperation with respect to fixation probabilities in multi-player games with random payoffs. Theoretical Population Biology 145, 1-21, 2022.
  * [27] Kroumi, D., Martin, É., Li, C., Lessard, S. Effect of Variability in Payoffs on Conditions for the Evolution of Cooperation in a Small Population. Dynamic Games and its Applications  11, 803-834, 2021.
  * [28] Ladret, V., Lessard, S. Evolutionary game dynamics in a finite asymmetric two-deme population and emergence of cooperation. Journal of Theoretical Biology 225, 137-151, 2008.
  * [29] Lambert, A. Probability of fixation under weak selection: a branching process unifying approach. Theoretical Population Biology 69, 419-441, 2006.
  * [30] Lande, R., Engen, S., Saether, B.E. Stochastic Population Dynamics in Ecology and Conservation. Oxford University Press, Oxford, 2003
  * [31] Lessard, S. An exact sampling formula for the Wright–Fisher model and a conjecture about the finite-island model. Genetics 177, 1249-1254, 2007.
  * [32] Lessard, S. Diffusion approximations for one-locus multi-allele kin selection, mutation and random drift in group-structured populations: a unifying approach to selection models in population genetics. Journal of Mathematical Biology 59, 659-696, 2009.
  * [33] Lessard, S., Ladret, V. The probability of fixation of a single mutant in an exchangeable selection model. Journal of Mathematical Biology 54, 721-744, 2007.
  * [34] Li, C., Lessard, S. Randomized matrix game in a finite population: effect of stochastic fluctuations in the payoffs on the evolution of cooperation. Theoretical Population Biology 134, 77-91, 2020.
  * [35] Li, C., Lessard, S. The effect of the opting-out strategy on conditions for selection to favor the evolution of cooperation in a finite population. Journal of Theoretical Biology 510, 110543, 2021.
  * [36] Maynard Smith, J. Evolution and the Theory of games. Cambridge University Press: Cambridge, 1982.
  * [37] Maynard Smith, J., Price, G.R. The logic of animal conflict. Nature 246, 15–18, 1973.
  * [38] May, R.M. Stability and Complexity in Model Ecosystems. Princeton University Press, Princeton, 1973.
  * [39] Moran, P.A.P. Random processes in genetics. Proceedings of the Cambridge Philosophical Society 54, 60-71, 1958.
  * [40] Nagylaki, T. The strong-migration limit in geographically structured populations. Journal of Mathematical Biology 9, 101-114, 1980.
  * [41] Nagylaki, T. Geographical invariance in population genetics. Journal of Theoretical Biology 99, 159-172, 1982.
  * [42] Nagylaki, T. Geographical invariance and the strong-migration limit in subdivided populations. Journal of Mathematical Biology 41, 123-142, 2000.
  * [43] Nowak, M.A., Sasaki, A., Taylor, S., Fudenberg, D. Emergence of cooperation and evolutionary stability in finite populations. Nature 428, 646-650, 2004.
  * [44] Nowak, M.A. Evolutionary Dynamics. The Belknap Press of Harvard University Press: Cambridge, MA, USA, 2006.
  * [45] Otto, S.P., Whitlock, M.C. The probability of fixation in populations of changing size. Genetics 146, 723-733, 1997.
  * [46] Parsons, T.L., Quince, C. Fixation in haploid populations exhibiting density dependence I: the non-neutral case. Theoretical Population Biology 72, 121-135, 2007.
  * [47] Parsons, T.L., Quince, C. Fixation in haploid populations exhibiting density dependence II: the quasi-neutral case. Theoretical Population Biology 72, 468-479, 2007.
  * [48] Roze, D., Rousset, F. Selection and drift in subdivided populations: a straightforward method for deriving diffusion approximations and applications involving dominance, selfing and local extinctions. Genetics 165, 2153-2166, 2003.
  * [49] Rousset, F. Inbreeding and relatedness coefficients: what do they measure? Heredity 88, 371-380, 2002.
  * [50] Slatkin, M. Fixation probabilities and fixation times in a subdivided population. Evolution  35, 477-488, 1981.
  * [51] Slatkin, M. Gene flow in natural populations. Annual Review of Ecology, Evolution, and Systematics 16, 393-430, 1985.
  * [52] Sugden, R. The Economics of Rights, Co-operation and Welfare. Blackwell, Oxford, 1986.
  * [53] Takahata, N. Genealogy of neutral genes and spreading of selected mutations in a geographically structured population. Genetics  129, 585-595, 1991.
  * [54] Taylor, P.D., Jonker, L.B. Evolutionary stable strategies and game dynamics. Mathematical Biosciences 40, 145-156, 1978.
  * [55] Uecker, H., Hermisson, J. On the fixation process of a beneficial mutation in a variable environment. Genetics  188, 915-930, 2011.
  * [56] Wakeley, J. Polymorphism and divergence for island-model species. Genetics 163, 411-420, 2007.
  * [57] Wakeley, J., Takahashi, T. The many-demes limit for selection and drift in a subdivided population. Theoretical Population Biology 66, 83-91, 2004
  * [58] Whitlock, M.C. Fixation probability and time in subdivided populations. Genetics 164, 767-779, 2003.
  * [59] Wild, G., Taylor, P.D., Fitness and evolutionary stability in game theoretic models of finite populations. Proceedings of the Royal Society B 271, 2345-2349, 2004
  * [60] Wright, S. Evolution in Mendelian populations. Genetics  16, 97-159, 1931.
  * [61] Zeeman, R.C. Populations dynamics from game theory. In Global Theory of Dynamical Systems; Nitecki, Z.H., Robinson, R.C., Editors; Springer: New York, 1980.
  * [62] Zhang, B.Y., Fan, S.J., Li, C., Zheng, X.D., Bao, J.Z., Cressman, R., Tao, Y. Opting out against defection leads to stable coexistence with cooperation. Scientific Reports 6, 35902, 2016.
  * [63] Zheng, X.D., Li, C., Lessard, S., Tao, Y., Evolutionary stability concepts in a stochastic environment, Physical Review E 96, 032414, 2017.
  * [64] Zheng, X.D., Li, C., Lessard, S., Tao, Y., Environmental noise could promote stochastic local stability of behavioral diversity evolution, Physical Review Letters 120, 218101, 2018