
%\documentclass[twocolumn,10pt]{article}

% algorithmic uses latex2e command, but jair style is latex209

\documentstyle[jair,twoside,11pt,theapa,epsfig,graphicx,algorithmic]{article}
\input epsf %% at top of file

\jairheading{23}{2005}{647-683}{7/04}{5/05}
\ShortHeadings{An Improved Search Algorithm for Optimal Multiple-Sequence Alignment} 
{Schroedl}
\firstpageno{647}


%% \def\newcommand#1{%
%%    \def\@commandname{#1}%
%%    \@getoptionalarg\@continuenewcommand
%% }%
%% \def\@continuenewcommand{%
%%    \@numparams = \ifx\@optionalarg\empty 0\else\@optionalarg \fi \relax
%%    \@newcommand
%% }%
%% \def\@newcommand#1{%
%%    \def\@startdef{\expandafter\edef\@commandname}%
%%    \ifnum\@numparams=0
%%       \let\@paramdef = \empty
%%    \else
%%       \ifnum\@numparams>9
%%          \errmessage{\the\@numparams\space is too many parameters}%
%%       \else
%%          \ifnum\@numparams<0
%%             \errmessage{\the\@numparams\space is too few parameters}%
%%          \else
%%             \edef\@paramdef{%
%%                \ifcase\@numparams
%%                   \empty  No arguments.
%%                \or ####1%
%%                \or ####1####2%
%%                \or ####1####2####3%
%%                \or ####1####2####3####4%
%%                \or ####1####2####3####4####5%
%%                \or ####1####2####3####4####5####6%
%%                \or ####1####2####3####4####5####6####7%
%%                \or ####1####2####3####4####5####6####7####8%
%%                \or ####1####2####3####4####5####6####7####8####9%
%%                \fi
%%             }%
%%          \fi
%%       \fi
%%    \fi
%%    \expandafter\@startdef\@paramdef{#1}%
%% }%



% correct bad hyphenation here
% \hyphenation{op-tical net-works semi-conduc-tor IEEEtran}

%\usepackage{algorithmic}

\def\pea{\emph{PEA}$^*$}
\def\ast{$A^*$}
\def\ida{\emph{IDA}$^*$}
\def\DP{\emph{DP}}
\def\iddp{\emph{IDDP}}
\def\bali{\mbox{\emph{BAliBASE}}}
%\newcommand{\GRAPHIC}[1]{\includegraphics[width=0.5\columnwidth,height=0.5\columnwidth]{#1}}
\def\GRAPHIC#1{\includegraphics[width=0.5\columnwidth,height=0.5\columnwidth]{#1}}

\begin{document}

\sloppy

% paper title
\title{An Improved Search Algorithm\\ for Optimal Multiple-Sequence Alignment} 


% author names and affiliations
% use a multiple column layout for up to three different
% affiliations
%\author{\authorblockN{Stefan Schroedl\\
%\authorblockA{848 14th St\\
%%San Francisco CA 94114\\
%stefan.schroedl@gmx.de}}}

%% \author{Stefan Schroedl\\
%% 848 14th St\\
%% San Francisco CA 94114\\
%% (415) 522-1148\\
%% stefan.schroedl@gmx.de}


\author{\name Stefan Schroedl \email stefan.schroedl@gmx.de\\
       \addr  848 14th St\\
 San Francisco CA 94114\\ +1 (415) 522-1148}
   


% make the title area
\maketitle

\begin{abstract}

\emph{Multiple sequence alignment (MSA)} is a ubiquitous problem in
computational biology. Although it is \emph{NP}-hard to find an
\emph{optimal} solution for an arbitrary number of sequences, due to
the importance of this problem researchers are trying to push the
limits of exact algorithms further. Since MSA can  be cast as a
classical path finding problem, it is attracting a growing number of
AI researchers interested in heuristic search algorithms
as a challenge with actual practical relevance.

In this paper, we first review two previous, complementary lines of
research. Based on \emph{Hirschberg's algorithm}, \emph{Dynamic
  Programming} needs  $O(kN^{k-1})$ space to store both the search
frontier and the nodes needed to reconstruct the
solution path, for $k$ sequences of length $N$. \emph{Best first
  search}, on the other hand, has the advantage of bounding the search
space that has to be explored  using a heuristic. However, it is
necessary to maintain all explored nodes up to the final solution in
order to prevent the search from  re-expanding them at higher
cost. Earlier approaches to reduce the \emph{Closed} list are either
incompatible with pruning methods for the \emph{Open} list, or must
retain at least the boundary of the \emph{Closed} list. 

In this article, we present an algorithm that attempts at combining
the respective advantages; like \ast\, it uses a 
heuristic for pruning the search space, but reduces both the maximum
\emph{Open} 
and \emph{Closed} size to  $O(kN^{k-1})$, as in Dynamic
Programming. 
The underlying idea is to conduct a series 
of searches with successively increasing upper bounds, but using the
DP ordering as the key for the \emph{Open} priority queue. With a suitable
choice of thresholds, in practice, a running
time below four times that of~\ast\  can be expected.

In our experiments we show that our algorithm  outperforms  one of the  
currently most successful algorithms for optimal multiple sequence
alignments, \emph{Partial Expansion} \ast, both in time
and memory. Moreover, we apply a \emph{refined heuristic} based on
optimal alignments not only of pairs of 
sequences, but of larger subsets. This  idea is not new; however, to
make it practically relevant we show that it is equally important to
bound the heuristic computation appropriately, or the overhead can
obliterate any possible gain. 

Furthermore, we discuss a number of improvements in time and space
efficiency with regard to practical implementations. 

Our algorithm, used in conjunction with higher-dimensional heuristics,
is able to calculate for the first time the optimal
alignment for almost all of the problems in Reference~1 of the
benchmark database \bali. 


%\tableofcontents

\end{abstract}

%\vspace{0.5cm}

%\noindent {\bf Keywords:} Multiple Sequence Alignment, Heuristic
%Search, Dynamic Programming

%\pagebreak
% no keywords


\section{Introduction: Multiple Sequence Alignment}

The \emph{multiple sequence alignment problem (MSA)} in computational
biology consists in aligning several sequences, e.g.\  related genes
from different organisms, in order to reveal similarities and
differences across the group. Either DNA can be directly compared, and 
the underlying alphabet $\Sigma$ consists of the   
 set $\{${\tt C,G,A,T}$\}$ for the four standard nucleotide bases
 cytosine, guanine, adenine and thymine; or we can compare proteins,
 in which case $\Sigma$ comprises the twenty amino  acids.

 Roughly speaking, we try to write the sequences one above the other 
such that the columns with matching 
letters are maximized; thereby gaps (denoted here by an additional
letter~``{\tt \_}'') may be inserted into either of
them in order to shift the remaining characters into better
corresponding positions. Different letters in the same column can be
interpreted as being caused by point mutations during the course of
evolution that substituted one amino acid by another one; 
gaps can be seen as insertions or deletions (since the direction of
change is often not known, they are also  collectively referred to as
\emph{indels}). Presumably, the alignment with the fewest mismatches
or indels constitutes the biologically most plausible explanation.

There is a host of applications of MSA within computational biology;
e.g., for determining the evolutionary relationship between
species, for detecting functionally active sites which tend to be
preserved best across homologous sequences, and for predicting 
three-dimensional protein structure.


Formally, one 
associates a cost with an alignment and tries to 
find the (mathematically) \emph{optimal} alignment, i.e., that one
with minimum cost. When designing a cost function, computational
efficiency and 
biological meaning have to be taken into account. The most widely-used definition is the 
\emph{sum-of-pairs} cost function. First, we are given a
symmetric $(|\Sigma|+1)^2$ matrix containing penalties (scores) for
substituting a letter 
with another one (or a gap). In the simplest case, this could be one
for a mismatch and zero for a match, but more biologically relevant
scores have been developed. \fullciteA{day} have
proposed a model of molecular evolution where they estimate the
exchange probabilities of amino acids for different amounts of
evolutionary divergence; this gives rise to the so-called \emph{PAM
  matrices}, where PAM250 is generally the most widely used; \fullciteA{pet91} refined the statistics based on a larger body of   
experimental data. Based on
such a substitution matrix, the sum-of-pairs cost of an alignment is
defined as the sum of penalties between all letter pairs in
corresponding column positions. 

A pairwise alignment can be conveniently depicted as a path between two
opposite corners in a two-dimensional grid (Needleman and Wunsch, 1981): one sequence
is placed on 
the horizontal axis from left to right, the other one on the vertical
axis from top to bottom. If there is no gap in either string, the path
moves diagonally down and right; a gap in the vertical (horizontal) string is
represented as a horizontal (vertical) move right (down), since a
letter is consumed in only one of the strings. The alignment graph is
directed and acyclic, where a (non-border) vertex has incoming edges from the
left, top, and top-left adjacent vertices, and outgoing edges to the
right, bottom, and bottom-right vertices.

Pairwise alignment can be readily generalized to the simultaneous alignment
of multiple sequences, by considering higher-dimensional lattices. For
example, an alignment of three sequences can be visualized as a path
in a cube. Fig.~\ref{fig:cube} illustrates an example for the
strings {\tt ABCB}, {\tt BCD}, and {\tt DB}. It also shows
the computation of the sum-of-pairs cost, for a hypothetical
substitution matrix. A real example (problem {\tt 2trx} of \bali, see
Sec.~\ref{sec:bali}) is given in Fig.~\ref{fig:2trx}.

\begin{figure}[tb]
\begin{center}

%\centering{
%\begin{minipage}{0cm}
\begin{verbatim}
       Alignment:       Substitution matrix:
       
       A B C _ B           A B C D _
       _ B C D _         A 0 2 4 2 3
       _ _ _ D B         B   1 3 3 3
                         C     2 2 3
 Cost: 6+7+8+7+7 = 35    D       1 3
                         _         0

\end{verbatim}
%\end{minipage}
%}

%\epsfxsize7cm
%\epsfysize7cm
%\GRAPHIC{3sa}
\includegraphics[width=0.5\columnwidth]{3sa}
%\vspace{1cm}

\end{center}
\caption{Fictitious alignment problem: Column representation, cost
  matrix, three-dimensional visualization of the alignment path through the cube.\label{fig:cube}} 
\end{figure}


\begin{figure}%[hbt]
\begin{center}
{%\fontsize{5}{12}
%\tiny
\begin{verbatim}
  1thx    _aeqpvlvyfwaswcgpcqlmsplinlaantysdrlkvvkleidpnpttvkkyk______vegvpal
  1grx    __mqtvi__fgrsgcpysvrakdlaeklsnerdd_fqyqyvdiraegitkedlqqkagkpvetvp__
  1erv    agdklvvvdfsatwcgpckmikpffhslsekysn_viflevdvddcqdvasece______vksmptf
  2trcP   _kvttivvniyedgvrgcdalnssleclaaeypm_vkfckira_sntgagdrfs______sdvlptl
 
  1thx    rlvkgeqildstegvis__kdkllsf_ldthln_________
  1grx    qifvdqqhiggytdfaawvken_____lda____________
  1erv    qffkkgqkvgefsgan___kek_____leatine__lv____
  2trcP   lvykggelisnfisvaeqfaedffaadvesflneygllper_

\end{verbatim}
}

\end{center}

\caption{Alignment of problem {\tt 2trx} of \bali, computed with
  algorithm settings as described in Sec.~\ref{sec:bali}.\label{fig:2trx}} 
\end{figure}



A number of improvements can be integrated into the sum-of-pairs cost,
like associating weights with sequences, and using different
substitution matrices for sequences of varying evolutionary distance.  
A major issue in multiple sequence alignment algorithms is their
ability to handle gaps. Gap penalties can be made dependent on the
neighbor letters. Moreover, it has been found \cite{altgap}
that assigning a fixed score for each indel sometimes does not produce
the biologically most plausible alignment. Since the insertion of a
sequence of $x$ letters is more likely than $x$ separate insertions of
a single letter, gap cost 
functions have been introduced that depend on the length of a gap. A
useful approximation are \emph{affine gap costs}, which distinguish between
opening and extension of a gap and charge $a + b *x$ for
a gap of length $x$, for appropriate $a$ and $b$. Another frequently
used modification is to waive the penalties for gaps at the beginning
or end of a sequence.

Technically, in order to deal with affine gap costs we can no longer
identify nodes in the search graph with lattice vertices, since the
cost associated with an edge depends on the preceding edge in the
path. Therefore, it is more suitable to store lattice edges in the
priority queue, and let the transition costs for $u\rightarrow v, v
\rightarrow w$ be the sum-of-pairs substitution costs for using one
character from each sequence or a gap, plus the incurred gap penalties
for $v\rightarrow w $ followed by $u\rightarrow v$. This representation was
adopted in the program \emph{MSA} \fullcite{gupta96}. Note that
the state space in this representation grows by a factor of $2^k$. 
 An example of how successor costs are calculated, with the cost
 matrix of Fig.~\ref{fig:cube} and a gap opening penalty of 4, is
 shown in Fig.~\ref{fig:edgecost}.

\begin{figure}
\begin{center}
%\epsfxsize6cm
%\epsfysize7cm
\includegraphics[width=0.35\columnwidth]{edgecost}
%\GRAPHIC{edgecost}
\end{center}
\caption{Example of computing path costs with affine gap function;
 the substitution matrix of Fig.~\ref{fig:cube} and a gap opening
 penalty of 4 is used.\label{fig:edgecost}} 
\end{figure}

For convenience of terminology in the sequel we will still refer to
nodes when dealing with the search algorithm.

\section{Overview}

\fullciteA{wang} have shown that the optimal
multiple sequence alignment problem is \emph{NP}-hard; therefore, we
cannot hope to achieve an efficient algorithm for an arbitrary number
of sequences. As a consequence, alignment tools most widely used in
practice sacrifice the sound theoretical basis of exact algorithms,
and are heuristic in nature \fullcite{chan}. A wide variety of techniques
has been developed. \emph{Progressive} methods build up the alignment
gradually, starting with the closest sequences and successively adding
more distant ones. \emph{Iterative} strategies refine an initial alignment
through a sequence of improvement steps. 



Despite their limitation to moderate number of sequences, however, the
research into exact algorithms is still going on, trying to push the
practical boundaries further. They still
form the building block of heuristic techniques, and incorporating
them into existing tools could improve them. For example, an algorithm
iteratively aligning two groups of sequences at a time could do this
with three or more, to better avoid local minima. Moreover, it is
theoretically important to have the ``gold standard'' available for
evaluation and comparison, even if not for all problems. 

Since MSA can be cast as a minimum-cost path finding problem, it turns
out that it is amenable to heuristic search algorithms developed in
the AI community; these are actually among the currently best
approaches. Therefore, while many researchers in this area have 
often used puzzles and games in the past to study heuristic search
algorithms, recently there has been a rising interest in MSA as a
testbed with practical relevance, e.g., 
\fullcite{korf-dcbds,korf-dcfs,PE,zhou03b};
%(Korf, 1999; Korf and Zhang, 2000; Yoshizumi, Miura, and Ishida, 2000; Zhou and Hansen, 2003); 
its study has also led to 
major improvements of general search techniques.

It should be pointed out that the definition of the MSA problem as
given above is not the only one; it competes with other attempts
at formalizing biological meaning, which is often imprecise
or depends on the type of question the biologist investigator is
pursuing. E.g., in this paper we are only concerned with \emph{global 
  alignment} methods, which find an alignment of entire
sequences. Local methods, in contrast, are geared towards finding
maximally similar partial sequences, possibly ignoring the
remainder. 

In the next section, we briefly review previous
approaches, based on dynamic programming and incorporating lower and
upper bounds. In Sec.~\ref{sec:iddp}, we describe a new algorithm that
combines and extends some of these 
ideas, and allows to reduce the storage of \emph{Closed} nodes by
partially recomputing the solution path at the end
(Sec.~\ref{sec:sparse}). Moreover, it turns out that our algorithm's
iterative deepening strategy can be transferred to find a good balance
between the computation of \emph{improved heuristics} and the main
search (Sec.~\ref{sec:ndimheu}), an
issue that has previously been a major obstacle for their practical
application. Sec.~\ref{sec:exp} presents an experimental
comparison with \emph{Partial Expansion} \ast\ \fullcite{PE}, 
one of the currently most successful approaches. We also solve all but
two problems of Reference 1 of the widely used benchmark database
\bali\  \cite{bali}. To the best of our knowledge, this has not been
achieved previously with an exact algorithm.


\section{Previous Work}

A number of exact algorithms have been developed previously that can
compute alignments of a moderate number of sequences. Some of them are
mostly constrained by available memory, some by the required
computation time, and some on both. We can roughly group them into two
categories: those based on the dynamic programming paradigm, which
proceed primarily in breadth-first fashion; and best-first search,
utilizing lower and upper bounds to prune the search space. Some recent
research, including our new algorithm introduced in
Sec.~\ref{sec:iddp}, attempts to beneficially combine these approaches.

\subsection{Dijkstra's Algorithm and Dynamic Programming}

\citeA{Dijkstra} presented a general algorithm for
finding the shortest (resp.\ minimum cost) path in a directed graph.
It uses a \emph{priority queue (heap)} to store nodes $v$ together with the
shortest found distance from the start node $s$ (i.e., the top-left corner of
the grid) to $v$ (also called the
$g$-value of $v$). Starting with only $s$ in the priority queue, 
in each step, an edge with the minimum $g$-value is removed from the priority 
queue; its \emph{expansion} consists in generating all of its successors
(vertices to the right and/or below) reachable in one step,
computing their respective $g$-value by adding the edge cost to the
previous $g$-value, and inserting them in turn into the priority queue
in case this newly found distance is smaller than their previous
$g$-value. By the time a node is expanded, the $g$-value is
guaranteed to be the minimal path cost from the start node,
$g^*(v)=d(s,v)$. The 
procedure runs until the priority queue becomes empty, or the target
node $t$ (the bottom-right corner of the grid) has been reached; its $g$-value
then constitutes the optimal solution cost $g^*(t) = d(s,t)$ of the alignment
problem.  In order to trace back the 
path corresponding to this cost, we move backwards to the start node choosing
predecessors with minimum cost. The nodes can either be stored in a
fixed matrix structure 
corresponding to the grid, or they can be dynamically generated; in
the latter case, we can explicitly store at each node a
backtrack-pointer to this optimal parent.  

For integer edge costs, the priority queue can 
be implemented as a bucket array pointing to doubly linked lists
\cite{Dial:CACM:69}, so that all operations can be performed in
constant time (To be precise, the \emph{DeleteMin}-operation also
needs a pointer that runs through all different $g$-values once;
however, we can neglect this in comparison to the
number of expansions). To expand a vertex, at most $2^k-1$
successor vertices have to be generated, since we have the choice
of introducing a gap in each sequence. Thus, Dijkstra's algorithm can
solve the multiple sequence alignment problem in $O(2^kN^k)$ time and
$O(N^k)$ space for $k$ sequences of length $\leq N$.


A means to reduce the number of nodes that have to be stored 
for path reconstruction is by associating a counter with each node
that maintains the number of children whose backtrack-pointer refers 
to them \cite{gupta96}. Since each node can be expanded at most once, 
after this the number of referring backtrack-pointers can only 
decrease, namely, whenever a cheaper path to one of its children is
found. If a node's reference count goes to zero, whether immediately 
after its expansion or when it later loses a child, it can be deleted
for good. This way, we only keep 
nodes in memory that have at least one descendant currently in the
priority queue. Moreover, auxiliary data structures for vertices and
coordinates are most efficiently stored in tries (prefix trees); they
can be equipped with reference counters as well and be freed
accordingly when no longer used by any edge. 


The same complexity as for Dijkstra's algorithm holds for
\emph{dynamic programming (DP)}; it 
differs from the former one in that it scans the nodes in a
fixed order that is known beforehand (hence, contrary to the name the
exploration scheme is actually static). The exact order of the scan can
vary  (e.g., row-wise or column-wise), as long as it is compatible
with the topological ordering 
of the graph (e.g., for two sequences that the cells left, top, and
diagonally top-left have been explored prior to a cell). One
particular such ordering is that of \emph{antidiagonals}, diagonals
running from upper-right to lower-left. The calculation of the
antidiagonal of a node merely amounts to summing up its $k$ coordinates.

Hirschberg (1975\nocite{hirsch}) noticed that in order to determine only
the cost of 
the optimal alignment $g^*(t)$, it would not be necessary to store the
whole matrix; 
instead, when proceeding e.g. by rows it suffices to keep track
of only $k$ of them at a time, deleting each row as soon as the next
one is completed. This reduces the space requirement by one dimension
from $O(N^k)$ to $O(k N^{k-1})$. In order to recover the solution path
at the end, re-computation 
of the lost cell values is needed. A
\emph{Divide-and-conquer}-strategy applies the algorithm twice to
half the grid each, once in forward and once in backward direction,
meeting at a fixed middle row. By adding the corresponding  
forward and backward distances in this middle row and finding the
minimum, one cell 
lying on an optimal path can be recovered. This cell essentially
splits the problem into two smaller subproblems, one from the upper
left corner to it, and the other one to the lower right corner; they
can be recursively solved using the same method. In two dimensions,
the computation time is at most doubled, and the overhead reduces even
more in higher dimensions.

The \emph{FastLSA} algorithm \cite{davidson01} further refines Hirschberg's
algorithm by exploiting additionally available memory to store more than
one node on an optimal path, thereby reducing the number of
re-computations. 

\subsection{Algorithms Utilizing Bounds}

While Dijkstra's algorithm and dynamic programming can be viewed as
variants of \emph{breadth-first search}, we achieve \emph{best first
  search} if we expand nodes $v$ in the order of an estimate (lower bound)
of the total cost of a path from $s$ to the $t$ passing
through $v$. Rather than using the $g$-value as in 
Dijkstra's algorithm, we use $f(v) := g(v)+h(v)$ as the heap key, where $h(v)$ is
a lower bound on the cost of an optimal path from $v$ to $t$.
If $h$ is indeed \emph{admissible}, then the first
solution found is guaranteed to be optimal \cite{Nilsson}. This is the
classical best-first search algorithm, the \ast\ \emph{algorithm},
well known in the artificial intelligence 
community. In this context, the priority queue maintaining the
generated nodes is often also called
the \emph{Open} list, while the nodes that have already been expanded and 
removed from it constitute the
\emph{Closed} list. Fig.~\ref{fig:layers} schematically depicts a
snapshot during a two-dimensional alignment problem,
where all nodes with $f$-value no larger than the current $f_{\min}$
have been expanded. Since the accuracy of the heuristic decreases with
the distance to the goal, the typical `onion-shaped' distribution
results, with the bulk being located closer to the start node, and
tapering out towards higher levels. 
\begin{figure}
\begin{center}
%\epsfxsize8cm
%\epsfysize7cm
\includegraphics[width=0.5\columnwidth]{layers}
\end{center}
\caption{Snapshot during best-first search in pairwise alignment (schematically).\label{fig:layers}}
\end{figure}

The \ast\ algorithm can significantly  reduce the total number of
expanded and generated nodes; therefore, in higher dimensions it is
clearly superior to dynamic programming. However, in contrast to the
Hirschberg algorithm, it still stores all of the explored nodes in the
\emph{Closed} list. Apart from keeping track of the solution path,
this is necessary to prevent the search from ``leaking back'', in the
following sense.

A heuristic $h$ is called \emph{consistent} if $h(x) \leq h(x') +
c(x,x')$, for any node $x$ and its child $x'$. A consistent heuristic
ensures that (as in the case of Dijkstra's algorithm) at the time a
node is expanded, its $g$-value is optimal, and hence it is never expanded
again. However, if we try to delete the \emph{Closed} nodes, then 
there can be topologically smaller nodes in \emph{Open} with a higher
$f$-value; when those are expanded at a later stage, they can lead to
the re-generation of the node at a non-optimal $g$-value, since the first
instantiation is no longer available for duplicate
checking. In Fig.~\ref{fig:layers}, nodes that might be subject to
spurious re-expansion are marked ``X''.


Researchers have tried to avoid these leaks, while retaining the
basic \ast\ search scheme. Korf proposed to store a list of forbidden
operators with each node, or to place the parents of a deleted node on
\emph{Open} with $f$-value infinity \cite{korf-dcbds,korf-dcfs}.   
However, as \fullciteA{zhou03} remark, it is hard to combine this
algorithm with techniques for reduction of the \emph{Open} list, and moreover
the storage of operators lets the size of the nodes grow exponentially
with the number of sequences. In their algorithm, they keep track of
the \emph{kernel} of the \emph{Closed} list, which is defined as the set of nodes
that have only \emph{Closed} nodes as parents; otherwise a \emph{Closed} node is
said to be in the \emph{boundary}. The key idea is that only the
boundary nodes have to be maintained, since they shield the kernel
from re-expansions. Only when the algorithm gets close to the memory
limit nodes from the kernel are deleted; the backtrack pointer of the
children is changed to the parents of the deleted nodes, which become
\emph{relay nodes} for them. For the final reconstruction of the
optimal solution path, the algorithm is called recursively for each
relay node to bridge the gap of missing edges.

In addition to the \emph{Closed} list, also the \emph{Open} list can grow
rapidly in sequence alignment problems. Particularly, since in the
original \ast\ algorithm the expansion of a node  
generates all of its children at once, those whose $f$-value is
larger than the optimal cost $g^*(t)$ are kept in the heap up to the
end, and waste much of the available space. 

If an upper bound $U$ on the
optimal solution cost $g^*(t)$ is known, then nodes $v$ with $f(v) > U$
can be pruned right away;  this idea is used in several articles
%(Spouge, 1989; Gupta et al., 1996).
\cite{spouge,gupta96}. 
%\nocite{spouge,gupta96}
One of the most successful approaches is
Yoshizumi et al.'s (2000\nocite{PE}) 
%\fullciteA{PE}'s
\emph{Partial Expansion}
\ast\  (\pea). Each node stores an additional value 
$F$, which is the minimum $f$-value of all of its yet ungenerated
children. In each step, only a node with minimum $F$-value is
expanded, and only 
those children with $f=F$ are generated. This algorithm clearly only
generates nodes with $f$ value no larger than the optimal cost, which
cannot be avoided altogether. However, the overhead in computation
time is considerable: in the straightforward implementation, if we
want to maintain nodes of constant size, generating one edge requires
determining the $f$-values of all successors, such that for an
interior node which eventually will be fully expanded the computation
time is of the order of the square of the number of successors, which
grows as $O(2^k)$ with the number of sequences $k$. As a remedy, in
the paper it is proposed to relax the condition by generating all
children with $f \leq F + C$, for some small $C$.

An alternative general search strategy to \ast\ that uses only linear space is
\emph{iterative deepening} \ast (\ida) \cite{korf_ida}. The basic
algorithm conducts a depth-first search up to a pre-determined
threshold for the $f$-value. During the search, it keeps track of the
smallest $f$-value of a generated successor that is larger than the
threshold. If no solution is found, this provides an increased
threshold to be used in the next search iteration.


\fullciteA{wah} suggested more liberal schemes for
determining the next threshold dynamically in order to minimize the
number of recomputations. \ida\ is most
efficient in tree structured search spaces. However, it is
difficult to detect duplicate expansions without additional memory;
Therefore, unfortunately  it is not applicable in 
lattice-structured graphs like in the sequence alignment problem due to the
combinatorially explosive number of paths between any two given nodes.
 
A different line of research tries to restrict the search space of the
breadth-first approaches by incorporating bounds. \citeA{ukkonen}
presented an algorithm for the pairwise alignment problem which is
particularly efficient for similar sequences; its computation time
scales as $O(d m)$, where $d$ is the optimal solution cost. First
consider the problem of deciding whether a solution exists whose cost
is less than some upper threshold $U$. We can restrict the evaluation of the
\DP\ matrix to a band of diagonals where the minimum number of
indels required to reach the diagonal, times the minimum indel cost,
does not exceed $U$. In general, starting with a minimum $U$ value, we
can successively double $G$ until the test returns a solution; the
increase of computation time due to the recomputations is then also
bounded by  a factor of 2.

Another  approach for multiple sequence alignment is to make use of
the lower bounds $h$ from $A^*$. The key idea is the following: Since
all nodes with an $f$-value lower than $g^*(t)$  have
to be expanded anyway in order to
guarantee optimality,  we might as well explore them in any reasonable
order, like that of Dijkstra's algorithm or \DP, if we only knew the
optimal cost. Even slightly higher upper bounds will still help
pruning. \citeA{spouge} proposed to bound \DP\ to 
vertices $v$ where $g(v)+h(v)$ is smaller than an upper bound for~$g^*(t)$.

\emph{Linear Bounded Diagonal Alignment (LBD-Align)} \cite{davidson01}
uses an upper bound in order to reduce the computation time and memory
in solving a pairwise alignment problem by dynamic programming. The algorithm
calculates the \DP\ matrix one antidiagonal at a time,
starting in the top left corner, and working down
towards bottom-right. While \ast\ would have to check the bound in
every expansion, \emph{LBD-Align} only checks the top and
bottom cell of each diagonal. If e.g.\ the top cell of a
diagonal has been pruned, all the remaining cells in that row can be
pruned as well, since they are only reachable through it; this means
that the pruning frontier on the next row can be shifted down by
one. Thus, the pruning overhead can be reduced from a quadratic to a 
linear amount in terms of the sequence length.



\subsection{Obtaining Heuristic Bounds}

Up to now we have assumed lower and upper bounds, without specifying
how to derive them. Obtaining an inaccurate \emph{upper bound} on $g^*(t)$ is
fairly easy, since we can use the cost of \emph{any} valid path
through the lattice. Better estimates are e.g.\ available from
heuristic linear-time alignment programs such as \emph{FASTA} and
\emph{BLAST} \cite{altblast}, 
which are a standard method for database searches. Davidson (2001) \nocite{davidson01}
employed a local beam search scheme. 

\citeA{gusfield93star} proposed an approximation
called the \emph{star-alignment}. Out of all the sequences to be
aligned, one \emph{consensus sequence} is chosen such that the sum of its
pairwise alignment costs to the rest of the sequences is minimal. Using
this ``best'' sequence as the center, the other ones are aligned using
the ``once a gap, always a gap'' rule. Gusfield showed that the cost
of the optimal alignment is greater or equal to the cost of this star
alignment, divided by $(2-2/k)$.

For use in heuristic estimates, lower bounds on the $k$-alignment are
often based on optimal  
alignments of subsets of $m<k$ sequences. In general, for a vertex
$v$ in $k$-space, we are looking for a lower bound for a path from $v$
to the target corner $t$. Consider first the case $m=2$. The cost of such
a path is, by definition, the sum of its edge costs, where each edge
cost in turn is the sum of all pairwise (replacement or gap)
penalties. Each multiple sequence alignment induces a pairwise
alignment for sequences $i$ and $j$, by simply copying rows $i$ and
$j$ and ignoring columns with a ``{\tt \_}'' in both rows. These
pairwise alignments can be visualized as the projection of an
alignment onto its \emph{faces}, cf.~Fig.~\ref{fig:cube}.

By interchange of the summation order, the sum-of-pairs cost is the sum
of all pairwise alignment costs of the respective paths projected on
a face, each of which cannot be smaller than the optimal pairwise
path cost. Thus, we can construct an admissible heuristic $h_{pair}$
by computing, 
for each pairwise alignment and for each cell in a pairwise problem,
the cheapest path cost to the goal node. 

The optimal solutions  to all pairwise alignment problems needed for
the lower bound $h$ values 
are usually computed prior to the main search in a preprocessing step
\cite{Ikeda94algorithms}. To this end, it suffices to apply the
ordinary \DP\ procedure; however, since this time we are 
interested in the lowest cost of a path from $v$ to $t$, it runs in
backward direction, proceeding from the lower right corner to
the upper left, expanding all possible parents of a vertex in each step.

Let $U$ be an upper bound on the cost of an optimal multiple sequence
alignment $G$. The sum of all optimal alignment costs $L_{ij}=d(s_{ij},t_{ij})$ for
 pairwise subproblems $i,j \in \{1, \dots, k\}, i < j$, call it
$L$, is a lower bound on $G$. 
\fullciteA{carr} pointed out that by the additivity of the sum-of-pairs
cost function, any pairwise alignment induced by the optimal multiple
sequence alignment can at most be $\delta = U-L$ larger than the
respective optimal pairwise alignment. This bound can be used to
restrict the number of values that have to be computed in the
preprocessing stage and have to be stored for the calculation of the
heuristic: for the pair of sequences $i,j$, only those nodes $v$ are
feasible such that a path from the start node $s_{i,j}$ to the goal
node $t_{i,j}$ exists with total cost no more than $L_{i,j} +
\delta$. To optimize the storage requirements, we can combine the
results of two searches. First, a forward pass determines for each
relevant node $v$ the minimum distance $d(s_{ij},v)$ from the start node. The
subsequent backward pass uses this distance like an 'exact heuristic'
and stores the distance $d(v,t_{ij})$ from the target node only for those
nodes with $d(s_{ij},v) + d(v,t_{ij}) \leq d(s,t) +
\delta$\footnote{A slight technical complication arises
  for affine gap costs: recall that \DP\
  implementations usually charge the gap opening penalty
  to the $g$-value of the edge $e$ starting the gap, while the 
  edge $e'$ ending the gap carries no extra penalty at all. However,
  since the sum of pairs heuristics $h$ is 
  computed in backward direction, using the same algorithm we would assign the
  penalty for the same path instead to $e'$. This 
  means that the heuristic $f=g+h$ would no longer be guaranteed to be
a lower bound, since it contains the penalty twice. As a remedy, it is
necessary to make the computation symmetric by charging both the
beginning and end of a gap with half the cost each. The case of the
beginning and end of the sequences can be handled most conveniently by
starting the search from a ``dummy'' diagonal edge $((-1, \dots, -1),
(0, \dots, 0))$, and defining the target edge to be the dummy diagonal
edge $((N,\dots, N), (N+1, \dots, N+1))$, similar to the arrows shown in
Fig.~\ref{fig:cube}.}.

Still, for larger alignment problems the required storage size can be
extensive. The program \emph{MSA} \cite{gupta96} allows the user to adjust
$\delta$ to values below the Carrillo-Lipman bound individually for
each pair of sequences. This makes it
possible to generate at least heuristic alignments if time or memory
doesn't allow for the complete solution; moreover, it can be recorded
during the search if the $\delta$-bound was actually reached. In the
negative case, optimality of the found solution is still guaranteed;
otherwise, the user can try to run the program again with slightly
increased bounds. 

The general idea of precomputing simplified problems and storing the
solutions for use as a heuristic has been explored under the name of
\emph{pattern databases} \cite{PatternDatabases}. However, these
approaches implicitly assume that the computational cost can be
amortized over many search instances to the same target. In contrast,
in the case of MSA, the heuristics are instance-specific, so that we
have to strike a balance. We will discuss this in greater depth in
Sec.~\ref{sec:tradeoff}. 
 
%\goodbreak
\section{Iterative-Deepening Dynamic Programming\label{sec:iddp}}

As we have seen, a fixed search order as in dynamic programming can
have several advantages over pure best-first selection.

\begin{itemize}
\item Since \emph{Closed} nodes can never be reached more than once
  during the search, it is safe to delete useless ones (those that
  are not   part   of any shortest path to the current \emph{Open} nodes) and to apply
  path compression schemes, such as the Hirschberg algorithm. No
  sophisticated schemes for avoiding 'back leaks' are required, such as
  the above-mentioned methods of core set maintenance and dummy node insertion
  into \emph{Open}. 
\item Besides the size of the \emph{Closed} list, the memory requirement of
  the \emph{Open} list is determined by the \emph{maximum} number of
  nodes that are open \emph{simultaneously at any time} while the algorithm
  is running. When the $f$-value is used as the key for the priority  
  queue, the \emph{Open} list usually contains all nodes with $f$-values in some
  range $(f_{min}, f_{min}+\delta)$; this set of nodes is generally spread
  across all over the search space, since $g$ (and accordingly $h=
  (f-g)$) can vary arbitrarily between $0$ and $f_{min}+\delta$. As
  opposed to that, if \DP\ 
  proceeds along levels of antidiagonals or rows, at any iteration at
  most $k$   levels have to be maintained at the same time, and hence
  the size of the \emph{Open} list can be controlled more effectively.
  In Fig.~\ref{fig:layers}, the pairwise alignment is partitioned
  into antidiagonals: the maximum number of open nodes in any two
  adjacent levels is four, while the total amounts to
  seventeen\footnote{Contrary to what the figure might suggest, \ast\
  can open more than two nodes per level in pairwise alignments, if
  the set of nodes no worse than some $f_{\min}$ contains ``holes''.}.

\item For practical purposes, the running time should not only be
  measured in terms of the number of node expansions, but one
  should also take into account the \emph{execution time} needed for
  an expansion. By arranging 
  the exploration order such that edges with the same head node (or
  more generally, those sharing a common coordinate prefix) are
  dealt with one after the other, much of the computation can be
  cached, and edge generation can be sped up significantly. We will
  come back to this point in Sec.~\ref{sec:ndimheu}.
\end{itemize}

The remaining issue of a static exploration scheme consists in
adequately bounding the search space using the $h$-values. \ast\ is
known to be minimal in 
terms of the number of node expansions. If we knew the cost $g^*(t)$ of a
cheapest solution path beforehand, we could simply proceed level by level
of the grid, however only immediately prune generated edges $e$ whenever
$f(e)>g^*(t)$. This would ensure that we only generate those edges that would
have been generated by algorithm $A^*$, as well. An upper threshold would
additionally help reduce the size of the \emph{Closed} list, since a
node  can be pruned if all of its children lie beyond the threshold;
additionally, if this node is the only child of its parent, this can
give rise to a propagating chain of ancestor deletions.

We propose to apply a search scheme that carries out a
series of searches with successively larger thresholds, until a
solution is found (or we run out of memory or patience). The use of
such an upper bound parallels that in the \ida\ algorithm.

The resulting algorithm, which we will refer to as \emph{Iterative-Deepening
Dynamic Programming (IDDP)}, is sketched in Fig.~\ref{alg:iddp}. The
outer loop initializes the threshold with a lower bound (e.g.,
$h(s)$), and, unless a solution is found, increases it up to an upper
bound. In the same
manner as in the  \ida\ algorithm, in order to make sure that
at least one additional edge is explored in each iteration the
threshold has to be increased correspondingly at least to the minimum
cost of a fringe edge that exceeded the previous threshold. This
fringe increment is maintained in the variable \emph{minNextThresh},
initially estimated as the upper bound, and repeatedly decreased in
the course of  the following expansions. 
\begin{figure}
{\small
\begin{algorithmic}
\STATE{{\bf procedure} IDDP(Edge startEdge, Edge targetEdge, int lowerBound, int upperBound) }
\STATE{int thresh = lowerBound}

\COMMENT{Outer loop: Iterative deepening phases}
\WHILE{(thresh $\leq$ upperBound)}
\STATE{Heap h = $\{$(startEdge, 0)$\}$}
\STATE{int minNextThresh = upperBound}
\STATE{}\COMMENT{Inner loop: Bounded dynamic programming}
\WHILE{({\bf not} h.IsEmpty())}
\STATE{Edge e = h.DeleteMin()}
\COMMENT{Find and remove an edge with minimum level}
\IF{(e == targetEdge)}
\STATE{}\COMMENT{Optimal alignment found}
\STATE{{\bf return} TraceBackPath(startEdge, targetEdge)}
%\STATE{}\COMMENT{Update current best solution}
%\STATE{print(``Alignment found'')}
%\STATE{{\bf return} GetPath(e)}
%\STATE{bestSolution = e}
%\STATE{{\bf continue}}
\ENDIF
\STATE{Expand(e, thresh, minNextThresh)}
\ENDWHILE
\STATE{int threshIncr = ComputeThreshIncr()}
\COMMENT{Compute search threshold for next iteration, see text}
\STATE{ thresh = max(thresh + threshIncr, minNextThresh)}
\ENDWHILE
\STATE{ print(``No alignment with cost at most upperBound found'')}
\end{algorithmic}
}
\caption{Algorithm \emph{Iterative-Deepening Dynamic Programming.}\label{alg:iddp}}
\end{figure}


\begin{figure}[tbh]
{\small
\begin{algorithmic}
\STATE{{\bf procedure} Expand(Edge e, int thresh, int minNextThresh) }
\FORALL{Edge child $\in$ Succ(e)}
\STATE{}
\COMMENT{Retrieve child or tentatively generate it if not yet
  existing, set boolean variable `created' accordingly}
%% \IF{(child.GetH() == {\bf nil})}
%% \STATE{ print(``Delta-bound on faces reached, final solution might not be optimal'')}
%% \STATE{{\bf continue}}
%% \ENDIF
\STATE{int newG = e.GetG() + GapCost(e, child)}
\STATE{\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \  +  child.GetCost()}
\STATE{int newF = newG + child.GetH()}
\IF{(newF $\leq$ thresh {\bf and} newG $<$ child.GetG())}
\STATE{}\COMMENT{Shorter path than current best
  found, estimate within threshold} 
\STATE{child.SetG(newG)}
\STATE{UpdateEdge(e, child, h)}
\COMMENT{Update search structures}
\ELSIF{(newF $>$ thresh)}
\STATE{minNextThresh =}
\STATE{\ \ \ min(minNextThresh, newF)} 
\STATE{}\COMMENT{Record minimum of pruned edges}
\IF{(created)}
\STATE{Delete(child)}
\COMMENT{Make sure only promising edges are stored}
\ENDIF
\ENDIF
\ENDFOR
\IF{(e.ref == 0)}
\STATE{DeleteRec(e)}
\COMMENT{No promising children could be inserted into the heap}
\ENDIF
\end{algorithmic}
}
\caption{Edge expansion in \iddp.\label{alg:expand}}
\end{figure}

In each step of the inner loop, we select and remove a node
from the priority queue whose level is minimal. As explained later in
Sec.~\ref{sec:ndimheu}, it is favorable to break ties according to the
lexicographic order of target nodes. Since the total number of
possible levels is comparatively small and known in advance, the
priority queue can be implemented using an array of linked lists
\cite{Dial:CACM:69}; this provides constant time operations for
insertion and deletion.

%% In contrast to \ast, it is not admissible to terminate the inner loop with the first  
%% solution found, since there might still be cheaper paths left in the
%% remaining incompletely explored levels. Instead, the algorithm always
%% runs until the heap becomes empty, and then returns the best found
%% solution recorded in variable \emph{bestSolution}.

The expansion of an edge $e$ is partial (Fig.~\ref{alg:expand}). A
child edge might already 
exist from an earlier expansion of an edge with the same head vertex;
we have to test if we can decrease the $g$-value. Otherwise, we
generate a new edge, if only temporarily for the sake of calculating its
$f$-value; 
that is, if its $f$-value exceeds the search threshold of the current
iteration, its memory is immediately reclaimed. Moreover, in this case
the fringe threshold \emph{minNextThresh} is updated. In a practical
implementation, we can prune unnecessary accesses to partial alignments
\emph{inside} the calculation of the heuristic \emph{e.GetH()} as soon as
as the search threshold has already been reached. 

The relaxation of a child edge within the threshold is performed by 
the subprocedure \emph{UpdateEdge} (cf.~Fig.~\ref{alg:update}). This
is similar to the corresponding relaxation step in $A^*$, updating the
child's $g$- and $f$ values, its parent pointers, and inserting it
into \emph{Open}, if not already contained. However, in contrast to 
  best-first search, it is inserted into the heap 
according to the antidiagonal level of its head vertex. Note that in
the event that the former parent loses its last child, propagation of
deletions   
(Fig.~\ref{alg:DeleteRec}) can ensure that only those \emph{Closed} 
nodes continue to be stored that belong to some solution
path. Edge deletions can also ensue deletion of dependent vertex and
coordinate data structures (not shown in the pseudocode). The other
situation that gives rise to deletions is if immediately after the
expansion of a node no children are pointing back to it (the children 
might either be reachable more cheaply from different nodes, or their
$f$-value might exceed the threshold). 

\begin{figure}[tbh]
{\small
\begin{algorithmic}
\STATE{ {\bf procedure} UpdateEdge(Edge parent, Edge child, Heap h)}
%\STATE{ h.DeleteIfContains(child,}
%\STATE{\ \ \  \mbox{child.GetHead().GetLevel()})}
%\COMMENT{DecreaseKey}
\STATE{parent.ref$++$}
\STATE{child.GetBacktrack().ref$--$}
\IF{(child.GetBacktrack().ref == 0)}
\STATE{DeleteRec(child.GetBacktrack())}
\COMMENT{The former parent has lost its last child and becomes useless}
\ENDIF
\STATE{child.SetBacktrack(parent)}
\IF{({\bf not} h.Contains(child))}
\STATE{h.Insert(child, child.GetHead().GetLevel())}
\ENDIF
\end{algorithmic}
}
\caption{Edge relaxation in \iddp.\label{alg:update}}
\end{figure}
\begin{figure}
{\small
\begin{algorithmic}
\STATE{ {\bf procedure} DeleteRec(Edge e)}
\IF{(e.GetBacktrack() $\neq$ {\bf nil})}
\STATE{ e.GetBacktrack().ref$--$}
\IF {(e.GetBacktrack().ref == 0)}
\STATE{DeleteRec(e.GetBacktrack())}
\ENDIF
\ENDIF
\STATE{Delete(e)}
\end{algorithmic}
}
\caption{Recursive deletion of edges that are no longer part of any solution path.\label{alg:DeleteRec}}
\end{figure}
\begin{figure}[th]
{\small
\begin{algorithmic}
\STATE{ {\bf procedure} TraceBack(Edge startEdge, Edge e)}
\IF{(e == startEdge)}
\STATE{{\bf return}}
\COMMENT{End of recursion}
\ENDIF
\IF{(e.GetBackTrack().GetTarget() $\neq$ e.GetSource())}
\STATE{}\COMMENT{Relay node: recursive path reconstruction}
\STATE{IDDP( e.GetBackTrack(), e, e.GetF(), e.GetF())}
\ENDIF
\STATE{OutputEdge(e)}
\STATE{TraceBack(startEdge, e.GetBackTrack())}
\end{algorithmic}
}
\caption{Divide-and-Conquer solution reconstruction in reverse order.\label{alg:traceBack}}
\end{figure}


The correctness of the algorithm can be shown analogously to the
soundness proof of \ast. If the threshold is smaller than $g^*(t)$, the
\DP\ search will terminate without encountering a solution; otherwise,
only nodes are pruned that cannot be part of an optimal path. The
invariant holds that there is always a node in each level which 
lies on an optimal path and is in the \emph{Open} list. Therefore, if the
algorithm terminates only when the heap runs empty, the best found
solution will indeed be optimal. 


The iterative deepening strategy results in an overhead computation time
due to re-expansions, and we are trying to restrict this overhead as
much as possible. More precisely, we want to minimize the ratio
$$\nu = \frac{n_{\mbox{\emph{\scriptsize IDDP}}}}{n_{A^*}},$$
where $n_{\mbox{\emph{\scriptsize IDDP}}}$ and $n_{A^*}$ denote the number of expansions in
\emph{IDDP} and $A^*$, respectively. One way to do so \cite{wah} is to
choose a threshold sequence $\theta_1, \theta_2, \dots$ such that the
number of expansions $n_i$ in stage $i$ satisfies 
$$n_i = r n_{i-1},$$ for some fixed ratio $r$. If we choose $r$ too
small, the number of re-expansions and hence the computation time will
grow rapidly, if we choose it too big, then the threshold of the last
iteration can exceed the optimal solution cost significantly, and we
will explore many irrelevant edges.  Suppose that $n_0 r^p < n_{A^*} \leq 
n_0 r^{p+1}$. Then the algorithm performs $p+1$ iterations. In the worst
case, the overshoot will be maximal if $A^*$ finds the optimal
solution just above the previous threshold, $n_{A^*} = n_0 r^p +
1$. The total number of expansions is $n_0 \sum_{i=0}^{p+1} r^i =
n_0 \frac{r(r^{p+1}-1)}{r-1}$, and the ratio $\nu$ becomes approximately
$\frac{r^2}{r-1}$. By setting the derivative of this expression to
zero, we find that the optimal value for $r$ is 2; the number of
expansions should double from one search stage to the next. If we
achieve doubling, we will expand at most four times as many nodes as
$A^*$.   

Like in Wah and Shang's \citeyear{wah} scheme, we dynamically adjust the threshold using runtime
information. Procedure \emph{ComputeThreshIncr} stores the sequence of
expansion numbers and thresholds from the previous search stages, and
then uses curve fitting for extrapolation (in the first few
iterations without sufficient data available, a very small
default threshold is applied). We found 
that the distribution of nodes $n(\theta)$ with $f$-value smaller or
equal to threshold $\theta$ can be modeled 
very accurately according to the exponential approach
$$n(\theta) = A \cdot B^{\theta}.$$
Consequently, in order to attempt to double the number of expansions, we choose
the next threshold according to $$\theta_{i+1}=\theta_i+
\frac{1}{\log_{2}B}.$$ 

%\goodbreak
\section{Sparse Representation of Solution Paths\label{sec:sparse}}

When the search progresses along antidiagonals, we do not have to fear
back leaks, and are free to prune \emph{Closed} nodes. Similarly as in Zhou and Hansen's \citeyear{zhou03} work, however, we only want to delete them lazily and incrementally
when being forced by the algorithm approaching the computer's
memory limit.  

When deleting an edge $e$, the backtrack-pointers of its child edges
that refer to it are redirected to the respective predecessor of $e$,
whose reference count is increased accordingly. In the resulting
\emph{sparse solution path} representation, backtrack pointers can
point to any optimal ancestors.

After termination of the main search, we trace back the pointers
starting with the goal edge; this is outlined in Procedure
\emph{TraceBack} (Fig.~\ref{alg:traceBack}), which prints out the
solution path in reverse order. Whenever an edge $e$ points back to an
ancestor $e'$ which is not its direct parent, we apply an
auxiliary search from start edge $e'$ to goal edge $e$ in order to
reconstruct the missing links of the optimal solution path. 
The search threshold can now be fixed at the known solution cost;
moreover, the auxiliary search can prune those edges that cannot be
ancestors of $e$ because they have some coordinate greater than the
corresponding coordinate in $e$. Since also the shortest distance 
between $e$ and $e'$ is known, we can stop at the first path that is
found at this cost. To improve the efficiency of the 
auxiliary search even further, the heuristic could be recomputed to
suit the new target. Therefore, the cost of restoring the solution path is usually marginal compared to that of the main search.

Which edges are we going to prune, in which order? For simplicity,
assume for the moment that the \emph{Closed} list consists of a single
solution path. According to the Hirschberg approach, we would keep only one
edge, preferably lying near the center of the search space (e.g., on
the longest anti-diagonal), in order to
minimize the complexity of the two auxiliary searches. With additional
available space allowing to store three relay edges, we would
divide the  search space into four subspaces of about equal size (e.g.,
additionally storing the antidiagonals half-way between the middle
antidiagonal and the start node resp.\ the target node). By extension,
in order to incrementally save space under diminishing resources we would
first keep only every other level, then every fourth, and so on, until 
only the start edge, the target edge, and one edge half-way on the
path would be  left. 

%We can decrease the density in $\lfloor \log_2 k
%N\rfloor$ steps. 

Since in general the \emph{Closed} list contains multiple solution paths
(more precisely, a tree of solution paths),
we would like to have about the same density of relay edges on each of
them. For the case of $k$ sequences, an edge reaching level
$l$  with its head node can originate with its tail node from level
$l-1, \dots, l-k$. Thus, not every solution path passes through each
level, and deleting every other level could result in leaving one path
completely intact, while extinguishing another totally. Thus, it is
better to consider contiguous \emph{bands} of $k$ levels each, instead
of individual levels. Bands of this size cannot be skipped by any path.
The total number of antidiagonals in an alignment problem of $k$
sequences of length $N$ is $k\cdot N - 1$; thus, we can decrease the
density in $\lfloor \log_2 N \rfloor$ steps. 

A technical implementation issue concerns the ability to enumerate all edges
that reference some given prunable edge, without explicitly storing
them in a list. However, the reference counting method described above
ensures that any \emph{Closed} edge can be reached by following a path
bottom-up from some edge in \emph{Open}. The procedure is sketched in
Fig.~\ref{alg:SparsifyClosed}. The variable \emph{sparse} denotes the
interval between level bands that are to be maintained in memory. In the
inner loop, all paths to \emph{Open} nodes are traversed in backward
direction; for each edge $e'$ that falls into a prunable band, the pointer
of the successor $e$ on the path is redirected to its respective backtrack
pointer. If $e$ was the last edge referencing $e'$, the latter one is
deleted, and the path traversal continues up to the start edge. When all \emph{Open} nodes have been
visited and the memory bound is still exceeded, the outer loop tries
to double the number of prunable bands by increasing \emph{sparse}.



\begin{figure}[tb]
{\small
\begin{algorithmic}
\STATE{{\bf procedure} SparsifyClosed()}
\FOR{(int sparse = 1 {\bf to} $\lfloor \log_2 N \rfloor$)}
%\STATE{}\COMMENT{Increase interval between stored level bands}
\WHILE{(UsedMemory() $>$ maxMemory {\bf and} {\bf exists} $\{$Edge e $\in$ \emph{Open} $|$ e.GetLastSparse() $<$ sparse$\})$}
\STATE{Edge pred = e.GetBacktrack()}
\STATE{}\COMMENT{Trace back solution path}
\WHILE{(pred $\neq$ {\bf nil} {\bf and} e.GetLastSparse() $<$ sparse)}
\STATE{e.SetLastSparse(sparse)}
\COMMENT{Mark to avoid repeated trace-back}
\IF{($\lfloor$pred.GetHead().GetLevel() / $k$$\rfloor$ {\bf mod}
  $2^{\mbox{\scriptsize sparse}}$ $\neq$ 0)}
\STATE{}\COMMENT{pred lies in prunable band: redirect pointer}
\STATE{e.SetBacktrack(pred.GetBacktrack())}
\STATE{e.GetBacktrack().ref$++$}
\STATE{pred.ref$--$}
\IF{(pred.ref == 0)}
\STATE{}\COMMENT{e is the last remaining edge referring to pred}
\STATE{DeleteRec(pred)}
\ENDIF
\ELSE
\STATE{}\COMMENT{Not in prunable band: continue traversal}
\STATE{e = e.GetBacktrack()}
\ENDIF
\STATE{pred = e.GetBacktrack()}
\ENDWHILE
\ENDWHILE
\ENDFOR
\end{algorithmic}
}
\caption{Sparsification of \emph{Closed} list under restricted memory.\label{alg:SparsifyClosed}}
\end{figure}


Procedure \emph{SparsifyClosed} is called regularly during the search,
e.g., after each expansion. However, a naive version as described
above would incur a huge overhead in computation time, particularly when the algorithm's memory consumption is close to the 
limit. Therefore, some optimizations are necessary. First, we avoid
tracing back the same solution path at the same (or lower)
\emph{sparse} interval by recording for each edge the interval when it
was traversed the last time (initially zero); only for an increased
variable \emph{sparse} there can be anything left for further
pruning. In the worst case, each edge will be inspected  $\lfloor \log_2
 N \rfloor$ times. 
Secondly, it would be very inefficient to actually inspect each \emph{Open}
node in the inner loop, just to find that its solution path has been
traversed previously, at the same or higher \emph{sparse} value;
however, with an appropriate bookkeeping strategy it is possible to
reduce the time for this search overhead to~$O(k)$.

%\goodbreak


\section{Use of Improved Heuristics\label{sec:ndimheu}}


As we have seen, the estimator $h_{pair}$, the sum of optimal
pairwise goal distances, gives a lower bound on the actual path
length. However, more powerful heuristics are also conceivable. 
While their computation will require more resources, the trade-off 
can prove itself worthwhile; the tighter the estimator is, the smaller 
is the space that the main search needs to explore.


\subsection{Beyond Pairwise Alignments}

\citeA{kobayashi-improvement} suggested to generalize  $h_{pair}$ by considering 
optimal solutions for subproblems of size $m>2$. They proved that the following heuristics are
admissible and more informed than the pairwise estimate. 
\begin{itemize}
\item $h_{all,m}$
is the sum of all $m$-dimensional optimal costs, divided by $ {k-2}
\choose {m-2}$. 
\item $h_{one, m}$ splits the sequences into two sets of sizes $m$
  and $k-m$; the heuristic is the sum of the optimal cost of the
first subset, plus that of the second one, plus the sum of all
2-dimensional optimal costs of all pairs of sequences in different
  subsets. Usually, $m$ is chosen close to $k/2$. 
\end{itemize}

These improved heuristics can reduce the main search effort by orders
of magnitudes. However, in contrast to pairwise sub-alignments, time
and space resources devoted to compute and store higher-dimensional
heuristics are in general no longer negligible compared to the main
search. 
\citeA{kobayashi-improvement} noticed that even for the case $m=3$ of triples of
sequences, it can be impractical
to compute the entire subheuristic $h_{all,m}$. As
one reduction, they show that it suffices to restrict oneself to
nodes where the path cost does not exceed the optimal path cost of the
subproblem by more than 

$$ \delta = {{k - 2} \choose { m - 2}} U - \sum_{i_1, \dots, i_{m}} d(s_{i_1, \dots, i_{m}}, t_{i_1, \dots, i_{m}});$$

this threshold can be seen as a generalization of the Carrillo-Lipman
bound. However, it can still incur excessive overhead in space and
computation time for the computation of the $k \choose m$ lower-dimensional
subproblems.  A drawback is that it
requires an upper bound $U$, on whose accuracy also the algorithm's
efficiency hinges. We could improve this bound by applying more
sophisticated heuristic methods, but it seems counterintuitive to spend
more time doing so which we would rather use to calculate the exact solution. 
In spite of its advantages for the main search, the expensiveness of
the heuristic calculation appears as a major obstacle.

\fullciteA{octtree} suggested to partition the  
heuristic into (hyper-) cubes using a hierarchical \emph{oct-tree}
data structure; in contrast to ``full'' cells, ``empty'' cells only
retain the values at their surface. When the main search tries to
use one of them, its interior values are recomputed on demand. Still,
this work assumes that each node in the entire heuristic is calculated
at least once using dynamic programming.


We see one cause of the dilemma in the implicit assumption that a
\emph{complete} computation is necessary. The bound $\delta$
above refers to the worst-case, and can generally include
many more nodes than actually required in the main
search. However, since we are only dealing with the heuristic, we can 
actually afford to miss some values occasionally; while this might slow down
the main search, it cannot compromise the optimality of the final
solution. Therefore, we propose to generate the 
heuristics with a much smaller bound $\delta$. Whenever the attempt to
retrieve a value of the $m$-dimensional subheuristic fails during the
main search, we simply revert to replacing it by the sum of the $m
\choose 2$  optimal pairwise goal distances it covers.

%% We 
%% also maintain a statistics of the frequency of these misses. When the
%% fault rate rises above some threshold, we know that the
%% $\delta$-bound was chosen too low, interrupt the current search and restart it with
%% an increased threshold $\delta$. Similarly as in the case of the search
%% threshold of \iddp, since the use of an improved heuristic can lead to
%% exponential savings in time and space, this restarting scheme can be
%% performed in a way so that its contribution to the overall execution time
%% is bounded. 

We believe that the \iddp\ algorithm lends itself well
to make productive use of higher-dimensional heuristics. Firstly and
most importantly, the strategy of searching to adaptively increasing
thresholds can be transferred to the $\delta$-bound as well; this will
be addressed in more detail in the next section. 

Secondly, as far as a
practical implementation is concerned, it is 
important to take into account not only how a higher-dimensional
heuristic affects the number of node expansions,
but also their time complexity. This time is dominated by the
number of accesses to sub-alignments. With
$k$ sequences, in the worst case an edge has $2^k-1$
successors, leading to a total of $$(2^k -1){{k} \choose {m}}$$
evaluations for $h_{all,m}$. One possible improvement is to enumerate all
edges emerging from a given vertex in lexicographic order, and to store
partial sums of heuristics of prefix subsets of sequences for later
re-use. In this way, if we allow for a cache of linear size, the
number of accesses is reduced to   
$$\sum_{i={m}}^{i=k} 2^i {{i-1}
  \choose {m-1}};$$ 
correspondingly, for a quadratic cache we only need 
$$\sum_{i=m}^{i=k} 2^i {{i-2} \choose {m-2}}$$
evaluations. For instance, in aligning 12 sequences using $h_{all,3}$,
a linear cache reduces the evaluations to about 37 percent within one
expansion. 
%For a notable speedup in computation time, the reduction
%in expansions achieved by a heuristic has to at least offset the
%increase in computation time for one expansion. 

As mentioned above, in contrast to \ast, \iddp\ gives us the freedom to
choose any particular expansion order of the edges within a given
level. Therefore, when we sort edges lexicographically according 
to the target nodes, much of the cached prefix information can 
be shared additionally across consecutively expanded
edges. The higher the dimension of the subalignments, the larger are
the savings. In our experiments, we experienced speedups of up to
eighty percent in the heuristic evaluation. 



\subsection{Trade-Off between Computation of Heuristic and Main
  Search\label{sec:tradeoff}}
%In contrast to pairwise sub-alignments, time and space resources
%devoted to compute and store higher-dimensional heuristics are in
%general no longer negligible compared to the main search. 

As we have seen, we can control the size of the precomputed
sub-alignments by choosing the bound 
$\delta$ up to which $f$-values of edges are generated
beyond the respective optimal solution cost. There is
obviously a trade-off between the auxiliary and main searches. It is
instructive to consider 
the \emph{heuristic miss ratio} $r$, i.e., the fraction of calculations of
the heuristic $h$ during the main search when a requested entry in a
partial MSA has not been precomputed. The optimum
for the main search is achieved if the
heuristic has been computed for every requested edge ($r=0$). 
Going beyond that point will generate an unnecessarily large heuristic
containing many entries that will never be actually used. On the other
hand, we are free to allocate less effort to the heuristic, resulting in $r>0$
and consequently decreasing performance of the main search. Generally,
the dependence has an S-shaped form, as exemplified in 
Fig.~\ref{fig:time_tradeoff} for the
case of problem {\tt 1tvxA} of \bali\ (cf.~next section). Here, the 
execution time of one iteration of the main search at a fixed
threshold of 45 above the lower bound is shown, which includes the
optimal solution.

\begin{figure}
\begin{center}
%\epsfxsize8cm
%\epsfysize7cm
\GRAPHIC{time_tradeoff}
\end{center}
\caption{Trade-off between heuristic and main search:
  Execution times for problem {\tt 1tvxA} as a function of heuristic miss ratio.\label{fig:time_tradeoff}}
\end{figure}

Fig.~\ref{fig:time_tradeoff} illustrates the overall time trade-off
between auxiliary and main search, if we fix $\delta$ at different
levels. The minimum total
execution time, which is the sum of auxiliary and main search, is attained at
about $r=0.15$ (5.86 seconds). The plot for the corresponding memory usage
trade-off has a very similar shape. 
\begin{figure}[tb]
\begin{center}
%\epsfxsize8cm
%\epsfysize7cm
\GRAPHIC{time_tradeoff2}
\end{center}
\caption{Time of the last iteration in the main search for problem
  {\tt 1tvxA} as a function of heuristic miss ratio.\label{fig:time_tradeoff2}}
\end{figure}

Unfortunately, in general we do not know in advance the right amount
of auxiliary search. As mentioned above, choosing $\delta$ according
to the Carrillo-Lipman bound will ensure that every requested
sub-alignment cost will have been precomputed; however, in general we will
considerably overestimate the necessary size of the
heuristic. 

As a remedy, our algorithm \iddp\ gives us the opportunity to recompute
the heuristic in each threshold iteration in the main search. In this
way, we can adaptively strike a balance between the two. 

When the currently experienced miss rate $r$ rises above some
threshold, we can suspend the 
current search, recompute the pairwise alignments with an increased threshold
$\delta$, and resume the main search with the improved
heuristics. 

Like for the main search, we can accurately predict the
auxiliary computation time and space at threshold $\delta$ using
exponential fitting. Due to the lower dimensionality,
it will generally increase less steeply; however, the constant factor
might be higher for the heuristic,  due to the combinatorial number of
$ k \choose m$ alignment problems to be solved.

A doubling scheme 
as explained above can bound the overhead to within a constant
factor of the effort in the last iteration. In this way, when also
limiting the heuristic computation time by a fixed fraction of the main
search, we can ensure as an expected upper bound that the overall execution
time stays within a constant factor of the search time that would be
required using only the pairwise heuristic. 

If we knew the exact relation between $\delta$, $r$, and
the speedup of the main search, an ideal strategy would double the
heuristic whenever the expected computation time is smaller than the
time saved in the 
main search. However, as illustrated in Fig.~\ref{fig:time_tradeoff2},
this dependence is more complex than simple exponential growth, it
varies with the search depth and specifics of the problem. Either we
would need a more elaborate model of the 
search space, or the algorithm would have to conduct exploratory
searches in order to estimate the relation. 
We leave this issue to future work, and restrict ourselves here to a
simplified, conservative heuristic: 
We hypothesize that the main search can be made twice as fast by a
heuristic doubling if the miss rate $r$ rises above 25 percent; in
our experiments, we found that this assumption is almost always true.
In this event, since the effective branching factor of the main search
is reduced by the improved heuristic, we also ignore the history of
main search times in the exponential extrapolation procedure for
subsequent iterations.

\section{Experimental Results\label{sec:exp}}

In the following, we compare \iddp\ to one of the currently most successful approaches,
\emph{Partial Expansion} $A^*$. We empirically explore the benefit of higher-dimensional 
heuristics; finally, we show its feasibility by means of the benchmark database \bali.


\subsection{Comparison to Partial Expansion~\boldmath{$A^*$}}

For the first series of evaluations, 
%we compared \iddp\ to
%\emph{Partial Expansion} $A^*$ 
we ran \iddp\ on the same set of sequences as chosen by
\citeA{PE} (elongation factors EF-TU and EF-1$\alpha$ from various
species, with a high degree of similarity). As in this work,
substitution costs were chosen according to the PAM-250 matrix.
 The  applied heuristic was the sum of optimal pairwise goal distances. The
expansion numbers do not completely match with their results, however,
since we applied the biologically more realistic affine gap costs:
gaps of length $x$ were charged $8 + 8\cdot x$, except at the 
beginning and end of a sequence, where the penalty was $8\cdot
x$.   

All of the following experiments were run under RedHat Linux 7.3 on an
Intel Xeon\texttrademark\  CPU with 3.06 GHz, and main memory of 2
Gigabytes; we used the gcc 2.96 compiler. 

The total space consumption of a search algorithm is determined by the 
peak number of \emph{Open} and \emph{Closed} edges over the entire
running time. Table~\ref{tab:comp} and Fig.~\ref{fig:mem} give these
values for the series of 
successively  larger sets of input sequences 
%(with the sequences numbered as defined in \cite{PE}) $1-4$, $1-5$,
\cite<with the sequences numbered as defined in>{PE} $1-4$, $1-5$,
$\dots$, $1-12$.

\begin{table}[tb!]
\begin{center}
{\small
\begin{tabular}{rrrrr}
& {\bf Num} & {\bf Time} & {\bf Max} & {\bf
  Max} \\
 & {\bf Exp} & {\bf[sec]}& {\bf Open} & {\bf Open +}\\
& & & & {\bf Closed}\\
\hline \\
&\multicolumn{2}{r}{\boldmath $A^*$}\\
\input{new_experiment/table_ast.tbl}
\hline \\
&\multicolumn{2}{r}{\bf \pea}\\
\input{new_experiment/table_pea.tbl}
\hline \\
&\multicolumn{2}{r}{\bf \iddp}\\
\input {new_experiment/table_norm.tbl}
\hline
\end{tabular}
}
\end{center}
\caption{Algorithm comparison for varying number of input sequences (elongation factors EF-TU and EF-1$\alpha$).\label{tab:comp}}
\end{table}
\begin{figure}
\begin{center}
%\epsfxsize8cm
%\epsfysize7cm
\GRAPHIC{new_experiment/oc}
\end{center}
\caption{Memory requirements for \ast, \iddp, and \pea (elongation factors EF-TU and EF-1$\alpha$).\label{fig:mem}}
\end{figure}

With our implementation, the basic $A^*$ algorithm could be carried
out only up to 9 sequences, before exhausting our computer's  main
memory.

Confirming the results of \citeA{PE}, Partial Expansion requires
only about one percent of this space. Interestingly, during the
iteration with the peak in total numbers of nodes held in memory, no
nodes are actually closed except in problem 6. This might be explained
with the high degree of similarity between sequences in this example.
Recall that \pea\ only closes a node if all of its successors have an
$f$-value of no more than the optimal solution cost; if the span
to the lower bound is small, each node can have at least
one ``bad'' successor that exceeds this difference.

\iddp\
reduces the memory requirements further by a factor of about 6. The
diagram also 
shows the maximum size of the \emph{Open} list alone. For few sequences, the
difference between the two is dominated by the linear length to store
the solution path. As the problem size increases, however, the
proportion of the \emph{Closed} list of the total memory drops to about only 12
percent for 12 sequences. The total number of expansions (including
all search stages) is slightly higher than in \pea; however, due to
optimizations made possible by the control of the expansion order, the
execution time at 12 sequences is reduced by about a third.


Since \pea\ does not prune edges, its maximum space usage is always the total
number of edges with $f$-value smaller than 
$g^*(t)$ (call these edges the \emph{relevant edges}, since they have to be
inspected by each admissible algorithm). In \iddp, on the other
hand, the \emph{Open} list can only comprise 
$k$ adjacent levels out of those edges (not counting the possible threshold
overshoot, which would contribute a factor of at most 2). Thus, the
improvement of \iddp\  over \pea\ will tend to increase with 
the overall number of levels (which is the sum of all string lengths),
divided by the number of sequences; in other words, with the average
sequence length. 

Moreover, the ratio depends on how well the heuristic suits the
particular problem. 
%% Typically, the number of explored nodes has
%% an 'onion-shaped' distribution with the bulk being located closer to
%% the start node, and  tapering out towards higher levels, since the
%% accuracy of the heuristic increases inversely with target distance.
Fig.~\ref{fig:level} shows the distribution of all
edges with $f$ value smaller or equal to $g^*(t)$, for the case of 9
of the example sequences. This problem is quite extreme as the bulk of
these edges is concentrated in a 
small level band between 1050 and 1150. As an example with a more even
distribution, Fig.~\ref{fig:level_1cpt} depicts the situation for
problem {\tt 1cpt} from  Reference~1 in the benchmark set \bali\
\cite{bali} with heuristic $h_{all,3}$. In this case, the proportion of
the overall 19492675 relevant edges that are maximal among all 4
adjacent levels amounts to only 0.2 percent. The maximum \emph{Open} size in
\iddp\ is 7196, while the total number of edges generated by \pea\ is
327259, an improvement by about a factor of 45.

{
\begin{figure}[tb]
\begin{center}
%\epsfxsize8cm
%\epsfysize7cm
\GRAPHIC{lev}
\end{center}
\caption{Distribution of relevant edges over levels (elongation
  factors EF-TU and EF-1$\alpha$); compare to the schematic projection
  in Fig.~\ref{fig:layers}.\label{fig:level}}
\end{figure}
\samepage
\begin{figure}[tb]
\begin{center}
%\epsfxsize8cm
%\epsfysize7cm
\GRAPHIC{lev_1cpt}
\end{center}
\caption{Distribution of relevant edges over levels, problem {\tt 1cpt} from \bali.\label{fig:level_1cpt}}
\end{figure}
}

\subsection{Multidimensional Heuristics\label{sec:exp_heu_cmp}}


On the same set of sequences, we compared different
improved heuristics in order to get an impression for their respective
potential. Specifically, we ran \iddp\ with heuristics $h_{pair}$,
$h_{all,3}$, $h_{all,4}$, and 
$h_{one,k/2}$ at various thresholds $\delta$. Fig.~\ref{fig:heu_time}
shows the total execution time for computing the heuristics, and
performing the main search. In each case, we manually selected a value
for $\delta$ which minimized this time. It can be seen that the times for
$h_{one,k/2}$ lie only a little bit below  $h_{pair}$; For
few sequences (less than six), the computation of the 
heuristics $h_{all,3}$ and $h_{all,4}$ dominates their overall time. 
With increasing dimensions, however, this investment
starts to yield growing returns, with $h_{all,3}$ being the fastest
algorithm, requiring only 5 percent of the time of $h_{pair}$ at 12 sequences. 
\begin{figure}[tb]
\begin{center}
%\epsfxsize8cm
%\epsfysize7cm
\GRAPHIC{new_experiment/heu_time}
\end{center}
\caption{Comparison of execution times (including calculation of
  heuristics), elongation factors EF-TU and EF-1$\alpha$.\label{fig:heu_time}}
\end{figure}

As far as memory is concerned, Fig.~\ref{fig:heu_mem} reveals that the maximum size of the \emph{Open} and \emph{Closed}
list, for the chosen $\delta$ values, is very similar for  $h_{pair}$
and $h_{one,k/2}$ on the one hand, and  $h_{all,3}$ and  $h_{all,4}$
on the other hand.

At 12 sequences, $h_{one,6}$ saves only about 60 percent
of edges, while  $h_{all,3}$ only needs 2.6 percent and $h_{all,4}$
only 0.4 percent of the space required by the pairwise heuristic.
%(not counting the space for the heuristics itself). 
Using \iddp, we never ran out of main memory; even larger test sets could be
aligned, the range of the shown diagrams was limited by our patience
to wait for the results for more than two days.

\begin{figure}
\begin{center}
%\epsfxsize8cm
%\epsfysize7cm
\GRAPHIC{new_experiment/heu_mem}
\end{center}
\caption{Combined maximum size of \emph{Open} and \emph{Closed}, for different heuristics (elongation factors EF-TU and EF-1$\alpha$).\label{fig:heu_mem}}
\end{figure}

Based on the experienced burden of computing the heuristic,
\citeA{kobayashi-improvement} concluded that $h_{one,m}$ should be preferred to
$h_{all,m}$. We do not quite agree with this judgment. We see that 
the heuristic $h_{all,m}$ is able to
reduce the search space of the main search considerably stronger than
$h_{one,m}$, so that it can be more beneficial with an appropriate
amount of heuristic computation.

\subsection{The Benchmark Database \bali\label{sec:bali}}

\bali\ \cite{bali} is a widely used database of manually-refined
multiple sequence alignments specifically designed for the evaluation
and comparison of multiple sequence alignment programs. 
%Core blocks are identified excluding non-superimposable regions. 
The alignments are classified into
8 reference sets. Reference~1 contains alignments of up to six
about equidistant  sequences. All the
sequences are of similar length; they are grouped into 9 classes,
indexed by sequence length and the percentage of identical amino
acids in the same columns. Note that many of these
problems are indeed much harder than the elongation factor examples from
the previous section; despite consisting of fewer
sequences, their dissimilarities are much more pronounced.

We applied our algorithm to Reference~1, with substitution costs
according to the PET91 matrix \cite{pet91} and affine gap costs of $9
\cdot x + 8$, except for leading and trailing gaps, where no gap
opening penalty was charged. For all 
instances, we precomputed the pairwise sub-alignments up to a fixed
bound of 300 above the optimal solution; the optimal solution was
found within this bound in all cases, and the effort is
generally marginal compared to the overall computation.
For all problems involving more than three sequences, the heuristic
$h_{all,3}$ was applied.
 
Out of the 82 alignment problems in Reference~1, our algorithm could
solve all but 2 problems (namely, \emph{1pamA} and \emph{gal4}) on our
computer. Detailed results are listed in 
Tables~\ref{tab:bali_1_1} through \ref{tab:bali_3_3}. 

%% %Fig.~\ref{fig:balistat} shows the so-called \emph{SP}-score in the core
%% %blocks, that is, the percentage of correctly aligned pairs within the
%% %biologically significant motifs. 

Thompson, Plewniak, and Poch (1999) compared a number of widely used heuristic alignment
tools using the so-called \emph{SP}-score; their software calculates
the percentage of correctly aligned pairs within the biologically
significant motifs. They found that all programs perform about equally
well for the sequences with medium and high amino acid identity;
differences only occurred for the case  of the more distant sequences
with less than 25 percent identity, the so-called ``twilight
zone''. Particularly challenging was the group of short
sequences. In this subgroup, the three highest scoring programs are 
PRRP, CLUSTALX, and SAGA, with respective median scores of 0.560,
0.687, and 0.529. The medium score for the alignments found in our
experiments amounts to 0.558; hence, it is about as good as PRRP, and
only beaten by CLUSTALX. While we focused in our experiments on
algorithmic feasibility rather than on solution quality, it would be
worthwhile to attempt to improve the alignments found by these program
using their more refined penalty functions. CLUSTALX, for example, uses
different PAM matrices depending on the evolutionary distance of
sequences; moreover, it assigns weights to sequences (based on a
phylogenetic tree), and gap penalties are made
position-specific. All of these  improvements can be easily integrated
into the basic sum-of-pairs 
cost function, so that we could attempt to compute an optimal
alignment with respect to these metrics. We leave this line of
research for future work.

%% \begin{table}[htb]
%% \begin{center}
%% \begin{tabular}{r|rrrr}
%% & {\bf Problem} & {\bf PRRP} & {\bf CLUSTALX} & {\bf
%%   SAGA} \\
%% \hline
%% 1aboA & 0.560 & 0.687 & 0.529 & \\
%% 1idy & 0.606 & 0.705 & \\
%% \end{tabular}
%% \end{center}
%% \caption{SP-score for hardest problems in Reference 1.\label{tab:V!}}
%% \end{table}

%% For sequences of medium and high
%% similarity, the median score 
%% lies between 93.8\% and 98.6\%. The more interesting case is that of
%% the more distant sequences with less than 25 percent identity, the
%% so-called twilight zone. The problems are particularly hard for short
%% sequences, where the median score reaches  42.5\%; for medium-size and
%% long sequences, the score amounts to 68.4\% and 63.3\%,
%% respectively. The overall median score is 95.2\%.

%% \begin{figure}
%% \begin{center}
%% %\epsfxsize8cm
%% \epsfysize7cm
%% \GRAPHIC{balistat}
%% \end{center}
%% \caption{Median \emph{SP} score for Reference~1 of \bali.\label{fig:balistat}}
%% \end{figure}

Fig.~\ref{fig:bali_open} shows the maximum number of edges that have
 to be stored in \emph{Open} during the search, in dependence of the
 search threshold in the final iteration. For better comparability, we
 only included those problems in the diagram that consist of 5
 sequences.  The logarithmic  scale emphasizes that the growth fits an
 exponential curve quite  well. Roughly speaking, an increase of the
 cost threshold by 50 leads  to a ten-fold increase in the space
 requirements. This relation is  similarly applicable to the number of
 expansions  (Fig.~\ref{fig:bali_iter}).
\begin{figure}
\begin{center}
%\epsfxsize8cm
%\epsfysize7cm
\GRAPHIC{bali_open}
\end{center}
\caption{Maximum size of \emph{Open} list, dependent on the final search
  threshold (\bali).\label{fig:bali_open}}
\end{figure}
\begin{figure}
\begin{center}
%\epsfxsize8cm
%\epsfysize7cm
\GRAPHIC{bali_iter}
\end{center}
\caption{Number of expansions in the final search iteration (\bali).\label{fig:bali_iter}}
\end{figure}

Fig.~\ref{fig:bali_open_closed}
depicts the proportion between the maximum \emph{Open} list size and
the combined maximum size of \emph{Open} and \emph{Closed}. It is
clearly visible that due to the pruning of edges outside of possible
solution paths, the  \emph{Closed} list contributes less and less to the
overall space requirements the more difficult the problems become.
\begin{figure}
\begin{center}
%\epsfxsize8cm
%\epsfysize7cm
\GRAPHIC{bali_open_closed}
\end{center}
\caption{Maximum number of \emph{Open} edges, divided by combined maximum of
  \emph{Open} and \emph{Closed} (\bali).\label{fig:bali_open_closed}}
\end{figure}

Finally, we estimate the reduction in the size of the \emph{Open} list
compared to all relevant edges by the ratio of the maximum \emph{Open} size in the
last iteration of \iddp\ to the total number of expansions in this
stage, which is equal to the number of edges with $f$-value less or
equal to the threshold. Considering possible overshoot of \iddp,
algorithm \pea\ would expand 
at least half of these nodes. The proportion ranges between 0.5 to 5
percent (cf.~Fig.~\ref{fig:bali_band}). Its considerable scatter
indicates the dependence on individual problem properties; however, a
slight average decrease can be noticed for the more difficult problems.
\begin{figure}
\begin{center}
%\epsfxsize8cm
%\epsfysize7cm
\GRAPHIC{bali_band}
\end{center}
\caption{Percentage of reduction in \emph{Open} size (\bali).\label{fig:bali_band}}
\end{figure}

\clearpage
\goodbreak

\section{Conclusion and Discussion}

We have presented a new search algorithm for optimal multiple sequence
alignment that combines the effective use of a heuristic bound as in
best-first search with the ability of the dynamic programming approach
to reduce the maximum size of the \emph{Open} and \emph{Closed} lists
by up to one order 
of magnitude of the sequence length. The algorithm performs a
series of searches with successively increasing bounds that explore
the search space in \DP\ 
order; the thresholds are chosen adaptively so that the expected overhead
in recomputations is bounded by a constant factor. 

We have demonstrated
that the algorithm can outperform one of the  currently most successful
algorithms for optimal multiple sequence alignments, \emph{Partial
  Expansion} \ast, both in terms of 
computation time and memory consumption. Moreover, the
iterative-deepening strategy alleviates the use of partially computed
higher-dimensional heuristics. To the best of our
knowledge, the algorithm is the first one that is  
able to solve standard benchmark alignment problems in \bali\  with a
biologically realistic cost function including affine gap costs
without end gap penalties. The quality of the alignment is in the 
range of the best heuristic programs; while we have concentrated on
algorithmic feasibility, we deem it worthwhile to incorporate their
refined cost metrics for better results; we will study this question
in future work.

%In the final version of this paper, we will
%be able to provide the source code of the program.

Recently, we learned about related approaches developed simultaneously 
and independently by \citeA{zhou03b,zhou04}. \emph{Sweep}\ast\ explores a search graph according to 
layers in a partial order, but still uses the $f$-value for selecting
nodes within one layer. \emph{Breadth-First Heuristic Search}
implicitly defines the layers in a graph with uniform costs according
to the breadth-first traversal. Both algorithms incorporate upper
bounds on the optimal solution cost for pruning; however, the idea of
adaptive threshold determination to limit re-expansion overhead to a
constant factor is not described. Moreover, they do not consider the
flexible use of additional memory to minimize the divide-and-conquer
solution reconstruction phase.  

Although we described our algorithm entirely within the framework of
the MSA problem, it is straightforward to transfer it to any domain
in which the state space graph is directed and acyclic. Natural
candidates include applications where such an ordering is imposed
by time or space coordinates, e.g., finding the most likely
path in a Markov model.

Two of the \bali\ benchmark problems could still not be solved by our
algorithm within the computer's main memory limit. Future work will
include the integration of techniques exploiting secondary 
memory. We expect that the level-wise exploration scheme of our
algorithm lends itself naturally to \emph{external search algorithms},
another currently very active research topic in Artificial
Intelligence and theoretical computer science.

\acks{The author would like to thank the reviewers of this article whose comments 
have helped in significantly improving it.}

\newpage
{\samepage
\onecolumn
\section*{Appendix A}

\nopagebreak
\begin{table}[h]
\caption{Results for \bali\, Reference~1, group of short sequences with low
  amino acid identity. The columns denote: S --- number of aligned
  sequences; $\delta$ --- upper bound for precomputing optimal solutions
  for partial problems in last iteration of main search;  $g^*(t)$ --- optimal
  solution cost; $h(s)$ --- lower bound for
  solution cost, using heuristics; \#Exp --- total number of
  expansions in all iterations of the main search;
  \#Op --- peak number of edges in \emph{Open} list over the course 
  of the search; \#Op+Cl --- peak combined number of edges in either
  \emph{Open} or \emph{Closed} list during search; \#Heu --- peak
  number of sub---alignment edge costs stored as heuristic; Time: total
  running time including auxiliary and main search, in seconds; Mem ---
  peak total memory usage for face alignments, heuristic, and main
  search, in KB.\label{tab:bali_1_1}} 
\begin{center}
\small
\begin{tabular}{lrrrrrrrrrrr}
 &  {\bf S} &  {\boldmath $\delta$} &  {\boldmath
  $g^*(t)$} & {\boldmath $h(s)$}  &  {\bf \#Exp} & {\bf\#Op} &  {\bf \#Op+Cl} & {\bf \#Heu} &
{\bf Time} & {\bf Mem}\\ 
\input{adap_heu/results_1_1.tex}
\end{tabular}
\end{center}
\end{table}
}

\begin{table}[h]
\caption{Short sequences, medium similarity.\label{tab:bali_1_2}}
\begin{center}
\small
\begin{tabular}{lrrrrrrrrrrr}
 &  {\bf S} &  {\boldmath $\delta$} &  {\boldmath
  $g^*(t)$} & {\boldmath $h(s)$}  &  {\bf \#Exp} & {\bf\#Op} &  {\bf \#Op+Cl} & {\bf \#Heu} &
{\bf Time} & {\bf Mem}\\ 
\input{adap_heu/results_1_2.tex}
\end{tabular}
\end{center}
\end{table}

\begin{table}[h]
\caption{Short sequences, high similarity.\label{tab:bali_1_3}}
\begin{center}
\small
\begin{tabular}{lrrrrrrrrrrr}
 &  {\bf S} &  {\boldmath $\delta$} &  {\boldmath
  $g^*(t)$} & {\boldmath $h(s)$}  &  {\bf \#Exp} & {\bf\#Op} &  {\bf \#Op+Cl} & {\bf \#Heu} &
{\bf Time} & {\bf Mem}\\ 
\input{adap_heu/results_1_3.tex}
\end{tabular}
\end{center}
\end{table}

\begin{table}[h]
\caption{Medium-length sequences, low similarity.\label{tab:bali_2_1}}
\begin{center}
\small
\begin{tabular}{lrrrrrrrrrrr}
 &  {\bf S} &  {\boldmath $\delta$} &  {\boldmath
  $g^*(t)$} & {\boldmath $h(s)$}  &  {\bf \#Exp} & {\bf\#Op} &  {\bf \#Op+Cl} & {\bf \#Heu} &
{\bf Time} & {\bf Mem}\\ 
\input{adap_heu/results_2_1.tex}
\end{tabular}
\end{center}
\end{table}

\begin{table}[h]
\caption{Medium-length sequences, medium similarity.\label{tab:bali_2_2}}
\begin{center}
\small
\begin{tabular}{lrrrrrrrrrrr}
 &  {\bf S} &  {\boldmath $\delta$} &  {\boldmath
  $g^*(t)$} & {\boldmath $h(s)$}  &  {\bf \#Exp} & {\bf\#Op} &  {\bf \#Op+Cl} & {\bf \#Heu} &
{\bf Time} & {\bf Mem}\\ 
\input{adap_heu/results_2_2.tex}
\end{tabular}
\end{center}
\end{table}

\begin{table}[h]
\caption{Medium-length sequences, high similarity.\label{tab:bali_2_3}}
\begin{center}
\small
\begin{tabular}{lrrrrrrrrrrr}
 &  {\bf S} &  {\boldmath $\delta$} &  {\boldmath
  $g^*(t)$} & {\boldmath $h(s)$}  &  {\bf \#Exp} & {\bf\#Op} &  {\bf \#Op+Cl} & {\bf \#Heu} &
{\bf Time} & {\bf Mem}\\ 
\input{adap_heu/results_2_3.tex}
\end{tabular}
\end{center}
\end{table}

\begin{table}[h]
\caption{Long sequences, low similarity.\label{tab:bali_3_1}}
\begin{center}
\small
\begin{tabular}{lrrrrrrrrrrr}
 &  {\bf S} &  {\boldmath $\delta$} &  {\boldmath
  $g^*(t)$} & {\boldmath $h(s)$}  &  {\bf \#Exp} & {\bf\#Op} &  {\bf \#Op+Cl} & {\bf \#Heu} &
{\bf Time} & {\bf Mem}\\ 
\input{adap_heu/results_3_1.tex}
\end{tabular}
\end{center}
\end{table}

\begin{table}[h]
\caption{Long sequences, medium similarity.\label{tab:bali_3_2}}
\begin{center}
\small
\begin{tabular}{lrrrrrrrrrrr}
 &  {\bf S} &  {\boldmath $\delta$} &  {\boldmath
  $g^*(t)$} & {\boldmath $h(s)$}  &  {\bf \#Exp} & {\bf\#Op} &  {\bf \#Op+Cl} & {\bf \#Heu} &
{\bf Time} & {\bf Mem}\\ 
\input{adap_heu/results_3_2.tex}
\end{tabular}
\end{center}
\end{table}

\begin{table}[h]
\
\caption{Long sequences, high similarity.\label{tab:bali_3_3}}
\begin{center}
\small
\begin{tabular}{lrrrrrrrrrrr}
 &  {\bf S} &  {\boldmath $\delta$} &  {\boldmath
  $g^*(t)$} & {\boldmath $h(s)$}  &  {\bf \#Exp} & {\bf\#Op} &  {\bf \#Op+Cl} & {\bf \#Heu} &
{\bf Time} & {\bf Mem}\\ 
\input{adap_heu/results_3_3.tex}
\end{tabular}
\end{center}
\end{table}

\clearpage
%\newpage
%\bibliographystyle{msa}
\bibliographystyle{theapa}
% argument is your BibTeX string definitions and bibliography database(s)
\bibliography{align}
%

% that's all folks
\end{document}


