Fitness dependence of the fixation-time distribution for evolutionary dynamics on graphs

TL;DR

This paper analyzes how mutant fitness influences the distribution of fixation times in evolutionary graph models, revealing distinct behaviors on different network structures and identifying conditions for universal and skewed distributions.

Contribution

It provides exact solutions for fixation-time distributions under varying fitness levels on complete graphs and rings, highlighting their dependence and discontinuities.

Findings

Fixation times are normally distributed on rings with non-neutral fitness.

On complete graphs, fixation times follow a fitness-weighted convolution of Gumbel distributions.

Neutral fitness leads to highly skewed fixation-time distributions on both networks.

Abstract

Evolutionary graph theory models the effects of natural selection and random drift on structured populations of competing mutant and non-mutant individuals. Recent studies have found that fixation times in such systems often have right-skewed distributions. Little is known, however, about how these distributions and their skew depend on mutant fitness. Here we calculate the fitness dependence of the fixation-time distribution for the Moran Birth-death process in populations modeled by two extreme networks: the complete graph and the one-dimensional ring lattice, obtaining exact solutions in the limit of large network size. We find that with non-neutral fitness, the Moran process on the ring has normally distributed fixation times, independent of the relative fitness of mutants and non-mutants. In contrast, on the complete graph, the fixation-time distribution is a fitness-weighted convolution of two Gumbel distributions. When fitness is neutral the fixation-time distribution jumps discontinuously and becomes highly skewed on both the complete graph and the ring. Even on these simple networks, the fixation-time distribution exhibits rich fitness dependence, with discontinuities and regions of universality. Extensions of our results to two-fitness Moran models, times to partial fixation, and evolution on random networks are discussed.

Tables

Table 1. Table 1: Asymptotic fixation-time statistics for the Moran Birth-death and Death-birth processes on the complete graph and the 1D lattice. Together with the mean and variance, the standardized distributions give a complete statistical description of the fixation time. The mean and variance given are to leading order in N 𝑁 N for each case.

		Asymptotic Fixation-Time Statistics
Network	Fitness Level	Mean	Variance	Standardized Distribution
1D Lattice	$r=1$	$\frac{1}{6}{N}^3$	$\frac{1}{90}{N}^6$	Highly Skewed [Eqs. (11) & (12)]
	$r>1$	$\frac{r+1}{2\left(r-1\right)}{N}^2$	$\frac{\left(r+1\right)\left(r^2+r+1\right)}{3{\left(r-1\right)}^3}{N}^3$	$𝒩(0,1)$
Complete Graph	$r=1$	$N^2$	$\left(\frac{{\pi }^2}{3}-3\right)N^4$	Highly Skewed [Eqs. (IV.1) & (24)]
	$r>1$	$\frac{r+1}{r-1}N\log N$	$\frac{{\pi }^2{\left(r+1\right)}^2}{6{\left(r-1\right)}^2}{N}^2$	$G+r{G}^{\prime }$

Equations

$$\dot{\mathbf{p}}(t)=\Omega \cdot \mathbf{p}(t),$$

$${\Omega }_{mn}=b_n{\delta }_{m,n+1}+d_n{\delta }_{m,n-1}-(b_n+d_n){\delta }_{m,n},$$

$$f_T(t)=\sum _{j=1}^{N-1}\left(\prod _{k=1,k\ne j}^{N-1}\frac{{\lambda }_k}{{\lambda }_k-{\lambda }_j}\right){\lambda }_j{e}^{-{\lambda }_{j}t}.$$

$$\kappa_{n}(N)=(n-1)!\left(\sum_{m=1}^{N-1}\frac{1}{\lambda_{m}^{n}}\right)\Bigg{/}\left(\sum_{m=1}^{N-1}\frac{1}{\lambda_{m}^{2}}\right)^{n/2},$$

$$[{\Omega }_{\text{RW}}{]}_{nm}=\frac{r}{1+r}{\delta }_{m,n+1}+\frac{1}{1+r}{\delta }_{m,n-1}-{\delta }_{m,n}.$$

$$\varphi (\omega ):=E[\exp (i\omega T)]=1{\tilde{\Omega }}_{\text{tr}}(i\omega +{\tilde{\Omega }}_{\text{tr}}{)}^{-1}{\mathbf{p}}_{\text{tr}}(0).$$

$$E[T^n]=(-1{)}^{n}n!1{\tilde{\Omega }}_{\text{tr}}^{-n}{\mathbf{p}}_{\text{tr}}(0),$$

$${\kappa }_n(N)=\frac{\sum _{i_1,i_2,\dots ,i_n=1}^{N-1}\frac{w_{i_1{i}_2\cdots i_n}^n(r,N|m_0)}{(b_{i_1}+d_{i_1})(b_{i_2}+d_{i_2})\cdots (b_{i_n}+d_{i_n})}}{{\left(\sum _{i,j=1}^{N-1}\frac{w_{ij}^2(r,N|m_0)}{(b_i+d_i)(b_j+d_j)}\right)}^{n/2}},$$

$$b_m=\frac{r}{rm+N-m},d_m=\frac{1}{rm+N-m},$$

$$\lambda_{m}=\frac{2}{N}-\frac{2}{N}\cos\bigg{(}\frac{m\pi}{N}\bigg{)},\quad m=1,2,\dots,N-1.$$

$${\kappa }_n=(n-1)!\frac{\zeta (2n)}{\zeta (4{)}^{n/2}},$$

$$\varphi (\omega )=e^{-\sqrt{\frac{5}{2}}\omega }\Gamma \left(1-\frac{90^{1/4}\sqrt{\omega }}{\pi }\right)\Gamma \left(1+\frac{90^{1/4}\sqrt{\omega }}{\pi }\right).$$

$$\lambda_{m}=1-\frac{2\sqrt{r}}{1+r}\cos\bigg{(}\frac{m\pi}{N}\bigg{)},\quad m=1,2,\dots,N-1.$$

$$S_n\approx \frac{N}{\pi }{\int }_0^{\pi }\frac{1}{(1-2\sqrt{r}/(1+r)\cos x{)}^n}\text{d}x.$$

$${\kappa }_n\sim \frac{1}{N^{(n-2)/2}}\stackrel{N\to \infty }{\to }0,n\ge 3.$$

$${\kappa }_3\approx \frac{2+2r(r+4)}{(r+1)\sqrt{(r^2-1)}}\frac{1}{\sqrt{N}}.$$

$$\begin{split}\Delta&=\bigg{|}\frac{(1-2\sqrt{r}/(1+r))^{n}}{(1-2\sqrt{r}/(1+r)\cos(\pi/N))^{n}}-1\bigg{|}\\ &=\frac{\sqrt{r}n\pi^{2}}{(\sqrt{r}-1)^{2}N^{2}}+\mathcal{O}(1/N^{4}).\end{split}$$

$$\kappa_{3}^{\infty}=2\left(\sum_{m=1}^{N-1}m^{3}\right)\bigg{/}\left(\sum_{m=1}^{N-1}m^{2}\right)^{3/2}\approx\frac{3\sqrt{3}}{2}\frac{1}{\sqrt{N}},$$

$${\tilde{\kappa }}_3={\kappa }_3\left[r^{-q}+\frac{3\sqrt{3}}{4}(1-r^{-q})\right]$$

$$b_m=\frac{rm}{rm+N-m}\cdot \frac{N-m}{N-1},$$

$$d_m=\frac{N-m}{rm+N-m}\cdot \frac{m}{N-1}.$$

$${\lambda }_m=\frac{m(m+1)}{N(N-1)},m=1,2,\dots N-1.$$

$${\kappa }_n$$

$$\begin{array}{rl}\sigma f_T(\sigma t+\mu )=c_{\sigma }\sum _{j=1}^{\infty }(-& 1{)}^{j+1}j(j+1)(2j+1)\\ & \times \exp \left[j(j+1)(c_{\sigma }t+1)\right],\end{array}$$

$$f(t)={\beta }^{-1}{e}^{-(t-\alpha )/\beta }\exp (-e^{-(t-\alpha )/\beta }).$$

$$\frac{T-\mu }{\sigma }\stackrel{d}{\to }\frac{G+r{G}^{\prime }}{\sqrt{1+r^2}},$$

$$T\approx \sum _{m=1}^{N-1}\mathcal{E}(b_m+d_m).$$

$${\kappa }_n=\frac{1+r^n}{(1+r^2{)}^{n/2}}\times \frac{(n-1)!\zeta (n)}{\zeta (2{)}^{n/2}},$$

$$[{\Omega }_{\text{RW}}{]}_{nm}=\frac{r}{1+r}{\delta }_{m,n+1}+\frac{1}{1+r}{\delta }_{m,n-1}-{\delta }_{m,n}.$$

$${\tilde{\Omega }}_{\text{tr}}=S{\Omega }_{\text{RW}}{S}^{-1}D,$$

Full text

I Introduction

Reproducing populations undergo evolutionary dynamics. Mutations can endow individuals with a fitness advantage, allowing them to reproduce more quickly and outcompete non-mutant individuals Nowak (2006). Two natural questions arise: If a single mutant individual is introduced into a population, what is the probability that the mutant lineage will spread and ultimately take over the population (an outcome known as fixation)? And if fixation occurs, how much time does it take?

These questions have been addressed, in part, by evolutionary graph theory, which studies evolutionary dynamics in structured populations. Thanks to this approach, fixation probabilities are now well understood for various models on various networks Maruyama (1974a, b); Slatkin (1981); Lieberman et al. (2005); Antal and Scheuring (2006); Houchmandzadeh and Vallade (2011); Díaz et al. (2014); Kaveh et al. (2015); Jamieson-Lane and Hauert (2015); Altrock et al. (2017); Tkadlec et al. (2019). Less is known about fixation times. Given a model of evolutionary dynamics, one would like to predict the mean, variance, and ideally the full distribution of its fixation times.

Of these quantities, the mean is the best understood. Numerical and analytical results exist for mean fixation times on both deterministic Kimura (1980); Slatkin (1981); Antal and Scheuring (2006); Altrock and Traulsen (2009); Frean et al. (2013); Askari and Samani (2015); Altrock et al. (2017); Askari et al. (2017); Tkadlec et al. (2019) and random Askari and Samani (2015); Askari et al. (2017); Farhang-Sardroodi et al. (2017); Hajihashemi and Aghababaei Samani (2019) networks. Yet although mean fixation times are important to study, the information they provide can be misleading, because fixation-time distributions tend to be broad and skewed and hence are not well characterized by their means alone Dingli et al. (2007); Altrock et al. (2011); Ashcroft et al. (2015); Altrock et al. (2017); Ying et al. (2018). Initial analytical results have determined the asymptotic fixation-time distribution for several simple networks, but only when the relative fitness of the mutants is infinite Aldous (2013); Ottino-Löffler et al. (2017, 2017). For other values of the relative fitness, almost nothing is known. Preliminary results suggest that at neutral fitness (when mutants and non-mutants are equally fit), the fixation-time distribution becomes highly right-skewed Ottino-Löffler et al. (2017).

In this paper we investigate the full fitness dependence of fixation-time distributions for the Moran process Moran (1958a, b), a simple model of evolutionary dynamics. In the limit of large network size, we derive asymptotically exact results for the fixation-time distribution and its skew for two network structures at opposite ends of the connectivity spectrum: the complete graph, in which every individual interacts with every other individual; and the one-dimensional ring lattice, in which each individual interacts only with its nearest neighbors on a ring.

The specific model we consider is the Moran Birth-death (Bd) process 111We use the convention that capital letters designate a fitness dependent step in the Moran process (e.g. for the Bd process nodes give birth at a rate proportional to their fitness, but die with uniform probability). See Ref. (Ottino-Löffler et al., 2017, Box 2) for a detailed explanation of this nomenclature., defined as follows. On each node of the network there is an individual, either mutant or non-mutant. The mutants have a fitness level $r$, which designates their relative reproduction rate compared to non-mutants. When $r>1$, the mutants have a fitness advantage, whereas when $r=1$ they have neutral fitness. At each time step we choose a node at random, with probability proportional to its fitness, and choose one of its neighbors with uniform probability. The first individual gives birth to an offspring of the same type. That offspring replaces the neighbor, which dies. The model population is updated until either the mutant lineage takes over (in which case fixation occurs) or the mutant lineage goes extinct (a case not considered here).

As mentioned above, the distribution of fixation times is often skewed. The skew emerges from the stochastic competition between mutants and non-mutants through multiple mechanisms. For instance, when the mutants have neutral fitness the process resembles an unbiased random walk. We find that the asymptotic fixation-time distribution for a simple random walk is only skewed when the walk is unbiased. The lack of bias allows for occasional long recurrent excursions (that are suppressed in biased walks) during successful runs to fixation. The fixation-time distribution is strongly skewed because there are many ways to execute such walks that are much longer than usual, but comparably few ways for mutants to sweep through the population much faster than usual.

Depending on network structure, the fixation-time skew can also come from a second, completely separate mechanism, which involves characteristic slowdowns that arise because individuals do not discriminate between mutants and non-mutants during the replacement step of the Moran process. For example, when very few non-mutants remain, the mutants can waste time replacing each other. These slowdowns are reminiscent of those seen in a classic problem from probability theory, the coupon collector’s problem, which asks: How long does it take to complete a collection of $N$ distinct coupons if a random coupon is received at each time step? The intuition for the long slowdowns is clear: when nearly all the coupons have been collected, it can take an exasperatingly long time to collect the final few, because one keeps acquiring coupons that one already has. The problem was first solved by Erdős and Rényi, who proved that for large $N$, the time to complete the collection has a Gumbel distribution Erdős and Rényi (1961). In fact, for evolutionary processes with infinite fitness there exists an exact mapping onto coupon collection Ottino-Löffler et al. (2017, 2017). Remarkably, while this correspondence breaks down for finite fitness, the coupon collection heuristic still allows us to predict correct asymptotic fixation-time distributions for non-neutral fitness.

In the following sections we show that for $N\gg 1$, the neutral-fitness Moran process on the complete graph and the one-dimensional ring lattice has highly skewed fixation-time distributions, and we solve for their cumulants exactly. For non-neutral fitness the fixation-time distribution is normal on the lattice and a weighted convolution of Gumbel distributions on the complete graph. These results are novel; apart from the infinite fitness limit and some partial results at neutral fitness (noted below), the fitness dependence of these distributions was previously unknown.

We begin by developing a general framework for computing fixation-time distributions and cumulants of birth-death Markov chains, and then apply it to the Moran process to prove the results above. We also consider the effects of truncation on the process and examine how long it takes to reach partial, rather than complete, fixation. The fixation-time distributions have rich dependence on the fitness level and the degree of truncation, with both discontinuities and regions of universality. To conclude, we discuss extensions of our results to two-fitness Moran models and to more complicated network topologies.

II General Theory for Birth-Death Markov Processes

For simplicity, we restrict attention to network topologies and initial mutant populations for which the probability of adding or removing a mutant in a given time step depends only on the number of existing mutants, not on where the mutants are located on the network. The state of the system can therefore be defined in terms of the number of mutants, $m=0,1,\dots,N$, where $N$ is the total number of nodes on the network. The Moran process is then a birth-death Markov chain with $N+1$ states, transition probabilities $b_{m}$ and $d_{m}$ determined by the network structure, and absorbing boundaries at $m=0$ and $m=N$. In this section we review several general analytical results for absorbing birth-death Markov chains, explaining how they apply to fixation times in evolutionary dynamics. We also develop an approach, which we call visit statistics, that enables analytical estimation of the asymptotic fixation time cumulants.

On more complicated networks, the probability of adding or removing a mutant depends on the configuration of existing mutants. For some of these networks, however, the transition probabilities can be accurately estimated using a mean-field approximation Hajihashemi and Aghababaei Samani (2019); Ying et al. (2018); Ottino-Löffler et al. (2017, 2017). Then, to a good approximation, the results below apply to such networks as well.

II.1 Eigendecomposition of the birth-death process

Assuming a continuous-time process, the state of the Markov chain described above evolves according to the master equation,

[TABLE]

where $\mathbf{p}(t)$ is the probability of occupying each state of the system at time $t$ and $\Omega$ is the transition rate matrix, with columns summing to zero. In terms of the transition probabilities $b_{m}$ and $d_{m}$, the entries of $\Omega$ are

[TABLE]

where $m$ and $n$ run from [math] to $N$, $\delta_{m,n}$ is the Kronecker delta, and $b_{0}=d_{0}=b_{N}=d_{N}=0$. The final condition guarantees the system has absorbing boundaries with stationary states $p_{m}=\delta_{m,0}$ and $p_{m}=\delta_{m,N}$ when the population is homogeneous. Thus we can decompose the transition matrix into stationary and transient parts, defining the transient part $\Omega_{\text{tr}}$ as in Eq. (2), but with $m,\,n=1,\dots,N-1$. The transient transition matrix acts on the transient states of the system, denoted $\mathbf{p}_{\text{tr}}(t)$. The eigenvalues of $\Omega_{\text{tr}}$ are real and strictly negative, since probability flows away from these states toward the absorbing boundaries. To ease notation in the following discussion and later applications, we shall refer to the positive eigenvalues of $-\Omega_{\text{tr}}$ as the eigenvalues of the transition matrix, denoted $\lambda_{m}$, where $m=1,\dots,N-1$.

From the perspective of Markov chains, the fixation time $T$ is the time required for first passage to state $m=N$, given $m_{0}$ initial mutants, $p_{m}(0)=\delta_{m,m_{0}}$. At time $t$, the probability that state $N$ has been reached (i.e., the cumulative distribution function for the first passage times) is simply $\varphi_{m_{0}}^{-1}p_{N}(t)$, where $\varphi_{m_{0}}$ is the fixation probability given $m_{0}$ initial mutants. The distribution of first passage times is therefore $\varphi_{m_{0}}^{-1}\dot{p}_{N}(t)=\varphi_{m_{0}}^{-1}b_{N-1}p_{N-1}(t)$. Since we normalize by the fixation probability, this is precisely the fixation-time distribution conditioned on reaching $N$.

The solution to the transient master equation is the matrix exponential $\mathbf{p}_{\text{tr}}(t)=\exp(\Omega_{\text{tr}}t)\cdot\mathbf{p}_{\text{tr}}(0)$, yielding a fixation-time distribution $\varphi_{m_{0}}^{-1}b_{N-1}[\exp(\Omega_{\text{tr}}t)\cdot\mathbf{p}_{\text{tr}}(0)]_{N-1}$ Asmussen (2003). If we assume one initial mutant $m_{0}=1$ this becomes $\varphi_{1}^{-1}b_{N-1}[\exp(\Omega_{\text{tr}}t)]_{N-1,1}$. The matrix exponential can be evaluated in terms of the eigenvalues $\lambda_{m}$ by taking a Fourier (or Laplace) transform (for details, see Ref. Ashcroft et al. (2015)). For a single initial mutant, the result is that the fixation time $T$ has a distribution $f_{T}(t)$ given by

[TABLE]

This formula holds as long as the eigenvalues $\lambda_{m}$ are distinct, which for birth-death Markov chains occurs when $b_{m}$ and $d_{m}$ are non-zero (except at the absorbing boundaries) Keilson (1979). Generalizations of this result for arbitrarily many initial mutants have also recently been derived, in terms of eigenvalues of the transition matrix and certain sub-matrices Ashcroft et al. (2015).

The distribution in Eq. (3) is exactly that corresponding to a sum of exponential random variables with rate parameters $\lambda_{m}$. The corresponding cumulants equal $(n-1)!\sum_{m=1}^{N-1}(\lambda_{m})^{-n}$. As our primary interest is the asymptotic shape of the distribution, we normalize $T$ to zero mean and unit variance and study $(T-\mu)/\sigma$, where $\mu$ and $\sigma$ denote the mean and standard deviation of $T$. The standardized distribution is then given by $\sigma f_{T}(\sigma t+\mu)$. The rescaled fixation time has cumulants

[TABLE]

which, for many systems including those considered below, are finite as $N\rightarrow\infty$. When the limit exists, we define the asymptotic cumulants by $\kappa_{n}=\lim_{N\rightarrow\infty}\kappa_{n}(N)$. In particular, because we have standardized our distribution, the third cumulant $\kappa_{3}$ is the skew. In practice the limit $N\rightarrow\infty$ is taken by computing the leading asymptotic behavior of both the numerator and denominator in Eq. (4). As we will see below the scaling of these terms with $N$ depends on both the population network structure and the mutant fitness (see also asymptotic analysis in Supplemental Material, Sections S3 & S4 SM ). This approach allows us to characterize the asymptotic shape of the fixation-time distribution in terms of the constants $\kappa_{n}$. Since $\lambda_{m}>0$, it is clear from this expression that, for finite $N$, the skew and all higher order cumulants must be positive, in agreement with results for random walks with non-uniform bias Bakhtin (2018). As $N\rightarrow\infty$ this is not necessarily true; in some cases the cumulants vanish.

The eigendecomposition gives the fixation-time distribution and cumulants in terms of the non-zero eigenvalues of the transition matrix. In general the eigenvalues must be found numerically, but in cases where they have a closed form expression the asymptotic form of the cumulants and distribution can often be obtained exactly.

II.2 Analytical cumulant calculation: Visit statistics

In this section we develop machinery to compute the cumulants of the fixation time analytically without relying on matrix eigenvalues. For this analysis, we specialize to cases where $b_{m}/d_{m}=r$ for all $m$, relevant for the Moran processes considered below. These processes can be thought of as biased random walks overlaid with non-constant waiting times at each state.

It is helpful to consider the Markov chain conditioned on hitting $N$, with new transition probabilities $\tilde{b}_{m}$ and $\tilde{d}_{m}$ so that the fixation probability $\varphi_{m_{0}}=1$. If $X_{t}$ is the state of the system at time $t$, then $\tilde{b}_{m}=\mathcal{P}(X_{t}=m\rightarrow X_{t+1}=m+1|X_{\infty}=N)$ with $\tilde{d}_{m}$ defined analogously. We derive explicit expressions for $\tilde{b}_{m}$ and $\tilde{d}_{m}$ in Supplemental Material, Section S1 SM . Conditioning is equivalent to a similarity transformation on the transient part of the transition matrix: $\tilde{\Omega}_{\text{tr}}=S\,\Omega_{\text{tr}}\,S^{-1}$, where $S$ is diagonal with $S_{mm}=1-1/r^{m}$. Furthermore, since $b_{m}/d_{m}=r$, we can decompose $\Omega_{\text{tr}}=\Omega_{\text{RW}}D$, where $D$ is a diagonal matrix, $D_{mm}=b_{m}+d_{m}$, that encodes the time spent in each state and $\Omega_{\text{RW}}$ is the transition matrix for a random walk with uniform bias,

[TABLE]

Applying the results of the previous section and using the fact that the columns of $\Omega$ sum to zero, we can write there fixation-time distribution of the conditioned Markov chain as $f_{T}(t)=-\mathbf{1}\tilde{\Omega}_{\text{tr}}\exp(\tilde{\Omega}_{\text{tr}}t)\mathbf{p}_{\text{tr}}(0)$, where $\mathbf{1}$ is the row vector containing all ones. This distribution has characteristic function Asmussen (2003)

[TABLE]

and the derivatives $(-i)^{n}\phi^{(n)}(0)$ give the moments of $T$

[TABLE]

in terms of $\tilde{\Omega}_{\text{tr}}^{-1}=D^{-1}S\Omega_{\text{RW}}^{-1}S^{-1}$. This inverse has a nice analytical form because $S$ and $D$ are diagonal and $\Omega_{\text{RW}}$ is tridiagonal Toeplitz. We call this approach visit statistics because the elements $V_{ij}$ of $V=-S\Omega_{\text{RW}}^{-1}S^{-1}$ encode the average number of visits to state $i$ starting from state $j$.

Each power of $\tilde{\Omega}_{\text{tr}}$ in Eq. (7) produces products of $(b_{i}+d_{i})$ that arise in linear combinations determined by the visit numbers $V_{ij}$. Therefore, the cumulants of the fixation time have the general form

[TABLE]

where $w_{i_{1}i_{2}\cdots i_{n}}^{n}(r,N|m_{0})$ are weighting factors based on the visit statistics of the biased random walk, given the initial number of mutants $m_{0}$. In what follows, we always assume $m_{0}=1$ and suppress the dependence of the weighting factors on initial condition, writing $w_{i_{1}i_{2}\cdots i_{n}}^{n}(r,N)$ instead. A detailed derivation of Eq. (8) and explicit expressions for $w_{ij}^{2}(r,N)$ and $w_{ijk}^{3}(r,N)$ are given in the Appendix below.

To the best of our knowledge this representation of the fixation-time cumulants has not been previously derived, although a similar approach was recently used to compute mean fixation times for evolutionary dynamics on complex networks Hajihashemi and Aghababaei Samani (2019). This expression is equivalent to the well-known recurrence relations for absorption-time moments of birth-death processes Altrock et al. (2011); Goel and Richter-Dyn (1974) but is easier to handle asymptotically, and can be useful even without explicit expressions for $w_{i_{1}i_{2}\cdots i_{n}}^{n}(r,N)$. Estimating the sums in Eq. (8) allows us to compute the asymptotic fixation time cumulants exactly.

II.3 Recurrence relation for fixation-time moments

Evaluation of the eigenvalues of the transition matrix for large systems can be computationally expensive, with the best algorithms having run times quadratic in matrix size. Numerical evaluation of the expression given in Eq. (8) is even worse, as it requires summing $\mathcal{O}(N^{n})$ elements. If only a finite number of fixation time cumulants (and not the full distribution) are desired, there are better numerical approaches. Using standard methods from probability theory Keilson (1965), we derive a recurrence relation that allows numerical moment computation with run time linear in system size $N$. For completeness we provide the full derivation of the reccurence for the fixation-time skew in Supplemental Material, Section S2 SM .

II.4 Equivalence between advantageous and disadvantageous mutations

In the following applications, we will generally speak of the mutants as having a fitness advantage, designated by the parameter $r>1$. Our results, however, can be immediately extended to disadvantageous mutations. In particular, the fixation-time distributions (conditioned on fixation occurring) for mutants of fitness $r$ and $1/r$ are identical. When a mutant with fitness $1/r$ is introduced into the population (and eventually reaches fixation), the non-mutants are $r$ times as fit as the mutants. Therefore, this system is equivalent to another system that starts with $N-1$ fitness $r$ mutants which eventually die out (the mutants in the former system are the non-mutants in the latter). It has been shown that the times to go from one initial mutant to fixation ($m=1\rightarrow m=N$) and from $N-1$ initial mutants to extinction ($m=N-1\rightarrow m=0$) have identical distribution Ashcroft et al. (2015). Thus indeed, the conditioned fixation-time distributions are identical for mutants of fitness $r$ and $1/r$. Of course the fixation probability is very different in the two cases: for the disadvantageous mutations it approaches 0 for large $N$ Lieberman et al. (2005).

III One-Dimensional Lattice

We now specialize to Moran Birth-death (Bd) processes, starting with the one-dimensional (1D) lattice. We assume periodic boundary conditions, so that the $N$ nodes form a ring. The mutants have relative fitness $r$, meaning they give birth $r$ times faster, on average, than non-mutants do.

Starting from one mutant, suppose that at some later time $m$ of the $N$ nodes are mutants. On the 1D lattice, the population of mutants always forms a connected arc, with two mutants at the endpoints of the arc. Therefore, the probability $b_{m}$ of increasing the mutant population by one in the next time step is the probability of choosing a mutant node at an endpoint to give birth, namely $2r/(rm+N-m)$, times the probability $1/2$ that the neighboring node to be replaced is not itself a mutant. (The latter probability equals 1/2 because there are two neighbors to choose for replacement: a mutant neighbor on the interior of the arc and a non-mutant neighbor on the exterior. Only the second of these choices produces an increase in the number of mutants.) Multiplying these probabilities together we obtain

[TABLE]

where the probability $d_{m}$ of decreasing the mutant population by one is found by similar reasoning. Note that this derivation fails for $m=1$ ($m=N-1$) when the arc of mutants (non-mutants) contains only one node, but one can check Eq. (9) still holds for these cases. These quantities play the role of transition probabilities in a Markov transition matrix. The next step is to find the eigenvalues of that matrix.

III.1 Neutral fitness

First we work out the eigenvalues for the case of neutral fitness, $r=1$. In this case, the transition probabilities are equal, $b_{m}=d_{m}=1/N$, and independent of $m$. Therefore, the Moran process is simply a random walk, with events occurring at a rate of $2/N$ per time step. The associated transition matrix is tridiagonal Toeplitz, which has eigenvalues given by

[TABLE]

Applying Eq. (4) and computing the leading asymptotic form of the sums $S_{n}=\sum_{m=1}^{N-1}(\lambda_{m})^{-n}$ (see Supplemental Material, Section S3 SM ), we find that as $N\rightarrow\infty$, the fixation-time distribution has cumulants

[TABLE]

where $\zeta$ denotes the Riemann zeta function. In particular, the skew $\kappa_{3}=4\sqrt{10}/7\approx 1.807$, as previously calculated by Ottino-Löffler et al. Ottino-Löffler et al. (2017) via martingale methods. The other cumulants (and characteristic function below) haven’t previously been computed for the Bd process on the 1D lattice. The largeness of the skew stems from the recurrent property of the random walk. As $N\rightarrow\infty$, long walks with large fixation times become common and the system revisits each state infinitely often Durrett (2010).

Knowledge of the cumulants allows us to obtain the exact characteristic function of the fixation-time distribution:

[TABLE]

Although we cannot find a simple expression for the distribution itself, we can efficiently evaluate it by taking the inverse Fourier transform of the characteristic function numerically. Figure 1(a) shows that the predicted fixation-time distribution agrees well with simulations.

III.2 Non-neutral fitness

Next, consider $r\neq 1$ with the transition probabilities given by Eq. (9). Then the eigenvalues of the transition matrix are no longer expressible in closed form. If $r$ is not too large, however, the probabilities $b_{m}$ and $d_{m}$ do not vary dramatically with $m$, the number of mutants. In particular, $b_{m}\sim 1/N$ for all $m$ when $N$ is large. Therefore, as a first approximation we treat the Bd process on a 1D lattice as a biased random walk with $b_{m}=r/(1+r)$ and $d_{m}=1/(1+r)$. The eigenvalues of the corresponding transition matrix are

[TABLE]

The cumulants again involve sums $S_{n}=\sum_{m=1}^{N-1}(\lambda_{m})^{-n}$, which can be approximated in the limit $N\rightarrow\infty$ by,

[TABLE]

Since the integral is independent of $N$ and converges for $r\neq 1$, each of the sums scales linearly: $S_{n}\sim N$. Thus, using Eq. (4), we see that all cumulants past second order approach 0,

[TABLE]

Hence the fixation-time distribution is asymptotically normal, independent of fitness level.

By evaluating the integrals in Eq. (14), we can more precisely compute the scaling of the cumulants. For the skew we find

[TABLE]

The integral approximation becomes accurate when the first term in the sums $S_{n}$ becomes close to the value of the integrand evaluated at the lower bound ($x=0$). The fractional difference between these quantities is

[TABLE]

Then we have $\Delta\ll 1$ when $N\gg N_{c}$ where $N_{c}\approx 2\pi\sqrt{n}/(r-1)$ (assuming $r$ is near 1). For the skew, we require the sums with $n=2$ and $3$, giving $N_{c}\approx 10/(r-1)$.

The above calculation fails for $r\gg 1$, because when $r=\infty$ the transition probabilities $b_{m}=1/m$ have different asymptotic behavior as $N\rightarrow\infty$. In particular, more time is spent waiting at states with large $m$. The process still has normally distributed fixation times Ottino-Löffler et al. (2017), but the skew becomes

[TABLE]

for large $N$. Notice that the coefficient is different from that given by the infinite-$r$ limit of Eq. (16), $\kappa_{3}\approx 2/\sqrt{N}$. We conjecture that there is a smooth crossover between these two scaling laws with the true skew given approximately by

[TABLE]

for some exponent $q$, where $\kappa_{3}$ is the skew given in Eq. (16). For small $r$ this ansatz has skew similar to that of a random walk, but captures the correct large-$r$ limit. We do not have precise theoretical motivation for this ansatz, but as discussed below, it works quite well.

Numerical calculation of the skew for the 1D lattice was performed using the recurrence relation method discussed in Section II.3. The results are shown in Fig. 2 for a few values of $r$. This calculation confirms our initial hypothesis, near neutral fitness the waiting times are uniform enough that the process is well approximated by a biased random walk and the skew approaches 0, scaling in excellent agreement with Eq. (16). When $N\ll N_{c}$, the bias is not sufficient to give the mutants a substantial advantage: the process is dominated by drift and the fixation-time distribution has large skew $\kappa_{3}\approx 1.807$, as found in the preceding section. For $N\gg N_{c}$, selection takes over, the cumulants approach 0, and the distribution becomes normal. A similar crossover appears in the study of the fixation probability, where a transition from $\varphi_{1}\sim 1/N$ to $\varphi_{1}\sim 1-1/r$ is seen when $N$ passes a critical system size (that is slightly different than $N_{c}$). For large fitness $r\gg 1$, the ansatz Eq. (19) captures the scaling behavior if we use an exponent $q=1/2$. Direct numerical simulations of the process confirm that, for any $r>1$, the fixation time on the 1D lattice has an asymptotically normal distribution [Fig. 1(b)].

The random walk approximation allows us to find the asymptotic scaling of the fixation-time cumulants, but ignores the heterogeneity of waiting times present in the Moran process. Using visit statistics we can compute the cumulants exactly with Eq. (8) and rigorously prove they vanish as $N\rightarrow\infty$, verifying that the waiting times have no influence on the asymptotic form of the distribution. For details, see Supplemental Material, Section S3 SM .

Our analysis of the 1D lattice reveals an intriguing universality property of its fixation-time distribution. For any value of relative fitness $r$ other than $r=1$, the fixation-time distribution approaches a normal distribution as $N\rightarrow\infty$. Thus, for $r\neq 1$ the asymptotic shape of the distribution is universal and independent of $r$ (though bear in mind, its mean and variance do depend on $r$).

When $r=1$, corresponding to precisely neutral fitness, the unbiased random walk yields a qualitatively different distribution with considerably larger skew. This qualitative change as $r$ passes through unity leads to a discontinuous jump in the skew at $r=1$.

As one might expect, the discontinuity stems from passage to the infinite-$N$ limit. For finite but large $N$, the distribution varies continuously with $r$, though our numerical results indicate that the sharp increase in skew still occurs very close to $r=1$. We will see in the next section that the discontinuity and highly skewed distribution at neutral fitness persist when we alter the network structure from a locally connected 1D lattice to a fully connected complete graph.

IV Complete Graph

Next we consider the Moran process on a complete graph, useful for modeling well-mixed populations in which all individuals interact. By similar reasoning to above, given $m$ mutants the probability of adding a mutant in the next time step is

[TABLE]

while the probability of subtracting a mutant is

[TABLE]

Interestingly, as we will see in this section, these transition probabilities give rise to a fitness dependent fixation-time distribution, in stark contrast to the universality of the normal distribution observed on the 1D lattice.

IV.1 Neutral Fitness

Again we begin with neutral fitness $r=1$. Now $b_{m}=d_{m}=(Nm-m^{2})/(N^{2}-N)$. The eigenvalues of this transition matrix also have a nice analytical form:

[TABLE]

The asymptotic form of the sums $S_{n}=\sum_{m=1}^{N-1}(\lambda_{m})^{-n}$, can be found by taking the partial fraction decomposition of $(\lambda_{m})^{-n}$ and evaluating each term individually. The resulting cumulants are

[TABLE]

Our knowledge of the eigenvalues also allows us to obtain a series expression for the asymptotic distribution using Eq. (3). For $N\rightarrow\infty$ the standardized distribution is,

[TABLE]

where to leading order in $N$ the mean and standard deviation are $\mu=N^{2}$ and $\sigma=c_{\sigma}N^{2}$, with $c_{\sigma}=\sqrt{\pi^{2}/3-3}$. This distribution was previously found using a different approach by Kimura, who also computed the first few fixation-time moments Kimura (1970). We have extended these results, obtaining the cumulants to all orders. Figure 3 shows that the predicted asymptotic distribution agrees well with numerical experiments.

The numerical value of the fixation-time skew for the Birth-death process on the complete graph is $\kappa_{3}=6\sqrt{3}(10-\pi^{2})/(\pi^{2}-9)^{3/2}\approx 1.6711$, slightly less than that for the 1D lattice. This decrease is the result of two competing effects contributing to the skew. First, since the birth and death transition probabilities are the same, the process is a random walk, which has a highly skewed fixation-time distribution, as shown above. The average time spent in each state, however, varies with $m$. For instance, when $m=1\text{ or }N-1$, $b_{m}\rightarrow 0$ for large $N$. But if $m=\alpha N$ for some constant $0<\alpha<1$ independent of $N$, then $b_{m}$ approaches a constant.

Intuitively, the beginning and end of the mutation-spreading process are very slow because the transition probabilities are exceedingly small. To start, the single mutant must be selected by chance to give birth from the $N$ available nodes, a selection problem like finding a needle in a haystack. Similarly, near fixation the reproducing mutant must find and replace one of the few remaining non-mutants, again choosing it by chance from an enormous population.

The characteristic slowing down at certain states is reminiscent of “coupon collection”, as discussed earlier. Erdős and Rényi proved that for large $N$, the normalized time to complete the coupon collection follows a Gumbel distribution Erdős and Rényi (1961), which we denote by $\text{Gumbel}(\alpha,\beta)$ with density

[TABLE]

For the Moran process, each slow region is produced by long waits for the random selection of rare types of individuals: either mutants near the beginning of the process or non-mutants near the end. In the next section we show that the two coupon collection regions of the Bd process on a complete graph lead to fixation-time distributions that are convolutions of two Gumbel distributions. In the case of neutral fitness, these Gumbel distributions combine with the random walk to produce a new highly skewed distribution with cumulants given by Eq. (IV.1).

IV.2 Non-neutral fitness

We saw in Section III.2 that when the average time spent in each state is constant or slowly varying the fixation-time distribution is asymptotically normal. Birth-death dynamics on the complete graph, however, exhibit coupon collection regions at the beginning and end of the process, where the transition probabilities vanish. We begin this section with a heuristic argument that correctly gives the asymptotic fixation-time distribution in terms of independent iterations of coupon collection.

Differentiating $b_{m}$ with respect to $m$, we find the slope near $m=0$ is $(r+1)/N$, while the slope near $m=N$ has magnitude $(r+1)/(rN)$ for $N\gg 1$. The transition rates approach zero at each of these points, so we expect behavior similar to coupon collection giving rise to two Gumbel distributions. Since the slope is greater for $m$ near 0 than for $m$ near $N$, the Moran process completes its coupon collection faster near the beginning of the process than near fixation.

This heuristic suggests that the asymptotic fixation time should be equal in distribution to the sum of two Gumbel random variables, one weighted by $r$, which is the ratio of the slopes in the coupon collection regions. Specifically, if $T$ is the fixation time with mean $\mu$ and variance $\sigma^{2}$, we expect

[TABLE]

where $\xrightarrow{d}$ means convergence in distribution for large $N$. Here $G$ and $G^{\prime}$ denote independent and identically distributed Gumbel random variables with zero mean and unit variance. It is easy to check that the correct distribution is $\text{Gumbel}(-\gamma\sqrt{6}/\pi,\sqrt{6}/\pi)$, where $\gamma\approx 0.5772$ is the Euler-Mascheroni constant.

Let us make this argument more rigorous. Previous theoretical analysis showed that in the infinite fitness limit, the fixation time has an asymptotically Gumbel distribution Ottino-Löffler et al. (2017). This result can be recovered within our framework, since when $r=\infty$ it follows that $d_{m}=0$, so the eigenvalues of the transition matrix are just $\lambda_{m}=b_{m}=(N-m)/(N-1)$ and the cumulants can be directly calculated using Eq. (4).

For large (but not infinite) fitness, the number of mutants is monotonically increasing, to good approximation, since the probability that the next change in state increases the mutant population is $r/(1+r)\approx 1$. The time spent waiting in each state, however, changes dramatically, especially near $m=1$. Here, $b_{1}\rightarrow 0$ for large $N$, in stark contrast to the infinite fitness system where $b_{1}\rightarrow 1$. The time spent at each state, $t_{m}$ is an exponential random variable, $\mathcal{E}(b_{m}+d_{m})$. In this approximation each state is visited exactly once, so the total fixation time is a sum of these waiting times:

[TABLE]

But this sum of exponential random variables has density given by Eq. (3), with the substitution $\lambda_{m}\rightarrow b_{m}+d_{m}$. Thus, the cumulants of $(T-\mu)/\sigma$ are

[TABLE]

which are exactly the cumulants corresponding to the sum of Gumbel random variables given in Eq. (26). In the limit $r\rightarrow\infty$, the first term in Eq. (28) becomes $1$, and the cumulants are those for a single Gumbel distribution, in agreement with previous results Ottino-Löffler et al. (2017).

Remarkably, these cumulants are exact for any $r>1$, not just in the large-$r$ limit. We can see this directly for the skew $\kappa_{3}$ using the visit statistics approach, computing the asymptotic form of Eq. (8) with the complete graph transition probabilities, Eqs. (20) and (21). Details of the asymptotic analysis are provided in Supplemental Material, Section S4 SM . Numerical simulations of the Moran process corroborate our theoretical results. As shown in Fig. 4, for $r=1.1$ and $r=5$ the agreement between simulated fixation times and the predicted convolution of Gumbel distributions is excellent, at least when $N$ is sufficiently large. Again, our calculations show a discontinuity in the fixation-time distribution at $r=1$. In particular, the $r\rightarrow 1$ limit of the cumulants for non-neutral fitness in Eq. (28) is not the same as the cumulants for neutral fitness found in the preceding section [Eq. (IV.1)].

For smaller networks, it is fascinating to see how the results converge to the asymptotic predictions as $N$ grows. Figure 5 shows how the skew of the fixation-time distribution depends on $r$ and $N$ for the complete graph. As discussed in Section II.4, the fixation-time distributions for these systems are invariant under $r\rightarrow 1/r$. Therefore we show the skew for all $r>0$, to emphasize the intriguing behavior near neutral fitness, where $r=1$. We find that non-uniform convergence of the fixation-time skew leads to the discontinuity predicted at $r=1$. For finite $N$, the skew is a non-monotonic function of $r$ and has a minimum value at some fitness $r_{\text{min}}(N)$. Furthermore, at fixed $r$, the convergence to the $N=\infty$ limit is itself non-monotone. Though beyond the scope of the current study, further investigation of this finite-$N$ behavior would be worth pursuing.

V Partial fixation times

In many applications, we may be interested in the time to partial fixation of the network. For instance, considering cancer progression Williams and Bjerknes (1972); Sottoriva et al. (2013); Bozic et al. (2013) or the incubation of infectious diseases Ottino-Löffler et al. (2017), symptoms can appear in a patient even when a relatively small proportion of cells are malignant or infected. We therefore consider $T_{\alpha}$, the total time to first reach $\alpha N$ mutants on the network, where $0<\alpha<1$. The methods developed in Section II apply to these processes as well. For the eigendecomposition approach we instead use the sub-matrix of $\Omega_{\text{tr}}$ containing the first $\alpha N$ rows and columns. In calculations involving the numerical recurrence relations or visit statistics, we simply cut the sums off at $\alpha N$ instead of $N$ and for the latter, replace $w_{i_{1}i_{2}\cdots i_{n}}(r,N)$ with $w_{i_{1}i_{2}\cdots i_{n}}(r,\alpha N)$.

V.1 One-dimensional lattice

Truncating the Moran Bd process on the 1D lattice by a factor $\alpha$ has no effect on the asymptotic shape of the fixation-time distributions. In both the neutral fitness system and the random walk approximation to the non-neutral fitness system, the transition matrix has no explicit dependence on the state or system size [aside from proportionality factors that cancel in Eq. (4)]. Thus, the eigenvalues are identical to those calculated previously, but correspond to a smaller effective system size $\alpha N$. Taking the limit $N\rightarrow\infty$ therefore yields the same asymptotic distributions found in Section III.

V.2 Complete graph: truncating coupon collection

The complete graph exhibits more interesting dependence on truncation. Since the transition probabilities have state dependence, the eigenvalues change with truncation (they don’t correspond to the same system with smaller effective $N$). Our intuition from coupon collection, however, lets us predict the resulting distribution.

First consider non-neutral fitness. Then there are two coupon collection stages, one near the beginning and another near the end of the process, and together they generate a fixation-time distribution that is a weighted convolution of two Gumbel distributions. The effect of truncating the process near its end should now become clear: it simply removes the second coupon collection. The truncated process stops before the mutants have to laboriously find and replace the last remaining non-mutants. Therefore, we intuitively expect the fixation time for non-neutral fitness to be distributed according to a single Gumbel distribution, regardless of fitness level.

The only exception occurs if $r=\infty$; then no coupon collection occurs at the beginning of the process either, as the lone mutant is guaranteed to be selected to give birth in the first time step, thanks to its infinite fitness advantage. Thus, when fitness is infinite and the process is truncated at the end, both coupon collection phases are removed and the fixation times are normally distributed.

Similar reasoning applies to the Birth-death process with neutral fitness. It also has two coupon collection regions, one of which is removed by truncation. In this case, however, the random walk mechanism contributes to the skew of the overall fixation-time distribution, combining non-trivially with the coupon collection-like process. We find that the skew of the fixation time depends on the truncation factor $\alpha$, varying between $6\sqrt{3}(10-\pi^{2})/(\pi^{2}-9)^{3/2}\approx 1.6711$ when $\alpha=1$, and $\sqrt{3}\approx 1.732$ when $\alpha=0$. A derivation of this $\alpha\rightarrow 0$ limit is given in Supplemental Material, Section S4 SM .

V.3 Summary of main results

The main results from Sections III–V are summarized in Fig. 6, which shows the asymptotic fitness dependence of fixation-time skew for each network considered in this paper. We again show the skew for all $r>0$ (not just $r>1$) to emphasize the discontinuities at zero, neutral, and infinite fitness. On the 1D lattice, independent of the truncation factor $\alpha$, the Bd process has normally distributed fixation times, except at neutral fitness where the distribution is highly skewed. The complete graph fixation-time distributions are the weighted convolution of two Gumbel distributions for $r\neq 1$, again with a highly skewed distribution at $r=1$. With truncation by a factor $\alpha<1$, the distribution for the complete graph is Gumbel for $1<r<\infty$, and normal for $r=\infty$. With neutral fitness the fixation distribution is again highly skewed, with skew dependent on the truncation factor $\alpha$.

VI Extensions

It is natural to ask whether our results are generic; do the same fixation-time distributions appear in other models of evolutionary dynamics? Here we explore the robustness of our results to various changes in the model update dynamics and the network topology. The main finding is that our results are insensitive to these changes, at least qualitatively. The distributions typically remain right-skewed and even follow the same functional forms derived above.

VI.1 Other update dynamics

VI.1.1 Two-fitness Moran process

The Moran Bd processes considered above require a single fitness level, designating the relative reproduction rates between mutants and non-mutants. Another common model is the Moran Birth-Death (BD) process Note (1), which has a second fitness level $\tilde{r}$ measuring the resilience of mutants versus non-mutants during the replacement step Kaveh et al. (2015). Taking this into account, when a mutant or non-mutant is trying to replace its neighbors, mutants are replaced with probability proportional to $1/\tilde{r}$. Taking $\tilde{r}=1$ returns to the model used throughout the preceding sections. The two-fitness model may better capture the complexity of real-world evolutionary systems but does not generally give rise to qualitatively different fixation-time distributions. For brevity, we simply discuss the resulting fixation-time distributions for the BD model. Details supporting the results quoted below are provided in the Supplemental Material, Section S5 SM .

Writing down the transition probabilities for the Moran BD process, we find that $b_{m}/d_{m}\rightarrow r\tilde{r}$ as $N\rightarrow\infty$. This motivates the definition of an effective fitness level, $r_{\text{eff}}=r\tilde{r}$. When $r_{\text{eff}}\neq 1$ our results from above translate to this model. On the 1D lattice the fixation times are normally distributed, while on the complete graph the fixation time distribution is a weighted convolution of Gumbel distributions $G+(r/\tilde{r})G^{\prime}$, with relative weighting $r/\tilde{r}$ (instead of $r$). When $r_{\text{eff}}=1$, the process is asymptotically unbiased and we expect a highly skewed fixation-time distribution. This is indeed the case, although numerical calculations indicate there is an entire family of distributions, dependent on $r=1/\tilde{r}$.

It is interesting to contrast the above observations with a result in evolutionary dynamics known as the isothermal theorem. The theorem states that for $\tilde{r}=1$, the Moran process on a large class of networks, known as isothermal graphs, has fixation probability identical to the complete graph Lieberman et al. (2005). Recent work has shown that this breaks down if $\tilde{r}\neq 1$; the fixation probability develops new network dependence Kaveh et al. (2015). In contrast, even isothermal graphs (including the complete graph and 1D lattice) have fixation-time distributions that depend on network structure. The two-fitness BD model breaks the universality in fixation probabilities predicted by the isothermal theorem, but leads to the same family of fixation distributions that arise due to network structure.

VI.1.2 The Death-Birth Moran process

A two-fitness Death-Birth (DB) Moran process Note (1) is also frequently used to study evolutionary dynamics. In this model, the birth and death events are reversed in order. At each time step a node is chosen at random, with probability proportional to $1/\tilde{r}$, and one of its neighbors is chosen with probability proportional to $r$. The first individual dies and is replaced by an offspring of the same type as the neighbor. The process continues until the mutation either reaches fixation or goes extinct.

The BD and DB processes obey a duality property Kaveh et al. (2015). Starting from the BD transition probabilities, if we swap the two fitness levels $r\leftrightarrow\tilde{r}$ and substitute $m\rightarrow N-m$ (which swaps mutants and non mutants), we obtain the DB transition probabilities. Therefore, the transition matrix for the DB model is identical to that for the corresponding dual BD process, but has the main-, super-, and sub-diagonal entries reversed in order. This leaves the matrix eigenvalues unchanged, so that the DB process has identical fixation-time distributions to those given in the preceding section for the dual BD process.

In principle, the correspondence between DB and BD fixation times could break down for the truncated process considered in Section V. In practice, however, the results are again generally identical. For the truncated DB process, the fixation times on the 1D lattice remain normally distributed. On the complete graph, one of two coupon collection regions is removed by truncation leading to fixation-times following a single Gumbel distribution.

One exception, where the dual models yield different results under truncation, is at infinite fitness. As in Section V, at infinite fitness ($r\rightarrow\infty$) the BD model performs a single coupon collection near fixation, which is cut off by truncation, leading to a normal fixation-time distribution. In contrast, in the dual infinite-fitness DB model ($\tilde{r}\rightarrow\infty$) the coupon collection occurs at the beginning of the process and even under truncation the Gumbel fixation-time distribution is preserved. This effect was previously observed by Ottino-Löffler et al. Ottino-Löffler et al. (2017).

VI.2 Other networks: Approximate results via mean-field transition probabilities

While the 1D lattice and complete graph provide illustrative exactly solvable models of the fitness dependence of fixation-time distributions, other networks may be more realistic. On more complicated networks the analytical tools developed here fail because the transition probabilities (the probability of adding or subtracting a mutant given the current state) depend on the full configuration of mutants, not just the number of mutants. Such systems can still be modeled as a Markov process, but the state space becomes prohibitively large. Fortunately, for certain networks the effect of different configurations can be averaged over, giving a mean-field approximation to the transition probabilities. This approach has been used on a variety of networks to calculate fixation times Hajihashemi and Aghababaei Samani (2019); Ying et al. (2018); Ottino-Löffler et al. (2017, 2017). In this section we discuss how such mean-field approaches can be used to calculate fixation-time distributions for evolution on several different networks.

VI.2.1 Erdős-Rényi random graph

We start with the Erdős-Rényi random graph, for which the mean-field transition probabilities were recently estimated Hajihashemi and Aghababaei Samani (2019). The result is identical to the complete graph probabilities [Eqs. (20)-(21)] up to a constant factor $1-2/Np$, which depends on the edge probability $p$ for the network. This correction is important for computing the mean fixation time, but does not affect the shape of the fixation-time distribution, since proportionality factors cancel in Eq. (4). Therefore we expect the asymptotic fixation-time distribution will be a weighted sum of two Gumbel distributions. This prediction holds for infinite fitness, where the fixation time on an Erdős-Rényi network has a Gumbel distribution Ottino-Löffler et al. (2017).

Preliminary simulations show that the Erdős-Rényi network has the expected fixation-time distributions for $p=1/4$ and $r=2$ (see Figure 7). Further investigation is required to determine the range of fitness and edge probabilities for which this result holds asymptotically (as $N\rightarrow\infty$). For constant $p$, the average degree is proportional to the system size $\langle k\rangle=pN$, similar to the complete graph. It may be, however, that for some $p$ and $r$ the mean-field approximation is not sufficient to capture the higher-order moments determining the shape of the distribution. It is also traditional to consider $N$-dependent edge probabilities with $p(N)$ chosen, for example, to fix $\langle k\rangle$. It is unclear whether such graphs will behave like the ring (due to their sparsity), like the complete graph (due to their short average path length), or somewhere in between these extreme cases. In the same vein, which other networks admit accurate mean-field approximations to the transition probabilities? Do many complex networks have fixation-time distributions identical to the complete graph?

VI.2.2 Stars and superstars: evolutionary amplifiers

Another nice approximation maps the Moran process on a star graph, a simple amplifier of selection, onto a birth-death Markov chain Frean et al. (2013). The resulting transition probabilities exhibit coupon collection regions, similar to the complete graph. The ratio of slopes near these regions (few mutants or non-mutants), however, is $r^{2}$. Our heuristic predicts the fixation-time distribution on the star is $G+r^{2}G^{\prime}$. In addition to amplifying fixation probability, the star increases fixation-time skew. This raises a broader question: do evolutionary amplifiers also amplify fixation-time skew? Computing fixation times for evolutionary dynamics on superstars (which more strongly amplify selection Lieberman et al. (2005)) remains an open problem.

VI.2.3 Growth of cancerous tumors: evolutionary dynamics on $d$-dimensional lattices

Mean-field arguments have also been applied to $d$-dimensional lattices in the infinite-fitness limit Ottino-Löffler et al. (2017, 2017). In this limit the mutant population grows in an approximately spherical shape near the beginning of the process and the population of non-mutants is approximately spherical near fixation. The surface area to volume ratio of the $d$-dimensional sphere gives the probability of adding a mutant. With finite fitness, non-mutants can now replace their counterparts and the surface of the sphere of growing mutants roughens Williams and Bjerknes (1972). For near-neutral fitness, the configuration of mutants resembles the shape of real cancerous tumors. Perhaps mean-field approaches can draw connections between the fitness-dependent roughness of growing mutant populations and fixation-time distributions for evolution on lattices.

VII Summary

In this paper we have obtained the first closed-form solutions for the fitness dependence of fixation-time distributions of the Moran Birth-death process on the 1D lattice and complete graph. Previous analyses were restricted to the limit of infinite fitness, with some partial results for neutral fitness. To reiterate our new results: There is a dichotomy between neutral and non-neutral fitness. When fitness is neutral, the distribution always exhibits a discontinuity; whether the graph is complete or a 1D lattice, the skew jumps up discontinuously in either case. On the other hand, when fitness is non-neutral but otherwise arbitrary, the results depend strongly on network topology. Specifically, on the complete graph the fixation-time distribution is a fitness-weighted convolution of Gumbel distributions and hence is always skewed, whereas on the 1D lattice the distribution is normal and hence is never skewed.

Together with the mean and variance, the distributions derived here give a complete statistical description of the asymptotic fixation time (see Table 1). Our analysis revealed that these results are robust in the sense that similar distributions arise under truncation, in some other models, and in some other network structures, including the Erdős-Rényi random graph.

VIII Future Directions

Though the model we have focused on here (the Moran Birth-death model) is deliberately simplified, we expect our results will be useful in applications. For instance, the theory should allow a more refined analysis of the rate of evolution, by extending the seminal work by Kimura, whose neutral theory of evolution predicted a molecular clock Kimura (1968). In his model, neutral mutations become fixed at a constant rate, independent of population size. This result, with some refinements, is now used widely in estimating evolutionary time scales Kumar (2005). The fixation-time distributions discussed here should allow one to go beyond Kimura’s classic analysis to capture the full range of evolutionary outcomes, by providing information about the expected deviations from the constant-rate molecular clock, as well as how this prediction is affected by population structure. More generally, it would be interesting to study the implications of these distributions for rates of evolution at various fitness levels.

Furthermore, our results provide concrete predictions that are testable via bacterial evolution experiments. Does the same fitness and network structure dependence of fixation-time distributions arise in real systems?

Future theoretical studies could analyze random networks and lattices more deeply, as well as stars and superstars, the prototypical evolutionary amplifiers Lieberman et al. (2005). More sophisticated models involving evolutionary games are also of interest. These have skewed fixation-time distributions Ashcroft et al. (2015) whose asymptotic form remains unknown. Finally, we hope that methods developed here will prove useful in other areas, such as epidemiology Doering et al. (2007), ecology Doering et al. (2005), and protein folding Zwanzig et al. (1992), where stochastic dynamics may similarly give rise to skewed first-passage times.

Acknowledgments

We thank Bertrand Ottino-Löffler and Jacob Scott for helpful discussions. Research supported by the Cornell Presidential Life Science Fellowship and NSF Graduate Research fellowship grant DGE-1650441 to D. H. and by NSF grants DMS-1513179 and CCF-1522054 to S.H.S.

Appendix A Visit statistics

In this Appendix we formulate the visit statistics approach. We first provide further details in the derivation of the series expression for the fixation-time cumulants given in Eq. (8), and then explicitly compute the weighting factors that appear in this expression to third order. This result requires constant selection, $b_{m}/d_{m}=r$, as is the case for the Moran process. Under constant selection the transient transition matrix can be written as $\Omega_{\text{tr}}=\Omega_{\text{RW}}D$, where $D$ is diagonal with elements $D_{mm}=b_{m}+d_{m}$ and $\Omega_{\text{RW}}$ is the transition matrix for a random walk,

[TABLE]

Since we are interested in the fixation-time distribution, we condition on fixation occurring. As discussed in Section II.2 (see also Supplemental Material, Section I SM ), the conditioned transition matrix $\tilde{\Omega}_{\text{tr}}=S\,\Omega_{\text{tr}}\,S^{-1}$, where $S$ is diagonal with $S_{mm}=1-1/r^{m}$. Combining these results, we have that

[TABLE]

where we have used the fact that both $D$ and $S$ are diagonal matrices, and therefore commute.

We found in Section II.2 that the moments of the fixation time $T$ can be expressed as,

[TABLE]

where $\mathbf{1}$ is a row vector of ones and $\mathbf{p}_{\text{tr}}(0)$ is the initial state of the system, with $[p_{\text{tr}}(0)]_{m}=\delta_{m,m_{0}}$ for $m_{0}$ initial mutants. To compute these moments, we need the inverse $\tilde{\Omega}_{\text{tr}}^{-1}=D^{-1}S\,\Omega_{RW}^{-1}S^{-1}$. Since $\Omega_{RW}$ is a tridiagonal Toeplitz matrix, its inverse has a well-known form da Fonseca and Petronilho (2001):

[TABLE]

Hence the matrix $V=-S\,\Omega_{\text{RW}}^{-1}S^{-1}$ has elements

[TABLE]

The matrix $V$, sometimes called the fundamental matrix, encodes the visit statistics of the conditioned random walk: $V_{ij}$ is the mean number of visits to state $i$ from state $j$ before hitting the absorbing state $N$ Kemeny and Snell (1983). The Moran process has the same visit statistics, but on average spends a different amount of time, designated by $(b_{i}+d_{i})^{-1}$, waiting in each state.

While one could now compute the moments $\mu_{n}$ in Eq. (31) directly, we find that the cumulants yield nicer expressions. Furthermore, the normal and Gumbel fixation-time distributions, predicted by our simulations and approximate calculations, are more simply described in terms of their cumulants. The non-standardized cumulants $\kappa^{\prime}_{n}$ are linear combinations involving products of moments whose orders sum to $n$. Thus each term in the cumulants has $n$ powers of $D$ producing $n$ factors of $(b_{i}+d_{i})^{-1}$ with a weight designated by the visit statistics. With this observation, it is clear the standardized cumulants $\kappa_{n}=\kappa^{\prime}_{n}/(\kappa^{\prime}_{2})^{n/2}$ have the form given in Eq. (8),

[TABLE]

where $w_{i_{1}i_{2}\cdots i_{n}}^{n}(r,N|m_{0})$ are the weighting factors coming entirely from the visit statistics of a biased random walk (starting from $m_{0}$ initial mutants). As in the main text, we take the initial state to be a single mutant $m_{0}=1$, and will suppress the dependence of the weighting factors on initial condition, writing $w_{i_{1}i_{2}\cdots i_{n}}^{n}(r,N)$ instead. Generalizations to other cases are straightforward and are discussed briefly below.

We emphasize that even without explicit knowledge of the factors $w_{i_{1}i_{2}\cdots i_{n}}^{n}(r,N)$, this formulation can be extremely useful. For instance when $b_{i}+d_{i}$ is constant, these are just the cumulants for the (possibly biased) random walk, which were computed in Section III to approximate the Moran process on the 1D lattice. In particular, the sums over weighting factors obtained from setting $b_{i}+d_{i}=1$ in Eq. (34) have leading asymptotic form given by Eq. (14). This fact can be used to bound the cumulants even when $b_{i}+d_{i}\neq 1$, which in some cases is sufficient to determine the leading asymptotic behavior. When this is not possible, the weighting factors must be computed explicitly. We now turn our focus to derviing $w_{ij}^{2}(r,N)$ and $w_{ijk}^{3}(r,N)$.

We can compute the weighting factors by writing out the matrix multiplication of $\tilde{\Omega}_{\text{tr}}^{-1}$. First note that

[TABLE]

Then the first three moments of the fixation time are,

[TABLE]

The corresponding non-standardized cumulants are given by the usual formulas, $\kappa^{\prime}_{2}=\mu_{2}-\mu_{1}^{2}$ and $\kappa^{\prime}_{3}=\mu_{3}-3\mu_{2}\mu_{1}+2\mu_{1}^{3}$. In terms of the visit numbers the non-standardized cumulants become

[TABLE]

From here we can read off the weighting factors accordingly. For convenience, we can choose $w_{ij}^{2}(r,N)$ and $w_{ijk}^{3}(r,N)$ to be symmetric by averaging the numerators in Eq. (37) over the permutations of the indices. Then,

[TABLE]

where $\Pi_{2}$ is the set of permutations of $\{i,j\}$ and $\Pi_{3}$ are the permutations of $\{i,j,k\}$. We note that these expressions also hold for general initial condition by replacing the subscript $1$ with $m_{0}$. Plugging Eq. (33) into this expression for $w_{ij}^{2}$ we obtain, after some algebra,

[TABLE]

for $i\geq j$. Since we have constructed $w_{ij}^{2}(r,N)$ to be symmetric, when $j>i$ the formula is identical with $i$ and $j$ exchanged. Similarly, using Eq. (33) together with the expression for $w_{ijk}^{3}$ in Eq. (38) leads to

[TABLE]

for $i\geq j\geq k$. Again, the formula for different orderings of the indices $i,j,k$ is the same with the indices permuted appropriately, so that $w_{ijk}^{3}$ is perfectly symmetric.

This completes the derivation of the visit statistics expression for the fixation-time cumulants. Together, Eqs. (34), (39) and (40) give a closed form expression for the fixation-time skew which is manageable for the purpose of asymptotic approximations. The diagonal terms in the higher-order weighting factors are also particularly simple, $w_{ii\cdots i}^{n}(r,N)=(n-1)!V_{ii}^{n}$. While we will not explicitly compute them, the off diagonal weights $w_{i_{1}i_{2}\cdots i_{n}}^{n}(r,N)$ can be found by a straightforward generalization of the above procedure. Example applications of this approach are given in Supplemental Material, Sections III and IV SM , where we show that all cumulants of the fixation time vanish for the Moran process on the 1D lattice and compute the asymptotic skew for the Moran process on the complete graph.