Structure-based virtual screening for drug discovery: principles, applications and recent advances.

Structure-based drug discovery (SBDD) is becoming an essential tool in assisting fast and cost-efficient lead discovery and optimization. The application of rational, structure-based drug design is proven to be more efficient than the traditional way of drug discovery since it aims to understand the molecular basis of a disease and utilizes the knowledge of the three-dimensional structure of the biological target in the process. In this review, we focus on the principles and applications of Virtual Screening (VS) within the context of SBDD and examine different procedures ranging from the initial stages of the process that include receptor and library pre-processing, to docking, scoring and post-processing of topscoring hits. Recent improvements in structure-based virtual screening (SBVS) efficiency through ensemble docking, induced fit and consensus docking are also discussed. The review highlights advances in the field within the framework of several success studies that have led to nM inhibition directly from VS and provides recent trends in library design as well as discusses limitations of the method. Applications of SBVS in the design of substrates for engineered proteins that enable the discovery of new metabolic and signal transduction pathways and the design of inhibitors of multifunctional proteins are also reviewed. Finally, we contribute two promising VS protocols recently developed by us that aim to increase inhibitor selectivity. In the first protocol, we describe the discovery of micromolar inhibitors through SBVS designed to inhibit the mutant H1047R PI3Kα kinase. Second, we discuss a strategy for the identification of selective binders for the RXRα nuclear receptor. In this protocol, a set of target structures is constructed for ensemble docking based on binding site shape characterization and clustering, aiming to enhance the hit rate of selective inhibitors for the desired protein target through the SBVS process.

INTRODUCTION

The identification of lead compounds showing pharmacological activity against a biological target and the progressive optimization of the pharmacological properties and potency of these compounds are the focal points of early-stage drug discovery. To this end, the pharmaceutical industry has adopted the experimental screening of large libraries of chemicals against a therapeutically-relevant target (high-throughput screening or HTS) as a means to identify new lead compounds. Through HTS, active compounds, antibodies or genes, which modulate a particular biomolecular pathway, may be identified; these provide starting points for drug discovery and for understanding the role of a particular biochemical process in biology. Although HTS remains the method of choice for drug discovery in the pharma industry, the various drawbacks of this method, namely the high cost, the time-demanding character of the process as well as the uncertainty of the mechanism of action of the active ingredient have led to the increasing employment of rational, structure-based drug design (SBDD) with the use of computational methods.

SBDD is nowadays central to the efficient development of therapeutic agents and to the understanding of metabolic processes. SBDD is proven to be more efficient than the traditional way of drug discovery since it aims to understand the molecular basis of a disease and utilizes the knowledge of the three-dimensional (3D) structure of the biological target in the process. By using computational methods and the 3D structural information of the protein target, we are now able to investigate the underlying molecular interactions involved in ligand-protein binding and thus interpret experimental results in atomic-level detail. The use of computers in drug discovery bears the additional advantage of delivering new drug candidates more quickly and cost-efficiently.

State of the art structure-based drug design methods include virtual screening (VS) and de novo drug design; these serve as an efficient, alternative approach to HTS. In virtual screening, large libraries of drug-like compounds that are commercially available are computationally screened against targets of known structure, and those that are predicted to bind well are experimentally tested [1, 2]. However, database screening does not provide molecules that are structurally “novel” as these molecules have been previously synthesized by commercial vendors. Existing molecules can only be patented with a “method of use” patent covering their use for a unique application and not their chemical structure. In the de novo drug design approach, the 3D structure of the receptor is used to design structurally novel molecules that have never been synthesized before using ligand-growing programs and the intuition of the medicinal chemist [3].

Computer-aided drug discovery has recently had important successes: new biologically-active compounds have been predicted along with their receptor-bound structures and in several cases the achieved hit rates (ligands discovered per molecules tested) have been significantly greater than with HTS [1, 4-6]. Moreover, while it is rare to deliver lead candidates in the nM regime through VS, several reports in the recent literature describe the identification of nM leads directly from VS; these strategies will be discussed herein [7-9]. Therefore, computational methods play a prominent role in the drug design and discovery process within the context of pharmaceutical research.

In this review, we focus on the principles and applications of VS in the SBDD framework, starting from the initial stages of the process that include receptor and library pre-processing, to docking, scoring, and post-processing of top-scoring hits. We also highlight several successful studies and protocols that led to nM leads, discuss novel applications of Structure-Based VS (SBVS) such as substrate identification for the discovery of novel metabolic pathways, and provide recent trends in library design. Limitations of SBVS are also examined. Finally, we present two developed VS protocols that aim to enhance inhibitor selectivity for the target protein structure.

VIRTUAL SCREENING IN STRUCTURE-BASED DRUG DISCOVERY

The general scheme of a SBVS strategy is shown in Fig. () [1, 2, 5]. SBVS starts with processing the 3D target structural information of interest. The target structure may be derived from experimental data (X-ray, NMR or neutron scattering spectroscopy), homology modeling, or from Molecular Dynamics (MD) simulations. There are numerous fundamental issues that should be examined when considering a biological target for SBVS; for example, the druggability of the receptor, the choice of binding site, the selection of the most relevant protein structure, incorporating receptor flexibility, suitable assignment of protonation states, and consideration of water molecules in a binding site, to name a few. In fact, the identification of ligand binding sites on biological targets is becoming increasingly important. The need for novel modulators of protein/gene function has recently directed the scientific community to pursue druggable allosteric binding pockets. Another consideration for SBVS includes the careful choice of the compound library to be screened in the VS exercise according to the target in question, and the preprocessing of libraries in order to assign the proper stereochemistry, tautomeric, and protonation states.

Following library and receptor preparation, each compound in the library is virtually docked into the target binding site with a docking program. Docking aims to predict the ligand-protein complex structure by exploring the conformational space of the ligands within the binding site of the protein. A scoring function is then utilized to approximate the free energy of binding between the protein and the ligand in each docking pose. Docking and scoring produce ranked compounds, which are then post-processed by examining calculated binding scores, validity of generated pose, undesirable chemical moieties, metabolic liabilities, desired physicochemical properties, lead-likeness, and chemical diversity. Post-processing results in a small number of selected compounds, which proceed to experimental assaying [1, 2, 5, 10].

Protein Preparation Schemes for SBVS

The success of a SBVS campaign largely depends on reasonable starting structures for both the protein and the ligand. A typical PDB structure file consists only of heavy atoms (if the input is an X-ray structure) and may contain water molecules, cofactors, activators, ligands, and metal ions as well as several protein subunits. Moreover, the structure has generally no information on bond orders, topologies, or formal atomic charges. Terminal amide groups and asparagine residues may be misarranged as X-ray structures cannot unambiguously distinguish between O and NH2 groups. Ionization and tautomeric states are also unassigned in most cases, residue side chains or larger loops may be missing because of low resolution of a particular protein area, and steric clashes may exist.

To efficiently address the above-mentioned structural issues, several protein preparation schemes have been proposed [5, 11, 12]. The general proposed strategy is to first determine the protonation states of the aminoacids in the protein using available software. Popular freely available software include PROPKA [13], H++ [14], SPORES [15]. The next step is to assign hydrogen atoms and optimize protein hydrogen bonds according to an optimal hydrogen bond network. A widely-used software for these tasks is the PDB2PQR software [16]. The next steps are assignment of partial charges, capping of residues, treating metals, filling in missing loops and missing side chains, and minimizing the protein structure to relieve steric clashes. In addition, a decision needs to be made regarding whether water molecules will be left in or removed from the binding site. Several methods have been developed to tackle this challenging problem, such as 3D RISM [17-19], SZMAP [20], JAWS [21], WaterMap [18, 22], which are utilized in commercial software [23-25]. In case the protein is co-crystallized with substrates, cofactors, etc., these ligands must be prepared to create 3D geometries, assign proper bond orders, and generate accessible tautomer and ionization states prior to VS. These tasks may be performed within available free software packages such as the “Protein Preparation Wizard” of Maestro [26] or scripts that use a variety of different programs to prepare a protein structure for SBVS, e.g. WebPDB (see Table for a full list of freely-available programs).

The importance of protein preparation in docking performance has been recently reported [11]. Docking performance is significantly improved when employing a best-practice preparation scheme over using minimally-prepared structures from the PDB structure. In the only systematic study available in the literature regarding the influence of the protein preparation procedure in the success of a SBVS workflow, Sastry et al. have explored each of the steps involved in preparing a system for VS [11]. In their study, a large number of parameters in protein preparation is examined using the Glide validation set of 36 crystal structures and 1,000 decoys. Several protocols are considered and applied on the Directory of Useful Decoys (DUD) database [27, 28], showing that database enrichment is improved after proper receptor preparation and that neglecting certain steps of the preparation process produces a systematic degradation in enrichments, which may be large for some targets.

Binding Site Identification

Binding site identification is often an additional prerequisite for performing SBVS, when the binding site is not known or when new, allosteric modulators of protein function are sought. Ideally, the target binding site is a pocket, typically a concave, having a variety of probable hydrogen bond donors and acceptors and hydrophobic characteristics. Currently four approaches in the identification of putative binding sites exist in the literature: 1) Static approaches, where computational solvent mapping with chemical probes (small organic molecules) is utilized to identify binding hot spots on a 3D structure (from Xray, MD, etc). These approaches determine the hot spot druggability and provide information for drug design. Examples of such an approach include SiteMap [29], FTMap [30], Fpocket [31], MDpocket [32], QsiteFinder [33], MED-SUMO [34] and SiteHound-web [35]. Selecting a different chemical probe each time may result in the identification of a different type of binding site. 2) Dynamic approaches, where the probes and the protein evolve dynamically in time and more than one probe may be used in one simulation. Analysis of these simulations provides direct access to interaction free energies between the protein and small organic molecules, which can then be used to detect binding sites and predict the maximal affinity that a drug-like fragment could attain. Programs following this approach are MDMix [36], SILCS [37], MixMD [38]. In this context, studies that aim to understand ligand binding mechanisms using microsecond MD simulations or enhanced sampling techniques [39-41] may also be used for revealing new binding pockets using small organic molecules as probes, although this approach is not computationally-efficient. 3) In the mixed approach, the detection of ligand binding sites is performed with a chemical probe as in the static approach and in addition, the putative binding site is evaluated in terms of flexibility. Flexibility can be determined by performing normal mode analysis or by identifying flexible residues within the binding site and probing residue alternative conformations using a rotamer library. FTFlex, which is an extension of the FTMap server, uses the latter approach [42]. 4) The last approach uses water as a probe to identify putative binding sites on proteins. Approaches following the inhomogeneous solvation theory [43] include JAWS [21], WaterMap [22], and WATMD (Novartis), but also other methods have been developed such as 3D-RISM [17] and SZMAP [20].

Compound Database Preparation

The construction of compound databases is the next important step in the SBVS process. Databases for SBVS contain drug-like small molecules, often freely available or available via purchase or synthesis, which possess desirable characteristics such as stability and solubility in aqueous media, existence of appropriate functional groups to interact with biological targets and absence of toxic and undesirable moieties. Several rules have been applied to ensure ‘drug-likeness’, with the most popular being the “Lipinski Rule of Five" [44], which states that drug-like compounds should have molecular weight lower than 500, lipophilicity (logP) lower than 5, less than five hydrogen bond donors, and less than 10 hydrogen bond acceptors. As increasingly more compounds break some of these rules and enter the market (e.g. natural product drugs as well as 50% of marketed drugs do not comply with the “Rule of Five”) [45], attempts to improve the predictions of druglikeness have spawned many extensions to the Lipinski Rule of Five. Molecules may be

considered as drug-like within a range of log P in −0.4 to +5.6, molar refractivity from 40 to 130, molecular weight from 180 to 500, number of atoms from 20 to 70 (includes H-bond donors (e.g.; OH's and NH's) and H-bond acceptors (e.g.; N's and O's), polar surface area no greater than 140 Å2, and/or fewer than ten rotatable bonds [46]. Moreover, the “Rule of Three” is an extension to the Rule of Five specifically for fragments (molecular weight < 300, logP < 3, number of hydrogen bond donors and acceptors < 3, number of rotatable bonds < 3) [47]. However, in a recent study by Klebe and coworkers, it was shown that only four of their 11 discovered fragments were compliant with the Rule of Three [48]. Strict compliance to this rule would have strongly limited the variety of chemotypes among the fragment hits, indicating that it is advisable to follow these rules with a certain degree of flexibility, i.e. allow one violation to the rule each time [49]. The overly rigid application of strict cut-off points may introduce artificial distinctions between similar compounds and runs the risk of missing valuable opportunities [50]. A more recent rule accounting for the physico-chemical properties relating to preclinical toxicology outcomes is the “Pfizer's Rule of 3/75", which is based on the values of calculated partition coefficient (ClogP) and topological polar surface area (TPSA). According to this rule, compounds with ClogP lower than 3 and TPSA higher than 75 are approximately 2.5 times more likely to be safe in in vivo assays [51]. The “Jorgensen Rule-of-Three” is another widely-followed rule for lead like properties and states that the aqueous solubility measured as logS should be greater than -5.7, the apparent Caco-2 cell permeability should be faster than 22 nm/s, and the number of primary metabolites should be less than 7; these limits are based on the properties of 90% of 1700 oral drugs [52]. Bickerton et al. have proposed a measure of drug-likeness based on the concept of desirability called “the quantitative estimate of drug-likeness” (QED) [53]. The QED functions rank chemical structures by their merit relative to the properties of oral drugs. Moreover, QED functions have the additional advantage that they are based on the underlying distribution data of drug properties and may identify cases in which a generally unfavorable property may be tolerated when the other parameters are close to ideal. Other metrics of promise include ligand efficiency (LE = potency/heavy atom count) and lipophilic LE (LLE = potency - cLogP), which is a measure of lipophilicity that leads to the observed potency. These metrics show how well a molecule engages in receptor interactions; taking these rules into account provides excellent insights for lead optimization [49].

It should be noted that the majority of commercial compounds found in chemical libraries have a larger molecular weight and higher hydrophobicity compared to orally available drugs. This shift to higher hydrophobicity and lower solubility is also reflected in compounds with significant biological activity found in the ChEMBL database, indicating a trend towards larger, more hydrophobic compounds. In particular, drugs that have been recently approved were found to possess much lower median logP values than those previously suggested as optimal [54].

It has been demonstrated that docking huge compound collections such as PubChem (30M compounds) and ChemSpider (26M compounds), in other words ‘blind docking’, not only is time-consuming and computationally demanding, but also results in redundant possibilities and poses a great burden for compound selection. Apart from filtering for lead-like properties, one should exclude known toxicophores or metabolically liable moieties [55]. Moreover, several novel filters have been recently developed, aiming at the quality enhancement of database content. The Pan Assay Interference Compounds (PAINS) [56] and the ALARM-NMR [57] filters, contain compounds found to be chemically reactive and assay-interfering, appear as frequent hitters, and are not identified by toxicophoric filters. Therefore, a combination of filtering for desired pharmacological and adsorption, distribution, metabolism, excretion, and toxicological (ADMET) properties is advisable early in the drug design process. To efficiently filter a library of compounds against such criteria, several online tools have been developed.

Chembioserver is a publicly available online application specializing in filtering and selection of small molecules [58]. The objective of this application is to facilitate compound preparation prior to (or after) VS computations by utilizing its many sections, such as (i) basic search, (ii) filtering (steric clashes and toxicity), (iii) advanced filtering based on custom chosen physicochemical properties, (iv) clustering (according to structure and compound physicochemical properties providing representative compounds for each cluster), (v) customized pipeline and (vi) visualization of compound' properties through property graphs and thus, increase the efficiency and the quality of compounds that proceed to in vitro assaying. The FAF-Drugs2 server is a public tool for computationally filtering compounds by considering ADMET and physicochemical properties, and identifying key functional or undesirable moieties [59]. Moreover, users have the opportunity to select filtering thresholds and rules among 23 physicochemical and 204 substructure searching rules. Last but not least, this tool provides numerous distribution diagrams of the properties of the filtered compounds in a web-server version [60].

Finally, compound datasets should be preprocessed in realistic 3D representations. Therefore, the compound set to be used for SBVS should have realistic bond lengths and angles as these may not change during docking. Also, it must be devoid of accompanying fragments such as counter-ions, metals, and solvent molecules, all compounds need to have assigned bond order and filled valences, partial charges, an appropriate protonation state at physiological pH or at the pH of the interest and proper tautomeric states [61, 62]. However, it was shown that consideration of all possible tautomers (stable and unstable) yields slightly poorer results than including only the most stable form in water [63]. Docking programs usually include ligand preparation software such as Autodock Tools [64], LigPrep within the Schrödinger suite [65], MOE [23] and the MAPS platform [66]. Moreover, stand-alone programs exist for ligand preparation such as DISI [67], Pipeline pilot [68], or Hyperchem [69]. Grinter et al. have created a workflow through custom-made script files for preparing large libraries for docking, including scripts for ligand preparation, which are offered freely to the academic community [70].

Library Design

Although numerous drug-like compound libraries are freely available online [1, 5], users may often need to create a custom-made library. Libraries may be divided in a) generic virtual high throughput (vHTS) libraries, which contain large sets of compounds, b) diversity-oriented, containing highly chemically diverse compounds, c) target-oriented, which are designed with a specific target in mind, d) molecular property diversity libraries, which are designed with specific molecular property profiles (i.e. solubility, lipophilicity etc.), e) natural product libraries.

Online tools for library design include CLEVER [71] that supports chemical library manipulation, combinatorial chemical library enumeration using user-specified chemical components, chemical format conversion, as well as chemical compound analysis and filtration with respect to drug-likeness, lead-likeness, and fragment-likeness based on the physicochemical properties computed from the derived molecules. Also provided is an integrated property-based graphing component that visually depicts the diversity, coverage and distribution of selected compound collections. Another tool for library design is e-LEA3D, which performs combinatorial library design that is based on a user-drawn scaffold and reactants coming, for example, from a chemical supplier [72]. More freely-available as well as commercial software packages are listed in Table . A review on how to design target-focused compound libraries, e.g. for kinases, ion channels, GPCRs, protein-protein interfaces can be found in Ref. [73]. In another approach, by using a genetic algorithm, Lisurek et al. have identified substructures typically occurring in bioactive compounds using the World Drug Index; subsequently available compounds containing the selected substructures were collected from vendor libraries [74].

In order to build a compound library with a range of possible replacements for a single substructure, Wirth et al. have developed a useful platform for the applicability of Matched Molecular Pairs (MMPs) to optimize the binding affinity of a small-molecule ligand to a target protein, namely the SwissBioisostere database [75]. In this scheme, MMPs are defined as two molecules, which differ in one particular substituent and exhibit different properties. The underlying assumption of MMPs is that the difference in properties can be extrapolated to another pair of molecules exhibiting the same substitution pattern. Weber et al. have improved the bridging between MMPs and structural data by developing the Virtually Aligned Matched Molecular Pairs Including Receptor Environment (VAMMPIRE) database [76]. In their work, they have proposed that the involved atoms in the ligand-protein complexes play a critical role in the exchange of the substituent and therefore in the binding affinity change, providing the possibility to extrapolate activities from one biological system to another. In seek of designing diverse libraries to enrich the screening collection, Vainio et al., have constructed “Virtual Library”, a virtual automatic system that uses validated synthetic protocols and available starting materials to generates a large number of virtual compound libraries [77].

An informative analysis and novel insights on preparing a high-quality database of lead-like compounds for SBVS are presented in Ref. [78]. The application of functional group, lead-like, PAINS, and Tanimoto similarity exclusion filters is discussed resulting in a library of ~350,000 compounds that represent 80% of available leadlike space based on vendor information in what is proposed to be the globe’s highest quality collection of available screening compounds. The article is a nice review of filtering process steps also providing detailed information on the filters and pointing out several issues as resupply, minimum concentration to be tested, etc., which is of practical use to medicinal chemists.

Constructing libraries with inherent chemical diversity such as in Ref. [78] may be desired when investigating underexplored targets with few known ligands. However, the construction of targeted libraries based on properties of known ligands may also be required in order to identify hits with improved potency. Xing et al., designed combinatorial libraries to search for novel soluble epoxide hydrolase (sEH) inhibitors based on a benzoxazole template forming conserved hydrogen bonds with the catalytic machinery of sEH [79]. Consequently, screening of these libraries resulted to 90% hit rate and more than 300 submicromolar sEH inhibitors were finally discovered.

Docking & Scoring

A large number of docking programs have been developed recently, including AutoDock [64], Dock [80], FlexX [81], Glide [82], Gold [83], Surflex [84], ICM [85], LigandFit [86], Drugster [87], and eHiTS [88] (also reviewed in [5]). Docking entails predicting the protein-ligand complex structure and is followed by scoring in SBVS in order to rank the compounds. Docking programs utilize various methods of conformational search in order to explore the ligand conformational space; these are categorized as following: a) Systematic methods, which place ligands in the predicted binding site after considering all degrees of freedom, b) Random or stochastic torsional searches about rotatable bonds, such as Monte Carlo and genetic algorithms to “evolve” new low energy conformers, (c) Molecular Dynamics simulation methods and energy minimization for exploring the energy landscape of a molecule [reviewed in Ref. 89].

In order to rank compounds, docking programs utilize scoring functions that aim to estimate the free energy of binding of a ligand to a specific target based on a generated docked pose after docking different ligands of a database. A variety of scoring functions have hitherto been developed (reviewed in [1, 5]). Commonly-used scoring functions can be categorized as follows: (a) Force field-based functions that estimate the binding free energy by summing the strength of intermolecular van der Waals, electrostatic interactions and hydrogen bonding between all atoms of the two binding partners in the complex. Solvation and entropy contributions are also taken into account. (b) Empirical scoring functions that are based on counting the number of various types of interactions between the two binding partners, i.e. hydrophobic contacts, number of hydrogen bonds and number of rotatable bonds immobilized in complex formation. These functions have proven to be successful for many protein-ligand complexes. LUDI, FlexX, F-Score, ChemScore and Fresno are a few examples of functions that are contained in this group. (c) Knowledge-based functions that use statistical observations of intermolecular contacts in receptor-ligand complexes with known structural conformations. For instance, some functions of this category are Potential of Mean Force (PMF), DrugScore and SMoG (Small Molecule Growth). Furthermore, there exist scoring functions that combine two or more of the above mentioned scoring function categories.

It is commonly accepted that while docking results to successful binding pose prediction, scoring usually fails to correctly rank different compounds with the difficulty increasing in congeneric series [95, 96]. Hence, identifying correctly the right binding pose as the top-ranked one still remains a challenge [1]. Addressing the correspondence between predicted and experimental binding affinities in empirical scoring functions still remains unresolved. Based on that critical issue, Zilian et al. have developed a new scoring function, SFCscoreRF, with significantly improved performance in comparison to previously developed SFCscore functions, by using SFCscore descriptors and a PDBbind set of 1005 complexes as a training set in combination with random forest for regression [97]. HYDE is another novel scoring function for protein-ligand complexes which is based on HYdration and DEsolvation terms. It consistently describes hydrogen bonds, the hydrophobic effect and desolvation. Recently, Schneider et al. have validated the HYDE scoring function through large-scale docking experiments, which resulted to successful prediction of the correct binding mode in 93% of complexes when checked with the Astex diverse set [98]. Alternatively, Ravindranathan et al. have defined a physics based scoring function and more precisely, a variable dielectric model based on residue types for better description of protein–ligand electrostatics in MM-GBSA scoring, which results in a higher correlation with affinity data [99]. Certainly, while much progress has been made in delivering more accurate scoring functions, further improvements in this direction are desirable.

The assessment of the efficiency and accuracy of scoring functions has been addressed in a number of studies. Cross et al., have assessed six molecular docking programs (DOCK, FlexX, GLIDE, ICM, PhDOCK, and Surflex) in order to evaluate docking pose, scoring and thus VS accuracy [90]. The results revealed that ICM, GLIDE, and Surflex reproduced accurately the X-ray poses, while GLIDE and Surflex showed superior receiver operating characteristic (ROC) areas under the curve (AUC) and ROC enrichments after performing a VS exercise on the 40 protein targets in the Directory of Useful Decoys (DUD). McGann et al. have compared the performance of FRED and HYBRID docking programs on known datasets, where the second utilizes a modified version of FRED that uses both ligand- and structure-based information for docking [91]. The results have shown that HYBRID scoring function increases identified actives. Especially, when using multiple crystal structures considering protein flexibility, HYBRID efficiency increases even more without increasing significantly the docking time (~ 15 %). Liebeschuetz et al. have tested the four GOLD scoring functions in docking and scoring computations and demonstrated that the ChemPLP scoring function comes first for both pose and affinity prediction [92].

In an attempt to improve docking and scoring performance, Spitzer et al. used Surflex-Dock to incorporate protein flexibility in the docking process. Ligand-protein optimization before VS, resulted in large differences in the VS performance, indicating that ligand re-docking for pose prediction assessment may have limitations [93]. Docking protocols, which incorporate protein flexibility, performed better in terms of pose prediction. Applying various VS protocols (docking, 2D molecular similarity, and 3D molecular similarity) resulted in achieving better hit rates. Neves et al. have investigated the ICM flexible docking and scoring benchmarks for 85 co-crystal structures of the modified Astex data set and derived significant improvements, compared to single rigid docking [94]. Reparsky et al. also have obtained enrichment results for WScore, a new scoring function and sampling methodology incorporating WaterMap results into Glide screening [27]. In this study, WScore results show early enrichment and outerperform the DUD best-practices Glide SP results.

Improving Pose/Compound Selection After Docking (Post-processing)

A rate-limiting step in SBVS is often the need for a computational chemist expert to post-process compounds that result from a VS/docking exercise before selecting the ones that will pass to the experimental test phase. The implementation of the simplified scoring functions and sometimes the inadequate sampling of the conformational space for the ligand may lead to unrealistic poses, intra-ligand steric clashes, twisted amides, E/Z esters, imperfect hydrogen-bonding network, and poses based on shape complementarity. These poses may result to an unreasonably high score and need to be discarded. Therefore, visual inspection of thousands of docking poses is normally needed by the medicinal chemist in order to select the appropriate compound set for assaying. To this end, significant efforts have been dedicated to increase the efficiency and the quality of compound selection [58, 106].

In the process of post-processing the results, one may resort to rescoring. To examine the reliability of VS results without post-processing, Malmstrom et al. have investigated a set of known binders and non-binders of L99A T4 lysozyme using both free energy of binding (FEB) rescoring and other empirical docking scores and concluded that FEB rescoring may assist lead identification resulting from a SBVS exercise [107]. Ding et al. have introduced a directional approach, MIECs (molecular interaction energy components), which successfully predicted the binding affinities of protein-ligand complexes or distinguished binding from non-binding ligands [108]. This methodology was extended to SBVS. HIV-1 protease was used as a model system to assess the performance of the method on ranking docking poses and docked ligands and a support vector machine (SVM) was trained in evaluating the energetic and geometric characteristics that discriminate binding ligands from the non-binding ones. As previously-mentioned, several tools exist for the pre-processing of ligand libraries prior to the VS procedure. These tools may be also used for post-processing of SBVS results are readily available. The ChemBioServer for example, uses vdW filtering to remove compounds with steric clashes. Poses that are far from the energy minimum are unlikely to be adopted in nature and hence, should be discarded. Compounds that pass vdW filtering may be then subjected to more stringent physicochemical property filtering compared to the initial compound selection for SBVS. Subsequently, hierarchical clustering may be performedin order to group compounds with similar structures/physicochemical properties and derive subsets with maximal chemical diversity [58].

Improvements in SBVS Efficiency

Ensemble Docking (ED)

Selecting the ideal crystal structure of a receptor target is the first step in the SBVS process. Unfortunately, crystal structures provide a single conformation of the protein, which is influenced by the crystallization conditions, and do not include any information on protein dynamics. Furthermore, crystal structures are highly affected by ligand binding, which leads to conformational changes in both the protein and the ligand. Proteins can exhibit induced fit effects upon binding of a ligand, in which the protein conformation changes significantly. Hence, although crystal structures are a preferable starting point for SBDD, in some cases they provide misleading information. Therefore, many efforts have been made to include receptor flexibility in docking programs, as it confers a more realistic depiction of the modeled system [100]. In order to avoid biasing towards one protein conformation and implicitly including protein flexibility, one may use an ensemble of representative structures to dock candidate ligands. ED comprises docking a single ligand library against multiple rigid receptor conformations, contrary to the standard single rigid receptor docking methods. In this way, ED mimics the realistic structural variation of the protein-target [101].

Currently, a huge amount of protein structures that can be used as biological targets for specific diseases is available. Therefore, the challenge in SBVS has changed from “Is there a promising structure available?” to “How do I appropriately select a single or multiple protein structures for my screening exercise?”. Consequently, standardized protocols that specialize in the selection of suitable conformations for ensemble-based virtual screening and finally result to potent leads of the target are needed. For instance, Rueda et al. conclude that the largest binding sites in X-ray structures lead to the highest enrichments [102]. Osguthorpe et al. apply a shape-based characterization of the binding site using molecular dynamics (MD) to construct ensembles with structural diversity [103]. The ensembles for the performance of ED may be selected through X-ray crystallography, NMR spectroscopy or a combination of both and through computational techniques (molecular dynamics simulations, homology modeling) [100].

There are certain advantages that distinguish ED among other flexible-receptor docking methods, i.e. the capability of accounting for any type of protein motion. The major disadvantage of ED is the computational effort, which scales linearly with the number of receptor conformations in the ensemble. As a result, the ED calculations have a significantly longer duration comparing to standard docking protocols. Moreover, possible inaccurate results being output from ED computations may be due to the possible imperfections in scoring functions, which may lead to inaccurate predictions of certain protein-ligand interactions [101, 103].

Consensus Induced-Fit Docking

Consensus Induced Fit Docking (cIFD) is an alternative docking approach, which allows the protein binding site to adapt to multiple ligands during SBVS. This methodology has been practiced by Kalid et al., firstly by validating the cIFD protocol on COX-2, the estrogen receptor, and HIV reverse transcriptase, which were previously shown to be challenging for docking programs and then by proving the utility of cIFD in discovering novel irreversible Crm1 inhibitors [104]. Houston et al. have introduced Consensus Docking for improving the probability of identifying accurately docked poses [105]. In cIFD, a similar approach to consensus scoring schemes is applied; this methodology combines information about predicted binding modes rather than predicted binding affinities. It results in 82% success rate of pose prediction by using more than one docking program, while the success rate was found to be lower for each docking program alone, e.g. 55% for Autodock, 58% for DOCK, and 64% for Vina.

PROPOSED WORKFLOWS-PROTOCOLS

Nasr et al. have introduced a set of “binding site properties-based” guidelines that aim to optimize future prospective drug discovery protocols in presence of several available receptor structures [109]. Several trends in the properties of binding sites, including hydrophobicity, volume, and opening, were identified and tested in SBVS with several programs (Surflex-dock or ICM). Planesas et al. have proposed a three-step VS protocol to improve actives from decoys selection using VEGFR-2 (Kdr-kinase) inhibitors. This protocol includes the conventional docking step, a pharmacophore post-filter step, and a similarity search post-process [110]. When this protocol is retrospectively applied to VEGFR2 inhibitors, both the overall and early VS performance are improved. Osguthorpe et al. have generated a SBVS protocol through which a set of target structures is constructed for ED based on binding site shape characterization and clustering, in order to account for protein flexibility into the VS process [103]. The protocol’s steps are: a) receptor structures selection from X-ray crystallography and MD, b) binding site shape characterization, c) pairwise volume overlap computation, d) hierarchical clustering, e) docking calculations of known actives and decoys to the shape-diverse binding sites, and f) comparison of ED results with the performance of a single rigid receptor.

One approach to achieve nM hits directly from VS is to bias the compound selection toward known chemotypes of high potency. This high potency may be achieved through the design of target-focused libraries for virtual screening. Schröder et al. discovered covalent reversible inhibitors of cathepsin K through virtual screening of a focused library built for this purpose, which included compounds bearing electrophilic groups that act as warheads in mediating covalent reversible inhibition [8]. Substantial effort was invested in developing appropriate electrophilic warheads in their corresponding transition state. Covalent docking of this focused library resulted in identifying 21 inhibitors of cathepsin K, including 3 nanomolar leads, out of the 44 initially-selected putative binders.

In another SBVS for β2-adrenergic receptor ligands, virtual screening was performed on the lead-like subset of the ZINC database using the available crystal structure [9]. 0.05% of the top-ranked compounds (500 molecules) were visualized, and 25 molecules were finally selected based on their chemical diversity, commercial availability, and overall balance between polar and non-polar complementarity to the binding site. These efforts resulted in the identification of a9 nM inverse agonist of the β2-adrenergic receptor. Although the virtual screening procedure was unbiased, it was later found that the lead-like subset of the ZINC database is biased towards GPCR ligands.

A novel approach targeting GPCR inhibition, has unraveled nanomolar leads for the melanin-concentrating hormone-1 receptor (MCH-1R) [7]. This approach integrates GPCR modeling, prediction of antagonist binding site, design, synthesis, in combination with utilizing a focused library for screening. A primary hit compound from a pyranose-based VAST library (pyranose scaffolds mimic active GPCR peptides) was initially utilized to construct a high quality MCH-1R model. Subsequently, the model was validated by a virtual enrichment experiment and the model-driven structure-based expansion of the original hit was used to identify key interactions in the binding site of the protein. Following the model validation, a structure-based virtual screen of a library with ≤ 0.7 Tanimoto similarity to known MCH-1R ligands, provided a 14% hit rate and 10 novel chemotypes of potent MCH-1R antagonists, including two nanomolar leads. Ligand selection also took into consideration AMBER interaction energies in the range of the top 10% of known MCH-1R antagonists.

In all the above-mentioned approaches, a high hit rate and high affinity binders were achieved by using a focused library to identify lead candidates.

Apart from identifying novel lead candidates, other applications of SBVS may be envisaged such as the screening for novel substrates for genetically-engineered proteins. Recently, it has been shown that new specific substrates of genetically-engineered proteins can be used for identifying the direct protein targets of enzymes [111-113]. In this context, SBVS could be applied for the identification of substrates specifically-designed to target the mutated enzymes that are then used to elucidate signal transduction pathways. In this direction, the group of Matt Jacobson has recently unveiled a strategy for assigning valid biological functions to proteins identified in genome projects using SBVS [114]. Computer-aided strategies for functional discovery with ‘metabolite docking’ were used to predict substrate specificities of several enzymes encoded by a bacterial gene cluster. This work led to the correct prediction of the in vitro activity of a structurally characterized enzyme of unknown function (PDB 2PMQ), and also to the correct identification of the catabolic pathway in which Hyp-B 2-epimerase participates. These studies establish the utility of structure-guided functional predictions to enable the discovery of new metabolic and signal transduction pathways.

Another emerging application of SBVS is the search for inhibitors of multifunctional proteins. Stumpfe et al. [115] report a virtual screening approach for finding inhibitors of the cytohesins, a family of cytoplasmic proteins with multiple known functions. Two virtual screening approaches, fingerprint similarity searching, and support vector machine modeling were combined to identify structurally diverse compounds capable of inhibiting the different functions of the cytohesins. Specific inhibitors of guanine nucleotide exchange, insulin signaling, and leukocyte adhesion were identified from the screen.

Protocol for the Discovery of Mutant-Specific PI3Kα Inhibitors

The PIK3CA gene is one of the most frequently mutated oncogenes in human cancers. It encodes p110α, the catalytic subunit of phosphatidylinositol 3-kinase, α (PI3Kα) isoform. The PI3Kα protein is implicated in signaling cascades, which lead to cell proliferation, survival, and cell growth. One of the most frequent mutations in PI3Kα is a histidine changed to arginine in exon 20 (H1047R) [116]. Here, we describe a VS workflow to identify novel inhibitors of the H1047R mutant form of PI3Kα [95]. The workflow incorporates the use of ChemBioServer [58], a free web-based application, which enables the efficient final selection of compounds to be tested experimentally.

Initially, the crystal structure of the mutant H1047R PI3Kα (PDB ID: 3HIZ) was constructed using a combination of homology and loop modeling and the atomistic model of the full-length H1047R mutant was created using Modeller 9v8 [117]. The resulting model was solvated in water and employed in Molecular Dynamics (MD) simulations using the NAMD package [118] for ~70ns. The final snapshot of the simulation was extracted and used for binding site identification calculations through SiteMap module of Schrödinger v2.4. A non ATP binding site close to the H1047R mutation was identified to be among the top-ranked and was used in the present study. After binding site identification, we performed VS using the docking program Glide 5.7 [82]. The drug-like subset of the HitFinder collection from the Maybridge database (www.maybridge.com) was used for the VS exercise [119]. All structures were docked and scored using the Glide standard precision (SP) mode. The 10,000 top-ranked structures from the SP filter were redocked and rescored using the Glide extra precision (XP) mode. The complexes for the top-ranked 1,000 compounds resulting from the XP processing were submitted to further post-processing with the ChemBioServer. In order to filter out poses far from the energy minimum, we used the van der Waals filter of the server to remove compounds with steric clashes. The remaining compounds were then subjected to physicochemical property filtering based on the “Jorgensen rule of 3” [52] as well as toxicity filtering based on a database available in ChemBioServer, which contains known toxic moieties. Subsequently, a hierarchical clustering was performed for the remaining compounds using the Tanimoto coefficient and the Ward Clustering Linkage in order to maximize chemical diversity but at the same time minimize visual inspection efforts. Finally, the resulting cluster representatives were visually inspected and the ten most promising compounds were purchased and submitted to in vitro assay testing. Four out of the ten purchased molecules inhibited the PI3Kα protein activity in vitro in μm concentrations, indicating that the workflow described herein can be successfully applied in drug discovery [120]. The process is described in Fig. ().

A Proposed Method to Increase the Probability of Virtual Screening Hits: Application on the RXRα Nuclear Receptor

The availability of multiple solved protein structures of the same protein may offer the opportunity to increase the hit rate in biological assays of compounds selected in silico. The case of receptors in particular, which represent the majority of targets for pharmaceutical intervention, could be good examples for such approaches since several crystal structures of a particular receptor with agonists or antagonists are available. We propose to use such available structures in combination with the protocol developed below in order to increase the rate of compound identification that can indeed modulate the target in question.

Here, we use the nuclear receptor Retinoid X receptor alpha (RXRα) as an example for a selective SBVS screening protocol, which is based on ensemble docking in order to discover novel selective agonists. Retinoid X receptor (RXR) is a heterodimer partner for about one third of the 48 human nuclear receptor superfamily members, thus, its ligand binding pocket may be altered upon binding to a different heterodimerization partner. Heterodimerization greatly enhances DNA binding and transcriptional activation [121] and RXRα overexpression has been reported to enhance the transcriptional response to ligand binding. RXR-selective retinoids can potentiate the antiproliferative and apoptotic responses of breast cancer cell lines to PPAR ligands. Recently, RXR-selective ligands were discovered that inhibited proliferation of all-trans RA resistant breast cancer cells in vitro and caused regression of the disease in animal models. Therefore, RXRα is potentially a therapeutic target for cancer [121, 122].

Based on the fact that multiple RXRα receptor conformations exist depending on binding of the receptor to different heterodimerization partners, we seek to increase the specificity of identified binders for a given heterodimer partner by carefully choosing the RXRα receptor structure for SBVS. In our method, the choice of RXRα conformations for the SBVS is based on four criteria: a) Pairwise comparison of the receptor conformations according to RMSD calculations, b) analysis and clustering of RXRα structures comparing the binding-site shape and volume using SiteMap (Fig. ), c) docking of a small-database of known actives for a specific heterodimer partner to the resulting shape-diverse subset of binding sites resulting from (a) and (b) using Glide 5.8 SP and XP [82], d) retrieving representative protein conformations for the structure of interest from MD simulations using GROMACS [123]. VS performed on three different subsets of RXRα receptor conformations, which arise by binding to different heterodimerization partners, selected as mentioned above, may enhance the success rate of the process (Fig. ).

Docking compound databases for this SBVS exercise entails pre-processing to assign protonation and tautomeric states as described in [119]. Compounds for assying are selected as follows: Molecules that score high when docked in the RXRα protein ensemble that binds to the heterodimer partner of interest and at the same time score low for RXRα structures that bind to heterodimer partners of no interest, may be selected in order to achieve selectivity. Finally, a post-processing step is imposed to the top-scoring compounds by using Chembioserver [58] and FAF-Drugs2 [60] filtering tools as well as pharmacological property prediction with the QikProp software [124]. The workflow of this protocol is shown in Fig. ().

CONCLUSION

For decades, drug discovery was carried out using trial and error experimental techniques for screening large libraries of chemicals against a biological target. Recent advances in computer-aided drug design allow the tailored design of drugs for a target protein, shortening the development cycle of new drugs. The advent of Structure Based Virtual Screening has undoubtedly changed and improved the drug discovery process and has been established as one of the most promising in silico techniques for drug design. The focal point of this review is the detailed description of the SBVS steps in drug design, in combination with the presentation of recent advances and the introduction of various protocols that facilitate the identification of inhibitors with nM potency.

Furthermore, we also propose herein two novel VS protocols that aim to enhance inhibitor selectivity for the target protein against close homologs. First, a SBVS workflow, which was utilized to discover novel inhibitors of the H1047R mutant form of PI3Kα is presented. In this workflow, MD simulations in aqueous solution were carried out for both WT and mutant PI3Ka proteins. Binding site analysis in the kinase domain identified cavities in the vicinity of the H1047R mutation and the membrane binding regions. SBVS was performed in several binding pockets and top-ranked compounds in terms of predicted binding affinity were carefully post-processed to ensure the validity of docked poses, chemical diversity, desirable physicochemical properties and the absences of toxic or metabolically liable moieties. Finally, ten promising compounds were selected for in vitro assaying, four of which emerged as µM inhibitors of the H1047R mutant PI3Kα protein validating our approach [120]. A second SBVS protocol is contributed for the identification of binders for the RXRα nuclear receptor that are selective depending on the protein’s heterodimer partner. A structural RXRα ensemble was created for this purpose by selecting different RXRα crystal structures based on (i) RMSD calculations, (ii) binding-site shape and volume, (iii) docking of a small database of known actives, and (iv) choosing representative structures from MD simulations. SBVS was then performed on three different subsets of RXRα arising from different heterodimer complexes of RXRα and/or its binding to different ligands. Candidate binders of RXRα were selected for purchase with an eye on their different orientation at the binding site of the various structures and different interactions with specific surrounding residues in order to maximize their selectivity potential. In vitro assying of these compounds is still pending experimental testing.

Although SBVS is widely used nowadays by multiple academic and industrial research groups in the drug discovery process, it suffers from limitations that restrict its effectiveness (reviewed in Table ). Significant breakthroughs are required in order to address fundamental challenges such as, for example, scoring, target flexibility and appropriately treating water molecules. Such challenges ultimately lead to the query: Is SBVS an indispensable tool for modern drug design? The significant reduction in time and cost compared to the high-throughput screening process, the continuous efforts in improving the efficiency of docking programs and scoring functions, and the plentiful successful case studies that have led to low nM leads are only a few representative examples showing that SBVS is here to stay. However, prospective users of the method should remember that SBVS is not as simple as running a computer program. Careful choices need to be made; appropriately selecting the structural ensemble for the screening exercise, cautious preparation of the biological target and the database to be used, treatment of water molecules in the cavity, and careful post-processing of SBVS results are of utmost importance that ultimately result in enhancing lead identification and selectivity rates.

CONFLICT OF INTEREST

The authors confirm that this article content has no conflict of interest.

Table 1.

Free available software and tools for performing a SBVS workflow.

Software	Free Academic License	Website
GLARE [130]	Yes	http://glare.sourceforge.net/
CDK-Taverna [131]	Yes	http://cdktaverna.wordpress.com/
CLEVER [71]	Yes	http://datam.i2r.a-star.edu.sg/clever/
e-LEA3D [72]	Yes	http://chemoinfo.ipmc.cnrs.fr/
SmiLib v2.0 [132]	Yes	http://gecco.org.chemie.uni-frankfurt.de/smilib/
Library synthesizer (Tripod) [133]	Yes	http://tripod.nih.gov/?p=370
Swissbioisostere [75]	Yes	http://www.swissbioisostere.ch/
VAMMPIRE [76]	Yes	http://vammpire.pharmchem.uni-frankfurt.de/#!home
Virtual Library [77]	Yes	Contact authors
Pipeline pilot (Accelrys) [68]	No	http://www.accelrys.com
Reactor [134](ChemAxon)	No	http://www.chemaxon.com
OELib library enumeration [135]	No	http://www.eyesopen.com
CombiLibMaker and Legion (Tripos) [136]	No	http://www.tripos.com
QuaSAR-CombiGen [137]	No	http://www.chemcomp.com
ChemOffice CombiChem [138]	No	http://www.cambridgesoft.com/
ICM-Chemist [139]	No	http://www.molsoft.com
LUCIA [140]	No	http://www.eidogen-sertanty.com/

Table 2.

Free and commercial software packages for library design.

Software	Free Academic License	Website
GLARE [130]	Yes	http://glare.sourceforge.net/
CDK-Taverna [131]	Yes	http://cdktaverna.wordpress.com/
CLEVER [71]	Yes	http://datam.i2r.a-star.edu.sg/clever/
e-LEA3D [72]	Yes	http://chemoinfo.ipmc.cnrs.fr/
SmiLib v2.0 [132]	Yes	http://gecco.org.chemie.uni-frankfurt.de/smilib/
Library synthesizer (Tripod) [133]	Yes	http://tripod.nih.gov/?p=370
Swissbioisostere [75]	Yes	http://www.swissbioisostere.ch/
VAMMPIRE [76]	Yes	http://vammpire.pharmchem.uni-frankfurt.de/#!home
Virtual Library [77]	Yes	Contact authors
Pipeline pilot (Accelrys) [68]	No	http://www.accelrys.com
Reactor [134](ChemAxon)	No	www.chemaxon.com
OELib library enumeration [135]	No	http://www.eyesopen.com
CombiLibMaker and Legion (Tripos) [136]	No	http://www.tripos.com
QuaSAR-CombiGen [137]	No	http://www.chemcomp.com
ChemOffice CombiChem [138]	No	http://www.cambridgesoft.com/
ICM-Chemist [139]	No	http://www.molsoft.com
LUCIA [140]	No	http://www.eidogen-sertanty.com/

Table 3.

Advantages and Drawbacks of SBVS.

Virtual Screening
Advantages	Limitations
Time and cost reduction of screening process of millions of small molecules, compared to HTS	Many VS tools are applicable and successful to specific case studies (based on the training set) and not in general cases.
There is no need for physically existing compounds to perform the screening process, unlike HTS.	Compounds being identified by HTS are usually more bioactive than compounds identified by VS.
Different approaches of VS have been created for lead discovery depending each time on the availability of experimental information (SBVS Ligand-Based VS, Fragment-Based VS,etc.)	Weakness in perfect inclusion of receptor structural flexibility and of water in docking computations due to computational-cost and high complexity of its modeling
Several successful examples of identifying low nM leads that show the intended biological activity	Very potent leads (i.e. low nM) are rarely identified through VS.
A large number of docking programs and scoring functions	Scoring is still challenging in predicting accurately the correct binding pose and ranking of the compounds due to the difficulty in parameterizing the complexity of the ligand-receptor binding interactions and the approximations in calculating desolvation and entropic terms.
VS can use as input a desirable target structure complexed with a specific ligand even if there are no experimental data, through molecular modeling.	Predicted protein structures from homology modeling and predicted protein-ligand complexes may result to increased rates of false positive/negative results.
	Does not perform in congeneric series.