RESEARCH ARTICLE The Molecular Clock of Neutral Evolution Can Be Accelerated or Slowed by Asymmetric Spatial Structure Benjamin Allen1,2,3*, Christine Sample1, Yulia Dementieva1, Ruben C. Medeiros1, Christopher Paoletti1, Martin A. Nowak2,4 1 Department of Mathematics, Emmanuel College, Boston, Massachusetts, United States of America, 2 Program for Evolutionary Dynamics, Harvard University, Cambridge, Massachusetts, United States of America, 3 Center for Mathematical Sciences and Applications, Harvard University, Cambridge, Massachusetts, United States of America, 4 Department of Mathematics, Department of Organismic and Evolutionary Biology, Harvard University, Cambridge, Massachusetts, United States of America * benjcallen@gmail.com a11111 Abstract OPEN ACCESS Citation: Allen B, Sample C, Dementieva Y, Medeiros RC, Paoletti C, Nowak MA (2015) The Molecular Clock of Neutral Evolution Can Be Accelerated or Slowed by Asymmetric Spatial Structure. PLoS Comput Biol 11(2): e1004108. doi:10.1371/journal.pcbi.1004108 Editor: Carl T. Bergstrom, University of Washington, UNITED STATES Received: September 11, 2014 Accepted: January 2, 2015 Published: February 26, 2015 Copyright: © 2015 Allen et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Data Availability Statement: All relevant data are within the paper. Funding: BA and MAN are supported by the Foundational Questions in Evolutionary Biology initiative of the John Templeton Foundation, RFP-1202. http://www.templeton.org/node/1521. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Competing Interests: The authors have declared that no competing interests exist. 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 Ku. Our results apply to a general class of haploid, asexually reproducing, spatially structured populations. 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 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 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 PLOS Computational Biology | DOI:10.1371/journal.pcbi.1004108 February 26, 2015 1 / 32 Asymmetric Spatial Structure Alters the Molecular Clock Rate 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 single 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 ρ—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 [3, 4, 8–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–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ρ 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 ρ = 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 Fig. 1): PLOS Computational Biology | DOI:10.1371/journal.pcbi.1004108 February 26, 2015 2 / 32 Asymmetric Spatial Structure Alters the Molecular Clock Rate Fig 1. How spatial structure affects 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 bi and death rates di 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 u; Result 3). doi:10.1371/journal.pcbi.1004108.g001 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 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 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, . . ., 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. Replacement events are chosen according to a fixed probability distribution—called the replacement rule—specific to the model in question. Since we consider only neutral mutations that have no phenotypic effect, the probabilities of replacement events do not depend on the current population state. PLOS Computational Biology | DOI:10.1371/journal.pcbi.1004108 February 26, 2015 3 / 32 Asymmetric Spatial Structure Alters the Molecular Clock Rate This class includes many established evolutionary models. One important subclass is spatial Moran processes [23, 33–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 R & {1, . . ., N} of individuals may be replaced in a given replacement event. Parentage in a replacement event is recorded in an offspring-to-parent map α:R ! {1, . . ., N} (see Methods, [32, 38]), which ensures that each offspring has exactly one parent and allows us to trace lineages over time. For a given model in this class, we let eij denote the (marginal) 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 PN offspring of site i over a single time-step is bi ¼ j¼1 eij . The probability that node i dies (i.e., is PN replaced) in a time-step is di ¼ j¼1 eji . The death rate di can also be regarded as the rate of turnover at site i. The total expected number of offspring per time-step is denoted PN P PN B ¼ i¼1 bi ¼ i¼1 di ¼ i;j eij . 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 ρi denote the probability that a new mutation arising at site i becomes fixed. (ρi can also be understood as the reproductive value of site i [39].) We show in the Methods that the fixation probabilities ρi are the unique solution to the system of equations di ri ¼ N X eij rj j¼1 for i ¼ 1; . . . ; N; ð1Þ N X ri ¼ 1: i¼1 ð2Þ Equation (2) arises because ρ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 di, because each new offspring provides an independent chance of mutation. Specifically, if mutation occurs at rate u ( 1 per reproduction, then new mutations arise at site i at rate Nudi/B per generation [32]. Thus the fraction of mutations that arise at site i is di/B. The overall fixation probability ρ of new mutations, taking into account all possible initial sites, is therefore r¼ N 1X dr: B i¼1 i i ð3Þ The molecular clock rate K is obtained by multiplying the ﬁxation probability ρ by the total rate of mutation per generation: K ¼ Nur ¼ N Nu X dr: B i¼1 i i ð4Þ 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 formula is PN K ¼ Bur ¼ u i¼1 di ri . PLOS Computational Biology | DOI:10.1371/journal.pcbi.1004108 February 26, 2015 4 / 32 Asymmetric Spatial Structure Alters the Molecular Clock Rate Fig 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. doi:10.1371/journal.pcbi.1004108.g002 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 eij is constant over all i and j (Fig. 2a). In this case, the unique solution to Eqs. (1)–(2) is ρ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 Fig. 2b, this symmetry implies that ρi = 1/N for all i and K = u as in the well-mixed case. 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 Fig. 3. From Eqs. (2) and (4) we can see that K > u is equivalent to the condition dr > dr, where the bars indicate averages over i = 1, . . ., N. This means that the molecular clock is accelerated if and only if di and ρi are positively correlated over sites; that is, if and only if there is a positive spatial correlation between the arrival of new mutations and the success they enjoy. These results led us to seek general conditions on the spatial structure leading to faster, slower, or the same molecular clock rates as a well-mixed population. We first find Result 1. If the death rates di are constant over all sites i = 1, . . ., N, then ρ = 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 di are constant over i, then since PN i¼1 di ¼ B, it follows that di = B/N for each i. Substituting in Eq. (3) yields ρ = 1/N. Another condition leading to K = u is the following: Result 2. If the birth rate equals the death rate at each site (bi = di for all i = 1, . . ., N), then ρ = 1/N, and consequently K = u. Moreover, bi = di for all i = 1, . . ., N if and only if the fixation probability is the same from each site (ρi = 1/N for all i = 1, . . ., N). PLOS Computational Biology | DOI:10.1371/journal.pcbi.1004108 February 26, 2015 5 / 32 Asymmetric Spatial Structure Alters the Molecular Clock Rate Fig 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 (bi = 1 for all i). Solving Eqs. (1)–(2), we obtain site-specific fixation probabilities of ρH = 1/(1+n2) and ρL = n/(1+n2) for the hub and each leaf, respectively. From Eq. (4), the molecular clock rate is 2n K ¼ 1þn2 u, which equals u for n = 1 and is less than u for n ! 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 8 4 2 1 with self-replacement only in site 1. Solving Eqs. (1)–(2) we find r1 ¼ 15, r2 ¼ 15, r3 ¼ 15 and r4 ¼ 15. (The powers of two arise because there is twice as much gene flow in one direction as the other.) From Eq. (4), the molecular clock rate is K ¼ 16 u > u, thus the molecular clock is accelerated in this case. 13 doi:10.1371/journal.pcbi.1004108.g003 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 eij = eji for all i and j then K = u. Result 2 can be obtained by substituting ρi = 1/N for all i into Eq. (1) and simplifying to obtain bi = di 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 bi are constant over all sites i = 1, . . ., N, then ρ 1/N, and consequently K u, with equality if and only if the death rates di 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 Fig. 3a. Intuitively, the sites at which mutations occur most frequently are those with high death rates di; 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 eij) the fixation probability can take any value 0 ρ < 1, and consequently, the molecular clock can take any rate 0 K < Nu. PLOS Computational Biology | DOI:10.1371/journal.pcbi.1004108 February 26, 2015 6 / 32 Asymmetric Spatial Structure Alters the Molecular Clock Rate Fig 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 di is uniform across sites, as in this example (di = 0.2 for all i). (b) Result 2 asserts that ρi = 1/N for all i—again implying K = u—if and only if each site has birth rate equal to death rate, bi = di for all i, as in this example. Nodes are colored according to their rates of turnover di. doi:10.1371/journal.pcbi.1004108.g004 Fig 5. Result 4 shows that K can achieve any value 0 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 (ρH = 1) and impossible for those arising in leaves (ρL = 0). The overall fixation probability is equal by Eq. (3) to the rate of turnover at the hub: ρ = dH = 1−(N−1)a. The molecular clock rate is therefore K = Nuρ = 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 ρ 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 ! 0, we have K ! Nu. doi:10.1371/journal.pcbi.1004108.g005 PLOS Computational Biology | DOI:10.1371/journal.pcbi.1004108 February 26, 2015 7 / 32 Asymmetric Spatial Structure Alters the Molecular Clock Rate 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 Fig. 5. Any non-negative value of ρ less than 1 can be obtained by an appropriate choice of parameters (details in Methods). 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” (Fig. 6). The sizes of these subpopulations are labeled N" and N#, respectively. Each subpopulation is well-mixed, with replacement probabilities e" for each pair of upstream sites and e# for each pair of downstream sites. Dispersal between the subpopulations is represented by the replacement probabilities e! from each upstream site to each downstream site, and e from each downstream site to each upstream site. We assume there is net gene flow downstream, so that e! > e . Solving Eqs. (1)–(2), we find that the fixation probabilities from each upstream site and each downstream site, respectively, are r" ¼ r# ¼ e! ; N" e! þ N# e e : N" e! þ N# e ð5Þ These ﬁxation probabilities were previously discovered for a different model of a subdivided population [40]. Substituting these ﬁxation probabilities into Eq. (4) yields the molecular clock rate: K¼ N N" d" e! þ N# d# e u: B N" e ! þ N# e ð6Þ Fig 6. A model of a population divided into upstream and downstream subpopulations. Each subpopulation is well-mixed. The replacement probability eij equals e" if sites i and j are both upstream, e# if i and j are both downstream, e! if i is upstream and j is downstream, and e if i is downstream and j is upstream. We suppose there is net gene flow downstream, so that e! > e . 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" > d#. doi:10.1371/journal.pcbi.1004108.g006 PLOS Computational Biology | DOI:10.1371/journal.pcbi.1004108 February 26, 2015 8 / 32 Asymmetric Spatial Structure Alters the Molecular Clock Rate Above, d" and d# are the turnover rates in the upstream and downstream populations, respectively, and B = N" d"+N# d# is the total birth rate per time-step. In Methods, we show that K > u if and only if d" > d#; that is, the molecular clock is accelerated if and only if there is more turnover in the upstream population than in the downstream population. In the case of unidirectional gene flow, e = 0, the molecular clock rate is simply K = (d"/B) Nu. The quantity d"/B represents the relative rate of turnover in the upstream population, and can take any value in the range 0 d"/B < 1/N"; thus K takes values in the range 0 K < (N/ N")u. We note that the upper bound on K is inversely proportional to the size N" of the upstream population. The largest possible values of K are achieved when N" = 1, in which case K can come arbitrarily close to Nu. These bounds also hold if there are multiple downstream subpopulations, since for unidirectional gene flow, the spatial arrangement of downstream sites does not affect the molecular clock rate. In particular, if the hub and leaves in Fig. 5 are each replaced by well-mixed subpopulations, then K is bounded above by (N/NH)u, where NH is the size of the hub subpopulation. 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 [41]. 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 (Fig. 7; [42–44]). We focus on the proximal small intestine in mice, for which recent studies [41, 45] suggest there are * 5 active stem cells per crypt, each replaced * 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/ day between neighbors, so that di = 0.1/day for all i. 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, ρi = 1/5 for each of the five stem cell sites. The molecular clock rate is therefore P5 K ¼ u i¼1 di ri ¼ 0:1u substitutions per day. This accords with the empirical finding that, for a neutral genetic marker with mutation rate u % 1.1×10−4, substitutions accumulate at a rate ~ K % 1:1 Â 10À5 per crypt per day [41]. Does crypt architecture limit the rate of genetic change in intestinal tissue? Intestinal crypts in mice contain * 250 cells and replace all their cells about once per day [46]. If each crypt ~ were a well-mixed population, the molecular clock rate would be K ¼ Bu=N % 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 ideas Our results can also be applied to ideas that spread by imitation on social networks. In this setting, a mutation corresponds to a new idea that could potentially replace an established one. Neutrality means that all ideas are equally likely to be imitated. To investigate whether human social networks accelerate or slow the rate of idea substitution, we analyzed 973 Twitter networks from the Stanford Large Network Dataset Collection [47]. 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 PLOS Computational Biology | DOI:10.1371/journal.pcbi.1004108 February 26, 2015 9 / 32 Asymmetric Spatial Structure Alters the Molecular Clock Rate Fig 7. A simple model of cell replacement structure in epithelial crypts of the small intestine, based on results of [41] and [45]. A small number of stem cells (N" * 5) residing at the bottom of the intestinal crypt and are replaced at rate d" * 0.1 per stem cell per day. Empirical results [41, 45] suggest a cycle structure for stem cells. To achieve the correct replacement rate we set eij = 0.05/day for each neighboring pair. Stem cells in an individual crypt replace a much larger number of progenitor and differentiated cells (* 250; [46]). These downstream progenitor and differentiated cells are replaced about every day [46]. The hierarchical organization of intestinal crypts, combined with the low turnover rate of stem cells, limits the rate of neutral ~ genetic substitutions (K~ % 0:1u substitutions per day), since only mutations that arise in stem cells can fix. doi:10.1371/journal.pcbi.1004108.g007 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. To model the spread of ideas on these networks, we set eij = 1/L if j follows i and zero otherwise, where L is the total number of links. This can be instantiated by supposing that at each time-step, one followee-follower link is chosen with uniform probability. The follower either adopts the idea of the followee, with probability 1−u, or innovates upon it to create a new idea, with probability u, where u ( 1. With these assumptions, the resulting rate of idea substitution (as a multiple of u) depends only on the network topology and not on any other parameters. 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 di, thus K = u follows from Result 1 for these networks. We found a weak but statistically significant negative relationship between the PLOS Computational Biology | DOI:10.1371/journal.pcbi.1004108 February 26, 2015 10 / 32 Asymmetric Spatial Structure Alters the Molecular Clock Rate Fig 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 (slope % −0.00164 with 95% confidence interval (−0.0023, −0.001) based on the bootstrap method; R % −0.45). 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. doi:10.1371/journal.pcbi.1004108.g008 PLOS Computational Biology | DOI:10.1371/journal.pcbi.1004108 February 26, 2015 11 / 32 Asymmetric Spatial Structure Alters the Molecular Clock Rate Fig 9. Illustration of idea generation and spreading on Twitter “ego networks”. Network structure accelerates idea substitution (K > u) if and only if there is a positive spatial correlation between the generation of new ideas (which for our model occurs proportionally to the rate di of incoming ideas) and the probability of fixation ρi. Panels (a) and (b) show the networks with the slowest (K % 0.667u) and fastest (K % 1.085u) rates of idea substitution, respectively, among networks of size 13. The coloration of nodes corresponds to their rate of turnover di, with warmer colors indicating more rapid turnover. The size of nodes corresponds to their fixation probability ρi. doi:10.1371/journal.pcbi.1004108.g009 network size N and value K/u (slope % −0.00164 with 95% confidence interval (−0.0023, −0.001) based on the bootstrap method; R % −0.45). This negative relationship persists even if small networks with less than 10 nodes are removed (slope % −0.00156 with 95% confidence interval (−0.0023, −0.0009); R % −0.43). In summary, while some Twitter ego-networks accelerate the substitution of neutral innovations, the vast majority slow this rate (Figs. 8 and 9). One possible explanation for the rarity of networks that accelerate idea substitution has to do with the intrinsic relationship between the turnover rates di and the site-specific fixation probabilities ρi. From Eq. (1), we see that ρi can be written as the product P N ð1=di Þ Â eij rj , where the first factor can be interpreted as the “attention span” of node i j¼1 and the second can be interpreted as its influence. While these two factors are not strictly independent, we would not necessarily expect a systematic relationship between them in our Twitter network ensemble. In the absence of such a relationship, ρi is inversely related to di, which implies K < u. In other words, the most fertile nodes (in terms of generating new ideas) are also the most fickle (in terms of adopting the ideas of others); thus many new ideas are abandoned as soon as they are generated. This heuristic argument suggests that K > u, while possible, might be an uncommon occurrence in networks drawn from statistical or probabilistic ensembles. Discussion The spatial structure of a population affects its evolution in many ways, for example by promoting cooperative behaviors [36, 48–54], genetic variation [55–58], and speciation [59–61]. Asymmetric spatial structure in particular is known to have important consequences for adaptation [23, 27, 28, 35, 62–65] and for genetic diversity [66–68]. 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 probabilities of turnover (di) and net offspring production (bi−di) across sites. If both di and bi−di differ across sites, the molecular clock rate will in general differ PLOS Computational Biology | DOI:10.1371/journal.pcbi.1004108 February 26, 2015 12 / 32 Asymmetric Spatial Structure Alters the Molecular Clock Rate from that of a well-mixed population (Result 4). If additionally bi is constant across sites, Result 3 guarantees that neutral substitution is slowed relative to the well-mixed case. Our Result 4 shows that the rate of neutral substitution in a population can come arbitrarily close to Nu. However, this result depends on the existence of a single “hub” individual seeding the rest of the population as in Fig. 5. A similar but more plausible scenario (especially in sexually reproducing populations) involves a well-mixed hub subpopulation seeding one or more “leaf” subpopulations. In this case, the upper bound on the molecular clock rate is (N/NH)u, where NH is the size of the hub subpopulation. Conditions leading to altered molecular clock rate may occur frequently in natural populations. Asymmetric dispersal may result, for example, from prevailing winds [69, 70], water currents [71, 72], or differences in elevation [73, 74]. Differences in habitat quality may lead to variance in birth and death rates across sites (nonuniform bi and di). It is known that such asymmetries in spatial structure have important consequences for adaptation [23, 62–64] and for genetic diversity [66–68]; 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. One important assumption made in our work is that mutation occurs with a constant probability per reproduction. Alternatively, one might suppose that heritable mutations accrue at a constant rate per individual (e.g. due to germline cell divisions). With this alternate assumption, mutations would arise at uniform rates over sites, resulting in a molecular clock rate of K = u for all spatial structures. The applicability of our results thus depends on the mutation process of the population in question, which may in many cases lie between these extremes. In humans, for example, recent evidence suggests that maternal mutations occur with constant probability per reproduction, whereas paternal mutations (which are more frequent) increase in probability with the father’s age [75–77]. In general, we expect deviations from K = u to scale with the extent to which mutations in a lineage depend on the number of generations rather than chronological time. Many somatic cell populations have strongly asymmetric patterns of replacement, with a small number of stem cell pools feeding a much larger number of progenitor and differentiated cells. The rate at which mutations accumulate in these populations has significant implications for the onset of cancer [78–80] and the likelihood of successful cancer therapy [81–84]. It is well-known that this rate is proportional to the rate of stem cell division, since only mutations that arise in stem cells can persist [45, 85–88]. Our work shows this how principle arises from a general analysis of neutral evolution in spatially structured populations. 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 [41] and/or allowing differentiated cells to revert to stem cells [89]. This restructuring may, in turn, alter the rate of genetic substitution, with further ramifications for cancer progression. The influence of social network topology on the spread of ideas and behaviors is a question of both theoretical and practical interest [90–97]. The neutral substitution rate K on social networks describes how innovations spread when they are equally likely to be imitated as an existing convention. Our finding that most Twitter ego networks have rates less than those of wellmixed populations contrasts with results from epidemiological models, which generally find that the heterogeneity of real-world social networks accelerates contagion [90, 98, 99]. We note, however, that this finding is sensitive to the assumption that individuals generate new ideas in proportion to the rate of incoming ideas (which itself is proportional to the number of individuals followed). Since the success of an idea varies according to the node at which it arises, the overall rate of substitution depends on the distribution of new ideas among nodes. PLOS Computational Biology | DOI:10.1371/journal.pcbi.1004108 February 26, 2015 13 / 32 Asymmetric Spatial Structure Alters the Molecular Clock Rate For example, if one were to instead assume that each individual generates new ideas at an equal rate, one would find that the network topology has no effect on the rate of substitution. 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 structure changes over time—due, for instance, to changes in climate, tumorogenesis, or social network dynamics [100–103]—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ρ) 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 [22, 30, 58, 104], leading to further ramifications for the molecular clock rate [105]. 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 frequencydependent selection [65, 106]. Re-analyzing existing models using this new mutation convention may reveal further surprises about how spatial structure affects evolution. Fig 10. 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 α is given by α(2) = 3, α(4) = 4, α(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 s = (M, R, R, M, R) to s0 = (M, R, R, M, M). doi:10.1371/journal.pcbi.1004108.g010 PLOS Computational Biology | DOI:10.1371/journal.pcbi.1004108 February 26, 2015 14 / 32 Asymmetric Spatial Structure Alters the Molecular Clock Rate Methods Class of Models In the class of models we consider, there are N sites indexed i = 1, . . ., 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, α), where R & {1, . . ., N} is the set of sites whose occupants are replaced by new offspring, and α:R ! {1, . . ., 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 Fig. 10. 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, α) in this distribution will be denoted p(R, α). Neutrality is represented by independence of the probabilities p(R, α) 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 2 {1, . . ., N}, a positive integer n and a finite sequence fðRk ; ak Þgn of replacement events such that k¼1 • p(Rk, αk) > 0 for all k, and • For all individuals j 2 {1, . . ., N}, ak1 ak2 Á Á Á akm ðjÞ ¼ i; where k1, . . . < km is the maximal subsequence of 1, . . ., n such that the compositions in Eq. (7) are well-deﬁned. We observe that Eq. (7) 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. It is also equivalent to the statement that the weighted digraph with edge weights eij is out-connected from at least one vertex. 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, α), the sites that are replaced by the offspring of i is the given by the preimage α−1(i) & R, i.e., the set of indices that map to i under α. The number of offspring of site i is equal to jα−1(i)j, 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 X pðR; aÞ jaÀ1 ðiÞj: bi ¼ E½jaÀ1 ðiÞj ¼ ðR;aÞ ð7Þ The death rate (probability of replacement) of site i is equal to P di ¼ Pr ½i 2 R ¼ pðR; aÞ: ðR; aÞ i2R The probability that the offspring of i displaces the occupant of j in a replacement event is P pðR; aÞ: eij ¼ Pr ½ j 2 R and aðjÞ ¼ i ¼ ðR; aÞ j2R aðjÞ ¼ i We observe that bi ¼ PN j¼1 eij and di ¼ PN j¼1 eji . PLOS Computational Biology | DOI:10.1371/journal.pcbi.1004108 February 26, 2015 15 / 32 Asymmetric Spatial Structure Alters the Molecular Clock Rate 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 si 2 {M, R}. The overall state of the process can be recorded as a string s = (s1, . . ., sN) of length N with alphabet {M, 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 s = (s1, . . ., sN) and the replacement event (R, α) occurs, then the new state s0 ¼ ðs01 ; . . . ; s0N Þ is given by & i2 R = si s0i ¼ saðiÞ i 2 R: The above assumptions describe a Markov chain on the set of strings of length N with alphabet {M, 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: (M, . . ., M) or (R, . . ., 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 A. Instead of two types, the ancestral Markov chain has N types, labeled 1, . . ., N, which correspond to the N members of a “founding generation” of the population. Evolution proceeds according to the given replacement rule, as described above. The states of the ancestral Markov chain are strings of length N with alphabet {1, . . ., N}. The ancestral Markov chain has a canonical initial state a0 = (1, . . ., 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 a0— in our notation, (A, a0)—has the useful feature that at any time t ! 0, the state a(t) = (a1(t), . . ., aN(t)) indicates the site occupied by each individual’s founding ancestor. In other words, if aj(t) = i, then the current occupant of site j is descended from the founder that occupied site i. To relate the evolutionary Markov chains A and M, consider any set mapping γ:{1, . . ., N} ! {M, R}. We think of γ(i) as giving the genetic type (M or R) of each member i = 1, . . ., N of the founding generation. The mapping γ induces a mapping ~ from states of A to states of M, g defined by ~ða1 ; . . . ; aN Þ ¼ ðgða1 Þ; . . . ; gðaN ÞÞ. For any state a(t) of (A, a0), the string g ~ðaðtÞÞ 2 f1; . . . ; NgN indicates the genetic type of each individual’s ancestor in the founding g generation. Since genetic types are inherited faithfully, it follows that ~ðaðtÞÞ gives the current g genetic type of each individual. Thus if M and A follow the same replacement rule, we have that for any such mapping γ and any string s 2 {M, R}N, g PðA;a0 Þ ½~ðaðtÞÞ ¼ s ¼ PðM;~ða0 ÞÞ ½sðtÞ ¼ s: g ð8Þ Basic results on fixation probabilities We define the fixation probability from site i, ρ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: ri ¼ lim PðM;mi Þ ½sðtÞ ¼ ðM; . . . ; MÞ: t!1 ð9Þ Above, mi denotes the initial state consisting of an M in position i and R’s elsewhere. The PLOS Computational Biology | DOI:10.1371/journal.pcbi.1004108 February 26, 2015 16 / 32 Asymmetric Spatial Structure Alters the Molecular Clock Rate ordered pair (M, mi) refers to the Markov chain M with initial state mi. P(M, mi)[s(t) = s] denotes the probability that the state of (M, mi) is s at time t ! 0. We can use the relationship between the ancestral Markov chain A and the evolutionary Markov chain M to obtain an alternate expression for the site-specific fixation probability ρi: Theorem 1. ρi = limt ! 1 P(A, a0) [a(t) = (i, . . ., i)]. In words, the site-specific fixation probability ρi equals the probability that founding individual i becomes the eventual ancestor of the whole population. Proof. For any i = 1, . . ., N, define the set mapping γi:{1, . . ., N} ! {M, R} by ( gi ðjÞ ¼ M R if j ¼ i otherwise: (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 ~i ða0 Þ ¼ mi . Comg bining Eq. (8) with the deﬁnition of ρi, we obtain ri ¼ lim PðM;mi Þ ½sðtÞ ¼ ðM; . . . ; MÞ t!1 ¼ lim PðA;a0 Þ ½~i ðaðtÞÞ ¼ ðM; . . . ; MÞ g t!1 ¼ lim PðA;a0 Þ ½aðtÞ ¼ ði; . . . ; iÞ: t!1 More generally, we can consider the probability of fixation from an arbitrary set of sites. For any set S & {1, . . ., N}, we let ρS denote the probability that the mutant type becomes fixed, given the initial state mS with mutants occupying the sites specified by S and residents occupying all other sites: rS ¼ lim PðM;mS Þ ½sðtÞ ¼ ðM; . . . ; MÞ: t!1 Site-specific fixation probabilities are additive in the following sense: Theorem 2. For any set S & {1, . . .N} of sites, rS ¼ In particular, N X ri ¼ 1: i¼1 X i2S ri : ð10Þ ð11Þ This result has previously been obtained for a number of specific evolutionary processes on graphs [39, 107, 108]. Intuitively, the probability of fixation from the initial state described by S equals the probability that one of the individuals in a site i 2 S becomes the eventual ancestor of the population. Since this cannot be true of more than one site in S, the overall probability ρS is obtained by summing over all i 2 S the individual probabilities ρi that site i contains the eventual ancestor. PLOS Computational Biology | DOI:10.1371/journal.pcbi.1004108 February 26, 2015 17 / 32 Asymmetric Spatial Structure Alters the Molecular Clock Rate Proof. Suppose that we are given a set S & {1, . . ., N} of sites initially occupied by mutants. This situation is described by the set mapping γS:{1, . . ., N} ! {M, R}, given by & M if i 2 S gS ðiÞ ¼ R otherwise: Invoking the relationship between A and M—in particular, Eq. (8) and Theorem 1—we have rS ¼ lim PðM;g~S ða0 ÞÞ ½sðtÞ ¼ ðM; . . . ; MÞ t!1 ~ ¼ lim PðA;a0 Þ ½gS ðaðtÞÞ ¼ ðM; . . . ; MÞ t!1 ¼ lim PðA;a0 Þ ½aðtÞ ¼ ði; . . . ; iÞ for some i 2 S t!1 X ¼ lim PðA;a0 Þ ½aðtÞ ¼ ði; . . . ; iÞ ¼ X i2S i2S t!1 ri : This proves Eq. (10). Eq. (2) now follows from letting S = {1, . . ., N}, and noting that, in this case, rS ¼ lim PðM;ðM;...;MÞÞ ½sðtÞ ¼ ðM; . . . ; MÞ ¼ 1: t!1 We now derive Eq. (1), which allows the fixation probabilities ρi to be calculated from the replacement rates eij. Theorem 3. For each i = 1, . . ., N, di ri ¼ N X j¼1 eij rj : Proof. Considering the change that can occur over a single time-step, we have the following recurrence relation: X ri ¼ ð1 À di Þri þ pðR; aÞ raÀ1 ðiÞ : ðR;aÞ The ﬁrst 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 −di)ρi from both sides and applying Theorem 2 with S = ρα−1(i) yields ! X X pðR; aÞ rj : di ri ¼ ðR;aÞ j2aÀ1 ðiÞ Now interchanging the summations on the right-hand side yields 0 1 C N X BX B C B rj B pðR; aÞC: di ri ¼ C @ ðR; aÞ A j¼1 j2R aðjÞ ¼ i PLOS Computational Biology | DOI:10.1371/journal.pcbi.1004108 February 26, 2015 18 / 32 Asymmetric Spatial Structure Alters the Molecular Clock Rate By deﬁnition, X ðR; aÞ j2R aðjÞ ¼ i pðR; aÞ ¼ eij : This completes the proof. Proof of Result 2 Theorem 4 (Result 2). The fixation probabilities from each site are equal to 1/N (ρi = 1/N for all i = 1, . . ., N) if and only if each site has birth rate equal to death rate (bi = di for all i = 1, . . ., N). Proof. Assume first that the fixation probabilities from each site are all equal to 1/N. Substituting ρi = 1/N for all i into Eq. (1) yields di = bi. This proves the “only if” direction. Next assume that the birth rate is equal to the death rate at each site (bi = di for all i = 1, . . ., N). Eq. (1) can then be rewritten as N X j¼1 eij ri ¼ N X eij rj : j¼1 PN Clearly, ρi = 1/N for all i satisﬁes the above equation for all i, and also satisﬁes i¼1 ri ¼ 1. Assumption 1 guarantees that the solution to these equations is unique. Therefore bi = di for all i implies ρ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, bi = 1/N for all i = 1, . . ., N, then ρ 1/N (and consequently K u) with equality if and only if the death rates are also constant over sites. Proof. We separate our proof into six 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 PN N N points of the function r ¼ i¼1 di ri with respect to the variables feij gi;j¼1 and fri gi¼1 and the constraints 8 N 1 >X > > eij ¼ > > N > j¼1 > > > N ! < X N X eji ri ¼ eij rj > j¼1 > j¼1 > > N >X > > > > ri ¼ 1: : i¼1 i 2 f1; . . . ; Ng; i 2 f1; . . . ; Ng; ð12Þ We obtain that the critical points are precisely those for which di = 1/N for all i. In the next three steps, we use partial derivatives of Eqs. (1), (2) and (3) to form a second-order Taylor expansion of ρ in the variables {eij}i 6¼ j around the critical points. Finally, we show that the critical points are maxima, 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 ∑i, j eij, is one. We will demonstrate this by showing that, if a PLOS Computational Biology | DOI:10.1371/journal.pcbi.1004108 February 26, 2015 19 / 32 Asymmetric Spatial Structure Alters the Molecular Clock Rate maximum of ρ is achieved by a combination of replacement rates feÃ gi;j , this maximum is also ij achieved using the normalized rates feij0 gi;j with e0ij ¼ eÃ =BÃ , BÃ ∑i, j eÃ . ij ij First, we see that the equations diÃ ri ¼ N X eÃ rj ij j¼1 N X j¼1 and di0 ri ¼ e0ij rj are equivalent, differing only by a factor of B. (Above, we have introduced the notation diÃ ¼ PN Ã PN 0 N 0 j¼1 eij and di ¼ j¼1 eij .) Thus the node-speciﬁc ﬁxation probabilities fri gi¼1 resulting from feij0 gi;j are the same as those resulting from feÃ gi;j . ij Turning now to the overall fixation probabilities ρÃ and ρ0 corresponding to feÃ gi;j and ij feij0 gi;j respectively, we see that rÃ ¼ N N N X dÃ X 1X Ã i di ri ¼ ri ¼ di0 ri ¼ r0 : BÃ BÃ i¼1 i¼1 i¼1 Thus the same maximum is achieved by the normalized replacement rates feij0 gi;j , which satisfy P 0 i;j eij ¼ 1. We conclude that there is no loss of generality in assuming B = 1. Together with uniform birth rates, this implies the constraint N X j¼1 eij ¼ 1 : N ð13Þ Step 2: Determine critical points PN We use the method of Lagrange multipliers to find critical points of the function r ¼ i¼1 di ri N N with respect to the variables feij gi;j¼1 and fri gi¼1 . Eqs. (13), (1), and (2) give the respective constraint equations gi N X 1 eij À ¼ 0; for i 2 f1; . . . ; Ng N j¼1 ! N N X X ri eji À eij rj ¼ 0; for i 2 f1; . . . ; Ng hi q j¼1 N i¼1 X ri À 1 ¼ 0: i¼1 Critical points with respect to the above constraints are solutions to rr ¼ N N X X li rgi þ mi rhi þ srq; i¼1 i¼1 where the λi’s, μi’s, and σ are the Lagrange multipliers. The individual components of the gradient are obtained by taking the partial derivative with respect to ekl, yielding rl ¼ lk þ rl ml À rl mk ; ð14Þ PLOS Computational Biology | DOI:10.1371/journal.pcbi.1004108 February 26, 2015 20 / 32 Asymmetric Spatial Structure Alters the Molecular Clock Rate for k, l 2 {1, . . ., N}, and the partial derivative with respect to ρk, yielding N X j¼1 ejk ¼ mk N X j¼1 ejk À N X mj ejk þ s; j¼1 ð15Þ for k = 1, . . ., N. Solving Eq. (14) for ρl gives rl ¼ lk : 1 À ml þ mk ð16Þ Note that in the case l = k, Eq. (16) becomes ρk = λk. Therefore, we replace λk with ρk in Eq. (16): rl ¼ Interchanging k and l, we obtain rk ¼ rl : 1 À mk þ ml ð18Þ rk : 1 À ml þ mk ð17Þ Combining Eqs. (17) and (18) yields μl = μk and ρl = ρk for all k, l 2 {1, . . ., N}. That is, all node-specific probabilities ρi are equal. It follows from Eq. (2) that ρi = 1/N for all i = 1, . . ., N. Furthermore, given that μl = μk, Eq. (15) yields dk = σ for all k = 1, . . ., N. That is, all death rates di are equal. Therefore di = 1/N for i = 1, . . ., N because we have assumed (without loss of genPN erality) that i¼1 di ¼ 1. In summary, we have shown that the overall fixation probability ρ has a critical point whenever all node-specific fixation probabilities and death rates are constant (di = ρi = 1/N for all i = 1, . . ., N). Consequently, ρ = 1/N at all critical points. It still remains to prove that this critical value of ρ 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 feÃ g i; j . Then, viewing ρ as a function of the independent variables feij g i; j , we form a ij i 6¼ j i 6¼ j second-order Taylor expansion of ρ around the critical point. (Since we are operating under PN the constraint j¼1 eij ¼ 1=N for all i, this set of variables suffices to determine the value of ρ.) 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 ρ about the critical point feij g i; j can be written as i 6¼ j X @r À 1 X @2r 1 3Á Dekl þ Dekl Demn þ O jD~ ; ej þ r¼ @ekl 2 k; l; m; n @ekl @emn N k; l k 6¼ l m 6¼ n k 6¼ l ð19Þ where all derivatives are taken at feÃ g i; j and Dekl ¼ ekl À eÃ . More simply, we write this expanij kl i 6¼ j sion as r¼ À 1 3Á þ rð1Þ þ rð2Þ þ O jD~ ; ej N PLOS Computational Biology | DOI:10.1371/journal.pcbi.1004108 February 26, 2015 21 / 32 Asymmetric Spatial Structure Alters the Molecular Clock Rate with ρ(1) and ρ(2) representing the ﬁrst- and second-order terms, respectively, in Eq. (19). The ﬁrst-order term ρ(1) is zero since feÃ g i; j is a critical point. Our goal will be to show that the secij i 6¼ j ond-order term, ρ(2) is always negative or zero. We can find an alternative expression for ρ(2) using the definition of ρ in Eq. (3) as follows. We introduce the notation Ddi ¼ di À diÃ ¼ di À 1=N. With this notation, Eq. (3) becomes r¼ N N X 1X ri þ Ddi ri N i¼1 i¼1 N 1 X ¼ þ Ddi ri : N i¼1 ð20Þ We substitute the ﬁrst-order Taylor series expansion for ρi, À 1 2Á þ rð1Þ þ O jD~ ; ej i N PN into Eq. (20), noting that Δdi = O(jΔe and i¼1 Ddi ¼ 0. This yields ~j) N À 1 X 1 2Á þ rð1Þ þ O jD~ Ddi ej r ¼ þ i N N i¼1 ri ¼ ¼ N À 1 X 3Á þ Ddi rð1Þ þ O jD~ : ej i N i¼1 Thus, the second-order term of the overall ﬁxation probability can be written as rð2Þ ¼ N X Ddi rð1Þ : i i¼1 ð21Þ Step 4: Obtain recursion for first-order term of the site-specific fixation probability We now investigate the properties of rð1Þ . We begin by rewriting Eq. (1) as i X eij rj : di ri ¼ eii ri þ j6¼i 1 Replacing eii with N À P j6¼i ij e and simplifying, we obtain X eij rj À ri : Ddi ri ¼ j6¼i Next we take the partial derivative of both sides with respect to ekl, where k 6¼ l, X @r X @eij @ðDdi Þ @r @r j ri þ Ddi i ¼ rj À ri þ eij À i ; @ekl @ekl @ekl @ekl @ekl i; j i; j j 6¼ i j 6¼ i 1 and evaluate at the critical point (rÃ ¼ N and Δdi = 0): i 1 @ðDdi Þ X Ã @rj @ri ¼ eij À : N @ekl @ekl @ekl j6¼i PLOS Computational Biology | DOI:10.1371/journal.pcbi.1004108 February 26, 2015 22 / 32 Asymmetric Spatial Structure Alters the Molecular Clock Rate Multiplying both sides by Δekl and then summing over all k, l 2 {1, . . ., N} for k 6¼ l yields X @rj X @r 1 X @ðDdi Þ Dekl ¼ eÃ Dekl À eÃ i Dekl : ij ij @ekl @ekl N k; l @ekl j; k; l j; k; l k 6¼ l k 6¼ l j 6¼ i k 6¼ l j 6¼ i ð22Þ We observe that rð1Þ ¼ i P @ri k; l k 6¼ l @ekl of ρi about the critical point, and ekl. Thus Eq. (22) becomes P Dekl is the first order term of the Taylor series expansion k; l k 6¼ l @ðDdi Þ @ekl Dekl ¼ Ddi , since Δdi is a linear function of the X X ð1Þ 1 Ddi ¼ eÃ rð1Þ À eÃ ri : ij j ij N j6¼i j6¼i The restriction j 6¼ i in the sums on the right-hand side can be removed, since the entire righthand side is zero when j = i. Rearranging further, we have Nrð1Þ i N N X ð1Þ X eÃ ¼ ÀDdi þ N eÃ rj : ij ij j¼1 j¼1 Since the birth rate is uniform over sites ðbi ¼ ﬁes to PN j¼1 1 eÃ ¼ N for all i), the above equation simpliij rð1Þ ¼ ÀDdi þ N i N X ð1Þ eÃ r j : ij j¼1 ð23Þ This equation provides a recurrence relation for the ﬁrst-order term of the Taylor expansion of ρi about the critical point. Step 5: Show second-order term of the overall fixation probability is nonpositive We now combine Eq. (23) with Eq. (21) to show that ρ(2) 0, with equality if and only if Δdi = 0 for all i, according to the following argument. We first multiply both sides of Eq. (23) by rð1Þ and sum over i = 1, . . . , N: i N X i¼1 rð1Þ i 2 ¼À PN N X X Ddi rð1Þ þ N eÃ rð1Þ rð1Þ : i j ij i i¼1 i;j By Eq. (21) we can substitute ρ(2) for rð2Þ ¼ À i¼1 Ddi rð1Þ . Then, solving for ρ(2), we obtain i 2 þN X i;j N X i¼1 rð1Þ i eÃ rð1Þ rð1Þ : j ij i ð24Þ We now make use of the fact that the product rð1Þ rð1Þ can be written as a difference of i j squares: 2 2 2 ! 1 ð1Þ : ri À rð1Þ À rð1Þ À rð1Þ rð1Þ rð1Þ ¼ À i j j i j 2 PLOS Computational Biology | DOI:10.1371/journal.pcbi.1004108 February 26, 2015 23 / 32 Asymmetric Spatial Structure Alters the Molecular Clock Rate Substituting this identity into Eq. (24) yields rð2Þ ¼ À N 2 N X 2 X ð1Þ 2 N X ri À eÃ rð1Þ À rð1Þ þ eÃ rð1Þ i j ij 2 i;j 2 i;j ij i i¼1 N X Ã ð1Þ 2 þ e r 2 i;j ij j N N N 2 N X 2 X X ð1Þ 2 N X ð1Þ ¼À ri À eÃ ri À rð1Þ þ rð1Þ eÃ j i ij ij 2 i;j 2 i¼1 i¼1 j¼1 N N N X ð1Þ 2 X Ã þ rj eij : 2 j¼1 i¼1 eÃ ¼ 1=N for each i. Furthermore, since feÃ gi;j is a critiij ij PN Ã cal point, death rates are uniform as well; thus i¼1 eij ¼ 1=N for each j. Making these substitutions yields Uniform birth rates implies that j¼1 PN rð2Þ ¼ À N N 2 1 X 2 1 X 2 N X Ã ð1Þ eij ri À rð1Þ þ rð1Þ þ rð1Þ j 2 i;j 2 i¼1 i 2 j¼1 j i¼1 2 N X Ã ð1Þ ¼À eij ri À rð1Þ : j 2 i;j N X rð1Þ i 2 À ð25Þ This shows that ρ(2) 0. We now consider the case that ρ(2) = 0. By Eq. (25), this occurs if and only if rð1Þ ¼ rð1Þ for i j any pair i, j with eÃ 6¼ 0. It follows that each term eÃ rð1Þ in the sum on the right-hand side of ij ij j PN Ã ð1Þ Eq. (23) can be replaced by eij ri . Making these substitutions and using j¼1 eÃ ¼ 1=N, we ij find that Δdi = 0. Thus ρ(2) = 0 implies that Δdi = 0 for each i. Conversely, if Δdi = 0 for each i, then Eq. (21) implies that ρ(2) = 0. Overall, we have shown that ρ(2) 0 with equality if and only if Δdi = 0 for each i. Step 6: Show that the critical points are maxima We now adopt a vector notation, in which e ¼ feij g i; j is an arbitrary combination of replacei 6¼ j ment probabilities satisfying the conditions of the theorem, and eÃ ¼ feÃ g i; j is an arbitrary ij i 6¼ j critical point of ρ. By Taylor’s theorem, we can rewrite the second-order expansion (19) of ρ around eÃ as rðeÞ ¼ 1 1 T T þ ðe À eÃ Þ HeÃ ðe À eÃ Þ þ ðe À eÃ Þ Re;eÃ ðe À eÃ Þ; N 2 ð26Þ where HeÃ is the Hessian matrix of ρ at eÃ and {Re;eÃ } is a family of symmetric matrices with the property that Re;eÃ converges to the zero matrix as e ! eÃ . Since HeÃ and Re;eÃ are symmetric, all eigenvalues of these matrices are real, and the eigenspaces corresponding to these eigenvalues are orthogonal. We showed in Step 5 that ρ(2) 0, which implies that HeÃ is negative semidefinite. We also showed that ρ(2) = 0 if and only if Δdi = 0 for all i; thus the null space of HeÃ is the vector space D consisting of all displacement vectors for which Δdi = 0 for all i. It follows that HeÃ is negative definite as an operator on the orthogonal complement D? of D. PLOS Computational Biology | DOI:10.1371/journal.pcbi.1004108 February 26, 2015 24 / 32 Asymmetric Spatial Structure Alters the Molecular Clock Rate We now fix a particular critical point eÃ , and consider a neighborhood NðeÃ ; eÞ of the form 0 0 È Ã É NðeÃ ; eÞ ¼ e0 þ v þ w j v 2 D; jjvjj < e; w 2 D? ; jjwjj < e : 0 For each point eÃ þ v þ w 2 NðeÃ ; eÞ, we observe that eÃ þ v is also a critical point of ρ, since 0 0 0 the addition of v 2 D leaves the value of di unchanged for each i. Now, applying the Taylor expansion (26) of rðeÃ þ v þ wÞ about eÃ þ v, yields 0 0 rðeÃ þ v þ wÞ ¼ 0 1 1 T þ w ðHeÃ þv Þw þ wT ðReÃ þvþw;eÃ þv Þw: 0 0 0 N 2 0 ð27Þ Let λH(v) denote the largest nonzero eigenvalue of HeÃ þv , and let λR(v, w) denote the largest eigenvalue of ReÃ þvþw;eÃ þv . Since HeÃ þv is negative semidefinite, λH(v) < 0. Further, since HeÃ þv is negative definite as an operator on D? and since w 2 D?, it follows that wT HeÃ þv w 0 0 0 0 0 lH ðvÞjjwjj : 2 ð28Þ We also have that, for each eÃ þ v þ w 2 NðeÃ ; eÞ, 0 0 wT ReÃ þvþw;eÃ þv w 0 0 0 0 lR ðv; wÞjjwjj : 2 ð29Þ Furthermore, since ReÃ þvþw;eÃ þv converges to the zero matrix as w ! 0, we have limw ! 0 λR(v, w) = 0 for each ﬁxed v. It follows that, for sufﬁciently small ε > 0, all points eÃ þ v þ w 2 0 NðeÃ ; eÞ satisfy 0 1 l ðvÞ þ lR ðv; wÞ < 0: 2 H ð30Þ Combining Eqs. (27)–(30), we have that for all eÃ þ v þ w 2 NðeÃ ; eÞ with ε sufﬁciently 0 0 small, 1 2 Ã rðe0 þ v þ wÞ 1=N þ l ðvÞ þ lR ðv; wÞ jjwjj 1=N: ð31Þ 2 H Equality is achieved in Eq. (31) if and only if w = 0, in which case eÃ þ v þ w is another critical 0 point. Thus ρ has a local maximum of 1/N at every critical point. It follows that these points are global maxima as well. We showed in Step 2 that the critical points of ρ are precisely those for which di = 1/N for all i. This completes the proof. 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 Fig. 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 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 dH = 1−(N−1)a for the hub and dL = a for each leaf. Solving Eqs. (1)–(2) we obtain the site-specific fixation probabilities ρL = 0 for each leaf and ρH = 1 for the hub. By Eq. (3), the overall fixation probability is equal to the rate of turnover at the hub: r ¼ dH ¼ 1 À ðN À 1Þa: Any value 0 ρ < 1 can be obtained by an appropriate choice of a, with 0 < a specifically, a = (1−ρ)/(N−1). This proves: 1/(N−1), PLOS Computational Biology | DOI:10.1371/journal.pcbi.1004108 February 26, 2015 25 / 32 Asymmetric Spatial Structure Alters the Molecular Clock Rate Theorem 6 (Result 4). For arbitrary spatial population structure (no constraints on eij) the fixation probability can take any value 0 ρ < 1, and consequently, the molecular clock can take any rate 0 K < Nu. We observe that if a = 1/N then the death rate is 1/N at each site, and therefore ρ = 1/N by Result 1. If a = 1/(N−1) then there is no turnover at the hub and thus ρ = 0. At the other extreme, as a approach zero, ρ comes arbitrarily close to unity. Exact expressions for ρ with N 3 To complement the above arguments, which apply to general population size N, we here derive exact expressions for the overall fixation probability ρ 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 [40]. Population size N = 2 We solve Eqs. (1)–(2) to obtain expressions for ρ1 and ρ2: r1 ¼ e12 ; e12 þ e21 e21 : e12 þ e21 ð32Þ r2 ¼ ð33Þ P2 We first consider uniform birth rates (bi ¼ j¼1 eij ¼ 1 for i = 1, 2). Substituting Eqs. (32) 2 and (33) into Eq. (3) (B = 1 in this case), gives an expression for the overall fixation probability: 1 ðe À e21 Þ r ¼ À 12 : e12 þ e21 2 2 Thus, in accordance with Result 3, ρ is less than or equal to 1 and consequently K u. Equality 2 occurs if and only if e12 = e21, which corresponds to equal node-speciﬁc ﬁxation probabilities (ρ1 = ρ2 = 1/N from Eqs. (32) and (33)) and equal death rates (d1 = d2 = 1/N). When the condition of uniform birth rates is relaxed, we find, by substituting Eqs. (32) and (33) into Eq. (3), the following expression for the overall fixation probability: ðe12 À e21 Þ B À d1 2 : Bðe12 þ e21 Þ r¼ 1 À 2 We discover that r ¼ 1, and consequently K = u, if (i) death rates are equal, d1 ¼ d2 ¼ B (Result 2 2 1), or (ii) node-speciﬁc ﬁxation probabilities are equal, r1 ¼ r2 ¼ 1 (Result 2). When death 2 rates are unequal and node-speciﬁc ﬁxation probabilities are unequal, there are cases in which r < 1 (K < u) and cases in which r > 1 (K > u). In particular, suppose that d1 > B. Then ρ is 2 2 2 less than 1 if e12 < e21 and greater than 1 if e12 > e21. 2 2 PLOS Computational Biology | DOI:10.1371/journal.pcbi.1004108 February 26, 2015 26 / 32 Asymmetric Spatial Structure Alters the Molecular Clock Rate Population size N = 3 We solve Eqs. (1) and (2) to obtain expressions for ρ1, ρ2 and ρ3: r1 ¼ e12 e13 þ e12 e23 þ e13 e32 ; e12 ðe13 þ e23 þ e31 Þ þ e13 ðe21 þ e32 Þ þ ðe21 þ e31 Þðe23 þ e32 Þ e21 e13 þ e21 e23 þ e23 e31 ; e12 ðe13 þ e23 þ e31 Þ þ e13 ðe21 þ e32 Þ þ ðe21 þ e31 Þðe23 þ e32 Þ e12 e31 þ e21 e32 þ e31 e32 : e12 ðe13 þ e23 þ e31 Þ þ e13 ðe21 þ e32 Þ þ ðe21 þ e31 Þðe23 þ e32 Þ ð34Þ r2 ¼ ð35Þ r3 ¼ ð36Þ Substituting Eqs. (34)–(36) into Eq. (3) gives an expression for the overall ﬁxation probability num r ¼ denom where num denom ¼ d1 ðe12 e13 þ e12 e23 þ e13 e32 Þ þ d2 ðe21 e13 þ e21 e23 þ e23 e31 Þ þ d3 ðe12 e31 þ e21 e32 þ e31 e32 Þ ¼ Bðe12 ðe13 þ e23 þ e31 Þ þ e13 ðe21 þ e32 Þ þ ðe21 þ e31 Þðe23 þ e32 ÞÞ: We factor to obtain an expression for num(Δρ) in terms of Ddi ¼ di À B for i 2 {1, 2, 3}, 3 ~ ~ ~ ~ numðDrÞ ¼ ðe12 Dd1 þ e21 Dd2 Þðd3 À b3 Þ À e32 Dd1 ðd1 À b1 Þ À e31 Dd2 ðd2 À b2 Þ: ð37Þ From Eq. (37) we see that r ¼ 1, and consequently K = u, in the case of equal death rates 3 (d1 ¼ d2 ¼ d3 ¼ B) and in the case of equivalent birth and death rate at each site (bi = di for i 2 3 {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 ~ Dd1 ¼ 0; Dd2 > 0; d2 > b2 ; and d3 < b3 then Eq. (37) yields num(Δρ) < 0 and consequently, r < 1 and K < u. On the other hand, if Dd1 ¼ 0; b2 > d2 > B ; and d3 > b3 then Eq. (37) yields 3 3 num(Δρ) > 0 and therefore, r > 1 and K > u. 3 Upstream-downstream populations We now turn to the upstream-downstream model introduced in the Results and in Fig. 6. Theorem 7. In the upstream-downstream model, ρ > 1/N, and consequently K > u, if and only if d" > d#. Proof. From Eq. (3) we obtain the following expression for the overall fixation probability r Á 1À N d r þ N # d# r# B " " " ! 1 1 1 1 þ N# d# r# À : ¼ þ N" d" r" À N B N N ¼ ð38Þ Since ﬁxation probabilities sum to one, and since the total population size is N = N"+N#, we have 1 1 N# r# À ¼ ÀN" r" À : ð39Þ N N PLOS Computational Biology | DOI:10.1371/journal.pcbi.1004108 February 26, 2015 27 / 32 Asymmetric Spatial Structure Alters the Molecular Clock Rate Substituting in Eq. (38) yields N" 1 1 : r" À r ¼ þ ðd" À d# Þ N N B ð40Þ It follows from Eq. (5) and from e! > e that ρ" > ρ#. Moreover, Eq. (39) implies that ρ"−1/N and ρ#−1/N have opposite signs, and since ρ" > ρ#, it follows that ρ"−1/N must be positive. Thus the second term on the right-hand side of Eq. (40) has the sign of d"−d#. We therefore conclude that ρ > 1/N, and consequently the molecular clock is accelerated relative to the well-mixed case (K > u), if and only if d" > d#. Acknowledgments We are grateful to Joshua B. Plotkin and David McCandlish for helpful comments. Author Contributions Conceived and designed the experiments: BA CS YD MAN. Performed the experiments: BA CS YD RCM CP. Analyzed the data: BA CS YD RCM CP. Wrote the paper: BA CS YD RCM CP MAN. References 1. 2. 3. 4. 5. 6. Zuckerkandl E, Pauling L (1965) Evolutionary divergence and convergence in proteins. Evolving genes and proteins 97: 97–166. doi: 10.1016/B978-1-4832-2734-4.50017-6 Kimura M (1968) Evolutionary rate at the molecular level. Nature 217: 624–626. doi: 10.1038/ 217624a0 PMID: 5637732 Ayala FJ (1999) Molecular clock mirages. BioEssays 21: 71–75. doi: 10.1002/(SICI)1521-1878 (199901)21:1%3C71::AID-BIES9%3E3.3.CO;2-2 PMID: 10070256 Bromham L, Penny D (2003) The modern molecular clock. Nature Reviews Genetics 4: 216–224. doi: 10.1038/nrg1020 PMID: 12610526 Pollak E (1982) The rate of mutant substitution in populations with overlapping generations. Genetical Research 40: 89–94. doi: 10.1017/S0016672300018930 PMID: 7141224 Balloux F, Lehmann L (2012) Substitution rates at neutral genes depend on population size under fluctuating demography and overlapping generations. Evolution 66: 605–611. doi: 10.1111/j.15585646.2011.01458.x PMID: 22276552 Lehmann L (2014) Stochastic demography and the neutral substitution rate in class-structured populations. Genetics 197: 351–360. doi: 10.1534/genetics.114.163345 PMID: 24594520 Ohta T, Kimura M (1971) On the constancy of the evolutionary rate of cistrons. Journal of Molecular Evolution 1: 18–25. doi: 10.1007/BF01659391 PMID: 4377445 Kimura M, Ohta T (1974) On some principles governing molecular evolution. Proceedings of the National Academy of Sciences 71: 2848–2852. doi: 10.1073/pnas.71.7.2848 Langley CH, Fitch WM (1974) An examination of the constancy of the rate of molecular evolution. Journal of Molecular Evolution 3: 161–177. doi: 10.1007/BF01797451 PMID: 4368400 Kimura M (1984) The Neutral Theory of Molecular Evolution. Cambridge University Press. Ayala FJ (1986) On the virtues and pitfalls of the molecular evolutionary clock. Journal of Heredity 77: 226–235. PMID: 3020121 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. doi: 10.1007/BF02603118 PMID: 3118047 Kumar S (2005) Molecular clocks: four decades of evolution. Nature Reviews Genetics 6: 654–662. doi: 10.1038/nrg1659 PMID: 16136655 Maruyama T (1970) On the fixation probability of mutant genes in a subdivided population. Genetical Research 15: 221–225. doi: 10.1017/S0016672300001543 PMID: 5480754 Lande R (1979) Effective deme sizes during long-term evolution estimated from rates of chromosomal rearrangement. Evolution 33: 234–251. doi: 10.2307/2407380 7. 8. 9. 10. 11. 12. 13. 14. 15. 16. PLOS Computational Biology | DOI:10.1371/journal.pcbi.1004108 February 26, 2015 28 / 32 Asymmetric Spatial Structure Alters the Molecular Clock Rate 17. 18. 19. 20. 21. Slatkin M (1981) Fixation probabilities and fixation times in a subdivided population. Evolution: 477– 488. doi: 10.2307/2408196 Nagylaki T (1982) Geographical invariance in population genetics. Journal of Theoretical Biology 99: 159–172. doi: 10.1016/0022-5193(82)90396-4 PMID: 7169795 Tachida H, Iizuka M (1991) Fixation probability in spatially changing environments. Genetical Research 58: 243–251. doi: 10.1017/S0016672300029992 PMID: 1802806 Barton N (1993) The probability of fixation of a favoured allele in a subdivided population. Genetical Research 62: 149–149. doi: 10.1017/S0016672300031748 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. PMID: 14704194 Whitlock MC (2003) Fixation probability and time in subdivided populations. Genetics 164: 767–779. PMID: 12807795 Lieberman E, Hauert C, Nowak MA (2005) Evolutionary dynamics on graphs. Nature 433: 312–316. doi: 10.1038/nature03204 PMID: 15662424 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. doi: 10.1098/rspa.2008.0058 Patwa Z, Wahl LM (2008) The fixation probability of beneficial mutations. Journal of The Royal Society Interface 5: 1279–1289. doi: 10.1098/rsif.2008.0248 Masuda N, Ohtsuki H (2009) Evolutionary dynamics and fixation probabilities in directed networks. New Journal of Physics 11: 033012. doi: 10.1088/1367-2630/11/3/033012 Houchmandzadeh B, Vallade M (2011) The fixation probability of a beneficial mutation in a geographically structured population. New Journal of Physics 13: 073020. doi: 10.1088/1367-2630/13/7/ 073020 Shakarian P, Roos P, Johnson A (2012) A review of evolutionary graph theory with applications to game theory. Biosystems 107: 66–80. doi: 10.1016/j.biosystems.2011.09.006 PMID: 22020107 Voorhees B (2013) Birth–death fixation probabilities for structured populations. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Science 469: 20120248. doi: 10.1098/ rspa.2012.0248 Constable GWA, McKane AJ (2014) Population genetics on islands connected by an arbitrary network: An analytic approach. Journal of Theoretical Biology 358: 149–165. doi: 10.1016/j.jtbi.2014.05. 033 PMID: 24882790 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. doi: 10.1098/rspa.2013.0730 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. doi: 10.1007/s00285-0120622-x PMID: 23179131 Liggett TM (1985) Particle Systems. Springer. Nowak MA, Bonhoeffer S, May RM (1994) More spatial games. International Journal of Bifurcation and Chaos 4: 33–56. doi: 10.1142/S0218127494000046 Antal T, Redner S, Sood V (2006) Evolutionary dynamics on degree-heterogeneous graphs. Physical Review Letters 96: 188104. doi: 10.1103/PhysRevLett.96.188104 PMID: 16712402 Nowak MA, May RM (1992) Evolutionary games and spatial chaos. Nature 359: 826–829. doi: 10. 1038/359826a0 Taylor P, Lillicrap T, Cownden D (2011) Inclusive fitness analysis on mathematical groups. Evolution 65: 849–859. doi: 10.1111/j.1558-5646.2010.01162.x PMID: 21044061 Wakano JY, Ohtsuki H, Kobayashi Y (2012) A mathematical description of the inclusive fitness theory. Theoretical Population Biology. Maciejewski W (2014) Reproductive value in graph-structured populations. Journal of Theoretical Biology 340: 285–293. doi: 10.1016/j.jtbi.2013.09.032 PMID: 24096097 Lundy IJ, Possingham HP (1998) Fixation probability of an allele in a subdivided population with asymmetric migration. Genetical Research 71: 237–245. doi: 10.1017/S0016672398003139 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. doi: 10.1016/j.stem.2013.08.001 PMID: 24035355 22. 23. 24. 25. 26. 27. 28. 29. 30. 31. 32. 33. 34. 35. 36. 37. 38. 39. 40. 41. PLOS Computational Biology | DOI:10.1371/journal.pcbi.1004108 February 26, 2015 29 / 32 Asymmetric Spatial Structure Alters the Molecular Clock Rate 42. 43. Lopez-Garcia C, Klein AM, Simons BD, Winton DJ (2010) Intestinal stem cell replacement follows a pattern of neutral drift. Science 330: 822–825. doi: 10.1126/science.1196236 PMID: 20929733 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. doi: 10.1016/j.cell.2010.09.016 PMID: 20887898 Vermeulen L, Snippert HJ (2014) Stem cell dynamics in homeostasis and cancer of the intestine. Nature Reviews Cancer 14: 468–480. doi: 10.1038/nrc3744 PMID: 24920463 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. doi: 10.1126/science. 1243148 PMID: 24264992 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. doi: 10.1016/j.stem.2012.09.009 PMID: 23040474 Leskovec J (2011). Stanford large network dataset collection. http://snap.stanford.edu/data/index. html URL http://snap.stanford.edu/data/index.html Durrett R, Levin S (1994) The importance of being discrete (and spatial). Theoretical Population Biology 46: 363–394. doi: 10.1006/tpbi.1994.1032 Nakamaru M, Matsuda H, Iwasa Y (1997) The evolution of cooperation in a lattice-structured population. Journal of Theoretical Biology 184: 65–81. doi: 10.1006/jtbi.1996.0243 PMID: 9039401 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. doi: 10.1006/jtbi.1998.0730 PMID: 9750181 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. doi: 10.1038/nature04605 PMID: 16724065 Nowak MA, Tarnita CE, Antal T (2010) Evolutionary dynamics in structured populations. Philosophical Transactions of the Royal Society B: Biological Sciences 365: 19–30. doi: 10.1098/rstb.2009.0215 Allen B, Nowak MA (2014) Games on graphs. EMS Surveys in Mathematical Sciences 1: 113–151. doi: 10.4171/EMSS/3 Débarre F, Hauert C, Doebeli M (2014) Social evolution in structured populations. Nature Communications 5: 4409. doi: 10.1038/ncomms4409 PMID: 24598979 Felsenstein J (1976) The theoretical population genetics of variable selection and migration. Annual Review of Genetics 10: 253–280. doi: 10.1146/annurev.ge.10.120176.001345 PMID: 797310 Hedrick PW, Ginevan ME, Ewing EP (1976) Genetic polymorphism in heterogeneous environments. Annual Review of Ecology and Systematics 7: 1–32. doi: 10.1146/annurev.es.07.110176.000245 Slatkin M (1987) Gene flow and the geographic structure of natural populations. Science 236: 787– 792. doi: 10.1126/science.3576198 PMID: 3576198 Cherry JL, Wakeley J (2003) A diffusion approximation for selection and drift in a subdivided population. Genetics 163: 421–428. PMID: 12586727 MacArthur RH, Wilson EO (1967) The Theory of Island Biogeography. Princeton, NJ: Princeton University Press. Doebeli M, Dieckmann U (2003) Speciation along environmental gradients. Nature 421: 259–264. doi: 10.1038/nature01274 PMID: 12529641 De Aguiar M, Baranger M, Baptestini E, Kaufman L, Bar-Yam Y (2009) Global patterns of speciation and diversity. Nature 460: 384–387. doi: 10.1038/nature08168 PMID: 19606148 Kawecki TJ, Holt RD (2002) Evolutionary consequences of asymmetric dispersal rates. The American Naturalist 160: 333–347. doi: 10.1086/341519 PMID: 18707443 Lenormand T (2002) Gene flow and the limits to natural selection. Trends in Ecology & Evolution 17: 183–189. doi: 10.1016/S0169-5347(02)02497-7 Garant D, Forde SE, Hendry AP (2007) The multifarious effects of dispersal and gene flow on contemporary adaptation. Functional Ecology 21: 434–443. doi: 10.1111/j.1365-2435.2006.01228.x Maciejewski W, Fu F, Hauert C (2014) Evolutionary game dynamics in populations with heterogenous structures. PLoS Computational Biology 10: e1003567. doi: 10.1371/journal.pcbi.1003567 PMID: 24762474 Nagylaki T (1978) Clines with asymmetric migration. Genetics 88: 813–827. PMID: 17248820 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. doi: 10.1016/j.tpb.2006.11.001 PMID: 17275866 44. 45. 46. 47. 48. 49. 50. 51. 52. 53. 54. 55. 56. 57. 58. 59. 60. 61. 62. 63. 64. 65. 66. 67. PLOS Computational Biology | DOI:10.1371/journal.pcbi.1004108 February 26, 2015 30 / 32 Asymmetric Spatial Structure Alters the Molecular Clock Rate 68. 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. doi: 10.1086/648311 PMID: 19860541 Tackenberg O, Poschlod P, Bonn S (2003) Assessment of wind dispersal potential in plant species. Ecological Monographs 73: 191–205. doi: 10.1890/0012-9615(2003)073%5B0191:AOWDPI%5D2. 0.CO;2 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. doi: 10.1126/science.1095210 PMID: 15155945 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. doi: 10.1111/j.1365-294X.2010.04842.x PMID: 20875065 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. doi: 10.1073/pnas.1100473108 PMID: 21876126 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. doi: 10.1073/pnas.0901652106 PMID: 19805178 O’keefe K, Ramakrishnan U, Van Tuinen M, Hadly EA (2009) Source–sink dynamics structure a common montane mammal. Molecular Ecology 18: 4775–4789. doi: 10.1111/j.1365-294X.2009.04382.x PMID: 19863718 Sun JX, Helgason A, Masson G, Ebenesersdóttir SS, Li H, et al. (2012) A direct characterization of human mutation based on microsatellites. Nature genetics 44: 1161–1165. doi: 10.1038/ng.2398 PMID: 22922873 Kong A, Frigge ML, Masson G, Besenbacher S, Sulem P, et al. (2012) Rate of de novo mutations and the importance of father/’s age to disease risk. Nature 488: 471–475. doi: 10.1038/nature11396 PMID: 22914163 Campbell CD, Eichler EE (2013) Properties and rates of germline mutations in humans. Trends in Genetics 29: 575–584. doi: 10.1016/j.tig.2013.04.005 PMID: 23684843 Knudson AG (2001) Two genetic hits (more or less) to cancer. Nature Reviews Cancer 1: 157–162. doi: 10.1038/35101031 PMID: 11905807 Attolini CSO, Michor F (2009) Evolutionary theory of cancer. Annals of the New York Academy of Sciences 1168: 23–51. doi: 10.1111/j.1749-6632.2009.04880.x PMID: 19566702 Vogelstein B, Papadopoulos N, Velculescu VE, Zhou S, Diaz LA, et al. (2013) Cancer genome landscapes. Science 339: 1546–1558. doi: 10.1126/science.1235122 PMID: 23539594 Goldie JH, Coldman AJ (1998) Drug resistance in cancer: mechanisms and models. Cambridge University Press. Komarova NL, Wodarz D (2005) Drug resistance in cancer: principles of emergence and prevention. Proceedings of the National Academy of Sciences 102: 9714–9719. doi: 10.1073/pnas.0501870102 PMID: 15980154 Michor F, Nowak MA, Iwasa Y (2006) Evolution of resistance to cancer therapy. Current Pharmaceutical Design 12: 261–271. doi: 10.2174/138161206775201956 PMID: 16454743 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. doi: 10.7554/eLife.00747 PMID: 23805382 Nowak MA, Michor F, Iwasa Y (2003) The linear process of somatic evolution. Proceedings of the National Academy of Sciences 100: 14966–14969. doi: 10.1073/pnas.2535419100 PMID: 14657359 Dingli D, Traulsen A, Pacheco JM (2007) Stochastic dynamics of hematopoietic tumor stem cells. Cell Cycle 6: 461–466. doi: 10.4161/cc.6.4.3853 PMID: 17329969 Bozic I, Nowak MA (2013) Unwanted evolution. Science 342: 938–939. doi: 10.1126/science. 1247887 PMID: 24264980 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. doi: 10.1002/embr.201337799 PMID: 24355609 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. doi: 10.1016/j.cell.2011. 08.030 PMID: 21854987 May RM, Lloyd AL (2001) Infection dynamics on scale-free networks. Physical Review E 64: 066112. doi: 10.1103/PhysRevE.64.066112 69. 70. 71. 72. 73. 74. 75. 76. 77. 78. 79. 80. 81. 82. 83. 84. 85. 86. 87. 88. 89. 90. PLOS Computational Biology | DOI:10.1371/journal.pcbi.1004108 February 26, 2015 31 / 32 Asymmetric Spatial Structure Alters the Molecular Clock Rate 91. 92. 93. Watts DJ (2002) A simple model of global cascades on random networks. Proceedings of the National Academy of Sciences 99: 5766–5771. doi: 10.1073/pnas.082090499 PMID: 16578874 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. doi: 10.1056/NEJMsa066082 PMID: 17652652 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. doi: 10.1073/pnas. 0706851105 PMID: 18216267 Castellano C, Fortunato S, Loreto V (2009) Statistical physics of social dynamics. Reviews of Modern Physics 81: 591. doi: 10.1103/RevModPhys.81.591 Centola D (2010) The spread of behavior in an online social network experiment. Science 329: 1194– 1197. doi: 10.1126/science.1185231 PMID: 20813952 Hill AL, Rand DG, Nowak MA, Christakis NA (2010) Infectious disease modeling of social contagion in networks. PLoS Computational Biology 6: e1000968. doi: 10.1371/journal.pcbi.1000968 PMID: 21079667 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. doi: 10. 1145/2187836.2187907 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. doi: 10.1007/s10051-002-8996-y 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. doi: 10.1016/j.physa. 2006.07.017 Skyrms B, Pemantle R (2009) A dynamic model of social network formation. In: Adaptive Networks, Springer. pp. 231–251. Rosvall M, Bergstrom CT (2010) Mapping change in large networks. PloS ONE 5: e8694. doi: 10. 1371/journal.pone.0008694 PMID: 20111700 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. doi: 10. 1073/pnas.1108243108 PMID: 22084103 Wardil L, Hauert C (2014) Origin and structure of dynamic cooperative networks. Scientific Reports 4. doi: 10.1038/srep05725 PMID: 25030202 Baxter GJ, Blythe RA, McKane AJ (2008) Fixation and consensus times on a network: A unified approach. Physical Review Letters 101: 258701. doi: 10.1103/PhysRevLett.101.258701 PMID: 19113759 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. doi: 10.1098/rspb.2013.0211 PMID: 23677339 Tarnita CE, Taylor PD (2014) Measures of relative fitness of social behaviors in finite structured population models. The American Naturalist 184: 477–488. doi: 10.1086/677924 PMID: 25226183 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. doi: 10.1098/rspa.2010.0067 Shakarian P, Roos P, Moores G (2013) A novel analytical method for evolutionary graph theory problems. Biosystems 111: 136–144. doi: 10.1016/j.biosystems.2013.01.006 PMID: 23353025 94. 95. 96. 97. 98. 99. 100. 101. 102. 103. 104. 105. 106. 107. 108. PLOS Computational Biology | DOI:10.1371/journal.pcbi.1004108 February 26, 2015 32 / 32