\section{Introduction}

Computational chemists working in the area of 
structured-based drug design consider 
both the chemical and the geometric properties of the interacting
molecules when developing new pharmaceutical 
drugs~\cite{balbes,itai,lengauer-gen}.  The underlying assumption 
is that drug activity, or pharmacophoric activity, is obtained through
the molecular recognition and binding of one molecule (ligand) to a
pocket of another, usually larger, molecule (receptor). 
This assumption is supported by a 
number of experimental observations of 
molecules with geometric and chemical
complementarity in their active, or binding, conformations~\cite{}.

When the three dimensional structure of the receptor is known,
docking methods~\cite{} exploit both the geometric and the chemical information
available. However, the geometric structure of relatively few molecules 
has been obtained via X-ray crystallography or NMR techniques.
As a result, computational chemists often try to 
develop pharmaceutical drugs for receptors whose structure is
unknown~\cite{barnum,martin-pha}. 
The starting point in this case 
is a collection of ligands that have been experimentally found
to interact with the considered receptor. 
By examining the chemical properties and the possible shapes of these
ligands, chemists seek to identify  
a set of features embedded in $\Re^3$ that is contained
in all, or most, of the active conformations of all ligands.
This is called the pharmacophore and it is considered responsible for 
the observed drug activity.
The features of the pharmacophore interact with features of the receptor,
while the rest of the ligand acts as a scaffold. Once a
pharmacophore has been isolated, it can be used to 
further improve the activity  of pharmaceutical drugs~\cite{glen}.

This paper considers the following problem: Given a set of ligands that
interact with the same receptor, find geometric invariants of these ligands,
i.e., a set of features embedded in $\Re^3$ that is present in one or more
valid conformations of the ligands.  
We refer to this problem as the 
pharmacophore identification problem since the pharmacophore is such an  
invariant. Solving this problem requires dealing 
efficiently with large amounts of spatial data and shape information.
Ligand molecules are very flexible and can assume 
many distinct potentially valid conformations. 
A valid conformation is a rigid spatial realization of the atoms of a
molecule whose energy is below a predefined threshold \cite{}.
It is not known, a priori, which of the many possible valid conformations of a
ligand is the one that contains the desired pharmacophore. 

Besides providing starting points for  designing  effective pharmaceutical 
drugs,  
there are several important applications where geometric
invariant identification is useful.
For example, the geometric invariants identified by our system
can help in the formulation of database queries that may retrieve
functionally equivalent, but structurally novel molecules
from molecular databases \cite{willet}.
Geometric invariants can also suggest alignments of molecules
for input to CoMFA (Comparative Molecular Field Analysis) and
other 3D QSAR (Quantitative Structure-Activity Relationship) methods
\cite{barnum,itai}. These are methods that predict the activity of
hypothetical compounds on the assayed activity of previously
synthesized ones.

This paper describes our efforts to prototype an 
an integrated software system, called RAPID
(RAndomized Pharmacophore Identification for Drug design),
for addressing the pharmacophore identification problem.  
We present briefly the overall 
structure of RAPID in Section \ref{sec:approach}, 
and outline related work in Section~\ref{sec:related}.
The two main modules of RAPID, conformational
search and identification of invariants, are described in
Sections~\ref{sec:conf} and~\ref{sec:inv}.
Our system attempts to deal with
the numerous complications that arise in real examples. We 
report in Section~\ref{sec:exp} some preliminary experimental
results. Finally, in 
Section~\ref{sec:disc}, we conclude with a discussion of our short-term
goals in further developing our system, as well some open issues that
merit consideration in future work.
 
\section{Overview of RAPID}
\label{sec:approach}

\begin{figure}
\centerline{
\setlength{\unitlength}{1mm}
\begin{picture}(150,40)
\put(10,5){\psfig{figure=2mols_new.eps,height=3.5cm}}
\put(100,5){\psfig{figure=points.xfig.eps,height=2.5cm}}
\put(35,0){(a)}
\put(65,0){(b)}
\put(117,0){(c)}
\end{picture}
}
\caption{1TMN: (a) space-filling, (b) stick, and (c) set-of-points models}
\label{fig:model}
\end{figure}

RAPID tries to identify geometric invariants among a collection of 
small ligand molecules like the one shown in Figure\ref{fig:model}.
This ligand is called 1TMN and it is one of the several 
inhibitors of the protease 
thermolysin. Figure\ref{fig:model}(a) shows  
the space filling model of 1TMN. In this model a 
van der Walls sphere~\cite{} is considered around each atom center.
Figure\ref{fig:model}(b) shows the corresponding stick model in which only  
chemical bonds are drawn. The degrees of freedom of ligands   
include bond lengths, bond angles
(angles between two consecutive bonds), and dihedral or torsional angles
(angles formed by the first and
third of three consecutive bonds, viewed along the
axis of the second bond). In practice, only the torsional degrees of
freedom are considered since these are the ones that exhibit large
variations in their values \cite{}. Figure\ref{fig:model}(c) shows
the conformation of 1TMN of Figure\ref{fig:model}(a) as a set of points in
$\Re^3$. These points may represent atom centers
or groups of atoms aggregated to one point endowed with
a feature common to all these atoms (e.g., a rigid benzene ring) \cite{}.
%The representation of a
%conformation by points is consistent with the projected
%uses of our results (candidate pharmacophores, database queries)
%since bond information is not relevant in these applications.
We assume that once a conformation is given, one can automatically
transform it to a unique collection of points.
For the purposes of this paper, a conformation $C$
may represent a valid geometric embedding of a molecule where
bonds are retained, or a set of points without bond information. 

In RAPID the identification of geometric invariants in a collection of
flexible ligands denoted by $M = \{M_1, M_2, \ldots, M_N\}$ is treated 
as a two-stage process \footnote{A third module of RAPID, currently
under development, involves the computation of molecular surfaces (see
Finn et al~\cite{finn-et-al}).} addressing  the two following problems:

\begin{problem} {\bf (Conformational Search)}
Given a collection of ligands $M = \{M_1, M_2, \ldots, M_N\}$,
the degrees of freedom of each of them, and an energy 
function $E$, find for
each $M_i$ a set of conformations $\{C_{i1}, C_{i2}, \ldots, C_{ik_i}\}$,
such that for each ligand $i$, $i= 1,\ldots, N$,  
\[ E(C_{ij}) \leq \mbox{THRESHOLD}, \mbox{\ and\ } \mbox{DRMS}(C_{ij}, C_{il}) >
\mbox{TOLERANCE\ }
j,l=1,\ldots, k_i, j\neq l\]
where THRESHOLD and TOLERANCE are pre-specified values, and DRMS is a
distance function.
\end{problem}

\begin{problem} {\bf (Invariant Identification)} 
Given a set of molecules $M = \{M_1, M_2, \ldots, M_n\}$, 
where each $M_i$ is a collection of
conformations
$\{C_{i1}, C_{i2}, \ldots, C_{ik_i}\}$,
determine a
set of labeled points $S$ in $\Re^3$ with the property that for all $i \in
\{1,\ldots,n\}$, there exists some $C_{ij} \in {\cal C}(M_i)$, such
that $S$ is congruent to some subset of $C_{ij}$. The solution $S$, if
it exists, is called an {\em invariant} of $M$.
\end{problem}

\vspace{0.2cm}
Although at this stage the two modules of RAPID work independently, we
plan to support their interaction as the system develops. For example,
the invariant identification module  
will be able to request from the conformational search module 
%more valid conformations of a specific
%molecule, or 
conformations that contain certain features
of the molecule in specified spatial positions.
RAPID also needs to deal with a host of complicating factors.
In practice, the input may contain ligands that do not bind firmly to a
receptor and do not contain a pharmacophore. This requires us to consider a
relaxation of problem 2 above, where a geometric invariant 
need only be present in
conformations of some $K$ of the $N$ molecules.
Furthermore, an invariant may be
perturbed in its spatial layout in distinct ligands' conformations,
or there may be a large number of spurious invariants of differing
sizes that do not correspond to a pharmacophore. Each of these
contributes to a combinatorial explosion in the search process for
the invariants. We defer a detailed discussion of such issues to
the final version of this paper.

%\section{Molecular Modeling}
%\label{sec:problem}
%
%The degrees of freedom of a
%molecule include bond lengths, bond angles
%(angles between two consecutive bonds), and dihedral or torsional angles
%(angles formed by the first and
%third of three consecutive bonds, viewed along the
%axis of the second bond). Bond lengths and bond angles
%do not seem to exhibit large
%changes in their values and are frequently modeled by
%standard geometries calculated from statistical analysis
%of structural data in X-ray databases
%\cite{boyd,guida}.  Torsional angles, however,
%vary significantly and this affects the conformation of a
%molecule.
%
%Valid conformations of a ligand have energy that is
%below a predefined threshold value. For energy calculations,
%empirical force fields are used in practice instead of
%more detailed methods like quantum mechanics.
%While empirical force fields
%offer only an approximation of the energy of the molecule,
%they have nevertheless given correct predictions in a number of cases
%\cite{blaney,boyd,guida}. Force fields include
%terms for bond-stretch, bond-angle, and torsional-angle
%deformations, as well as terms for
%van der Waals and Coulomb potentials \cite{lybrand}.
%Frequently, terms that model solvation effects
%are also present, but since these are difficult to model,
%the molecule is often assumed to be in vacuum \cite{blaney,boyd}.
%
%A conformation of a molecule can be abstracted as a set
%of points in $\Re^3$ which represent atom centers.
%Each point has a label that describes the type of the atom that
%gave rise to it or a feature of that
%atom (e.g., positive charge). It is also possible that
%groups of atoms are aggregated to one point endowed with
%a feature common to all these atoms (e.g., a rigid benzene ring
%may be aggregated to a point). The representation of a
%conformation by points is consistent with the projected
%uses of our results (candidate pharmacophores, database queries)
%since bond information is not relevant in these applications.
%We assume that once a conformation is given, one can automatically
%transform it to a unique collection of points.
%For the purposes of this paper, a conformation $C$
%may represent a valid geometric embedding of a molecule where
%bonds are retained, or a set of points without
%bond information.
%
%In our system, the identification of geometric invariants in a
%collection of flexible ligands is viewed as a two-stage process:
%{\bf (a)} conformational search for each of the ligands under
%consideration; and, {\bf (b)} identification of invariants from the
%point sets arising from the conformations produced. Note that the
%invariant need only be present in at least one conformation
%of most, but not necessarily all, of the ligands; some of the
%ligands under consideration may not bind firmly to a receptor
%and may not contain a pharmacophore at all. We 
%emphasize again that any real system has
%to deal with a host of complicating factors: an invariant may be
%perturbed in its spatial layout in distinct ligands' conformations, 
%or there may be a large number of spurious invariants of differing
%sizes that do not correspond to a pharmacophore. Each of these
%contributes to a combinatorial explosion in the search process for
%the invariants.
%

\section{Related Work}
\label{sec:related}

As far as conformational search is concerned, both systematic and
randomized techniques are being investigated~\cite{blaney,leach,lengauer-gen}.
Systematic search methods sample each torsional
degree of freedom of the ligand at regularly spaced intervals and
minimize all conformations produced~\cite{lipton}.
These techniques can be prohibitively expensive~\cite{leach}.
Several heuristics are employed to quickly prune away
conformations that are ``close'' to previously generated
conformations 
%in an effort to enhance the diversity of the
%sample~
\cite{smellie}. Randomized methods work as follows: 
conformations are obtained by applying
random increments to torsional degrees of freedom of the molecule
starting with a user-specified initial conformation~\cite{ghose},
or with a previously discovered low-energy conformation~\cite{chang}.
Recent articles, which attempt to compare different methods, emphasize the
superior quality of
the results obtained with randomized methods \cite{ghose}.
Other considered methods produce low energy conformations 
obeying distance constraints using  
inverse kinematics \cite{manocha}, algebraic methods \cite{go70}, or 
distance geometry techniques\cite{crippen}.
Let us note here that the protein folding problem is also 
a conformational search problem but its large size 
prohibits the use of the above techniques
(see Dill~\cite{dill}).

Invariant identification is related to well-studied problems in geometric
optimization~\cite{BP,MIS,KMY}. While determining the congruity of two
point-sets in $\Re^3$ is tractable~\cite{BP} in absence of complications
such as noise, invariant identification is more closely related to the
problem of identifying the {\em largest common point set} (LCP).
Unfortunately, the LCP problem turns out to exceedingly difficult; in
fact, even for $m$ collections of $n$ points on the real line, the LCP
cannot be approximated to within an $n^\epsilon$ factor unless $P=NP$,
and only weak positive results are known~\cite{MIS,KMY}. Of course, there
are constraints in the molecular structures, but there are also various
complications as discussed earlier.

In the computational chemistry literature, the most popular algorithms
are based on clique-detection. For instance, the DISCO
system~\cite{martin-pha} initially considers a pair of
conformations $c_1$ and $c_2$ belonging to different molecules and
constructs a ``correspondence graph'' $G$.
The nodes of $G$ are all node pairs, and an edge is created if the
pairs in each of the connected nodes can be matched simultaneously.
The Bron-Kerbosch clique detection
algorithm~\cite{bron} is then used to find cliques in $G$. These
correspond to invariants in $c_1$ and $c_2$, and thus to
candidate  pharmacophores. Maximum clique detection is
NP-hard~\cite{garey}, but this algorithm seems to work
well in practice~\cite{martin-pha,willet}. This approach has been 
generalized  
to $n$ conformations by choosing a reference conformation and comparing
it with the other $n-1$ conformations. 
%Common parts of all pairwise invariants need to be computed at the end. 
Unfortunately, if a large number of
conformations per molecule are to be considered, then there is a
tremendous blow-up in the number of primitive operations performed by
such algorithms~\cite{barnum}, rendering them impractical. 
%A major
%issue with the proposed generalization is that it treats each distinct
%conformation of a molecule as a separate molecule. 
%\mnote{Needs more work}
Given this situation,
several new approaches are under development. For instance, one idea is
to start with small invariants (2-3 features) and gradually expand
them~\cite{barnum}. Our approach is fairly different and involves heavy
use of randomized search techniques~\cite{rand}.

\section{Conformational Search}
\label{sec:conf}

%A randomized technique is used for exploring the conformational
%space of small flexible molecules.
%The method borrows ideas from the success of randomized techniques
%for planning in high-dimensional configuration
%spaces \cite{kavraki-th}.

The goal of conformational search 
as defined in Problem 1, Section \ref{sec:approach} 
is to produce a number of distinct geometric conformations of a ligand that are
below an energy threshold.  

We proceed as follows. Initially a large number 
of conformations are generated at random. 
In contrast with previous randomized search methods, 
we obtain a random conformation by selecting each
degree of freedom from its allowed
range according to a user specified distribution.
This distribution is frequently the uniform distribution.
However, if some a priori information is
available about the preferred values of a particular
degree of freedom~\cite{klebe}, then the corresponding values are selected
according to a distribution that reflects the a priori 
information (i.e. a gaussian distribution).
An efficient minimizer~\cite{brent,pierre} is then used  to obtain
conformations at local energy minima. Minimization is the most
time-consuming step during conformational search, so we have carefully
optimized this procedure.

To obtain a representative set of conformations from or sample,  
we partition it into sets that
reflect geometric similarity as captured by the distance measure DRMS.
We define DRMS($C_i$,$C_j$) as the root mean square distances of the
corresponding atoms of $C_i$ and $C_j$, after $C_i$ is transformed to $C_j$.
This transformation is computed by using a basis of three predefined
atoms $a_1$, $a_2$, and $a_3$ as described in ~\ref{sec:inv}.
We perform a greedy clustering of the conformations by placing a given
conformation in an existing cluster if its distance from the ``center''
of that cluster is less than the predefined value TOLERANCE. If no such cluster is
found, a new cluster is created. The center of a cluster is the
conformation with the lowest energy in the cluster.
%[replace with furthest point heuristic].

We omit further details of the clustering, but note that
the implementation employs heuristics for improving the
quality of output, e.g., it is ensured the the cluster center 
is close to the average conformation in the cluster and
a post-processing step checks that all cluster centers
are sufficiently far apart. 
While REDUCE produces reasonable experimental results,
we are investigating further optimizations via incremental clustering
algorithms \cite{?}.

