Spherical Harmonics Coefficients for Ligand-Based Virtual Screening of Cyclooxygenase Inhibitors

PLOS ONE, Jul 2011

Background Molecular descriptors are essential for many applications in computational chemistry, such as ligand-based similarity searching. Spherical harmonics have previously been suggested as comprehensive descriptors of molecular structure and properties. We investigate a spherical harmonics descriptor for shape-based virtual screening. Methodology/Principal Findings We introduce and validate a partially rotation-invariant three-dimensional molecular shape descriptor based on the norm of spherical harmonics expansion coefficients. Using this molecular representation, we parameterize molecular surfaces, i.e., isosurfaces of spatial molecular property distributions. We validate the shape descriptor in a comprehensive retrospective virtual screening experiment. In a prospective study, we virtually screen a large compound library for cyclooxygenase inhibitors, using a self-organizing map as a pre-filter and the shape descriptor for candidate prioritization. Conclusions/Significance 12 compounds were tested in vitro for direct enzyme inhibition and in a whole blood assay. Active compounds containing a triazole scaffold were identified as direct cyclooxygenase-1 inhibitors. This outcome corroborates the usefulness of spherical harmonics for representation of molecular shape in virtual screening of large compound collections. The combination of pharmacophore and shape-based filtering of screening candidates proved to be a straightforward approach to finding novel bioactive chemotypes with minimal experimental effort.

A PDF file should load here. If you do not see its contents the file may be temporarily unavailable at the journal website or you do not have a PDF plug-in installed and enabled in your browser.

Alternatively, you can download the file locally and open with any standalone PDF reader:


Spherical Harmonics Coefficients for Ligand-Based Virtual Screening of Cyclooxygenase Inhibitors

et al. (2011) Spherical Harmonics Coefficients for Ligand-Based Virtual Screening of Cyclooxygenase Inhibitors. PLoS ONE 6(7): e21554. doi:10.1371/journal.pone.0021554 Spherical Harmonics Coefficients for Ligand-Based Virtual Screening of Cyclooxygenase Inhibitors Quan Wang 0 Kerstin Birod 0 Carlo Angioni 0 Sabine Gro sch 0 Tim Geppert 0 Petra Schneider 0 Matthias 0 Rupp 0 Gisbert Schneider 0 Paul Wrede, Charite-Universitatsmedizin Berlin, Germany 0 1 Frankfurt Institute for Advanced Studies (FIAS), Goethe University , Frankfurt, Germany , 2 Institute for Clinical Pharmacology, Goethe University , Frankfurt, Germany , 3 Institute of Pharmaceutical Sciences, Swiss Federal Institute of Technology (ETH), Zu rich, Switzerland, 4 Machine Learning Group, Technical University , Berlin , Germany Background: Molecular descriptors are essential for many applications in computational chemistry, such as ligand-based similarity searching. Spherical harmonics have previously been suggested as comprehensive descriptors of molecular structure and properties. We investigate a spherical harmonics descriptor for shape-based virtual screening. Methodology/Principal Findings: We introduce and validate a partially rotation-invariant three-dimensional molecular shape descriptor based on the norm of spherical harmonics expansion coefficients. Using this molecular representation, we parameterize molecular surfaces, i.e., isosurfaces of spatial molecular property distributions. We validate the shape descriptor in a comprehensive retrospective virtual screening experiment. In a prospective study, we virtually screen a large compound library for cyclooxygenase inhibitors, using a self-organizing map as a pre-filter and the shape descriptor for candidate prioritization. Conclusions/Significance: 12 compounds were tested in vitro for direct enzyme inhibition and in a whole blood assay. Active compounds containing a triazole scaffold were identified as direct cyclooxygenase-1 inhibitors. This outcome corroborates the usefulness of spherical harmonics for representation of molecular shape in virtual screening of large compound collections. The combination of pharmacophore and shape-based filtering of screening candidates proved to be a straightforward approach to finding novel bioactive chemotypes with minimal experimental effort. - Funding: The authors received an academic MOE software license from Chemical Computing Group Inc., Montreal, Canada. M.R. acknowledges partial support from DFG grant MU 987/4-2 and the FP7-ICT programme of the European Community under the PASCAL2 network of excellence, ICT-216886. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Competing Interests: The authors received an academic MOE software license from Chemical Computing Group Inc., Montreal, Canada. There are no patents, products in development or marketed products to declare. This does not alter the authors adherence to all the PLoS ONE policies on sharing data and materials, as detailed online in the guide for authors. . These authors contributed equally to this work. Current address: Institute for Pure and Applied Mathematics (IPAM), University of California, Los Angeles, Los Angeles, California, United States of America Ligand-based virtual screening [1,2], quantitative structureproperty and structure-activity relationships [3,4], and other concepts in computational medicinal chemistry are based on the similarity principle [5], which states that (structurally) similar compounds generally exhibit similar properties. Such methods require quantitative representations of molecules, usually in the form of chemical descriptors, i. e., computable numerical attributes in vector form [6]. Numerous molecular 3D-descriptors and alignment methods have been proposed. Examples include CoMFA (comparative molecular field analysis) [7], Randic molecular profiles [8], 3DMoRSE code (3D-molecule representation of structures based on electron diffraction) [9], invariant moments and radial scanning and integration [10], radial distribution function descriptors [11], WHIM (weighted holistic invariant molecular descriptors) [12], length-to-breadth ratios [13], USR (ultrafast shape recognition, based on statistical moments) [14], ROCS (rapid overlay of chemical structures, based on Gaussian densities) [15], VolSurf (volumes and surfaces of 3D molecular fields) [16], GETAWAY (geometry, topology, and atom weights assembly) [17], and shrinkwrap surfaces [18], to name just a few prominent representatives. In computer graphics, several methods exist for the more general problem of comparing arbitrary 3D objects [19,20], including distribution-based shape histograms [21], the D2 shape descriptor [22], and, the scaling index method [23]; the viewbased methods of extended Gaussian images [24], and the light field descriptor [25]; the surface decomposition-based methods of Zernike moments [26], REXT (radialized spherical extent function) [27], and spherical harmonics descriptors [28]. Spherical harmonics have been used in cheminformatics as a global feature-based parametrization method of molecular shape [2838]. Their attractive properties with regard to rotations make them an intuitive and convenient choice as basis functions when searching in a rotational space [31]. A review article by Venkatraman et al. [38] highlights applications of spherical harmonics to protein structure comparison, ligand binding site similarity, protein-protein docking, and virtual screening. Jakobi et al. [37] use spherical harmonics in their ParaFrag approach to derive 3D pharmacophores of molecular fragments. Recently, Ritchie and co-workers have applied the ParaSurf and ParaFit methodologies [32,33] (Cepos InSilico Ltd., Erlangen, Germany) in a virtual screening study on the directory of useful decoys (DUD) data set [39], which motivates 3D shape-property combinations specifically for flexible ligands [40]. The DUD data set was also used in a comparative analysis of the performance of various shape descriptors alone and in combination with property and pharmacophore features [41]. See the section on related methods for further discussion of spherical harmonics approaches. In this work, we introduce a partially rotation-invariant descriptor of molecular shape based on spherical harmonics decomposition coefficients. The idea is to decompose the molecular surface using spherical harmonics and to use the norm of the decomposition coefficients as a description of molecular shape. In this, we take advantage of the fact that the norm of the coefficients does not change under rotation around the z-axis, which we align to the primary axis of the molecule. We retrospectively evaluate our descriptor, and prospectively apply it to screen for novel inhibitors of the enzymes cyclooxygenase-1 (COX-1) and cyclooxygenase-2 (COX-2). Particular focus is on the practical application of the virtual screening technique as an evaluation of its actual suitability for early-phase drug discovery. Materials and Methods Spherical harmonics Let l[f0,1,2, . . .g, let jmjl, let (h,w) indicate spherical coordinates, and let Plm denote the Legendre polynomials [42]. The spherical harmonics [43] Ylm of order l (frequency, angular quantum number) and degree m (azimuthal quantum number), form an orthonormal (with respect to integration over the unit sphere) and complete set of basis functions (Fig. 1). They are solutions to Laplaces differential equation +2Y~0 in spherical coordinates [44]. Any square-integrable spherical function f (h,w) can be decomposed as frequency detail. The coefficients are unique, and can therefore be used as feature vectors for shape description. Rotation of a molecule (its shape function f ) changes the coefficients. A conventional solution is to define a canonical orientation of the molecule. For the purpose of shape comparisons, this implies an alignment of the compared molecules, with all associated problems and computational requirements. As an alternative, we use a partial orientation in conjunction with certain rotational invariance properties of the coefficients. Descriptor definition Let X~(~x1, . . . ,~xn)T[Rn|3 denote the Cartesian coordinates of points ~x1, . . . ,~xn[R3 sampled from a molecular surface. We assume that the surface is star-like (single-valued) in the sense that rays radiating outward from the molecules origin intersect the surface only once (this is more of an issue for proteins; as argued elsewhere [35], small molecules are little, if at all, affected). Let B[Rn|k denote the spherical harmonics basis functions Ylm evaluated at ~x1, . . . ,~xn, 0lL, {lml, with L[N the maximum order used and k~ PlL~0 Plm~{l 1~(Lz1)2 the number of basis functions. The sampled molecular surface X can be reconstructed using a matrix C[Rk|3 of coefficients as X~BC. The coefficient matrix is given by C~B{X, where B{ denotes the pseudo-inverse [46] of B. Lemma. The p-norm of the rows of C does not change under rotation around the z-axis (polar axis, change in w). Proof It is sufficient to consider a single coefficient, i.e., n~k~1, and ~c~b{1~x~(x,y,z)=b. Here, ~x~(x,y,z)[R3 is a sampled surface point, ~c[R3 are coefficients, and b[R is the spherical harmonics basis function Ylm. From Eq. 1, it is clear that if w changes, only the part exp (imw) of the spherical harmonic basis function Ylm(h,w) changes, while the rest of Ylm stays constant. Thus, b~ exp (imw)d for some constant d[R depending only on m,l,h, but not on w. Since j exp (imw)j~1, and thus jbj~jdj, with complex coefficients clm. The spherical harmonics decomposition can be viewed as a generalization of the Fourier decomposition to three dimensions [45]. The coefficients clm of an harmonic expansion can be found using the orthonormality property. Multiplying each side of Eq. 2 by the complex conjugate Y lm(h,w) and integrating over the sphere yields f (h,w)Y lm(h,w) sin (h)dwdh: Small values of l correspond to low frequencies, and describe the overall low-resolution shape; higher values of l add finer, highz p 1=p jj~c0jjp~jdj{1 jj~x0jjp. Since the rotation jj~xjjp~jj~x 0jjp, and it follows that jj~cjjp~jj~c 0jjp. Our spherical harmonics descriptor is the k-component vector (jj~c 0jjp,jj~c1{1jjp,jj~c1jjp,jj~c 1jjp,jj~c2{2jjp, . . . ,jj~c LLjjp) 0 0 1 of the norms of the coefficients. It is a description of molecular shape that is invariant to rotations around the z-axis. Before spherical harmonics decomposition, we place molecules into a common frame of reference by translating their center of gravity to the coordinate system origin and by aligning their first principle component (the direction of maximum variance as given by principle component analysis [47]) with the z-axis. In other words, we align molecules according to their longest spatial extent, and then apply our descriptor which is invariant to rotations around the z-axis. Descriptor computation Gaussian contact surfaces [48] of all compounds were computed using MOE (Molecular Operating Environment, version 2009.10, Chemical Computing Group Inc., Montreal, Canada, www. chemcomp.com ). Spherical harmonics decomposition was then carried out on the vertices of these surfaces, giving approximate coefficients [49]. To limit computational expense, we truncated spherical harmonics expansions after order L~9. The resulting k~(9z1)2~100 decomposition coefficients were sufficient to represent fine molecular detail and approximately reconstruct the original molecular surfaces (Fig. 2). The partial rotational invariance of the coefficient norms jjclmjj is demonstrated in Fig. 3. Computation was done in Matlab (The MathWorks, version R2007a, www.mathworks.com ), partly based on code by Dr. Andrew Hanna (University of East Anglia, United Kingdom, www.cmp.uea.ac.uk/,aih ). Average computing time was v3 seconds per compound, which is acceptable for medium-sized libraries but will require speed-up for high-throughput virtual screening. Related methods Spherical harmonics have been widely used in cheminformatics as a global feature-based parametrization method of molecular shape [2838]. Most current approaches, including ours, use the center of gravity as the center of the spherical harmonics decomposition. Molecular surface sampling can be done by sampling iso-probability surfaces of molecular property densities. One aspect in which methods differ is the way they deal with rotations in 3D space. Ritchie and Kemp [31] apply the rotational property of spherical harmonics (a rotation of the surface can be simulated by rotating the expansion coefficients) to maximize the pairwise superposition of two molecules. The software ParaSurf superposes molecules using a brute-force rotational search over the three Euler rotation angles [50]. In a recent publication, Cai et al. [36] use a similar approach to obtain the minimal root-mean-square distance between a ligand molecule and a target protein. In these related studies, molecular surfaces were rotated by transforming their expansion coefficients. Standard orientation of compounds prior to spherical harmonics decomposition was proposed by Morris et al. [34]. Their work registered molecules and binding pockets in a standard frame by translating their center of mass to the coordinate origin and aligning their variance-covariance matrix to the axes of the coordinate system. They then use the coefficients of a real spherical harmonics expansion to describe and compare the molecular shape of binding pockets and ligands. This approach aligns molecules to minimize rotation-dependent differences in the coefficients. Rotation-invariant spherical harmonics descriptors were applied by Kazhdan et al. [28] and Mavridis et al. [35,51], using the fact that expansion coefficients of the same order l transform among themselves to construct rotationally invariant spherical harmonics coefficients jj Plm~{l clmjj2. In their approach, coefficients of the same order l are binned together, thereby losing information contained in the individual degrees m, but gaining complete rotational invariance. In this work, we combine partial orientation of the molecules with the magnitude of the expansion coefficients as a partially rotation-invariant shape descriptor. Our proposed descriptor retains more information than the spherical harmonics descriptors by Kazhdan et al. [28] and Mavridis et al. [35,51] in the sense that coefficients within the same order are not summed up, but kept. Compared with standard orientation methods, our descriptor is potentially less susceptible to problems in the orientation step than most others because only the first (and most stable) principle component is used for orientation. Retrospective evaluation For retrospective validation, we ranked the compounds in a database according to their similarity to a reference compound, as measured by Euclidean distance and our descriptor. Two conceptually different collections of reference data were used, the DUD data set (release 2, from http://dud.docking.org/r2 , unmodified data) [39], and the COBRA data set (version 10.3, 11 244 compounds annotated with activity on a total of 677 individual macromolecular targets) [52]. COBRA 10.3 contains 168 COX-2 inhibitors. Gaussian contact surfaces were generated with the MOE 2009.10 (Molecular Operating Environment, Chemical Computing Group Inc., Montreal, Canada, www.chemcomp.com ) GaussianSurface function, with parameter pos set to aPos a, rad dock_aRadius a, nearpos aPos a, neardist 5, maxMb 1, and fuzzy 0. All other parameters were kept at their default values. Virtual screening experiments in COBRA were carried out using a single conformation generated by CORINA (version 2007, Molecular Networks GmbH, Erlangen, Germany, www.molecular-networks.com ). We used the selective COX-2 inhibitor SC-558 and the nonselective inhibitor indomethacin as queries for ligand-based similarity searching, with the conformations extracted from the crystal structure (protein data bank [53] identifiers (PDB ID) 6cox [54] and 4cox [54]). Enrichment factors [55], receiver operating characteristic curves (ROC curves [56]), and the area under these curves (ROC AUC) were used as performance measures. Prospective virtual screening We screened the ChemBridge compound pool (457 226 compounds, ChemBridge Corp., San Diego, USA, www.chembridge.com) for potential COX ligands using a single CORINA conformer query as in the retrospective screening. The database was preprocessed using the washing procedure in MOE (protonation of strong bases and de-protonation of strong acids; all other parameters were kept at their default values). To reduce computational effort and allow for pharmacophore feature-based compound ranking, the screening compound pool was pre-filtered using a self-organizing map (SOM [57]) trained on the ChemBridge collection and 275 COX-1 and COX-2 inhibitors from the COBRA database. SOM topology was toroidal with 20|20 neurons (1 200 molecules per neuron on average); compounds were represented using the 150-dimensional CATS2D topological pharmacophore descriptor [58,59] and compared using the Manhattan distance. The initial width of the Gaussian neighborhood function was set to 5; training was terminated after 5:106 steps (using each compound 10 times on average). We used the MOLMAP software tool for SOM generation [60]. After pre-filtering, 21 950 compounds of the ChemBridge database that were similar to the COX inhibitors from the COBRA database were retained for virtual screening using our spherical harmonics shape descriptor. Two potent COX inhibitors served as reference molecules (queries; Fig. 4). All parameters were set to the values used in retrospective virtual screening. The spherical harmonics descriptor was calculated for the 21 950 retained molecules and for the two reference molecules. Enzyme inhibition assay Inhibition of COX-1 (ovine) and COX-2 (human recombinant) activity was measured using a COX inhibitor screening assay kit (Cayman Chemicals, Ann Arbor, MI, USA, www.caymanchem. com ), according to the manufacturers protocol. SC-560, a selective COX-1 inhibitor, and celecoxib, a selective COX-2 inhibitor, served as positive controls. The COX inhibitor screening assay directly measures the amount of prostaglandins PGE2, PGD2 and PGF2a produced by SnCl2 reduction of COXderived PGH2. In addition to this protocol, the amounts of prostaglandins were quantified by LC-MS/MS analysis as described previously [61]. Whole blood assay COX-1 whole blood assay. One-milliliter heparinized human blood samples were incubated with 4 ml test substance (in DMSO) or 4 ml DMSO (control) for 10 min at 37 0 C . After this, thrombocyte aggregation was stimulated by addition of calcium ionophore A23187 (50 mM) for 30 min at 37 0C . Plasma was separated by centrifugation for 20 min at 2000 rpm, 4 0 C and kept at {80 0 C until assayed for TXB2 by LC-MSMS (see below). COX-2 whole blood assay. For the determination of COX2 activity, 1 ml of heparinized human blood was incubated at 37 0 C with 10 ml of acetylsalicylic acid (1 mg ml{1 in PBS), 4 ml DMSO or 4 ml inhibitor (in DMSO) for 10 min. After this, 2 ml of LPS (5 mg ml{1 in DMSO) was added and incubated for 24 hr at 37 0 C . The reaction was terminated by quickly chilling on ice. Plasma was separated by centrifuging (20 min, 2 000 rpm, 4 0 C ), stored at {80 0 C until analysis of prostaglandins by LC-MS/MS within two weeks. Extraction procedure. 250 ml plasma was incubated with 600 ml 45 mM H3PO4, 100 ml 0:15 M EDTA, 10 ml BHT (butylated hydroxytoluene, 2 mg ml{1), 20 ml MeOH, 20 ml internal standard PGE2 {d4 (25 ng ml{1), PGD2 {d4 (25 ng ml{1), PGF 2a {d4 (10 ng ml{1), 6k PGF {1 a {d4 (10 ng ml{1), TXB 2 {d4 (25 ng ml{1) for 1 min, and passed through a ABS ELUT-Nexus cartridge (Varian, Darmstadt, Germany) preconditioned with methanol (1 ml), followed by distilled water (1 ml). The cartridge was washed with distilled water (1 ml) and 30 % MeOH (1 ml). PGE2, PGF2a, 6{keto{PGF la and TXB2 were eluted with hexaneethylacetate-isopropranolol (30:65:5, v/v, 1 ml). After vaporating the solvent under nitrogen atmosphere, the residue was reconstituted in 50 ml acetonitrile / formic acid. PGE2 concentrations were quantified by means of a validated LC-MS/ MS assay described previously [61]. The lower limit of quantification was 10 pg ml{1. Results and Discussion We validated our spherical harmonics (SpH) descriptor in a retrospective setting (statistical validation on known data), and in a prospective study to obtain biochemical confirmation of our model. Retrospective evaluation As a first analysis, we used the DUD compound collection for a preliminary comparison of selected shape- and structure-based virtual screening methods. ROC AUC [62] values were computed average + std.dev. aValues reproduced from the original study by Vainio et al. [68], with ShaEP and ROCS run in shape-only mode. DUD COX-2 data: 426 actives, 13 289 decoys. Query SC558 has PDB ID 6cox; indomethacin has PDB ID 4cox. n.d. = not determined. doi:10.1371/journal.pone.0021554.t001 for each of the methods compared. ROC AUC values lie in the interval 0,1 , with values closer to 1 indicating higher enrichment of actives in a ranked list of compounds. The analysis was limited to the original COX-2 data from DUD (426 actives, 13 289 decoys). We did not perform exhaustive comparative analyses of virtual screening performance or focus on early recognition of actives [63,64], as the primary purpose of this study was to determine whether our SpH descriptor might be a useful shapebased filtering criterion for COX inhibitors. Retrospective screening was restricted to COX-2, our original target. Table 1 summarizes the results obtained for CATS2D (topological pharmacophore descriptor [58]), LIQUID (threedimensional pharmacophore descriptor using Gaussian feature points; v1: hydrogen-bond donors, hydrogen-bond acceptors, lipophilic [65]; v2: additional aromatic, positive and negative charge features (manuscript in preparation)), PRPS (Gaussian pseudoreceptor model [66,67]), ShaEP (field-based subgraph matching [68]), and ROCS (Gaussian shape model [15,69]). For the DUD COX-2 data, ROC AUC values indicate better than random performance for all methods. SpH yielded an average of 0.86, which compares to Ritchies ParaFit spherical harmonics descriptor (note that the ParaFit ROC AUC value is not given in the original publication; we estimated it from graphical material provided in the articles supplementary material [41]). Among the tested methods, SpH performed best for the selective COX-2 inhibitor SC-558 (Fig. 4) yielding a ROC AUC = 0.91. Notably, high values were also obtained for indomethacin (Fig. 4), a nonselective COX inhibitor (COX-1 IC 50 = 18 nM; COX-2 IC 50 = 26 nM) [70]. Apparently, only the PRPS pseudoreceptor model distinguished between the selective (ROC AUC = 0.83) and the non-selective (ROC AUC = 0.15) query. In contrast to DUD (unmodified data), the COBRA database contains only druglike bioactive compounds. Ranking of the COBRA database with SC-558 as query resulted in an enrichment factor (computed for the first percentile) of 23. We compared this result to those obtained by the shapelets [71] method from our group, using the same version of the COBRA database and the same reference structure. The shapelets shape-only virtual screening method achieved a comparable enrichment factor of 24. ROC curves are presented in Fig. 5 (numbers for shapelets refer to COBRA version 8.4 containing 8 311 compounds including 136 COX-2 inhibitors). In summary, our spherical harmonics coefficients-based approach SpH achieves notable enrichment of actives and seems suitable for COX-2 inhibitor retrieval. This outcome is in agreement with the study of shape-based virtual screening approaches by Ritchie et al. [41], who report high hit rates for COX-2 using shape descriptors. We conclude that spherical harmonics-based decomposition of molecular shape captures structural features that are relevant for virtual screening. Due to the limited number of published prospective applications [72], it seems premature to render any conclusion regarding certain implementation preferences or best-in-class spherical harmonics methods. To further assess our SpH approach, we performed a prospective study using SpH in a virtual screening cascade with the aim to identify new COX inhibitors. Prospective virtual screening Virtual screening. We used a SOM to pre-select potential COX inhibitors from the screening compound pool. The SOM (Fig. 6) of COX activity islands contains six neurons with more than three ligands (neurons (1,16), (1,14), (1,15), (7,18), (18,14), (10,14) with 49, 25, 15, 14, 12, 11 ligands, respectively). We selected all compounds from the ChemBridge database contained in these neurons, 21 950 in total (5 % of the pool). In the second virtual screening step, SpH was used for shapebased filtering. Two reference molecules (SC-558 and indomethacin; Fig. 4) resulted in two ranked lists of the pre-filtered ChemBridge compounds. 10 duplicates were found among the 50 top-ranking compounds from the two lists (20% overlap). In total, 12 compounds were selected by visual inspection, preferring potentially new scaffolds (cherry-picking, Fig. 7), and submitted for activity determination in a direct enzyme inhibition and a whole blood assay. Assay results We determined the COX-inhibitory activity of 12 compounds by performing a commercially available competitive COXinhibition assay using purified COX-1 (ovine) and COX-2 (human recombinant) enzymes. Compounds 5 and 9 inhibit COX-1 in a concentration dependent manner (Fig. 8 and Table 2). At 100 mM compounds 5 and 9 inhibit COX-1 activity to 72 % and 89 %, respectively. Both compounds have only marginal effects on COX-2-activity at concentrations up to 100 mM. All other substances have no effect on COX-1 or COX-2 activity in this in vitro assay. While this outcome supports our general virtual screening approach, we failed to retrieve COX-2 inhibitors. This might be a consequence of using the selective COX-2 inhibitor SC-558 in combination with the non-selective COX inhibitor indomethacin as queries for the spherical harmonics shape filter. Apparently, the COX activity island on the SOM and SpH consensus filtering eliminated COX-2 specific features. It is also possible that there were no hitherto unidentified COX-2 ligands in the compound pool. In the whole blood assay (Fig. 9, Table 3), compounds 5 and 9 are less effective, with maximum COX-1 inhibition of about 30 % and no COX-2 inhibitory efficacy. Interestingly, in this assay, compounds 6, 10, 2 and 8 inhibit TXB 2 production in a concentration dependent manner up to 72 %, 72 % and 61 % at 100 mM, respectively. Compounds 6 and 10 have only marginal inhibitory potency on PGE 2 production, which points to selective COX-1 inhibitors in vivo. Compound 2 also inhibits PGE 2 production comparable to TXB 2, indicating that this compound is a COX-unselective inhibitor. In contrast, substance 8 increases the PGE 2 amount in a concentration dependent manner, which argues for an activator of PGE 2 production in the cellular context. All other compounds show only very weak or no effect on PGE 2 production. The inhibitory data obtained from the whole blood assay might be meaningful for further hit optimization. Compounds that are active in this assay are not snatched away by binding to serum albumin, but cross the cell membrane and overcome possible interactions with cellular substances or enzymes. This could explain why compounds 5 and 9 are active in the enzyme assay, but inactive in the whole blood assay. In contrast, compounds 6, 10, 2 and 8, which were more active in the whole blood assay, possibly interact with the arachidonic acid pathway in other ways than direct inhibition of COX-1 or COX-2. Also, these compounds might be metabolized by cellular enzymes to more active derivatives, but this hypothesis needs to be tested by further experiments. Compound 8 is of special interest, as it induces PGE 2 production up to 322 %. This increase could be due to an activation of enzyme activity, possibly by binding to the inactive monomer of the COX-homodimer complex [73,74], or, due to an enhancement of COX-2 protein, either by transcriptional or posttranscriptional mechanisms. As a preliminary novelty check, similarity searches were performed using SciFinder Web (2010-10-21) for data retrieval from the CAS database (Chemical Abstracts Service, Columbus, Ohio, USA; www.cas.org). For none of the actives any reference to COX inhibition was found, and only for compound 9 substructure Figure 8. COX-2 inhibition in vitro assay results. Shown are COX-1 (blue) and COX-2 (red) inhibition. Celecoxib and SC-560 are known inhibitors selective for COX-2 and COX-1, respectively. doi:10.1371/journal.pone.0021554.g008 1 2 ? 3 4 5 6? 7 8 ? 9 10? 11 12 2.1 9.7 0.7 3.2 8.1 68:3{ 9.4 2.6 0.7 5.5 2.9 0.1 Primed (0) compounds are shown in Fig. 8, starred (?) compounds are shown in Fig. 9. Discussed compounds are shown in bold face. {For these three measurements, high standard deviations were observed. This could be due to solubility problems, impurities, protein degradation, or other unspecific effects. The corresponding measurements should be treated carefully. doi:10.1371/journal.pone.0021554.t002 matches (lacking the meta methyl group) were retrieved with regard to bioactivities other than COX inhibition. It is therefore reasonable to conclude that COX inhibition by compounds 5 and 9 represents a novel finding resulting from our study. We did not perform additional analytical investigations of compound integrity and purity other than those provided by the compound supplier. Therefore, we cannot exclude that the activities measured in the assays might be partially owed to decomposition or oxidation products. Analog compound design and testing will be mandatory. TXB2 is indicative of COX-1 activity, PGE2 of COX-2 activity relative to DMSO control. Primed (0) compounds are shown in Fig. 8, starred (?) compounds are shown in Fig. 9. Discussed compounds are shown in bold face. doi:10.1371/journal.pone.0021554.t003 Conclusions We presented a favorable retrospective evaluation of the SpH approach using COX-2 data from the DUD collection, and in a first prospective application demonstrated the usefulness of the descriptor in combination with a self-organizing map for retrieving bioactive ligands from a large compound pool. Although we did not retrieve a potent COX-2 inhibitor, which is likely owed to the setup of the virtual screening cascade, two novel COX-1 inhibitors were discovered. Future research will have to focus on mathematical descriptions of molecular shape that allow for enzyme subtype-selective ligand screening. We introduced the magnitude of spherical harmonics coefficients as a partially rotation-invariant descriptor of molecular shape. In retrospective validation on the DUD dataset, the performance (as estimated by ROC AUC) of our shape-only method was comparable to other shape-based similarity searching methods. Results show that the magnitude of spherical harmonics decomposition coefficients can be used to describe molecular shape in a partially rotation-invariant way, resulting in a notable enrichment of active compounds in virtual and real screening studies. The combination of pharmacophore filtering by a selforganizing map and shape-filtering by spherical harmonics descriptors might be a useful two-step virtual screening protocol for hit retrieval from large screening compound collections. Conceived and designed the experiments: QW SG MR GS. Performed the experiments: QW KB CA TG MR GS. Analyzed the data: QW SG PS MR GS. Contributed reagents/materials/analysis tools: PS SG. Wrote the paper: QW MR GS. 1. Bohm HJ , Schneider G, eds. ( 2000 ) Virtual Screening for Bioactive Molecules . Weinheim, Germany: Wiley-VCH. 2. Douguet D ( 2008 ) Ligand-based approaches in virtual screening . Curr Comput Aided Drug Des 4 : 180 - 190 . 3. Jurs P ( 2003 ) Quantitative structure-property relationships . In: Gasteiger J, ed. Handbook of Chemoinformatics: From Data to Knowledge, Wiley, volume 3 , chapter 1 .2. pp 1314 - 1135 . 4. Kubinyi H ( 2003 ) QSAR in drug design . In: Gasteiger J, ed. Handbook of Chemoinformatics: From Data to Knowledge, Wiley, volume 4 , chapter 4 .2. pp 1532 - 1554 . 5. Johnson M , Maggiora G, eds. Concepts and Applications of Molecular Similarity . New York : Wiley. 6. Rupp M , Schneider P , Schneider G ( 2009 ) Distance phenomena in highdimensional chemical descriptor spaces: Consequences for similarity-based approaches . J Comput Chem 30 : 2285 - 2296 . 7. Cramer III R , Patterson D , Bunce J ( 1988 ) Comparative molecular field analysis (CoMFA). 1. Effect of shape on binding of steroids to carrier proteins . J Am Chem Soc 110 : 5959 - 5967 . 8. Randic M , Kleiner A , Alba LD ( 1994 ) Distance/distance matrices . J Chem Inf Comput Sci 34 : 277 - 286 . 9. Schuur J , Gasteiger J ( 1997 ) Infrared spectra simulation of substituted benzene derivatives on the basis of a 3D structure representation . Anal Chem 69 : 2398 - 2405 . 10. Robinson D , Barlow T , Richards G ( 1997 ) The utilization of reduced dimensional representations of molecular structure for rapid molecular similarity calculations . J Chem Inf Comput Sci 37 : 943 - 950 . 11. Hemmer M , Steinhauer V , Gasteiger J ( 1999 ) Deriving the 3D structure of organic molecules from their infrared spectra . Vib Spectros 19 : 151 - 164 . 12. Gramatica P , Corradi M , Consonni V ( 2000 ) Modelling and prediction of soil sorption coefficients of non-ionic organic pesticides by molecular descriptors . Chemosphere 41 : 763 - 777 . 13. Todeschini R , Consonni V ( 2000 ) Handbook of Molecular Descriptors . Weinheim, Germany: Wiley- VCH, first edition. 14. Ballester PJ , Richards WG ( 2007 ) Ultrafast shape recognition to search compound databases for similar molecular shapes . J Comput Chem 28 : 1711 - 1723 . 15. Grant A , Gallardo A , Pickup B ( 1996 ) A fast method of molecular shape comparison: A simple application of a Gaussian description of molecular shape . J Comput Chem 17 : 1653 - 1666 . 16. Cruciani G , Crivori P , Carrupt PA , Testa B ( 2000 ) Molecular fields in quantitative structurepermeation relationships: the VolSurf approach . J Mol Struct 503 : 17 - 30 . 17. Consonni V , Todeschini R , Pavan M , Gramatica P ( 2002 ) Structure / response correlations and similarity / diversity analysis by GETAWAY descriptors. 2. Application of the novel 3D molecular descriptors to QSAR/QSPR studies . J Chem Inf Comput Sci 42 : 693 - 705 . 18. Van Drie JH ( 1997 ) ' 'shrink-wrap'' surfaces: A new method for incorporating shape into pharmacophoric 3D database searching . J Chem Inf Comput Sci 37 : 38 - 42 . 19. Shilane P , Min P , Kazhdan M , Funkhouser T ( 2004 ) The Princeton shape benchmark . In: Proceedings of the International Conference on Shape Modeling & Applications (SMI 2004 ), Genova, Italy, June 7- 9 IEEE Computer Society. pp 167 - 178 . 20. Iyer N , Jayanti S , Lou K , Kalyanaraman Y , Ramani K ( 2005 ) Threedimensional shape searching: State-of-the-art review and future trends . Comput Aided Des 37 : 509 - 530 . 21. Ankerst M , Kastenmuller G , Kriegel HP , Seidl T ( 1999 ) 3D shape histograms for similarity search and classification in spatial databases . In: Gu ting R , Papadias D , Lochovsky F, eds. Proceedings of the 6th International Symposium on Spatial Databases (SSD 1999 ), Hong Kong , China, July 20 -23 Springer. pp 207 - 228 . 22. Osada R , Funkhouser T , Chazelle B , Dobkin D ( 2001 ) Matching 3D models with shape distributions . In: Proceedings of the 3rd International Conference on Shape Modeling & Applications (SMI 2001 ), Genova, Italy, May 7- 11 IEEE Computer Society. pp 154 - 166 . 23. Jamitzky F , Stark R , Bunk W , Thalhammer S , Rath C , et al. ( 2001 ) Scalingindex method as an image processing tool in scanning-probe microscopy . Ultramicroscopy 86 : 241 - 246 . 24. Horn B ( 1984 ) Extended Gaussian images . Proc IEEE 72 : 1671 - 1686 . 25. Chen DY , Tian XP , Shen YT , Ouhyoung M ( 2003 ) On visual similarity based 3D model retrieval . Comput Graph Forum 22 : 223 - 232 . 26. Novotni M , Klein R ( 2003 ) 3D Zernike descriptors for content based shape retrieval . In: Proceedings of the 8th ACM Symposium on Solid modeling and Applications (SM 2003 ), Seattle, Washington, USA, June 16 -20 Association for Computing Machinery. pp 216 - 225 . 27. Vranic D ( 2003 ) An improvement of rotation invariant 3D-shape based on functions on concentric spheres . In: Proceedings of the International Conference on Image Processing (ICIP 2003 ), Barcelona, Spain, September 14-17 IEEE Computer Society , volume 3 . pp 757 - 760 . 28. Kazhdan M , Funkhouser T , Rusinkiewicz S ( 2003 ) Rotation invariant spherical harmonic representation of 3D shape descriptors . In: Kobbelt L, Schroder P , Hoppe H, eds. Proceedings of the 1st Eurographics Symposium on Geometry Processing (SGP 2003 ), Aachen, Germany, June 23 -25 Association for Computing Machinery. pp 167 - 175 . 29. Max N , Getzoff E ( 1988 ) Spherical harmonic molecular surfaces . IEEE Comput Graph Appl 8 : 42 - 50 . 30. Duncan B , Olson A ( 1993 ) Approximation and characterization of molecular surfaces . Biopolymers 33 : 219 - 229 . 31. Ritchie D , Kemp G ( 1999 ) Fast computation, rotation, and comparison of low resolution spherical harmonic molecular surfaces . J Comput Chem 20 : 383 - 395 . 32. Ritchie DW , Kemp GJL ( 2000 ) Protein docking using spherical polar fourier correlations . Proteins: Struct Funct Bioinf 39 : 178 - 194 . 33. Lin JH , Clark T ( 2005 ) An analytical, variable resolution, complete description of static molecules and their intermolecular binding properties . J Chem Inf Model 45 : 1010 - 1016 . 34. Morris R , Najmanovich R , Kahraman A , Thornton J ( 2005 ) Real spherical harmonic expansion coefficients as 3D shape descriptors for protein binding pocket and ligand comparisons . Bioinformatics 21 : 2347 - 2355 . 35. Mavridis L , Hudson B , Ritchie D ( 2007 ) Toward high throughput 3D virtual screening using spherical harmonic surface representations . J Chem Inf Model 47 : 1787 - 1796 . 36. Cai W , Xu J , Shao X , Leroux V , Beautrait A , et al. ( 2008 ) SHEF: A vHTS geometrical filter using coefficients of spherical harmonic molecular surfaces . J Mol Model 14 : 393 - 401 . 37. Jakobi AJ , Mauser H , Clark T ( 2008 ) Parafrag-an approach for surface-based similarity comparison of molecular fragments . J Mol Model 14 : 547 - 558 . 38. Venkatraman V , Sael L , Kihara D ( 2009 ) Potential for protein surface shape analysis using spherical harmonics and 3D Zernike descriptors . Cell Biochem Biophys 54 : 23 - 32 . 39. Huang N , Shoichet B , Irwin J ( 2006 ) Benchmarking sets for molecular docking . J Med Chem 49 : 6789 - 6801 . 40. Perez-Nueno VI , Venkatraman V , Mavridis L , Clark T , Ritchie DW ( 2011 ) Using spherical harmonic surface property representations for ligand-based virtual screening . Mol Inf 30 : 151 - 159 . 41. Venkatraman V , Perez-Nueno VI , Mavridis L , Ritchie DW ( 2010 ) Comprehensive comparison of ligand-based virtual screening tools against the DUD data set reveals limitations of current 3D methods . J Chem Inf Model 50 : 2079 - 2093 . 42. Abramowitz M , Stegun I ( 1972 ) Handbook of Mathematical Functions . New York : Dover. 43. Press W , Teukolsky S , Vetterling W , Flannery B ( 2007 ) Numerical Recipes . The Art of Scientific Computing. Cambridge : Cambridge University Press , third edition. 44. Vilenkin NY ( 1968 ) Special Functions and the Theory of Group Representations , volume 22 of Translations of Mathematical Monographs. Washington DC: American Mathematical Society. 45. Funkhouser T , Min P , Kazhdan M , Chen J , Halderman A , et al. ( 2003 ) A search engine for 3D models . ACM Trans Graph 22 : 83 - 105 . 46. Ben-Israel A , Greville T ( 2003 ) Generalized Inverses . Theory and Applications Springer, second edition. 47. Jolliffe I ( 2004 ) Principle Component Analysis . New York : Springer, second edition. 48. Grant A , Pickup B , Nicholls A ( 2001 ) A smooth permittivity function for PoissonBoltzmann solvation methods . J Comput Chem 22 : 608 - 640 . 49. Brechbu hler C , Gerig G , Ku bler O ( 1995 ) Parametrization of closed surfaces for 3-D shape description . Comput Vis Image Understand 61 : 154 - 170 . 50. Clark T ( 2010 ) ParaSurf 10 User Manual CePos InSilico Ltd ., The Old Vicarage, 132 Bedford Road, Kempston, United Kingdom . 51. Mavridis L , Ritchie DW ( 2010 ) 3D-blast: 3D protein structure alignment, comparison, and classification using spherical polar fourier correlations . In: Proceedings of the 15th Pacific Symposium on Biocomputing (PSB 2010 ), Maui, Hawaii, USA, January 3- 7 . pp 281 - 292 . 52. Schneider P , Schneider G ( 2003 ) Collection of bioactive reference compounds for focused library design . QSAR Comb Sci 22 : 713 - 718 . 53. Berman H , Westbrook J , Feng Z , Gilliland G , Bhat T , et al. ( 2000 ) The protein data bank . Nucleic Acids Res 28 : 235 - 242 . 54. Kurumbail R , Stevens A , Gierse J , McDonald J , Stegeman R , et al. ( 1996 ) Structural basis for selective inhibition of cyclooxygenase-2 by anti-inammatory agents . Nature 384 : 644 - 648 . 55. Hawkins P , Warren G , Skillman G , Nicholls A ( 2008 ) How to do an evaluation: pitfalls and traps . J Comput Aided Mol Des 22 : 179 - 190 . 56. Jain A , Nicholls A ( 2008 ) Recommendations for evaluation of computational methods . J Comput Aided Mol Des 22 : 133 - 139 . 57. Kohonen T ( 2001 ) Self-Organizing Maps . New York : Springer, third edition. 58. Schneider G , Neidhart W , Giller T , Schmid G ( 1999 ) ''Scaffold-hopping'' by topological pharmacophore search: A contribution to virtual screening . Angew Chem Int Ed 38 : 2894 - 2896 . 59. Fechner U , Schneider G ( 2004 ) Optimization of a pharmacophore-based correlation vector descriptor for similarity searching . QSAR Comb Sci 23 : 19 - 22 . 60. Schneider G , Wrede P ( 1998 ) Artificial neural networks for computer-based molecular design . Progr Biophys Mol Biol 70 : 175 - 222 . 61. Schmidt R , Coste O , Geisslinger G ( 2005 ) LC-MS/MS-analysis of prostaglandin E2 and D2 in microdialysis samples of rats . J Chrom B 826 : 188 - 197 . 62. Fawcett T ( 2006 ) An introduction to ROC analysis . Pattern Recogn Lett 27 : 861 - 874 . 63. Mackey MD , Melville JL ( 2009 ) Better than random? The chemotype enrichment problem . J Chem Inf Model 49 : 1154 - 1162 . 64. Truchon JF , Bayly C ( 2007 ) Evaluating virtual screening methods: Good and bad metrics for the ''early recognition'' problem . J Chem Inf Model 47 : 488 - 508 . 65. Tanrikulu Y , Nietert M , Scheffer U , Proschak E , Grabowski K , et al. ( 2007 ) Scaffold hopping by''fuzzy'' pharmacophores and its application to RNA targets . Chem Bio Chem 8 : 1932 - 1936 . 66. Tanrikulu Y , Schneider G ( 2008 ) Pseudoreceptor models in drug design: Bridging ligand- and receptor-based virtual screening . Nat Rev Drug Discov 7 : 667 - 677 . 67. Tanrikulu Y , Proschak E , Werner T , Geppert T , Todoroff N , et al. ( 2009 ) Homology model adjustment and ligand screening with a pseudoreceptor of the human histamine H4 receptor . Chem Med-Chem 4 : 820 - 827 . 68. Vainio MJ , Puranen JS , Johnson MS ( 2009 ) ShaEP: Molecular overlay based on shape and electrostatic potential . J Chem Inf Model 49 : 492 - 502 . 69. Rush III TS , Grant JA , Mosyak L , Nicholls A ( 2005 ) A shape-based 3-d scaffold hopping method and its application to a bacterial protein-protein interaction . J Med Chem 48 : 1489 - 1495 . 70. Riendeau D , Percival MD , Boyce S , Brideau C , Charleson S , et al. ( 1997 ) Biochemical and pharmacological profile of a tetrasubstituted furanone as a highly selective COX-2 inhibitor . Br J Pharmacol 121 : 105 - 117 . 71. Proschak E , Rupp M , Derksen S , Schneider G ( 2008 ) Shapelets: Possibilities and limitations of shape-based virtual screening . J Comput Chem 29 : 108 - 114 . 72. Ripphausen P , Nisius B , Bajorath J ( 2011 ) State-of-the-art in ligand-based virtual screening . Drug Discov Today 16 : 372 - 376 . 73. Yuan C , Rieke CJ , Rimon G , Wingerd BA , Smith WL ( 2006 ) Partnering between monomers of cyclooxygenase-2 homodimers . Proc Natl Acad Sci USA 103 : 6142 - 6147 . 74. Vecchio AJ , Simmons DM , Malkowski MG ( 2010 ) Structural basis of fatty acid substrate binding to cyclooxygenase -2. J Biol Chem 285 : 22152 - 22163 .

This is a preview of a remote PDF: http://journals.plos.org/plosone/article/file?id=10.1371/journal.pone.0021554&type=printable

Quan Wang, Kerstin Birod, Carlo Angioni, Sabine Grösch, Tim Geppert, Petra Schneider, Matthias Rupp, Gisbert Schneider. Spherical Harmonics Coefficients for Ligand-Based Virtual Screening of Cyclooxygenase Inhibitors, PLOS ONE, 2011, DOI: 10.1371/journal.pone.0021554