\documentstyle[times,psfig,epsf,11pt]{article}

%% USE THIS TO GET MORE INTER-LINE SPACING, IF WE HAVE EXTRA SPACE
\renewcommand{\baselinestretch}{1.24}

\setlength{\topmargin}{0in}
\setlength{\textheight}{9in}
\setlength{\textwidth}{6.5in}
\setlength{\evensidemargin}{0in}
\setlength{\oddsidemargin}{0in}
\setlength{\footheight}{0.0in}
\setlength{\headheight}{0.0in}
\setlength{\headsep}{0.0in}

\newcommand{\omt}[1]{}
\newcommand{\mnote}[1]{{\tiny\marginpar[#1]{#1}}}
\newcommand{\bigvec}[1]{{\stackrel{\rightarrow}{#1}}}
\newsavebox{\theorembox}
\newsavebox{\lemmabox}
\newsavebox{\claimbox}
\newsavebox{\factbox}
\newsavebox{\defnbox}
\newsavebox{\corollarybox}
\newsavebox{\propositionbox}
\newsavebox{\examplebox}
\newsavebox{\exercisebox}
\newsavebox{\problembox}
\newsavebox{\conjecturebox}
\newsavebox{\oqbox}
\newsavebox{\remarkbox}
\savebox{\theorembox}{\noindent\bf Theorem}
\savebox{\lemmabox}{\noindent\bf Lemma}
\savebox{\claimbox}{\noindent\bf Claim}
\savebox{\factbox}{\noindent\bf Fact}
\savebox{\defnbox}{\noindent\bf Definition}
\savebox{\corollarybox}{\noindent\bf Corollary}
\savebox{\propositionbox}{\noindent\bf Proposition}
\savebox{\examplebox}{\noindent\bf Example}
\savebox{\exercisebox}{\noindent\bf Exercise}
\savebox{\problembox}{\noindent\bf Problem}
\savebox{\conjecturebox}{\noindent\bf Conjecture}
\savebox{\remarkbox}{\noindent\bf Remark}
\savebox{\remarkbox}{\noindent\bf Assumption}
\savebox{\oqbox}{\noindent\bf Open Question}
\newtheorem{theorem}{\usebox{\theorembox}}
\newtheorem{lemma}{\usebox{\lemmabox}}
\newtheorem{claim}{\usebox{\claimbox}}
\newtheorem{fact}{\usebox{\factbox}}
\newtheorem{remark}{\usebox{\remarkbox}}
\newtheorem{assumption}{\usebox{\remarkbox}}
\newtheorem{corollary}{\usebox{\corollarybox}}
\newtheorem{proposition}[theorem]{\usebox{\propositionbox}}
\newtheorem{exercise}{\usebox{\exercisebox}}
\newtheorem{problem}{\usebox{\problembox}}
\newtheorem{conjecture}{\usebox{\conjecturebox}}
\newtheorem{oq}{\usebox{\oqbox}}
\newtheorem{defn}{\usebox{\defnbox}}
\def\blackslug{\hbox{\hskip 1pt \vrule width 4pt height 8pt depth 1.5pt
\hskip 1pt}}
\newtheorem{exprime}{Example}%[section]
%\newenvironment{example}{\noindent\exprime\em {\vspace{-0.05 in}}}{
%\hfill $\blackslug$ \par \addvspace{\bigskipamount}{\vspace{0in}}}
\newenvironment{example}{\noindent\exprime\em {\vspace{-0.05 in}}}{
\hfill $\blackslug$ \par \addvspace{\bigskipamount}{\vspace{0in}}}

%Miscellaneous commands for proofs
\newcommand{\qed}{\mbox{}\hspace*{\fill}\nolinebreak\mbox{$\rule{0.7em}{0.7em}$}}
\newenvironment{proof}{\par{\bf Proof:}}{\(\qed\) \par}
\newenvironment{proofsketch}{\par{\bf Proof Sketch:}}{\(\qed\) \par}
\newenvironment{defproof}[1]{\par{\bf Proof #1:}}{\(\qed\) \par}
\newcommand{\case}[1]{{\small \sc Case #1:}}
\newcommand{\subcase}[1]{{\em Case #1:}\/}

\renewcommand{\hat}[1]{\widehat{#1}}
\newcommand{\opt}{\mbox{\footnotesize\sc opt}}
\newcommand{\parent}{\mbox{\sc parent}}
\newcommand{\precstar}{\stackrel{\mbox{\raisebox{-0.5ex}[0pt][0pt]{\tiny$\star$}}}{\prec}}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% ALGORITHM ENVIRONMENT STUFF

\newcommand{\step}[1]{\item[{\bf #1}]}	%macro for steps

  % Boxing command (especially for algorithms)
\long\gdef\boxit#1{\vspace{5mm}\begingroup\vbox{\hrule\hbox{\vrule\kern3pt
      \vbox{\kern3pt#1\kern3pt}\kern3pt\vrule}\hrule}\endgroup}

  %Rajeev's environment for outermost block in an algorithm
\newenvironment{algorithm}[3]%
  { \begin{flushleft}%
   \begin{minipage}[b]{0.95\textwidth}%   %minipage keeps algorithm on same page
      %To allow it to break, use \leftmargin and \rightmargin below.
   \noindent{\bf The #1 Algorithm\/}: \vspace*{-0.1in}%
   \begin{description}
    {\setlength{\listparindent}{-2.5mm}%
     %\addtolength{\leftmargin}{0.1in}%
     \setlength{\parskip}{0ex}%
     \setlength{\itemsep}{0in}
     \setlength{\parsep}{\parskip}%
     \setlength{\topsep}{-0.45in}%
     \setlength{\partopsep}{0in}}%
     \item[{\bf Input:}] #2 \vspace*{-0.1in}
     \item[{\bf Output:}] #3 \vspace*{-0.1in}
   \end{description}
   \begin{description}{}%
    {\setlength{\listparindent}{-2.5mm}%
     %\addtolength{\leftmargin}{0.1in}%
     \setlength{\parskip}{0ex}%
     \setlength{\itemsep}{0in}
     \setlength{\parsep}{\parskip}%
     \setlength{\topsep}{-0.45in}%
     \setlength{\partopsep}{0in}}%
   %\item[]%
  }%
  {\end{description}%
   \end{minipage}%
   \end{flushleft}%
  }

  %Rajeev's environment for indented block in an algorithm
\newenvironment{block}%
  {\begin{list}{}%
     {\addtolength{\leftmargin}{0.1in}%
     \setlength{\listparindent}{-2.5mm}%
     \setlength{\parskip}{0ex}%
     \setlength{\parsep}{\parskip}%
     \setlength{\itemsep}{0in}%
     \setlength{\partopsep}{0in}%
     \setlength{\topsep}{-0.1in}}%
     \vspace*{-0.1in}
     \item[]%
  }%
  {\end{list}%\vspace*{-0.1in}%
  }

  %Algorithm keywords
\newcommand{\Algorithm}{{\bf algorithm\ }}
\newcommand{\Procedure}{{\bf procedure\ }}
\newcommand{\Function}{{\bf function\ }}
\newcommand{\Var}{{\bf var\ }}
\newcommand{\Comment}{{\bf comment:\ }}
\newcommand{\Begin}{{\bf begin\ }}
\newcommand{\End}{{\bf end\ }}
\newcommand{\If}{{\bf if\ }}
\newcommand{\Or}{{\bf \ or\ }}
\newcommand{\And}{{\bf \ and\ }}
\newcommand{\Not}{{\bf not\ }}
\newcommand{\Then}{{\bf then\ }}
\newcommand{\Else}{{\bf else\ }}
\newcommand{\Repeat}{{\bf repeat\ }}
\newcommand{\Times}{{\bf times\ }}
\newcommand{\Forever}{{\bf forever\ }}
\newcommand{\For}{{\bf for\ }}
\newcommand{\While}{{\bf while\ }}
\newcommand{\Until}{{\bf until\ }}
\newcommand{\From}{{\bf from\ }}
\newcommand{\To}{{\bf to\ }}
\newcommand{\By}{{\bf by\ }}
\newcommand{\All}{{\bf all\ }}
\newcommand{\Such}{{\bf such\ }}
\newcommand{\That}{{\bf that\ }}
\newcommand{\Do}{{\bf do\ }}
\newcommand{\Existentially}{{\bf existentially\ }}
\newcommand{\Universally}{{\bf universally\ }}
\newcommand{\Choose}{{\bf choose\ }}
\newcommand{\Guess}{{\bf guess\ }}
\newcommand{\Accept}{{\bf accept}}
\newcommand{\Reject}{{\bf reject}}
\newcommand{\Halt}{{\bf halt}}
\newcommand{\Return}{{\bf return\ }}
\newcommand{\Send}{{\bf send\ }}
\newcommand{\Receive}{{\bf receive}}
\newcommand{\Returns}{{\bf returns\ }}
\newcommand{\Case}{{\bf case }}
\newcommand{\Of}{{\bf of }}
\newcommand{\Let}{{\bf let }}
\newcommand{\Goto}{{\bf go to }}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\renewcommand{\v}[1]{{\bf #1}}


\title{\LARGE \bf RAPID: Randomized Pharmacophore Identification \\ for Drug
Design\thanks{This research is supported by a grant from Pfizer Central
Research.} \\ \mbox{\ }}


\author{
{\Large \sc P.~Finn}\thanks{Pfizer Central Research. {\sf E-mail:}
{\tt finn@sandwich.pfizer.com}.}
\and
{\Large \sc  L. E.~Kavraki}\thanks{Partially supported by startup funds
from Rice University. Currently at Department of 
Computer Science, Rice University, Houston, TX 77005. 
{\sf E-mail:} {\tt kavraki@cs.rice.edu }%\mbox{\ \ \ \ \ \ \ \ \ }
}
\and
{\Large \sc J.-C.~Latombe%\mbox{\ \ \ \ \ \ \ \ \ }
}
\and
{\Large \sc R.~Motwani}\thanks{Partially supported
by an Alfred P.~Sloan Research Fellowship, an IBM Faculty
Partnership Award, an ARO MURI Grant DAAH04-96-1-0007, and NSF Young
Investigator Award CCR-9357849, with matching funds from IBM,
Mitsubishi, Schlumberger Foundation, Shell Foundation, and Xerox
Corporation.}
\and
{\Large \sc C.~Shelton%\mbox{\ \ \ \ \ \ \ \ \ }
}
\and
{\Large \sc S.~Venkatasubramanian%\mbox{\ \ \ \ \ }
}
\and
{\Large \sc A.~Yao}\thanks{On leave from Princeton University. {\sf E-mail:}
{\tt yao@cs.princeton.edu}}
}

\date{
\begin{center}
{\large Department of Computer Science} \\%[0.15cm]
{\large Stanford University} \\%[0.15cm]
{\large Stanford, CA 94305} \\[0.15cm]
{\large \tt
\{latombe,motwani,cshelton,suresh,yao\}@cs.stanford.edu}\\[0.15cm]
\end{center}
}

\begin{document}

\begin{titlepage}

\def\thepage{}

\maketitle

\begin{abstract}
This paper describes a randomized approach for finding invariants
in a set of flexible ligands (drug molecules) that underlies an integrated
software system called RAPID currently under development. An invariant
is a collection of features embedded in $\Re^3$ which is present in one
or more of the possible low-energy conformations of each ligand.
Such invariants of chemically distinct molecules are useful for
computational chemists since they may represent candidate
pharmacophores. A pharmacophore contains the parts of the ligand
that are primarily responsible for its interaction and binding with
a specific receptor. It is regarded as an inverse image of a receptor
and is used as a template for building better pharmaceutical drugs.
The identification of pharmacophores is crucial in drug design
since the structure of the targeted receptor is frequently unknown,
but a number of molecules that interact with the receptor have been discovered
by experiments. It is expected that our techniques and the
results produced by our system will prove useful in other applications
such as molecular database screening and comparative molecular field analysis.
\end{abstract}
\end{titlepage}

\section{Introduction}

Computational chemists working in the area of 
structured-based drug design consider 
both chemical and 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 results showing  
molecules with geometric and chemical
complementarity in their active, or binding, conformations~\cite{boyd}.

When the three dimensional structure of the receptor is known,
docking methods~\cite{blaney} 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 discovered 
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 some active conformation of each (or most) of the 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{boyd}.
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 
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}.
In this abstract, we focus on the latter module and describe in some
detail a randomized approach for solving an abstract geometric problem
lying at the core of pharmacophore identification -- the search for
common subsets in a collection of point sets embedded in $\Re^3$.
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}
}
\vspace{-0.3cm}
\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{boyd} 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{boyd}. 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{martin-pha}.
%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$. For each $M_i$, find a set of conformations
$\{C_{i1}, C_{i2}, \ldots, C_{ik_i}\}$, such that for each
$i= 1,\ldots, N$ and $j = 1,\ldots, k_i$,
\[ E(C_{ij}) \leq \mbox{THRESHOLD}, \mbox{\ and\ } \mbox{d}(C_{ij},
C_{il}) > \mbox{TOLERANCE\ } \mbox{for\ } l \neq j
\]
where THRESHOLD and TOLERANCE are pre-specified values, and $d(.,.)$ is a
distance function.
\end{problem}

\begin{problem} {\bf (Invariant Identification)} 
Given a collection of ligands $M = \{M_1, M_2, \ldots, M_N\}$, 
where each $M_i$ is a set of conformations
${\cal C}(M_i) = \{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}$. A 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 the 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} and 
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, 
which obey distance constraints, using  
inverse kinematics \cite{manocha}, algebraic methods \cite{go70}, or 
distance geometry techniques~\cite{crippen}.
Note that the {\em 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 the well-studied problem
in geometric optimization~\cite{BP,ATT,MIS,KMY} of finding
common point sets. While determining the congruence of two
point-sets in $\Re^3$ is tractable~\cite{BP,ATT} {\em in the 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 the
problem is polynomially solvable when $m=2$~\cite{ATT} but that is not
good enough for our purposes. Also, 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, DISCO
\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.
A 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}. The 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, 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
low-energy conformations of a given ligand.
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 our 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 square root of the  mean of the 
squared 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$ (see Section~\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 cluster center 
is close to the average cluster conformation and
a post-processing step checks that all cluster centers
are sufficiently far apart. 
While we obtain reasonable experimental results,
we are investigating optimizations via incremental clustering
algorithms \cite{MCC}.


\section{Identification of Invariants}
\label{sec:inv}

The set of cluster centers, denoted by 
${\cal C}(M) = {\cal C}(M_1) \cup \ldots \cup {\cal C}(M_N) $
is the input for the invariant identification module. 
%Let ${\cal C}(M)$ denote the set of low-energy conformations for a molecule $M$. 
%In particular, ${\cal C}(M)$ is represented as a collection
Each conformation in ${\cal C}(M)$ is now represented as 
a set of labeled points in $\Re^3$ (see Section \ref{sec:approach}).
%   This is obtained by extracting the 
%labels from each conformation.  
We wish to determine a structure $S$ that is
congruent to a substructure of some conformation in every molecule.
The congruence relation is with respect to 3-D rotations and
translations that preserve equality of labels.

%\begin{problem}
%Given a set of molecules $M = \{M_1, M_2, \ldots, M_n\}$, 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}

%Note that the congruence relation is with respect to 3-D rotations and
%translations that preserve equality of labels.
Our formulation of the pharmacophore identification problem in
Problem 2 (Section~\ref{sec:approach})
ignores the issue of error in the data. In
practice, atom positions are fuzzy, and it may not be possible to align
points exactly. Therefore, we adopt the convention that 
two points $p_1$ and $p_2$
are said to {\em match} when $\mid p_1 - p_2 \mid \leq \epsilon$,
where $\epsilon$ is the {\em point location error}.

The pharmacophore identification 
problem is a variant of the {\em largest common point set problem
(LCP)} in three dimensions, which is the following. Given $s$ point
sets $P_1, P_2, \ldots, P_s$ in $\Re^3$, determine the point set of
maximum cardinality congruent to some subset of each point set.
For convenience, we assume that each point set $P_i$ has cardinality
exactly $n$. For arbitrary $s$ and dimension $d$,
LCP is hard to approximate within a factor of $n^\epsilon$,
for some $\epsilon \ge 0$~\cite{MIS}. In the sequel, we consider the following
version of LCP, called LCP-$\alpha$: determine a point set $S$ of size $|S| \geq
\alpha n$ congruent to some subset of each $P_i, 1 \leq i \leq s$.
%where $n = \min_{i=1}^s |P_i|$.
The motivation for focusing on this subproblem is that it better
captures the application where it is known that a pharmacophore has
to be a reasonably large invariant.
We begin by focusing on the case $s = 2$, and later indicate how we
generalize to arbitrary values of $s$.


\subsection{Phase 1: Pairwise Matching}
In this section, we focus on the invariant identification problem for
two point sets, denoted by MATCH. This problem has been studied extensively
in the literature~\cite{ragh,pr}. For the case when $\alpha = 1$, i.e.,
determining the congruence of two point sets, Atkinson~\cite{ATK} obtains
a solution in time $O(n \log n)$ in two and
three dimensions. However, for general $\alpha$, the problem is much
harder, relating to combinatorial distance properties of $\Re^3$ and is
examined by Akutsu, Tamaki, and Tokuyana~\cite{ATT}, in 
which a worst-case bound of $O(n^{4.6})$ is
obtained for unknown $\alpha$, and $O(n^{2.6}/\alpha^2)$ randomized when
$\alpha$ is known. For two dimensions, the corresponding bounds are
$O(n^{3.2})$ and $O(n^{2.2}/\alpha)$.

In this section, we describe two random-sampling schemes for solving
LCP-$\alpha$. The first is quite natural, and the second uses more careful
random sampling. Under certain assumptions about the distribution on
distances between points in a set, we prove loose bounds on the running
time of these schemes, illustrating the superiority of the second
approach. Our results are confirmed by experimental evidence in
Section~\ref{sec:exp}.

For a set of points $S$, let $D_S$ be the multi-set of all inter-point
distances in $S$. For a given {\em cluster width} $\delta$, we can
partition $D_S$ into $L$ {\em buckets}, where  bucket $i$ contains all
entries $d \in D_S \mbox{\ such that \ } i\delta \le d < (i+1)\delta$.
For the purposes of theoretical analysis, we assume that given a point set
$P$, each such distance bucket is of constant size, for fixed $\delta$;
our experimental results will later justify this assumption.
The notation $g(n) = \tilde{O}(f(n))$, where $f$ and $g$ are functions, implies
that $g(n) = O(f(n) \log n)$. Note that in three dimensions, a unique
transformation $T$ (upto reflection) between two point sets $P_1$ and $P_2$
is determined by matching three points $p, q$, and $r$ in $P_1$ with
three points $s, t$, and $u$ in $P_2$. Furthermore, is it known that~\cite{DL}:

%\begin{fact}
%In three dimensions, a unique transformation $T$ (upto reflection) between two
%point sets $P_1, P_2$ is determined by three points $p, q, r \in P_1$, and
%$s,t,u \in P_2$. 
%\end{fact}

\begin{proposition}
Given a transformation $T$ from $P_1$ to $P_2$, in time $\tilde{O}(n)$ we
can determine for each point $p \in P_1$ a corresponding point $p' \in P_2$
such that $|T(p)-p'|$ is minimized.
\end{proposition}

The basic random sampling method is as follows.

\paragraph{BASIC-SAMPLE:} 
For some constant $c$, perform $c/\alpha^3$ iterations of the 
following sampling process: sample a triplet of points $<p_1, p_2, p_3>$
randomly from  $P_1$; determine three points in $P_2$ congruent to this set;
compute the resulting induced transformation and determine the number of
points in $P_1$ matching corresponding points in $P_2$; and, if the number
exceeds $\alpha n$, declare SUCCESS.

\begin{theorem}
Given a common subset $S$ of size $|S| \ge \alpha n$, the probability that
BASIC-SAMPLE fails to declare SUCCESS is at most $1/n$.
\end{theorem}

The proof of this theorem is straightforward. Observe that the probability
of picking three points that belong to the common substructure is at least
$1/\alpha^3$.

\begin{theorem}
BASIC-SAMPLE runs in time $\tilde{O}(n/\alpha^3)$ using space $O(n^2)$.
\end{theorem}

The running time can be written as 
$T = n \sum_{i=1}^{1/\alpha^3} \Delta_i $ 
where $\Delta_i$ is the number of candidate transformations generated
during iteration $i$. Using the notation of \cite{ATT}, let $H_P(t)$ denote
the number of triangles congruent to $t$ in point set $P$. It is the
case that  $\Delta_i = H_{P_1}(t_i)*H_{P_2}(t_i)$, where $t_i$ is the
triangle induced by the triplet chosen at iteration $i$.
>From our assumption on the size of distance buckets, we produce a
constant number of transformations in each stage, and the above bound
holds.  

%Figure~\ref{fig:search} shows experimental results indicating that the
%number of candidate transformations generated indeed behaves as
%$O(1/\alpha^3)$. 

Besides BASIC-SAMPLE, we consider the following more sophisticated sampling
scheme.

\paragraph{PARTITION-SAMPLE:}
For some constant $c$, perform $c$ iterations of the following sampling
process: randomly select two subsets $A$ and $B$ of size $1/\alpha$ from
$P_1$; also select a subset $C$ of size $1/\alpha$ from $P_2$; store all
distances $d(p,q)$, for all $p \in C$ and $q \in P_2-C$, in a hash table; 
for every triangle $(a, b, p)$ with $a \in A, b \in B$, and $p \in
P-(A \cup B)$, probe for $d(p,a)$ and $d(p,b)$ in the hash table to
determine all matching triplets $(c, p_1, p_2)$ with $c \in C$ and
$p_1, p_2 \in P_2-C$; finally, as before, if the resulting transformation
induces a match of more than $\alpha n$ points, declare SUCCESS.

\begin{theorem}
Given a common subset $S$ of size $|S| \ge \alpha n$, the probability
that PARTITION-SAMPLE fails to declare SUCCESS is at most $e^{-c\alpha n}$.
\end{theorem}

We defer the proof of the preceeding theorem to the final version of the paper. 

\begin{theorem}
PARTITION-SAMPLE runs in time $\tilde{O}(n^2/\alpha^2-n/\alpha^3)$ using space
$O(n/\alpha)$. 
\end{theorem}

As before, the running time of the algorithm can be written as 
$ T = n \sum_{i=1}^{n/\alpha^2-1/\alpha^3} \Delta_i$. If our assumption
uniform hashing assumption is valid, then the bound follows.
Figure~\ref{fig:search} illustrates experimental indications that
this is the case. Intuitively, this partitioning favors ``large''
solutions over small ones. Given a solution $S$ of size $K$ in $P_1$,
the probability that each of $A, B$, and $P_1-A-B$ contain at least one
point of $S$ increases with $K$.

%In Section\ref{sec:exp}, we show experimental results supporting 
%our claims. For example,  
%Figure~\ref{fig:search} indicates that the  
%number of candidate transformations generated by BASIC-SAMPLE indeed behaves as
%$O(1/\alpha^3)$, while Figure~\ref{fig:compare} shows that 
%that PARTITION-SAMPLE
%consistently outperforms BASIC-SEARCH.


\paragraph{Eliminating Redundant Solutions.}
Any MATCH algorithm could return many solutions that satisfy a containment
relationship with respect to each other. To prevent such redundant
solutions from propagating through the complete search procedure, we need
to {\em prune} them. Given invariants $S_1$ and $S_2$, we can check if $S_1
\subseteq S_2$ by invoking MATCH on these two sets with $\alpha$ set to
1. Clearly, since $\alpha$ is 1, a solution will be produced only if
$S_1$ is contained in $S_2$. This pruning is done in a pair-wise fashion on the
set of invariants before these invariants are used as input to
the next phase.

\paragraph{Inaccuracies in Computing the Transformation.}
The module at the core of the substructure search algorithm is the
following: Given points $p_1, p_2, p_3  \in P_1$, and points $p'_1, p'_2,
p'_3 \in P_2$, compute a transformation $T$ from $P_1$ to $P_2$ that maps
$p_i \to p'_i$, for $i = 1,2,3$. For our 
purposes we need to compute a transformation that,
given this initial set of correspondences, matches the largest number of
points from $P \mbox{\ to\ }Q$. Such a transformation can be computed in
many different ways. The simplest such way is: {\bf (i)} align $p_1$ and
$q_1$; {\bf (ii)} align the vectors $\bigvec{p_1p_2}$ and $\bigvec{q_1q_2}$;
and, {\bf (iii)} align $p_3$ and $q_3$ by a rotation about $\bigvec{p_1p_2}$.

This approach over-constrains the fit between the two triangles.
In general, it is sufficient that these pairs match
within distance $\epsilon$, at such an orientation that the largest number
of points 
from $P_1$ can be made to match points in $P_2$. However, we do not know these
other points in advance, and as a result of over-constraining, it is possible
that we may not be able to force as many points to match as we should.

If we knew the correspondences between points in $P_1 \mbox{\ and\ }P_2$,
then a more ``balanced'' transformation could be obtained by computing the
transformation minimizing $\mbox{DRMS}(P'_1, P'_2), P'_1\mbox{\ and\ } P'_2
$ being the subsets of $P_1$ and $P_2$ that correspond~\cite{A87} (see
Section \ref{sec:conf} for the definition of DRMS). This suggests
the obvious refinement: Use the procedure above to generate a plausible
correspondence list, and then use a DRMS minimizing procedure to compute a
refined transformation.

However, this approach is more time-consuming, and it is possible that DRMS
minimization merely replaces one form of over-constraining by another. In
general, if we consider the six-dimensional space of all rigid
transformations in $\Re^3$, the set of transformations satisfying the given
correspondences is actually an enclosed volume. Choosing any
transformation from within this volume will produce an acceptable
transformation between the two given point sets. However, when we
generate candidate solutions and attempt to match them against other point
sets, the ``right'' choice of transformation could be crucial in
preventing points from being lost. Our implementation employs various
heuristics to try and minimize the effect of this problem. One such
heuristic is the following.
Determine a {\em seed} transformation $T$ by any of the
above methods, and then sample three random pairs from the set of
correspondences that $T$ induces, using these pairs to construct a new
transformation. Clearly, in a perfect world, we will obtain $T$
again. However, given the inaccuracies in point location, it turns out that
some choices of triplets may yield a larger number of correspondences than
before. We defer experimental results in this regard to a final version
of the paper.

\subsection{Phase 2: Multiple Matching}

The previous phase  yields methods for solving LCP on two sets. 
In Phase 2, candidate solutions obtained from the previous phase
are tested against the remaining molecules to determine the invariant.
Each MATCH call operates on two conformations. Since each molecule is
represented by many conformations, we extend MATCH to two molecules by
doing all pair-wise matches between the sets of conformations.

There are various strategies one could use to process multiple
molecules. A simple strategy  does
the obvious linear merging. We take each solution and compare it with the
next molecule. We do this for all current solutions, concatenate and prune
the results, and repeat with a new molecule.
Another approach, which is in some ways less order-sensitive, is a
tree-based merge. Here, we run the two-molecule algorithm on distinct
pairs of molecules, and recursively combine the results in a binary
tree-like fashion.
Both of these approaches are simple and are correct, if point location
error is not permitted. 

In addition to this, we may wish to find an invariant that does not exist
in {\em all} 
the molecules, but in some fixed number of them. Our current strategy
would fail to do this, because all candidate solutions are compared
against every molecule. We use a modified merging strategy here to  
keep  track of the number of times an invariant fails to match against a
molecule, and only reject those which exceed the maximum allowed number
of failures.

If an invariant $S$ is contained in a point set $P$, then a call to
MATCH($S, P$) should return a set of solutions containing
$S$. This can easily be checked, which therefore means that we can verify
whether an invariant matches against a conformation. Our
marking algorithm associates {\em marks} with each candidate
solution. All original conformations are initialized with zero marks.
All invariants that result from a MATCH acquire the sum of the marks of the
input structures.
Every time a solution fails to match against a new molecule
(which means that it fails to match with {\em any} conformation of that
molecule), we add a mark to it and continue to propagate it. 
We reject any solutions for which
the number of marks exceeds the prescribed maximum.

%For the sake of brevity, we omit detailed experimental analysis of these
%%strategies, Table \ref{CRS-COMPARE} uses a linear search to compute
%invariants, with no failures in matching invariants permitted.


%%The
%%problem is that the first two molecules decide the complete set of
%%candidate solutions which are then pruned and refined. A sequence that
%puts the most similar molecules first will yield solutions that are
%large and can survive subsequent refinement without being discarded
%completely. However, if two less similar molecules are placed first, the
%solutions produced are smaller, and more prone to subsequent removal. Note
%that in a perfect world, all molecules are equally similar (with
%respect to an invariant) to each other, and the order is immaterial.

%The tree search method, on the other hand, is less perturbed by the order
%because of the more global nature of the search, which yields many more
%candidate solutions initially. However, the approach does work better if
%more similar molecules are paired together initially.

%%In the original linear merge, the list of current invariants is updated
%%as follows:
%%\[ INV_{i+1} = \mbox{MMATCH}(INV_i, M_{i+1}).\]
%%The modified linear merge procedure does the following:
%%\[ INV_{i+1} = \mbox{MMATCH}(INV_i, M_{i+1}) \cup \mbox{MMATCH}(M_i,
%%M_{i+1}) \cup INV_i \oplus 1 \]
%%where $INV_i \oplus k$ is defined as
%%$\forall I \in INV_i, \mbox{mark}(I) \leftarrow \mbox{mark}(I)+k.$
%%Every invariant produced by $\mbox{MMATCH}(M_i, M_{i+1})$ is
%%initialized with $i-1$ marks.

%%Analogously, for the tree merge let the $i$th set of invariants produced
%at level $k$ of the merge-tree be denoted as $INV^k_i$. If we define $INV^0_i = M_i, \forall i
%\in \{1,\ldots, n\}$, then the list of current invariants is updated as
%follows:
%\[ INV^k_{i} = \mbox{MMATCH}(INV^{k-1}_{2i}, INV^{k-1}_{2i-1}).\]
%Let $N(INV^k_i)$ be the number of molecules that are leaves of the
%merge-tree whose root is $INV^k_i$. We can write the procedure for
%updating the marked list of current invariants as follows:
%\[ INV^k_{i} = \mbox{MMATCH}(INV^{k-1}_{2i}, INV^{k-1}_{2i-1}) \cup
%	       (INV^{k-1}_{2i} \oplus N(INV^{k-1}_{2i-1})) \cup
%	       (INV^{k-1}_{2i-1} \oplus N(INV^{k-1}_{2i})).\]

%%\begin{algorithm}{MARK-SEARCH}{$M_1, M_2, \ldots, M_n, K$}{Invariants that
%%occur in at least $K$ of the molecules}

%%\begin{enumerate}
%%\item \Let ${\cal S}$ = MMATCH($M_1, M_2$)
%%\item $\forall S \in {\cal S}$, Mark($S$) $\leftarrow 0$
%%\item \For $i = 3 \ldots n $
%%\item 	${\cal R} =$ MMATCH($M_{i-1}, M_i$)
%%\item	$\forall R \in {\cal R}$, Mark($R$) $\leftarrow i-2$
%%\item	${\cal S} \leftarrow {\cal S} \cup {\cal R}$
%%\item	${\cal S}' \leftarrow$ MMATCH(${\cal S}, M_i$)
%%\item	$\forall S \in {\cal S}$,
%%\item	    \If $S$ fails to match with $M_i$,
%%\item	         \Then Mark($S$)++;
%%\item	               \If Mark($S$) $\leq K$,
%%\item	                   \Then add $S$ to ${\cal S}'$
%%\item 	${\cal S} \leftarrow {\cal S}'$
%%\end{enumerate}
%%\end{algorithm}

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

The section reports preliminary experimental results; more 
detailed results and analysis will be included in the final 
version of the paper. All reported timings are on an SGI Indigo2
with a 175 MHz MIPS R8000 processor and 384MB RAM. Code was written in
C/C++, and compiled using SGI CC. 

\vspace{-0.2cm}
\paragraph{Input.} 
In Figure~\ref{fig:mols}, we show four different inhibitors of the
protease thermolysin. These molecules fit into the same cavity of
thermolysin and by their presence inhibit the activity associated
with that cavity. This example was chosen because all the 
inhibitors have been crystallized with thermolysin and their 
active conformations are known and recorded in the PDB database. 
Note that 1TLP has 69 atoms and 10 torsional
degrees of freedom, 4TMN has 68 atoms and 15 degrees of freedom, 5TMN
has 64 atoms and 13 degrees of freedom, and 6TMN has 63 atoms and 12
degrees of freedom.

\begin{figure}
\centerline{
\psfig{figure=4mols.ps,height=3cm}
}
\vspace{-0.4cm}
\caption{1TLP, 4TMN, 5TMN, and 6TMN are inhibitors of thermolysin.}
\label{fig:mols}
\end{figure}

\begin{figure}
\centerline{
\setlength{\unitlength}{1mm}
\begin{picture}(100,38)
\put(-30,0){\psfig{figure=c1.ps,height=3.4cm}}
\put(19,0){\psfig{figure=c2.ps,height=3.4cm}}
\put(55,0){\psfig{figure=c4.ps,height=3.4cm}}
\put(85,0){\psfig{figure=c5.ps,height=3.4cm}}
%\put(100,5){\psfig{figure=c5.ps,bbllx=11mm,bblly=240mm,bburx=50mm,bbury=270mm}}
\end{picture}
}
\vspace{-0.3cm}
\caption{Different clusters of 1TLP.}
\label{fig:clusters}
\end{figure}



\vspace{-0.2cm}
\paragraph{Conformational Search.}
Each of the molecules in Figure~\ref{fig:mols} was run through
our conformational search software. A THRESHOLD value of 20 Kcals/mol was used
for the energy of the valid conformations. It took 4.3 hours to 
to produce 4000 valid conformations for each molecule and reduce
these to a set of representatives.
The timings reported are very reasonable if one takes into 
account that 10, 15, 13, and 12 dimensional spaces are searched.
The number of representatives produced was 850, 1192, 1024, and 955 for 
1TLP, 4TMN, 5TMN, and 6TMN, respectively.
The TOLERANCE value was set to 1.23 $\AA$ for all runs. 
A few clusters of 1TLP are shown in Figure~\ref{fig:clusters}. 
Since we know the active conformations of all the molecules, we 
performed a consistency check at the end of conformational search and 
we confirmed that the active conformations were 
close to one of the produced representatives.

An important issue in conformational search is to decide the number
of valid conformations to produce. At this stage, this number is 
determined experimentally. We plan to interleave the generation of 
conformations and their clustering in the future and  
stop producing new conformations when these fall into existing
clusters and do not increase significantly the overall 
number of clusters. 

\vspace{-0.2cm}
\paragraph{Identification of Invariants.}

%For the sake of brevity, we omit detailed experimental analysis of these
%%strategies, Table \ref{CRS-COMPARE} uses a linear search to compute
%invariants, with no failures in matching invariants permitted.

To find the invariants in these four molecules, the conformations produced by 
the first stage were given as input to the search algorithm described in 
Section~\ref{sec:inv}. For these tests, we considered each non-hydrogen
atom to be a separate feature or point. Thus, each conformation had approximately
$30$ features drawn from $6$ feature classes. We define the ``solution''
to be the overlapping portions of the molecules when aligned as shown in
Figure~\ref{fig:plot}. This is the upper right hand
T-shaped portion of this diagram.  The entire invariant
consists of roughly $7$ atoms and an additional $7$ atoms of
``scaffolding,'' or connecting atoms with no pharmacophore functionality.

Table~\ref{CRS-PARM} details the parameters used for invariant identification.
The parameters $\epsilon$, $\delta$, and $\alpha$ are as described in
Section~\ref{sec:inv}. We set at 5 the number of times the MATCH algorithm
is run per comparison, or the maximum number of invariants found per MATCH.
%The value $n$ is the number of times the MATCH algorithm is run per
%comparison, or the maximum number of invariants found per MATCH.
These parameters
were chosen by experimentation. For example, $\delta$ was selected to
be much lower than $\epsilon$. We found that higher values of $\delta$
increased the running time without improving the quality of the solution.


In Section \ref{sec:inv} we assumed that the distance distribution for
molecules satisfies certain properties, and showed certain performance
bounds based on this. Here, we show experimental validation of these
bounds.
For the sake of brevity, we omit detailed experimental analysis of 
variants of Multiple Search and use a linear search to compute invariants, 
with no failures in matching invariants permitted.

Figure~\ref{fig:search}(a) plots the number of candidate transformations
produced ($\sum_{i=1}^{n/\alpha^2} \Delta_i$) by using the procedure 
PARTITION-SAMPLE 
against
$\alpha$. Superimposed on this is a plot of $n/\alpha^2-1/\alpha^3$.
We do the same for BASIC-SAMPLE in Figure \ref{fig:search}(b). 
Here we superimpose the plot of $1/\alpha^3$.
For both
algorithms, the experimental results follow the predicted behavior.
In Figure~\ref{fig:compare}, we compare the two schemes experimentally.
The comparison 
clearly shows that PARTITION-SAMPLE performs 
better experimentally, although
its asymptotic performance seems worse.

This discrepancy is partly explained in Figure~\ref{fig:4mols}. 
For each molecule in the input set, we plot its distance
distribution function (which BASIC-SAMPLE samples from), as well the
distribution function induced by the subset of distances that PARTITION-SAMPLE
samples (for varying values of $\alpha$). 

We finally present in Table \ref{CRS-COMPARE} a
comparison between the two schemes on large data sets, 
where we vary the number of conformations per molecule. We include
the active conformation in our sets. In all cases, 
the quality of solutions (in terms of the largest solution found) is
comparable, and PARTITION-SAMPLE consistently runs faster than BASIC-SAMPLE.
When the number of conformations increases, more invariants
are produced  because some of the added
conformations had additional ``scaffolding'' which also could be matched.
Yet, in every case, the invariant of the largest size was the ``correct
solution.'' Many non-trivial solutions of size $8-14$ were also produced.

%produced are given in Table \ref{CRS-COMPARE}.
%In all cases, the correct invariant was found.  In trials with more
%conformations, larger sized solutions were produced because some of the added
%conformations had additional ``scaffolding'' which also could be matched.
%Yet, in every case, the invariant of the largest size was the ``correct
%solution.'' Many non-trivial solutions of size $8-14$ were produced, as
%well a large number of trivial solutions of size $4$, 

\vspace{-0.2cm}
\section{Discussion}
\label{sec:disc}

We briefly discuss some enhancements currently underway that we hope
to report in the final version of the paper. 
%We are performing more
%extended experiments to better evaluate the 
%strengths and deficiencies of RAPID. 
Our goal is to optimize the modules of the system  
to perform experiments that 
involve 5-20 ligands and a large number of conformations per ligand.
As far as conformational search is concerned, 
we are investigating dynamic clustering techniques.
By looking at the relative changes in 
the number of clusters found, we hope that we will be able 
to automatically end the conformational search stage. 
For the invariant identification stage, we will try to 
reduce the dependence of the implementation on the ordering of 
the input, a phenomenon observed in the current system. 
Last but not the least, our short-term plans include the  
implementation of an algorithm that will extract features from conformations 
in the way the computational chemists understand them. Currently, 
each atom gives rise to a feature that bears as a label the
corresponding atom type. Computational chemists, however, have rules to 
group atoms into features that describe the chemical behavior 
of parts of the molecule, e.g., hydrophobic parts and positively  
charged parts. Features will reduce the number of points we
consider per conformation and will hopefully  
help us deal with larger examples. 


\newpage

\renewcommand{\baselinestretch}{1}
\normalsize
\pagenumbering{roman}

\begin{figure}
\centerline{
\psfig{figure=active.ps,height=3.4cm} }
\vspace{-0.3cm}
\caption{1TLP, 4TMN, 5TMN, and 6TMN overlapped in their active
conformations.}
\vspace{1.5cm}
\label{fig:plot}
\end{figure}


\begin{table}
\begin{center}
\vspace{-0.2cm}
\begin{tabular}{|c|l|l|} \hline
Parameter	& Search Value	& Prune Value 	\\ \hline
$\epsilon$	& $1.3$		& $3.5$		\\ 
$\delta$	& $0.5$		& $1.0$		\\ 
$\alpha$	& $0.3$		& $1.0$		\\ 
%$n$		& $5$		& $5$		\\ 
\hline
\end{tabular}
\end{center}
\vspace{-0.4cm}
\caption{Parameters for invariant search trials.}
\vspace{1.5cm}
\label{CRS-PARM}
\end{table}

\begin{figure}
\centerline{
\setlength{\unitlength}{1mm}
\begin{picture}(100,58)
\put(-30,5){\psfig{figure=hs.ps,height=5cm}}
\put(55,5){\psfig{figure=rs.ps,height=5cm}}
\put(10,0){(a)}
\put(93,0){(b)}
\end{picture}
}
\vspace{-0.3cm}
\caption{Experimental results vs. theoretical predictions.}
\label{fig:search}
\end{figure}

\newpage 

\begin{figure}
\centerline{
\psfig{figure=hr.ps,height=5.5cm} }
\vspace{-0.3cm}
\caption{Relative performance of BASIC-SAMPLE and PARTITION-SAMPLE.}
\label{fig:compare}
\end{figure}

\begin{figure}
\centerline{
\setlength{\unitlength}{1mm}
\begin{picture}(100,115)
\put(-30,5){\psfig{figure=1tlp.ps,height=5cm}}
\put(10,0){(c)}
\put(55,5){\psfig{figure=4tmn.ps,height=5cm}}
\put(93,0){(d)}
\put(-30,65){\psfig{figure=5tmn.ps,height=5cm}}
\put(10,60){(a)}
\put(55,65){\psfig{figure=6tmn.ps,height=5cm}}
\put(93,60){(b)}
\end{picture}
}
\vspace{-0.3cm}
\caption{Distance distribution functions of 
(a) 5TMN, (b) 6TMN, (c) 1TLP, (d) 4TMN.}
\label{fig:4mols}
\end{figure}

\newpage 

{\small
\begin{table}
\begin{center} \begin{tabular}{|c|c|l|l|l|l|l|l|l|l|l|l|l|l|l|} \hline
Number of       & Time to       &
\multicolumn{13}{|c|}{Histogram of Solution Sizes}      \\  \hline
Conformations   & Completion (in s)
& 4 & 5 & 6 & 7 & 8 & 9 & 10 & 11 & 12 & 13 & 14 & 15 & total\\ \hline \hline
\multicolumn{15}{|c|}{BASIC-SAMPLE} \\ \hline
1& 4.62 &   & 2  & 1  &   &   &   &    &  & & & 1 &  &    \\ \hline
11& 2127.66 & 90  & 34  & 20  & 5 & 6 & 4  & 1 & 1 & 1 &  & & &  \\ \hline
21& 24175.42 & 289  & 110  & 46  & 35  & 42  &18 & 3 &6 & 2 & 2
&  1& &  \\ \hline
\multicolumn{15}{|c|}{PARTITION-SAMPLE} \\ \hline
1& 2.85 &   &  & 1  & & 1 & & & & & 1 &  &  &   \\ \hline
11& 656.47  & 20  & 9  & 8 & 11  &  11  & 3  & & & &  & & &\\ \hline
21& 8060.00  & 73 & 55  & 53  & 30  & 21  & 25 & 7 & 8 & 1 &
1 & 4 & 1 &  \\ \hline
\end{tabular}
\end{center}
\vspace{-0.3cm}
\caption{Comparison of BASIC-SAMPLE and PARTITION-SAMPLE}.
\label{CRS-COMPARE}
\end{table}
}


\iffalse
\section*{Acknowledgment}
This research is supported by a grant from Pfizer Central Research.
%Dr.~Paul Finn of Pfizer Central Research has provided invaluable help
%in identifying and understanding the problems discussed in this paper
%as well as many insightful comments on this work.

Rajeev Motwani's research is partially supported
by an Alfred P.~Sloan Research Fellowship, an IBM Faculty
Partnership Award, an ARO MURI Grant DAAH04-96-1-0007, and NSF Young
Investigator Award CCR-9357849, with matching funds from IBM,
Mitsubishi, Schlumberger Foundation, Shell Foundation, and Xerox
Corporation. Lydia Kavraki is partially supported by startup funds from 
Rice University.
\fi



{\small
%\bibliography{bio}
%\bibliographystyle{abbrv}

\begin{thebibliography}{10}

\bibitem{BP}
T.~Akutsu.
\newblock On determining the congruence of point sets in higher dimensions.
\newblock {\em Lecture Notes in Computer Science 834}, page~38.

\bibitem{MIS}
T.~Akutsu and M.~{Halld\'{o}rsson}.
\newblock On the approximation of largest common subtrees and largest common
  point sets.
\newblock {\em Lecture Notes in Computer Science 834}, pages 405--413.

\bibitem{ATT}
T.~Akutsu, H.~Tamaki, and T.~Tokuyama.
\newblock Distribution of distances and triangles in a point set and
algorithms for computing the largest common point set.
\newblock Manuscript, 1996.

\bibitem{A87}
K.~S.~Arun, T.~S.~Huang, and S.~D.~Blostein.
\newblock Least-squares Fitting of Two 3-D Point Sets.
\newblock {\em IEEE Trans. Pattern. Anal. Machine Intell.}, volume~PAMI-9,
No.~5, September 1987.

\bibitem{ATK}
M.~D.~Atkinson.
\newblock An optimal algorithm for geometric congruence.
\newblock {\em J. Algorithms} 8 (1987). pages 134--183.

\bibitem{balbes}
L.~Balbes, S.~Mascarella, and D.~Boyd.
\newblock A perspective of modern methods in computer-aided drug design.
\newblock In K.~Lipkowitz and D.~B. Boyd, editors, {\em Reviews in
  Computational Chemistry}, volume~5, pages 337--370. VCH Publishers, 1994.

\bibitem{barnum}
D.~Barnum, J.~Greene, A.~Smellie, and P.~Sprague.
\newblock Identification of common functional components among molecules.
\newblock To appear in {\em J. Chem. Inf. Comput. Sci.}, 1996.

\bibitem{blaney}
J.~Blaney and S.~Dixon.
\newblock A good ligand is hard to find: Automated docking methods.
\newblock {\em Perspectives in Drug Discovery and Design}, 1:301--319, 1993.

\bibitem{boyd}
D.~B. Boyd.
\newblock Aspects of molecular modeling.
\newblock In K.~Lipkowitz and D.~B. Boyd, editors, {\em Reviews in
  Computational Chemistry}, volume~1, pages 321--351. VCH Publishers, 1990.

\bibitem{brent}
R.~Brent.
\newblock {\em Algorithms for finding zeros and extrema of functions without
  calculating derivatives}.
\newblock PhD thesis, Stanford University, 1971.

\bibitem{bron}
C.~Bron and J.~Kerbosch.
\newblock Finding all cliques of an undirected subgraph.
\newblock {\em Commun. ACM}, 16:575--577, 1973.

\bibitem{chang}
G.~Chang, W.~Guida, and W.~Still.
\newblock An internal coordinate monte-carlo method for searching
  conformational space.
\newblock {\em J. Am. Chem. Soc.}, 111:4379--4386, 1989.

\bibitem{MCC}
M.~Charikar, C.~Chekuri, T.~Feder, and R.~Motwani.
\newblock Incremental Clustering.
\newblock Manuscript. December 1996.

\bibitem{crippen}
G.~Crippen and T.~Havel.
\newblock {\em Distance Geometry and Molecular Conformation}.
\newblock Research Studies Press, Letchworth, U.K., 1988.

\bibitem{dill}
K.~Dill.
\newblock Folding proteins: Finding a needle in a haystack.
\newblock {\em Current Opinion in Structural Biology}, 3:99--103, 1993.

\bibitem{DL}
D.~Dobkin and R.~Lipton.
\newblock Multidimensional search problems.
\newblock {\em SIAM J. Computing} 5:181--186 (1976)
\bibitem{finn-et-al}
P.~Finn, D.~Halperin, L.~Kavraki, J.-C. Latombe, R.~Motwani, C.~Shelton, and
  S.~Venkatasubramanian.
\newblock Geometric manipulation of flexible ligands.
\newblock In M.~Lin and D.~Manocha, editors, {\em LNCS Series - 1996 ACM
  Workshop on Applied Computational Geometry}. Springer-Verlag, 1996.

\bibitem{garey}
M.~Garey and D.~Johnson.
\newblock {\em Computers and Intractability: A Guide to the Theory of
  NP-Completeness}.
\newblock Freeman, San Francisco, 1980.

\bibitem{ghose}
A.~Ghose, J.~Kowalczyk, M.~Peterson, and A.~Treasurywala.
\newblock Conformational searching methods for small molecules: I. study of the
  sybyl search method.
\newblock {\em J. of Computational Chemistry}, 14(9):1050--1065, 1993.

\bibitem{glen}
R.~Glen, G.~Martin, A.~Hill, R.~Hyde, P.~Wollard, J.~Salmon, J.~Buckingham, and
  A.~Robertson.
\newblock Computer-aided design and synthesis of 5-substituted tryptamines and
  their pharmacology at the {$5-HT_{10}$} receptor: Discovery of compounds with
  potential anti-migraine properties.
\newblock {\em J. of Medical Chemistry}, 38:3566--3580, 1995.

\bibitem{go70}
N.~Go and H.~Scherga.
\newblock Ring closure and local conformational deformations of chain
  molecules.
\newblock {\em Macromolecules}, 3(2):178--187, 1970.

\bibitem{guida}
W.~Guida, R.~Bohacek, and M.~Erion.
\newblock Probing the conformational space available to inhibitors in the
  thermolysin active site using monter carlo/energy minimization techniques.
\newblock {\em J. of Comp. Chem.}, 13(2):214--228, 1992.

\bibitem{itai}
A.~Itai, N.~Tomioka, M.~Yamada, A.~Inoue, and Y.~Kato.
\newblock Molecular superposition for rational drug design.
\newblock In Kubingi, editor, {\em {3D QSAR} in Drug Design: Theory, Methods,
  and Applications}. ESCOM, 1995.

\bibitem{kavraki-th}
L.~Kavraki.
\newblock {\em Random Networks in Configuration Space for Fast Path Planning}.
\newblock PhD thesis, Stanford University, 1995.

\bibitem{KMY}
S.~Khanna, R.~Motwani, and F.~F. Yao.
\newblock Approximation algorithms for the largest common subtree problem.
\newblock Technical Report STAN-CS-95-1545, Department of Computer Science,
  Stanford University, 1995.

\bibitem{klebe}
G.~Klebe and T.~Mietzener.
\newblock A fast and efficient method to generate biologically relevant
  conformations.
\newblock {\em J. of Computer-Aided Molecular Design}, 8:583--606, 1994.

\bibitem{leach}
A.~Leach.
\newblock A survey of methods for searching the conformational space of small
  and medium sized molecules.
\newblock In K.~Lipkowitz and D.~Boyd, editors, {\em Reviews in Computational
  Chemistry}, volume~2, pages 1--47. VCH Publishers, 1991.

\bibitem{lengauer-gen}
T.~Lengauer.
\newblock Algorithmic research problems in molecular bioinformatics.
\newblock In {\em IEEE Proc. of the 2nd Israeli Symposium on the Theory of
  Computing and Systems}, pages 177--192, 1993.

\bibitem{lipton}
M.~Lipton and W.~Still.
\newblock The multiple minimum problem in molecular modeling: Tree searching
  internal coordinate conformational space.
\newblock {\em J. of Computational Chemistry}, 9(4):343--355, 1988.

\bibitem{lybrand}
T.~Lybrand.
\newblock Computer simulation of biomolecular systems using molecular dynamics
  and free energy perturbation methods.
\newblock In K.~Lipkowitz and D.~B. Boyd, editors, {\em Reviews in
  Computational Chemistry}, volume~1, pages 295--320. VCH Publishers, 1990.

\bibitem{manocha}
D.~Manocha, Y.~Zhu, and W.~Wright.
\newblock Conformational analysis of molecular chains using nano-kinematics.
\newblock {\em Computer Application of Biological Sciences (CABIOS)},
  11(1):71--86, 1995.

\bibitem{martin-pha}
Y.~Martin, M.~Bures, E.~Danaher, J.~DeLazzer, and I.~Lico.
\newblock A fast new approach to pharmacophore mapping and its application to
  dopaminergic and benzodiazepine agonists.
\newblock In {\em J. of Computer-Aided Molecular Design}, volume~7, pages
  83--102, 1993.

\bibitem{martin-db}
Y.~C. Martin, M.~G. Bures, and P.~Willet.
\newblock Searching databases of three-dimensional structures.
\newblock In K.~Lipkowitz and D.~B. Boyd, editors, {\em Reviews in
  Computational Chemistry}, volume~1, pages 213--256. VCH Publishers, 1990.

\bibitem{rand}
R.~Motwani and P.~Raghavan.
\newblock {\em Randomized Algorithms.}
\newblock Cambridge University Press, 1995.

\bibitem{pierre}
D.~Pierre.
\newblock {\em Optimization theory with applications}.
\newblock Dover, NY, 1986.

\bibitem{ragh}
P.~Raghavan and S.~Irani.
\newblock Combinatorial and experimental results for randomized point
matching algorithms.
\newblock {\em Proc. 12th. ACM. SCG}, pp. 68-77, 1996.

\bibitem{pr}
P.~J.~Rezende and D.~T.~Lee.
\newblock Point set pattern matching in $d$-dimensions.
\newblock {\em Algorithmica,} 13, 387-404, 1995.

\bibitem{smellie}
A.~Smellie, S.~Kahn, and S.~Teig.
\newblock Analysis of conformational coverage: 1. validation and estimation of
  coverage.
\newblock {\em J. Chem. Inf. Comput. Sci.}, 35:285--294, 1995.

\bibitem{willet}
P.~Willet.
\newblock Searching for pharmacophoric patterns in databases of
  three-dimensional chemical structures.
\newblock {\em J. of Molecular Recognition}, 8:290--303, 1995.

\end{thebibliography}
}


\end{document}






