Species Trees Forcing Parsimony to Fail

TL;DR

This paper demonstrates that the parsimony method for inferring species trees can fail under specific conditions involving three consecutive inner edges with certain lengths, especially as species trees grow larger.

Contribution

It proves that no matter the true 9-taxon or larger species tree topology, parsimony fails if it contains three specific consecutive inner edges with particular length proportions.

Findings

Parsimony fails with certain edge length configurations.

Failure probability increases with species tree size.

Applicable only when such edge length conditions are not met.

Abstract

To the known fact that Parsimony method sometimes fails on the problem of inferring species trees from gene trees, here we proved that no mater of what topology the true 9-taxon and greater species tree is the only thing one needs to break down Parsimony is to have in this species tree three consecutive inner edges not going through the root but perhaps ending on it with lengths T1, T2, T3 of some proportions. Obviously, the probability to meet these lengths is growing in general with the size of species tree. Therefore, Parsimony may be applied only when the described lengths of edges cannot be met in the tree.

Tables

Table 1. Table 1: 15 topological 5-taxon trees.

${\tau }_1$: $((B,C),A,(D,E))$,   ${\tau }_2$: $((C,D),A,(B,E))$,   ${\tau }_3$: $((C,E),A,(B,D))$,

${\tau }_4$: $((A,E),B,(C,D))$,   ${\tau }_5$: $((A,D),B,(C,E))$,   ${\tau }_6$: $((A,C),B,(D,E))$,

${\tau }_7$: $((A,B),C,(D,E))$,   ${\tau }_8$: $((A,D),C,(B,E))$,   ${\tau }_9$: $((A,E),C,(B,D))$,

${\tau }_{10}$: $((A,B),D,(C,E))$,   ${\tau }_{11}$: $((A,C),D,(B,E))$,   ${\tau }_{12}$: $((A,E),D,(B,C))$,

${\tau }_{13}$: $((A,B),E,(C,D))$,   ${\tau }_{14}$: $((A,C),E,(B,D))$,   ${\tau }_{15}$: $((A,D),E,(B,C))$.

Table 2. Table 2: The parsimony scores p ​ a ​ r ​ s τ j ​ ( σ i ) ≡ m i ​ j 𝑝 𝑎 𝑟 subscript 𝑠 subscript 𝜏 𝑗 subscript 𝜎 𝑖 subscript 𝑚 𝑖 𝑗 pars_{\tau_{j}}(\sigma_{i})\equiv m_{ij} for all 15 possible output trees σ 𝜎 \sigma with respect to the matrix representation of all 15 possible input trees τ 𝜏 \tau .

	${\tau }_1$	${\tau }_2$	${\tau }_3$	${\tau }_4$	${\tau }_5$	${\tau }_6$	${\tau }_7$	${\tau }_8$	${\tau }_9$	${\tau }_{10}$	${\tau }_{11}$	${\tau }_{12}$	${\tau }_{13}$	${\tau }_{14}$	${\tau }_{15}$
${\sigma }_1$	2	4	4	4	4	3	3	4	4	4	4	3	4	4	3
${\sigma }_2$	4	2	4	3	4	4	4	3	4	4	3	4	3	4	4
${\sigma }_3$	4	4	2	4	3	4	4	4	3	3	4	4	4	3	4
${\sigma }_4$	4	3	4	2	4	4	4	4	3	4	4	3	3	4	4
${\sigma }_5$	4	4	3	4	2	4	4	3	4	3	4	4	4	4	3
${\sigma }_6$	3	4	4	4	4	2	3	4	4	4	3	4	4	3	4
${\sigma }_7$	3	4	4	4	4	3	2	4	4	3	4	4	3	4	4
${\sigma }_8$	4	3	4	4	3	4	4	2	4	4	3	4	4	4	3
${\sigma }_9$	4	4	3	3	4	4	4	4	2	4	4	3	4	3	4
${\sigma }_{10}$	4	4	3	4	3	4	3	4	4	2	4	4	3	4	4
${\sigma }_{11}$	4	3	4	4	4	3	4	3	4	4	2	4	4	3	4
${\sigma }_{12}$	3	4	4	3	4	4	4	4	3	4	4	2	4	4	3
${\sigma }_{13}$	4	3	4	3	4	4	3	4	4	3	4	4	2	4	4
${\sigma }_{14}$	4	4	3	4	4	3	4	4	3	4	3	4	4	2	4
${\sigma }_{15}$	3	4	4	4	3	4	4	3	4	4	4	3	4	4	2

Table 3. Table 3: The probabilities   p i = p ​ ( τ i | Σ 1 ) subscript 𝑝 𝑖 𝑝 conditional subscript 𝜏 𝑖 subscript Σ 1 p_{i}=p(\tau_{i}|\Sigma_{1})  for   i = 1 , 15 ¯ 𝑖 ¯ 1 15 i=\overline{1,15} ,  where   X = e − T 1 𝑋 superscript 𝑒 subscript 𝑇 1 X=e^{-T_{1}} , Y = e − T 2 , Z = e − T 3 formulae-sequence 𝑌 superscript 𝑒 subscript 𝑇 2 𝑍 superscript 𝑒 subscript 𝑇 3 Y=e^{-T_{2}},\ Z=e^{-T_{3}} .

$p_1$	$X/3-XY/3+X{Y}^3/18+X{Y}^3{Z}^6/90$
$p_2$	$X{Y}^3/18+X{Y}^3{Z}^6/90$
$p_3$	$X{Y}^3/18+X{Y}^3{Z}^6/90$
$p_4$	$X{Y}^3/18+X{Y}^3{Z}^6/90$
$p_5$	$X{Y}^3/18+X{Y}^3{Z}^6/90$
$p_6$	$X/3-XY/3+X{Y}^3/18+X{Y}^3{Z}^6/90$
$p_7$	$1-2X/3-2Y/3+XY/3+X{Y}^3/18+X{Y}^3{Z}^6/90$
$p_8$	$X{Y}^3/18+X{Y}^3{Z}^6/90$
$p_9$	$X{Y}^3/18+X{Y}^3{Z}^6/90$
$p_{10}$	$Y/3-XY/6-X{Y}^3/9+X{Y}^3{Z}^6/90$
$p_{11}$	$XY/6-X{Y}^3/9+X{Y}^3{Z}^6/90$
$p_{12}$	$XY/6-X{Y}^3/9+X{Y}^3{Z}^6/90$
$p_{13}$	$Y/3-XY/6-X{Y}^3/18-2X{Y}^3{Z}^6/45$
$p_{14}$	$XY/6-X{Y}^3/18-2X{Y}^3{Z}^6/45$
$p_{15}$	$XY/6-X{Y}^3/18-2X{Y}^3{Z}^6/45$

Table 4. Table 4: The probabilities p i p ​ c = p ​ ( τ i | Σ 2 ) subscript superscript 𝑝 𝑝 𝑐 𝑖 𝑝 conditional subscript 𝜏 𝑖 subscript Σ 2 p^{pc}_{i}=p(\tau_{i}|\Sigma_{2}) for i = 1 , 15 ¯ 𝑖 ¯ 1 15 i=\overline{1,15} .

$p_1^{pc}$	$=p_3^{pc}=p_5^{pc}=p_6^{pc}=p_8^{pc}=p_9^{pc}=p_{11}^{pc}=p_{12}^{pc}=\left(1/18\right)XY+XY{Z}^6/90,$
$p_2^{pc}$	$=p_4^{pc}=p_7^{pc}=p_{10}^{pc}=XY{Z}^6/90+Y/3-\left(5/18\right)XY$,
$p_{13}^{pc}$	$=1-\left(2/45\right)XY{Z}^6+\left(4/9\right)XY-\left(2/3\right)X-\left(2/3\right)Y,$
$p_{14}^{pc}$	$=p_{15}^{pc}=-\left(2/45\right)XY{Z}^6+XY/9.$

Table 5. Table 5: The probabilities p i p ​ b = p ​ ( τ i | σ 7 ) subscript superscript 𝑝 𝑝 𝑏 𝑖 𝑝 conditional subscript 𝜏 𝑖 subscript 𝜎 7 p^{pb}_{i}=p(\tau_{i}|\sigma_{7}) for i = 1 , 15 ¯ 𝑖 ¯ 1 15 i=\overline{1,15} .

$p_1^{pb}$	$=p_6^{pb}=-\left(1/3\right)ZXY+\left(1/15\right)X{Y}^{3}Z+\left(1/3\right)X,$
$p_2^{pb}$	$=p_3^{pb}=p_4^{pb}=p_5^{pb}=\left(1/15\right)X{Y}^{3}Z,$
$p_7^{pb}$	$=1+\left(1/3\right)ZXY-\left(2/3\right)X+\left(1/15\right)X{Y}^{3}Z-\left(2/3\right)ZY,$
$p_8^{pb}$	$=p_9^{pb}=\left(1/15\right)X{Y}^{3}Z,$
$p_{10}^{pb}$	$=p_{13}^{pb}=\left(1/3\right)ZY-\left(1/6\right)ZXY-\left(1/10\right)X{Y}^{3}Z,$
$p_{11}^{pb}$	$=p_{12}^{pb}=p_{14}^{pb}=p_{15}^{pb}=\left(1/6\right)ZXY-\left(1/10\right)X{Y}^{3}Z.$

Equations

$$g_{ij}(T)=\sum _{k=j}^i{e}^{-k(k-1)T/2}\frac{(2k-1)(-1{)}^{k-j}{j}_{(k-1)}{i}_{[k]}}{j!(k-j)!i_{(k)}},$$

$${\mathbf{s}}^{cat}:=M\mathbf{p},$$

$$\mathbf{s}^{cat}:=M\mathbf{p}=\left(\begin{array}[]{l}3-X/3+2Y/3+XY/3-XY^{3}/18-XY^{3}Z^{6}/90\\ 4-Y/3-XY^{3}/18-XY^{3}Z^{6}/90\\ 4-Y/3-XY^{3}/18-XY^{3}Z^{6}/90\\ 4-Y/3-XY^{3}/18-XY^{3}Z^{6}/90\\ 4-Y/3-XY^{3}/18-XY^{3}Z^{6}/90\\ 3-X/3+2Y/3+XY/3-XY^{3}/18-XY^{3}Z^{6}/90\\ 2+2X/3+2Y/3+XY/3-XY^{3}/18-XY^{3}Z^{6}/90\\ 4-XY/3-XY^{3}/18-XY^{3}Z^{6}/90\\ 4-XY/3-XY^{3}/18-XY^{3}Z^{6}/90\\ 3+2X/3-Y/3+XY/6+XY^{3}/9-XY^{3}Z^{6}/90\\ 4-X/3-XY/6+XY^{3}/9-XY^{3}Z^{6}/90\\ 4-X/3-XY/6+XY^{3}/9-XY^{3}Z^{6}/90\\ 3+2X/3-Y/3+XY/6+XY^{3}/18+2XY^{3}Z^{6}/45\\ 4-X/3-XY/6+XY^{3}/18+2XY^{3}Z^{6}/45\\ 4-X/3-XY/6+XY^{3}/18+2XY^{3}Z^{6}/45\end{array}\right).$$

$${\gamma }^{cat}<{\zeta }^{cat},{\eta }^{cat}<{\theta }^{cat},{\gamma }^{cat}<{\alpha }^{cat}<{\beta }^{cat},{\gamma }^{cat}<{\delta }^{cat},{\gamma }^{cat}<{ϵ}^{cat}.$$

$${\gamma }^{cat}-{\delta }^{cat}=s_7^{cat}-s_8^{cat}-=-2+(2/3)XY+(2/3)X+(2/3)Y.$$

$${\gamma }^{cat}-{ϵ}^{cat}=s_7^{cat}-s_{10}^{cat}=-1+XY/6-X{Y}^3/6+Y=-(Y-1)(X{Y}^2+XY-6)/6<0.$$

$$=-2+XY(1-Y^2/3)/2+X+2Y/3<-2+(1/2)XY(2/3)+X+2Y/3<-2+2/6+1+2/3=0.$$

$$\eta^{cat}-\theta^{cat}=s^{cat}_{13}-s^{cat}_{14}=-1+XY/3+X-Y/3=(Y+3)(X-1)/3<0,\ {\rm i.e.}\ \eta^{cat}<\theta^{cat}.\ \$$

$$\mathbf{s}^{pc}:=M\mathbf{p}^{pc}=\left(\begin{array}[]{l}4-XYZ^{6}/90-Y/3-XY/18\\ 3-XYZ^{6}/90-X/3+2Y/3+5XY/18\\ 4-XYZ^{6}/90-Y/3-XY/18\\ 3-XYZ^{6}/90-X/3+2Y/3+5XY/18\\ 4-XYZ^{6}/90-Y/3-XY/18\\ 4-XYZ^{6}/90-Y/3-XY/18\\ 3-XYZ^{6}/90+2X/3-Y/3+5XY/18\\ 4-XYZ^{6}/90-X/3-XY/18\\ 4-XYZ^{6}/90-X/3-XY/18\\ 3-XYZ^{6}/90+2X/3-Y/3+5XY/18\\ 4-XYZ^{6}/90-X/3-XY/18\\ 4-XYZ^{6}/90-X/3-XY/18\\ 2+2XYZ^{6}/45+2X/3+2Y/3+2XY/9\\ 4+2XYZ^{6}/45-4XY/9\\ 4+2XYZ^{6}/45-4XY/9\end{array}\right).$$

$${\zeta }^{pc}<{\alpha }^{pc},{\zeta }^{pc}<{\beta }^{pc},{\zeta }^{pc}<{\gamma }^{pc},{\zeta }^{pc}<{\delta }^{pc},{\zeta }^{pc}<{\eta }^{pc}.$$

$${\alpha }^{pc}-{\zeta }^{pc}=s_1^{pc}-s_{13}^{pc}=2-XY{Z}^6/18-Y-5XY/18-2X/3.$$

$$1-XY{Z}^6/18-X+XY/18=1-X[Y{Z}^6/18+1-Y/18]>$$

$$1-Y{Z}^6/18-1+Y/18=Y[-Z^6+1]/18>0.$$

$$1-XY{Z}^6/18-Y+XY/18=1-Y[X{Z}^6/18+1-X/18]>$$

$$>-X{Z}^6/18+X/18=X[-Z^6+1]/18>0.$$

$${\delta }^{pc}-{\zeta }^{pc}=s_8^{pc}-s_{13}^{pc}=2-(1/18)XY{Z}^6-X-(5/18)XY-(2/3)Y>0\mathrm{and}$$

$${\eta }^{pc}-{\zeta }^{pc}=s_{14}^{pc}-s_{13}^{pc}=2-(2/3)XY-(2/3)X-(2/3)Y>0,\mathrm{so}{\delta }^{pc}>{\zeta }^{pc}$$

$$\mathbf{s}^{pb}:=M\mathbf{p}^{pb}=\left(\begin{array}[]{l}3+(1/3)ZXY-(1/3)X-(1/15)XY^{3}Z+(2/3)ZY\\ 4-(1/15)XY^{3}Z-(1/3)ZY\\ 4-(1/15)XY^{3}Z-(1/3)ZY\\ 4-(1/15)XY^{3}Z-(1/3)ZY\\ 4-(1/15)XY^{3}Z-(1/3)ZY\\ 3+(1/3)ZXY-(1/3)X-(1/15)XY^{3}Z+(2/3)ZY\\ 2+(1/3)ZXY+(2/3)X-(1/15)XY^{3}Z+(2/3)ZY\\ 4-(1/3)ZXY-(1/15)XY^{3}Z\\ 4-(1/3)ZXY-(1/15)XY^{3}Z\\ 3+(1/6)ZXY+(2/3)X+(1/10)XY^{3}Z-(1/3)ZY\\ 4-(1/6)ZXY-(1/3)X+(1/10)XY^{3}Z\\ 4-(1/6)ZXY-(1/3)X+(1/10)XY^{3}Z\\ 3+(1/6)ZXY+(2/3)X+(1/10)XY^{3}Z-(1/3)ZY\\ 4-(1/6)ZXY-(1/3)X+(1/10)XY^{3}Z\\ 4-(1/6)ZXY-(1/3)X+(1/10)XY^{3}Z\end{array}\right).$$

$${\alpha }^{pb}>{\gamma }^{pb},{\beta }^{pb}>{\gamma }^{pb},{\delta }^{pb}>{\gamma }^{pb},{\zeta }^{pb}>{\gamma }^{pb},{\eta }^{pb}>{\gamma }^{pb}.$$

$${\beta }^{pb}-{\gamma }^{pb}=s_2^{pb}-s_7^{pb}=2-ZY-(1/3)ZXY-(2/3)X>0$$

$${\delta }^{pb}-{\gamma }^{pb}=s_8^{pb}-s_7^{pb}=2-(2/3)ZXY-(2/3)X-(2/3)ZY>0.$$

$$>1-Y(1+X(1-Y^2)/6)>1-Y(1+(1-Y^2)/6)=1-7Y/6+Y^3/6=:h$$

$$(2-x)-ZYX/2-X{Y}^2/6+2/3)>1-ZY(2/3+X(1-Y^2/3)/2)$$

$$>1-Y(2/3+X(1/2-Y^2/6))>1-Y(2/3+1/2-Y^2/6)=1-7Y/6+Y^3/6=h>0.$$

$$r_3=2r_2-1.$$

Full text

Key Words: MRP; parsimony; coalescence theory; evolutionary trees; caterpillar tree; caterpillar subtree; caterpillar measure; caterpillar score; revolution number.

1 Introduction

This article is an extended version of [Miheevs] (http://ceur-ws.org/Vol-2254/10000206.pdf).

Constructing evolutionary species trees is one of the most interesting problems in biology. It means finding the relations and mutual ancestors of existing and extinct species and also the time of formation of new species. Paleontology itself gives very poor information of species trees structure and time lengths of their edges. More precisely species trees can be built by analysis of genomes of species. One considers the set of species $\{a_{i}\}_{1}^{N}$ and their gene groups $\{G^{j}\}_{1}^{K}$, where $G^{j}=\{A^{j}_{i}\}_{i=1}^{N}$ is a set of some functionally relative to each other genes. For example, one group can be responsible for hemoglobin production, another one can define the eye color and so on. Let the gene $A^{j}_{i}$ be discovered in the species $a_{i}$. Then in each $j$-th functional group one can establish the relations between the genes in the form of unrooted tree, where the leaves are the set of elements of $j$-th group. The structure of such trees for different groups can coincide (i.e. be topologically identical) or do not coincide. The number of these coincidences defines the frequency of the corresponding gene tree. These frequencies are the base of parsimony method to construct evolutionary trees which sometimes gives wrong results.

If it is known that a method being applied to some type of problems may fail, why would anyone still use it in this area? Well, in phylogenetics most methods give probabilistic answers. Therefore, getting sometimes wrong answers doesn’t necessarily imply that method is bad. To make a final conclusion about the quality of the method one could estimate how often the wrong results appear. Then one would compare the obtained frequency with the frequencies of other methods’ failures. Having this information on the table, a researcher can decide if the method is acceptable. However, we did better than this. We have found the set of all combinations of parameters of $5$-taxon species tree when Parsimony is guaranteed to fail. Why tree with 5 taxa? It is known that Parsimony always gives right answers on $k$-taxon species tree for $k=3,4$. So, considering $k=5$ is quite logical from computational point of view. Also the smaller $k$ when things go bad, the louder the warning.

The phylogenetics society has intuitive tendency to use Parsimony less and less. Nevertheless, many biologists still do it because of simplicity of the method. They should not be judged for that since simplicity is a strong argument. The results of this paper will show them the danger of Parsimony. However, forewarned is forearmed. If a researcher is sure that their resulting species tree doesn’t have the combination of parameters we showed to be bad, they can safely use fast and simple Parsimony on 5-taxon trees.

For a specific rooted species tree with the known time lengths of the edges, using Coalescence method one can obtain the probabilities of gene trees. That allows to find analytically the rooted species tree and region of lengths of its edges, when parsimony fails.

2 Preliminaries

We consider an evolution tree $T$ (here and after we mean binary tree) of 5 species $a$, $b$, $c$, $d$ and $e$ with some parameters $T_{1}$, $T_{2}$, $T_{3}$ – time in coalescence units between the branching points (see an example in fig. 1).

Based on the Coalescent model [Rosenberg], [SemSte], [Wakeley], [Baum] the program COAL

[DegnanSalter, WangDegnan] yields the probabilities of all 15 possible unrooted gene trees for 5 genes $A$, $B$, $C$, $D$, $E$, such that $A$ is discovered in species $a$, $B$ is discovered in the species $b$ and so on.

All these species have related genes $A$, $B$, $C$, $D$, $E$, respectively. The genes originated from each other or from mutual ancestor. The branching in gene and species trees may not coincide. That creates a problem of inferring species trees from gene trees.

Also, for 5 species there are 15 different unrooted species trees or 105 rooted ones. In each species tree, one can fit any of 15 different gene trees. However, the amount of mutations needed for this fitting will differ generally from one gene tree to another. So, for each gene tree from these 15 and each species tree from the same 15 species trees one can correspond some non-negative integer number of mutations (parsimony score). These can be written in a $15\times 15$ matrix $M$, where rows correspond to species trees and columns correspond to gene trees.

If one knows the frequencies of different gene trees $P=(p_{1},...,p_{15})^{T}$, then the mathematical expectation of the number of mutations for each of 15 possible species trees can be calculated by multiplying the matrix $M$ by the vector-column $P$.

One can assume that the most probable species tree for the given sample of gene trees (or 15 gene trees with assigned probabilities) corresponds to the minimal expectation of mutations. This is the idea of Parsimony method [AllmanRhodes], its application to the problem of inferring species trees from gene trees is called Matrix Representation with Parsimony (MRP) [WangDegnan].

3 Numerical Experiments. Performance of Unrooted MRP for 5-Taxon Species Tree Inference.

Any 5 genes can be joined in one unrooted tree in 15 ways as given in the following list:

Each unrooted tree may be transformed into a rooted tree by introducing a root to an edge. As a result, we have 7 rooted versions for each unrooted tree.

Step 1. Compute parsimony scores by Fitch-Hartigan [Hartigan]

Thus, if $\vec{N}=(N_{1},N_{2},...,N_{15})^{\rm\tiny T}$ is the vector of counts of 15 topological trees in the input, $M$ is the matrix of entries in Table 2 and vector-column  $\vec{S}=(pars(\sigma_{1}),\ pars(\sigma_{2}),...,pars(\sigma_{15}))^{\rm\tiny T}$ then $\vec{S}=M\vec{N}$. Here $(pars(\sigma_{i})$ is parsimony score of species tree on the collection of gene trees $\tau_{1},...,\tau_{15}$.

Step 2. Pick the smallest entry or entries in $\vec{S}$ to determine the most parsimonious tree(s).

To study the 5-taxon case further we need to use Coalescent Theory. The coalescent model, introduced by Kingman in [Kingmn], describes the coalescence of lineages as we move backwards in time within a single species (Note that in biology the understanding of the word ‘species’ may vary. Here we use this word in the same meaning as ‘population’). By “gluing” together such species or populations to form a tree, one gets the Multi-species Coalescent Model, which describes the production of gene trees within species trees.

Let $g_{ij}(T)$ denote the probability that $i$ lineages (genes) since time [math] have coalesced to exactly $j$ lineages at time $T$ under the coalescent model.

General formulas for the $g_{ij}(T)$ were derived in [Tavaré]:

[TABLE]

where  $a_{(k)}=a(a+1)\cdots(a+k-1)$  for  $k\geq 1$  with  $a_{(0)}=1$  (the partial permutation); and$a_{[k]}=a(a-1)\cdots(a-k+1)$  for  $k\geq 1$ with $a_{[0]}=1$.

Some of these formulas for small indexes are

$g_{11}(T)=1,\qquad g_{21}(T)=1-e^{-T},\qquad g_{22}(T)=e^{-T},$

$g_{31}(T)=1-(3/2)e^{-T}+(1/2)e^{-3T},\qquad g_{32}(T)=(3/2)e^{-T}-(3/2)e^{-3T},\qquad g_{33}(T)=e^{-3T},$

$g_{41}(T)=1-(9/5)e^{-T}+e^{-3T}-(1/5)e^{-6T},\qquad g_{42}(T)=(9/5)e^{-T}-3e^{-3T}+(6/5)e^{-6T},$

$g_{43}(T)=2e^{-3T}-2e^{-6T},\qquad g_{44}(T)=e^{-6T}$

Let 3 rooted species tree be  $\Sigma_{1}:=((((a,b):T_{1},c):T_{2},d):T_{3},e)$,  $\Sigma_{2}:=(((a,b):T_{1},(c,d):T_{2}):T_{3},e)$ and  $\Sigma_{3}:=(((a,b):T_{1},c):T_{2},(d,e):T_{3})$.  They are rooted versions of $\sigma_{7},\sigma_{13}$ and $\sigma_{7}$ again, respectively (see Fig. 2).

4 An experiment with caterpillar

Using the program COAL [DegnanSalter, WangDegnan] to get probabilities of rooted gene trees and the formulas (1) for $g_{ij}$, we calculate the probabilities $p_{i}=p(\tau_{i}|\Sigma_{1})$ for $i=\overline{1,15}$, which are listed in Table 3 ($X:=e^{-T_{1}},\ Y:=e^{-T_{2}},\ Z:=e^{-T_{3}},$) after simplification in Maple 15.

Let vector-column $\mathbf{p}$ be $(p_{1},p_{2},\cdots,p_{15})^{\tiny T}$.

We consider the product

[TABLE]

where each entry is the expectation of parsimony score of a possible output tree for MRP (see (2)).

[TABLE]

We discover that in $\mathbf{s}^{cat}$ some entries are equal. Let’s denote them as following

$\alpha^{cat}:=s^{cat}_{1}=s^{cat}_{6}=3-X/3+2Y/3+XY/3-XY^{3}/18-XY^{3}Z^{6}/90$,

$\beta^{cat}:=s^{cat}_{2}=s^{cat}_{3}=s^{cat}_{4}=s^{cat}_{5}=4-Y/3-XY^{3}/18-XY^{3}Z^{6}/90$,

$\gamma^{cat}:=s^{cat}_{7}:=2+2X/3+2Y/3+XY/3-XY^{3}/18-XY^{3}Z^{6}/90$,

$\delta^{cat}:=s^{cat}_{8}=s^{cat}_{9}=4-XY/3-XY^{3}/18-XY^{3}Z^{6}/90,$

$\epsilon^{cat}:=s^{cat}_{10}:=3+2X/3-Y/3+XY/6+XY^{3}/9-XY^{3}Z^{6}/90$,

$\zeta^{cat}:=s^{cat}_{11}=s^{cat}_{12}=4-X/3-XY/6+XY^{3}/9-XY^{3}Z^{6}/90,$

$\eta^{cat}:=s^{cat}_{13}:=3+2X/3-Y/3+XY/6+XY^{3}/18+2XY^{3}Z^{6}/45$,

$\theta^{cat}:=s^{cat}_{14}=s^{cat}_{15}=4-X/3-XY/6+XY^{3}/18+2XY^{3}Z^{6}/45$.

The analytical comparison of these values we form in the following

**Proposition 1.**For any $X,Y,Z\in(0,1)$, the following inequalities hold:

[TABLE]

Proof.

Consider the difference  $\alpha^{cat}-\gamma^{cat}=s^{cat}_{1}-s^{cat}_{7}=-X+1.$ Since  $X\in(0,1)$,  $\gamma^{cat}<\alpha^{cat}$.  Similarly, $\alpha^{cat}-\beta^{cat}=s^{cat}_{1}-s^{cat}_{2}=-1+XY/3-X/3+Y=(Y-1)(X+3)/3<0$ implies $\alpha^{cat}<\beta^{cat}$. Let’s look at

[TABLE]

Since $X,Y$ and $X,Y$ have upper limit 1, $\gamma^{cat}-\delta^{cat}<0$. Observe that

[TABLE]

Thus, $\gamma^{cat}<\epsilon^{cat}$.  Consider  $\gamma^{cat}-\zeta^{cat}=s^{cat}_{7}-s^{cat}_{11}=-2+XY/2-XY^{3}/6+X+2Y/3=$

[TABLE]

Thus,   $\gamma^{cat}<\zeta^{cat}$.  Finally,

[TABLE]

$\square$

However, the 3D-graphs on Figures 3 show that the expressions $\eta^{cat}$ and $\gamma^{cat}$ can not be put in one order for all $X,Y,Z\in(0,1)$.

There is a large region were $\gamma^{cat}<\eta^{cat}$ but, nevertheless, there is also a region where $\eta^{cat}<\gamma^{cat}$. The last defines the parameters $T_{1},T_{2},T_{3}$ where MRP will fail to recover the tree topology of the true species tree producing the gene tree distribution, even when given an arbitrary large sample of gene trees. Figure 3.left shows that provided $Y$ is not too large, regardless of $X,Z$, MRP will return the correct species tree.

To determine this cutoff for $Y$, we set $\eta^{cat}(1,Y,0)=\gamma^{cat}(1,Y,0)$ and solve to get $Y=0.935...$ (the solutions of  $(1/3)Y^{3}-(7/2)Y+3=0$ are  $Y_{1,2,3}\approx 2.670...,-3.605...,0.935...$).

5 An experiment with pseudo-caterpillar

Let the pseudo-caterpillar species tree be $\Sigma_{2}=(((a,b):T_{1},(c,d):T_{2}):T_{3},e)$. We use COAL, the formulas for $g_{ij}$ and Maple 15 to calculate the probabilities $p^{pc}_{i}=p(\tau_{i}|\Sigma_{2})$ for $i=\overline{1,15}$. Their list is shown in Table 4 after simplifications.

In the table 4: $X=e^{-T_{1}}$,  $Y=e^{-T_{2}},\ Z=e^{-T_{3}}$

Now assuming  ${\mathbf{p}}^{pc}=(p_{1}^{pc},p_{2}^{pc},\cdots,p_{15}^{pc})^{\tiny T}$  we see that the product  $\mathbf{s}^{pc}:=M\mathbf{p}^{pc}$ (see (3)) is the vector of expected parsimony scores of possible output trees with a pseudo-caterpillar species tree.

[TABLE]

We discover that in $\mathbf{s}^{pc}$ some entries are equal. Let’s denote them as following

$\alpha^{pc}:=s^{pc}_{1}=s^{pc}_{3}=s^{pc}_{5}=s^{pc}_{6}=4-XYZ^{6}/90-Y/3-XY/18$,

$\beta^{pc}:=s^{pc}_{2}=s^{pc}_{4}=3-XYZ^{6}/90-X/3+2Y/3+5XY/18$,

$\gamma^{pc}:=s^{pc}_{7}=s^{pc}_{10}=3-XYZ^{6}/90+2X/3-Y/3+5XY/18$,

$\delta^{pc}:=s^{pc}_{8}=s^{pc}_{9}=s^{pc}_{11}=s^{pc}_{12}=4-XYZ^{6}/90-X/3-XY/18$,

$\zeta^{pc}:=s^{pc}_{13}=2+2XYZ^{6}/45+2X/3+2Y/3+2XY/9$,

$\eta^{pc}:=s^{pc}_{14}=s^{pc}_{15}=4+2XYZ^{6}/45-4XY/9$.

Let’s compare these expressions.

Proposition 2.  For any $X,Y,Z\in(0,1)$, the following inequalities hold:

[TABLE]

Proof. Consider the difference

[TABLE]

It is easy to see that, since $X,Y,Z\in(0,1)$,  $\alpha^{pc}-\zeta^{pc}>0$, i.e. $\alpha^{pc}>\zeta^{pc}$.

Observe that  $\beta^{pc}-\zeta^{pc}=s^{pc}_{2}-s^{pc}_{13}=$

[TABLE]

[TABLE]

Thus, $\beta^{pc}>\zeta^{pc}$. Similarly,  $\gamma^{pc}-\zeta^{pc}=s^{pc}_{7}-s^{pc}_{13}=$

[TABLE]

[TABLE]

Thus, $\gamma^{pc}>\zeta^{pc}$. Observe

[TABLE]

[TABLE]

and  $\eta^{pc}>\zeta^{pc}.$    $\square$

This implies that if the true species tree is 5-taxon pseudo-caterpillar, MRP, for a sufficiently large data set, will give with probability 1 the unrooted species tree topology for all $T_{1},T_{2},T_{2}\in(0,1)$.

6 An experiment with psuedo-balanced

The last tree we need to consider (since all other are just permutations of taxon names) is the pseudo-balanced species tree $\Sigma_{3}=(((a,b):T_{1},c):T_{2},(d,e):T_{3})$.

The same chain of actions give us the probabilities $p^{pb}_{i}$ (Table 5).

Again, multiplying the matrix $M$ by the vector $\mathbf{p}^{pb}:=(p^{pb}_{1},p^{pb}_{2},\cdots,p^{pb}_{15})$, gives the vector of expected parsimony scores of possible output trees (4) with a pseudo-balanced species tree $\Sigma_{3}$ as input.

[TABLE]

Let’s denote $\alpha^{pb}:=s^{pb}_{1}=s^{pb}_{6}$,   $\beta^{pb}:=s^{pb}_{2}=s^{pb}_{3}=s^{pb}_{4}=s^{pb}_{5}$,   $\gamma^{pb}:=s^{pb}_{7}$,   $\delta^{pb}:=s^{pb}_{8}=s^{pb}_{9}$,  $\zeta^{pb}:=s^{pb}_{10}=s^{pb}_{13}$   and   $\eta^{pb}:=s^{pb}_{11}=s^{pb}_{12}=s^{pb}_{14}=s^{pb}_{15}$.

Proposition 3. For any $X,Y,Z\in(0,1)$, the following inequalities hold:

[TABLE]

Proof. Since $X\in(0,1)$, we immediately have that the difference $\alpha^{pb}-\gamma^{pb}$ equals to $s^{pb}_{1}-s^{pb}_{7}=1-X>0$.  So,  $\alpha^{pb}>\gamma^{pb}$. Also

[TABLE]

and

[TABLE]

Thus,  $\beta^{pb}>\gamma^{pb}$  and  $\delta^{pb}>\gamma^{pb}$.

Consider  $\zeta^{pb}-\gamma^{pb}=s^{pb}_{10}-s^{pb}_{7}=1-ZY(1+(1/6)X(1-Y^{2}))$

[TABLE]

Since  $h^{\prime}=-7/6+(1/2)Y^{2}<0$  and  $h(0)<0$,  $h>0$ for all $Y\in(0,1)$.  Thus,  $\zeta^{pb}>\gamma^{pb}$.

Finally, observe that  $\eta^{pb}-\gamma^{pb}=s^{pb}_{11}-s^{pb}_{7}=$

[TABLE]

[TABLE]

So, $\eta^{pb}>\gamma^{pb}$.    $\square$

So, if the true species tree is 5-taxon pseudo-balanced $\Sigma_{3}$, MRP, for a sufficiently large data set, will give with probability 1 the correct unrooted species tree topology for all $T_{1},T_{2},T_{2}\in(0,1)$.

7 Generalization of results.

7.1 Caterpillar Subtree.

Definition 2. There is a rooted tree $T$ with number of taxa equal to  $|T|=:n$. Let $T_{cat}(T)$ be a caterpillar subtree of this tree. The number $\displaystyle Cat(T):=\max_{T_{cat}(T)\subset T}|T_{cat}(T)|$ for a particular tree $T$ is called caterpillar score for the tree $T$. The number $\displaystyle cat(n):=\min_{|T|=n}\ Cat(T)$ is called caterpillar measure.

It is clear that $cat(n)$ is an increasing function with respect to $n$. There are may be a few consecutive numbers $n$ such that $cat(n)=k$ for some given natural $k$.

Definition 3. Let’s call number $\displaystyle r_{k}:=\min_{cat(n)=k}n$ the revolution number.

Observe that for the caterpillar lengths $1,2,3$ their revolution numbers are $r_{1}=1,r_{2}=2,r_{3}=3$, because these trees are caterpillar themselves. Note, that the third and the second revolution numbers are connected by

[TABLE]

**Theorem 1.***The revolution numbers $r_{k}$ may be calculated recursively  $r_{k}=2r_{k-1}-1$, $k=4,5...$. *

Proof. Let $T$ be a $k$-taxon tree. Since $T$ is a binary tree, we can think of $T$ as two subtrees $T^{\prime},T^{\prime\prime}$ glued together only by two edges at the root (see Figure 4). Observe that for any caterpillar subtree of $T^{\prime}$ one of these $T^{\prime},T^{\prime\prime}$ being transformed properly brings only one edge to the caterpillar.

On the other hand,  $r_{k}$ is an increasing function. So, the first bifurcation in the root will be the worst in the sense of caterpillar score for the tree $T$ when this bifurcation divides the tree into two subtrees $T^{\prime}$ and $T^{\prime\prime}$ such that $|T^{\prime}|=|T^{\prime\prime}|$ for even $|T|$, and \bigl{|}|T^{\prime}|-|T^{\prime\prime}|\bigr{|}=1 for odd $|T|$. Further, we use mathematical induction.

*Base of induction. * The formula (5) is the case $k=2$.

Assumption of induction. Let $r_{3},...,r_{k}$ calculated recursively be revolution numbers.

Inductive step. We need to prove that $r_{k+1}=2r_{k}-1$ is the revolution number. Let us take a tree $T$ with $r_{k+1}$ taxa and consider the worst bifurcation in its root. As mentioned above, the worst bifurcation in the root of $T$ forms two subtrees $T^{\prime}$ and $T^{\prime\prime}$ such that $|T^{\prime}|=r_{k}-1,\ |T^{\prime\prime}|=r_{k}$.

Theorem 1 allows to continue the sequence of the revolution numbers for the caterpillar measures $3,4,5,6,7,...$ as $r_{k}=3,5,9,17,33,...$, respectively. For trees with number of taxa in $[r_{k},2r_{k}-1]$, the caterpillar measure is $k$.

7.2 MRP on trees with the number of taxa greater than 5

Theorem 2.* If a true species rooted tree $\hat{G}$ contains 5-taxon caterpillar subtree, then MRP may fail to obtain the unrooted version of $\hat{G}$ from the set of gene trees generated by Coalescent model from $\hat{G}$. *

Proof. For the number of species greater than 5 and the same number of genes one can make the following construction.

Take a caterpillar tree $\Gamma$ of 5 species $a$, $b$, $c$, $d$, $e$ with $T_{1}$, $T_{2}$, $T_{3}$ and root $\rho$, such that parsimony fails (fig. 3). Take an arbitrary tree $G$, where every edge is very close to 0 ($T_{i}\approx 0$,  $i>4$). Connect $\Gamma$ to $G$ through its root $\rho$ and the edge $\epsilon$ with the length $T_{4}$. Make $T_{4}$ big enough so the genes $A$, $B$, $C$, $D$, $E$ coalesce in $\epsilon$ if they didn’t in $\Gamma$. No matter what is on the upper end of $\epsilon$, root of entire tree or inner node created by $\epsilon$ on some edge of $G$ .

The numeration of $n$ possible gene trees we do in the following way: First 15 trees will have the same subtree $G$ and different topology or permutation of $A$, $B$, $C$, $D$, $E$. Other $n-16$ trees can be numerated in any order, and we set their probabilities $p_{16},...,p_{n}$ equal to zero, since $T_{i}$,  $i\in\{5,...,n\}$, can be taken infinitively small. Therefore, the set of gene trees is numerated and the probabilities of them are presented by vector-column $p=(p_{1},...,p_{15},p_{16},...,p_{n})^{T}$, where all the entries below 15-th equal zero and the first 15 are the same that obtained ones for 5-taxon experiment. The matrix $M_{n\times n}$ has dimension $n\times n$, but only submatrix $M_{n\times 15}$ does participate in calculation of expectations of gene mutations due to $p_{16},...,p_{n}=0$.

Moreover, the parsimony incorrect choice may be shown on submatrix $M_{15\times 15}$ in upper corner. Observe that the elements of $M_{15\times 15}$ are the sums of the elements of $M$ obtained earlier in performance of MRP for 5-taxon trees (see table 2) and some constant number generated by the constant subtree $G$ with coalesced gene $A+B+C+D+E$. This means that each of $s^{cat}_{1}$, $s^{cat}_{2}$, … , $s^{cat}_{15}$ from Section 2 must be increased by the constant value   $V\sum_{1}^{15}p_{i}$,  to be a new mathematical expectation for the new big tree. So, the minimum among the first 15 rows must be achieved in the same index. Therefore, being wrong for 5-taxon caterpillar species tree the parsimony becomes wrong for the constructed tree $\hat{G}$ as well.     $\square$

Corollary 1. *MRP on a set of gene trees with 5 taxa or more may yield wrong result. If one applies MRP on set of gene trees with 9 taxa or more, MRP may fail even more probably, since 9-taxon species tree always has a caterpillar subtree, which may have unfortunate lengths of inner edges from the small region in Fig. 3. *

We have established what these unfortunate lengths are. But how to find these caterpillar subtrees? One may use the following

Theorem 3. *If a tree has three consecutive inner edges not contaning the root betwen them but perhaps one of these edges ending on it then the tree has a 5-taxon caterpillar subtree, which contains these three inner edges. *

Proof. At Fig. 5 we can see three consecutive edges denoted 1,2 and 3 between nodes $a$,$b$, $c$ and $d$, respectively and the third edge is closest to the root. Since these edges are inner, all the nodes must be bifurcation points and so each of $b$, $c$, $d$ has one more edge running towards a taxon or a clade opposite to the root and $a$ has two more such edges. Let $d$ be the root then simply contracting each of the mentioned clades to one of its taxa we get a 5-taxon caterpillar subtree. If $d$ is not the root, throw away one of the edges at $d$ to make $d$ the root of the 5-taxon caterpillar subtree.

**Corollary 2.***It is enough to know that in a species tree with amount of taxa 5 or more there are three consecutive inner edges not going through the root but perhaps ending on it with lengths $T_{1},\ T_{2}$ and $T_{3}$ from the small region of cube in Figure 3 to conclude that Parsimony is guaranteed to fail on this tree. *

8 Conclusions and Future work

The fact that Parsimony may fail is not new. However, here we proved that no mater of what topology the true 9-taxon and greater species tree is the only condition to fail Parsimony is to have in this tree three consecutive inner edges not going through the root but perhaps ending on it with lengths $T1,T2,T3$ (which are times in coalescence units between the branching points) of some proportions. Obviously, the probability to meet these lengths is growing in general with the size of species tree. So, if one wants to safely use MRP on a set of $n$-taxon gene trees, it is need to know somehow that the resulting $n$-taxon species tree cannot have any of “bad” topologies and edge lengths from “bad” regions. This paper makes it possible for $n\leq 5$. Also, one may apply MRP on a set of 5-taxon gene trees and if the result is the caterpillar tree or topology $\Sigma_{3}$ from Fig. 2, it is true.

One may consider 6-taxon species trees the way we did in this paper and check the existence of 6-taxon topology which forces Parsimony to fail when lengths of inner edges have some proportions. Then prove, perhaps following our ideas, that every tree with some number of taxa greater than 6 has this 6-taxon topology subgraph. Then do the same for 7 taxa and so on.

If one studies all “bad” parameters’ regions of all “bad” topologies for all trees with the amount of taxa less or equal some $n$, it becomes theoretically possible to check either MRP can be applied for a set of gene trees of $n$ taxa (if, of course, the researcher knows enough information about possible results). However, taking into account the factorial growth of the amount of binary trees with respect to the amount of taxa, the problem to find “safe zone” for MRP becomes extremely hard. Unfortunately, we don’t see any way around besides doing that scheme for each $k<n$. So, someone has to be very motivated to use MRP to go through with the research. Nowadays, there are good coalescence based methods, for example, [YufengWu], [RanYan] and [EmmsKelly]

It could be interesting to study stability questions of the coalescent model with uncertainty in data applying the thoughts from [Zubov].