In this section, the basic dynamic programming method for sequence
alignment is described \cite{SalBlu90}. 
This method forms the core of the pairwise and
multiple sequence and structure comparisons as well as of the sequence
database searching.


\subsection{Pairwise comparison}
\label{SECTION:dp}

The residue by residue scores $W_{ij}$ can be used directly in the
sequence alignment algorithm of Needleman \& Wunsch \cite{NeeWun70}
to obtain the comparison of two protein sequences or structures. The only
difference between the two types of comparison is in the type
of the comparison matrix. In the case of sequence, the amino acid 
substitution matrix is used. In the case of 3D structure, the Euclidean 
distance (or some function of it) between two equivalent atoms in 
the current optimal superposition is used \cite{SalBlu90}.

The problem of the optimal alignment of two sequences as addressed by
the algorithm of Needleman \& Wunsch is as follows. We are given two
sequences of elements and an $M$ times $N$ score matrix $\cal W$ where
$M$ and $N$ are the numbers of elements in the first and second
sequence. The scoring matrix is composed of scores $W_{ij}$ describing
differences between elements $i$ and $j$ from the first and second
sequence respectively. The goal is to obtain an optimal set of
equivalences that match elements of the first sequence to the elements
of the second sequence. The equivalence assignments are subject to the
following ``progression rule'': for elements $i$ and $k$ from the
first sequence and elements $j$ and $l$ from the second sequence, if
element $i$ is equivalenced to element $j$, if element $k$ is
equivalenced to element $l$ and if $k$ is greater than $i$, $l$
must also be greater than $j$. The optimal set of equivalences is the
one with the smallest alignment score. The alignment score is a sum of
scores corresponding to matched elements, also increased for
occurrences of non-equivalenced elements ({\it ie} gaps). For a detailed
discussion of this and related problems see \cite{SanKru83}.

We summarize the dynamic programming formulae used by \MODELLER\ to
obtain the optimal alignment since they differ slightly from those
already published \cite{SelAlg74,Got82}. The recursive dynamic
programming formulae that give a matrix $\cal D$ are:

\begin{equation}
  \begin{array}{ccl}
  D_{i,j} & = & \min \left\{ \begin{array}{l}
                          P_{i,j}                 \\
                          D_{i-1,j-1} + W_{i,j}   \\
                          Q_{i,j}
                 \end{array} \right.      \vspace{0.2in} \\
  P_{i,j} & = & \min \left\{ \begin{array}{l}
                         D_{i-1,j} + g(1)         \\
                         P_{i-1,j} + v
                 \end{array} \right.      \vspace{0.2in} \\
  Q_{i,j} & = & \min \left\{ \begin{array}{l}
                         D_{i,j-1} + g(1)         \\
                         Q_{i,j-1} + v
                 \end{array} \right.      \vspace{0.7in}
  \end{array}
  \label{eq:29}
\end{equation}
where $g(l)$ is a linear gap penalty function:
\begin{equation}
   g(l) = u + v \cdot l \, .
   \label{eq:30}
\end{equation}
Note that only a vector is needed for the storage of $P$ and $Q$.  The
uppermost formula in Eq.~\ref{eq:29} is calculated for $i=M$ and
$j=N$. Variable $l$ is a gap length and parameters $u$ and $v$ are
gap-penalty constants.

The arrays $\cal D$, $\cal P$ and $\cal Q$ are initialized as follows:
\begin{equation}
 \begin{array}{ccl}
  D_{i,0} & = & \left\{ \begin{array}{lc}
			 0,    &       i \leq e   \\
			 g(i-e), & e  <  i \leq N
                 \end{array}  \right.         \vspace{0.2in}  \\
  D_{0,j} & = & \left\{ \begin{array}{lc}
			 0,    &       j \leq e   \\

			 g(j-e), & e  <  j \leq N
                 \end{array}  \right.        \vspace{0.2in}   \\
  P_{i,0} & = & Q_{i,0} = \infty , \qquad i = 1,2,\ldots,M  \vspace{0.1in} \\
  P_{0,j} & = & Q_{0,j} = \infty , \qquad j = 1,2,\ldots,N
 \end{array}
 \label{eq:31}
\end{equation}
where parameter $e$ is the maximal number of elements at sequence
termini which are not penalized with a gap-penalty if not
equivalenced. A segment at the terminus of length $e$ is termed an
``overhang''. Note a difference from \cite{Got82} in the
initialization of the $\cal P$ and $\cal Q$ arrays. Also note that
only vectors $Q_i$ and $P_j$ need to be stored in computer, not the
whole arrays.

The minimal score $d_{M,N}$ is obtained from
\begin{equation}
  d_{M,N} = \min(D_{i,N}, D_{M,j})
  \label{eq:32}
\end{equation}
where $i = M,M-1,\ldots,M-e \quad \hbox{and} \quad j =
N,N-1,\ldots,N-e$ to allow for the overhangs. The equivalence
assignments are obtained by backtracking in matrix $\cal D$.
Backtracking starts from the element $D_{i,j} = d_{M,N}$.


\subsection{Variable gap penalty}
\label{SECTION:variablegap}

This work is still in progress and is not described here.

% The recursive equation for calculating the $D$ array when using the variable gap
% penalty function is:
% \begin{equation}
%   \begin{array}{ccl}
%   D_{i,j} & = & \min \left\{ \begin{array}{l}
%                           \min_{1 \le l \le L, \; 1 \le l \le L} 
%                           D_{i-l,j} + g_i(l) \\
%                           D_{i-1,j-1} + W_{i,j}     \\
%                       \end{array} \right.  \\[8mm]
%   G_{i,j} & = & \min_{1 \le l \le L, \; 1 \le k \le L}
%                          D_{l,k} + g_j(1) \\
%                  \end{array} \right.      \\[7mm]
%   g_i(l) & = & u f_2 + v l \qquad \mbox{gap opening or inserting in sequence `j'}  \\[2mm]
%   g_j(l) & = & u f_1 + v l \qquad \mbox{gap opening or inserting in structure `i'}  \\
%   \end{array}
%   \label{eq:29b}
% \end{equation}
% See Section~\ref{SECTION:dp} for the description of recursion, $D$, $Q$, $W$, 
% $u$, $v$, array initialization, and backtracking.
% See Section~\ref{COMMAND:align2d} for the definition of functions $f_1$
% and $f_2$. Functions $f_1$ and and $f_2$ modulate the structure dependent, 
% variable gap opening penalty. The gap extension penalty, $v$, is a 
% constant, as is the other factor in the gap opening penalty, $u$.
% The only difference in the recursion from the linear gap penalty case
% (Eq.~\ref{eq:29}) is that a slightly slower procedure for finding
% the optimal gap lengths must be used for gap openings in sequence
% `$j$' because of the penalty dependence on the distance between the
% two spanning \Calpha\ atoms in structure `$i$'. The CPU time is saved
% by limiting the minimization over $l$ to values of $l$ smaller than
% $L$; this is equivalent to limiting the maximal length of a gap to 
% $L = \K{OFF\_DIAGONAL}$ positions.  
% 
% The only \MODELLER\ command using the variable gap penalty is \C{ALIGN2D}.
% 
% It seems that the same algorithm can also be used for more complicated
% length dependencies of gap penalty. For example, the exponential length
% dependence could be used, following Gonnet and Benner. Another
% possible improvement is to determine the weights in $f_1$ and $f_2$ by
% analyzing reference structure--structure alignments. The algorithm
% could also be used to align structures with structures, in which case
% it would be symmetric.

\subsection{Local \versus\ global alignment}

The Kruskal and Sankoff version of the local alignment is implemented
\cite{SanKru83}; this is very similar to the \cite{SmiWat81} method.
All the routines for the local alignment are exactly the same as the
routines for the global alignment except that during the construction
of matrix $D$ the alignment is restarted each time the score becomes
higher than a cutoff. The second difference is that the backtracking
starts from the lowest element in the matrix, wherever it is.

\subsection{Similarity \versus\ distance scores}

Each scoring matrix contains a flag determining whether it is a
distance or similarity matrix. An appropriate optimization is used
automatically. This is achieved by using exactly the same code except
that one side of comparisons is multiplied by $-1$ when dealing with
similarities as opposed to distances.


\subsection{Multiple comparisons}
   \label{sec:multiple}

In the discussion of the previous section, we have assumed that the
sequences or structures would be compared in a pairwise manner.
However, such pairwise comparisons of several related proteins may not
be self consistent, {\it ie} the following transitivity rule can be
broken: If residue $a$ from protein $A$ is equivalent to residue $b$
in protein $B$ which in turn is equivalent to residue $c$ in protein
$C$ then the residue $a$ from protein $A$ must also be equivalent to
residue $c$ from protein $C$. This property is not always attained in
the set of usual pairwise comparisons relating a group of similar
proteins. For this reason we proceed by simultaneously aligning all
proteins. This is achieved by aligning the second sequence with the
first one, the third sequence with the alignment of the first two, 
\etc. A more general tree-like growth of the multiple alignment is
not yet implemented.

% Recently, several methods for efficient multiple sequence alignment
% have been described (Hogeweg \& Hesper, J.Mol.Evol.
% 20, 175--186, 1984;  Barton \& Sternberg,  J. Mol. Biol. 198, 327--337,
% 1987; Feng \& Doolittle, J.Mol.Evol. 25, 351--360, 1987). The 
% essence of these approaches is in
% the way of overcoming the computationally too exhaustive dynamic
% programming algorithm for finding the global minimum of the multiple
% alignment. We adopt here a simplification of the `tree' alignment by
% successively adding the proteins to the growing alignment, as they
% occur in the input file; thus, the tree has a very simple and
% pre-determined structure and the final alignment may depend on the
% order of proteins in the input file. The following description assumes
% a general tree.
% 
% The alignment procedure follows the tree structure, from the tips of
% the branches to the root. First, the two most similar proteins are
% aligned using the pairwise procedure described above to give the first
% sub-alignment. Then either the third protein is compared with the
% first sub-alignment or the second most similar pair of proteins is
% compared in the usual pairwise way.  Which of the two alternatives
% will be realized depends on the actual tree topology. Suppose the
% latter is the case and the tree then implies that the two already
% aligned pairs have to be joined next.  This is again accomplished with
% the standard pairwise alignment technique
% except that the scoring matrix $\cal W$ is defined by averaging the
% differences between the two groups of already aligned proteins. In
% general, when two sub-alignments have to be compared, the scoring
% matrix element $W_{ij}$ is defined by
% \begin{equation}
%   W_{ij} = {1\over{n_1 \cdot n_2}} \, \sum_{l=1}^{n_1} \sum_{k=1}^{n_2}
%                \left\{ \begin{array}{ll}
%                W_{i'j'}^{lk} & \mbox{if two residues are compared} \\
%                G             & \mbox{if a residue and a gap are compared} \\
%                0             & \mbox{if two gaps are compared}
%                \end{array}  \right. \, .
%   \label{eq:34}
% \end{equation}
% Indices $i$ and $j$ are indices of the two positions compared, for the
% first and second sub-alignment respectively. Subscripts $i'$ and $j'$
% are true indices for residues from proteins $l$ and $k$, occupying
% sub-alignment positions $i$ and $j$, respectively. $W^{lk}_{i'j'}$ is
% an element of a familiar scoring matrix $\cal W$
% for a pairwise comparison of proteins $l$ and $k$.
% Parameter $G$ is the corresponding score from the scoring matrix (gap
% residue type is 21). If the scoring matrix does not explicitly specify
% this score, it is 0.  Parameters $n_1$ and $n_2$ are the numbers of
% proteins in the first and second sub-alignment, respectively.

If the number of all proteins is $N$, $N-1$ alignments must be made to
obtain the final multiple comparison. It is noted that once an
equivalence or gap is introduced it is not changed in later stages.

% The tree-like addition of proteins (Hogeweg \& Hesper, 1984; Feng \&
% Doolittle, 1987) implemented here is contrasted to its special case
% --- a simple sequential addition (Barton \& Sternberg, 1987). The
% rationale behind this choice is that the most similar proteins, and
% therefore the ones that are compared most accurately, are aligned
% first.  In addition to a better alignment of similar proteins, this
% also results in improved performance in later stages, since superior
% sub-alignments are used. Our experience is that the tree-like addition
% performs slightly better than the simple sequential addition ---
% especially when one aligns a group of proteins with wide range of
% similarities.
