\documentclass{article}
\usepackage[affil-it]{authblk}
\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
\usepackage{url}
\usepackage{hyperref}
\hypersetup{colorlinks=false,pdfborder={0 0 0}}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[english]{babel}

\journal{Journal of Molecular Graphics and Modelling}



\begin{document}

\title{Fragment Oriented Molecular Shapes}


\author{Ethan Hain}
\affil{Affiliation not available}

  
\author{Awaiting Activation}
\affil{Affiliation not available}
  
\author{David Koes}
\affil{Affiliation not available}
  



\date{\today}

\maketitle




% introduction and background
\hyphenation{mol-e-cules}

Molecular shape is a fundamental concept in medicinal chemistry \cite{Nicholls2010}, and shape-based virtual screens have successfully identified novel inhibitors \cite{RushIII2005,McMasters2009,Muchmore2006,Naylor_2009}.
Shape is used both to assess the similarity of a candidate molecule to a set of known actives and to evaluate the complementarity of a molecular shape to the shape of the binding site on the target receptor.  Here we describe a novel method for \textit{exactly} specifying shape constraints for querying a large library of molecular shapes.  The shape constraints are derived from both the receptor, which describes where a molecule cannot be, and from the shape of a known binder, which describes where the molecule should be.
Our method is unique in relying on a manually positioned ligand \textit{anchor fragment}. This anchor fragment requirement makes our approach particularly applicable to a fragment-based drug discovery workflow \cite{Rees2004,Congreve2008} as shape constraints are a natural way to search for compounds that extend an identified fragment structure while remaining complementary to the receptor.
Additionally, the anchor fragment, by defining a fixed coordinate system, enables the indexing of large libraries of molecular shapes.  This indexing allows search times to scale sub-linearly with the size of the library, resulting in search performance that is on an interactive time scale. 

Shape-based virtual screening typically attempts to identify the most similar molecules in a virtual library to known active molecules or to a pseudo-ligand that is derived from the desired binding site \cite{Ebalunode2008}.  Shape similarity is usually assessed either through alignment methods, which seek to maximize the three dimensional overlap of two shapes, or through feature vector methods, which transform shapes into a low-dimension vector of features that can be efficiently compared.  As part of the similarity calculation, molecular shapes may be further annotated with electrostatic or pharmacophore features \cite{Vainio2009,Cheeseright2006,Thorner1996,Tervo2005,Marin2008,Sastry2011}.

Alignment methods try to find the optimal overlay of two molecules to either maximize the overlapping volume or the correspondence between feature points, such as molecular field extrema \cite{Vainio2009,Cheeseright2006}.  The predominant method of maximizing volume overlap is to represent the molecular shapes as a collection of Gaussians \cite{Good1993,Grant1996}, sample several starting points, and use numerical optimization to find a local maximum \cite{Hawkins_2007}. An alternative method is to represent the molecular shape by discrete features and use point correspondence algorithms to generate the alignment.  Potential features include pharmacophore features \cite{Sastry2011}, field points \cite{Thorner1996,Vainio2009,Cheeseright2006}, or hyperbolical paraboloid representations of patches of molecular surface \cite{Proschak2008}. A number of performance optimizations to the alignment process have been described \cite{Grant1996,RushIII2005,Sastry2011,Fontaine2007}, but alignment methods remain computationally expensive relative to feature vector methods, unless shapes can be pre-aligned to a canonical coordinate system, as is done with Volumetric Aligned Molecular Shapes (VAMS) \cite{VAMS}.

Feature vector methods reduce molecular shapes to a simple vector of Boolean or numerical features.  Shape similarity is then determined by comparing these vectors using a metric such as Tanimoto or Euclidean distance.  The numerical features can be computed using geometric moments \cite{Ballester2007,Schreyer_2012}, ray-tracing histograms \cite{Zauhar2003}, or a small set of reference shapes \cite{Haigh2005,Putta2002}.
Feature vectors enable computationally efficient screening (millions of shape comparisons per a second) \cite{Ballester2007}, but lack the accuracy and interpret-ability of alignment methods \cite{Nicholls2010}.  Critically, a feature vector similarity does not generate a molecular overlay suitable for visual inspection and analysis.

In our approach, fragment oriented molecular shapes (FOMS), we eliminate the computational burden of alignment by requiring the presence of a common anchor fragment.  Molecules are trivially aligned by a direct superposition of anchor fragments, and the fragment defines a standard coordinate system for describing the shape of the molecule.  Prepositioned molecular fragments are a common component of \textit{de novo} drug design \cite{Schneider2005} where ligands are `grown' from a prepositioned fragment to fit the binding site.  Prepositioned fragments have also been successfully used in structure-based design to identify high-affinity inhibitors for a multitude of targets \cite{Kick1997,Murray1997,Li1998,Liebeschuetz2002}. Notably, our AnchorQuery web service \cite{Koes_2012} provides interactive pharmacophore virtual screening of billions of custom compounds for protein-protein interaction (PPI) inhibitors by pre-aligning compounds to amino acid side-chain motifs corresponding to anchor residues \cite{Rajamani_2004, Meireles_2010} in the PPI interface.
Fragment-based drug discovery workflows \cite{Rees2004,Congreve2008} can provide a physical basis for the selection and positioning of an appropriate anchor fragment.  Alternatively, virtual docking methods may be used \cite{Brenke2009}.

Anchor fragments present a different modality for shape-based screening: the user is required to identify a fragment structure with a meaningful binding mode and the search space is limited to compounds that contain the specified fragment. These requirements enable a new type of search language that supports explicit shape constraints. In essence, a \textit{partial similarity search} \cite{Bronstein2009} can be performed, where instead of optimizing similarity with the entirety of a query shape, the shape constraints specify only part of the shape in detail (e.g., within the binding site) while leaving other parts unspecified (e.g., interactions with solvent).  In addition, the use of anchor fragments enable a new mechanism of search.  Instead of evaluating the query against every molecule in the virtual screening library, the molecular shapes of the library can be \textit{indexed} so that searches need only evaluate a fraction of the library.  This allows large libraries of millions of shapes to be searched on an interactive time scale of a few seconds.


Here we describe the retrospective virtual screening performance of FOMS and explore the potential for our explicit shape constraints to define highly specific filters for the creation of significantly enriched subsets of large virtual libraries.


\section*{Methods}

We describe our representation of molecular shapes, how we define and construct shape constraints for searching, and briefly discuss how these shapes can be indexed to support efficient searches. We also describe our benchmarking methodology for performing a retrospective evaluation of FOMS.\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.42\columnwidth]{figures/rholig/rholig}
\caption{{%startfig 0.8
%nocap%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.42\columnwidth]{figures/voxellig/voxellig}
\caption{{%endfig 0.8
%bigcap
\label{voxelshape}
 The Rho-Kinase inhibitor from PDB 2H9V and  the voxelization of its molecular, solvent-excluded volume.
Voxlized image generated with Sproxel \protect\cite{sproxel}.%
}}
\end{center}
\end{figure}


\subsection*{Shape Representation}

A molecular shape in FOMS is a discretized solvent-excluded volume that is calculated from the heavy atoms of a molecular conformation using a water probe with radius 1.4{\AA}.  The volume is discretized into 0.5{\AA$^3$} voxels (three dimensional pixels) and stored as an oct-tree, an efficient data structure for representing volumetric data \cite{Zhang2009}.  Most oct-tree operations take time proportional to the surface area of a shape ($\approx n^2$) instead of the volume ($\approx n^3$). An example of a voxelized shape is shown in Figure~\ref{voxelshape}.
The grid coordinates of a molecular shape are all computed relative to the identified anchor fragment which is aligned to the axes at the origin.  If a molecule contains multiple instances of an anchor fragment or if there are inherent symmetries in the fragment, then a distinct voxelized shape is generated for every valid alignment of the fragment(s).
For example, a fragment consisting of ring with a single exit vector will map to four distinct alignments on a ligand that contains a ring with two exit vectors.
We grow/shrink molecular shapes by adding/removing the appropriate number of surface voxels.\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/rhoall/rhoall}
\caption{{%nocap
%startfig 0.48%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/rhorecligi/rhorecligi}
\caption{{%endfig 0.48
%bigcap
\label{iptsshape} (Left) The Rho-Kinase inhibitor (sticks) and receptor (surface) from PDB 2H9V shown with the identified interaction points (spheres).  (Right) Shape constraints derived from this structure.  The interaction points become a minimum shape constraint (yellow) by applying a radius of 1{\AA}.  Here the receptor volume has been shrunk by 1{\AA} (two layers of surface voxels were removed) to form a maximum shape constraint (inverse shown in cyan).
Voxlized image generated with Sproxel \protect\cite{sproxel}.%
}}
\end{center}
\end{figure}

\subsection*{ Shape Constraints}

Since we assume all shapes are registered to a common coordinate system defined by the anchor fragment, it is possible to exactly specify regions of space within this coordinate system that a molecule should and should not occupy.  Following the nomenclature of VAMS \cite{VAMS}, we refer to these constraints as minimum and maximum shape constraints. These constraints combine to form an expressive and exacting specification for a desired target shape.
A \textit{minimum shape constraint} sets a strict lower bound on the volumetric shape of a target shape.  Every voxel within the minimum shape constraint must be contained within the target shape.  A minimum shape constraint can be used to require that a target shape has a specific binding mode and minimum bulkiness.
A \textit{maximum shape constraint} sets a strict upper bound on the volumetric shape of a target shape.  Every voxel of the target shape must be contained within the maximum shape constraint.  The maximum shape can be used to constrain the total volume of the target shape and prevent the target shape from overlapping undesirable areas, such as space filled by a receptor. 
Shape constraints are distinct from shape similarity. Unlike shape similarity, which produces a continuous ranking of similarity with respect to a query shape, shape constraints are binary filters: a shape either matches the constraints or does not.

The structure of a receptor-ligand complex provides a natural starting point for the development of shape constraints.  The inverse of the receptor shape provides a starting point for defining a maximum shape constraint.  The constraint can be expanded within the binding site as needed to incorporate tolerance for small steric clashes and receptor flexibility while the shape can be shrunk to avoid excessive exposure of the target ligand to solvent. 
For visualization purposes we display the \textit{inverse} of the maximum shape, which corresponds to the region of space the desired shape is excluded from, to more clearly visualize the minimum shape, which is, by definition, fully subsumed in the maximum shape.

Similarily, the ligand shape can be used to develop a minimum shape constraint.  However, constructing a minimum shape constraint directly from the ligand shape, even if the shape is shrunk or skeletonized, results in a highly specific shape query that may be limited in its ability to retrieve novel chemical scaffolds.
As an alternative to directly using the shape of a bound ligand, we derive \textit{interaction points} from the ligand-receptor complex. Interaction points are points at the center of clusters of ligand atoms that are interacting with the receptor.  Interacting ligand atoms are defined to be those no more than 6{\AA} away from a receptor atom, as measured between atom centers.  These atoms are then clustered into groups of three or more that are no more than 4{\AA} across.  The center of the cluster is an interaction point.  An example of interaction points identified from a ligand-receptor complex is shown with the corresponding minimum and maximum shape constraints in Figure~\ref{iptsshape}.  Interaction points provide a more general specification of the binding mode of a ligand that is less dependent on the ligand chemical scaffold than the full ligand shape and ignores non-interacting components of the ligand.

\subsection*{Shape Indexing}

As with VAMS\cite{VAMS}, we use a matching and packing\cite{Koes_2014} bulk-loading algorithm to initialize an efficient data structure for volumetric shape constraint searches.
  Briefly, shapes are stored in a GSS-tree\cite{keim1999} where each leaf of the tree is a single molecular shape and each internal node includes a maximum included volume (MIV) and minimum surrounding volume (MSV).  The MSV is the union of all the molecular shapes beneath the node in the tree while the MIV is the intersection.  By appropriately applying minimum and maximum shape constraints to the MIV and MSV, it can be determined if any of the shapes lower in the tree have the potential to match the constraints.  If not, the entire subtree of molecular shapes can be eliminated from consideration, resulting in a sub-linear running time.

\subsection*{Shape Similarity}
We use the shape Tanimoto \cite{RushIII2005} to compute the similarity of two voxelized shapes A and B:
$$\delta(A,B) = \frac{A \cap B}{A \cup B}$$
where a larger score indicates a greater degree of similarity. For shape similarity evaluations we consider only similarity with the query ligand.





\subsection*{Dataset}

The specific modality of fragment-oriented molecular shapes requires the creation of a custom benchmark for assessing the virtual screening performance of the method.  In order to construct the shape constraints, a receptor-ligand structure is required, and, in order to screen a library, all compounds in the library must contain the desired anchor fragment.  

We use the Maximum Unbiased Validation (MUV) benchmark \cite{Rohrer2009} as a starting point.  MUV includes sets of 30 active and 15000 property-matched decoy compounds for each of 17 targets.  Compounds are selected from PubChem bioactivity data using a methodology that both reduces the similarity of actives (to avoid analogue bias \cite{Good2008}) and increases the similarity between actives and decoys (which helps prevent artificial enrichment \cite{Verdonk2004}). MUV is also noteworthy in that the decoys were all assessed to be inactive in the initial high-throughput screen against the target providing some measure of experimental evidence that the negatives are true negatives (as opposed to schemes that generate decoys through random sampling \cite{Mysinger_2012}).  The construction of the MUV dataset makes it particularly challenging for ligand similarity approaches. In an evaluation of several screening protocols (including molecular shape), poor enrichments were obtained for all protocols for the majority of query molecules \cite{Tiikkainen_2009}. The use of such a challenging benchmark allows us to critically evaluate FOMS and shape constraints in the context of a dataset where traditional approaches often fail.

Of the 17 targets in the MUV dataset, we identified 10 that had a receptor-ligand structure in the Protein Data Bank (PDB) where the ligand had sub-$\mu$M affinity.  The interaction diagrams of these structures and their PDB codes are shown in Figure~\ref{targets}.  For each of these structures we manually identified interacting fragments that could potentially serve as anchor fragments.  For each target we selected relatively generic functional groups (at most 7 atoms) that were sufficiently common among both the actives and decoys to yield meaningful results and that were clearly forming interactions with the receptor. An exit vector atom was included in generic ring fragments to limit the number of innate symmetries in the fragment and reduce the number for aligned poses. Chosen fragments are shown in Table~\ref{fragtable} which includes the number of matching actives and decoys for each benchmark.  A second fragment is chosen for the PKA target to illustrate the effect of fragment specificity.\selectlanguage{english}
\begin{table*} \centering
    \begin{tabular}{l  c c c c}
        Target & SMARTS & Confirmatory Assay ID & Actives & Decoys \\ 
        \hline
        CathG & c1ccccc1[!H] & 832 & 24 & 12847 \\ % $_3$
        ER$\alpha$ agonist & c1ccccc1[!H] & 737 & 29 & 13390 \\ %$_2$
        ER$\alpha$ & c1ccccc1[!H] & 713 & 22 & 13063 \\ %_3
        ER$\beta$ & c1ccccc1[!H] & 733 & 23 & 12839 \\  %_3
        FAK & a1aaaaa1[!H] & 810 & 29 & 13891 \\ % $_3$
        FXIa & a1aaaa1[!H] & 846 & 24 & 7691 \\ %$_4$
        HIVrt & a1aaaa1[!H] & 652 & 18 & 9535 \\ %$_4$
        HSP90 & c1ccccc1[!H] & 712 & 29 & 12915 \\ % $_2$
        PKA & c1[c,n]cccn1 & 466 & 25 & 5582 \\ % $_1$
        PKA$_\mathrm{alt}$ & a1aaaaa1[!H] &  466 & 29 & 13985 \\ % $_5$
        Rho & c1[c,n]cccn1 & 644 & 18 & 4593 \\ % $_1$
    \end{tabular} 
    \caption{{\label{fragtable}
SMARTS expressions for anchor fragments evaluated in each MUV benchmark shown with the number of active and decoy compounds from the benchmark that contain that fragment. Every MUV benchmark starts with 30 active compounds and 15000 decoys. }} 
\end{table*}\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=1.00\columnwidth]{figures/all2d/all}
\caption{{%fullfig
\label{targets} Graphical depictions of the ten receptor-ligand structures used to define shape constraints. The fragment of the ligand used for alignment is highlighted.  Diagrams generated with PoseView \protect\cite{Stierand_2006}.%
}}
\end{center}
\end{figure}

\subsection*{Virtual Screening Evaluation}

In order to investigate the utility of the FOMS approach we consider both the ability of shape constraint filters to generate enriched subsets and the quality of shape similarity rankings generated using fragment aligned molecules.  We compare to VAMS, which aligns all molecules to a canonical reference system based on their moments of inertia, and rdMolAlign, the shape alignment module of RDKit \cite{rdkit} which aligns molecules using Open3DALIGN \cite{Tosco_2011}. Open3DALIGN uses shape and property information to perform an atom-based alignment that attempts to match the largest number of atoms having similar partial charges and MMFF94 atom types. 
The input conformers for RDKit alignment were the same  aligned poses used with VAMS, and the computed similarity score is the Tanimoto coefficient. 

For each target, conformers of the active and decoy compounds were generated using RDKit \cite{rdkit}.  A maximum of 100 conformers with a minimum RMSD difference of 0.7{\AA} and a maximum energy difference of 10 kcal/mol were generated for each compound.  For each fragment considered, the corresponding conformers were extracted into fragment-specific subsets.  The subsets were then preprocessed to create VAMS and FOMS search databases.  The VAMS database stores a single pose for each conformation aligned along its moments of inertia. For FOMS, if a compound contains multiple instances of the anchor fragment or the fragment contains symmetries, multiple poses per a conformation are stored to account for the multiple fragments/symmetries.

For each target, a single reference structure was identified by searching BindingDB \cite{Liu_2007}, PDBbind \cite{Wang_2005}, and Binding MOAD \cite{Hu_2005} for the complex with the best binding affinity.  The structures were visually inspected to identify buried anchor fragments that made key contacts with the protein, as shown in Figure~\ref{targets} and Table~\ref{fragtable}. All shape queries were constructed using the reference complexes.

We evaluated both shape constraints and shape similarity using the single reference complex. For the shape constraint search we evaluated multiple maximum shape constraints by shrinking the receptor shape by 0, 0.5, 1.0, 1.5, and 2.0 {\AA}ngstroms.  For the minimum shape constraints we used interaction points with minimal radius (one voxel). We evaluated every possible query of interaction points.  That is, if $n$ interaction points were identified, we evaluated $2^n$ minimum shapes. Selecting subsets of interaction points reduces the dependency of the query on the full shape of the reference molecule and more closely mimics the use case where the shape constraints are manually constructed by an expert to select only the most important interactions. As each shape constraint query is a filter, it does not return a total ranking of database compounds, so we report the true positive and false positive rate of the selected subset. We calculate a p-value for the returned subset (relative to the null hypothesis of random selection) using a hypergeometric test.  Reported p-values are corrected for multiple comparisons using the Bonferroni correction (the calculated p-value is multiplied by the number of comparisons).

Shape similarity results are reported using the area under the curve (AUC) of the receiver operating characteristic (ROC) curve. This metric is not biased by the number of compounds and has well-developed statistical properties \cite{Jain_2007}. The AUC is equivalent to the probability that a randomly selected active and randomly selected inactive compound will be correctly ranked by a method.  An AUC of 0.5 corresponds to random chance while a perfect predictor exhibits an AUC of 1.0. 


\section*{Results}

Consistent with previous studies \cite{Tiikkainen_2009}, we find the MUV dataset to be a challenging target for shape-based screening with few targets demonstrating AUCs far from random performance.  Overall FOMS either matched or exceeded the virtual screening performance of VAMS and RDKit while retaining most of the benefits of pre-aligned molecules. Specifically, it is orders of magnitude faster than the optimizing alignment of RDKit suggesting, at a minimum, that FOMS is a viable method for rapidly pre-screening large libraries. In general, the Pareto frontier of the shape constraint queries delineates virtual screen performance equivalent or better to a full shape similarity search while performing queries orders of magnitude faster.

For each target, we plot the ROC curves for FOMS, VAMS, and RDKit similarity search as well as the results for every interaction point shape constraint query. Shape constraint results along the Pareto frontier are highlighted as solid circles and the most statistically significant result is annotated with its Bonferroni-corrected $p$-value. We also plot the total time required for the FOMS, VAMS, and RDKit similarity search and provide a box plot of the distribution of times for the shape constraint searches. All time plots are plotted on a logarithmic scale due to the orders of magnitude difference between the methods. Shape constraint search time varies considerably depending on the number of hits retrieved.  Most queries are highly selective, return few or no results, and take less than a hundredth of a second. We do not include the time needed to generate the shape constraints in the query time since the receptor-based constraints are generated once and cached.  We plot times for only those queries that generate results (queries that were too selective to return any matches were nearly instantaneous).  For these queries the median time, as shown in the box plots, remains below a tenth of a second.  A few queries are non-selective and return subsets of compounds that are comparable in size to the full database. These queries approach or exceed the running time of a full similarity search. The statistically most significant shape constraint query, annotated with a p-value in the ROC plot, has its running time plotted as a circle.  These informative queries typically exceed the median time, but still take less than a second.\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/cathg.f3/cathg.f3}
\caption{{%startfig 0.48
%nocap%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/cathg.f3.time/cathg.f3.time}
\caption{{%endfig 0.48
%bigcap
\label{cathg} Cathepsin G virtual screening performance and search time for various methods.%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/ip/ip}
\caption{{\label{cathgip} Cathepsin G interaction points that defined the shape constraint query with the most significant enrichment. The interaction points (magenta spheres) and fragment (magenta sticks) delineate only half the molecule and are located near a charge-charge interaction and aromatic interaction.  Unlike molecular similarity, shape constraints are capable of selecting molecules that sterically match these features while ignoring the remainder of the query ligand shape and avoiding significant clashes with the receptor shape (for this query the receptor was reduced by 1.5{\AA}).%
}}
\end{center}
\end{figure}

\subsection*{Cathepsin G}
 Virtual screening results for cathepsin G are shown in Figure~\ref{cathg}. All three shape similarity methods exhibited poor performance with this query with AUCs near 0.5, but many of the interaction point shape constraint queries were significantly better than random.  The statistically most significant interaction point query is shown in Figure~\ref{cathgip}.  The fragment for the query is deeply buried within the S1 pocket and the interaction points match the S2 pocket and phosphonate moiety of the ligand. However, there are no interaction points selected in the S3/S4 pocket.  Interactions in this pocket can increase affinity, but are not necessary for binding \cite{de_Garavilla_2005}.  Unlike whole-ligand similarity methods, shape constraint search can ignore a significant part of the geometry of the ligand (the naphthalene in the S3 pocket) and select for only specific parts of the query shape.  In this case, this selectivity results in substantially better virtual screening results that is two orders of magnitude faster than they pre-aligned similarity search methods.\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/eralpha-pot.f2/eralpha-pot.f2}
\caption{{%startfig 0.48
%nocap
%eralpha pot%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/eralpha-pot.f2.time/eralpha-pot.f2.time}
\caption{{%endfig 0.48
%bigcap
\label{eralphapot} ER$\alpha$ agonists virtual screening performance and search time for various methods.%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/eralpha.f3/eralpha.f3}
\caption{{%startfig 0.48
%nocap%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/eralpha.f3.time/eralpha.f3.time}
\caption{{%endfig 0.48
%bigcap
\label{eralpha}  ER$\alpha$ inhibitors virtual screening performance and search time for various methods.%
}}
\end{center}
\end{figure}

\subsection*{Estrogen Receptor $\alpha$}

Virtual screening results for estrogen receptor $\alpha$ agonists and inhibitors are shown in Figures~\ref{eralphapot} and \ref{eralpha}. RDKit was marginally successful at identifying inhibitors with an AUC of 0.59 but the remaining methods were less successful with AUCs ranging from 0.43 to 0.52 and none of the shape-based methods exhibited early enrichment. The shape constraint queries lacked statistical significance and the most significant queries identified a large percentage of the total compounds in the database, resulting in running times that sometimes exceeded a second.\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/erbeta.f3/erbeta.f3}
\caption{{%startfig 0.48
%nocap
%erbeta%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/erbeta.f3.time/erbeta.f3.time}
\caption{{%endfig 0.48
%bigcap
\label{erbeta} ER$\beta$  virtual screening performance and search times for various methods.%
}}
\end{center}
\end{figure}

\subsection*{Estrogen Receptor $\beta$}

Virtual screening results for estrogen receptor $\beta$ are shown in Figure~\ref{erbeta}. None of the methods outperform random selection and all methods fail to exhibit early enrichment.  Although shape constraints are capable of improving on the performance of the similarity methods, they still lack significance.\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/fak.f1/fak.f3}
\caption{{%startfig 0.48
%nocap
%fak1%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/fak.f2.time/fak.f3.time}
\caption{{%endfig 0.48
%bigcap
\label{fak} FAK virtual screening performance and search times for various methods.%
}}
\end{center}
\end{figure}

\subsection*{Focal Adhesion Kinase}

Virtual screening results for focal adhesion kinase (FAK) are shown in Figure~\ref{fak}. Similar to ER$\beta$, none of the methods outperform random selection.  However, there is at least one shape constraint query capable of generating an enriched subset with statistical significance. Note that unlike the other kinases, Rho and PKA, for FAK a nitrogen containing ring pattern was not chosen for fragment matching since fewer than half of the actives in the benchmark had such a pattern.  Unlike the other kinases, using such a pattern for FAK does not improve the virtual screening results (data not shown).\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/fxia.f4/fxia.f4}
\caption{{%startfig 0.48
%nocap
%fxia%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/fxia.f4.time/fxia.f4.time}
\caption{{%endfig 0.48
%bigcap
\label{fxia} Factor XIa virtual screening performance and search times for various methods.%
}}
\end{center}
\end{figure}

\subsection*{Factor XIa}

Virtual screening results for factor XIa (FXIA) are shown in Figure~\ref{fxia}. All three shape similarity methods perform worse than random with AUCs ranging from 0.31 (VAMS) to 0.42 (RDKit and FOMS), with both VAMS and FOMS failing to identify any actives early in the screen.  Shape constraints do not select enriched subsets and the pre-alignment methods exhibit no early enrichment behavior.\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/hivrt.f4/hivrt.f4}
\caption{{%startfig 0.48
%nocap
%hivrt%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/hivrt.f4.time/hivrt.f4.time}
\caption{{%endfig 0.48
%bigcap
\label{hivrt} HIV reverse transcriptase  virtual screening performance and search times for various methods.%
}}
\end{center}
\end{figure}

\subsection*{HIV Reverse Transcriptase}
Virtual screening results for HIV reverse transcriptase (HIV-RT) are shown in Figure~\ref{hivrt}. Again, none of the methods outperform random selection. However, FOMS displays the best early enrichment up to a true positive rate of 17\% with RDKit also achieving some early enrichment. The shape constraint queries are statistically significant and can identify enriched subsets with greater proficiency than the shape similarity methods.

%include ip figure?\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/hsp90.f2/hsp90.f2}
\caption{{%startfig 0.48
%nocap
%hsp90%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/hsp90.f2.time/hsp90.f2.time}
\caption{{%endfig 0.48
%bigcap
\label{hsp90} Heat shock protein 90 virtual screening performance and search times for various methods.%
}}
\end{center}
\end{figure}

\subsection*{Heat Shock Protein 90}

Virtual screening results for heat shock protein 90 (HSP90) are shown in Figure~\ref{hsp90}. All three shape similarity methods perform poorly with AUCs from 0.38 (VAMS) to 0.49 (RDKit). The shape constraint queries can produce statistically significant enrichments and the Pareto frontier matches or exceeds the ROC curve of RDKit. 


 %   The best SMARTS expression and approach for Hsp90 performed significantly better than random performance for FOMS, but not for VAMS (cite table). Hsp 90 inhibitors interact with the ATP binding site. The most important interactions are with Asp93, Gly97, Lys58 , and Asn51 residues, with the study reporting the crystal structure highlighting Asp93 as the most important (cite HSP). Only the Asp93 interaction is covered by our SMARTS expressions and fragment pre-alignments, while the other interactions are also key for Hsp90 inhibition, with the Gly backbone being second in importance. As with other targets in the MUV dataset, the poor performance of FOMS relative to random performance can be partially attributed to the large number of interactions that are not covered by our fragment selection due to the nature of the ligand binding mode.
 %   Docking with 2VCI in smina for the active compounds from the MUV dataset were performed. This revealed 4 potential hydrophobic pockets, with only one covered by our fragment pre-alignment. There were also many active conformers hydrogen-bonding with Asp93 and at least one other residue as well. These docking results indicate the importance of inclusion of every hydrogen-bonding interaction and more than one hydrophobic pocket for Hsp90 inhibition. Our SMARTS expression and fragment pre-alignment only covered one of these interactions and one of the pockets, which most likely led to FOMS not having high performance. 
%More specific SMARTS expressions led to higher performance for FOMS and VAMS (cite table). This makes sense as the more specific fragment choice would help both methods select compounds that can make the hydrogen-bond with Asp93.\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/pka.f/pka.f}
\caption{{%startfig 0.48
%nocap
%pka%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/pka.f.time/pka.f.time}
\caption{{\textit{}%endfig 0.48
%bigcap
\label{pka} Protein kinase A virtual screening performance and search times for various methods.%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/pka.f5/pka.f5}
\caption{{%startfig 0.48
%nocap
%pka5%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/pka.f5.time/pka.f5.time}
\caption{{%endfig 0.48
%bigcap
\label{pka5} Protein kinase A virtual screening performance and search times for various methods with a less selective fragment choice.%
}}
\end{center}
\end{figure}

\subsection*{Protein Kinase A}

Virtual screening results for protein kinase A (PKA) are shown in Figure~\ref{pka}. For this kinase FOMS substantially outperforms both RDKit (AUC 0.64) and VAMS (AUC 0.55) with an AUC of 0.90.  Shape constraints can perform even better than the full similarity search.  Much of this performance is due to the use of a more selective fragment that properly orients the nitrogen in the ring fragment.  Figure~\ref{pka5} shows the result when a more generic fragment, an arbitrary aromatic ring (PKA$_\mathrm{alt}$ in Table~\ref{fragtable}), is used. The enhanced performance of FOMS disappears, and the ability of shape constraints to generated enriched subsets is reduced. In comparison, the optimizing shape similarity of RDKit performs about the same (AUC 0.62). The use of this fragment pulls in thousands more decoys and introduces many more valid fragment alignments for both actives and decoys.  This implies that the good performance of FOMS and shape constraints with this target is likely due the ability of the fragment to properly position compounds to make conserved interactions with the hinge region of the kinase, as shown in Figure~\ref{targets}.  Compounds that contain the specified fragment but cannot position the fragment to make these interactions while maintaining the steric compatibility of the whole ligand with the query ligand and receptor are filtered out.

%    The best SMARTS expression and approach for PKA performed significantly better than random performance for both FOMS and VAMS (cite stats table). The FOMS approach was significantly better than the VAMS approach for PKA (cite stats table).
%Like many other targets in this dataset, PKA inhibitors interact with the ATP binding site. These inhibitors are relatively smaller in size compared to other inhibitors and have few electrostatic interactions. In the example of 1Q8U and its inhibitor, there is only one hydrogen bond with Val123, where the rest of the ligand makes multiple Van der Waals contacts with other parts of the target (cite PKA). SMARTS expressions and fragment pre-alignment covers the only electrostatic interaction and most of the hydrophobic interactions involved in PKA inhibitors, which most likely led to the high performance of FOMS (cite table).
%More specific SMARTS expressions led to higher performance for FOMS (cite table). This makes sense as the more specific fragment choice would cover the only necessary hydrogen-bond interaction with the receptor. Less specific SMARTS expressions led to higher performance for the VAMS method (cite table). (explanation for why this is).\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/rho.f/rho.f}
\caption{{%startfig 0.48
%nocap
%rho%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/rho.f.time/rho.f.time}
\caption{{%endfig 0.48
%bigcap
\label{rho} Rho-kinase 2 virtual screening performance and search times for various methods.%
}}
\end{center}
\end{figure}

\subsection*{Rho-Kinase2}
Virtual screening results Rho-Kinase2 are shown in Figure~\ref{rho} and are similar to those of PKA. FOMS achieves an AUC of 0.94, VAMS an AUC of 0.56, and RDKit an AUC of 0.66. The best shape constraints match or exceed the performance of FOMS. As with PKA, if a less selective fragment is used (a generic aromatic ring), virtual performance is substantially reduced (data not shown).



%    The best SMARTS expression and approach for Rho-Kinase2 performed significantly better than random performance for both FOMS and VAMS (cite stats table). The FOMS approach was significantly better than the VAMS approach for Rho-Kinase2 (cite stats table). 
%Similar to PKA, Rho-Kinase2 inhibitors are relatively small compared to other inhibtors and have few electrostatic interactions. In the case of 2H9V and its inhibitor, there are three hydrogen-bonding interactions with conserved Met, Asn, and Asp. There are also important hydrophobic interactions with other residues inside the pocket (cite Rho). Our SMARTS expression and fragment pre-alignment covers the hydrogen-bonding interaction with the Met and important hydrophobic interactions with conserved Val and Leu residues. This most likely led to the high performance of FOMS for Rho-Kinase2 (cite table).

%More specific SMARTS expressions led to lower performance for both the FOMS and VAMS methods. (explanation for why this is ).\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/aucs-1.00/aucs-1.00}
\caption{{\label{aucs} AUCs for shape similarity methods across all targets shown with their 95\% confidence intervals.  Confidence intervals were calculated using 2000 iterations of bootstrapping.%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{table} \centering
    \begin{tabular}{l c c}
        Target & Best $p$-value & Sig. Queries \\ 
        \hline
CathG & 4.23e-06 & 1\% \\
ER$\alpha$ & 0.0505 & 0\% \\
ER$\alpha$ agonist & 0.0106 & 0\% \\
ER$\beta$ & 0.0328 & 0\% \\
FAK & 0.000289 & 4\% \\
FXIa & 0.387 & 0\% \\
HIVrt & 0.000951 & 5\% \\
HSP90 & 0.000755 & 2\% \\
PKA & 1.14e-23 & 80\% \\
Rho & 8.41e-15 & 80\% \\
    \end{tabular} 
    \caption{{\label{pvaltable} Significance of the best shape constraint queries and the percentage of all evaluated queries with a Bonferroni-corrected $p$-value below 0.01.}} 
\end{table}\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/confs/confs}
\caption{{%startfig 0.48
%confs
%nocap%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[h!]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/cdist/cdist}
\caption{{%endfig 0.48
%bigcap
\label{confs}  (Left) The effect of the maximum number of conformers generated for each molecule on screening performance, and (Right) the change in distribution of similarity scores for active and decoy compounds using FOMS for different numbers of conformers.%
}}
\end{center}
\end{figure}

\section*{Discussion}

The MUV dataset, with its focus on eliminating analogue bias, is particularly resistant to single-query shape-based virtual screens \cite{Tiikkainen_2009}.  This is reflected in our overall results, shown in Figure~\ref{aucs}, where only two targets (Rho and PKA) achieve AUCs where the 95\% confidence interval does not overlap with 0.5 (random performance). The remaining targets likely lack meaningful whole-molecule shape complementary between the query ligand and the active compounds of the benchmark.  One exception may be HIV-rt, where there is clear early enrichment which indicates that a subset of the actives may be compatible with the query molecule.
FOMS dramatically outperformed other methods for the Rho and PKA targets due to correct positioning of a fragment with key, conserved interactions.  For comparison, Figure~\ref{aucs} also shows the performance of a 2D fingerprint, the OpenBabel FP2 \cite{O_Boyle_2011} path-based fingerprint.  2D information is more successful for three targets (ER$\beta$, FXIa, HSP90), and has comparable or worse performance for the remaining targets, illustrating the orthogonality of 2D and 3D approaches.

Our assumption in the design of this study was that, since the evaluated methods utilize rigid conformers, it would be necessary to generate a large sampling of conformers to ensure the biologically relevant conformer was screened and was correctly registered as a hit.  Surprisingly, when we tested this assumption with the Rho Kinase target (results were similar for PKA), we found that the virtual screening performance for all three shape methods was relatively insensitive to the number of conformers used (Figure~\ref{confs}).  Interestingly, this is not because the top ranked conformer is the most representative conformer.  Instead, as demonstrated in Figure~\ref{confs}, reducing the number of conformers does result in a reduction in similarity scores for active compounds.  However, a compensating reduction in scores is observed in the decoy set as the number of conformers sampled is decreased, resulting in similar virtual screening results.  This suggests that sampling a large number of conformations will be beneficial in terms of producing meaningful poses with high similarity to the query ligand, a key advantage of shape-based methods, but, counter-intuitively, may not be strictly necessary for virtual screening purposes.

FOMS essentially provides a rapid means of template docking \cite{Ruiz_Carmona_2014,abagyan2015icm,Koes_2012} using shape-based scoring. The disadvantage of fragment-oriented approaches is they are critically dependent on the choice of fragment and its proper positioning in defining the query. Provided these requirements can be met, there are several advantages to shape-based fragment alignment search.  By enforcing the fragment alignment, key interaction are guaranteed to be conserved. Previous studies have demonstrated the importance of adding pharmacophoric properties (or `color') to shape similarity \cite{Hawkins_2007}.  Fragment alignment introduces a hard bias toward matching a key portion of the query molecule without introducing any additional computation, as required by more general methods.  In fact, as we have shown, pre-alignment substantially reduces the computational overhead. 

Prealignment, whether to fragments (FOMS) or canonical internal coordinates (VAMS) is orders of magnitude faster than methods that dynamically optimize the alignment.  This holds true even if the cost to create the search database is taken into account.  The time to create the databases scales with the number of molecular shapes (about 10 shapes a second on our system) and compares favorable with RDKit search (2 molecules a second).
As the common case is for a fragment oriented database to be re-used to for queries by multiple users investigating multiple targets, in practice the cost of database creation gets amortized into insignificance.

A major advantage of fragment alignment is that it enables the use of shape constraints.
Shape constraint search generally tracked or improved upon the performance of FOMS similarity ranking (e.g. Figure~\ref{cathg}).  As shown in Table~\ref{pvaltable}, shape constraints were able to generate statistically significant ($p < 0.01$) enriched subsets for six of the ten targets.
Unlike whole-molecule shape similarity, shape constraints can select for a \textit{subshape} of the query ligand (i.e., perform partial shape matching) and specifically filter out potential clashes with the receptor.  Furthermore, shape constraints support an indexing search algorithm that scales sub-linearly.  In our evaluation, shape constraint searches were orders of magnitude faster than even VAMS similarity search and generally took well under a second. 

Shape constraints provide a novel and unique method for specifying molecular shape queries.  Although we have investigated automatically generated shape constraints, our assumption is that intelligently designed constraints created by human experts will substantially outperform these interaction point based constraints.  The speed of shape constraint searches enables interactive analysis and refinement of shape queries where the expert user sculpts their desired molecule guided by the results of searches of full-sized databases. We anticipate that fragment-oriented shape constraints will prove a useful tool in both fragment-based drug discovery workflows, which are naturally centered around a fragment, and lead optimization workflows, where the core scaffold of the lead compound can serve as a fragment.



\subsection*{Acknowledgements}

This work was supported by the National Institute of Health [R01GM097082, R01GM108340].

\selectlanguage{english}
\FloatBarrier
\bibliographystyle{plain}
\bibliography{bibliography/converted_to_latex.bib%
}

\end{CJK}\end{document}

