<figure><img src="image_1.png" /></figure>

HAL
open science

# Algorithms to reconstruct past indels: the deletion-only parsimony problem

Jordan Moutet, Eric Rivals, Fabio Pardi

<figure><img src="image_2.png" /></figure>

## To cite this version:

Jordan Moutet, Eric Rivals, Fabio Pardi. Algorithms to reconstruct past indels: the deletion-only parsimony problem. PLoS Computational Biology, 2025, 21 (7), pp.e1012585. 10.1371/journal.pcbi.1012585. lirmm-05212128

HAL Id: lirmm-05212128

https://hal-lirmm.ccsd.cnrs.fr/lirmm-05212128v1

Submitted on 17 Aug 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Algorithms to reconstruct past indels:
the deletion-only parsimony problem

Jordan Moutet¹, Eric Rivals¹, Fabio Pardi¹*

¹LIRMM, Université de Montpellier, CNRS, Montpellier, France

*pardi@lirmm.fr

## Abstract

Ancestral sequence reconstruction is an important task in bioinformatics, with applications ranging from protein engineering to the study of genome evolution. When sequences can only undergo substitutions, optimal reconstructions can be efficiently computed using well-known algorithms. However, accounting for indels in ancestral reconstructions is much harder. First, for biologically-relevant problem formulations, no polynomial-time exact algorithms are available. Second, multiple reconstructions are often equally parsimonious or likely, making it crucial to correctly display uncertainty in the results.

Here, we consider a parsimony approach where only deletions are allowed, while addressing the aforementioned limitations. First, we describe an exact algorithm to obtain all the optimal solutions. The algorithm runs in polynomial time if only one solution is sought. Second, we show that all possible optimal reconstructions for a fixed node can be represented using a graph computable in polynomial time. While previous studies have proposed graph-based representations of ancestral reconstructions, this result is the first to offer a solid mathematical justification for this approach. Finally we provide arguments for the relevance of the deletion-only case for the general case.

**Author summary:** An exciting frontier in evolutionary biology is the ability to reconstruct DNA or protein sequences from species that lived in the distant past. By analyzing sequences from present-day species, we aim to infer the sequences of their common ancestors —a process known as ancestral sequence reconstruction. This task has far-reaching applications, such as resurrecting ancient proteins and studying the biology of extinct organisms. However, a significant challenge remains: the lack of well-established methods for inferring past deletions and insertions —mutations that remove or add segments of genetic code. In this paper, we present algorithms that lay the groundwork for addressing this gap. We show that finding the reconstructions involving only deletion events, while minimizing their number, can be done efficiently. Additionally, we show that all optimal solutions can be represented using specialized graphs. While previous studies have proposed graph-based representations of ancestral reconstructions, we are the first to provide a rigorous mathematical foundation for the use of these graphs.
# 1 Introduction

Biological sequences undergo mutations during evolution, and the most common ones are substitutions (when a nucleotide or amino acid is replaced by another), insertions (when a sequence is added inside another larger sequence), and deletions (when a part of a sequence is lost). Ancestral Sequence Reconstruction (ASR) is the task that consists in recovering the evolutionary history of a set of sequences related by a known phylogeny, that is the mutation events that explain the differences between these sequences. Finding where exactly these events occurred in the phylogeny amounts to determining the sequences that were present at the internal nodes of the phylogeny, that is, the ancestral sequences. ASR is a key task in bioinformatics, with wide-ranging applications, from protein engineering in biomedicine [1] and biotechnology [2], to protein resurrection for vaccine development [3] and paleontology [4], to the reconstruction of ancient genomes [5]. Additionally, ASR lies at the foundation of other tasks in bioinformatics, for example sequence alignment [6] and phylogenetic placement [7, 8], meaning that any improvements to the methodology for ASR are likely to benefit these other subjects.

Polynomial-time algorithms to infer past substitutions were introduced decades ago, and now form a solid basis for ASR [9, 10, 11, 12, 13, 14]. On the other hand, the inference of indels (insertions and deletions), despite multiple methodological attempts, does not have an established theoretical basis [6]. The lack of broadly adopted solutions for indels compared to substitutions can be explained by multiple factors. First, no exact algorithms that are polynomial both in the number of sequences and their length are known to solve the most biologically-relevant formulations of this problem. Thus, some algorithms are exact but have worst-case running times exponential either in the number of sequences [15, 16, 17] or their length [18, 19]. Usually, the available polynomial-time algorithms are heuristic [5, 16, 20], or they are based on simplified modeling approaches such as indel coding and its variants [21, 22, 23] or the assumption of site independence [24, 25]. Another challenge for the inference of past indels is that multiple reconstructions are often equally parsimonious/likely, meaning that it can be misleading to just return a single optimal solution. Instead, it is important to provide some insight about the other reconstructions. Exploring the diversity of solutions not only carries an additional computational cost, but it also poses the problem of how to correctly display uncertainty in the results. Even though this issue is also present for the inference of past substitutions, it is much more problematic for indels, as they can overlap each other. A graph-based solution to represent multiple reconstructions for a single node in the phylogeny, showing multiple alternative starting and ending positions for each putative gap, was recently put forward by Foley et al. [23]. Despite the problems outlined above, indel inference shows a rebirth of interest, with two parsimony-based works published in the past year [25, 26].

Similarly to substitutions, indels can be inferred using two broad approaches. The conceptually simpler one is maximum parsimony. Here, a number of formulations of the problem were already proposed two decades ago [18, 27, 15]. Although similar in spirit, these approaches differed in terms of some important modeling choices (e.g. on how to count indels between two sequences, and on whether to allow multiple introductions of the same character). The formu-
lation by Chindelevitch et al. [27] is arguably the most faithful to the biological context, and was later adapted into number of maximum-likelihood approaches for indel inference [28, 16, 29]. It leads to a problem called the Indel Parsimony Problem (IPP), which we formally define below (Def. 5). Although the IPP is NP-hard for trees of unbounded degree, its complexity remains unknown for trees with bounded degree, in particular binary trees, and under the unit-cost model, consisting of minimizing the number of indels irrespective of their size or place in the tree.

The second broad approach for indel inference is the statistical one. In practice one assumes a probabilistic model explicitly describing the evolution of indels (e.g. [30, 31, 32, 33, 34]) and then seeks the evolutionary history or histories with highest probability [19, 16, 29, 35]. This idea can be formulated either as an optimization problem or in a Bayesian context [36, 24]. For many tasks in phylogenetics, maximum parsimony has a computational advantage over statistical methods, which may rapidly become considerable as the number of sequences grows [37, 38]. It is reasonable to expect a similar advantage for indel inference. Moreover statistical approaches may benefit from the insights gained from studying the algorithmics of maximum parsimony (see, e.g., the similarities between Sankoff's and Felsenstein's algorithms [10, 39]).

In this paper, we focus on the specific case of parsimony-based ASR where solutions are constrained to only involve deletions. This problem, which we call the deletion-only parsimony problem (DPP), was first introduced by Chindelevitch et al. [27]. As we discuss more extensively in Section 7, the DPP is relevant to the solution of the general problem (the IPP mentioned above), as parts of the IPP can be optimally solved alongside this problem. The first contribution of this paper is the proposal of an exact algorithm to find all optimal solutions of the DPP. Since the number of solutions may grow exponentially in the size of the input alignment, our algorithm will also run in exponential time. However, it is easy to modify the algorithm so that it can find a single optimal solution in polynomial time. Naturally, we formally prove the correctness of this algorithm (i.e. that all, and only, the optimal solutions of the DPP are found). Secondly, we propose a polynomial alternative for the exponential-time part of the algorithm, constructing a graph for each internal node of the phylogenetic tree. Each of these graphs represents exactly the set of all reconstructions that can be found at the specified node in an optimal solution for the whole tree. These graphs are similar to the ones constructed by Foley et al. [23] and Tule et al. [26], also known as *partial order graphs* (POGs; [40]). We discuss the relation with these and other recent works in Section 8.

After this introduction, Section 2 settles the notations and definitions that are necessary to introduce the IPP and the DPP. Then, Section 3 describes our algorithm, and illustrates it step by step with a recurring example. Section 4 presents the theoretical results about the algorithm's correctness. Then, in Section 5, we describe and illustrate our graph-based representation of optimal solutions, and the algorithm to construct it. Afterwards, in Section 6, we analyze the asymptotic complexity of all the provided algorithms. Section 7 deals with the issue of the relevance of the DPP for the solution of the IPP, and, finally Section 8 concludes with some remarks about our original contributions, and the questions left open. The proofs of all our mathematical statements and a
number of other accessory results are given as Supporting Information.

# 2 Preliminaries

## 2.1 Basic notation

**Notation for phylogenetic trees.** We define a (phylogenetic) tree $T$ as a binary rooted tree, whose root, set of edges, set of nodes and set of leaves are denoted by $r(T)$, $E(T)$, $V(T)$ and $L(T)$, respectively. Given two nodes $u, v \in V(T)$, we say that $v$ is a *descendant* of $u$ and that $u$ is an *ancestor* of $v$ if the path from $r(T)$ to $v$ traverses $u$; in particular, if $v \neq u$ we say that $v$ is a *strict descendant* of $u$, and that $u$ is a *strict ancestor* of $v$. $T_u$ denotes the subtree of $T$ induced by the set of descendants of $u$. Finally, if $(p, w) \in E(T)$ and $w$ is a descendant of $p$, we say that $p$ is the *parent* of $w$, and that $w$ is a *child* of $p$.

**Notation and assumptions for alignments.** Given a set $X$, a (multiple sequence) alignment $A$ for $X$ is a table whose columns are labeled by a set of consecutive integers, called *sites*, and whose rows are in a one-to-one correspondence with the elements of $X$, called *taxa*. In particular, given $x \in X$ and site $i$:

* $A_x$ denotes the row of $A$ corresponding to $x$;

* $A[i]$ denotes the column $i$ of $A$;

* $A_x[i]$ denotes the element of $A$ at column $i$ and at the row corresponding to $x$.

Although in general $A_x[i]$ may be an element of any fixed alphabet, throughout this paper we assume $A$ has been preprocessed so as to be binary, i.e. $A_x[i] \in \{0, 1\}$, where $1$ and $0$ are used to represent presence and absence of a character at site $i$, respectively. For any row $A_x$ of an alignment, $\mathbb{1}(A_x)$ denotes the set of sites that are set to $1$ (i.e. filled), that is $\mathbb{1}(A_x) = \{k : A_x[k] = 1\}$.

For any alignment in input, we also assume that it has been preprocessed by adding a first column $A[0]$ and a last column $A[m+1]$, with both $A[0]$ and $A[m+1]$ only containing $1$s, and where $m$ is the number of columns in the original alignment. Finally we assume that no column $A[i]$ contains only $0$s.

**Notation about intervals and gaps.** Given two integers $i$ and $j$, the set of all integers $k$ with $i \le k \le j$ is denoted by $[i, j]$ and is called an *interval*. Two intervals $[i, j]$ and $[i', j']$ are said to *overlap* if $[i, j] \cap [i', j'] \neq \emptyset$; in particular, if $[i, j] \not\subseteq [i', j']$ and $[i, j] \not\supseteq [i', j']$, then $[i, j]$ and $[i', j']$ are said to *partially overlap*. Moreover, $A_x[i, j]$ denotes the sequence of characters at row $x$ in alignment $A$ between (and including) sites $i$ and $j$.

**Definition 1 (gap).** Let $[i, j] \subseteq [1, m]$. $A_x[i, j]$ is a gap if all characters in $A_x[i, j]$ are $0$s and $A_x[i-1] = A_x[j+1] = 1$.

## 2.2 Alignment extensions and indels

Alignments encountered in practice usually only contain information for the leaves of a phylogenetic tree and are thus alignments for $L(T)$ where $T$ is an
unknown tree. Given a putative tree $T$ and an alignment for $L(T)$, the goal of
ancestral indel reconstruction is to recover the pattern of gaps at the internal
nodes of $T$. This idea underlies the following definitions.

**Definition 2 (alignment extension).** Let $T$ be a tree and $A$ an alignment for $L(T)$. An alignment extension of $A$ for $T$ is an alignment $A^+$ for $V(T)$, such that $A_v^+ = A_v$ for any $v \in L(T)$.

**Definition 3 (phylogenetic correctness).** An alignment extension $A^+$ for $T$ is said to be phylogenetically correct if for any two nodes $u, v \in V(T)$ and any site $i$, such that $A_u^+[i] = A_v^+[i] = 1$, we have $A_x^+[i] = 1$ for each node $x$ in the path between $u$ and $v$.

The biological intuition behind phylogenetic correctness is the following. If two nodes $u, v$ satisfying $A_u^+[i] = A_v^+[i] = 1$ were separated by a node $x$ with $A_x^+[i] = 0$, that would imply that either the character at site $i$ has been introduced by two separate insertion events, or that it has first been lost with a deletion and then re-introduced with an insertion. Both of these scenarios contradict the phylogenetic interpretation of alignments, in which two characters should only be placed at the same site in an alignment if they derived from the same ancestral character. For these reasons, the problems that we consider in this paper (Definitions 5 and 6 below) constrain alignment extensions to be phylogenetically correct. This constraint is equivalent to the common interpretation of “Dollo’s law” in phylogenetics, which asserts that no (complex) character can be introduced more than once during the course of evolution [41]. Recent works impose phylogenetic correctness by presupposing Dollo’s law [25].

The last missing ingredient before we can introduce the problems to solve is a definition of insertions and deletions, and how to count them. In a parsimony context, this can be done as in Definition 4, which is equivalent to the one by Chindelevitch et al. [27].

**Definition 4 (insertion, deletion, unit-cost distance).** Let $T$ be a tree and $(u, v) \in E(T)$. Let $A$ be an alignment for $X \supseteq \{u, v\}$. We say that $A$ has a *deletion* from $u$ to $v$ at interval $[i, j]$ if the following statements are all verified

1. $(A_u[i], A_v[i]) = (1, 0)$

2. $(A_u[j], A_v[j]) = (1, 0)$

3. $(A_u[k], A_v[k]) \neq (1, 1), \forall k \in [i, j]$

4. $[i, j]$ is maximal among all intervals satisfying the three points above.

Similarly, we say that $A$ has an *insertion* from $u$ to $v$ at interval $[i, j]$ if there is a deletion from $v$ to $u$ at $[i, j]$. Together, insertions and deletions are called *indels*. Sites $i$ and $j$ are called the *endpoints* of the indel. The number of indels from $u$ to $v$ is also called the *unit-cost distance* between $u$ and $v$ and here will be noted as $d(A_u, A_v)$.

**Example 1.** Let $A_u$ and $A_v$ be defined as follows:

0 1 2 3 4 5 6 7
$u$: 1 1 0 0 1 1 1 1
$v$: 1 0 0 1 0 1 0 1
Here we have $d(A_u, A_v) = 3$, and more specifically a deletion at $[1, 4]$, an insertion at $\{3\} = [3, 3]$, and a deletion at $\{6\} = [6, 6]$. We use this example to illustrate the biological appropriateness of Definition 4, supported by a few informal arguments.

First, an interesting feature is that insertions and deletions can occur at overlapping intervals (e.g. $[1, 4]$ and $[3, 3]$ here). Indeed, other ways to define and count indels may result in considering the character losses at sites 1 and 4 as two separate deletions, implying a number of indels larger than 3 in this example [29, 15, 17]. However, this seems unjustified in a maximum parsimony context: we can imagine that at an intermediate node $x$ (subdividing edge $(u, v)$ into $(u, x)$ and $(x, v)$) we may have $A_x = 10000101$, which can be explained with a single deletion acting on both sites 1 and 4 simultaneously from $u$ to $x$ (this is possible because sites 1 and 4 contain consecutive characters in the sequence at $u$). Additionally $A_x$ implies a deletion at site 6 from $u$ to $x$ and an insertion from $x$ to $v$ at site 3. The existence of a scenario employing only 3 indels implies that, under maximum parsimony, $d(A_u, A_v) \le 3$.

Moreover, any definition resulting in a number of indels smaller than 3 is also biologically implausible: since from $u$ to $v$ some characters are lost (at sites 1, 4, 6) and some other characters are gained (site 3) there must be at least one deletion and one insertion. Let us first assume that this is possible and that the insertion occurs before the deletion. Then the single deletion would have to simultaneously erase sites 1, 4, 6 without erasing intermediate sites 3 and 5. This is biologically unrealistic, as a single deletion should only be allowed to erase a set of consecutive characters. Let us now assume the opposite scenario, where the insertion occurs after the deletion. Since the deletion must remove the characters at sites 1 and 6 and all characters that are present in $u$ between them, this would result in the loss of the character at site 5, followed by its reintroduction. As argued for phylogenetic correctness, this should not be allowed to happen under a biologically-relevant modelisation of indels (see paragraph following Def. 3).

## 2.3 Problem definition

Chindelevitch et al. [27] introduced the following problem and showed that it is NP-complete on trees of unbounded degree.

**Definition 5 (INDEL PARSIMONY PROBLEM).**

*Input:* A tree $T$ and an alignment $A$ for $L(T)$.

*Output:* $A^+$, a phylogenetically correct extension of $A$, minimizing the cost

$$c(A^+) = \sum_{(u,v) \in E(T)} d(A_u^+, A_v^+).$$

Here, we focus on a variant of this problem, where insertions are not allowed. We do this by specifying an additional constraint on $A^+$:

**Definition 6 (DELETION-ONLY PARSIMONY PROBLEM – DPP).**

*Input:* A tree $T$ and an alignment $A$ for $L(T)$.

*Output:* $A^+$, an alignment extension of $A$ minimizing $c(A^+)$ under the constraint that, for all $(p, w) \in E(T)$ where $p$ is the parent of $w$, all indels from $p$ to $w$ must be deletions.
We will refer to the extensions $A^+$ that satisfy the constraint of the DPP as the *candidate solutions* of the DPP. Fig. 1 and 2 show an example of the input and of a candidate solution of the DPP, respectively. We note that the candidate solution in Fig. 2 is not optimal for the DPP, as we will show that there exist several solutions with a strictly lower $c(A^+)$.

**Observation 1.** Let $A^+$ be a candidate solution of the DPP. Then, for every node $w$ and site $k$:

1. If $A_w^+[k] = 0$ then $A_u^+[k] = 0, \forall u \in V(T_w)$.

2. If $A_w^+[k] = 1$ then $A_u^+[k] = 1, \forall u$ in the path between $r(T)$ and $w$.

3. $A^+$ is phylogenetically correct.

4. $A_{r(T)}^+[k] = 1$.

Points 1, 2 in Observation 1 are due to the fact that, because insertions are not allowed, no site $k$ can be such that $A_p^+[k] = 0$ and $A_w^+[k] = 1$ for any $p, w$ where $p$ is the parent of $w$. Point 3 follows from the fact that any violation of phylogenetic correctness would also imply the existence of one such site. Finally, point 4 is due to the fact that if $A_{r(T)}^+[k] = 0$, then a column of the input alignment $A$ would only consist of gaps (0s).

<figure><img src="image_4.png" /><figcaption>Figure 1: Example input of the DPP: a phylogenetic tree $T$ and an alignment $A$ for $L(T)$.</figcaption></figure>

<table><thead><tr><th></th><th>0</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th><th>10</th><th>11</th><th>12</th></tr></thead><tbody><tr><th>r</th><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td></tr><tr><th>a</th><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td></tr><tr><th>x</th><td>1</td><td>1</td><td>0</td><td>1</td><td>1</td><td>1</td><td>0</td><td>1</td><td>1</td><td>0</td><td>0</td><td>1</td><td>1</td></tr><tr><th>b</th><td>1</td><td>1</td><td>0</td><td>1</td><td>1</td><td>1</td><td>0</td><td>1</td><td>1</td><td>0</td><td>0</td><td>1</td><td>1</td></tr><tr><th>y</th><td>1</td><td>1</td><td>0</td><td>0</td><td>1</td><td>0</td><td>0</td><td>0</td><td>1</td><td>1</td><td>0</td><td>0</td><td>1</td></tr><tr><th>c</th><td>1</td><td>1</td><td>0</td><td>0</td><td>1</td><td>0</td><td>0</td><td>0</td><td>1</td><td>1</td><td>0</td><td>0</td><td>1</td></tr><tr><th>d</th><td>1</td><td>0</td><td>0</td><td>0</td><td>1</td><td>0</td><td>0</td><td>0</td><td>0</td><td>1</td><td>0</td><td>0</td><td>1</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>Figure 2: (Right:) A candidate solution $A^+$ for the example in Fig. 1. The rows shaded in grey are those that have been added when going from the input $A$ to the extension $A^+$. They correspond to the internal nodes of $T$. (Left:) The $c(A^+) = 9$ deletions in $A^+$ are here depicted as black boxes appearing on the edges where the deletions occur, along with their endpoints. The relative order of these boxes on their edge is uninformative.</figcaption></figure>

<table><thead><tr><th></th><th>0</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th><th>10</th><th>11</th><th>12</th></tr></thead><tbody><tr><th>r</th><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td></tr><tr><th>a</th><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td></tr><tr><th>x</th><td>1</td><td>1</td><td>0</td><td>1</td><td>1</td><td>1</td><td>0</td><td>1</td><td>1</td><td>1</td><td>0</td><td>1</td><td>1</td></tr><tr><th>b</th><td>1</td><td>1</td><td>0</td><td>1</td><td>1</td><td>1</td><td>0</td><td>1</td><td>1</td><td>0</td><td>0</td><td>1</td><td>1</td></tr><tr><th>y</th><td>1</td><td>1</td><td>0</td><td>0</td><td>1</td><td>0</td><td>0</td><td>0</td><td>1</td><td>1</td><td>0</td><td>0</td><td>1</td></tr><tr><th>c</th><td>1</td><td>1</td><td>0</td><td>0</td><td>1</td><td>0</td><td>0</td><td>0</td><td>1</td><td>1</td><td>0</td><td>0</td><td>1</td></tr><tr><th>d</th><td>1</td><td>0</td><td>0</td><td>0</td><td>1</td><td>0</td><td>0</td><td>0</td><td>0</td><td>1</td><td>0</td><td>0</td><td>1</td></tr></tbody></table>

# 3 The deletion-only algorithm

We now describe an algorithm that solves the DPP (Def. 6). The algorithm is designed to find not just one optimal solution, but all of them.
**High-level description** Our algorithm for the DPP consist of two phases. The first phase is a bottom-up procedure, starting from the leaves and proceeding up to the root of the tree, conveying information about the presence of gaps in the descendants of each node. By doing so, we can determine some gaps that are necessarily present in every optimal solution (called 0-gaps, see next paragraph), and take note of the parts that can only be determined later on. The second phase is a top-down procedure, from the root to the leaves. It allows us to resolve the parts undetermined during the first phase, sometimes with multiple choices, providing a way to generate every possible optimal solution. Algorithm 1 provides the pseudo code for our algorithm, where procedures BottomUp and TopDown are defined further below (Algorithms 2 and 3 respectively).

**Algorithm-specific notation.** For each node we store a set of labeled gaps defined as intervals $[i, j] \subseteq [1, m]$ labeled with one of three symbols: 0, C or P. If an interval is labeled by 0, C, P we call it 0-gap, C-gap or P-gap, respectively. We explain the choice of this naming later below. If a 0-gap $[i, j]$ is stored for node $u$, we also say that $u$ has a 0-gap $[i, j]$ or that $[i, j]$ is a 0-gap at $u$. The same applies to C-gaps and P-gaps.

Remark that labeled gaps are those stored by the algorithm during its execution, and should not be confused with the gaps *tout-court* which are those that can be observed in an alignment, such as the input $A$ or the output $A^+$ (Definition 1).

**Bottom-up phase.** We now describe how labeled gaps are created by the algorithm for every node $w$, in a bottom-up fashion. Algorithm 2 provides a recursive implementation of the bottom-up phase. If $w$ is a leaf, its labeled gaps are set to the maximal intervals of sites filled with 0s in $A_w$ and labeled as 0-gaps. In the example of Fig. 1 leaf $d$ has three 0-gaps: [1, 3], [5, 8] and [10, 11].

Now let $w$ be an internal node and $u$ and $v$ its children. For every pair of overlapping labeled gaps $[i_u, j_u]$, $[i_v, j_v]$ at $u$ and $v$, respectively, we store the labeled gap $[i, j] = [i_u, j_u] \cap [i_v, j_v]$ for node $w$ and set its label with one of the following rules.

**R1.1** If $[i_u, j_u] = [i_v, j_v]$ and neither of them is a P-gap at $u$ and $v$, then $[i, j]$ is labeled as a 0-gap at $w$.

<caption>**Algorithm 1:** DeletionOnly($T$, $A$)</caption>

**Data:** $T$ a binary tree; $A$ an alignment for $L(T)$
**Result:** The set of all alignment extensions optimal for the DPP
1. Let $S$ be a data structure containing, $\forall u \in V(T)$, a set $S_u$ of labeled gaps at $u$
   // initially $S_u = \emptyset, \forall u \in V(T)$
2. BottomUp($r(T)$, $A$, $S$)
   // this fills $S_u, \forall u \in V(T)$
3. $Sol \leftarrow \{S\}$
4. TopDown($r(T)$, $Sol$)
   // this fills $Sol$
   // at this point $\forall S \in Sol, \forall u \in V(T)$, all labeled gaps in $S_u$ are 0-gaps
5. **Return** $Sol$
<table><caption>**Algorithm 2:** BottomUp($w$, $A$, $S$)</caption><thead><tr><th colspan="2">Data:</th><td>$w$, a node of the binary tree $T$; $A$, an alignment for $L(T)$; $S$, a data structure containing a set $S_u$ of labeled gaps at $u$, $\forall u \in V(T)$</td></tr><tr><th colspan="2">Result:</th><td>No return value, but the sets $S_z$ are filled $\forall z \in V(T_w)$</td></tr></thead><tbody><tr><td>1</td><td>if $w$ is a leaf then</td><td></td></tr><tr><td>2</td><td>S<sub>w</sub> &larr; gaps in A<sub>w</sub></td><td></td></tr><tr><td>3</td><td>for [i, j] &isin; S<sub>w</sub> do</td><td></td></tr><tr><td>4</td><td>S<sub>w</sub> &larr; {[i<sub>u</sub>, j<sub>u</sub>] &cap; [i<sub>v</sub>, j<sub>v</sub>]} : ([i<sub>u</sub>, j<sub>u</sub>], [i<sub>v</sub>, j<sub>v</sub>]) &isin; S<sub>u</sub> &times; S<sub>v</sub>} &setminus; {0}</td><td></td></tr><tr><td>5</td><td>else</td><td></td></tr><tr><td>6</td><td>Let u, v be the two children of w</td><td></td></tr><tr><td>7</td><td>BottomUp(u, A, S)</td><td></td></tr><tr><td>8</td><td>BottomUp(v, A, S)</td><td></td></tr><tr><td>9</td><td>S<sub>w</sub> &larr; {[i<sub>u</sub>, j<sub>u</sub>] &cap; [i<sub>v</sub>, j<sub>v</sub>]} : ([i<sub>u</sub>, j<sub>u</sub>], [i<sub>v</sub>, j<sub>v</sub>]) &isin; S<sub>u</sub> &times; S<sub>v</sub>} &setminus; {0}</td><td></td></tr><tr><td>10</td><td>for [i, j] &isin; S<sub>w</sub> do</td><td></td></tr><tr><td>11</td><td>assign label of [i, j] considering rules R1.1 to R1.6</td><td></td></tr></tbody></table>

**R1.2** If $[i_u, j_u] = [i_v, j_v]$ and only one of them is a P-gap at $u$ or $v$, then $[i, j]$ is labeled as a C-gap at $w$.

**R1.3** If $[i_u, j_u] = [i_v, j_v]$ and both are P-gaps at $u$ and $v$, then $[i, j]$ is labeled as a P-gap at $w$.

**R1.4** If $[i_u, j_u] \subsetneq [i_v, j_v]$ and $[i_u, j_u]$ is not a P-gap at $u$, then $[i, j]$ is labeled as a C-gap at $w$. The same applies if the roles of $u$ and $v$ are reversed.

**R1.5** If $[i_u, j_u] \subsetneq [i_v, j_v]$ and $[i_u, j_u]$ is a P-gap at $u$, then $[i, j]$ is labeled as a P-gap at $w$. The same applies if the roles of $u$ and $v$ are reversed.

**R1.6** If $[i_u, j_u]$ and $[i_v, j_v]$ partially overlap, then $[i, j]$ is labeled as a P-gap at $w$.

These rules cover every case and are exclusive. In fact, R1.1, R1.2, R1.3 obviously cover all equality cases, R1.4 and R1.5 cover the strict inclusion cases, and finally R1.6 covers the partial overlap case. Table 1 displays these labeling rules as a function of the labels of $[i_u, j_u]$, $[i_v, j_v]$, and of whether the two gaps are identical, strictly included in one another, or partially overlapping.

<table><caption>Table 1: Depiction of how the label of $[i, j] = [i_u, j_u] \cap [i_v, j_v]$ is assigned as a function of whether $[i_u, j_u]$, $[i_v, j_v]$ are identical (row =), nested (rows &sub;, &sup;), or partially overlapping (row "else"), and as a function of the labels of $[i_u, j_u]$ (u columns) and of $[i_v, j_v]$ (v columns).</caption><thead><tr><th></th><th>u</th><th>v</th><th>u</th><th>v</th><th>u</th><th>v</th></tr><tr><th></th><th>0,C</th><th>0,C</th><th>0,C</th><th>P</th><th>P</th><th>P</th></tr></thead><tbody><tr><th>=</th><td>0</td><td></td><td>C</td><td></td><td>P</td><td></td></tr><tr><th>&sub;</th><td>C</td><td></td><td>C</td><td></td><td>P</td><td></td></tr><tr><th>&sup;</th><td>C</td><td></td><td>P</td><td></td><td>P</td><td></td></tr><tr><th>else</th><td>P</td><td></td><td>P</td><td></td><td>P</td><td></td></tr></tbody></table>

The rationale for these rules is to guarantee the following properties. If $[i, j]$ is a 0-gap at $w$, it means that $w$ is either a leaf with this gap, or if we let $u, v$
be the children of $w$, that there exists at least one leaf in $T_u$ and at least one leaf in $T_v$ such that $[i, j]$ is a gap at those leaves. If $[i, j]$ is a C-gap, it means that only one (but not the other) of $T_u$ and $T_v$ has some leaf with the same gap. Finally, if $[i, j]$ is a P-gap at $w$, it means that none of the leaves that descend from $w$ have this same gap. See Lemma 3 in S9 for a proof.

Labeled gaps also have implications about an optimal solution $A^+$ of the DPP, as shown by a few propositions stated in the next section: Proposition 1 shows that if a site $k$ does not belong to any labeled gap at $w$, then we must have $A_w^+[k] = 1$. Proposition 2 shows that if $[i, j]$ is a 0-gap at $w$, then $A_w^+[i, j]$ must be a gap; Proposition 3 shows that if $[i, j]$ is a P-gap at $w$, then $A_w^+[i, j] = A_p^+[i, j]$ where $p$ is $w$'s parent (hence the use of letter P for "parent"); finally Proposition 4 shows that if $[i, j]$ is a C-gap at $w$, then either $A_w^+[i, j]$ must be a gap, or $A_w^+[i, j] = A_p^+[i, j]$ (hence letter C for "choice").

<figure><img src="image_6.png" /><figcaption>Figure 3: A depiction of the labeled gaps constructed by the bottom-up phase on the input shown in Fig. 1. Note that Algorithm 2 does not store labeled gaps in this form.</figcaption></figure>

**Example for the bottom-up phase.** At the end of the bottom-up phase, the example shown in Fig. 1 produces the following labeled gaps at the internal nodes of the input tree.

Node $y$ has two C-gaps, [2, 3] and [5, 7], both obtained by application of rule R1.4 on two 0-gaps at $c$ and $d$ (the children of $y$) with one 0-gap strictly contained in the other. Node $y$ also has a 0-gap [10, 11] obtained by application of rule R1.1 on two identical 0-gaps at $c$ and $d$.

Node $x$ has two C-gaps, [2, 2] and [6, 6], both obtained by application of rule R1.4 on a 0-gap at $b$ strictly contained in a C-gap at $y$. Node $x$ also has a P-gap [10, 10] obtained by application of rule R1.6 on two partially overlapping 0-gaps at $b$ and $y$.

Finally, $r$ has no labeled gap, as no two labeled gaps at $a$ and $x$ overlap. The labeled gaps listed above are depicted in Fig. 3, although we remark that the algorithm does not need to store labeled gaps in table form.

**Top-down phase.** The goal of the second phase is to remove C-gaps and P-gaps and replace them with one or more 0-gaps. These 0-gaps, together with those already determined during the bottom-up phase, will constitute the actual gaps in the output. As we explain below, sometimes (namely when applying rule R2.2), multiple optimal choices are possible, in which case we duplicate the current solution and carry on.

Algorithm 3 provides a recursive implementation of this phase in full detail, where intermediate solutions (those still containing C-gaps and P-gaps) are represented using a data structure $S$ (e.g. a map) that allows to retrieve a set $S_w$ of labeled gaps for each node $w$. We visit every node of the tree in pre-order. It is easy to see that the root $r(T)$ cannot have any labeled gaps,
because otherwise the input alignment would have one or more columns only containing 0s. (See Corollary 2 in S9.) Thus, no action is needed for $r(T)$. For every other node $w$, we process each of its C-gaps and P-gaps using one of the rules R2.1 to R2.3 below, where $p$ denotes the parent of $w$.

**R2.1** If $[i, j]$ is a C-gap at $w$, and $[i, j]$ is a 0-gap at $p$, then $[i, j]$ is reset as a 0-gap at $w$.

**R2.2** If $[i, j]$ is a C-gap at $w$, and $[i, j]$ is not a 0-gap at $p$, we duplicate the current solution. In the original solution we reset $[i, j]$ as a 0-gap at $w$. In the duplicate, we copy at $w$ all 0-gaps that are in $[i, j]$ at $p$. That is, for each 0-gap $[i', j'] \subseteq [i, j]$ at $p$, we set $[i', j']$ as a 0-gap at $w$, and $[i, j]$ is removed from the set of labeled gaps at $w$.

**R2.3** If $[i, j]$ is a P-gap at $w$, then we copy at $w$ all 0-gaps that are in $[i, j]$ at $p$. That is, for each 0-gap $[i', j'] \subseteq [i, j]$ at $p$, we set $[i', j']$ as a 0-gap at $w$. Moreover, $[i, j]$ is removed from the set of labeled gaps at $w$.

As a result of the repeated application of these rules, in the end every stored solution is encoded as a collection of 0-gaps. To produce an alignment extension $A^+$ it suffices to set every site within a 0-gap to 0 and every other site to 1.

We remark that every time we copy at $w$ all 0-gaps that are in $[i, j]$ at $w$'s parent $p$, this has the effect of setting $A_w^+[i, j] = A_p^+[i, j]$. This is a consequence of the fact that no two labeled gaps $[i_w, j_w]$, $[i_p, j_p]$ at $w$ and $p$ respectively, can be partially overlapping (see Lemma 1 in S9).

<table><caption>**Algorithm 3: TopDown($w, Sol$)**</caption><thead><tr><th></th><th></th></tr></thead><tbody><tr><td>1</td><td>if w is not a leaf then</td></tr><tr><td>2</td><td>for S &isin; Sol do</td></tr><tr><td>3</td><td>&nbsp;&nbsp;&nbsp;&nbsp;Let S<sub>w</sub><sup>C</sup> be the subset of S<sub>w</sub> containing the C-gaps</td></tr><tr><td>4</td><td>&nbsp;&nbsp;&nbsp;&nbsp;Let S<sub>w</sub><sup>P</sup> be the subset of S<sub>w</sub> containing the P-gaps</td></tr><tr><td>5</td><td>&nbsp;&nbsp;&nbsp;&nbsp;for [i, j] &isin; S<sub>w</sub><sup>C</sup> &cup; S<sub>w</sub><sup>P</sup> do</td></tr><tr><td>6</td><td>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Apply one of rules R2.1 to R2.3 // moves [i, j] out of S<sub>w</sub><sup>C</sup> &cup; S<sub>w</sub><sup>P</sup></td></tr><tr><td>7</td><td>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;if R2.2 is applied then</td></tr><tr><td>8</td><td>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Let S' be the duplicated solution</td></tr><tr><td>9</td><td>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Sol &leftarrow; Sol &cup; {S'}</td></tr><tr><td></td><td>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;// S' added to the end of Sol; will be visited later on</td></tr><tr><td></td><td>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;// At this point, every S &isin; Sol only has 0-gaps in node w</td></tr><tr><td>10</td><td>&nbsp;&nbsp;&nbsp;&nbsp;Let u and v the children of w</td></tr><tr><td>11</td><td>&nbsp;&nbsp;&nbsp;&nbsp;TopDown(u, Sol)</td></tr><tr><td>12</td><td>&nbsp;&nbsp;&nbsp;&nbsp;TopDown(v, Sol)</td></tr></tbody></table>
**Example for the top-down phase.** The top-down phase starting from the labeled gaps shown in Fig. 3 produces 16 different optimal solutions. (Because there are four C-gaps and each C-gap can be processed in two possible ways by rule R2.2.) Fig. 4 shows how the top-down phase constructs one of these optimal solutions.

At $x$, both C-gaps are reset as 0-gaps by one of the two choices for rule R2.2 (steps 1,2 in Fig. 4), and the P-gap is set to 1 by rule R2.3, as $r$ has no 0-gap at the interval $[10, 10]$ (step 3). At $y$, we only apply rule R2.2: C-gap $[2, 3]$ is replaced by 0-gap $[2, 2]$, which is a 0-gap at $x$ (in practice this sets $A_y^+[2, 3]$ as a copy of $A_x^+[2, 3]$; see step 4 in Fig. 4); the C-gap $[5, 7]$ at $y$ is instead reset as a 0-gap (step 5).

The optimal solution produced in the end is depicted in its entirety in Fig. 5.

<table><thead><tr><th>step</th><th></th><th>0</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th><th>10</th><th>11</th><th>12</th></tr></thead><tbody><tr><td rowspan="3">0</td><td>r</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td></tr><tr><td>x</td><td>1</td><td>1</td><td>C</td><td>1</td><td>1</td><td>C</td><td>1</td><td>1</td><td>P</td><td>1</td><td>1</td><td>1</td><td>1</td></tr><tr><td>y</td><td>1</td><td>1</td><td>C</td><td>C</td><td>1</td><td>C</td><td>C</td><td>1</td><td>1</td><td>0</td><td>0</td><td>1</td><td>1</td></tr><tr><td rowspan="3">1</td><td>r</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td></tr><tr><td>x</td><td>1</td><td>1</td><td>0</td><td>1</td><td>1</td><td>C</td><td>1</td><td>1</td><td>1</td><td>P</td><td>1</td><td>1</td><td>1</td></tr><tr><td>y</td><td>1</td><td>1</td><td>C</td><td>C</td><td>1</td><td>C</td><td>C</td><td>1</td><td>1</td><td>0</td><td>0</td><td>1</td><td>1</td></tr><tr><td rowspan="3">2</td><td>r</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td></tr><tr><td>x</td><td>1</td><td>1</td><td>0</td><td>1</td><td>1</td><td>1</td><td>0</td><td>1</td><td>1</td><td>P</td><td>1</td><td>1</td><td>1</td></tr><tr><td>y</td><td>1</td><td>1</td><td>C</td><td>C</td><td>1</td><td>C</td><td>C</td><td>1</td><td>1</td><td>0</td><td>0</td><td>1</td><td>1</td></tr><tr><td rowspan="3">3</td><td>r</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td></tr><tr><td>x</td><td>1</td><td>1</td><td>0</td><td>1</td><td>1</td><td>1</td><td>0</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td></tr><tr><td>y</td><td>1</td><td>1</td><td>C</td><td>C</td><td>1</td><td>C</td><td>C</td><td>1</td><td>1</td><td>0</td><td>0</td><td>1</td><td>1</td></tr><tr><td rowspan="3">4</td><td>r</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td></tr><tr><td>x</td><td>1</td><td>1</td><td>0</td><td>1</td><td>1</td><td>1</td><td>0</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td></tr><tr><td>y</td><td>1</td><td>1</td><td>0</td><td>1</td><td>1</td><td>C</td><td>C</td><td>1</td><td>1</td><td>0</td><td>0</td><td>1</td><td>1</td></tr><tr><td rowspan="3">5</td><td>r</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td></tr><tr><td>x</td><td>1</td><td>1</td><td>0</td><td>1</td><td>1</td><td>1</td><td>0</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td></tr><tr><td>y</td><td>1</td><td>1</td><td>0</td><td>1</td><td>1</td><td>0</td><td>0</td><td>0</td><td>1</td><td>1</td><td>0</td><td>0</td><td>1</td></tr></tbody></table>

Figure 4: A sequence of intermediate solutions produced by the top-down phase on the input of Fig. 1, leading to an optimal solution (step 5). Each solution is represented as an alignment restricted to the internal nodes $r$, $x$ and $y$, even though the algorithm does not store labeled gaps in this form. Colored rectangles indicate sites covered or previously covered (lighter color) by labeled gaps. At step 0, no labeled gap has been processed and the intermediate solution coincides with that produced by the bottom-up phase (Fig. 3). Each of the subsequent steps processes exactly one C- or P-gap. At steps 1, 2, 4 and 5, rule R2.2 is employed and a different choice from the one depicted could have been made.

<figure><img src="image_7.png" /></figure>

<table><thead><tr><th></th><th>0</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th><th>10</th><th>11</th><th>12</th></tr></thead><tbody><tr><th>r</th><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td></tr><tr><th>a</th><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td></tr><tr><th>x</th><td>1</td><td>1</td><td>0</td><td>1</td><td>1</td><td>1</td><td>0</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td></tr><tr><th>b</th><td>1</td><td>1</td><td>0</td><td>1</td><td>1</td><td>1</td><td>0</td><td>1</td><td>1</td><td>0</td><td>0</td><td>1</td><td>1</td></tr><tr><th>y</th><td>1</td><td>1</td><td>0</td><td>1</td><td>1</td><td>0</td><td>0</td><td>0</td><td>1</td><td>0</td><td>0</td><td>1</td><td>1</td></tr><tr><th>c</th><td>1</td><td>1</td><td>0</td><td>0</td><td>1</td><td>0</td><td>0</td><td>0</td><td>1</td><td>1</td><td>0</td><td>0</td><td>1</td></tr><tr><th>d</th><td>1</td><td>0</td><td>0</td><td>0</td><td>1</td><td>0</td><td>0</td><td>0</td><td>0</td><td>1</td><td>0</td><td>0</td><td>1</td></tr></tbody></table>

Figure 5: An optimal solution $A^+$ for the example in Fig. 1, found following the steps detailed in Fig. 4. (Left:) Tree representation of $A^+$, where each of the 8 deletions in $A^+$ is depicted as a black box, as in Fig. 2. (Right:) Table representation of $A^+$. Colored rectangles indicate sites previously covered by labeled gaps.
# 4 Correctness of the algorithm

In order to prove that our algorithm finds all (and only) the optimal solutions of the DPP, we need the following key results.

**Proposition 1.** At the end of the bottom-up phase, site $k$ is not in a labeled gap at $w$ if and only if $A_w^+[k] = 1$ for any candidate solution $A^+$ of the DPP.

**Proposition 2.** If $[i, j]$ is set as a 0-gap at $w$ during the bottom-up phase then, in every optimal solution $A^+$ of the DPP, $A_w^+[i, j]$ is a gap.

**Proposition 3.** If $[i, j]$ is set as a P-gap at $w$ in the bottom-up phase then, in every optimal solution of the DPP, $A_w^+[i, j] = A_p^+[i, j]$ with $p$ being the parent of $w$.

**Proposition 4.** If $[i, j]$ is set as a C-gap at $w$ in the bottom-up phase then, in every optimal solution $A^+$ of the DPP, $A_w^+[i, j]$ is a gap or $A_w^+[i, j] = A_p^+[i, j]$ with $p$ being the parent of $w$.

The four propositions above, together with a few other auxiliary results, allow us to prove the correctness of the algorithm. All proofs can be found in supplement S9.

**Theorem 1.** $A^+$ is among the alignment extensions constructed by Algorithm 1 if and only if $A^+$ is an optimal solution of the DPP.

# 5 Graph-based representation of solutions

Since the top-down phase can generate a number of solutions exponential in the number of C-gaps, we now describe a way to derive a compact representation of the set of all DPP-optimal reconstructions for any given node. It consists of a collection of graphs, computable in polynomial time.

For each internal node $w \in V(T) \setminus L(T)$, we define a directed graph $G_w = (V_w, E_w)$, where $V_w = [0, m+1]$ is the same for every $w$, and coincides with the set of sites. As for the set of arcs $E_w$, it is defined by the following algorithm, which should be executed after the bottom-up phase, and replaces the top-down phase of the algorithm described in Sec. 3. We call this the *graph construction phase*.

The graph construction phase again visits the nodes of the tree with a pre-order traversal. Thus, for any non-root node $w$ with parent $p$, by the time $w$ is traversed $E_p$ has already been constructed. Moreover the set $S_w$ of labeled gaps at $w$ is also available.

Initially set $E_w$ as follows:

**R3.0** $E_w = \{(k, k+1) : k \in [0, m]\}$ and neither $k$ nor $k+1$ are in a labeled gap at $w$

Then, for each labeled gap $[i, j]$ at $w$:

**R3.1** If $[i, j]$ is a 0-gap, then add arc $(i-1, j+1)$ to $E_w$.

**R3.2** If $[i, j]$ is a P-gap, then for each arc $(h, k) \in E_p$ with $h, k \in [i-1, j+1]$, add $(h, k)$ to $E_w$.
**R3.3** If $[i, j]$ is a C-gap, then first add arc $(i-1, j+1)$ to $E_w$; next, for each arc $(h, k) \in E_p$ with $h, k \in [i-1, j+1]$, add $(h, k)$ to $E_w$.

In other words, if we only consider the set of arcs in $E_w$ that are between two vertices in $[i-1, j+1]$, which we denote with $E_w[[i-1, j+1]]$, we have:

$$E_w[[i-1, j+1]] = \begin{cases} \{(i-1, j+1)\} & \text{if } [i, j] \text{ is a 0-gap,} \\ E_p[[i-1, j+1]] & \text{if } [i, j] \text{ is a P-gap,} \\ E_p[[i-1, j+1]] \cup \{(i-1, j+1)\} & \text{if } [i, j] \text{ is a C-gap.} \end{cases}$$

Note that since the root $r(T)$ has no labeled gap (see Corollary 2 in S9), the algorithm sets $E_r(T) = \{(k, k+1) : k \in [0, m]\}$.

**Correctness of the graph construction phase.** Theorem 2 below states that each directed path in $G_w$ from the first to the last site corresponds to one possible way of setting $A_w^+$ optimally, and inversely every optimal $A_w^+$ corresponds to one such path. Thus $G_w$ provides a compact representation of all the optimal ancestral reconstructions for node $w$.

More formally, the graph $G_w$ verifies the following property, where we recall that $\mathbb{1}(A_w^+) = \{k \in [0, m+1] : A_w^+[k] = 1\}$. The proof can be found in supplement S10.

**Theorem 2.** Let $V \subseteq [0, m+1]$ and $w \in V(T) \setminus L(T)$. $V$ is the vertex set of a directed path from 0 to $m+1$ in $G_w$ if and only if $\exists A^+$ optimal solution of the DPP such that $V = \mathbb{1}(A_w^+)$.

**Example for the graph construction phase.** Suppose we executed the bottom-up phase on the input of Fig. 1, producing the labeled gaps in Fig. 3. The graph construction phase starts by constructing $G_r$, which (as noted above) is the trivial graph where two vertices are only adjacent if they are consecutive.

Next, $G_x = ([0, 12], E_x)$ is set to the graph shown in Fig. 6. In detail, rule R3.0 initially sets $E_x = \{(0, 1), (3, 4), (4, 5), (7, 8), (8, 9), (11, 12)\}$ (as sites 2, 6 and 10 are in labeled gaps). Then rule R3.3 is applied to C-gap [2, 2], which leads to adding first arc (1, 3) and then arcs (1, 2), (2, 3) (copied from $G_r$). Next, rule R3.3 is applied to C-gap [6, 6], which similarly leads to adding arcs (5, 7), (5, 6), (6, 7) to $E_x$. Finally, rule R3.2 is applied to P-gap [10, 10],

<figure><img src="image_8.png" /><figcaption>Figure 6: Graph representation of all possible $A_x^+$, where $A^+$ is constrained to be an optimal solution for the example in Fig. 1. (Top:) The graph $G_x$ obtained by the graph construction phase, with one path highlighted. (Bottom:) The four possible $A_x^+$ when $A^+$ is an optimal solution. By Theorem 2, there is a one-to-one correspondence between these $A_x^+$ and the paths from vertex 0 to vertex 12 in $G_x$. The path highlighted in $G_x$ corresponds to the highlighted $A_x^+$.</figcaption></figure>
which results in adding arcs (9, 10), (10, 11) (copied from $G_r$). This completes the construction of $G_x$.

Next, $G_y = ([0, 12], E_y)$ is set to the graph shown in Fig. 7. In detail, rule R3.0 initially sets $E_y = \{(0, 1), (8, 9)\}$, as the labeled gaps at $y$ are [2, 3], [5, 7], [10, 11]. Then rule R3.3 is applied to C-gap [2, 3], which leads to adding to $E_y$ arc (1, 4) and then all the arcs in $E_x[[1, 4]]$. Next, rule R3.3 is applied to C-gap [5, 7], which leads to adding to $E_y$ arc (4, 8) and then all the arcs in $E_x[[4, 8]]$. Finally, rule R3.1 is applied to 0-gap [10, 11], which results in adding arc (9, 12) to $E_y$. This completes the construction of $G_y$, and the graph construction phase terminates.

<figure><img src="image_9.png" /><figcaption>Figure 7: Graph representation of all possible $A_y^+$, where $A^+$ is constrained to be an optimal solution for the example in Fig. 1. (Top:) The graph $G_y$ obtained by the graph construction phase, with one path highlighted. (Bottom:) The nine possible $A_y^+$ when $A^+$ is an optimal solution. The path highlighted in $G_y$ corresponds to the highlighted $A_y^+$.</figcaption></figure>

To illustrate Theorem 2, we focus on $G_y$ (Fig. 7). Theorem 2 implies that the set of paths in $G_y$ from 0 to 12 is in a one-to-one correspondence with the set of binary strings that can be obtained as $A_y^+$ for some optimal solution $A^+$.

To verify this, we can imagine deriving all 16 optimal solutions $A^+$ and then taking the row $A_y^+$ in each of them. Doing so produces 9 different values for $A_y^+$, the binary rows depicted in Fig. 7. (Sometimes pairs of optimal solutions that have different $A_x^+$ may have the same $A_y^+$, which explains why from 16 optimal solutions we only get 9 different $A_y^+$.) It is easy to check that each of these 9 values for $A_y^+$, when translated as $\mathbb{1}(A_y^+)$, gives the vertex set of a path from 0 to 12 in $G_y$. Conversely, each path in $G_y$, when converted into the binary string that has a 1 at site $k$ if the path traverses $k$, corresponds to one of the 9 rows shown in Fig. 7.
# 6 Complexity

Here we state the main results about the complexity of our algorithms for the DPP. Full proofs can be found in supplement S11. The asymptotic complexity will be expressed as a function of the following parameters:

* $n$, the number of leaves of $T$ or equivalently the number of rows of $A$;

* $s$, the number of optimal solutions of the DPP;

* $b$, the number of boundaries in $A$, where a *boundary* is defined as a pair of consecutive sites $(k, k+1)$, such that $A[k] \neq A[k+1]$, i.e., their corresponding columns differ;

* $m$, the number of columns in $A$ (excluding the first and last column which we only added for convenience of notation).

For the analysis of Algorithm 1, we assume that both the input alignment $A$ and the output alignment $A^+$ are encoded in *gap form*, instead of the table form presented in the Preliminaries (sec. 2.1). The gap form of an alignment $A$ consists of encoding each row of $A$ as the set of gaps it contains. For example, the input alignment of Fig. 1 can be encoded as $A_a = \emptyset, A_b = \{[2, 2], [6, 6], [9, 10]\}, A_c = \{[2, 3], [5, 7], [10, 11]\}, A_d = \{[1, 3], [5, 8], [10, 11]\}$.

We assume that alignments are represented in gap form because, as a result, running times and memory requirements only depend on parameter $b$ rather than $m$. In practice, for many alignments, we can expect $b$ to be much smaller than $m$. We can now state the first of the two main results of this section.

**Theorem 3.** Assume that the input alignment $A$ and every solution $A^+$ in output is represented in gap form. Then Algorithm 1 runs in $O(nbs)$ time.

Note that an obvious bound is $s \le 2^c$, where $c$ denotes the number of C-gaps at the end of the bottom-up phase. Although the worst-case complexity of Algorithm 1 is exponential in $c$, we can avoid the exponential factor in two ways. First, we can simply keep only one solution, which can be achieved by removing the duplication in rule R2.2 and skipping lines 7–9 in Algorithm 3. The complexity then becomes $O ngu)$.

Alternatively the top-down phase can be replaced by the graph construction phase described in Sec. 5. We now analyze the complexity of this algorithm. A key observation here is given by the following proposition.

**Proposition 5.** For any internal node $w$, $G_w$ is a planar graph.

Proposition 5 can be used to obtain the other main complexity result:

**Theorem 4.** The algorithm that first runs the bottom-up phase and then the graph construction phase runs in $O(nm)$ time.

Note that the complexity in Theorem 4 depends on $m$ rather than on $b$. This is only because, for simplicity, we have chosen to describe the graph construction phase so that the graphs have $[0, m+1]$ as vertex set. It is possible to see that it is feasible to modify the graph construction phase so that each graph has a vertex set of size $O(b)$. Doing so would allow us to obtain an algorithm running in $O ngu)$ time. We leave this as future work.
# 7 Practical relevance of the deletion-only case

Although this may not be evident at first sight, the relevance of deletion-only case extends beyond the solution of the DPP. In fact, algorithms to solve the DPP can be used to solve (parts of) many instances of the general problem (the IPP; Definition 5). This is due to three observations, which first we state informally and then illustrate with an example below.

1. Many instances of the IPP can be solved with a divide-and-conquer approach, by combining solutions for subproblems of the original problem. A way to do this was shown by Chindelevitch et al. [27], and another one is implicit in the work by Fredslund, Hein, and Scharling [18].

2. The IPP is *time-reversible*, in the sense that if the input tree is rerooted on a node different from the original root, then every candidate solution retains its cost. As a consequence, optimal solutions of the IPP (excluding the reconstruction for the root) are independent of the position of the root.

3. Recall that phylogenetic correctness may constrain an internal node *x* to have many of its sites equal to 1. Now suppose that there exists a node *x* such that $A_x^+$ must only consist of 1s, for any candidate solution $A^+$ of the IPP. Then if we re-root the tree in *x*, it is easy to see that none of the indels in any candidate solution $A^+$ can be an insertion. But then we can solve the IPP by solving the DPP on the re-rooted tree with one of our algorithms.

<figure><img src="image_10.png" /><figcaption>Figure 8: An instance $(T, A)$ of the IPP that can be solved by applying the algorithms for the DPP.</figcaption></figure>

The following example shows how these three ideas can be combined. Consider instance $(T, A)$ of the IPP in Fig. 8, where for simplicity we do not show the dummy columns $A[0]$ and $A[5]$. While at first sight there is no reason why the deletion-only case would be relevant to solve this example, it turns out that we can indeed use our algorithms for the DPP to solve it.

First, this is a case where we can decompose the problem in two subproblems, i.e., an example of point 1 above. Namely, we can split the input alignment in two sub-alignments $A_1$ and $A_2$, where $A_1$ consists of columns 1 and 2 of the original alignment, while $A_2$ contains columns 3 and 4. The solutions for the instance in Fig. 8 can be obtained by combining the solutions for the IPP on $(T, A_1)$ and those for the IPP on $(T, A_2)$. This can be proven by applying previous results [27, 18], but the main idea here is that, as can be easily verified, every indel must be within one of the two sub-alignments $A_1$ or $A_2$.

Because $(T, A_1)$ and $(T, A_2)$ are symmetric, we can just concentrate on the instance $(T, A_1)$ of the IPP. Here phylogenetic correctness implies that every
sequence on the path between *a* and *b* must have 1s at both sites 1 and 2. Thus, if we create a new node *r'* subdividing edge (*x*, *a*), we know that we must have $A_{1r'}^{+}[1, 2] = 11$. We then re-root the tree in *r'*, which gives us a new rooted tree $T'$ shown in Fig. 9. For the reasons explained in points 2 and 3 above, solving the IPP on (*T*, $A_1$) is equivalent to solving the DPP on ($T'$, $A_1$). (Excluding the ancestral reconstruction for the root, which we discuss below.)

<figure><img src="image_11.png" /><figcaption>Figure 9: An instance ($T'$, $A_1$) of the DPP that provides a solution of the first half of the IPP problem in Fig. 8. This instance has two optimal solutions with $c(A_1^+) = 2$. Both solutions have $A_{1r'}^+ = A_{1x}^+ = 11$, but they differ at node *y* where they have $A_{1y}^+ \in \{11, 01\}$.</figcaption></figure>

It is now easy to check that, for the instance in Fig. 9, our algorithm gives the two solutions with 11 at node *x* and with one of {11, 01} at node *y*. By symmetry the IPP on (*T*, $A_2$) has the two solutions with one of {11, 10} at node *x* and with 11 at node *y*. Combining each solution for one sub-problem with each solution for the other sub-problem gives four solutions for the original IPP on (*T*, *A*) (Fig. 8). They are the ones with $A_x^+ \in \{1111, 1110\}$ and $A_y^+ \in \{1111, 0111\}$. As for the root, it is easy to see that setting it equal to one of the reconstructions for its children always leads to an optimal solution. Thus, here we can set $A_r^+ = A_x^+$ or $A_r^+ = A_y^+$.

**Applicability of the DPP to a real-world dataset.** To investigate whether the strategy illustrated above can be applied frequently, we empirically evaluated it on OrthoMaM, a large collection of curated alignments of orthologous coding sequences [42]. Specifically, we downloaded OrthoMaM v12a [43], which includes 15,868 alignments from 190 complete mammalian genomes. Each alignment is available in four versions (nucleotide / amino acid, filtered / unfiltered) and is accompanied by a maximum-likelihood tree inferred from it. The filtered alignments are obtained by applying state-of-the-art alignment cleaning methods (implemented in MACSE, HMMcleaner and PhylteR; see [43] for details). For each alignment, both deletions and insertions are likely responsible for the gaps present in it.

We implemented the construction of the graph data structure introduced by [18]. In short, a consequence of this graph is that the set of gap characters in the input alignment (i.e. the set of pairs (*x*, *k*) for which $A_x[k] = 0$) is partitioned into “components” so that the resolution of the IPP can be performed independently for each of the components. For each of the components, our code then tests if the tree can be re-rooted so that resolution of the IPP is equivalent to that of the DPP, as described in the example above. If this is the case, the component is classified as solvable, and not solvable otherwise. To facilitate interpretation of these results, we then classify alignment sites, rather than components, as not solvable, partially solvable or entirely solvable: a site *k* with gaps
is *not solvable* if all of the gap characters at that site belong to components that are not solvable. A site *k* is *partially solvable* if it contains some gap characters belonging to a solvable component and some other gap characters belonging to a component that is not solvable. A site is *entirely solvable* if all of its gap characters belong to solvable components. The code to produce these statistics for any input alignment is available at https://gite.lirmm.fr/rivals/asr_indel.

<table><caption>Table 2: Experiment on the two OrthoMaM amino acid alignments datasets. For each dataset, we show the total number of sites with gaps (across 15,868 alignments), and the percentages of these sites that are entirely, partially, or not solvable by reduction to the DPP.</caption><thead><tr><th>Dataset</th><th>gapped sites</th><th>entirely s.</th><th>partially s.</th><th>not solvable</th></tr></thead><tbody><tr><td>filtered</td><td>8,439,079</td><td>26.80%</td><td>23.00%</td><td>50.20%</td></tr><tr><td>unfiltered</td><td>12,070,178</td><td>15.07%</td><td>15.78%</td><td>69.15%</td></tr></tbody></table>

Table 2 shows the results of the site classification procedure described above on the totality of the amino acid alignments, for the filtered and unfiltered OrthoMaM datasets. As shown under “gapped sites”, the filtering step results in the removal of about 30% of the sites with gaps across the entirety of the 15,868 alignments. In the filtered dataset, out of the sites with gaps about a quarter are entirely solvable, a quarter are partially solvable and half are not solvable at all. The proportion of sites with gaps that can be solved, partially or entirely, drops to about 30% for the unfiltered dataset.

Aside from showing that the DPP can be applied relatively frequently, the fact that the percentages of solvability drop for the unfiltered alignments shows that the parts of the alignments that are filtered out tend to pose more problems to our approach. This is not unexpected, as alignment cleaning methods tend to filter out the low-quality portions of multiple alignments, which are also the ones that are most likely to violate the assumptions necessary for the reduction to the DPP (presence of an internal node that can be used as a root, phylogenetic correctness, etc.).

# 8 Discussion

**Novel contributions and related works.** We have focused on a version of the Indel Parsimony Problem where only deletions are allowed (the DPP; Definition 6). We have described an exact algorithm that constructs all the optimal solutions for the DPP, and proved it is correct. While the algorithm's running time may be exponential in the size of the input, this is only due to the exponential number of solutions. The algorithm can be easily changed to only find a single optimal solution in polynomial time (e.g. skip lines 7 – 9 in Algorithm 3). We note that an algorithm to find an optimal solution of the DPP was also proposed by Chindelevitch et al. [27]. Unfortunately, the pseudo-code provided in that paper appears to be incorrect, although the main underlying ideas are sound. See supplement S12, where we copied the pseudo-code and provide a simple example where it fails to return an optimal solution. Importantly, the main difference with our work is that we seek the whole set of optimal solutions, a goal which significantly complicates the task.
Two recently published parsimony-based works also deserve some discussion [25, 26]. Iglhaut et al. [25] proposed an extension of Fitch's algorithm [9] that handles insertions and deletions. The proposed algorithm is not intended to solve any explicit formulation of the parsimony principle for indels. Importantly, like in the original algorithm by Fitch, ancestral reconstruction is executed one site at a time, independently. Second, Tule et al. [26] present a Mixed-Integer Programming (MIP) formulation of indel inference, which provides an exact solution, albeit not polynomial in the worst case. However, their formulation is radically different from that of the Indel Parsimony Problem [27]. For example, phylogenetic correctness is not imposed, and the number of indel events between the reconstructions for two adjacent nodes can be very different (e.g. Tule et al. [26] count three separate indel events between 101010101 and 100000001 whereas Def. 4 only one).

Back to our paper, since there is no reason to favor one optimal scenario over an equally-parsimonious one, we consider it important to conceive ways to represent compactly the diversity of the reconstructed solutions. Our other main contribution is that we have shown that, in the context of the DPP, it is possible to efficiently construct a graph $G_w$ representing the set of optimal reconstructions $A_w^+$ for any internal node $w$ of the phylogeny. While the software GRASP does indeed construct similar graphs [23], we note that there is no proven global optimality guarantee for the reconstructions represented by the graphs of GRASP. We believe that our result (i.e. Theorem 2) is the first to provide a mathematical justification for the use of partial order graphs [40] in the context of ancestral sequence reconstruction. More recently, Tule et al. [26] also infer a graph for each internal node, but each of these graphs only has a single path from the first to the last site, meaning that it cannot be used to represent multiple reconstructions.

Naturally, we have derived the asymptotic complexity for all our algorithms. These algorithms are asymptotically optimal, because their running time is linear in the size of their inputs, or in the size on memory of its intended output (the set of all optimal solutions, or the graphs representing them).

Finally, we have shown that restricting the model to deletions is not as significant a limitation as it may initially seem. In an experiment on 15,868 alignments from a real-world dataset, we found that about 30% to 50% of the sites (for unfiltered and filtered alignments, respectively) can be at least partially solved by using our algorithms, even though the optimal solutions involve both deletions and insertions.

**Beyond unit costs and binary trees.** In both the IPP and the DPP (Definitions 5 and 6) each indel has a constant unit cost. However, it may be biologically relevant to consider more general indel cost functions. For example, it may be interesting to penalize longer indels by setting their cost with the affine function $\alpha+\beta\ell$, where $\ell$ is the length of the indel [18, 27, 26]. Additionally the cost of insertions may be different from that of deletions (in the IPP), or the cost of an indel may depend on the edge where it occurs, allowing to model edge lengths [27]. These observations naturally raise the question of whether our algorithms can be extended to solve these more general formulations. While we leave the treatment of affine and asymmetric costs as open questions, it is interesting to note that the DPP with edge-dependent costs can also be solved
using an extension of our techniques.

In supplement S13, we describe how our algorithms can be extended to solve the DPP where for each edge $(u, v) \in E(T)$ we specify a cost $c_{uv}$ for the deletions occurring from $u$ to $v$. The only thing that needs to be changed is the way labels 0, P and C are assigned to labeled gaps during the bottom-up phase. All other aspects of our algorithms remain unchanged, including the fact that we can represent all optimal solutions using the graph construction phase of Section 5. To define the label of each labeled gap, we compute by dynamic programming two numerical scores describing the consequences of taking a deletion at that gap, and those of not taking it. Depending on the difference between these two scores, it may always be convenient to take a deletion at that labeled gap—in which case we assign label 0— never convenient to take a deletion —implying label P— or equally convenient —implying label C. See supplement S13 for full details, including a pseudocode that replaces Algorithm 2.

We remark that the algorithm of supplement S13 can also be used to solve the DPP on a non-binary tree $T$: to that end it suffices to take a binary resolution $T'$ of $T$ and assign an infinite (or sufficiently large) cost to each edge of $T'$ with no corresponding edge in $T$.

**Open problems.** The computational complexity result given by Chindelevitch et al. [27] implied the NP-hardness of two versions of the IPP:

* the IPP on trees with unbounded degree,

* the IPP on binary trees with branch-dependent indel costs.

This notably leaves open whether it is NP-hard to solve the IPP on binary trees and unit costs. In fact, we have chosen to tackle the DPP for binary trees and unit costs because of its potential impact on this open problem. While it may be interesting to clarify whether it is possible to extend our algorithms to more general indel costs—as discussed above—we believe that the most important avenue is to investigate whether our approach can be leveraged to also solve the general IPP on binary trees and unit costs. Should this prove impossible, it may be interesting to investigate more extensively its limitations on real-world data—for example with datasets covering a wider phylogenetic range than OrthoMaM—and to provide heuristics for the parts of the input alignments that cannot be solved exactly with our algorithms.

## References

[1] Hendrikse, N. M., Holmberg Larsson, A., Svensson Gelius, S., Kuprin, S., Nordling, E., and Syrén, P.-O. “Exploring the therapeutic potential of modern and ancestral phenylalanine/tyrosine ammonia-lyases as supplementary treatment of hereditary tyrosinemia”. In: *Scientific Reports* 10.1, 1315 (2020). QC 20200220.

[2] Ladics, G. S., Han, K.-H., Jang, M. S., Park, H., Marshall, V., Dersjant-Li, Y., and Sewalt, V. J. “Safety evaluation of a novel variant of consensus bacterial phytase”. In: *Toxicology Reports* 7 (2020), pp. 844–851.
[3] Ducatez, M. F., Bahl, J., Griffin, Y., Stigger-Rosser, E., Franks, J., Barman, S., Vijaykrishna, D., Webb, A., Guan, Y., Webster, R. G., et al. "Feasibility of reconstructed ancestral H5N1 influenza viruses for cross-clade protective vaccine development". In: *Proceedings of the National Academy of Sciences* 108.1 (2011), pp. 349–354.

[4] Chang, B. S. W., Jönsson, K., Kazmi, M. A., Donoghue, M. J., and Sakmar, T. P. "Recreating a Functional Ancestral Archosaur Visual Pigment". In: *Molecular Biology and Evolution* 19.9 (Sept. 2002), pp. 1483–1489.

[5] Blanchette, M., Green, E. D., Miller, W., and Haussler, D. "Reconstructing large regions of an ancestral mammalian genome in silico". In: *Genome Research* 14.12 (2004), pp. 2412–2423.

[6] Redelings, B. D., Holmes, I., Lunter, G., Pupko, T., and Anisimova, M. "Insertions and Deletions: Computational Methods, Evolutionary Dynamics, and Biological Applications". In: *Molecular Biology and Evolution* 41.9 (Aug. 2024), msae177.

[7] Linard, B., Swenson, K., and Pardi, F. "Rapid alignment-free phylogenetic identification of metagenomic sequences". In: *Bioinformatics* 35.18 (Jan. 2019), pp. 3303–3312.

[8] Romashchenko, N., Linard, B., Pardi, F., and Rivals, E. "EPIK: precise and scalable evolutionary placement with informative k-mers". In: *Bioinformatics* 39.12 (2023), btad692.

[9] Fitch, W. M. "Toward Defining the Course of Evolution: Minimum Change for a Specific Tree Topology". In: *Systematic Biology* 20.4 (Dec. 1971), pp. 406–416.

[10] Sankoff, D. "Minimal Mutation Trees of Sequences". In: *SIAM Journal on Applied Mathematics* 28.1 (1975), pp. 35–42.

[11] Yang, Z., Kumar, S., and Nei, M. "A new method of inference of ancestral nucleotide and amino acid sequences." In: *Genetics* 141.4 (1995), pp. 1641–1650.

[12] Koshi, J. M. and Goldstein, R. A. "Probabilistic reconstruction of ancestral protein sequences". In: *Journal of Molecular Evolution* 42 (1996), pp. 313–320.

[13] Schlüter, D., Price, T., Mooers, A. Ø., and Ludwig, D. "Likelihood of ancestor states in adaptive radiation". In: *Evolution* 51.6 (1997), pp. 1699–1711.

[14] Pupko, T., Pe, I., Shamir, R., and Graur, D. "A fast algorithm for joint reconstruction of ancestral amino acid sequences". In: *Molecular Biology and Evolution* 17.6 (2000), pp. 890–896.

[15] Snir, S. and Pachter, L. "Phylogenetic profiling of insertions and deletions in vertebrate genomes". In: *Annual International Conference on Research in Computational Molecular Biology*. Springer. 2006, pp. 265–280.

[16] Diallo, A. B., Makarenkov, V., and Blanchette, M. "Exact and heuristic algorithms for the indel maximum likelihood problem". In: *Journal of Computational Biology* 14.4 (2007), pp. 446–461.
[17] Snir, S. and Pachter, L. "Tracing the most parsimonious indel history". In: *Journal of Computational Biology* 18.8 (2011), pp. 967–986.

[18] Fredslund, J., Hein, J., and Scharling, T. "A large version of the small parsimony problem". In: *Algorithms in Bioinformatics: Third International Workshop, WABI 2003, Budapest, Hungary, September 15-20, 2003. Proceedings* 3. Springer. 2003, pp. 417–432.

[19] Kim, J. and Sinha, S. "Indelign: a probabilistic framework for annotation of insertions and deletions in a multiple alignment". In: *Bioinformatics* 23.3 (Nov. 2006), pp. 289–297.

[20] Blanchette, M., Diallo, A. B., Green, E. D., Miller, W., and Haussler, D. "Computational reconstruction of ancestral DNA sequences". In: *Phylogenomics* (2008), pp. 171–184.

[21] Simmons, M. P. and Ochoterena, H. "Gaps as characters in sequence-based phylogenetic analyses". In: *Systematic Biology* 49.2 (2000), pp. 369–381.

[22] Ashkenazy, H., Penn, O., Doron-Faigenboim, A., Cohen, O., Cannarozzi, G., Zomer, O., and Pupko, T. "FastML: a web server for probabilistic reconstruction of ancestral sequences". In: *Nucleic Acids Research* 40.W1 (May 2012), W580–W584.

[23] Foley, G., Mora, A., Ross, C. M., Bottoms, S., Sützl, L., Lamprecht, M. L., Zaugg, J., Essebier, A., Balderson, B., Newell, R., et al. "Engineering indel and substitution variants of diverse and ancient enzymes using Graphical Representation of Ancestral Sequence Predictions (GRASP)". In: *PLOS Computational Biology* 18.10 (2022), e1010633.

[24] Jowkar, G., Pečerska, J., Maiolo, M., Gil, M., and Anisimova, M. "ARPIP: Ancestral Sequence Reconstruction with Insertions and Deletions under the Poisson Indel Process". In: *Systematic Biology* 72.2 (July 2023), pp. 307–318.

[25] Iglhaut, C., Pečerska, J., Gil, M., and Anisimova, M. "Please Mind the Gap: Indel-Aware Parsimony for Fast and Accurate Ancestral Sequence Reconstruction and Multiple Sequence Alignment Including Long Indels". In: *Molecular Biology and Evolution* 41.7 (June 2024), msae109.

[26] Tule, S., Foley, G., Zhao, C., Forbes, M., and Bodén, M. "Optimal phylogenetic reconstruction of insertion and deletion events". In: *Bioinformatics* 40.Supplement 1 (June 2024), pp. i277–i286.

[27] Chindelevitch, L., Li, Z., Blais, E., and Blanchette, M. "On the inference of parsimonious indel evolutionary scenarios". In: *Journal of Bioinformatics and Computational Biology* 04.03 (2006), pp. 721–744.

[28] Diallo, A. B., Makarenkov, V., and Blanchette, M. "Finding maximum likelihood indel scenarios". In: *Comparative Genomics: RECOMB 2006 International Workshop, RCG 2006 Montreal, Canada, September 24-26, 2006 Proceedings* 4. Springer. 2006, pp. 171–185.

[29] Diallo, A. B. "Inference of Insertion and Deletion Scenarios for Ancestral Genome Reconstruction and Phylogenetic Analyses: Algorithms and Biological Applications". PhD thesis. 2009.
[30] Thorne, J. L., Kishino, H., and Felsenstein, J. "An evolutionary model for maximum likelihood alignment of DNA sequences". In: *Journal of Molecular Evolution* 33 (1991), pp. 114–124.

[31] Thorne, J. L., Kishino, H., and Felsenstein, J. "Inching toward reality: an improved likelihood model of sequence evolution". In: *Journal of Molecular Evolution* 34 (1992), pp. 3–16.

[32] Miklós, I., Lunter, G. A., and Holmes, I. "A "long indel" model for evolutionary sequence alignment". In: *Molecular Biology and Evolution* 21.3 (2004), pp. 529–540.

[33] Bouchard-Côté, A. and Jordan, M. I. "Evolutionary inference via the Poisson indel process". In: *Proceedings of the National Academy of Sciences* 110.4 (2013), pp. 1160–1166.

[34] De Maio, N. "The cumulative indel model: fast and accurate statistical evolutionary alignment". In: *Systematic Biology* 70.2 (2021), pp. 236–257.

[35] Holmes, I. "Historian: accurate reconstruction of ancestral sequences and evolutionary rates". In: *Bioinformatics* 33.8 (Jan. 2017), pp. 1227–1229.

[36] Westesson, O., Barquist, L., and Holmes, I. "HandAlign: Bayesian multiple sequence alignment, phylogeny and ancestral reconstruction". In: *Bioinformatics* 28.8 (Jan. 2012), pp. 1170–1171.

[37] Zhang, C., Huelsenbeck, J. P., and Ronquist, F. "Using parsimony-guided tree proposals to accelerate convergence in Bayesian phylogenetic inference". In: *Systematic Biology* 69.5 (2020), pp. 1016–1032.

[38] Kramer, A. M., Thornlow, B., Ye, C., De Maio, N., McBroome, J., Hinrichs, A. S., Lanfear, R., Turakhia, Y., and Corbett-Detig, R. "Online Phylogenetics with matOptimize Produces Equivalent Trees and is Dynamically More Efficient for Large SARS-CoV-2 Phylogenies than de novo and Maximum-Likelihood Implementations". In: *Systematic Biology* 72.5 (May 2023), pp. 1039–1051.

[39] Felsenstein, J. "Evolutionary trees from DNA sequences: a maximum likelihood approach". In: *Journal of Molecular Evolution* 17 (1981), pp. 368–376.

[40] Lee, C., Grasso, C., and Sharlow, M. F. "Multiple sequence alignment using partial order graphs". In: *Bioinformatics* 18.3 (Mar. 2002), pp. 452–464.

[41] Farris, J. S. "Phylogenetic analysis under Dollo's Law". In: *Systematic Biology* 26.1 (1977), pp. 77–88.

[42] Ranwez, V., Delsuc, F., Ranwez, S., Belkhir, K., Tilak, M.-K., and Douzery, E. J. "OrthoMaM: a database of orthologous genomic markers for placental mammal phylogenetics". In: *BMC evolutionary biology* 7 (2007), pp. 1–12.

[43] Allio, R., Delsuc, F., Belkhir, K., Douzery, E. J., Ranwez, V., and Scornavacca, C. "OrthoMaM v12: A database of curated single-copy ortholog alignments and trees to study mammalian evolutionary genomics". In: *Nucleic Acids Research* 52.D1 (2024), pp. D529–D535.
Supplementary Information to
Algorithms to reconstruct past indels:
the deletion-only parsimony problem.

This document contains all results not fully detailed in the main text, including
proofs of all mathematical statements.

S9 Correctness of the basic algorithm

All the propositions that follow are relative to a given instance (T, A) of the
Deletion-only Parsimony problem (DPP).

**Lemma 1.** Let $w$ be an internal node. If $[i, j]$ is a labeled gap at $w$ at the end of the bottom-up phase, then all nodes descending from $w$ must have a labeled gap that is a superset of $[i, j]$.

*Proof.* By induction. If $w$ is a leaf, then it only has itself as descendant and the statement is trivially true. Otherwise, let $u$ and $v$ be the children of $w$. The bottom-up phase must have obtained $[i, j]$ as $[i_u, j_u] \cap [i_v, j_v]$ where $[i_u, j_u]$ is a labeled gap at $u$ and $[i_v, j_v]$ is a labeled gap at $v$. By induction hypothesis, all nodes descending from $u$ must have a labeled gap that is a superset of $[i_u, j_u] \supseteq [i, j]$, and similarly all nodes descending from $v$ must have a labeled gap that is a superset of $[i_v, j_v] \supseteq [i, j]$. Thus the statement holds for $w$. $\square$

**Corollary 1.** Let $w$ be an internal node. If $[i, j]$ is a labeled gap at $w$ at the end of the bottom-up phase, then all leaves descending from $w$ must have a 0-gap that is a superset of $[i, j]$.

*Proof.* Trivial consequence of Lemma 1 and of the fact that every labeled gap at a leaf is a 0-gap. $\square$

**Corollary 2.** The root $r(T)$ does not have any labeled gap.

*Proof.* If $[i, j]$ was a labeled gap at $r(T)$, then by Corollary 1 every leaf of $T$ would have a 0-gap that is a superset of $[i, j]$. But then every site $k \in [i, j]$ would be such that $A_w[k] = 0, \forall w \in L(T)$, which violates the assumption that there are no columns only consisting of 0s (Sec. 2.1). $\square$

**Lemma 2.** At the end of the bottom-up phase, site $k$ is in a labeled gap at $w$ if and only if $A_x[k] = 0, \forall x \in L(T_w)$.

*Proof.* If a site $k$ is in a labeled gap at $w$, then $A_x[k] = 0, \forall x \in L(T_w)$ holds because of Corollary 1.

The opposite direction can be proven by induction on the height of $w$. If $w$ is a leaf, $A_x[k] = 0, \forall x \in L(T_w)$ implies $A_w[k] = 0$, which in turn implies that site $k$ is in a 0-gap at $w$. Now suppose that $w$ has children $u$ and $v$, and that $A_x[k] = 0, \forall x \in L(T_w)$. This implies $A_x[k] = 0, \forall x \in L(T_u)$ and $A_x[k] = 0, \forall x \in L(T_v)$. Then, by induction hypothesis, site $k$ must be both in a labeled gap $[i_u, j_u]$ at $u$, and in a labeled gap $[i_v, j_v]$ at $v$. But then site $k$ must also be in their intersection $[i, j] = [i_u, j_u] \cap [i_v, j_v]$, which is a labeled gap at $w$. $\square$
**Corollary 3.** Let $p$ be the parent node of $w$. At the end of the bottom-up phase, if a site $k$ is not in a labeled gap at $w$, then $k$ is also not in a labeled gap at $p$.

*Proof.* By Lemma 2, if $k$ is not in a labeled gap at $w$, then $\exists x \in L(T_w) \subseteq L(T_p)$ such that $A_x[k] = 1$, meaning that $k$ is not in a labeled gap at $p$. $\square$

**Lemma 3.** Let $w$ be a node. If $w$ is a leaf, and $[i, j]$ is a labeled gap at $w$, then $[i, j]$ is a 0-gap at $w$. If $w$ is an internal node, let $u$ and $v$ be its two children. At the end of the bottom-up phase, the following holds:

1. If $[i, j]$ is a 0-gap at $w$ then there exist $x_u$ leaf of $T_u$ and $x_v$ leaf of $T_v$, such that $[i, j]$ is a 0-gap at $x_u$ and $x_v$.

2. If $[i, j]$ is a C-gap at $w$ then exactly one of $T_u$ and $T_v$ has a leaf $x$ such that $[i, j]$ is a 0-gap at $x$.

3. If $[i, j]$ is a P-gap at $w$ then neither $T_u$ nor $T_v$ has a leaf $x$ such that $[i, j]$ is a 0-gap at $x$.

*Proof.* If $w$ is a leaf, the statement is trivially true. Now, suppose $w$ is an internal node, and that the lemma is verified for $u$ and $v$ (induction hypothesis). Let $[i_u, j_u]$ and $[i_v, j_v]$ be the labeled gaps in $u$ and $v$ respectively, such that $[i, j] = [i_u, j_u] \cap [i_v, j_v]$.

If $[i, j]$ is a 0-gap at $w$ then rule R1.1 must have been applied. But then $[i, j] = [i_u, j_u] = [i_v, j_v]$ and $[i, j]$ is not a P-gap at either $u$ or $v$. The induction hypothesis applied to $u$ implies that because $[i, j]$ is either a 0-gap or a C-gap at $u$, there exists at least one leaf in $T_u$ with $[i, j]$ as a 0-gap. The same argument holds for $v$. Therefore, $w$ has at least one leaf in each of its subtrees with $[i, j]$ as a 0-gap. This proves point 1 in the statement.

Second, if $[i, j]$ is a C-gap at $w$, then it has to be obtained via one of rules R1.2 or R1.4. If R1.2 is applied, then $[i, j] = [i_u, j_u] = [i_v, j_v]$ and without loss of generality, suppose that $[i, j]$ is a P-gap at $u$ but not at $v$. The induction hypothesis applied to $u$ and $v$ implies that there is no leaf $x_u$ in $T_u$ with a 0-gap $[i, j]$, but there is at least one leaf $x_v$ in $T_v$ with a 0-gap $[i, j]$. Thus, $w$ has such a leaf in exactly one of $T_u, T_v$. If R1.4 is applied, then without loss of generality, $[i, j] = [i_u, j_u] \subsetneq [i_v, j_v]$ and $[i, j]$ is not a P-gap at $u$. The induction hypothesis applied to $u$ implies that there is a leaf $x_u$ in $T_u$ with $[i, j]$ as a 0-gap. Moreover, by Corollary 1, we know that any leaf $x_v$ in $T_v$ has $[i', j'] \supseteq [i_v, j_v] \supsetneq [i, j]$ as a 0-gap, meaning that $[i, j]$ cannot be a 0-gap at $x_v$. Thus, exactly one of $T_u, T_v$ has a leaf having $[i, j]$ as a 0-gap. This proves point 2 in the statement.

Finally, if $[i, j]$ is a P-gap at $w$, then it can only have been obtained via one of rules R1.3, R1.5, or R1.6. If R1.3 is applied, then $[i, j] = [i_u, j_u] = [i_v, j_v]$ and $[i, j]$ is a P-gap at $u$ and $v$. The induction hypothesis implies that no leaf $x$ in $T_u$ or $T_v$ has $[i, j]$ as 0-gap. If R1.5 is applied, then, without loss of generality, $[i, j] = [i_u, j_u] \subsetneq [i_v, j_v]$ and $[i, j]$ is a P-gap at $u$. The induction hypothesis applied to $u$ implies that there is no leaf $x$ in $T_u$ such that $[i, j]$ is a 0-gap at $x$. And by Corollary 1, there is also no such leaf in $T_v$. If R1.6 is applied, then, without loss of generality, $[i_u, j_u] = [i - k, j]$ and $[i_v, j_v] = [i, j + l]$ with $k, l > 0$. By Corollary 1, it means that for every leaf $x$ in $T_u$, $x$ has a 0-gap $[i', j'] \supseteq [i - k, j] \supsetneq [i, j]$. Thus, it is impossible that $[i, j]$ is a 0-gap at any leaf of $T_u$. The same argument holds for $T_v$. This proves point 3 in the statement. $\square$
Recall that a candidate solution of the DPP is a phylogenetically correct alignment extension $A^+$ where only deletions occur. Optimal solutions are candidate solutions of minimum cost $c(A^+)$.

**Proposition 1.** At the end of the bottom-up phase, site $k$ is not in a labeled gap at $w$ if and only if $A_w^+[k] = 1$ for any candidate solution $A^+$ of the DPP.

*Proof.* First recall that site $k$ is in a labeled gap at $w$ if and only if $A_x[k] = 0, \forall x \in L(T_w)$ (Lemma 2).

Thus, if $k$ is not in a labeled gap at $w$, then at least one of its descendant leaves $x$ is such that $A_x[k] = 1$. But then any candidate solution $A^+$ must satisfy $A_w^+[k] = 1$, by Observation 1.2. This proves the forward implication in the statement.

To prove the other direction of the statement, assume $A_w^+[k] = 1$ for all candidate solutions of the DPP. By contradiction suppose that $k$ is in a labeled gap at $w$. Then, because $A_x[k] = 0, \forall x \in L(T_w)$, we can construct a candidate solution $A^+$ with a deletion occurring on an edge ancestral to $w$ and at an interval containing $k$. This implies that there exists a candidate solution with $A_w^+[k] = 0$, contradicting the initial assumption. □

**Lemma 4** ("Homogeneity of labeled gaps"). At the end of the bottom-up phase, if $[i, j]$ is a labeled gap at node $w$, then sites $i-1$ and $j+1$ are not in any labeled gap at $w$.

*Proof.* The statement is trivially true if $w$ is a leaf.

Suppose now that $w$ is an internal node and that the statement is true for $u$ and $v$, the two children of $w$. Let $[i_u, j_u]$ and $[i_v, j_v]$ be the labeled gaps at $u$ and $v$ respectively, such that $[i, j] = [i_u, j_u] \cap [i_v, j_v]$. By induction hypothesis, we know that sites $i_u - 1$ and $j_u + 1$ are not in any labeled gap at $u$. The same goes for $i_v - 1$ and $j_v + 1$ at $v$.

Because every labeled gap at $w$ is obtained as the intersection of one labeled gap at $u$ and one labeled gap at $v$, none of $\{i_u - 1, i_v - 1, j_u + 1, j_v + 1\}$ can be in a labeled gap at $w$. Because $\{i - 1, j + 1\} \subseteq \{i_u - 1, i_v - 1, j_u + 1, j_v + 1\}$, the statement follows. □

**Corollary 4.** If $[i, j]$ is set as a labeled gap at node $w$ in the bottom-up phase, then in every candidate solution $A^+$ of the DPP, $A_w^+[i-1] = A_w^+[j+1] = 1$.

*Proof.* Directly from Proposition 1 and Lemma 4. □

**Lemma 5.** Let $A$ be an alignment for set $X \supseteq \{x, y\}$. Suppose changing $A$ into $A'$ so that $A_x'[i, j]$ and $A_y'[i, j]$ are now only composed of 0s, while leaving the rest of $A_x$ and $A_y$ unchanged. This change

1. removes all indels between $x$ and $y$ acting at an interval contained in $[i, j]$;

2. cannot increase the number of the other indels between $x$ and $y$, i.e. those at an interval not contained in $[i, j]$.

*Proof.* Point 1 is trivially true. As for point 2, let us prove that the number of deletions at an interval not contained in $[i, j]$ cannot increase. The proof for insertions is symmetric.

Let $K$ be the set of sites $k \in [1, i-1] \cup [j+1, m]$ such that $(A_x[k], A_y[k]) = (1, 0)$. Define partition $P(A)$ of $K$ as the partition that puts two sites $k_1, k_2$ in
the same subset if and only if there is no site $h$ with $k_1 < h < k_2$ such that $(A_x[h], A_y[h]) = (1, 1)$. A different partition $P(A')$ of $K$ is obtained in the same way for $A'$ instead of $A$. Clearly $|P(A)|$ equals the number of deletions from $A_x$ to $A_y$ not contained in $[i, j]$, while $|P(A')|$ equals the number of deletions of from $A'_x$ to $A'_y$ not contained in $[i, j]$. Point 2 of the statement is then equivalent to proving $|P(A)| \ge |P(A')|$.

Now note that because $A'$ may remove from $A$ sites such that $(A_x[k], A_y[k]) = (1, 1)$, then $P(A')$ must be a coarser partition than $P(A)$, meaning that each element of $P(A)$ is a subset of some element of $P(A')$. This implies $|P(A)| \ge |P(A')|$. □

The following Observation introduces a way of partitioning the set of deletions in any $A^+$ candidate solution of the DPP. We will use this partition in the proofs of several lemmas below. (Namely Propositions 2, 3, 4 and Lemma 6.)

**Observation 2.** Let $A^+$ be a candidate solution of the DPP. Suppose $[i, j]$ is set as a labeled gap at node $w$ in the bottom-up phase. Let $p$ be the parent node of $w$. Define the following sets of deletions:

* $D_{p,w}$ contains the deletions from $p$ to $w$ at a subinterval of $[i, j]$,

* $D_{in}$ contains the deletions occurring in $T_w$ at a subinterval of $[i, j]$,

* $D_{out}$ contains the deletions occurring in $T_w$ at an interval that is not contained in $[i, j]$,

* $D_{else}$ contains the deletions occurring outside of $E(T_w) \cup \{(p, w)\}$ or in $(p, w)$ at an interval with no common site with $[i, j]$.

Then, $\{D_{p,w}, D_{in}, D_{out}, D_{else}\}$ is a partition of the set of deletions in $A^+$.

*Proof.* Any deletion in $A^+$ occurs on some edge $(x, y) \in E(T)$ at some interval $[i', j']$. Depending on whether $[i', j'] \subseteq [i, j]$, or $[i', j'] \cap [i, j] = \emptyset$ or the remaining case, and on whether $(x, y) = (p, w)$, or $(x, y) \in E(T_w)$ or $(x, y) \notin E(T_w) \cup \{(p, w)\}$, we have 9 possible cases, depicted in Table S1. The four sets defined above correspond to a partition of these 9 cases. The only case that they do not cover (bottom-left case in Table S1) is not possible: Corollary 4 implies $A_w^+[i-1] = A_w^+[j+1] = 1$, meaning that any deletion in $(p, w)$ is at an interval that is either contained in $[i, j]$ (i.e. the deletion is in $D_{p,w}$), or not overlapping with it (i.e. the deletion is in $D_{else}$). □

<table><caption>Table S1: Partitioning deletions with respect to a labeled gap $[i, j]$ at $w$. Suppose a deletion occurs on edge $(x, y)$ at interval $[i', j']$. The rows correspond to three different ways $[i', j']$ can be related to $[i, j]$. The columns correspond to the three cases $(x, y) = (p, w)$, $(x, y) \in E(T_w)$ and $(x, y) \notin E(T_w) \cup \{(p, w)\}$.</caption><thead><tr><th>[i', j']...</th><th>(p, w)</th><th>E(T<sub>w</sub>)</th><th>else</th></tr></thead><tbody><tr><td>... &sub; [i, j]</td><td>D<sub>p,w</sub></td><td>D<sub>in</sub></td><td>D<sub>else</sub></td></tr><tr><td>&cap;[i, j] = &empty;</td><td>D<sub>else</sub></td><td>D<sub>out</sub></td><td>D<sub>else</sub></td></tr><tr><td>else</td><td>&empty;</td><td>D<sub>out</sub></td><td>D<sub>else</sub></td></tr></tbody></table>

**Proposition 2.** If $[i, j]$ is set as a 0-gap at $w$ during the bottom-up phase then, in every optimal solution $A^+$ of the DPP, $A_w^+[i, j]$ is a gap.
**Proof.** If $w$ is a leaf, $A_w[i,j]$ must be a gap and therefore every extension $A^+$ of $A$ must be such that $A_w^+[i,j]$ is a gap. Now suppose $w$ is not a leaf. Since $[i,j]$ is set as a 0-gap for $w$ in the bottom-up phase, it means by Lemma 3 that there exist $x_u \in L(T_u)$ and $x_v \in L(T_v)$, with $u$ and $v$ children of $w$, such that $[i,j]$ is a 0-gap at both $x_u$ and $x_v$. And with Corollary 1, we know that every leaf $x$ in $T_w$ has a 0-gap $[ix, jx] \supseteq [i,j]$. We reason by contradiction and suppose there exists an optimal deletion-only solution $A^+$ where $[i,j]$ is not a gap for $w$. Let us partition the set of deletions in $A^+$ as described in Observation 2. We then have $c(A^+) = |D_{p,w}| + |D_{in}| + |D_{out}| + |D_{else}|$.

First note that $A_w^+[i-1] = A_w^+[j+1] = 1$, by Corollary 4. Since we assumed that $A_w^+[i,j]$ is not a gap, $\exists k \in [i,j]$ such that $A_w^+[k] = 1$. Note that at least one deletion at an interval including site $k$ must occur in the path from $w$ to $x_u$. This deletion must belong to $D_{in}$, because if it acted on an interval not contained in $[i,j]$ but overlapping with $[i,j]$, then this would contradict the fact that $[i,j]$ is a 0-gap at $x_u$. We can reason similarly for $x_v$, and thus conclude $|D_{in}| \ge 2$.

Now, let us take $B^+$, a copy of $A^+$, and modify $B^+$ by setting $B_x^+[h] = 0$, $\forall x$ descendant of $w$ and $\forall h \in [i,j]$. By doing so, we create a new deletion in $B^+$ between $w$ and its parent $p$, at an interval $[i_0, j_0] \subseteq [i,j]$ containing all sites $k \in [i,j]$ such that $A_w^+[k] = 1$. Note that $i_0 = \min\{k \in [i,j] : A_p^+[k] = 1\}$ and $j_0 = \max\{k \in [i,j] : A_p^+[k] = 1\}$. $B^+$ is clearly an extension of $A$ that only involves deletions, so it is a candidate solution of the DPP. Let $\{D'_{p,w}, D'_{in}, D'_{out}, D'_{else}\}$ be the partition of the deletions in $B^+$ defined in the same way as $\{D_{p,w}, D_{in}, D_{out}, D_{else}\}$ for $A^+$.

Now note that all deletions in $D_{in}$ are no longer deletions in $B^+$, by Lemma 5. The same holds for the deletions in $D_{p,w}$, which in $A^+$ act on intervals that are strictly contained in the interval $[i_0, j_0]$ of the new deletion from $p$ to $w$ in $B^+$. Thus, $|D'_{in}| = 0$ and $|D'_{p,w}| = 1$.

As for $D_{out}$, Lemma 5 implies that the number of deletions in $T_w$ at an interval not contained in $[i,j]$ cannot increase when switching from $A^+$ to $B^+$, meaning that $|D_{out}| \ge |D'_{out}|$.

Finally, we focus on $D_{else}$. First, it is easy to see that any deletion occurring outside $\{(p,w)\} \cup E(T_w)$ is unaffected by the changes between $A^+$ and $B^+$. Similarly, any deletion from $p$ to $w$ at an interval not overlapping with $[i,j]$ is unaffected by these changes because of $A_w^+[i-1] = A_w^+[j+1] = 1 = B_w^+[i-1] = B_w^+[j+1]$. Thus, $|D_{else}| = |D'_{else}|$.

In conclusion, $c(A^+) \ge 2+|D_{p,w}|+|D_{out}|+|D_{else}| > 0+1+|D'_{out}|+|D'_{else}| = c(B^+)$, which contradicts the assumption that $A^+$ is optimal. $\square$

**Proposition 3.** If $[i,j]$ is set as a P-gap at $w$ in the bottom-up phase then, in every optimal solution of the DPP, $A_w^+[i,j] = A_p^+[i,j]$ with $p$ being the parent of $w$.

**Proof.** We reason by contradiction and suppose there exists an optimal solution $A^+$ where $A_w^+[i,j] \ne A_p^+[i,j]$. We partition the set of deletions in $A^+$ as described in Observation 2. We then have $c(A^+) = |D_{p,w}|+|D_{in}|+|D_{out}|+|D_{else}|$. $A_w^+[i,j] \ne A_p^+[i,j]$ and $A_w^+[i-1] = A_w^+[j+1] = 1$ (consequence of Corollary 4), together imply $|D_{p,w}| \ge 1$.

Now set $B^+$ as a copy of $A^+$, and then $\forall x$ descendant of $w$, modify $B_x^+[i,j]$ as follows:
• if $A_x^+[i-1] = A_x^+[j+1] = 1$ then set $B_x^+[i,j] = A_p^+[i,j]$,

• otherwise set $B_x^+[h] = 0, \forall h \in [i, j]$.

Note that we necessarily set $B_w^+[i, j] = A_p^+[i, j]$, since $A_w^+[i-1] = A_w^+[j+1] = 1$ by Corollary 4. Moreover, by Lemma 3 the condition $A_x^+[i-1] = A_x^+[j+1] = 1$ cannot be true for any leaf in $T_w$. Thus, in any path between $w$ and one of the leaves in $T_w$, there exists an edge $(x, y)$ where that condition is true for $x$ but false for $y$, meaning that we set $B_x^+[i, j] = A_p^+[i, j]$ and $B_y^+[i, j] = 0 \dots 0$. Let $E$ be the set of such edges $(x, y)$. Because for any edge $(x, y) \in E$ at least one of $A_y^+[i-1]$ or $A_y^+[j+1]$ equals 0, in $A^+$ there exists a deletion from $x$ to $y$ at an interval containing at least one of $i-1$ or $j+1$. This deletion is necessarily an element of $D_{out}$. Moreover, $B^+$ is an extension of $A$ and only involves deletions, so it is a candidate solution of the DPP.

Let $\{D'_{p,w}, D'_{in}, D'_{out}, D'_{else}\}$ be the partition of the deletions in $B^+$ defined in the same way as $\{D_{p,w}, D_{in}, D_{out}, D_{else}\}$ for $A^+$. Note that $B_w^+[i, j] = A_p^+[i, j] = B_p^+[i, j]$ implies $|D'_{p,w}| = 0$. Moreover, for the same reasons as in the proof of Proposition 2, $|D_{else}| = |D'_{else}|$.

We now show that $|D_{out}| \ge |D'_{out}|$ by showing that for any $(x, y) \in E(T_w)$ there are at least as many deletions from $x$ to $y$ in $D_{out}$ than in $D'_{out}$. We consider three cases, depending on whether $(x, y) \in E$, or $(x, y)$ is a strict descendant of an element of $E$, or a strict ancestor. (We say that an edge $(x, y)$ is a strict descendant of another edge $(x', y')$ if $x$ and $y$ are strict descendants of $x'$ and $y'$, respectively. Strict ancestors are defined symmetrically.)

First, if $(x, y)$ is a strict descendant of an element of $E$, then $B_x^+[i, j] = B_y^+[i, j] = 0 \dots 0$ and, by Lemma 5, the number of deletions from $x$ to $y$ cannot be larger in $D'_{out}$ than in $D_{out}$.

Second, if $(x, y) \in E(T_w)$ is a strict ancestor of an element of $E$, then by construction all of $A_x^+[i-1], A_y^+[i-1], A_x^+[j+1], A_y^+[j+1]$ are equal to 1. But then, any deletion in $D_{out}$ and $D'_{out}$ from $x$ to $y$ must be at an interval that does not overlap with $[i-1, j+1]$. This means that the deletions from $x$ to $y$ are the same in $D_{out}$ and $D'_{out}$, and in particular that their number does not change.

Third, if $(x, y) \in E$, as we noted earlier this means that a deletion from $x$ to $y$ in $D_{out}$ occurs at an interval $[i', j']$ containing $i-1$ or $j+1$ (potentially both). But then, because of $B_y^+[i, j] = 0 \dots 0$, this deletion is replaced in $B^+$ by a potentially different deletion from $x$ to $y$ at an interval including $[i', j']$ and all sites $k \in [i, j]$ such that $B_x^+[k] = 1$. Clearly, this new deletion is an element of $D'_{out}$. Furthermore, any other deletion from $x$ to $y$ in $D'_{out}$ was already present at the same interval in $D_{out}$. Thus, for every deletion from $x$ to $y$ in $D'_{out}$ there exists at least one deletion in $D_{out}$ that is either identical to it, or occurring on a smaller interval, which implies that there are at least as many deletions from $x$ to $y$ in $D_{out}$ than in $D'_{out}$. This completes the proof of $|D_{out}| \ge |D'_{out}|$.

Finally, by considering the same three cases as we did to prove $|D_{out}| \ge |D'_{out}|$, it is easy to see that $|D'_{in}| = 0$. In conclusion, $c(B^+) = 0 + 0 + |D'_{out}| + |D'_{else}| < 1 + |D_{in}| + |D_{out}| + |D_{else}| \le c(A^+)$, which contradicts the optimality of $A^+$. $\square$

**Proposition 4.** If $[i, j]$ is set as a C-gap at $w$ in the bottom-up phase then, in every optimal solution $A^+$ of the DPP, $A_w^+[i, j]$ is a gap or $A_w^+[i, j] = A_p^+[i, j]$ with $p$ being the parent of $w$.
**Proof.** Since $[i, j]$ is set as a C-gap for $w$ in the bottom-up phase, by Lemma 3 one child of $w$, say $u$, has a descendant leaf $x_u$ such that $[i, j]$ is a gap at $x_u$, while the other child $v$ is such that all its descendant leaves have a gap that is a proper superset of $[i, j]$. We reason by contradiction and suppose there exists an optimal solution $A^+$ of DPP where $A_w^+[i, j]$ is not a gap and $A_w^+[i, j] \neq A_p^+[i, j]$. Let us partition the set of deletions in $A^+$ as described in Observation 2. We then have $c(A^+) = |D_{p,w}| + |D_{in}| + |D_{out}| + |D_{else}|$.

First note that $A_w^+[i-1] = A_w^+[j+1] = 1$, by Corollary 4. Since we assumed that $A_w^+[i, j]$ is not a gap, $\exists k \in [i, j]$ such that $A_w^+[k] = 1$. Note that at least one deletion at an interval including site $k$ must occur in the path from $w$ to $x_u$. This deletion must belong to $D_{in}$, because if it acted on an interval not contained in $[i, j]$ but overlapping with $[i, j]$, then this would contradict the fact that $[i, j]$ is a 0-gap at $x_u$. Thus, $|D_{in}| \ge 1$.

Moreover, the assumption that $A_w^+[i, j] \neq A_p^+[i, j]$ implies that there exists $k' \in [i, j]$ such that $A_p^+[k'] = 1$ and $A_w^+[k'] = 0$. Now note that $A_w^+[i-1] = A_w^+[j+1] = 1$ implies that the deletion from $p$ to $w$ that deletes site $k'$ must occur at a subinterval of $[i, j]$. Therefore, $|D_{p,w}| \ge 1$.

Now, define the candidate solution $B^+$ in the same way as in the proof of Proposition 2, that is, $B_x^+[h] = 0$, $\forall x$ descendant of $w$ and $\forall h \in [i, j]$, while $B^+$ is a copy of $A^+$ elsewhere. Let $\{D'_{p,w}, D'_{in}, D'_{out}, D'_{else}\}$ be the partition of the deletions in $B^+$ defined in the same way as $\{D_{p,w}, D_{in}, D_{out}, D_{else}\}$ for $A^+$. For the same reasons as those detailed in the proof of Proposition 2, we have $|D'_{in}| = 0, |D'_{p,w}| = 1, |D_{out}| \ge |D'_{out}|$ and $|D_{else}| = |D'_{else}|$.

Thus, $c(A^+) \ge 1+1+|D_{out}|+|D_{else}| > 0+1+|D'_{out}|+|D'_{else}| = c(B^+)$, which contradicts the assumption that $A^+$ is optimal.

□

**Lemma 6.** Let $w \in V(T) \setminus \{r(T)\}$ and let $p$ be the parent of $w$. If $[i, j]$ is set as a C-gap at $w$ in the bottom-up phase, then:

1. $\forall A^+$ optimal solution of DPP, such that $A_w^+[i, j]$ is a gap, $\exists B^+$ optimal solution of DPP, such that

* $B_w^+[i, j] = B_p^+[i, j]$ and

* $B_x^+[k] = A_x^+[k]$, whenever $x \in V(T) \setminus V(T_w)$ or $k \notin [i, j]$.

2. $\forall B^+$ optimal solution of DPP, such that $B_w^+[i, j] = B_p^+[i, j]$, $\exists A^+$ optimal solution of DPP, such that

* $A_w^+[i, j]$ is a gap and

* $A_x^+[k] = B_x^+[k]$, whenever $x \in V(T) \setminus V(T_w)$ or $k \notin [i, j]$.

**Proof.** Without loss of generality, suppose $u$ is the child of $w$ that has a descendant leaf $x_u$ such that $[i, j]$ is a gap at $x_u$, while the other child $v$ is such that all its descendant leaves have a gap that is a proper superset of $[i, j]$.

(Proof of 1.) Let $A^+$ be an optimal solution of DPP, such that $A_w^+[i, j]$ is a gap. If $A_p^+[i, j]$ is a gap, then we already have $A_w^+[i, j] = A_p^+[i, j]$, meaning that statement 1 is trivially true for $B^+ = A^+$. Now, suppose that $A_p^+[i, j]$ is not a gap. Let $K$ be the non-empty set of sites $k \in [i, j]$ such that $A_p^+[k] = 1$. Let us partition the set of deletions in $A^+$ as described in Observation 2. Thus, we have $c(A^+) = |D_{p,w}| + |D_{in}| + |D_{out}| + |D_{else}|$, with $|D_{p,w}| = 1, |D_{in}| = 0$. ($|D_{p,w}| = 1$
is because a deletion from *p* to *w* occurs at [min *K*, max *K*]; |*D*<sub>in</sub>| = 0 because $A_x^+[h] = 0, \forall x \in V(T_w), \forall h \in [i, j].$

Let us take $B^+$, a copy of $A^+$, and modify $B^+$ as follows: $\forall k \in K$ and $\forall x \in V(T_v) \cup \{w\}$, set $B_x^+[k] = 1$ if $A_x^+[i-1] = A_x^+[j+1] = 1$. This has the effect of setting $B_x^+[i,j] = B_p^+[i,j] (= A_p^+[i,j])$ for $x=w$ and possibly for several descendants of $v$. Up until the end of this proof, let $\{D_{p,w}', D_{in}', D_{out}', D_{else}'\}$ be the partition of the deletions in $B^+$ defined in the same way as $\{D_{p,w}, D_{in}, D_{out}, D_{else}\}$ for $A^+$. Clearly, $B_w^+[i,j] = B_p^+[i,j]$ implies $|D_{p,w}'| = 0$. Also note that $|D_{in}'| = 1$ because we created a deletion from $w$ to $u$ at interval [min $K$, max $K$]. Moreover, $|D_{else}'| = |D_{else}|$, because these deletions remain untouched.

Now let us focus on the “new” deletions in $B^+$, i.e. those in $B^+$ but not in $A^+$. Let $(x, y)$ be an edge where a new deletion occurs. It is easy to see that the only way for a new deletion to arise in $B^+$ at edge $(x, y)$ is that we have changed the sites in set $K$ from 0s in $A_x^+$ into 1s in $B_x^+$, but that these sites have remained 0s in $B_y^+$. This occurs for $(x, y) = (w, u)$, where the new deletion is the reason why $|D_{in}'| = 1$. The only other case when this occurs is for $(x, y) \in E(T_v) \cup \{(w, v)\}$ whenever $A_x^+[i-1] = A_x^+[j+1] = 1$ and $(A_y^+[i-1], A_y^+[j+1]) \neq (1, 1)$. This means that an “old” deletion from $x$ to $y$ in $D_{out}$ occurred at an interval $[i', j']$ containing $i-1$ or $j+1$ (potentially both). But then, it is easy to see that the new deletion in $B^+$ is identical to the old one, except that it occurs at the potentially larger interval $[\min(K \cup \{i'\}), \max(K \cup \{j'\})]$. This new deletion is clearly an element of $D_{out}'$.

In conclusion, aside from the new deletion in $D_{in}'$, any other new deletion in $B^+$ is an element of $D_{out}'$ that replaces an element of $D_{out}$ by extending its interval. This proves that $c(B^+) = 0 + 1 + |D_{out}'| + |D_{else}'| \le |D_{p,w}| + |D_{in}| + |D_{out}| + |D_{else}| = c(A^+)$. Thus, $B^+$ is an optimal solution of the DPP, such that $B_w^+[i, j] = B_p^+[i, j]$.

(Proof of 2.) Let $B^+$ be an optimal solution of DPP, such that $B_w^+[i, j] = B_p^+[i, j]$. If $B_p^+[i, j]$ is a gap, then since $B_w^+[i-1] = B_w^+[j+1] = 1$ (Corollary 4) $B_w^+[i, j]$ is also a gap. Thus statement 2 is trivially true for $A^+ = B^+$. Now, suppose that $B_p^+[i, j]$ is not a gap. Thus, we have $c(B^+) = |D_{p,w}'| + |D_{in}'| + |D_{out}'| + |D_{else}'|$, with $|D_{p,w}'| = 0$ and $|D_{in}'| \ge 1$. Let us take $A^+$, a copy of $B^+$, and modify $A_x^+$ by setting $A_x^+[h] = 0, \forall h \in [i, j]$ and $\forall x \in V(T_w)$.

For the same reasons detailed in the proof of Proposition 2 we have $|D_{p,w}| = 1$, $|D_{in}| = 0$, $D_{else}' = D_{else}$ and, by Lemma 5, $|D_{out}'| \ge |D_{out}|$. Therefore, we have $c(A^+) \le c(B^+)$. Thus $A^+$ is an optimal solution of DPP, such that $A_w^+[i, j]$ is a gap.

Finally note that both transformations described above ensure $A_x^+[k] = B_x^+[k]$ for $x \notin V(T_w)$ or $k \notin [i, j]$, thus completing the proof of the statement. □

In the following statement and in the proof of Theorem 2 we use the following notation.

**Definition 7 (concatenation).** Let $A$ and $B$ be two alignments over $X$, let $x \in X$ and $k \in [0, m]$. $A_x[0, k]B_x[k+1, m+1]$ denotes the array of $(m+2)$ characters obtained by concatenating $A_x[0, k]$ and $B_x[k+1, m+1]$.

**Lemma 7.** Suppose that site $k$ does not belong to any labeled gap at $w$, at the end of the bottom-up phase. For any two optimal solutions of the DPP $D^+$ and
$E^+$, there exists an optimal solution of the DPP $F^+$ such that

$$F_w^+ = D_w^+[0, k] E_w^+[k+1, m+1]. \qquad (S1)$$

*Proof.* By induction on the depth of $w$. If $w = r(T)$, then $D_w^+$ and $E_w^+$ are only composed of 1s by Observation 1 and the statement trivially holds.

Now suppose $w \neq r(T)$ and that the statement holds for $p$, the parent of $w$. That is, there exists an optimal $F^+$ such that

$$F_p^+ = D_p^+[0, k] E_p^+[k+1, m+1]. \qquad (S2)$$

$F_w^+$ is not guaranteed to satisfy Equation (S1), because of the possible presence of C-gaps at $w$ (as we show below). We will see, however, that $F^+$ guarantees the existence of another optimal solution that does satisfy Equation (S1).

First note that for every site $h$ outside labeled gaps at $w$, $F_w^+[h] = D_w^+[h] = E_w^+[h] = 1$ by Proposition 1, meaning that $F_w^+[h]$ is consistent with Equation (S1). Similarly, for every 0-gap $[i, j]$ at $w$, $F_w^+[i, j] = D_w^+[i, j] = E_w^+[i, j] = 0 \dots 0$, by Proposition 2, meaning that $F_w^+[i, j]$ is consistent with Equation (S1).

Moreover, if $[i, j]$ is a P-gap at $w$ with $i, j < k$ then $F_w^+[i, j] = D_p^+[i, j]$ (Proposition 3) and $F_p^+[i, j] = D_p^+[i, j]$ (Equation (S2)) imply $F_w^+[i, j] = D_w^+[i, j]$. Symmetrically, $F_w^+[i, j] = E_w^+[i, j]$ for every P-gap $[i, j]$ at $w$ with $i, j > k$. So, whenever $[i, j]$ is a P-gap, $F_w^+[i, j]$ is consistent with Equation (S1).

The arguments above show that, if there are no C-gaps at $w$, then $F_w^+$ satisfies Equation (S1). Now consider a C-gap $[i, j]$ at $w$, and suppose, without loss of generality, $i, j < k$. By Proposition 4 we know that $F_w^+[i, j] \in \{0 \dots 0, F_p^+[i, j]\}$ and $D_w^+[i, j] \in \{0 \dots 0, D_p^+[i, j]\}$, with $F_p^+[i, j] = D_p^+[i, j]$ (Equation (S2)). If $F_w^+[i, j] \neq D_w^+[i, j]$, Lemma 6 guarantees the existence of another optimal solution $G^+$ such that $G_w^+[i, j] = D_w^+[i, j]$. Moreover $G_w^+$ only differs from $F_w^+$ at the sites in $[i, j]$.

Now note that $G_w^+$ may still not satisfy Equation (S1) because of the presence of other C-gaps at $w$. In this case, repeated application of Lemma 6 in the way we described above leads to prove the existence of an optimal solution $H^+$ such that (a) $H_w^+[i, j] = D_w^+[i, j]$ for every C-gap $[i, j]$ with $i, j < k$, (b) $H_w^+[i, j] = E_w^+[i, j]$ for every C-gap $[i, j]$ with $i, j > k$, (c) $H_w^+$ coincides with $F_w^+$ for every site that is not in a C-gap. Therefore $H_w^+$ satisfies Equation (S1). $\square$

**Theorem 1.** $A^+$ is among the alignment extensions constructed by Algorithm 1 if and only if $A^+$ is an optimal solution of the DPP.

*Proof.* Recall that Algorithm 1 stores labeled gaps in a data structure $S$, which allows retrieving the set $S_u$ of labeled gaps at $u$, $\forall u \in V(T)$. At any point during execution, multiple versions of $S$ are stored in a set $Sol$. Each $S$ represents an “intermediate solution”, in the following sense. Let the “intermediate alignment extension” $I^+(S)$ be defined as follows:

* $I_w^+[k](S) = 0$ if $\exists$ a 0-gap $[i, j] \in S_w$ with $k \in [i, j]$;

* $I_w^+[k](S) = 1$ if, $\forall [i, j] \in S_w$, $k \notin [i, j]$;

* $I_w^+[k](S) = \text{undefined}$, if $\exists$ a C-gap or a P-gap $[i, j] \in S_w$ with $k \in [i, j]$.
In other words, $I^+(S)$ contains a 0 for every site that is in a 0-gap in $S$, a 1 for every site that is outside a labeled gap in $S$, and is undefined elsewhere. (Strictly speaking, $I^+(S)$ is not an alignment extension, because of the presence of undefined entries in it.) In the following we say that the intermediate solution $S$ is extendable to alignment extension $A^+$ if we can set all the undefined entries in $I^+(S)$ so as to make it equal to $A^+$ (or, equivalently, if $I^+(S) = A^+$ for all entries where $I^+(S)$ is defined).

With the definitions above in mind, it is possible to prove the correctness of Algorithm 1 by proving the following two invariant properties. As we will show, they hold at any point during the execution of Algorithm 1 after *Sol* has been initialized.

**Invariant 1** For any $S \in Sol$, $S$ is extendable to an optimal solution.

**Invariant 2** For any optimal solution $A^+$, there exists $S \in Sol$ extendable to $A^+$.

Before proving these two invariant properties, we note that they imply the Theorem's statement. This is because they must hold in particular at the very end of the execution of Algorithm 1. At this point, every $S \in Sol$ is such that all of the labeled gaps stored in it are 0-gaps, meaning that $I^+(S)$ has no undefined entries, and thus every $S \in Sol$ is only extendable to a single alignment extension $I^+(S)$. Thus Invariant 1 proves that, at termination, any $S \in Sol$ is such that $I^+(S)$ is an optimal solution. In other words, any returned solution is optimal. Invariant 2 proves that for any optimal solution $A^+$, there exists $S \in Sol$ such that $I^+(S) = A^+$. In other words, any optimal solution is among those returned.

Let us now prove the two invariant properties. First note that when *Sol* is initialized with *Sol* $\leftarrow \{S\}$, Proposition 1 says that the sites that are set to 1 in $I^+(S)$ are precisely those that equal 1 in all candidate solutions. Moreover by Proposition 2, all sites that are set to 0 in $I^+(S)$ are also equal to 0 in any optimal solution. These two observations mean that, initially, the only *S* in *Sol* is extendable to any optimal solution, which implies that both Invariant 1 and Invariant 2 are true at this stage.

As for later stages during execution, note that *Sol* can only be changed because of application of rules R2.1, R2.2, R2.3. Specifically, rules R2.1 and R2.3 replace an intermediate solution $S_0$ with another intermediate solution $S'_0$, while rule R2.2 replaces a intermediate solution $S_0$ with two intermediate solutions which we call $S_1$ and $S_2$ in this proof. To prove that the two invariant properties remain true after application of these rules, it suffices to prove that for any optimal solution $A^+$: (a) in rules R2.1 and R2.3, if $S_0$ is extendable to $A^+$, then $S'_0$ is also extendable to $A^+$; (b) in rule R2.2, if $S_0$ is extendable to $A^+$, then either $S_1$ or $S_2$ remains extendable to solution $A^+$, while the other is extendable to another optimal solution $B^+$.

Before considering each of rules R2.1, R2.2, R2.3, recall that they are applied to a labeled gap [*i*, *j*] at a node *w*. Since the root *r*(*T*) cannot have any labeled gaps (Corollary 2), we can assume *w* ≠ *r*(*T*). We call *p* the parent of *w*. Note that, by the time one of these rules is applied to a labeled gap at *w*, all of the C- and P-gaps at *p* have already been resolved, meaning that $I_p^+(S_0)$ does not contain any undefined site.
First, suppose rule R2.1 was applied, turning the C-gap $[i, j]$ at node $w$ into a 0-gap. This replaces intermediate solution $S_0$ with $S'_0$, so that while $I_w^+[i, j](S_0)$ was a stretch of undefined sites, $I_w^+[i, j](S'_0)$ is now a gap. Now suppose $S_0$ is extendable to optimal solution $A^+$. Rule R2.1 presupposes that $I_p^+[i, j](S_0)$ is a gap, which implies that also $A_p^+[i, j]$ is a gap. But then Proposition 4 implies that $A_w^+[i, j]$ is also a gap. Because $I_w^+[i, j](S'_0)$ and $A_w^+[i, j]$ are both gaps, and $I_w^+[i, j](S'_0)$ is all that changed in $I^+(S'_0)$ with respect to $I^+(S_0)$, $S'_0$ is also extendable to $A^+$, which proves point (a) above for rule R2.1.

Second, suppose rule R2.3 was applied to intermediate solution $S_0$, replacing it with $S'_0$. In terms of $I^+(S)$, while $I_w^+[i, j](S_0)$ was a stretch of undefined sites, rule R2.3 has the effect of setting $I_w^+[i, j](S'_0) = I_p^+[i, j](S_0)$. Now suppose $S_0$ is extendable to optimal solution $A^+$, which implies that $I_p^+[i, j](S_0) = A_p^+[i, j]$. Proposition 3 implies $A_w^+[i, j] = A_p^+[i, j]$, which together with $I_w^+[i, j](S'_0) = I_p^+[i, j](S_0) = A_p^+[i, j]$ proven above implies $A_w^+[i, j] = I_w^+[i, j](S'_0)$. This means that $S'_0$ is also extendable to $A^+$, which proves point (a) above for rule R2.3.

Finally, suppose rule R2.2 was applied to intermediate solution $S_0$, replacing it in $Sol$ with $S_1$ and $S_2$. While $I_w^+[i, j](S_0)$ was a stretch of undefined sites, $I_w^+[i, j](S_1)$ is now a gap, whereas $I_w^+[i, j](S_2) = I_p^+[i, j](S_0)$. Now suppose $S_0$ is extendable to an optimal solution $A^+$, which implies that $I_p^+[i, j](S_0) = A_p^+[i, j]$. Proposition 4 implies that either $A_w^+[i, j]$ is a gap, or $A_w^+[i, j] = A_p^+[i, j]$. If $A_w^+[i, j]$ is a gap, then $I_w^+[i, j](S_1) = A_w^+[i, j]$, meaning that $S_1$ is extendable to $A^+$. If $A_w^+[i, j] = A_p^+[i, j]$, then $I_w^+[i, j](S_2) = I_p^+[i, j](S_0) = A_p^+[i, j] = A_w^+[i, j]$, meaning that $S_2$ is extendable to $A^+$. Thus one of $S_1$ or $S_2$ remains extendable to $A^+$.

Now, without loss of generality we suppose $S_1$ is extendable to $A^+$ and prove that $S_2$ must be extendable to another optimal solution $B^+$. The proof that assuming $S_2$ is extendable to $A^+$ implies that $S_1$ is extendable to $B^+$ is symmetric.

Because $I_w^+[i, j](S_1)$ is a gap, then $A_w^+[i, j]$ must also be a gap. But then, by Lemma 6, there exists another optimal solution $B^+$ such that $B_w^+[i, j] = B_x^+[i, j]$ and such that $A_x^+[k] = B_x^+[k]$, whenever $x \notin V(T_w)$ or $k \notin [i, j]$. Because $I^+(S_2)$ is such that $I_w^+[i, j](S_2) = I_p^+[i, j](S_0)$, and $I_x^+[k](S_2) = I_x^+[k](S_1)$ whenever $x \neq w$ or $k \notin [i, j]$, it is easy to see that $S_2$ is extendable to $B^+$. In detail,

* for any node $x \notin V(T_w)$, $I_x^+(S_2) = I_x^+(S_1)$ and $B_x^+ = A_x^+$, so since here $S_1$ is extendable to $A^+$, then $S_2$ is extendable to $B^+$ for these nodes.

* for node $w$ and sites $k \notin [i, j]$, $I_w^+[k](S_2) = I_w^+[k](S_1)$ and $B_w^+[k] = A_w^+[k]$, so since here $S_1$ is extendable to $A^+$, then $S_2$ is extendable to $B^+$ for these sites.

* for node $w$ and sites in $[i, j]$, $I_w^+[i, j](S_2) = I_p^+[i, j](S_0) = A_p^+[i, j] = B_p^+[i, j] = B_w^+[i, j]$, so $S_2$ is identical (and thus extendable) to $B^+$ at these entries.

* for $x \in V(T_w) \setminus \{w\}$, because of the order with which the top-down phase (Algorithm 3) visits nodes, all the labeled gaps at $x$ in $S_2$ are the same that have been set by the bottom-up phase (Algorithm 2). We can thus apply Proposition 1, which implies that the sites that are set to 1 in $I_x^+(S_2)$ are precisely those that equal 1 in all candidate solutions. Moreover by
Proposition 2, all sites that are set to 0 in $I_x^+(S_2)$ are also equal to 0 in any optimal solution. These two observations imply that in the descendants of $w$, $S_2$ is extendable to any optimal solution, and in particular to $B^+$.

This concludes the proof of point (b) above, and completes the entire proof. □

# S10 Correctness of the graph construction phase

**Lemma 8.** Suppose that $k \in [0, m+1]$ is not in a labeled gap at $w$ at the end of the bottom-up phase. Then every path from 0 to $m+1$ in $G_w$ must traverse $k$.

*Proof.* If there existed a path from 0 to $m+1$ in $G_w$ not traversing $k$, then this path would have to traverse an arc $(i, j)$ with $i < k < j$. But the only way an arc $(i, j)$ between non-consecutive sites can be in $G_w$ is because of the application of one of rules R3.1 – R3.3 to a labeled gap containing $k$, which leads to a contradiction. □

**Theorem 2.** Let $V \subseteq [0, m+1]$ and $w \in V(T) \setminus L(T)$. $V$ is the vertex set of a directed path from 0 to $m+1$ in $G_w$ if and only if $\exists A^+$ optimal solution of the DPP such that $V = \mathbb{1}(A_w^+)$. □

*Proof.* For the sake of conciseness, in this proof we use the following notation:

$$ \mathcal{P}_w(i, j) = \{V \subseteq [i, j] : V \text{ is the vertex set of a path from } i \text{ to } j \text{ in } G_w\}. $$

Using the notation above, we prove a more general statement, which reduces to the Theorem’s statement for $k = 0, k' = m+1$:

**Claim.** Suppose $k, k' \in [0, m+1]$ ($k \le k'$) are not in any labeled gap at $w$ at the end of the bottom-up phase. Then, $\forall V \subseteq [k, k']$:

$$ V \in \mathcal{P}_w(k, k') \Leftrightarrow \exists A^+ \text{ optimal : } V = \mathbb{1}(A_w^+) \cap [k, k'.] $$

We prove this Claim by double induction, over the depth of $w$ and over the number of labeled gaps at $w$ between $k$ and $k'$. More precisely we prove the following statements which imply the Claim: (base case 1) the claim holds for $w = r(T)$; (base case 2) the claim holds for any node $w$ when there are no labeled gaps between $k$ and $k'$; finally, assuming that (induction hypothesis 1) the statement holds for $p$, the parent of $w$, and that (induction hypothesis 2) it holds for $w$ if there are up to $N$ labeled gaps between $k$ and $k'$, then (induction step) the statement holds for $w$ if there are $N+1$ labeled gaps at $w$ between $k$ and $k'$.

(Base case 2): If there are no labeled gaps between $k$ and $k'$, then by construction $G_w$ only has one path from $k$ to $k'$, the path that traverses every vertex in $[k, k']$. That is, $\mathcal{P}_w(k, k') = \{[k, k']\}$. Moreover, every solution is such that $A_w^+$ is only composed of 1s between $k$ and $k'$, by Proposition 1. That is, $\mathbb{1}(A_w^+) \cap [k, k'] = [k, k']$. Thus, $\forall V \subseteq [k, k']$, the following equivalences hold: $V \in \mathcal{P}_w(k, k') \Leftrightarrow V = [k, k'] \Leftrightarrow \exists A^+ \text{ optimal : } V = \mathbb{1}(A_w^+) \cap [k, k']$.

(Base case 1): If $w = r(T)$, then by Corollary 2 there can be no labeled gaps between $k$ and $k'$ and the proof of base case 1 applies.
Before proving the induction step, we need to prove 3 special cases, corresponding to the scenario where $k$ and $k'$ precisely delimit a labeled gap $[k+1, k'-1]$.

$([k+1, k'-1] \text{ 0-gap}):$

If $[k+1, k'-1]$ is a 0-gap at $w$, we can prove the Claim as follows:

$V \in \mathcal{P}_w(k, k') \Leftrightarrow^{(1)} V = \{k, k'\} \Leftrightarrow^{(2)} \exists A^+ \text{ optimal : } V = 1(A_w^+) \cap [k, k'].$

Step (1) holds because by construction of $G_w$, it only has one path from $k$ to $k'$ and $\mathcal{P}_w(k, k') = \{\{k, k'\}\}$. Step (2) uses the fact that, by Proposition 2, every optimal $A^+$ must be such that $A_w^+[k+1, k'-1]$ is a gap, meaning that $1(A_w^+) \cap [k, k'] = \{k, k'\}$.

$([k+1, k'-1] \text{ P-gap}):$

If $[k+1, k'-1]$ is a P-gap at $w$, then by construction, the subgraphs of $G_w$ and $G_p$ induced by $[k, k']$ are identical. We can prove the Claim as follows:

$V \in \mathcal{P}_w(k, k') \Leftrightarrow^{(1)} V \in \mathcal{P}_p(k, k') \Leftrightarrow^{(2)} \exists A^+ \text{ optimal : } V = 1(A_p^+) \cap [k, k'] \Leftrightarrow^{(3)} \exists A^+ \text{ optimal : } V = 1(A_w^+) \cap [k, k'].$

Step (1) holds because $\mathcal{P}_w(k, k') = \mathcal{P}_p(k, k')$ by construction of $G_w$. Step (2) holds because by Corollary 3, $k$ and $k'$ are not in any labeled gap at $p$, meaning that we can apply induction hypothesis 1. Step (3) holds because $A_w^+[k, k'] = A_p^+[k, k']$ by Proposition 3.

$([k+1, k'-1] \text{ C-gap}):$

If $[k+1, k'-1]$ is a C-gap at $w$, then by construction the subgraphs of $G_w$ and $G_p$ induced by $[k, k']$ are identical, except for the addition of arc $(k, k')$ to $G_w$. We can prove the Claim as follows:

$V \in \mathcal{P}_w(k, k') \Leftrightarrow^{(1)} V = \{k, k'\}$ or $V \in \mathcal{P}_p(k, k') \Leftrightarrow^{(2)}$

$V = \{k, k'\}$ or $\exists A^+$ optimal : $V = 1(A_p^+) \cap [k, k'] \Leftrightarrow^{(3)}$

$\exists A^+$ optimal : $V = 1(A_w^+) \cap [k, k'].$

Step (1) holds because $\mathcal{P}_w(k, k') = \mathcal{P}_p(k, k') \cup \{\{k, k'\}\}$ by construction of $G_w$. Step (2) holds again by Corollary 3 and the induction hypothesis 1. In step (3) the backward direction ($\Leftarrow^{(3)}$) holds because $1(A_w^+) \cap [k, k']$ must be equal to one of $\{k, k'\}$ or $1(A_p^+) \cap [k, k']$, by Proposition 4.

The forward direction ($\Rightarrow^{(3)}$) is slightly more complex. First suppose $V = \{k, k'\}$. In this case, since $\exists A^+$ optimal s.t. $A_w^+[k+1, k'-1]$ is a gap by Lemma 6, it means that $V = \{k, k'\} = 1(A_w^+) \cap [k, k']$. Now suppose $\exists A^+$ optimal s.t. $V = 1(A_p^+) \cap [k, k']$. In this case, since $\exists B^+$ optimal s.t. $B_w^+[k, k'] = B_p^+[k, k'] = A_p^+[k, k']$ by Lemma 6, it means that $V = 1(A_p^+) \cap [k, k'] = 1(B_w^+) \cap [k, k']$.

(Induction step):

Suppose there are $N+1$ labeled gaps at $w$ between $k$ and $k'$, with $N \ge 0$. Let $[i, j]$ be the rightmost labeled gap among these $N+1$ labeled gaps, meaning that we can decompose $[k, k']$ as: $[k, i-1] \cup [i, j] \cup [j+1, k']$ where $[k, i-1]$ contains $N$ labeled gaps, with both $k$ and $i-1$ outside of these labeled gaps, and $[j+1, k']$ does not contain, nor intersect, with any labeled gap. We can prove the induction step as follows:

$V \in \mathcal{P}_w(k, k') \Leftrightarrow^{(1)}$

$V \cap [k, i-1] \in \mathcal{P}_w(k, i-1)$ and

$V \cap [i-1, j+1] \in \mathcal{P}_w(i-1, j+1)$ and

$V \cap [j+1, k'] \in \mathcal{P}_w(j+1, k') \Leftrightarrow^{(2)}$

$\exists B^+$ optimal : $V \cap [k, i-1] = 1(B_w^+) \cap [k, i-1]$ and

$\exists C^+$ optimal : $V \cap [i-1, j+1] = 1(C_w^+) \cap [i-1, j+1]$ and

$\exists D^+$ optimal : $V \cap [j+1, k'] = 1(D_w^+) \cap [j+1, k'] \Leftrightarrow^{(3)}$
$\exists A^+$ optimal : $V \cap [k, k'] = 1(A_w^+) \cap [k, k']$.

Step (1) holds because every path from $k$ to $k'$ must traverse $i-1$ and $j+1$ by Lemma 8. Step (2) holds because:

$V \cap [k, i-1] \in \mathcal{P}_w(k, i-1) \Leftrightarrow \exists B^+$ optimal : $V \cap [k, i-1] = 1(B_w^+) \cap [k, i-1]$, by the induction hypothesis 2 (as there are $N$ labeled gaps between $k$ and $i-1$);

$V \cap [i-1, j+1] \in \mathcal{P}_w(i-1, j+1) \Leftrightarrow \exists C^+$ optimal : $V \cap [i-1, j+1] = 1(C_w^+) \cap [i-1, j+1]$, because as we have proven above the Claim is true when $k = i-1$ and $k' = j+1$ delimit a labeled gap;

$V \cap [j+1, k'] \in \mathcal{P}_w(j+1, k') \Leftrightarrow \exists D^+$ optimal : $V \cap [j+1, k'] = 1(D_w^+) \cap [j+1, k']$, by the base case 2 (as there are no labeled gaps between $j+1$ and $k'$).

Finally, while step (=$\Leftrightarrow$) is trivially true, step ($\Rightarrow$) can be proven by applying Lemma 7 twice: from $B^+$ and $C^+$ we can prove the existence of an optimal $E^+$ s.t. $E_w^+ = B_w^+[0, i-1]C_w^+[i, m+1]$, and from $E^+$ and $D^+$ we can prove the existence of an optimal $F^+$ s.t. $F_w^+ = E_w^+[0, j+1]D_w^+[j+2, m+1]$. (Recall Def. 7 for the concatenation of arrays of characters.) Noting that all solutions have 1s at sites $i-1$ and $j+1$ in $w$, all of the above implies

$$F_w^+[k, i-1] = B_w^+[k, i-1] \quad F_w^+[i-1, j+1] = C_w^+[i-1, j+1] \quad F_w^+[j+1, k'] = D_w^+[j+1, k']$$

Thus $F^+$ is such that

$V \cap [k, i-1] = 1(B_w^+) \cap [k, i-1] = 1(F_w^+) \cap [k, i-1]$ and

$V \cap [i-1, j+1] = 1(C_w^+) \cap [i-1, j+1] = 1(F_w^+) \cap [i-1, j+1]$ and

$V \cap [j+1, k'] = 1(D_w^+) \cap [j+1, k'] = 1(F_w^+) \cap [j+1, k']$

meaning that $V \cap [k, k'] = 1(F_w^+) \cap [k, k']$. Thus the Claim is verified with $A^+ = F^+$. This concludes the proof of step ($\Rightarrow$) and of the induction step. $\square$

# S11 Complexity analysis

## S11.1 Algorithm 1

**Remark 1.** In all the proofs that follow, we assume the following data structure for each stored $S$: $S$ is a map whose keys are the vertices in $V(T)$ and whose values $S_u$ are lists of labeled gaps. Moreover, for each $u \in V(T)$, the list $S_u$ is sorted by the indices of the labeled gaps, so that the labeled gaps with the smallest indices for their endpoints appear first.

Recall that $b$ denotes the number of boundaries in the input alignment $A$:

$$b = |\{k \in [0, m] : A[k] \neq A[k+1]\}|.$$

It is easy to see that $b$ is a very natural parameter to describe an alignment $A$. As we show in Lemma 9, the number $b$ also arises when counting the number of blocks, and of breakpoints in $A$, defined as follows.

Given an alignment $A$ for $X$, a block of $A$ is a maximal interval $[i, j]$ such that $A[k]$ is the same $\forall k \in [i, j]$. Given a gap $A_x[i, j]$, we say that its breakpoints are real numbers $i$ and $j+1$. The rationale behind the definition of breakpoint is that if we represent each site as a segments of length 1 over the real line (i.e., site 0 is the segment from 0.0 to 1.0, site 1 is the segment from 1.0 to 2.0, etc.) then a gap $[i, j]$ corresponds to the segment from $i$ to $j+1$. Similarly, if $[i, j]$ is a labeled gap at node $w$, its breakpoints are real numbers $i$ and $j+1$.
**Lemma 9 (Characterization of b).** Let $A$ be an alignment. The following propositions are equivalent:

1. A has $b$ boundaries.

2. A has $b+1$ blocks.

3. $|\{x \in \mathbb{R}: \text{ there exists a gap in } A \text{ which has } x \text{ as a breakpoint}\}| = b$

*Proof.* The equivalence between points 1 and 2 is trivial.

As for the equivalence between points 1 and 3, this can be shown by noting that $(k, k+1)$ is a boundary if and only if $k+1$ is a breakpoint for some gap in $A$. The forward direction (boundary $\Rightarrow$ breakpoint) can be proven as follows: if $k$ is such that $A[k] \neq A[k+1]$, then $\exists w \in L(T)$ such that $A_w[k] \neq A_w[k+1]$. But then two cases are possible: $A_w[k] = 1, A_w[k+1] = 0$, in which case $k+1$ is a breakpoint of a gap $[k+1, j]$ for some $j \ge k+1$; alternatively, $A_w[k] = 0, A_w[k+1] = 1$, in which case $k+1$ is a breakpoint of a gap $[i, k]$ for some $i \le k$. The converse implication (breakpoint $\Rightarrow$ boundary) can be shown analogously. $\square$

**Lemma 10.** Let $w \in V(T)$. At the end of the bottom-up phase, each breakpoint of a labeled gap at $w$ must be a breakpoint of some gap $A_x[i, j]$ for $x \in L(T_w)$.

*Proof.* By induction on the height of $w$. If $w$ is a leaf, then any labeled gap $[i, j]$ is a 0-gap and thus corresponds to a gap $A_w[i, j]$ with the same breakpoints as $[i, j]$.

If $w$ is an internal node, then the labeled gap $[i, j]$ must be obtained by the bottom-up phase as the intersection of labeled gaps $[i_u, j_u]$ at $u$ and $[i_v, j_v]$ at $v$, where $u$ and $v$ are the children of $w$. It is easy to see that a consequence of intersecting two intervals is that each breakpoint of $[i, j]$ is either a breakpoint of $[i_u, j_u]$ or of $[i_v, j_v]$. But by the inductive hypothesis the breakpoints of $[i_u, j_u]$ and of $[i_v, j_v]$ are breakpoints of some gap in $L(T_u) \cup L(T_v) = L(T_w)$. $\square$

**Lemma 11.** Let $w \in V(T)$. The number of labeled gaps at $w$ is at most $b/2$.

*Proof.* By Lemmas 9 and 10, any breakpoint of a given labeled gap must be one of the $b$ breakpoints of some gap in $A$. Moreover, considering that the labeled gaps at $w$ cannot overlap and cannot share a breakpoint (by Lemma 4), each breakpoint can be used at most once to create a labeled gap, and two breakpoints can create a labeled gap only if there is no other labeled gap between them. Thus, there are at most $b/2$ labeled gaps in $w$.

While in the argument above the number of labeled gaps is the one at the end of the bottom-up phase, the statement remains true for the number of 0-gaps at $w$ during the top-down phase. This is because 0-gaps inherit their breakpoints from labeled gaps at their ancestors, meaning that any breakpoint of a 0-gap must be one of the $b$ breakpoints of $A$. $\square$

**Lemma 12.** Assume that the input alignment is represented in gap form. Then the bottom-up phase runs in $O(np)$ time.

*Proof.* We are going to show here that the bottom-up phase processes a single node $w \in V(T)$ in $O(b)$ time, which yields the statement.

First, if $w \in L(T)$, each gap in $A_w$ is read and stored in $S_w$ as a 0-gap in constant time (lines 2 to 4 in Algorithm 2). Since there at most $O(b)$ such gaps (by Lemma 11), any leaf is processed in $O(b)$ time.
If $w$ is an internal node with children $u$ and $v$, $S_w$ can be efficiently computed from $S_u$ and $S_v$ in $O(b)$ time as follows (corresponding to line 9 in Algorithm 2). Initially $[i_u, j_u]$ is set to the first labeled gap in $S_u$ and $[i_v, j_v]$ to the first labeled gap in $S_v$. If $[i_u, j_u] \cap [i_v, j_v] \neq \emptyset$, then we can add this intersection to $S_w$. Next, if $j_u < j_v$ then $[i_u, j_u]$ is reset so that it represents the next labeled gap in $S_u$, if $j_v < j_u$ then $[i_v, j_v]$ is reset so that it represents the next labeled gap in $S_v$, and if $j_u = j_v$ then both $[i_u, j_u]$ and $[i_v, j_v]$ are reset to next labeled gaps in their respective lists. We iterate this procedure until the end of one of the two lists is reached. Clearly, this procedure requires $O(|S_u| + |S_v|) = O(b)$ time. Note that $|S_w| = O(b)$ for any node $w$, because of Lemma 11. Moreover, the label for each labeled gap (line 11) is computable in constant time, giving a total of $O(b)$ time to process node $w$. $\square$

**Theorem 3.** Assume that the input alignment $A$ and every solution $A^+$ in output is represented in gap form. Then Algorithm 1 runs in $O(nbs)$ time.

*Proof.* We start by analysing the running time required to duplicate solutions throughout the top-down phase (as specified by rule R2.2). Note that although $|Sol|$ can grow during execution, it is bounded by $s$. To reach the final number of solutions $s$, there must be exactly $s-1$ duplications throughout the top-down phase. Creating a new copy of $S$ can be done in $O(s)$ time, since $O(nb)$ is also the space required to store a solution $S$ (by Remark 1 and Lemma 11). Thus the total running time employed by duplications is $O(nbs)$.

We now consider the running time of the top-down phase (Algorithm 3 launched as TopDown($r(T)$, $Sol$)) *without* the duplication of solutions. We show that lines 3–9 inside the **for** $S \in Sol$ loop in Algorithm 1 can be executed in $O(b)$ time. For any fixed $S$, recall that by R2.2 some of its C-gaps at $w$ are reset as 0-gaps, while some other C-gaps are replaced by 0-gaps at $w$’s parent $p$. Resetting some of the C-gaps in $S_w^C$ as 0-gaps requires $O(|S_w^C|) = O(b)$ time. Moreover, it is easy to see that copying the 0-gaps at $p$ that fall within any P-gap at $w$, or within one of the remaining C-gaps in $S_w^C$, also requires $O(b)$ time in aggregate (because it can be done by scanning $S_p$ once, with $|S_p| = O(b)$ by Lemma 11).

Repeating the operations above once for each $S \in Sol$ and once for each internal node gives a total running time for the top-down phase of $O(nbs)$, excluding solution duplications.

Now recall that because the input alignment is in gap form, the bottom-up phase runs in $O nb$ time by Lemma 12. Finally note that returning $Sol$, whose elements are also in gap form, can also be done in $O(nbs)$ time. $\square$

## S11.2 Algorithm including the graph construction phase

**Proposition 5.** For any internal node $w$, $G_w$ is a planar graph.

*Proof.* Let $G_w = ([0, m+1], E_w)$. Below, we prove that for any pair of arcs $(h_1, k_1), (h_2, k_2) \in E_w$ one of the following two cases must be verified:

1. One of the two arcs occurs “to the left” of the other, meaning that we have either $h_1 < k_1 \le h_2 < k_2$ or $h_2 < k_2 \le h_1 < k_1$.

2. The two arcs are “nested”, meaning that either $[h_1, k_1] \subseteq [h_2, k_2]$ or $[h_1, k_1] \supseteq [h_2, k_2]$.
Clearly this implies the statement, as it means that drawing the graph with vertices $[0, m + 1]$ along a line in their natural ordering involves no crossing arcs.

We prove the necessity of points 1 or 2 by induction on the depth of $w$. If $w = r(T)$ then $E_w = \{(k, k+1) : k \in [0, m]\}$, which implies that case 1 is necessarily verified. Now suppose that points 1 or 2 are verified for node $p$, where $p$ is the parent of $w$. We show that one of them must also be true for $w$.

First consider the case where one of the two arcs connects consecutive sites, i.e. we can write $(h_1, k_1) = (i, i+1)$ without loss of generality. Let us consider the other arc. If $h_2 \ge i+1$ or $k_2 \le i$ then clearly point 1 is verified. Otherwise we must have $h_2 < i+1$ and $k_2 > i$, which implies point 2.

Now suppose that neither of $(h_1, k_1)$, $(h_2, k_2)$ connects consecutive sites. This implies that they both have been introduced by application of rules R3.1 – R3.3. Suppose that $(h_1, k_1)$ was introduced when processing labeled gap $[i_1, j_1]$ and that $(h_2, k_2)$ was introduced when processing labeled gap $[i_2, j_2]$. We now consider two cases: $[i_1, j_1] \ne [i_2, j_2]$ and $[i_1, j_1] = [i_2, j_2]$.

If $[i_1, j_1] \ne [i_2, j_2]$ then, by Lemma 4, these two labeled gaps cannot overlap, and they cannot cover consecutive sites (i.e. we cannot have $j_1 = i_2 - 1$ or $j_2 = i_1 - 1$). Without loss of generality let $[i_1, j_1]$ be to the left of $[i_2, j_2]$, i.e. $j_1 + 1 \le i_2 - 1$. But then, since the rules require $h_1, k_1 \in [i_1 - 1, j_1 + 1]$ and $h_2, k_2 \in [i_2 - 1, j_2 + 1]$, we must have $h_1, k_1 \le j_1 + 1 \le i_2 - 1 \le h_2, k_2$, which proves point 1 above.

Now suppose that both $(h_1, k_1)$ and $(h_2, k_2)$ were introduced when processing the same labeled gap $[i, j]$. If this is a 0-gap then by rule R3.1 we must have $(h_1, k_1) = (h_2, k_2) = (i-1, j+1)$ meaning that point 2 is verified. If $[i, j]$ is a P-gap, then $(h_1, k_1)$, $(h_2, k_2) \in E_p$ where $E_p$ is the arc set of $G_p$; but then either of points 1 or 2 must be verified, by the inductive hypothesis. Finally, if $[i, j]$ is a C-gap, then we can either have $(h_1, k_1)$, $(h_2, k_2) \in E_p$, a case that we have just treated, or at least one of the two arcs must be equal to $(i-1, j+1)$, which implies that the two arcs are nested, i.e. point 2.

This concludes the proof by induction that for every pair of arcs in $G_w$ either of points 1 or 2 holds, which implies that $G_w$ is planar. $\square$

**Corollary 5.** For any internal node $w$, the constructed graph $G_w = (V_w, E_w)$ is such that $|E_w| = O(m)$.

*Proof.* Because $G_w$ is planar (Proposition 5), we can deduce with Euler's formula for graphs that $|E_w| \le 3|V_w| - 6$. Thus, $G_w$ cannot have more than $3(m+2) - 6 = 3m$ arcs. $\square$

**Theorem 4.** The algorithm that first runs the bottom-up phase and then the graph construction phase runs in $O(nm)$ time.

*Proof.* By Lemma 12 the bottom-up phase runs in $O)$ time.

As for the graph construction phase, we now show that it takes $O(m)$ time to construct a single graph $G_w = (V_w, E_w)$. Clearly $V_w$ can be constructed in $O(m)$ time. Now consider the initialization of $E_w$ via rule R3.0, which depends on finding all sites that are not contained in any labeled gap of $S_w$. This can also be done in $O(m)$ time by considering each site in turn; for each site $k$ we only need to check whether $k$ belongs to one of two labeled gaps: either the labeled
gap [i, j] where we last found a previous site, or the labeled gap following [i, j] in $S_w$.

Next, we need to add to $E_w$ the arcs of rules R3.1, R3.2 and R3.3. It is easy to see that this can also be done in $O(m)$ time. One way to do this is to first consider all 0-gaps in $S_w$ and apply rule R3.1 to each. Then we consider all P-gaps in $S_w$. Here by rule R3.2, we need to find all arcs $(h, k) \in E_p$ such that $h, k \in [i-1, j+1]$ for some P-gap $[i, j] \in S_w$. To do this, for each $h \in [0, m]$ in ascending order, we consider all outgoing arcs $(h, k) \in E_p$ and for each arc we only need to check whether $h, k$ both belong to one of two P-gaps in $S_w$, similarly to what described in the last paragraph. Thus applying rule R3.2 to all P-gaps at once takes $O(|E_p|)$ time, which is $O(m)$ by Corollary 5. For the same arguments, applying rule R3.3 to all C-gaps takes $O(|E_p|) = O(m)$ time.

Thus, each graph is computed in $O(m)$ time. Since this must be repeated for all internal nodes, the graph construction phase requires $O(nm)$ time. Adding to this the bottom-up phase, the total complexity is $O(nb)+O(mn) = O(mn)$. $\square$

## S12 Chindelevitch et al.'s algorithm for DPP

Here, we show that the algorithm for DPP shown in Chindelevitch et al. [1], transcribed in Algorithm 4, may fail to return an optimal solution. Note that Chindelevitch et al. denote a deletion from $p$ to $u$ at interval $[i, j]$ as the pair $(u, [i, j])$. Moreover the edge set and vertex set of tree $T$ are denoted by $E_T$ and $V_T$, respectively.

<caption>**Algorithm 4:** DeletionParsimony (Alignment A, tree T)</caption>

```c
// Transcribed from Chindelevitch et al. [1]
1 Initialization:
2 for every leaf u do
3 | set $S_u = \{(i, j) | A_u(k) = 0 \forall i \le k \le j, A_u(i-1) = 1 = A_u(j+1)\}$
4 Recursion:
5 for each internal node $u \in V_T$, in a post-order traversal do
6 | Let $v$ and $w$ be the two children of $u$
7 | Set $R_v = \{(i, j) \in S_v | \exists(k, l) \in S_w$ with $(i, j) \subseteq (k, l)\}$
8 | Set $R_w = \{(i, j) \in S_w | \exists(k, l) \in S_v$ with $(i, j) \subseteq (k, l)\}$
9 | Set $S_u = (S_v - R_v) \cup (S_w - R_w)$
10 Termination:
11 Return $D = \bigcup_{u \in E_T} \bigcup_{(i,j) \in S_u - R_u} \{(u, [i, j])\}$.
```

Now consider the input depicted in Fig. S1, of a tree with 3 taxa, and an alignment of length 3 (+2 additional columns) of these taxa.

In the initialization, the algorithm computes a set $S_u$ for each leaf:

$$S_a = \{(1, 2)\}$$

$$S_b = \{(2, 3)\}$$

$$S_c = \emptyset$$
<figure><img src="image_12.png" /><figcaption>Figure S1: Input where the algorithm for DPP given by Chindelevitch et al. [1] fails to return an optimal solution. The unique optimal solution for this instance has two deletions: one from x to a at [1, 2] and the other from x to b at [2, 3].</figcaption></figure>

corresponding to the gaps at the leaves. Then, the recursion considers each internal node in post-order. When traversing node x, it sets

$$R_a = \{(i, j) \in S_a | \exists(k, l) \in S_b \text{ with } (i, j) \subseteq (k, l)\} = \emptyset$$

$$R_b = \{(i, j) \in S_b | \exists(k, l) \in S_a \text{ with } (i, j) \subseteq (k, l)\} = \emptyset$$

$$S_x = (S_a - R_a) \cup (S_b - R_b) = \{(1, 2); (2, 3)\}$$

Then, it considers the root r and sets

$$R_x = \{(i, j) \in S_x | \exists(k, l) \in S_c \text{ with } (i, j) \subseteq (k, l)\} = \emptyset$$

$$R_c = \{(i, j) \in S_c | \exists(k, l) \in S_x \text{ with } (i, j) \subseteq (k, l)\} = \emptyset$$

$$S_r = (S_x - R_x) \cup (S_c - R_c) = \{(1, 2); (2, 3)\}$$

Finally, let us consider the termination part, with the computing of D, that is the set of deletions found by the algorithm. Here line 11, as written in the published paper, sets

$$D = \bigcup_{u \in E_T} \bigcup_{(i,j) \in S_u - R_u} \{(u, [i, j])\}.$$

Note that $u \in E_T$ cannot be syntactically correct here, as sets $S_u, R_u$ and deletions $(u, [i, j])$ are defined for a node $u$, not an edge. We suppose that the authors mean that $u$ should be an internal node of $T$. Under this assumption,

$$D = \{(a, [1, 2]); (b, [2, 3]); (x, [1, 2]); (x, [2, 3])\}.$$

Clearly $D$ is suboptimal, as the (unique) optimal solution consist of deletions $\{(a, [1, 2]); (b, [2, 3])\}$.

We note that a way to fix Algorithm 4 for the particular input in Fig. S1 is to replace line 9 with

$$\text{Set } S_u = R_v \cup R_w.$$

However, it is easy to check that the resulting algorithm fails to return an optimal solution for the example in Fig. S2.

We also considered other ways to fix Algorithm 4. For example, replacing line 11 with

$$D = \bigcup_{(u,v) \in E_T} \bigcup_{(i,j) \in S_u - R_v} \{(u, [i, j])\}$$

also gives no consistent solution.
<figure><img src="image_13.png" /></figure>

<table><caption>Figure S2: Input where Algorithm 4 with line 9 replaced by “Set $S_u = R_v \cup R_w$” gives a suboptimal solution.</caption><thead><tr><th></th><th>0</th><th>1</th><th>2</th><th>3</th><th>4</th></tr></thead><tbody><tr><th>a</th><td>1</td><td>0</td><td>0</td><td>1</td><td>1</td></tr><tr><th>b</th><td>1</td><td>1</td><td>0</td><td>0</td><td>1</td></tr><tr><th>c</th><td>1</td><td>1</td><td>0</td><td>1</td><td>1</td></tr><tr><th>d</th><td>1</td><td>1</td><td>0</td><td>1</td><td>1</td></tr><tr><th>e</th><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td></tr></tbody></table>

# S13 Solving the DPP with edge-specific costs

Consider the following problem.

**Definition 8 (DPP WITH EDGE-SPECIFIC COSTS).**

*Input:* Tree $T$, alignment $A$ for $L(T)$, and, for each edge $(u, v) \in E(T)$, cost $c_{uv} > 0$.

*Output:* $A^+$, an extension of $A$, minimizing the cost

$$c(A^+) = \sum_{(u,v) \in E(T)} c_{uv} \cdot d(A_u^+, A_v^+)$$

under the constraint that for all $(p, w) \in E(T)$ where $p$ is the parent of $w$, all indels from $p$ to $w$ must be deletions.

In this problem each deletion has a cost that only depends on the edge where it occurs. This is useful to account for the fact that indel events are more likely to occur in some edges than others. Here we show how to modify the algorithm described in the main text, so as to solve this new problem.

The main idea is that the only thing that needs to be changed is the way labels 0, P and C are assigned to labeled gaps during the bottom-up phase. All other aspects of our work, including the graph-based representation of solutions, remain unchanged. In particular, if $[i, j]$ was a labeled gap at node $w$, $[i, j]$ remains a labeled gap in the new algorithm. Only its label may change. Later below we provide the pseudocode for the new bottom-up procedure (Algorithm 5) which replaces the old one (Algorithm 2). We do not provide a proof of correctness for this algorithm, but we describe the main ideas behind it and illustrate its functioning with an example.

## S13.1 Defining and computing $cost_w[i, j]($)

In order to assign the label of labeled gap $[i, j]$ at $w$, we need a few new definitions. Let $[i, j]$ be a labeled gap at node $w$. By Lemma 1, every leaf $x$ descending from $w$ must have a gap $A_x[i', j']$ such that $[i', j'] \supseteq [i, j]$. For brevity, we say that gap $A_x[i', j']$ *descends* from that labeled gap. Moreover, given a deletion from $u$ to $v$ at interval $[i, j]$, we say that it *covers* site $A_x[k]$ if $x \in L(T_v)$ and $k$ is one of the sites where $A_u^+[k] = 1$ and $A_v^+[k] = 0$. Finally, given any set of deletions, its *total cost* is trivially defined as the sum of the costs of the deletions in it.

**Definition 9.** Let $[i, j]$ be a labeled gap at $w$. Assume $A_w^+[k] = 1, \forall k \notin [i, j]$.
* $\text{cost}_w[i, j](\neg \emptyset)$ is the minimum total cost of a set of deletions in $T_w$ that cover all, and only, the sites in gaps descending from $[i, j]$.

* $\text{cost}_w[i, j](\emptyset)$ is the minimum total cost of a set of deletions in $T_w$ that cover all, and only, the sites $A_x[k]$ in gaps descending from $[i, j]$, such that $k \notin [i, j]$.

The intuition behind Definition 9 is the following: $\text{cost}_w[i, j](\emptyset)$ is the minimum cost to pay to set to 0 all gaps descending from $[i, j]$, under the assumption that $A_w^+[i, j]$ is already filled with 0s. The $\emptyset$ between the parentheses signifies that all sites in $[i, j]$ are assumed to be 0 at node $w$, and therefore also in any of its descendants $x \in V(T_w)$.

On the other hand, $\text{cost}_w[i, j](\neg \emptyset)$ is the minimum cost to pay to set to 0 all gaps descending from $[i, j]$, under the assumption that $A_w^+[i, j]$ is not filled with 0s. This explains the use of $\neg \emptyset$ between the parentheses. Note that, by definition, $\text{cost}_w[i, j](\emptyset) \le \text{cost}_w[i, j](\neg \emptyset)$.

Both quantities in Definition 9 can be computed recursively for all labeled gaps as follows. (We give a justification in the next subsection.) First, if $w$ is a leaf, for every labeled gap $[i, j]$ at $w$:

$$ \begin{aligned} \text{cost}_w[i, j](\emptyset) &= 0 \\ \text{cost}_w[i, j](\neg \emptyset) &= +\infty \end{aligned} \qquad (S3) $$

Now let $w$ be an internal node with children $u$ and $v$. Recall that any labeled gap $[i, j]$ at $w$ is obtained as the intersection of labeled gaps at its children, i.e. there exist $[i_u, j_u]$ and $[i_v, j_v]$, labeled gaps at $u$ and $v$ respectively, such that $[i, j] = [i_u, j_u] \cap [i_v, j_v]$. Then,

$$ \begin{aligned} \text{cost}_w[i, j](\neg \emptyset) = & \min \left\{ \begin{array}{l} c_{wu} + \text{cost}_u[i_u, j_u](\emptyset), \\ \text{cost}_u[i_u, j_u](\neg \emptyset) \end{array} \right\} + \\ & \min \left\{ \begin{array}{l} c_{wv} + \text{cost}_v[i_v, j_v](\emptyset), \\ \text{cost}_v[i_v, j_v](\neg \emptyset) \end{array} \right\} \end{aligned} \qquad (S4) $$

To compute $\text{cost}_w[i, j](\emptyset)$ we distinguish three cases. First, if $[i, j] = [i_u, j_u] = [i_v, j_v]$:

$$ \text{cost}_w[i, j](\emptyset) = \text{cost}_u[i, j](\emptyset) + \text{cost}_v[i, j](\emptyset). \qquad (S5) $$

Second, if the two labeled gaps in $w$'s children are strictly nested, that is, without loss of generality, if $[i, j] = [i_u, j_u] \subsetneq [i_v, j_v]$:

$$ \text{cost}_w[i, j](\emptyset) = \text{cost}_u[i, j](\emptyset) + \min \left\{ \begin{array}{l} c_{wv} + \text{cost}_v[i_v, j_v](\emptyset), \\ \text{cost}_v[i_v, j_v](\neg \emptyset) \end{array} \right\}. \qquad (S6) $$

Third, if $[i_u, j_u]$ and $[i_v, j_v]$ partially overlap:

$$ \text{cost}_w[i, j](\emptyset) = \text{cost}_w[i, j](\neg \emptyset). \qquad (S7) $$

Remark that $\text{cost}_w[i, j](\neg \emptyset)$ has already been computed with Eqn. S4.

## S13.2 Justification for the recurrences.

We focus on Eqn. S6, as it illustrates well all aspects behind equations S3 – S7. Similar arguments to the one below can be made for each of those equations.
Recall that Eqn. S6 assumes that $[i, j]$ is a labeled gap at $w$ and at its child $u$, and $[i_v, j_v] \supseteq [i, j]$ is a labeled gap at the other child $v$. We want to calculate $\text{cost}_w[i, j](\emptyset)$. Let $D_w$ be an optimal set of deletions satisfying the requirements in Definition 9. That is, $D_w$ is a minimum-cost set of deletions in $T_w$ covering all and only the sites $A_x[k]$ in gaps descending from $[i, j]$, but not the sites with $k \in [i, j]$ because these are already filled with 0s.

Because $A_w^+[i, j]$ is filled with 0s, there cannot be any deletion in $D_w$ from $w$ to $u$ at an interval inside $[i, j]$. Moreover, it is possible to show that any deletion from $w$ to $u$ at an interval outside $[i, j]$, i.e. not intersecting it or partially overlapping with it, would cover sites that are not in gaps that descend from the labeled gap $[i, j]$ at $w$. This is forbidden by Definition 9. So we can also exclude that $D_w$ contains any deletion from $w$ to $u$ at an interval outside $[i, j]$.

Because no deletion in $D_w$ can occur from $w$ to $u$, we can partition $D_w$ as

$$D_w = D_u \cup D_v \cup D_{wv}, \qquad (S8)$$

where $D_u$ contains the deletions of $D_w$ occurring in $T_u$, $D_v$ contains the deletions of $D_w$ occurring in $T_v$, and $D_{wv}$ contains the deletions occurring from $w$ to $v$.

Let us now focus on $D_{wv}$. Recall that $[i_v, j_v]$ is a labeled gap at $v$. Any deletion from $w$ to $v$ at an interval outside $[i_v, j_v]$ would cover sites that are not in gaps that descend from labeled gap $[i, j]$, so by Definition 9, $D_{wv}$ cannot contain any such deletion. Moreover there can be at most one deletion from $w$ to $v$ inside $[i_v, j_v]$, because two or more deletions at an interval within the same labeled gap would be suboptimal. All this means that $|D_{wv}| \le 1$.

Going back to the partition in (S8), we can then write $\text{cost}_w[i, j](\emptyset)$, which is the total cost of $D_w$, as the sum of the total costs of $D_u$, $D_v$ and $D_{wv}$, where the total cost of $D_{wv}$ is at most $c_{wv}$.

Recall that Eqn. S6 assumes that $A_w^+[i, j]$ is filled with 0s, which implies that the same holds for $A_u^+[i, j]$. Thus, $D_u$ is a set of deletions that satisfies all the requirements in the definition of $\text{cost}_u[i, j](\emptyset)$, and moreover it must have minimum total cost $\text{cost}_u[i, j](\emptyset)$ (because otherwise $D_w$ would not have minimum total cost itself). We have thus explained the first term in the right-hand side of Eqn. S6.

It remains to compute the total cost of $D_v \cup D_{wv}$. Here we must consider two cases: either $D_{wv}$ contains a deletion from $w$ to $v$ that sets to 0 all sites in $[i_v, j_v]$, in which case the total cost of $D_v \cup D_{wv}$ is $c_{wv} + \text{cost}_v[i_v, j_v](\emptyset)$, or $D_{wv} = \emptyset$, in which case we get $\text{cost}_v[i_v, j_v](\neg \emptyset)$. Naturally, because $D_v \cup D_{wv}$ must be optimal, its cost must be the minimum between these two options, i.e. the 2nd term of Eqn. S6: $\min\{c_{wv} + \text{cost}_v[i_v, j_v](\emptyset), \text{cost}_v[i_v, j_v](\neg \emptyset)\}$.

### S13.3 Assigning labels on the basis of $\text{cost}_w[i, j](\emptyset)$

For any labeled gap $[i, j]$ at $w$, once we have computed $\text{cost}_w[i, j](\emptyset)$ and $\text{cost}_w[i, j](\neg \emptyset)$, we can define its label as follows, where $p$ is the parent of $w$:

$$\text{label}_w[i, j] = \begin{cases} 0 & \text{if } c_{pw} + \text{cost}_w[i, j](\emptyset) < \text{cost}_w[i, j](\neg \emptyset), \\ P & \text{if } c_{pw} + \text{cost}_w[i, j](\emptyset) > \text{cost}_w[i, j](\neg \emptyset), \\ C & \text{if } c_{pw} + \text{cost}_w[i, j](\emptyset) = \text{cost}_w[i, j](\neg \emptyset). \end{cases} \qquad (S9)$$

Algorithm 5 shows how all the ideas presented here can be integrated in a new bottom-up phase.
### S13.4 A small example

Consider the small example in Fig. S3. We use it to illustrate: (1) the definition of $cost_w[i, j](X)$, (2) the recurrences, (3) the label assignment.

To illustrate the definition of $cost_w[i, j](X)$, consider labeled gap [3, 3] at node $x$. Its descendant gaps are $A_b[3, 3]$, $A_c[3, 3]$, $A_d[1, 3]$. Case $X = \emptyset$: if we assume that $A_x^+[3, 3] = 0$, then sites $A_b[3]$, $A_c[3]$ and $A_d[3]$ do not need to be covered; $cost_x[3, 3](\emptyset)$ is the minimum total cost of a set of deletions in $T_x$ that cover the remainder of the sites in the descendant gaps, that is, $A_d[1, 2]$. To cover them, we just need a single deletion from $y$ to $d$, which gives $cost_x[3, 3](\emptyset) = c_{yd} = 1$. Case $X = \neg \emptyset$: $cost_x[3, 3](\neg \emptyset)$ is the minimum total cost of a set of deletions in $T_x$ that cover exactly the sites in gaps $A_b[3, 3]$, $A_c[3, 3]$, $A_d[1, 3]$. To do this, it is optimal to take deletion [3, 3] from $x$ to $b$ (cost 2), deletion [3, 3] from $x$ to $y$ (cost 1), and deletion [1, 2] from $y$ to $d$ (cost 1). Hence, $cost_x[3, 3](\neg \emptyset) = 4$.

We now consider each labeled gap and show how $cost_w[i, j](X)$ is actually computed for each of them. First, Eqn. S4 sets

$$cost_y[1, 1](\neg \emptyset) = \min\{2 + 0, \infty\} + \min\{1 + 0, \infty\} = 3,$$

and Eqn. S6 sets

$$cost_y[1, 1](\emptyset) = 0 + \min\{1 + 0, \infty\} = 1.$$

The calculation for $cost_y[3, 3](X)$ is exactly the same as the one above for $cost_y[1, 1](X)$. Finally, Eqn. S4 sets

$$cost_x[3, 3](\neg \emptyset) = \min\{2 + 0, \infty\} + \min\{1 + 1, 3\} = 4,$$

and Eqn. S5 sets

$$cost_x[3, 3](\emptyset) = cost_b[3, 3](\emptyset) + cost_y[3, 3](\emptyset) = 0 + 1 = 1.$$

Finally, the labels of both labeled gaps at $y$ are set to 0, because $c_{xy} + 1 < 3$, and the label of [3, 3] at $x$ is set to C, because $c_{rx} + 1 = 4$. Note that this is very different to what we would obtain for the standard (constant costs) DPP, where we would assign label C to the labeled gaps at $y$ (by rule R1.4) and label 0 to [3, 3] at $x$ (by rule R1.1).

In the end, it is easy to check that the computed labels result in the two solutions such that $A_x^+[3, 3] \in \{0, 1\}$ and $A_y^+[1, 1] = A_y^+[3, 3] = 0$.

<figure><img src="image_14.png" /><figcaption>Figure S3: Resolution of a small instance of Problem 8. Left: input tree where cost $c_{uv}$ labels edge $(u, v)$. Right: table representation of $A^+$ showing labeled gaps as grey rectangles. Trivial labeled gaps (those at the leaves) are not shown for simplicity. For every other labeled gap, we report $cost_w[i, j](\emptyset)$ at the top of the rectangle, and $cost_w[i, j](\neg \emptyset)$ at the bottom. The assigned label appears to the right (after $\Rightarrow$).</figcaption></figure>
**Algorithm 5:** NewBottomUp($w$, $A$, $S$)

**Data:** $w$, a node of the binary tree $T$; $A$, an alignment for $L(T)$; $S$, a data structure containing a set $S_u$ of labeled gaps at $u$, $\forall u \in V(T)$; cost, a data structure to retrieve the values of $cost_u[i, j](\cdot)$, $\forall [i, j] \in S_u$

**Result:** $S_z$ and $cost_z$ are filled $\forall z \in V(T_w)$

// recall that here $S_u = \emptyset, \forall u \in V(T_w)$

if $w$ is a leaf then

$$ egin{array}{l} S_w \leftarrow \text{gaps in } A_w \\ \text{for } [i, j] \in S_w \text{ do} \\ \quad \text{cost}_w[i, j](\emptyset) = 0 \\ \quad \text{cost}_w[i, j](\neg \emptyset) = +\infty \\ \quad \text{assign label 0 to } [i, j] \end{array} $$

else

$$ egin{array}{l} \text{Let } u, v \text{ be the two children of } w \\ \text{BottomUp}(u, A, S) \\ \text{BottomUp}(v, A, S) \\ S_w \leftarrow \{[i_u, j_u] \cap [i_v, j_v] : ([i_u, j_u], [i_v, j_v]) \in S_u \times S_v\} \setminus \{\emptyset\} \\ \text{for } [i, j] \in S_w \text{ do} \\ \quad \text{Compute } cost_w[i, j](\neg \emptyset) \text{ using Eqn. S4} \\ \quad \text{Compute } cost_w[i, j](\emptyset) \text{ using one of Eqns. S5, S6, S7} \\ \quad \text{Let } p \text{ be the parent of } w \\ \quad \Delta = c_{pw} + cost_w[i, j](\emptyset) - cost_w[i, j](\neg \emptyset) \\ \quad L = \begin{cases} 0 & \text{if } \Delta < 0 \\ P & \text{if } \Delta > 0 \\ C & \text{if } \Delta = 0 \end{cases} \\ \quad \text{Assign label } L \text{ to } [i, j] \end{array} $$
### S13.5 Justification for Eqn. S9.

**Lemma 13.** Let $A^+$ be an optimal solution for Problem 8. Let $[i, j]$ be a labeled gap at node $w$. Let $D_w[i, j]$ be the set of deletions in $A^+$ occurring in $T_w$ and only covering sites in gaps descending from $[i, j]$.

$$ \text{The total cost of } D_w[i, j] = \begin{cases} \text{cost}_w[i, j](0) & \text{if } A_w^+[i, j] \text{ is a gap,} \\ \text{cost}_w[i, j](\neg 0) & \text{otherwise.} \end{cases} $$

We do not provide a proof for Lemma 13. The example in Figure S3 illustrates it well. Consider the optimal solution $A^+$ with $A_x^+[3, 3] = 1$ and $A_y^+[1, 1] = A_y^+[3, 3] = 0$. For this solution, $D_x[3, 3]$ contains three deletions: one from $x$ to $b$ at site 3, one from $x$ to $y$ at site 3, and one from $y$ to $d$ at site 2. Note that $D_x[3, 3]$ does not contain any deletion covering site 1. The total cost of $D_x[3, 3]$ equals 4, which is precisely $\text{cost}_x[3, 3](\neg 0)$, as predicted by Lemma 13. Now let $[i, j]$ be a labeled gap at node $w$ and let

* $D_w[i, j]$ be defined as in Lemma 13;

* $D_{pw}[i, j]$ be the set of deletions in $A^+$ from $p$ to $w$, at an interval contained in $[i, j]$.

For the sake of argument let us momentarily assume that $A_p^+[i, j]$ is not a gap. In this case, by Lemma 13, the total cost of $D_w[i, j] \cup D_{pw}[i, j]$ is:

$$ \begin{aligned} & c_{pw} + \text{cost}_w[i, j](0) && \text{if } A_w^+[i, j] \text{ is a gap,} \\ & = \text{cost}_w[i, j](\neg 0) && \text{if } A_w^+[i, j] = A_p^+[i, j], \end{aligned} \quad (S10) \\ & \ge c_{pw} + \text{cost}_w[i, j](\neg 0) && \text{otherwise.} \end{aligned} $$

Because $A^+$ is optimal, the total cost of $D_w[i, j] \cup D_{pw}[i, j]$ must be equal to the minimum across the three choices in Eqn. S10. Now let us consider three cases on the basis of the possible (in)equalities between $c_{pw} + \text{cost}_w[i, j](0)$ and $\text{cost}_w[i, j](\neg 0)$.

If $c_{pw} + \text{cost}_w[i, j](0) < \text{cost}_w[i, j](\neg 0)$, then the minimum of the three choices in Eqn. S10 is always the one in the first row. This means that, if $A_p^+[i, j]$ is not a gap, then $A_w^+[i, j]$ must be a gap. If, on the other hand, $A_p^+[i, j]$ is a gap, then by Observation 1, $A_w^+[i, j]$ is also a gap. Thus, if the $<$ relation holds in Eqn. S9, $A_w^+[i, j]$ must be a gap in any case, which justifies the fact that we assign the label 0 to it.

Similarly, if $c_{pw} + \text{cost}_w[i, j](0) > \text{cost}_w[i, j](\neg 0)$, then the minimum of the three choices in Eqn. S10 is always the one in the second row. This means that, if $A_p^+[i, j]$ is not a gap, then $A_w^+[i, j] = A_p^+[i, j]$. Clearly, this is also true if $A_p^+[i, j]$ is a gap. Since in any case $A_w^+[i, j] = A_p^+[i, j]$, this justifies the fact that if the $>$ relation holds in Eqn. S9, we assign label P.

Finally, it is intuitive that if $c_{pw} + \text{cost}_w[i, j](0) = \text{cost}_w[i, j](\neg 0)$, then both setting $A_w^+[i, j]$ as a gap or setting $A_w^+[i, j] = A_p^+[i, j]$ lead to an optimal solution, which justifies the assignment of label C in the third case of Eqn. S9. Note that the minimum of the three choices in Eqn. S10 can never be the one in the third row.
References

[1] Chindelevitch, L., Li, Z., Blais, E., and Blanchette, M. "On the inference of parsimonious indel evolutionary scenarios". In: *Journal of Bioinformatics and Computational Biology* 04.03 (2006), pp. 721-744.