An Adjacent-Swap Markov Chain on Coalescent Trees

The standard coalescent (Kingman:1982uj, ) is often used in evolutionary biology and population genetics to model the set of ancestral relationships among $n$ individuals as a rooted and timed binary tree called a genealogy (wakeley_coalescent_2008, ). In this manuscript, we are concerned with discrete tree topologies only and assume a unit interval length between consecutive coalescing (branching) events.

A labeled ranked tree shape $T^{L}\in\mathcal{T}^{L}_{n}$ with $n$ leaves is a rooted labeled and ranked binary tree. The leaves are labeled by the set $\{\ell_{1},\ell_{2},\ldots,\ell_{n}\}$ and internal nodes are labeled $1,\ldots,n-1$ in increasing order, with label $n-1$ at the root (Figure 1(A)). The cardinality of $\mathcal{T}^{L}_{n}$ is given by

$$|\mathcal{T}^{L}_{n}|=\frac{n!(n-1)!}{2^{n-1}}.$$

(1)

In many applications, the quantity of interest is the overall shape of the tree (Kirkpatrick1993, ; Frost2013, ; Maliet2018, ); the specific labels of the samples are un-important. In fact, a key feature of simple coalescent models, such as the Kingman coalescent, is the assumption that the leaves are exchangeable. The Tajima coalescent, a lumping of the Kingman coalescent that ignores leaf labels, was recently proposed for inferring past population size trajectories from binary molecular data (Palacios2019, ). This motivates the study of unlabeled ranked tree shapes with $n$ leaves $\mathcal{T}_{n}$ (Figure 1(B)). The appealing feature of modeling unlabeled ranked tree shapes is the reduction in the cardinality of the state space. The cardinality of $\mathcal{T}_{n}$ is given by the ($n-1$)-th Euler zig-zag number $E_{n-1}\ll|T^{L}_{n}|$ (kuznetsov1994increasing, ).

Other types of trees are used in the context of evolutionary models, such as cladeograms (or phylogenies). Cladeograms are rooted or unrooted binary trees with labeled leaves and no ranking of internal nodes. In diaconis1998matchings , the space of cladeograms with $l$ leaves is shown to be in bijection with the space of perfect matchings of $2n$ labels, with $n=l-1$. The bijection with perfect matchings is used to define a Markov chain on the space of cladeograms. Using the perfect matching perspective, the chain is analogous to a random transpositions chain on the space $S_{2n}$ of permutations of $2n$ elements. This representation allows the use of known results for the random transpositions chain to be applied to the chain on trees to determine sharp bounds on rates of convergence diaconis2002random . Motivated by this result, we propose a new representation of labeled and unlabeled ranked tree shapes as a type of matching. We use this representation to define an adjacent-swap chain on the spaces of labeled and unlabeled ranked tree shapes, and to study their mixing time. It is important to study the mixing time for chains on ranked trees because these results have implications for the design of MCMC algorithms to infer evolutionary parameters from molecular data in a Bayesian setting (drummond2012bayesian, ).

In the rest of this section we define a bijective representation of ranked trees as constrained ordered matchings and define the adjacent-swap chain. We then state our main results concerning mixing times and discuss related work. In Section 2, we state known results on Markov chains that will be used in the following sections. In section 3 we state important properties of the proposed adjacent-swap chain. The lower bound for the mixing time of the adjacent-swap chains is proven in section 4 and it is obtained by finding a specific function that gives a lower bound for the relaxation time of the chains.In section 5, we prove the upper bounds on the mixing time of the adjacent-swap chains via a coupling argument. In the discussion section we mention future directions and other related chains on labeled and un-labeled coalescent trees which could be defined using the matching representations, as well as discuss implications for MCMC methods applied in practice.

In this section, we review some results of Markov chains theory that will be used to bound the mixing time of the adjacent swap chains on $\mathcal{T}^{L}_{n}$ and $\mathcal{T}_{n}$. We use the coupling method to find the upper bounds. A coupling of Markov chains with transition matrix $P$ is a process $(X_{t},Y_{t})_{t\geq 0}$ such that both $(X_{t})$ and $(Y_{t})$ are marginally Markov chains with transition matrix $P$. The goal is to define a coupling of two chains from different starting distributions such that the two chains will quickly reach the same state, and once this occurs the two chains stay matched forever. This coupling time gives an upper bound on the mixing time, see Chapter 5 from LevinPeresWilmer2006 for more details.

To find a lower bound on the mixing time, we bound the relaxation time and use the following result:

The relaxation time is defined as the inverse of the spectral gap: $\tau_{n}=1/\gamma$, where $\gamma=1-|\lambda_{n,2}|$, one minus the absolute value of the second-largest eigenvalue of the transition matrix of the chain. We will bound the relaxation time using the following variational characterization:

We note that the transition matrix $P$ of the adjacent-swap chain on unlabeled ranked tree shapes is not symmetric. For example, consider the transition $x$ to $y$ of the type:

$$(0,a)^{k},(0,b)^{k+1}\to(0,0)^{k},(a,b)^{k+1},$$

(4)

for labels $a,b\neq 0$ (recall the $0$s represent leaf labels). The probability of this transition is $P(x,y)=1/(n-2)\cdot 1/4$, since once pairs $k$ and $k+1$ are chosen, there is a $1/4$ chance of choosing label $a$ from pair $k$ and label $0$ from pair $k+1$. However, the reverse transition:

$$(0,0)^{k},(a,b)^{k+1}\to(0,a)^{k},(0,b)^{k+1}$$

has probability $1/(n-2)\cdot 1/2$. It turns out that the stationary distribution $\pi$ for the chain is the Tajima coalescent distribution veber (also known as Yule distribution). For $T\in\mathcal{T}_{n}$

$$\pi(T)=\frac{2^{n-c(T)-1}}{(n-1)!},$$

where $c(T)$ is the number of cherries of $T$, i.e. the number of pairs of the type $(0,0)$ in the $COM_{n-1}(S)$ representation. Indeed, for transitions $x$ to $y$ of the type of (4), $y$ has one more cherry than $x$ and $\pi(x)P(x,y)=\pi(y)P(y,x)$; for other types of transitions $P(x,y)>0$ that do not affect the number of cherries, $P(x,y)=P(y,x)$, $\pi(x)=\pi(y)$ and the detailed balance equation is satisfied.

Another way of proving the Tajima coalescent distribution is the stationary distribution on unlabeled ranked tree shapes is to view the chain as a lumping of the chain on $\mathcal{T}^{L}_{n}$. This perspective also gives us an initial comparison of the relaxation times of the chains. The space $\mathcal{T}_{n}$ of unlabeled ranked-tree shapes can be considered as a set of equivalency classes of the trees $\mathcal{T}_{n}^{L}$. That is, for trees $x,y\in\mathcal{T}_{n}^{L}$, define the equivalence relation $x\sim y$ if the trees have the same ranked tree shape. From the $COM_{n-1}(S)$ perspective with $S=\{\ell_{1},\dots,\ell_{n},1,2,\dots,n-2\}$, two matchings are equivalent if all internal node labels $1,\ldots,n-2$ occur in the same pairs.

This equivalence relation induces a partition of $\mathcal{T}_{n}^{L}$ using equivalence classes. That is, we can write $\mathcal{T}_{n}^{L}$ as the disjoint union of sets $\Omega_{1},\dots,\Omega_{M}$, where $M=|\mathcal{T}_{n}|$, and all trees in $\Omega_{i}$ have the same ranked-tree shape.

In words, Lemma 3 says that probability of transitioning to a different ranked-tree shape is independent of the leaf configuration of the current tree. A well-known result (e.g. Lemma 2.5 in LevinPeresWilmer2006 ) is that the induced chain on the space of equivalence classes defined by $\tilde{P}([x],[y]):=P(x,[y])$ is a Markov chain; then observe that the induced chain $\tilde{P}$ is equivalent to adjacent-swap Markov chain $P$ on the space $\mathcal{T}_{n}$. This allows a comparison of the spectral gaps of the two chains. Lemma 12.9 in (LevinPeresWilmer2006, ) states the eigenvalues of the transition matrix of the lumped chain are eigenvalues of the transition matrix of the full chain, hence we have the following lemma.

From Lemma 3 we can immediately see that stationary distribution for $P$ is indeed the Tajima coalescent distribution. Suppose $T\in\mathcal{T}_{n}^{U}$ corresponds to the equivalence class $\Omega_{i}\subset\mathcal{T}_{n}^{L}$. Then,

$$\pi(T)=\sum_{y\in\Omega_{i}}\pi^{L}(y)=\frac{|\Omega_{i}|}{|\mathcal{T}_{n}^{L}|}.$$

A uniformly sampled labeled ranked tree shape with the leaf-labels erased gives an unlabeled ranked tree shape according to the Tajima distribution (veber, ). This implies that $|\Omega_{i}|/|\mathcal{T}_{n}^{L}|$ is equal to the Tajima distribution for the ranked tree shape corresponding to $\Omega_{i}$.

To prove the lower bounds in Theorem 1, we use the variational characterization of the relaxation time:

$$\tau_{n}=\sup_{f:\Omega\to\mathbb{R}}\frac{Var_{\pi}(f)}{\mathcal{E}(f)}.$$

(6)

We will use the internal tree length as a function $\varphi(x):\mathcal{T}_{n}\rightarrow\mathbb{R}$ to find a lower bound for $\tau_{n}$. Since the internal tree length is independent of the leaf labels, we get the same lower mixing time bound for the adjacent-swap chains on $\mathcal{T}_{n}$ and $\mathcal{T}^{L}_{n}$. For ease of exposition we will focus the following discussion on $\mathcal{T}_{n}$.

We will define the internal tree length function on the space $\mathcal{T}_{n}$ using the correspondence with $COM_{n-1}(S)$ and assuming unit length between consecutive coalescence events. For a label $j\in\{1,\dots,n-2\}$, let $I_{x}(j)$ denote the pair index $k$ such that $j\in\textbf{p}_{k}$. In the example of Figure 3, $I_{x}(2)=5$. Note that from the constraints (Eq. (2)), we have that

$$I_{x}(j)\in\{j+1,\dots,n-1\}.$$

We now define the internal tree length function on $\mathcal{T}_{n}$ as follows:

$$\varphi(x)=\sum_{j=1}^{n-2}(I_{x}(j)-j).$$

(7)

As an example consider the ranked unlabeled tree shape with $n=7$ depicted in Figure 3. Note that the function is minimized for the “caterpillar tree”: $(0,0)^{1},(0,1)^{2},(0,2)^{3},\ldots,(0,n-2)^{n-1}$, where $I_{x}(j)=j+1$ for all $j$, giving $\varphi(x)=n-2$.

This function $\varphi$ is useful in bounding the relaxation time because it is “local” with respect to the Markov chain. That is, suppose $P(x,y)>0$ and $x\neq y$. The change from $x$ to $y$ could have moved an interior node label to the left, in which case $\varphi(y)=\varphi(x)-1$. If it moved an interior node label to the right, then $\varphi(y)=\varphi(x)+1$. If two interior node labels were swapped, then $\varphi(y)=\varphi(x)$. Thus, $(\varphi(x)-\varphi(y))^{2}\leq 1$. The denominator of Equation (6) can then be upper-bounded by $1/2$.

To find the variance of the internal tree length, we note that Var${}_{\pi}[\varphi(X)]=$Var${}_{\pi^{L}}[\varphi(X)]$ since the internal tree length is independent of the leaf labels. We now re-state the standard coalescent of labeled ranked tree shapes as an urn process.

Consider an urn containing $n$ balls labeled $\ell_{1},\ldots,\ell_{n}$. At step k, draw two balls without replacement from the urn. Let $l_{1},l_{2}$ be the labels of the balls drawn, then set $\mathbf{p}_{k}=(l_{1},l_{2})^{k}$. Add in a new ball with label $k$. Repeat this process for $k=1,\dots,n-1$ until only a single ball labeled $(n-1)$ remains in the urn. The resulting sequence of pairs $\mathbf{p}=(\mathbf{p}_{1},\dots,\mathbf{p}_{n-1})$ corresponds to a labeled ranked tree shape $T\in\mathcal{T}^{L}_{n}$ drawn from the standard coalescent, i.e. $\mathbf{p}\sim\pi^{L}$.

We can now simplify the urn process as follows: start with $n$ white balls in the urn. At each step, draw two balls and add back in a single red ball (representing an interior node of the tree). Let $R_{0}=0$ and $R_{k}$ be the number of red balls in the urn after $k$ coalescence events. Note that after $k$ mergers (coalescence events), there are $n-k$ total balls left and $R_{n-1}=1$.

The simplified urn process is useful because the internal tree length can be computed by counting the number of red balls at each step and adding them all together (Figure 3), i.e.,

$$\varphi(\mathbf{p})=\sum_{k=1}^{n-2}R_{k}.$$

The quantities $\{R_{k}\}_{k=1}^{n-2}$ are fairly easy to analyze and have been studied before in various contexts. In janson2011 , the values are used to study the asymptotic behavior of the external tree length of the Kingman coalescent and it is shown that:

	$\displaystyle\mathbb{E}[R_{k}]=\frac{k(n-k)}{n-1},$
	$\displaystyle Cov(R_{k}\cdot R_{l})=\frac{k(k-1)(n-l)(n-l-1)}{(n-1)^{2}(n-2)},\,\,\,\,k\leq l.$

Using this, we get

$$\mathbb{E}_{\pi}[\varphi(\mathbf{p})]=\sum_{k=1}^{n-2}\frac{k(n-k)}{n-1}=\frac{1}{6}(n^{2}+n-6),$$

Then, we calculate

	$\displaystyle\sum_{k=1}^{n-2}\mathbb{E}[R_{k}^{2}]=\frac{1}{30}(n^{3}+3n^{2}+2n-30),$		(8)
	$\displaystyle 2\cdot\sum_{k=1}^{n-3}\sum_{l=k+1}^{n-2}\mathbb{E}[R_{k}R_{l}]=\frac{1}{180}(n-3)(5n^{3}+21n^{2}-14n-120).$		(9)

The second moment $\mathbb{E}_{\pi}[\varphi(\mathbf{p})^{2}]$ is the sum of these two lines (8) + (9), which gives

$$\text{Var}_{\pi}(\varphi(\mathbf{p}))=\mathbb{E}_{\pi}[\varphi(\mathbf{p})^{2}]-\mathbb{E}_{\pi}[\varphi(\mathbf{p})]^{2}=\frac{1}{90}n(n+1)(n-3).$$

The lower bound on the relaxation time is applicable to labeled and unlabeled ranked tree shapes, however Lemma 4 further allows us to state the first inequality on the left hand side of Theorem 5.

In this section we prove the upper bounds in Theorem 1 using a coupling argument. The coupling is similar to the one used by Aldous for analyzing the adjacent transpositions chain on $S_{n}$ aldous1983random .

We analyze a lazy version of the chain to make the coupling. That is, at each step, we generate a random coin flip $\theta\sim\text{Bernoulli}(1/2)$. If $\theta=1$, attempt a move of the chain, and if $\theta=0$ then make no change. Let $X_{t}=(\textbf{p}_{1},\dots,\textbf{p}_{n-1})$ and $Y_{t}=(\textbf{q}_{1},\dots,\textbf{q}_{n-1})$ be the two copies of the chain at time $t$, one started from $x$ and the other from $y$. We will first describe the coupling for the chain on $\mathcal{T}_{n}$.