The molecular clock of neutral evolution can be accelerated or slowed by
asymmetric spatial structure

Benjamin Allen1,2,3,∗, Christine Sample1, Yulia A. Dementieva1, Ruben C.
Medeiros1, Christopher Paoletti1, Martin A. Nowak2,4

1 Department of Mathematics, Emmanuel College, Boston, MA, USA
2 Program for Evolutionary Dynamics, Harvard University, Cambridge, MA, USA
3 Center for Mathematical Sciences and Applications, Harvard University,
Cambridge, MA, USA
4 Department of Mathematics, Department of Organismic and Evolutionary
Biology, Harvard University, Cambridge, MA, USA
$\ast$ E-mail<EMAIL_ADDRESS>
## Abstract

Over time, a population acquires neutral genetic substitutions as a
consequence of random drift. A famous result in population genetics asserts
that the rate, $K$, at which these substitutions accumulate in the population
coincides with the mutation rate, $u$, at which they arise in individuals:
$K=u$. This identity enables genetic sequence data to be used as a “molecular
clock” to estimate the timing of evolutionary events. While the molecular
clock is known to be perturbed by selection, it is thought that $K=u$ holds
very generally for neutral evolution. Here we show that asymmetric spatial
population structure can alter the molecular clock rate for neutral mutations,
leading to either $K<u$ or $K>u$. Deviations from $K=u$ occur because
mutations arise unequally at different sites and have different probabilities
of fixation depending on where they arise. If birth rates are uniform across
sites, then $K\leq u$. In general, $K$ can take any value between 0 and $Nu$.
Our model can be applied to a variety of population structures. In one
example, we investigate the accumulation of genetic mutations in the small
intestine. In another application, we analyze over 900 Twitter networks to
study the effect of network topology on the fixation of neutral innovations in
social evolution.

## Author Summary

Evolution is driven by genetic mutations. While some mutations affect an
organism’s ability to survive and reproduce, most are neutral and have no
effect. Neutral mutations play an important role in the study of evolution
because they generally accrue at a consistent rate over time. This result,
first discovered almost 50 years ago, allows neutral mutations to be used as a
“molecular clock” to estimate, for example, how long ago humans diverged from
chimpanzees and bonobos. We used mathematical modeling to study how the rates
of these molecular clocks are affected by the spatial arrangement of a
population in its habitat. We find that asymmetry in this spatial structure
can either slow down or speed up the rate at which neutral mutations accrue.
This effect could potentially skew our estimates of past events from genetic
data. It also has implications for a number of other fields. For example, we
show that the architecture of intestinal tissue can limit the rate of genetic
substitutions leading to cancer. We also show that the structure of social
networks affects the rate at which new ideas replace old ones. Surprisingly,
we find that most Twitter networks slow down the rate of idea replacement.

## Introduction

A half-century ago, Zuckerkandl and Pauling [1] discovered that amino acid
substitutions often occur with sufficient regularity as to constitute a
“molecular clock”. Theoretical support for this observation was provided by
Kimura [2], who argued that observed rates of amino acid substitution could
only be explained if the majority of substitutions are selectively neutral.
Under simple models of evolution, a neutral mutation has probability $1/N$ of
becoming fixed in a haploid population of size $N$. It follows that the rate
$K$ of neutral substitution per generation—given by the product of the
population size $N$, the mutation probability $u$ per reproduction, and the
fixation probability $\rho$—is simply equal to $u$. (A similar cancellation
occurs in diploids, leading again to $K=u$.) In other words, for any neutral
genetic marker, the rate of substitution at the population level equals the
rate of mutation at the individual level [2].

A number of factors can alter the rate of neutral substitution [3, 4],
including selection, changes in population demography over time, or mutation
rates that vary systematically by demographic classes [5, 6] or sex [7]. The
extent to which these factors compromise the applicability of the molecular
clock hypothesis to biological sequence data has been intensely debated [8, 9,
10, 11, 12, 13, 3, 4, 14].

However, it is generally thought that spatial structure alone (without spatial
variation in the mutation rate[7]) cannot alter the rate of neutral
substitution. This consensus is based on analyses [15, 16, 17, 18, 19, 20, 21,
22, 23, 24, 25, 26, 27, 28, 29, 30, 31] of a variety of models of spatially
structured populations. Each of these analyses found that the fixation
probability of a neutral mutation is $1/N$, and thus the rate of neutral
substitution $K=Nu\rho$ again simplifies to $u$.

Here we show that the absence of spatial effects on the rate of neutral
substitution in these models does not represent a general principle of
evolution. Rather, it is an artifact of two common modeling assumptions: (i)
all spatial locations are equivalent under symmetry, and (ii) mutations are
equally likely to arise in each location. While assumption (i) is relaxed in a
number of models, assumption (ii) is made almost universally. Either of these
assumptions alone is sufficient to guarantee $\rho=1/N$ and $K=u$, as we will
show. However, neither of these assumptions is necessarily satisfied in a
biological population. In particular, assumption (ii) is violated if some
spatial locations experience more turnover (births and deaths) than others.
Since each birth provides an independent opportunity for mutation, the rate at
which new mutations appear at a location is proportional to its turnover rate
[32]. Thus, even with a constant probability of mutation per birth, mutations
may appear with different frequency at different locations. If, in addition,
fixation probability depends on a mutant’s initial location, the rate of
neutral substitution is altered.

Our goal is to identify conditions under which the molecular clock rate is
maintained ($K=u$), accelerated ($K>u$), or slowed ($K<u$) by spatial
population structure. Our main results are as follows (see also Figure 1):

Figure 1: How spatial structure effects the molecular clock rate $K$. Relative
to the rate in a well-mixed population ($K=u$), spatial structure can either
accelerate ($K>u$) or slow ($K<u$) the accumulation of neutral substitutions,
depending on how birth rates $b_{i}$ and death rates $d_{i}$ vary across
sites. The rate is unchanged from that of a well-mixed population ($K=u$) if
either death rates are uniform across sites (Result 1), or the birth and death
rates are equal at each site (Result 2). Almost all previous studies of
neutral drift in spatially structured populations fall into one of these two
categories; thus the effects of spatial structure on the molecular clock rate
are unappreciated. We show that, in general, $K$ can take any non-negative
value less than $Nu$ (Result 4). If one adds the constraint that the birth
rate is the same at each site, then the molecular clock rate cannot exceed
that of a well-mixed population ($K\leq u$; Result 3).

  * Result 1 If the death rate is constant over sites, the molecular clock rate is identical to that of a well-mixed population: $K=u$.

  * Result 2 If the birth rate equals the death rate at each site, the molecular clock rate is identical to that of a well-mixed population: $K=u$.

  * Result 3 If the birth rate is constant over sites, the molecular clock rate is less than or equal to that of a well-mixed population: $K\leq u$.

  * Result 4 In general (with no constraints on birth or death rates), the molecular clock rate $K$ can take any non-negative value less than $Nu$.

## Results

### Obtaining the neutral substitution rate

Our investigations apply to a class of evolutionary models (formally described
in the Materials and Methods) in which reproduction is asexual and the
population size and spatial structure are fixed. Specifically, there are a
fixed number of sites, indexed $i=1,\ldots,N$. Each site is always occupied by
a single individual. At each time-step, a replacement event occurs, meaning
that the occupants of some sites are replaced by the offspring of others. The
replacement event is chosen according to a probability distribution—called the
replacement rule—which depends on the model in question. Since we consider
only neutral mutations that have no phenotypic effect, the probabilities of
replacement events are the same for every population state.

This class includes many established evolutionary models. One important
subclass is spatial Moran processes [33, 34, 23, 35], in which exactly one
reproduction occurs each time-step. This class also includes spatial Wright-
Fisher processes, in which the entire population is replaced each time-step
[36, 37]. In general, any subset of individuals may be replaced by the
offspring of any other subset in a given time-step.

For a given model in this class, we let $e_{ij}$ denote the probability that
the occupant of site $j$ is replaced by the offspring of site $i$ in a single
time-step. Thus the expected number of offspring of site $i$ over a single
time-step is $b_{i}=\sum_{j=1}^{N}e_{ij}$. The probability that node $i$ dies
(i.e., is replaced) in a time-step is $d_{i}=\sum_{j=1}^{N}e_{ji}$. The death
rate $d_{i}$ can also be regarded as the rate of turnover at site $i$. The
total expected number of offspring per time-step is denoted
$B=\sum_{i=1}^{N}b_{i}=\sum_{i=1}^{N}d_{i}=\sum_{i,j}e_{ij}$. We define a
generation to be $N/B$ time-steps, so that, on average, each site is replaced
once per generation.

We use this framework to study the fate of a single neutral mutation, as it
arises and either disappears or becomes fixed. The probability of fixation
depends on the spatial structure and the initial mutant’s location. We let
$\rho_{i}$ denote the probability that a new mutation arising at site $i$
becomes fixed. ($\rho_{i}$ can also be understood as the reproductive value of
site $i$ [38].) We show in the Materials and Methods that the fixation
probabilities $\rho_{i}$ are the unique solution to the system of equations

$\begin{cases}\displaystyle
d_{i}\rho_{i}=\sum_{j=1}^{N}e_{ij}\rho_{j}&\text{for $i=1,\ldots,N$},\\\
\displaystyle\sum_{i=1}^{N}\rho_{i}=1.\end{cases}$ (1)

The identity $\sum_{i=1}^{N}\rho_{i}=1$ arises because $\rho_{i}$ equals the
probability that the current occupant of site $i$ will become the eventual
ancestor of the population, which is true for exactly one of the $N$ sites.

To determine the overall rate of substitution, we must take into account the
likelihood of mutations arising at each site. The rate at which mutations
arise at site $i$ is proportional to the turnover rate $d_{i}$, because each
new offspring provides an independent chance of mutation. Specifically, if
mutation occurs at rate $u\ll 1$ per reproduction, then new mutations arise at
site $i$ at rate $Nud_{i}/B$ per generation [32]. Thus the fraction of
mutations that arise at site $i$ is $d_{i}/B$.

The overall fixation probability $\rho$ of new mutations, taking into account
all possible initial sites, is therefore

$\rho=\frac{1}{B}\sum_{i=1}^{N}d_{i}\rho_{i}.$ (2)

The molecular clock rate $K$ is obtained by multiplying the fixation
probability $\rho$ by the total rate of mutation per generation:

$K=Nu\rho=\frac{Nu}{B}\sum_{i=1}^{N}d_{i}\rho_{i}.$ (3)

The units of $K$ are substitutions per generation. Alternatively, the
molecular clock can be expressed in units of substitutions per time-step, in
which case the fomula is $\tilde{K}=Bu\rho=u\sum_{i=1}^{N}d_{i}\rho_{i}$.

### Effects of spatial structure

How does spatial structure affect the rate of neutral substitution? In a well-
mixed population, each individual’s offspring is equally likely to replace
each other individual, meaning that $e_{ij}$ is constant over all $i$ and $j$
(Fig. 2a). In this case, the unique solution to Eq. (1) is $\rho_{i}=1/N$ for
all $i$, and we recover Kimura’s [2] result $K=Nu(1/N)=u$. Moreover, if each
site is equivalent under symmetry, as in Figure 2b, this symmetry implies that
$\rho_{i}=1/N$ for all $i$ and $K=u$ as in the well-mixed case.

|
---|---
(a) | (b)
Figure 2: For spatial structures with symmetry, $K=u$. (a) For a well-mixed
population, represented by a complete graph with uniform edge weights, a
neutral mutation has a $1/N$ chance of fixation, where $N$ is the population
size. It follows that the rate $K$ of neutral substitution in the population
equals the rate $u$ of neutral mutation in individuals. (b) The same result
holds for spatial structures in which each site is _a priori_ identical, such
as the cycle with uniform edge weights.

However, asymmetric spatial structure can lead to faster ($K>u$) or slower
($K<u$) molecular clock rates than a well-mixed population, as shown in Figure
3. These results led us to seek general conditions for when spatial structure
lead to faster, slower, or the same molecular clock rates as a well-mixed
population.

|
---|---
(a) | (b)
Figure 3: Asymmetric spatial structure affects the rate of neutral
substitution. This is because the frequency of mutations and the probability
of fixation differ across sites. Turnover rates are indicated by coloration,
with red corresponding to frequent turnover and consequently frequent
mutation. (a) A star consists of a hub and $n$ leaves, so that the population
size is $N=n+1$. Edge weights are chosen so that the birth rates are uniform
($b_{i}=1$ for all $i$). Solving Eq. (1), we obtain site-specific fixation
probabilities of $\rho_{\mathrm{H}}=1/(1+n^{2})$ and
$\rho_{\mathrm{L}}=n/(1+n^{2})$ for the hub and each leaf, respectively. From
Eq. (3), the molecular clock rate is $K=\frac{2n}{1+n^{2}}u$, which equals $u$
for $n=1$ and is less than $u$ for $n\geq 2$. Thus the star structure slows
down the rate of neutral substitution, in accordance with Result 3.
Intuitively, the slowdown occurs because mutations are more likely to arise at
the hub, where their chances of fixation are reduced. (b) A one-dimensional
population with self-replacement only in site 1. Solving Eq. (1) we find
$\rho_{1}=\frac{8}{15}$ and $\rho_{2}=\frac{4}{15}$, $\rho_{3}=\frac{2}{15}$
and $\rho_{4}=\frac{1}{15}$. (The powers of two arise because there is twice
as much gene flow in one direction as the other). From Eq. (3), the molecular
clock rate is $K=\frac{16}{13}u>u$, thus the molecular clock is accelerated in
this case.

We first find

######  Result 1.

If the death rates $d_{i}$ are constant over all sites $i=1,\ldots,N$, then
$\rho=1/N$, and consequently $K=u$.

Thus the molecular clock rate is unaffected by spatial structure if each site
is replaced at the same rate (Fig. 4a). This result can be seen by noting that
if the $d_{i}$ are constant over $i$, then since $\sum_{i=1}^{N}d_{i}=B$, it
follows that $d_{i}=B/N$ for each $i$. Substituting in Eq. (2) yields
$\rho=1/N$.

|
---|---
(a) | (b)
Figure 4: Results 1 and 2 give conditions leading to $K=u$. (a) Our Result 1
states that the molecular clock has the same rate as in a well-mixed
population, $K=u$, if the rate of turnover $d_{i}$ is uniform across sites, as
in this example ($d_{i}=0.2$ for all $i$). (b) Result 2 asserts that
$\rho_{i}=1/N$ for all $i$—again implying $K=u$—if and only if each site has
birth rate equal to death rate, $b_{i}=d_{i}$ for all $i$, as in this example.
Nodes are colored according to their rates of turnover $d_{i}$.

Another condition leading to $K=u$ is the following:

######  Result 2.

If the birth rate equals the death rate at each site ($b_{i}=d_{i}$ for all
$i=1,\ldots,N$), then $\rho=1/N$, and consequently $K=u$. Moreover,
$b_{i}=d_{i}$ for all $i=1,\ldots,N$ if and only if the fixation probability
is the same from each site ($\rho_{i}=1/N$ for all $i=1,\ldots,N$).

Thus if births and deaths are balanced at each site, then all sites provide an
equal chance for mutant fixation (Fig. 4b). In this case the molecular clock
is again unchanged from the baseline value. In particular, if dispersal is
symmetric in the sense that $e_{ij}=e_{ji}$ for all $i$ and $j$ then $K=u$.
Result 2 can be obtained by substituting $\rho_{i}=1/N$ for all $i$ into Eq.
(1) and simplifying to obtain $b_{i}=d_{i}$ for all $i$ (details in Methods).
Alternatively, Result 2 can be obtained as a corollary to the Circulation
Theorem of Lieberman et al. [23].

Our third result reveals a “speed limit” to neutral evolution in the case of
constant birth rates:

######  Result 3.

If the birth rates $b_{i}$ are constant over all sites $i=1,\ldots,N$, then
$\rho\leq 1/N$, and consequently $K\leq u$, with equality if and only if the
death rates $d_{i}$ are also constant over sites.

In other words, a combination of uniform birth rates and nonuniform death
rates slows down the molecular clock. An instance of this slowdown in shown in
Figure 3a. Intuitively, the sites at which mutations occur most frequenly are
those with high death rates $d_{i}$; because of these high death rates, these
sites on the whole provide a reduced chance of fixation. The proof of this
result, however, is considerably more intricate than this intuition would
suggest (see Methods).

Finally, we investigate the full range of possible values for $K$ with no
constraints on birth and death rates. We find the following:

######  Result 4.

For arbitrary spatial population structure (no constraints on $e_{ij}$) the
fixation probability can take any value $0\leq\rho<1$, and consequently, the
molecular clock can take any rate $0\leq K<Nu$.

This result is especially surprising, in that it implies that the probability
of fixation of a new mutation can come arbitrarily close to unity. Result 4
can be proven by considering the hypothetical spatial structure illustrated in
Figure 5. Any non-negative value of $\rho$ less than 1 can be obtained by an
appropriate choice of parameters (details in Methods).

Figure 5: Result 4 shows that $K$ can achieve any value $0\leq K<Nu$. This is
proven by considering a population structure with unidirectional gene flow
from a hub (H) to $N-1$ leaves (L). Fixation is guaranteed for mutations
arising in the hub ($\rho_{\mathrm{H}}=1$) and impossible for those arising in
leaves ($\rho_{\mathrm{L}}=0$). The overall fixation probability is equal by
Eq. (2) to the rate of turnover at the hub: $\rho=d_{\mathrm{H}}=1-(N-1)a$.
The molecular clock rate is therefore $K=Nu\rho=N[1-(N-1)a]u$. It follows that
$K>u$ if and only if $a<1/N$. Intuitively, the molecular clock is accelerated
if the hub experiences more turnover (and hence more mutations) than the other
sites. Any value of $\rho$ greater than or equal to 0 and less than 1 can be
achieved through a corresponding positive choice of $a$ less than or equal to
$1/(N-1)$. For $a=1/(N-1)$ we have $K=0$, because mutations arise only at the
leaves where there is no chance of fixation. At the opposite extreme, in the
limit $a\to 0$, we have $K\to Nu$.

### Application to upstream-downstream populations

To illustrate the effects of asymmetric dispersal on the molecular clock, we
consider a hypothetical population with two subpopulations, labeled “upstream”
and “downstream” (Figure 6). The sizes of these subpopulations are labeled
$N_{\uparrow}$ and $N_{\downarrow}$, respectively. Each subpopulation is well-
mixed, with replacement probabilities $e_{\uparrow}$ for each pair of upstream
sites and $e_{\downarrow}$ for each pair of downstream sites. Dispersal
between the subpopulations is represented by the replacement probabilities
$e_{\rightarrow}$ from each upstream site to each downstream site, and
$e_{\leftarrow}$ from each downstream site to each upstream site. We assume
there is net gene flow downstream, so that $e_{\rightarrow}>e_{\leftarrow}$.

Figure 6: A model of a population divided into upstream and downstream
subpopulations. Each subpopulation is well-mixed. The replacement probability
$e_{ij}$ equals $e_{\uparrow}$ if sites $i$ and $j$ are both upstream,
$e_{\downarrow}$ if $i$ and $j$ are both downstream, $e_{\rightarrow}$ if $i$
is upstream and $j$ is downstream, and $e_{\leftarrow}$ if $i$ is downstream
and $j$ is upstream. We suppose there is net gene flow downstream, so that
$e_{\rightarrow}>e_{\leftarrow}$. We find that the molecular clock is
accelerated, relative to the well-mixed case, if and only if the upstream
subpopulation experiences more turnover than the downstream subpopulation:
$K>u$ if and only if $d_{\uparrow}>d_{\downarrow}$.

Solving Eq. (1), we find that the fixation probabilities from each upstream
site and each downstream site, respectively, are

$\begin{split}\rho_{\uparrow}=\frac{e_{\rightarrow}}{N_{\uparrow}e_{\rightarrow}+N_{\downarrow}e_{\leftarrow}},\\\
\rho_{\downarrow}=\frac{e_{\leftarrow}}{N_{\uparrow}e_{\rightarrow}+N_{\downarrow}e_{\leftarrow}}.\end{split}$
(4)

These fixation probabilities were previously discovered for a different model
of a subdivided population [39]. Substituting these fixation probabilities
into Eq. (3) yields the molecular clock rate:

$K=\frac{N}{B}\,\frac{N_{\uparrow}d_{\uparrow}e_{\rightarrow}+N_{\downarrow}d_{\downarrow}e_{\leftarrow}}{N_{\uparrow}e_{\rightarrow}+N_{\downarrow}e_{\leftarrow}}\,u.$
(5)

Above, $d_{\uparrow}$ and $d_{\downarrow}$ are the turnover rates in the
upstream and downstream populations, respectively, and
$B=N_{\uparrow}d_{\uparrow}+N_{\downarrow}d_{\downarrow}$ is the total birth
rate per time-step. In Methods, we show that $K>u$ if and only if
$d_{\uparrow}>d_{\downarrow}$; that is, the molecular clock is accelerated if
and only if there is more turnover in the upstream population than in the
downstream population.

### Application to epithelial cell populations

Our results are also applicable to somatic evolution in self-renewing cell
populations, such as the crypt-like structures of the intestine. Novel
labeling techniques have revealed that neutral mutations accumulate in
intestinal crypts at a constant rate over time [40]. The cell population in
each crypt is maintained by a small number of stem cells that reside at the
crypt bottom and continuously replace each other in stochastic manner (Figure
7; [41, 42, 43]). We focus on the proximal small intestine in mice, for which
recent studies [40, 44] suggest there are $\sim 5$ active stem cells per
crypt, each replaced $\sim 0.1$ times per day by one of its two neighbors. In
our framework, this corresponds to a cycle-structured population of size 5
with replacement rates $0.05/\mathrm{day}$ between neighbors, so that
$d_{i}=0.1/\mathrm{day}$ for all $i$.

Figure 7: A simple model of cell replacement structure in epithelial crypts of
the small intestine, based on results of [40] and [44]. A small number of stem
cells ($N_{\uparrow}\sim 5$) residing at the bottom of the intestinal crypt
and are replaced at rate $d_{\uparrow}\sim 0.1$ per stem cell per day.
Empirical results [40, 44] suggest a cycle structure for stem cells. To
achieve the correct replacement rate we set $e_{ij}=0.05/\mathrm{day}$ for
each neighboring pair. Stem cells in an individual crypt replace a much larger
number of progenitor and differentiated cells ($\sim 250$; [45]). These
downstream progenitor and differentiated cells are replaced about every day
[45]. The hierarchical organization of intestinal crypts, combined with the
low turnover rate of stem cells, limits the rate of neutral genetic
substitutions ($\tilde{K}\approx 0.1u$ substitutions per day), since only
mutations that arise in stem cells can fix.

Only mutations that arise in stem cells can become fixed within a crypt; thus
we need only consider the fixation probabilities and turnover rates among stem
cells. By symmetry among the stem cells, $\rho_{i}=1/5$ for each of the five
stem cell sites. The molecular clock rate is therefore
$\tilde{K}=u\sum_{i=1}^{5}d_{i}\rho_{i}=0.1u$ substitutions per day. This
accords with the empirical finding that, for a neutral genetic marker with
mutation rate $u\approx 1.1\times 10^{-4}$, substitutions accumulate at a rate
$\tilde{K}\approx 1.1\times 10^{-5}$ per crypt per day [40].

Does crypt architecture limit the rate of genetic change in intestinal tissue?
Intestinal crypts in mice contain $\sim 250$ cells and replace all their cells
about once per day [45]. If each crypt were a well-mixed population, the
molecular clock rate would be $\tilde{K}=Bu/N\approx u$ substitutions per day.
Thus the asymmetric structure of these epithelial crypts slows the rate of
neutral genetic substitution tenfold.

### Application to the spread of innovations

Our results can also be applied to innovations that spread by imitation on
social networks. In this setting, a mutation corresponds to an innovation that
could potentially replace an established convention. Neutrality means that the
innovation is equally likely to be imitated as the established convention.

To investigate whether human social networks promote or hinder the fixation of
selectively neutral innovations, we analyzed 973 Twitter networks from the
Stanford Large Network Dataset Collection [46]. Each of these “ego networks”
represents follower relationships among those followed by a single “ego”
individual (who is not herself included in the network). We oriented the links
in each network to point from followee to follower, corresponding to the
presumed direction of information flow. Self-loops were removed. To ensure
that fixation is possible, we eliminated individuals that could not be
reached, via outgoing links, from the node with greatest eigenvector
centrality. The resulting networks varied in size from 3 to 241 nodes (Fig.
8). All links were assigned equal weight, so that the molecular clock rate (as
a multiple of $u$) depends only on the network topology.

We found that the mean value of $K$ among these ego networks is $0.557u$, with
a standard deviation of $0.222u$. 19 of the 973 networks ($2\%$) have $K>u$.
Two networks have $K=u$ exactly; each of these has $N=3$ nodes and uniform in-
degree $d_{i}$, thus $K=u$ follows from Result 1 for these networks. We found
a weak but statistically significant negative relationship between $N$ and
$K/u$ ($\mathrm{slope}\approx-0.0016,\,R\approx-0.45,\,p<10^{-48}$). In
summary, while some Twitter ego-networks accelerate the substitution of
neutral innovations, the vast majority (especially the larger ones) slow this
rate.

|
---|---
Smallest $K$ | Largest $K$
$N=185,\;K\approx 0.088u$ | $N=66,\;K\approx 1.264u$
(a) | (b)
|
$N=3,\;K=u$ | $N=199,\;K\approx 1.065u$
(c) | (d)
|
$N=80,\;K\approx 0.473u$ |
(e) | (f)
Figure 8: The rate of neutral substitution $K$ on Twitter “ego networks”.
(a-e) Five of the 973 networks analyzed, including those with (a) the largest
value of $K$, (b) the smallest value of $K$, and (c) the fewest nodes. (f) A
scatter plot of $K/u$ versus $N$ reveals a weak negative correlation
($\mathrm{slope}\approx-0.0016,\,R\approx-0.45,\,p<10^{-48}$). The colored
dots on the scatter plot correspond to the networks shown in (a-e). The dashed
line corresponds to $K/u=1$, above which network topology accelerates neutral
substitution.

## Discussion

The spatial structure of a population affects its evolution in many ways, for
example by promoting cooperative behaviors [36, 47, 48, 49, 50, 51, 52, 53],
genetic variation [54, 55, 56, 57], and speciation [58, 59, 60]. Asymmetric
spatial structure in particular is known to have important consequences for
adaptation [61, 62, 23, 35, 63, 27, 28, 101] and for genetic diversity [64,
65, 66]. Our work shows that asymmetric spatial structure also affects the
rate of neutral substitution.

In light of Results 1 and 2, we see that the critical factors driving the
changes in molecular clock rate are differential rates of turnover ($d_{i}$)
and net offspring production ($b_{i}-d_{i}$) across sites. If both $d_{i}$ and
$b_{i}-d_{i}$ differ across sites, the molecular clock rate will in general
differ from that of a well-mixed population (Result 4). If additionally
$b_{i}$ is constant across sites, Result 3 guarantees that neutral
substitution is slowed relative to the well-mixed case.

While our class of models includes a number of limiting assumptions (e.g.,
constant population size), our qualitative results do not appear to depend on
these assumptions. We expect that differential rates of turnover and net
offspring production will also alter molecular clock rates in models outside
of the class considered here. However, the result that molecular clock rates
can come abritrarily close to $Nu$—which depends on the existence of a single
“hub” individual seeding the rest of the population as in Figure 5—may require
modification for sexually reproducing populations.

Conditions leading to altered molecular clock rate may occur frequently in
natural populations. Asymmetric dispersal may result, for example, from
prevailing winds [67, 68], water currents [69, 70], or differences in
elevation [71, 72]. Differences in habitat quality may lead to variance in
birth and death rates across sites (nonuniform $b_{i}$ and $d_{i}$). It is
known that such asymmetries in spatial structure have important consequences
for adaptation [61, 62, 23, 63] and for genetic diversity [64, 65, 66]; our
work shows that they also have consequences for the rate of neutral genetic
change. In particular, the molecular clock is accelerated if there is greater
turnover in “upstream” subpopulations.

Many somatic cell populations have strongly asymmetric patterns of
replacement, with small stem cell pools feeding much larger populations of
progenitor and differentiated cells. Our results support the idea that
infrequent turnover of stem cells suppresses the accumulation of somatic
mutations [73, 74, 44, 75, 76]. The rate at which these mutations accumulate
has significant implications for the onset of cancer [77, 78, 79] and the
likelihood of successful cancer therapy [80, 81, 82, 83]. It is important to
note, however, that cancer can alter the structure of cell hierarchies; for
example, by altering the number and replacement rate of stem cells [40] and/or
allowing differentiated cells to revert to stem cells [84]. This restructuring
may, in turn, alter the rate of genetic substitution.

The influence of social network topology on the spread of ideas and behaviors
is a question of both theoretical and practical interest [85, 86, 87, 88, 89,
90, 91, 92]. The neutral substition rate $K$ on social networks describes how
innovations spread when they are equally likely to be imitated as an existing
convention. Since it depends only on network topology, $K/u$ can be used as a
statistic to isolate how this topology affects the rate of social change. Our
finding that most Twitter ego networks have rates less than those of well-
mixed populations contrasts with results from epidemiological models, which
generally find that the heterogeneity of real-world social networks
accelerates contagion [85, 93, 94]. Further research is needed to determine
which kinds of real-world networks accelerate neutral substitution, and how
neutral substitution rates relate to other measures of diffusion and contagion
on networks.

If spatial structure remains constant over time, then the neutral substitution
rate $K$ is in all cases a constant multiple of the mutation rate $u$. In this
case, absent other complicating factors such as selection, neutral mutations
will accrue at a constant rate that can be inferred from genetic data.
However, if the spatial sturcture changes over time—due, for instance, to
changes in climate, tumorogenesis, or social network dynamics [95, 96, 97,
98]—the rate of neutral substitutions may change over time as well.

In our framework, the molecular clock rate is assumed to depend only on the
rate at which mutations arise and their probability of becoming fixed. This
approach assumes that the time to fixation is typically shorter than the
expected waiting time $1/(Nu\rho)$ for the next successful mutation. If this
is not the case, then substitution rates are also affected by fixation times.
These fixation times are themselves affected by spatial structure [57, 22, 99,
30], leading to further ramifications for the molecular clock rate [100].

The starting point of our analysis is that the convention, commonly assumed in
evolutionary models, that mutations arise with equal frequency at each site,
is not necessarily the most natural choice. If there is a constant probability
of mutation per birth, then mutations instead arise in proportion to the rate
of turnover at a site. Here we have applied this principle to study the rate
of neutral substitution. However, this principle also holds for advantageous
and deleterious mutations, as well as those whose effect varies with location.
It also applies to frequency-dependent selection [101, 102]. Re-analyzing
existing models using this new mutation convention may reveal further
surprises about how spatial structure affects evolution.

## Methods

### Class of Models

In the class of models we consider, there are $N$ sites indexed
$i=1,\ldots,N$, each always occupied by a single individual. At each time-
step, a replacement event occurs, in which the occupants of some positions are
replaced by the offspring of others. A replacement event is identified by a
pair $(R,\alpha)$, where $R\subset\\{1,\ldots,N\\}$ is the set of sites whose
occupants are replaced by new offspring, and $\alpha:R\to\\{1,\ldots,N\\}$ is
a set mapping indicating the parent of each new offspring. (This notation was
introduced in Ref. [32].) A sample replacement event is illustrated in Figure
9.

Figure 9: Illustration of a replacement event. In this case, the occupant of
site 3 produced one offspring, which displaced the occupant of site 2. The
occupant of site 4 produced two offspring, one remaining in site 4 and
displacing the parent, and the other displacing the occupant of site 5. Thus
the set of replaced positions is $R=\\{2,4,5\\}$ and the offspring-to-parent
map $\alpha$ is given by $\alpha(2)=3,\alpha(4)=4,\alpha(5)=4$. There is no
mutation in the example illustrated here, so each offspring inherits the type
of the parent. Thus the population transitions from state
$\mathbf{s}=(\mathrm{M},\mathrm{R},\mathrm{R},\mathrm{M},\mathrm{R})$ to
$\mathbf{s}^{\prime}=(\mathrm{M},\mathrm{R},\mathrm{R},\mathrm{M},\mathrm{M})$.

A model of neutral evolution is specified by a probability distribution over
the set of possible replacement events. We call this probability distribution
the replacement rule of the model. The probability of a replacement event
$(R,\alpha)$ in this distribution will be denoted $p(R,\alpha)$. Neutrality is
represented by independence of the probabilities $p(R,\alpha)$ from the state
of the evolutionary process.

The only assumption we place on the replacement rule is that it should be
possible for at least one site $i$ to contain the eventual ancestor of the
population:

######  Assumption 1.

There is an $i\in\\{1,\ldots,N\\}$, a positive integer $n$ and a finite
sequence $\\{(R_{k},\alpha_{k})\\}_{k=1}^{n}$ of replacement events such that

  * •

$p(R_{k},\alpha_{k})>0$ for all $k$, and

  * •

For all individuals $j\in\\{1,\ldots,N\\}$,

$\alpha_{k_{1}}\circ\alpha_{k_{2}}\circ\cdots\circ\alpha_{k_{m}}(j)=i,$ (6)

where $k_{1},\ldots<k_{m}$ is the maximal subsequence of $1,\ldots,n$ such
that the compositions in Eq. (6) are well-defined.

We observe that Eq. (6) traces the ancestors of $j$ backwards in time to $i$.
This assumption is equivalent to saying that there is at least one site $i$
such that mutations arising at site $i$ have nonzero fixation probability.
Assumption 1 precludes degenerate cases such as two completely separate
subpopulations with no gene flow between them, in which case fixation would be
impossible.

For a specific replacement event $(R,\alpha)$, the sites that are replaced by
the offspring of $i$ is the given by the preimage $\alpha^{-1}(i)\subset R$,
i.e., the set of indices that map to $i$ under $\alpha$. The number of
offspring of site $i$ is equal to $|\alpha^{-1}(i)|$, the cardinality (size)
of this preimage. Taking all possible replacement events into account, the
birth rate (expected offspring number) of site $i$ is given by

$b_{i}=\operatorname{E}\big{[}|\alpha^{-1}(i)|\big{]}=\sum_{(R,\alpha)}p(R,\alpha)\,|\alpha^{-1}(i)|.$

The death rate (probability of replacement) of site $i$ is equal to

$d_{i}=\Pr[i\in R]=\sum_{\begin{subarray}{c}(R,\alpha)\\\ i\in
R\end{subarray}}p(R,\alpha).$

The probability that the offspring of $i$ displaces the occupant of $j$ in a
replacement event is

$e_{ij}=\Pr[\text{$j\in R$ and
$\alpha(j)=i$}]=\sum_{\begin{subarray}{c}(R,\alpha)\\\ i\in R\\\
\alpha(j)=i\end{subarray}}p(R,\alpha).$

We observe that $b_{i}=\sum_{j=1}^{N}e_{ij}$ and $d_{i}=\sum_{j=1}^{N}e_{ji}$.

### The evolutionary Markov chain

To study the fixation of new mutations, we consider evolution with two genetic
types: mutant (M) and resident (R). The type occupying site $i$ in a given
state of the evolutionary process is denoted
$s_{i}\in\\{\mathrm{M},\mathrm{R}\\}$. The overall state of the process can be
recorded as a string $\mathbf{s}=(s_{1},\ldots,s_{N})$ of length $N$ with
alphabet $\\{\mathrm{M},\mathrm{R}\\}$.

We assume that there is no further mutation after an initial mutant appears;
thus offspring faithfully inherit the type of the parent. It follows that if
the current state is $\mathbf{s}=(s_{1},\ldots,s_{N})$ and the replacement
event $(R,\alpha)$ occurs, then the new state
$\mathbf{s}=(s_{1}^{\prime},\ldots,s_{N}^{\prime})$ is given by

$s_{i}^{\prime}=\begin{cases}s_{i}&i\notin R\\\ s_{\alpha(i)}&i\in
R.\end{cases}$

The above assumptions describe a Markov chain on the set of strings of length
$N$ with alphabet $\\{\mathrm{M},\mathrm{R}\\}$. We call this the evolutionary
Markov chain. It is straightforward to show that, from any initial state, this
Markov chain will eventually converge upon one of two absorbing states:
$(\mathrm{M},\ldots,\mathrm{M})$ or $(\mathrm{R},\ldots,\mathrm{R})$ [32]. In
the former case, we say that the mutant type has gone to fixation; in the
latter case we say that the mutant type has disappeared.

### The ancestral Markov chain

It is useful to consider a variation on the evolutionary Markov chain called
the _ancestral Markov chain_ , denoted $\mathcal{A}$. Instead of two types,
the ancestral Markov chain has $N$ types, labelled $1,\ldots,N$, which
correspond to the $N$ members of a “founding generation” of the population.
Evolution proceeds according to the given replacement rule, as described in
the Material and Methods. The states of the ancestral Markov chain are strings
of length $N$ with alphabet $\\{1,\ldots,N\\}$.

The ancestral Markov chain has a canonical initial state
$\mathbf{a}_{0}=(1,\ldots,N)$, in which the type of each individual
corresponds to its location. This initial state identifies the locations of
each founding $(t=0)$ member of the population. The ancestral Markov chain
with initial state $\mathbf{a}_{0}$—in our notation,
$(\mathcal{A},\mathbf{a}_{0})$—has the useful feature that at any time $t\geq
0$, the state $\mathbf{a}(t)=(a_{1}(t),\ldots,a_{N}(t))$ indicates the site
occupied by each individual’s founding ancestor. In other words, if
$a_{j}(t)=i$, then the current occupant of site $j$ is descended from the
founder that occupied site $i$.

To relate the evolutionary Markov chains $\mathcal{A}$ and $\mathcal{M}$,
consider any set mapping
$\gamma:\\{1,\ldots,N\\}\to\\{\mathrm{M},\mathrm{R}\\}$. We think of
$\gamma(i)$ as giving the genetic type (M or R) of each member $i=1,\ldots,N$
of the founding generation. The mapping $\gamma$ induces a mapping
$\tilde{\gamma}$ from states of $\mathcal{A}$ to states of $\mathcal{M}$,
defined by
$\tilde{\gamma}(a_{1},\ldots,a_{N})=(\gamma(a_{1}),\ldots,\gamma(a_{N}))$. For
any state $\mathbf{a}(t)$ of $(\mathcal{A},\mathbf{a}_{0})$, the string
$\tilde{\gamma}(\mathbf{a}(t))\in\\{1,\ldots,N\\}^{N}$ indicates the genetic
type of each individual’s ancestor in the founding generation. Since genetic
types are inherited faithfully, it follows that
$\tilde{\gamma}(\mathbf{a}(t))$ gives the current genetic type of each
individual. Thus if $\mathcal{M}$ and $\mathcal{A}$ follow the same
replacement rule, we have that for any such mapping $\gamma$ and any string
$\mathbf{s}\in\\{\mathrm{M},\mathrm{R}\\}^{N}$,

$\operatorname{P}_{(\mathcal{A},\mathbf{a}_{0})}\,[\tilde{\gamma}(\mathbf{a}(t))=\mathbf{s}]\quad=\quad\operatorname{P}_{(\mathcal{M},\tilde{\gamma}(\mathbf{a}_{0}))}\,[\mathbf{s}(t)=\mathbf{s}].$
(7)

### Basic results on fixation probabilities

We define the _fixation probability from site $i$_, $\rho_{i}$, to be the
probability that, from an initial state a mutant in site $i$ and residents in
all other sites, the mutant type goes to fixation:

$\rho_{i}=\lim_{t\to\infty}\operatorname{P}_{(\mathcal{M},\mathbf{m}_{i})}\,[\mathbf{s}(t)=(\mathrm{M},\ldots,\mathrm{M})]$
(8)

Above, $\mathbf{m}_{i}$ denotes the initial state consisting of an M in
position $i$ and R’s elsewhere. The ordered pair
$(\mathcal{M},\mathbf{m}_{i})$ refers to the Markov chain $\mathcal{M}$ with
initial state $\mathbf{m}_{i}$.
$\operatorname{P}_{(\mathcal{M},\mathbf{m}_{i})}[\mathbf{s}(t)=\mathbf{s}]$
denotes the probability that the state of $(\mathcal{M},\mathbf{m}_{i})$ is
$\mathbf{s}$ at time $t\geq 0$.

We can use the relationship between the ancestral Markov chain $\mathcal{A}$
and the evolutionary Markov chain $\mathcal{M}$ to obtain an alternate
expression for the site-specific fixation probability $\rho_{i}$:

######  Theorem 1.

$\rho_{i}=\lim_{t\to\infty}\operatorname{P}_{(\mathcal{A},\mathbf{a}_{0})}\,[\mathbf{a}(t)=(i,\ldots,i)]$.

In words, the site-specific fixation probability $\rho_{i}$ equals the
probability that founding individual $i$ becomes the eventual ancestor of the
whole population.

###### Proof.

For any $i=1,\ldots,N$, define the set mapping
$\gamma_{i}:\\{1,\ldots,N\\}\to\\{\mathrm{M},\mathrm{R}\\}$ by

$\gamma_{i}(j)=\begin{cases}\mathrm{M}&\text{if $j=i$}\\\
\mathrm{R}&\text{otherwise.}\end{cases}$

(Intuitively, this mapping describes the case that individual $i$ of the
founding generation is a mutant, and all others in the founding generation are
residents.) Note that $\tilde{\gamma}_{i}(\mathbf{a}_{0})=\mathbf{m}_{i}$.
Combining Eq. (7) with the definition of $\rho_{i}$, we obtain

$\displaystyle\rho_{i}$
$\displaystyle=\lim_{t\to\infty}\operatorname{P}_{(\mathcal{M},\mathbf{m}_{i})}\,[\mathbf{s}(t)=(\mathrm{M},\ldots,\mathrm{M})]$
$\displaystyle=\lim_{t\to\infty}\operatorname{P}_{(\mathcal{A},\mathbf{a}_{0})}\,[\tilde{\gamma_{i}}(\mathbf{a}(t))=(\mathrm{M},\ldots,\mathrm{M})]$
$\displaystyle=\lim_{t\to\infty}\operatorname{P}_{(\mathcal{A},\mathbf{a}_{0})}\,[\mathbf{a}(t)=(i,\ldots,i)].\qed$
(9)

More generally, we can consider the probability of fixation from an arbitrary
set of sites. For any set $S\subset\\{1,\ldots,N\\}$, we let $\rho_{S}$ denote
the probability that the mutant type becomes fixed, given the initial state
$\mathbf{m}_{S}$ with mutants occupying the sites specified by $S$ and
residents occupying all other sites:

$\rho_{S}=\lim_{t\to\infty}\operatorname{P}_{(\mathcal{M},\mathbf{m}_{S})}\,[\mathbf{s}(t)=(\mathrm{M},\ldots,\mathrm{M})]$

Site-specific fixation probabilities are additive in the following sense:

######  Theorem 2.

For any set $S\subset\\{1,\ldots N\\}$ of sites,

$\rho_{S}=\sum_{i\in S}\rho_{i}.$ (10)

In particular,

$\sum_{i=1}^{N}\rho_{i}=1.$ (11)

This result has previously been obtained for a number of specific evolutioanry
processes on graphs [103, 104, 38]. Intuitively, the probability of fixation
from the initial state described by $S$ equals the probability that one of the
individuals in a site $i\in S$ becomes the eventual ancestor of the
population. Since this cannot be true of more than one site in $S$, the
overall probability $\rho_{S}$ is obtained by summing over all $i\in S$ the
individual probabilities $\rho_{i}$ that site $i$ contains the eventual
ancestor.

###### Proof.

Suppose that we are given a set $S\subset\\{1,\ldots,N\\}$ of sites initially
occupied by mutants. This situation is described by the set mapping
$\gamma_{S}:\\{1,\ldots,N\\}\to\\{\mathrm{M},\mathrm{R}\\}$, given by

$\gamma_{S}(i)=\begin{cases}\mathrm{M}&\text{if $i\in S$}\\\
\mathrm{R}&\text{otherwise.}\end{cases}$

Invoking the relationship between $\mathcal{A}$ and $\mathcal{M}$—in
particular, Eq. (7) and Theorem 1—we have

$\displaystyle\rho_{S}$
$\displaystyle=\lim_{t\to\infty}\operatorname{P}_{(\mathcal{M},\tilde{\gamma_{S}}(\mathbf{a}_{0}))}\,[\mathbf{s}(t)=(\mathrm{M},\ldots,\mathrm{M})]$
$\displaystyle=\lim_{t\to\infty}\operatorname{P}_{(\mathcal{A},\mathbf{a}_{0})}\,[\tilde{\gamma_{S}}(\mathbf{a}(t))=(\mathrm{M},\ldots,\mathrm{M})]$
$\displaystyle=\lim_{t\to\infty}\operatorname{P}_{(\mathcal{A},\mathbf{a}_{0})}\,[\text{$\mathbf{a}(t)=(i,\ldots,i)$
for some $i\in S$}]$ $\displaystyle=\sum_{i\in
S}\lim_{t\to\infty}\operatorname{P}_{(\mathcal{A},\mathbf{a}_{0})}\,[\mathbf{a}(t)=(i,\ldots,i)]$
$\displaystyle=\sum_{i\in S}\rho_{i}.$

This proves Eq. (10). Eq. (11) now follows from letting $S=\\{1,\ldots,N\\}$,
and noting that, in this case,

$\rho_{S}=\lim_{t\to\infty}\operatorname{P}_{(\mathcal{M},(\mathrm{M},\ldots,\mathrm{M}))}\,[\mathbf{s}(t)=(\mathrm{M},\ldots,\mathrm{M})]=1.\qed$

We now derive Eq. (1), which allows the fixation probabilities $\rho_{i}$ to
be calculated from the replacement rates $e_{ij}$.

######  Theorem 3.

For each $i=1,\ldots,N$,

$d_{i}\rho_{i}=\sum_{j=1}^{N}e_{ij}\rho_{j}.$

###### Proof.

Considering the change that can occur over a single time-step, we have the
following recurrence relation:

$\rho_{i}=(1-d_{i})\rho_{i}+\sum_{(R,\alpha)}p(R,\alpha)\,\rho_{\alpha^{-1}(i)}.$

The first term above represents the case that the occupant of $i$ survives the
current time-step and becomes the eventual ancestor of the population, while
the second term represents the case that one of $i$’s offspring from the
current time-step is the eventual ancestor. Subtracting $(1-d_{i})\rho_{i}$
from both sides and applying Theorem 2 with $S=\rho_{\alpha^{-1}(i)}$ yields

$d_{i}\rho_{i}=\sum_{(R,\alpha)}p(R,\alpha)\left(\sum_{j\in\alpha^{-1}(i)}\rho_{j}\right).$

Now interchanging the summations on the right-hand side yields

$d_{i}\rho_{i}=\sum_{j=1}^{N}\rho_{j}\left(\sum_{\begin{subarray}{c}(R,\alpha)\\\
j\in R\\\ \alpha(j)=i\end{subarray}}p(R,\alpha)\right).$

By definition,

$\sum_{\begin{subarray}{c}(R,\alpha)\\\ j\in R\\\
\alpha(j)=i\end{subarray}}p(R,\alpha)=e_{ij}.$

This completes the proof. ∎

### Proof of Result 2

######  Theorem 4 (Result 2).

The fixation probabilities from each site are equal to $1/N$ ($\rho_{i}=1/N$
for all $i=1,\ldots,N$) if and only if each site has birth rate equal to death
rate ($b_{i}=d_{i}$ for all $i=1,\ldots,N$).

###### Proof.

Assume first that the fixation probabilities from each site are all equal to
$1/N$. Substituting $\rho_{i}=1/N$ for all $i$ into Eq. (1) yields
$d_{i}=b_{i}$. This proves the “only if” direction.

Next assume that the birth rate is equal to the death rate at each site
($b_{i}=d_{i}$ for all $i=1,\ldots,N$). Eq. (1) can then be rewritten as

$\sum_{j=1}^{N}e_{ij}\rho_{i}=\sum_{j=1}^{N}e_{ij}\rho_{j}.$

Clearly, $\rho_{i}=1/N$ for all $i$ satisfies the above equation for all $i$,
and also satisfies $\sum_{i=1}^{N}\rho_{i}=1$. Assumption 1 guarantees that
the solution to these equations is unique. Therefore $b_{i}=d_{i}$ for all $i$
implies $\rho_{i}=1/N$ for all $i$, proving the “if” direction. ∎

### Proof of Result 3

######  Theorem 5 (Result 3).

If birth rates are constant over sites, $b_{i}=1/N$ for all $i=1,\ldots,N$,
then $\rho\leq 1/N$ (and consequently $K\leq u$) with equality if and only if
the death rates are also constant over sites.

###### Proof.

We separate our proof into five steps. First, we show that there is no loss of
generality in assuming that $B=1$. Second, we use the method of Lagrange
multipliers to find the critical points of the function
$\rho=\sum_{i=1}^{N}d_{i}\rho_{i}$ with respect to the variables
$\\{e_{ij}\\}_{i,j=1}^{N}$ and $\\{\rho_{i}\\}_{i=1}^{N}$ and the constraints

$\begin{cases}\displaystyle\sum_{j=1}^{N}e_{ij}=\frac{1}{N}&i\in\\{1,\ldots,N\\},\\\
\displaystyle\left(\sum_{j=1}^{N}e_{ji}\right)\rho_{i}=\sum_{j=1}^{N}e_{ij}\rho_{j}&i\in\\{1,\ldots,N\\},\\\
\displaystyle\sum_{i=1}^{N}\rho_{i}=1.\end{cases}$ (12)

We obtain that the critical points are precisely those for which $d_{i}=1/N$
for all $i$. In the next three steps, we use partial derivatives of Eqs. (1)
and (2) to form a second-order Taylor expansion of $\rho$ in the variables
$\\{e_{ij}\\}_{i\neq j}$ around the critical points. We show that the second-
order part of this expansion is $\leq 0$ with equality only at the critical
points, completing the proof.

#### Step 1: Normalize the expected number of offspring per time-step

We first show that there is no loss of generality in assuming that the total
expected number of offspring per time-step, $B\equiv\sum_{i,j}e_{ij}$, is one.
We will demonstrate this by showing that, if a maximum of $\rho$ is achieved
by a combination of replacement rates $\\{e_{ij}^{*}\\}_{i,j}$, this maximum
is also achieved using the normalized rates $\\{e_{ij}^{\prime}\\}_{i,j}$ with
$e_{ij}^{\prime}=e_{ij}^{*}/B^{*}$, $B^{*}\equiv\sum_{i,j}e_{ij}$.

First, we see that the equations

$d_{i}^{*}\rho_{i}=\sum_{j=1}^{N}e_{ij}^{*}\rho_{j}\qquad\text{and}\qquad
d_{i}^{\prime}\rho_{i}=\sum_{j=1}^{N}e_{ij}^{\prime}\rho_{j}$

are equivalent, differing only by a factor of $B$. (Above, we have introduced
the notation $d_{i}^{*}=\sum_{j=1}^{N}e_{ij}^{*}$ and
$d_{i}^{\prime}=\sum_{j=1}^{N}e_{ij}^{\prime}$.) Thus the node-specific
fixation probabilities $\\{\rho_{i}\\}_{i=1}^{N}$ resulting from
$\\{e_{ij}^{\prime}\\}_{i,j}$ are the same as those resulting from
$\\{e_{ij}^{*}\\}_{i,j}$.

Turning now to the overall fixation probabilities $\rho^{*}$ and
$\rho^{\prime}$ corresponding to $\\{e_{ij}^{*}\\}_{i,j}$ and
$\\{e_{ij}^{\prime}\\}_{i,j}$ respectively, we see that

$\rho^{*}=\frac{1}{B^{*}}\sum_{i=1}^{N}d_{i}^{*}\rho_{i}=\sum_{i=1}^{N}\frac{d_{i}^{*}}{B^{*}}\rho_{i}=\sum_{i=1}^{N}d_{i}^{\prime}\rho_{i}=\rho^{\prime}.$

Thus the same maximum is achieved by the normalized replacement rates
$\\{e_{ij}^{\prime}\\}_{i,j}$, which satisfy $\sum_{i,j}e_{ij}^{\prime}=1$. We
conclude that there is no loss of generality in assuming $B=1$. Together with
uniform birth rates, this implies the constraint

$\sum_{j=1}^{N}e_{ij}=\frac{1}{N}.$ (13)

#### Step 2: Determine critical points

We use the method of Lagrange multipliers to find critical points of the
function $\rho=\sum_{i=1}^{N}d_{i}\rho_{i}$ with respect to the variables
$\\{e_{ij}\\}_{i,j=1}^{N}$ and $\\{\rho_{i}\\}_{i=1}^{N}$. Eqs. (13), (1), and
(11) give the respective constraint equations

$\displaystyle g_{i}$ $\displaystyle\equiv\sum_{j=1}^{N}e_{ij}-\frac{1}{N}=0,\
\text{for}\,i\in\\{1,\ldots,N\\}$ $\displaystyle h_{i}$
$\displaystyle\equiv\rho_{i}\left(\sum_{j=1}^{N}e_{ji}\right)-\sum_{i=1}^{N}e_{ij}\rho_{j}=0,\
\text{for}\,i\in\\{1,\ldots,N\\}$ $\displaystyle q$
$\displaystyle\equiv\sum_{i=1}^{N}\rho_{i}-1=0$

Critical points with respect to the above constraints are solutions to

$\nabla{\rho}=\sum_{i=1}^{N}\lambda_{i}\nabla{g_{i}}+\sum_{i=1}^{N}\mu_{i}\nabla{h_{i}}+\sigma\nabla{q},$

where the $\lambda_{i}$’s, $\mu_{i}$’s, and $\sigma$ are the Lagrange
multipliers. The individual components of the gradient are obtained by taking
the partial derivative with respect to $e_{kl}$, yielding

$\rho_{l}=\lambda_{k}+\rho_{l}\mu_{l}-\rho_{l}\mu_{k},$ (14)

for $k,l\in\\{1,\ldots,N\\}$, and the partial derivative with respect to
$\rho_{k}$, yielding

$\sum_{j=1}^{N}e_{jk}=\mu_{k}\sum_{j=1}^{N}e_{jk}-\sum_{j=1}^{N}\mu_{j}e_{jk}+\sigma,$
(15)

for $k=1,\ldots,N$. Solving Eq. (14) for $\rho_{l}$ gives

$\rho_{l}=\frac{\lambda_{k}}{1-\mu_{l}+\mu_{k}}.$ (16)

Note that in the case $l=k$, Eq. (16) becomes $\rho_{k}=\lambda_{k}$.
Therefore, we replace $\lambda_{k}$ with $\rho_{k}$ in Eq. (16):

$\rho_{l}=\frac{\rho_{k}}{1-\mu_{l}+\mu_{k}}.$ (17)

Interchanging $k$ and $l$, we obtain

$\rho_{k}=\frac{\rho_{l}}{1-\mu_{k}+\mu_{l}}.$ (18)

Combining Eqs. (17) and (18) yields $\mu_{l}=\mu_{k}$ and $\rho_{l}=\rho_{k}$
for all $k,l\in\\{1,\ldots,N\\}$. That is, all node-specific probabilities
$\rho_{i}$ are equal. It follows from Eq. (11) that $\rho_{i}=1/N$ for all
$i=1,\ldots,N$. Furthermore, given that $\mu_{l}=\mu_{k}$, Eq. (15) yields
$d_{k}=\sigma$ for all $k=1,\ldots,N$. That is, all death rates $d_{i}$ are
equal. Therefore $d_{i}=1/N$ for $i=1,\ldots,N$ because we have assumed
(without loss of generality) that $\sum_{i=1}^{N}d_{i}=1$.

In summary, we have shown that the overall fixation probability $\rho$ has a
critical point whenever all node-specific fixation probabilities and death
rates are constant ($d_{i}=\rho_{i}=1/N$ for all $i=1,\ldots,N$).
Consequently, $\rho=1/N$ at all critical points. It still remains to prove
that this critical value of $\rho$ is a maximum.

#### Step 3: Taylor series expansion

To prove that the overall fixation probability has an absolute maximum of
$1/N$, we pick a critical point $\\{e_{ij}^{*}\\}_{i,j}$. Then, viewing $\rho$
as a function of the independent variables
$\\{e_{ij}\\}_{\begin{subarray}{c}i,j\\\ i\neq j\end{subarray}}$, we form a
second-order Taylor expansion of $\rho$ around the critical point. (Since we
are operating under the constraint $\sum_{j=1}^{N}e_{ij}=1/N$ for all $i$,
this set of variables suffices to determine the value of $\rho$.) We will show
that the second-order term of this Taylor expansion is always less than or
equal to zero.

The second order multivariate Taylor series expansion for $\rho$ about the
critical point $\\{e_{ij}^{*}\\}_{\begin{subarray}{c}i,j\\\ i\neq
j\end{subarray}}$ can be written as

$\rho=\frac{1}{N}+\sum_{\begin{subarray}{c}k,l\\\ k\neq
l\end{subarray}}\frac{\partial{\rho}}{\partial{e_{kl}}}\Delta{e_{kl}}+\frac{1}{2}\sum_{\begin{subarray}{c}k,l,m,n\\\
m\neq n\\\ k\neq
l\end{subarray}}\frac{\partial^{2}{\rho}}{\partial{e_{kl}\partial{e_{mn}}}}\Delta{e_{kl}}\Delta{e_{mn}}+\mathcal{O}\left(|\Delta\vec{e}|^{3}\right),$
(19)

where all derivatives are taken at
$\\{e_{ij}^{*}\\}_{\begin{subarray}{c}i,j\\\ i\neq j\end{subarray}}$ and
$\Delta e_{kl}=e_{kl}-e_{kl}^{*}$. More simply, we write this expansion as

$\rho=\frac{1}{N}+\rho^{(1)}+\rho^{(2)}+\mathcal{O}\left(|\Delta\vec{e}|^{3}\right),$

with $\rho^{(1)}$ and $\rho^{(2)}$ representing the first- and second-order
terms, respectively, in Eq. (19). The first-order term $\rho^{(1)}$ is zero
since $\\{e_{ij}^{*}\\}_{\begin{subarray}{c}i,j\\\ i\neq j\end{subarray}}$ is
a critical point. Our goal will be to show that the second-order term,
$\rho^{(2)}$ is always negative or zero.

We can find an alternative expression for $\rho^{(2)}$ using the definition of
$\rho$ in Eq. (2) as follows. We introduce the notation $\Delta
d_{i}=d_{i}-d_{i}^{*}=d_{i}-1/N$. With this notation, Eq. (2) becomes

$\displaystyle\rho$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}\rho_{i}+\sum_{i=1}^{N}\Delta
d_{i}\rho_{i}$ $\displaystyle=\frac{1}{N}+\sum_{i=1}^{N}\Delta d_{i}\rho_{i}.$
(20)

We substitute the first-order Taylor series expansion for $\rho_{i}$,

$\rho_{i}=\frac{1}{N}+\rho_{i}^{(1)}+\mathcal{O}\left(|\Delta\vec{e}|^{2}\right),$

into Eq. (20), noting that $\Delta
d_{i}=\mathcal{O}\left(|\Delta\vec{e}|\right)$ and $\sum_{i=1}^{N}\Delta
d_{i}=0$. This yields

$\displaystyle\rho$ $\displaystyle=\frac{1}{N}+\sum_{i=1}^{N}\Delta
d_{i}\left(\frac{1}{N}+\rho_{i}^{(1)}+\mathcal{O}\left(|\Delta\vec{e}|^{2}\right)\right)$
$\displaystyle=\frac{1}{N}+\sum_{i=1}^{N}\Delta
d_{i}\rho_{i}^{(1)}+\mathcal{O}\left(|\Delta\vec{e}|^{3}\right).$

Thus, the second-order term of the overall fixation probility can be written
as

$\rho^{(2)}=\sum_{i=1}^{N}\Delta d_{i}\rho_{i}^{(1)}.$ (21)

#### Step 4: Determine first-order term of the site-specific fixation
probability

We now investigate the properties of $\rho_{i}^{(1)}$. We begin by rewriting
Eq. (1) as

$d_{i}\rho_{i}=e_{ii}\rho_{i}+\sum_{j\neq i}e_{ij}\rho_{j}.$

Replacing $e_{ii}$ with $\frac{1}{N}-\sum_{j\neq i}e_{ij}$ and simplifying, we
obtain

$\Delta d_{i}\rho_{i}=\sum_{j\neq i}e_{ij}\left(\rho_{j}-\rho_{i}\right).$

Next we take the partial derivative of both sides with respect to $e_{kl}$,
where $k\neq l$,

$\rho_{i}\frac{\partial(\Delta d_{i})}{\partial{e_{kl}}}+\Delta
d_{i}\frac{\partial{\rho_{i}}}{\partial{e_{kl}}}=\sum_{\begin{subarray}{c}i,j\\\
j\neq
i\end{subarray}}\frac{\partial{e_{ij}}}{\partial{e_{kl}}}\left(\rho_{j}-\rho_{i}\right)+\sum_{\begin{subarray}{c}i,j\\\
j\neq
i\end{subarray}}e_{ij}\left(\frac{\partial{\rho_{j}}}{\partial{e_{kl}}}-\frac{\partial{\rho_{i}}}{\partial{e_{kl}}}\right),$

and evaluate at the critical point ($\rho_{i}^{*}=\frac{1}{N}$ and $\Delta
d_{i}=0$):

$\frac{1}{N}\frac{\partial(\Delta d_{i})}{\partial{e_{kl}}}=\sum_{j\neq
i}e_{ij}^{*}\left(\frac{\partial{\rho_{j}}}{\partial{e_{kl}}}-\frac{\partial{\rho_{i}}}{\partial{e_{kl}}}\right).$

Multiplying both sides by $\Delta e_{kl}$ and then summing over all
$k,l\in\\{1,\ldots,N\\}$ for $k\neq l$ yields

$\frac{1}{N}\sum_{\begin{subarray}{c}k,l\\\ k\neq
l\end{subarray}}\frac{\partial(\Delta
d_{i})}{\partial{e_{kl}}}\Delta{e_{kl}}=\sum_{\begin{subarray}{c}j,k,l\\\
k\neq l\\\ j\neq
i\end{subarray}}e_{ij}^{*}\frac{\partial{\rho_{j}}}{\partial{e_{kl}}}\Delta{e_{kl}}-\sum_{\begin{subarray}{c}j,k,l\\\
k\neq l\\\ j\neq
i\end{subarray}}e_{ij}^{*}\frac{\partial{\rho_{i}}}{\partial{e_{kl}}}\Delta{e_{kl}}.$
(22)

We observe that $\rho_{i}^{(1)}=\sum_{\begin{subarray}{c}k,l\\\ k\neq
l\end{subarray}}\frac{\partial{\rho_{i}}}{\partial{e_{kl}}}\Delta{e_{kl}}$ is
the first order term of the Taylor series expansion of $\rho_{i}$ about the
critical point, and $\sum_{\begin{subarray}{c}k,l\\\ k\neq
l\end{subarray}}\frac{\partial(\Delta
d_{i})}{\partial{e_{kl}}}\Delta{e_{kl}}=\Delta d_{i}$, since $\Delta d_{i}$ is
a linear function of the $e_{kl}$. Thus Eq. (22) becomes

$\frac{1}{N}\Delta d_{i}=\sum_{j\neq i}e_{ij}^{*}\rho_{j}^{(1)}-\sum_{j\neq
i}e_{ij}^{*}\rho_{i}^{(1)}.$

The restriction $j\neq i$ in the sums on the right-hand side can be removed,
since the entire right-hand side is zero when $j=i$. Rearranging further, we
have

$N\rho_{i}^{(1)}\sum_{j=1}^{N}e^{*}_{ij}=-\Delta
d_{i}+N\sum_{j=1}^{N}e^{*}_{ij}\rho_{j}^{(1)}.$

Since the birth rate is uniform over sites
$(b_{i}=\sum_{j=1}^{N}e_{ij}^{*}=\frac{1}{N}$ for all $i$), the above equation
simplifies to

$\rho_{i}^{(1)}=-\Delta d_{i}+N\sum_{j=1}^{N}e_{ij}^{*}\rho_{j}^{(1)}.$ (23)

This equation provides a recurrence relation for the first-order term of the
Taylor expansion of $\rho_{i}$ about the critical point.

#### Step 5: Determine second-order term of the overall fixation probability

Finally, we combine Eq. (23) with Eq. (21) to show that $\rho^{(2)}\leq 0$.
The process is as follows. We first multiply both sides of Eq. (23) by
$\rho_{i}^{(1)}$ and sum over $i=1,\ldots,N$:

$\sum_{i=1}^{N}\left(\rho_{i}^{(1)}\right)^{2}=-\sum_{i=1}^{N}\Delta
d_{i}\rho_{i}^{(1)}+N\sum_{i,j}e_{ij}^{*}\rho_{i}^{(1)}\rho_{j}^{(1)}.$

By Eq. (21) we can substitute $\rho^{(2)}$ for $\sum_{i=1}^{N}\Delta
d_{i}\rho_{i}^{(1)}$. Then, solving for $\rho^{(2)}$, we obtain

$\displaystyle\rho^{(2)}=-\sum_{i=1}^{N}\left(\rho_{i}^{(1)}\right)^{2}+N\sum_{i,j}e_{ij}^{*}\rho_{i}^{(1)}\rho_{j}^{(1)}.$
(24)

We now make use of the fact that the product $\rho_{i}^{(1)}\rho_{j}^{(1)}$
can be written as a difference of squares:

$\rho_{i}^{(1)}\rho_{j}^{(1)}=-\frac{1}{2}\left[\left(\rho_{i}^{(1)}-\rho_{j}^{(1)}\right)^{2}-\left(\rho_{i}^{(1)}\right)^{2}-\left(\rho_{j}^{(1)}\right)^{2}\right].$

Substituting this identity into Eq. (24) yields

$\displaystyle\rho^{(2)}$
$\displaystyle=-\sum_{i=1}^{N}\left(\rho_{i}^{(1)}\right)^{2}-\frac{N}{2}\sum_{i,j}e_{ij}^{*}\left(\rho_{i}^{(1)}-\rho_{j}^{(1)}\right)^{2}+\frac{N}{2}\sum_{i,j}e_{ij}^{*}\left(\rho_{i}^{(1)}\right)^{2}$
$\displaystyle\qquad+\frac{N}{2}\sum_{i,j}e_{ij}^{*}\left(\rho_{j}^{(1)}\right)^{2}$
$\displaystyle=-\sum_{i=1}^{N}\left(\rho_{i}^{(1)}\right)^{2}-\frac{N}{2}\sum_{i,j}e_{ij}^{*}\left(\rho_{i}^{(1)}-\rho_{j}^{(1)}\right)^{2}+\frac{N}{2}\sum_{i=1}^{N}\left(\rho_{i}^{(1)}\right)^{2}\sum_{j=1}^{N}e_{ij}^{*}$
$\displaystyle\qquad+\frac{N}{2}\sum_{j=1}^{N}\left(\rho_{j}^{(1)}\right)^{2}\sum_{i=1}^{N}e_{ij}^{*}.$

Uniform birth rates implies that $\sum_{j=1}^{N}e^{*}_{ij}=1/N$ for each $i$.
Furthermore, since $\\{e_{ij}^{*}\\}_{i,j}$ is a critical point, death rates
are uniform as well; thus $\sum_{i=1}^{N}e^{*}_{ij}=1/N$ for each $j$. Making
these substitutions yields

$\displaystyle\rho^{(2)}$
$\displaystyle=-\sum_{i=1}^{N}\left(\rho_{i}^{(1)}\right)^{2}-\frac{N}{2}\sum_{i,j}e_{ij}^{*}\left(\rho_{i}^{(1)}-\rho_{j}^{(1)}\right)^{2}+\frac{1}{2}\sum_{i=1}^{N}\left(\rho_{i}^{(1)}\right)^{2}+\frac{1}{2}\sum_{j=1}^{N}\left(\rho_{j}^{(1)}\right)^{2}$
$\displaystyle=-\frac{N}{2}\sum_{i,j}e_{ij}^{*}\left(\rho_{i}^{(1)}-\rho_{j}^{(1)}\right)^{2}.$

We conclude that $\rho^{(2)}$ is always less than or equal to $0$.

Overall, we have shown that the critical points of $\rho$, subject to the
constraint $b_{i}=1/N$ for all $i$, are precisely those for which $d_{i}=1/N$
for all $i$, and that $\rho$ achieves its maximum value of $1/N$ at these
points. ∎

### Result 4

We now turn to the full range of possible values for $K$ with no constraints
on birth and death rates. Consider the example spatial structure illustrated
in Figure 5, which consists of a hub with outgoing edges to $n$ leaves, so
that the population size is $N=n+1$. There is an edge of weight $a$, $0\leq
a<1/(N-1)$, from the hub to each leaf. The hub also has a self-loop of weight
$1-(N-1)a$, so that $B=1$ births are expected per time-step. The death rates
are $d_{\mathrm{H}}=1-(N-1)a$ for the hub and $d_{\mathrm{L}}=a$ for each
leaf. Solving Eq. (1) we obtain the site-specific fixation probabilities
$\rho_{\mathrm{L}}=0$ for each leaf and $\rho_{\mathrm{H}}=1$ for the hub. By
Eq. (2), the overall fixation probability is equal to the rate of turnover at
the hub:

$\rho=d_{\mathrm{H}}=1-(N-1)a.$

Any value $0\leq\rho<1$ can be obtained by an appropriate choice of $a$, with
$0<a\leq 1/(N-1)$, specifically, $a=(1-\rho)/(N-1)$. This proves:

######  Theorem 6 (Result 4).

For arbitrary spatial population structure (no constraints on $e_{ij}$) the
fixation probability can take any value $0\leq\rho<1$, and consequently, the
molecular clock can take any rate $0\leq K<Nu$.

We observe that if $a=1/N$ then the death rate is $1/N$ at each site, and
therefore $\rho=1/N$ by Result 1. If $a=1/(N-1)$ then there is no turnover at
the hub and thus $\rho=0$. At the other extreme, as $a$ approach zero, $\rho$
comes arbitrarily close to unity.

### Exact expressions for $\rho$ with $N\leq 3$

To complement the above arguments, which apply to general population size $N$,
we here derive exact expressions for the overall fixation probability $\rho$
in the cases $N=2$ and $N=3$. Our results for node-specific fixation
probabilities coincide with those found previous for a different model of
evolution in a subdivided population [39].

#### Population size $N=2$

We solve Eq. (1) to obtain expressions for $\rho_{1}$ and $\rho_{2}$:

$\displaystyle\rho_{1}$ $\displaystyle=\frac{e_{12}}{e_{12}+e_{21}},$ (25)
$\displaystyle\rho_{2}$ $\displaystyle=\frac{e_{21}}{e_{12}+e_{21}}.$ (26)

We first consider uniform birth rates
($b_{i}=\sum_{j=1}^{2}e_{ij}=\frac{1}{2}$ for $i=1,2$). Substituting Eqs. (25)
and (26) into Eq. (2) ($B=1$ in this case), gives an expression for the
overall fixation probability:

$\rho=\frac{1}{2}-\frac{(e_{12}-e_{21})^{2}}{e_{12}+e_{21}}.$

Thus, in accordance with Result 3, $\rho$ is less than or equal to
$\frac{1}{2}$ and consequently $K\leq u$. Equality occurs if and only if
$e_{12}=e_{21}$, which corresponds to equal node-specific fixation
probabilities ($\rho_{1}=\rho_{2}=1/N$ from Eqs. (25) and (26)) and equal
death rates ($d_{1}=d_{2}=1/N$).

When the condition of uniform birth rates is relaxed, we find, by substituting
Eqs. (25) and (26) into Eq. (2), the following expression for the overall
fixation probability:

$\rho=\frac{1}{2}-\frac{(e_{12}-e_{21})\left(\frac{B}{2}-d_{1}\right)}{B(e_{12}+e_{21})}.$

We discover that $\rho=\frac{1}{2}$, and consequntly $K=u$, if (i) death rates
are equal, $d_{1}=d_{2}=\frac{B}{2}$ (Result 1), or (ii) node-specific
fixation probabilities are equal, $\rho_{1}=\rho_{2}=\frac{1}{2}$ (Result 2).
When death rates are unequal and node-specific fixation probabilities are
unequal, there are cases in which $\rho<\frac{1}{2}$ ($K<u)$ and cases in
which $\rho>\frac{1}{2}$ ($K>u$). In particular, suppose that
$d_{1}>\frac{B}{2}$. Then $\rho$ is less than $\frac{1}{2}$ if $e_{12}<e_{21}$
and greater than $\frac{1}{2}$ if $e_{12}>e_{21}$.

#### Population size $N=3$

We solve Eqs. (1) and (11) to obtain expressions for $\rho_{1}$, $\rho_{2}$
and $\rho_{3}$:

$\displaystyle\rho_{1}$
$\displaystyle=\frac{e_{12}e_{13}+e_{12}e_{23}+e_{13}e_{32}}{e_{12}\left(e_{13}+e_{23}+e_{31}\right)+e_{13}\left(e_{21}+e_{32}\right)+\left(e_{21}+e_{31}\right)\left(e_{23}+e_{32}\right)},$
(27) $\displaystyle\rho_{2}$
$\displaystyle=\frac{e_{21}e_{13}+e_{21}e_{23}+e_{23}e_{31}}{e_{12}\left(e_{13}+e_{23}+e_{31}\right)+e_{13}\left(e_{21}+e_{32}\right)+\left(e_{21}+e_{31}\right)\left(e_{23}+e_{32}\right)},$
(28) $\displaystyle\rho_{3}$
$\displaystyle=\frac{e_{12}e_{31}+e_{21}e_{32}+e_{31}e_{32}}{e_{12}\left(e_{13}+e_{23}+e_{31}\right)+e_{13}\left(e_{21}+e_{32}\right)+\left(e_{21}+e_{31}\right)\left(e_{23}+e_{32}\right)}.$
(29)

Substituting Eqs. (27)-(29) into Eq. (2) gives an expression for the overall
fixation probability $\rho=\frac{\text{num}}{\text{denom}}$ where

num
$\displaystyle=d_{1}\left(e_{12}e_{13}+e_{12}e_{23}+e_{13}e_{32}\right)+d_{2}\left(e_{21}e_{13}+e_{21}e_{23}+e_{23}e_{31}\right)$
$\displaystyle\quad+d_{3}\left(e_{12}e_{31}+e_{21}e_{32}+e_{31}e_{32}\right)$
denom
$\displaystyle=B\left(e_{12}\left(e_{13}+e_{23}+e_{31}\right)+e_{13}\left(e_{21}+e_{32}\right)+\left(e_{21}+e_{31}\right)\left(e_{23}+e_{32}\right)\right).$

We factor to obtain an expression for $\text{num}(\Delta\rho)$ in terms of
$\tilde{\Delta}d_{i}=d_{i}-\frac{B}{3}$ for $i\in\\{1,2,3\\}$,

$\displaystyle\text{num}(\Delta\rho)$
$\displaystyle=(e_{12}\tilde{\Delta}d_{1}+e_{21}\tilde{\Delta}d_{2})(d_{3}-b_{3})-e_{32}\tilde{\Delta}d_{1}(d_{1}-b_{1})-e_{31}\tilde{\Delta}d_{2}(d_{2}-b_{2})$
(30)

From Eq. (30) we see that $\rho=\frac{1}{3}$, and conseqently $K=u$, in the
case of equal death rates ($d_{1}=d_{2}=d_{3}=\frac{B}{3}$) and in the case of
equivalent birth and death rate at each site ($b_{i}=d_{i}$ for
$i\in\\{1,2,3\\}$). This conforms to Result 1 and 2. When death rates are
nonuniform and node-specific probabilities are also nonuniform, we find cases
in which $K>u$ and $K<u$. For example, if
$\tilde{\Delta}d_{1}=0,\tilde{\Delta}d_{2}>0,d_{2}>b_{2},$ and $d_{3}<b_{3}$
then Eq. (30) yields $\text{num}(\Delta\rho)<0$ and consequently,
$\rho<\frac{1}{3}$ and $K<u$. On the other hand, if
$\tilde{\Delta}d_{1}=0,b_{2}>d_{2}>\frac{B}{3},$ and $d_{3}>b_{3}$ then Eq.
(30) yields $\text{num}(\Delta\rho)>0$ and therefore, $\rho>\frac{1}{3}$ and
$K>u$.

### Upstream-downstream populations

We now turn to the upstream-downstream model introduced in the Results and in
Figure 6.

######  Theorem 7.

In the upstream-downstream model, $\rho>1/N$, and consequently $K>u$, if and
only if $d_{\uparrow}>d_{\downarrow}$.

###### Proof.

From Eq. (2) we obtain the following expression for the overall fixation
probability

$\displaystyle\rho$
$\displaystyle=\frac{1}{B}\left(N_{\uparrow}d_{\uparrow}\rho_{\uparrow}+N_{\downarrow}d_{\downarrow}\rho_{\downarrow}\right)$
$\displaystyle=\frac{1}{N}+\frac{1}{B}\left[N_{\uparrow}d_{\uparrow}\left(\rho_{\uparrow}-\frac{1}{N}\right)+N_{\downarrow}d_{\downarrow}\left(\rho_{\downarrow}-\frac{1}{N}\right)\right]$
(31)

Since fixation probabilities sum to one, and since the total population size
is $N=N_{\uparrow}+N_{\downarrow}$, we have

$N_{\downarrow}\left(\rho_{\downarrow}-\frac{1}{N}\right)=-N_{\uparrow}\left(\rho_{\uparrow}-\frac{1}{N}\right).$
(32)

Substituting in Eq. (31) yields

$\rho=\frac{1}{N}+(d_{\uparrow}-d_{\downarrow})\frac{N_{\uparrow}}{B}\left(\rho_{\uparrow}-\frac{1}{N}\right).$
(33)

It follows from Eq. (4) and from $e_{\rightarrow}>e_{\leftarrow}$ that
$\rho_{\uparrow}>\rho_{\downarrow}$. Moreover, Eq. (32) implies that
$\rho_{\uparrow}-1/N$ and $\rho_{\downarrow}-1/N$ have opposite signs, and
since $\rho_{\uparrow}>\rho_{\downarrow}$, it follows that
$\rho_{\uparrow}-1/N$ must be positive. Thus the second term on the right-hand
side of Eq. (33) has the sign of $d_{\uparrow}-d_{\downarrow}$. We therefore
conclude that $\rho>1/N$, and consequently the molecular clock is accelerated
relative to the well-mixed case ($K>u$), if and only if
$d_{\uparrow}>d_{\downarrow}$. ∎

## Acknowledgements

The Foundational Questions in Evolutionary Biology initiative at Harvard
University is sponsored by the John Templeton Foundation (RFP-12-02).

## References

  * 1. Zuckerkandl E, Pauling L (1965) Evolutionary divergence and convergence in proteins.  Evolving genes and proteins 97: 97–166.
  * 2. Kimura M (1968) Evolutionary rate at the molecular level.  Nature 217: 624–626.
  * 3. Ayala FJ (1999) Molecular clock mirages.  BioEssays 21: 71-75.
  * 4. Bromham L, Penny D (2003) The modern molecular clock.  Nature Reviews Genetics 4: 216–224.
  * 5. Pollak E (1982) The rate of mutant substitution in populations with overlapping generations.  Genetical Research 40: 89–94.
  * 6. Balloux F, Lehmann L (2012) Substitution rates at neutral genes depend on population size under fluctuating demography and overlapping generations.  Evolution 66: 605–611.
  * 7. Lehmann L (2014) Stochastic demography and the neutral substitution rate in class-structured populations.  Genetics 197: 351–360.
  * 8. Ohta T, Kimura M (1971) On the constancy of the evolutionary rate of cistrons.  Journal of Molecular Evolution 1: 18–25.
  * 9. Kimura M, Ohta T (1974) On some principles governing molecular evolution.  Proceedings of the National Academy of Sciences 71: 2848–2852.
  * 10. Langley CH, Fitch WM (1974) An examination of the constancy of the rate of molecular evolution.  Journal of Molecular Evolution 3: 161–177.
  * 11. Kimura M (1984) The Neutral Theory of Molecular Evolution.  Cambridge University Press.
  * 12. Ayala FJ (1986) On the virtues and pitfalls of the molecular evolutionary clock.  Journal of Heredity 77: 226–235.
  * 13. Li WH, Tanimura M, Sharp PM (1987) An evaluation of the molecular clock hypothesis using mammalian dna sequences.  Journal of Molecular Evolution 25: 330–342.
  * 14. Kumar S (2005) Molecular clocks: four decades of evolution.  Nature Reviews Genetics 6: 654–662.
  * 15. Maruyama T (1970) On the fixation probability of mutant genes in a subdivided population.  Genetical Research 15: 221–225.
  * 16. Lande R (1979) Effective deme sizes during long-term evolution estimated from rates of chromosomal rearrangement.  Evolution 33: 234–251.
  * 17. Slatkin M (1981) Fixation probabilities and fixation times in a subdivided population.  Evolution : 477–488.
  * 18. Nagylaki T (1982) Geographical invariance in population genetics.  Journal of Theoretical Biology 99: 159–172.
  * 19. Tachida H, Iizuka M (1991) Fixation probability in spatially changing environments.  Genetical Research 58: 243–251.
  * 20. Barton N (1993) The probability of fixation of a favoured allele in a subdivided population.  Genetical Research 62: 149–149.
  * 21. Roze D, Rousset F (2003) 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.
  * 22. Whitlock MC (2003) Fixation probability and time in subdivided populations.  Genetics 164: 767–779.
  * 23. Lieberman E, Hauert C, Nowak MA (2005) Evolutionary dynamics on graphs.  Nature 433: 312–316.
  * 24. Broom M, Rychtář J (2008) An analysis of the fixation probability of a mutant on special classes of non-directed graphs.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Science 464: 2609–2627.
  * 25. Patwa Z, Wahl LM (2008) The fixation probability of beneficial mutations.  Journal of The Royal Society Interface 5: 1279–1289.
  * 26. Masuda N, Ohtsuki H (2009) Evolutionary dynamics and fixation probabilities in directed networks.  New Journal of Physics 11: 033012.
  * 27. Houchmandzadeh B, Vallade M (2011) The fixation probability of a beneficial mutation in a geographically structured population.  New Journal of Physics 13: 073020.
  * 28. Shakarian P, Roos P, Johnson A (2012) A review of evolutionary graph theory with applications to game theory.  Biosystems 107: 66–80.
  * 29. Voorhees B (2013) Birth–death fixation probabilities for structured populations.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Science 469: 20120248.
  * 30. Constable GWA, McKane AJ (2014) Population genetics on islands connected by an arbitrary network: An analytic approach.  Journal of Theoretical Biology 358: 149–165.
  * 31. Monk T, Green P, Paulin M (2014) Martingales and fixation probabilities of evolutionary graphs.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Science 470: 20130730.
  * 32. Allen B, Tarnita CE (2014) Measures of success in a class of evolutionary models with fixed population size and structure.  Journal of Mathematical Biology 68: 109–143.
  * 33. Liggett TM (1985) Particle Systems.  Springer.
  * 34. Nowak MA, Bonhoeffer S, May RM (1994) More spatial games.  International Journal of Bifurcation and Chaos 4: 33–56.
  * 35. Antal T, Redner S, Sood V (2006) Evolutionary dynamics on degree-heterogeneous graphs.  Physical Review Letters 96: 188104.
  * 36. Nowak MA, May RM (1992) Evolutionary games and spatial chaos.  Nature 359: 826–829.
  * 37. Taylor P, Lillicrap T, Cownden D (2011) Inclusive fitness analysis on mathematical groups.  Evolution 65: 849–859.
  * 38. Maciejewski W (2014) Reproductive value in graph-structured populations.  Journal of Theoretical Biology 340: 285–293.
  * 39. Lundy IJ, Possingham HP (1998) Fixation probability of an allele in a subdivided population with asymmetric migration.  Genetical Research 71: 237–245.
  * 40. Kozar S, Morrissey E, Nicholson AM, van der Heijden M, Zecchini HI, et al. (2013) Continuous clonal labeling reveals small numbers of functional stem cells in intestinal crypts and adenomas.  Cell Stem Cell 13: 626–633.
  * 41. Lopez-Garcia C, Klein AM, Simons BD, Winton DJ (2010) Intestinal stem cell replacement follows a pattern of neutral drift.  Science 330: 822–825.
  * 42. Snippert HJ, van der Flier LG, Sato T, van Es JH, van den Born M, et al. (2010) Intestinal crypt homeostasis results from neutral competition between symmetrically dividing lgr5 stem cells.  Cell 143: 134–144.
  * 43. Vermeulen L, Snippert HJ (2014) Stem cell dynamics in homeostasis and cancer of the intestine.  Nature Reviews Cancer 14: 468–480.
  * 44. Vermeulen L, Morrissey E, van der Heijden M, Nicholson AM, Sottoriva A, et al. (2013) Defining stem cell dynamics in models of intestinal tumor initiation.  Science 342: 995–998.
  * 45. Barker N, van Oudenaarden A, Clevers H (2012) Identifying the stem cell of the intestinal crypt: strategies and pitfalls.  Cell Stem Cell 11: 452–460.
  * 46. Leskovec J (2011).  Stanford large network dataset collection.  URL http://snap.stanford.edu/data/index.html.
  * 47. Durrett R, Levin S (1994) The importance of being discrete (and spatial).  Theoretical Population Biology 46: 363–394.
  * 48. Nakamaru M, Matsuda H, Iwasa Y (1997) The evolution of cooperation in a lattice-structured population.  Journal of Theoretical Biology 184: 65–81.
  * 49. van Baalen M, Rand DA (1998) The unit of selection in viscous populations and the evolution of altruism.  Journal of Theoretical Biology 193: 631–648.
  * 50. Ohtsuki H, Hauert C, Lieberman E, Nowak MA (2006) A simple rule for the evolution of cooperation on graphs and social networks.  Nature 441: 502–505.
  * 51. Nowak MA, Tarnita CE, Antal T (2010) Evolutionary dynamics in structured populations.  Philosophical Transactions of the Royal Society B: Biological Sciences 365: 19–30.
  * 52. Allen B, Nowak MA (2014) Games on graphs.  EMS Surveys in Mathematical Sciences 1: 113–151.
  * 53. Débarre F, Hauert C, Doebeli M (2014) Social evolution in structured populations.  Nature Communications 5: 4409.
  * 54. Felsenstein J (1976) The theoretical population genetics of variable selection and migration.  Annual Review of Genetics 10: 253–280.
  * 55. Hedrick PW, Ginevan ME, Ewing EP (1976) Genetic polymorphism in heterogeneous environments.  Annual Review of Ecology and Systematics 7: 1–32.
  * 56. Slatkin M (1987) Gene flow and the geographic structure of natural populations.  Science 236: 787–792.
  * 57. Cherry JL, Wakeley J (2003) A diffusion approximation for selection and drift in a subdivided population.  Genetics 163: 421–428.
  * 58. MacArthur RH, Wilson EO (1967) The Theory of Island Biogeography.  Princeton, NJ: Princeton University Press.
  * 59. Doebeli M, Dieckmann U (2003) Speciation along environmental gradients.  Nature 421: 259–264.
  * 60. De Aguiar M, Baranger M, Baptestini E, Kaufman L, Bar-Yam Y (2009) Global patterns of speciation and diversity.  Nature 460: 384–387.
  * 61. Kawecki TJ, Holt RD (2002) Evolutionary consequences of asymmetric dispersal rates.  The American Naturalist 160: 333–347.
  * 62. Lenormand T (2002) Gene flow and the limits to natural selection.  Trends in Ecology & Evolution 17: 183–189.
  * 63. Garant D, Forde SE, Hendry AP (2007) The multifarious effects of dispersal and gene flow on contemporary adaptation.  Functional Ecology 21: 434–443.
  * 64. Nagylaki T (1978) Clines with asymmetric migration.  Genetics 88: 813–827.
  * 65. Vuilleumier S, Wilcox C, Cairns BJ, Possingham HP (2007) How patch configuration affects the impact of disturbances on metapopulation persistence.  Theoretical Population Biology 72: 77–85.
  * 66. Morrissey MB, de Kerckhove DT (2009) The maintenance of genetic variation due to asymmetric gene flow in dendritic metapopulations.  The American Naturalist 174: 875–889.
  * 67. Tackenberg O, Poschlod P, Bonn S (2003) Assessment of wind dispersal potential in plant species.  Ecological Monographs 73: 191–205.
  * 68. Muñoz J, Felicísimo ÁM, Cabezas F, Burgaz AR, Martínez I (2004) Wind as a long-distance dispersal vehicle in the southern hemisphere.  Science 304: 1144–1147.
  * 69. Collins CJ, Fraser CI, Ashcroft A, Waters JM (2010) Asymmetric dispersal of southern bull-kelp (durvillaea antarctica) adults in coastal new zealand: testing an oceanographic hypothesis.  Molecular Ecology 19: 4572–4580.
  * 70. Pringle JM, Blakeslee AM, Byers JE, Roman J (2011) Asymmetric dispersal allows an upstream region to control population structure throughout a speciesí range.  Proceedings of the National Academy of Sciences 108: 15288–15293.
  * 71. Angert AL (2009) The niche, limits to species’ distributions, and spatiotemporal variation in demography across the elevation ranges of two monkeyflowers.  Proceedings of the National Academy of Sciences 106: 19693–19698.
  * 72. O’keefe K, Ramakrishnan U, Van Tuinen M, Hadly EA (2009) Source–sink dynamics structure a common montane mammal.  Molecular Ecology 18: 4775–4789.
  * 73. Nowak MA, Michor F, Iwasa Y (2003) The linear process of somatic evolution.  Proceedings of the National Academy of Sciences 100: 14966–14969.
  * 74. Dingli D, Traulsen A, Pacheco JM (2007) Stochastic dynamics of hematopoietic tumor stem cells.  Cell Cycle 6: 461–466.
  * 75. Bozic I, Nowak MA (2013) Unwanted evolution.  Science 342: 938–939.
  * 76. Snippert HJ, Schepers AG, van Es JH, Simons BD, Clevers H (2014) Biased competition between lgr5 intestinal stem cells driven by oncogenic mutation induces clonal expansion.  EMBO Rep 15: 62–69.
  * 77. Knudson AG (2001) Two genetic hits (more or less) to cancer.  Nature Reviews Cancer 1: 157–162.
  * 78. Attolini CSO, Michor F (2009) Evolutionary theory of cancer.  Annals of the New York Academy of Sciences 1168: 23–51.
  * 79. Vogelstein B, Papadopoulos N, Velculescu VE, Zhou S, Diaz LA, et al. (2013) Cancer genome landscapes.  Science 339: 1546–1558.
  * 80. Goldie JH, Coldman AJ (1998) Drug resistance in cancer: mechanisms and models.  Cambridge University Press.
  * 81. Komarova NL, Wodarz D (2005) Drug resistance in cancer: principles of emergence and prevention.  Proceedings of the National Academy of Sciences 102: 9714–9719.
  * 82. Michor F, Nowak MA, Iwasa Y (2006) Evolution of resistance to cancer therapy.  Current Pharmaceutical Design 12: 261–271.
  * 83. Bozic I, Reiter JG, Allen B, Antal T, Chatterjee K, et al. (2013) Evolutionary dynamics of cancer in response to targeted combination therapy.  Elife 2.
  * 84. Gupta PB, Fillmore CM, Jiang G, Shapira SD, Tao K, et al. (2011) Stochastic state transitions give rise to phenotypic equilibrium in populations of cancer cells.  Cell 146: 633–644.
  * 85. May RM, Lloyd AL (2001) Infection dynamics on scale-free networks.  Physical Review E 64: 066112.
  * 86. Watts DJ (2002) A simple model of global cascades on random networks.  Proceedings of the National Academy of Sciences 99: 5766–5771.
  * 87. Christakis NA, Fowler JH (2007) The spread of obesity in a large social network over 32 years.  New England Journal of Medicine 357: 370–379.
  * 88. Rosvall M, Bergstrom CT (2008) Maps of random walks on complex networks reveal community structure.  Proceedings of the National Academy of Sciences 105: 1118–1123.
  * 89. Castellano C, Fortunato S, Loreto V (2009) Statistical physics of social dynamics.  Reviews of Modern Physics 81: 591.
  * 90. Centola D (2010) The spread of behavior in an online social network experiment.  Science 329: 1194–1197.
  * 91. Hill AL, Rand DG, Nowak MA, Christakis NA (2010) Infectious disease modeling of social contagion in networks.  PLoS Computational Biology 6: e1000968.
  * 92. Bakshy E, Rosenn I, Marlow C, Adamic L (2012) The role of social networks in information diffusion.  In: Proceedings of the 21st international conference on World Wide Web. ACM, pp. 519–528.
  * 93. Moreno Y, Pastor-Satorras R, Vespignani A (2002) Epidemic outbreaks in complex heterogeneous networks.  The European Physical Journal B-Condensed Matter and Complex Systems 26: 521–529.
  * 94. Nekovee M, Moreno Y, Bianconi G, Marsili M (2007) Theory of rumour spreading in complex social networks.  Physica A: Statistical Mechanics and its Applications 374: 457–470.
  * 95. Skyrms B, Pemantle R (2009) A dynamic model of social network formation.  In: Adaptive Networks, Springer. pp. 231–251.
  * 96. Rosvall M, Bergstrom CT (2010) Mapping change in large networks.  PloS ONE 5: e8694.
  * 97. Rand DG, Arbesman S, Christakis NA (2011) Dynamic social networks promote cooperation in experiments with humans.  Proceedings of the National Academy of Sciences 108: 19193–19198.
  * 98. Wardil L, Hauert C (2014) Origin and structure of dynamic cooperative networks.  Scientific Reports 4.
  * 99. Baxter GJ, Blythe RA, McKane AJ (2008) Fixation and consensus times on a network: A unified approach.  Physical Review Letters 101: 258701.
  * 100. Frean M, Rainey PB, Traulsen A (2013) The effect of population structure on the rate of evolution.  Proceedings of the Royal Society B: Biological Sciences 280.
  * 101. Maciejewski W, Fu F, Hauert C (2014) Evolutionary game dynamics in populations with heterogenous structures.  PLoS Computational Biology 10: e1003567.
  * 102. Tarnita CE, Taylor PD (2014) Measures of relative fitness of social behaviors in finite structured population models.  American Naturalist 184 (published online ahead of print).
  * 103. Broom M, Hadjichrysanthou C, Rychtář J, Stadler B (2010) Two results on evolutionary processes on general non-directed graphs.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Science 466: 2795–2798.
  * 104. Shakarian P, Roos P, Moores G (2013) A novel analytical method for evolutionary graph theory problems.  Biosystems 111: 136–144.