Hierarchical virtual screening approaches in small molecule drug discovery.

Virtual screening has played a significant role in the discovery of small molecule inhibitors of therapeutic targets in last two decades. Various ligand and structure-based virtual screening approaches are employed to identify small molecule ligands for proteins of interest. These approaches are often combined in either hierarchical or parallel manner to take advantage of the strength and avoid the limitations associated with individual methods. Hierarchical combination of ligand and structure-based virtual screening approaches has received noteworthy success in numerous drug discovery campaigns. In hierarchical virtual screening, several filters using ligand and structure-based approaches are sequentially applied to reduce a large screening library to a number small enough for experimental testing. In this review, we focus on different hierarchical virtual screening strategies and their application in the discovery of small molecule modulators of important drug targets. Several virtual screening studies are discussed to demonstrate the successful application of hierarchical virtual screening in small molecule drug discovery.

Introduction

Modern drug discovery process starts with the identification of initial hits that are further optimized to improve the potency, selectivity, metabolic stability and oral bioavailability. Among many ways of identifying initial hits in drug discovery, high-throughput screening (HTS) and virtual screening (VS) are most common. The VS was originally developed to bring down the cost of discovering new molecules using HTS. In the last two decades, advances in computational programs and processing power have made VS an important tool to identify starting points, inhibitors and chemical probes in various drug discovery campaigns [1], [2], [3], [4].

In light of the immense potential of VS methodologies in the identification of initial hits, various structure and ligand-based approaches were developed [1], [5], [6], [7], [8]. Structure-based methods rely on the structural information of the protein target and typically include methods such as molecular docking [9], [10], structure-based pharmacophores [11], [12] and de novo design [13], [14]. Ligand-based methods can work in the absence of structural information for the protein target. These methods include two or three-dimensional (2D or 3D) similarity searches [15], [16], ligand-based pharmacophore screenings [5], [17], machine learning approaches [18], [19], quantitative structure activity relationships (QSAR) [5], [20] among others. Ligand-based methods, however, require the availability of at least one known active molecule. Although utilizing these structure and ligand-based methods individually have demonstrated immense potential in retrieving initial hits, these methods are unable to fulfill all the practical requirements of drug discovery alone. Furthermore, with ever-increasing screening library size [21] and computational cost associated with some VS approaches especially flexible molecular docking [22], [23], [24], [25], it is indispensable to integrate different VS approaches to filter compounds. Ligand and structure-based methods can be combined in a sequential or parallel manner (Fig. 1 ). The most common way of combining these methods is to use them in a sequential funnel like manner commonly known as hierarchical VS (HLVS). In HLVS, a large small molecule library is reduced to a number of compounds that is small enough for biological assay by applying a series of filters (generally two or three) sequentially (Fig. 1A). In contrast to HLVS, there is parallel virtual screening (PVS) where several complementary methods are run in parallel and the best hits ranked according to each method are selected for biological testing (Fig. 1B). Although the retrospective analysis of literature data has shown the successful application of PVS [26], [27], [28], [29], [30], only a few applications in real world scenario could be found [27], [31], [32].

Fig. 1

Integration of ligand and structure-based approaches. (A) Hierarchical virtual screening (HLVS): series of filters (here similarity search, pharmacophore and molecular docking) are sequentially applied to bring down the number of compounds to be cherry-picked for biological assay. (B) Parallel virtual screening (PVS): ligand and structure-based filters are performed independently on the same or similar number of compounds.

In this paper, we review the current status of commonly used HLVS approaches and try to understand their utility in drug discovery campaigns of important therapeutic targets. Although the scientific literature is inundated by hierarchical computational approaches and protocols, we have restricted our review on only those studies that were validated by experimental assays. In the following sections we will outline the types of HLVS approaches utilized in various small molecule discovery campaigns. Later, we will discuss recent cases of successful hit identification utilizing HLVS protocols. Finally, we will describe the usefulness of HLVS in discovering inhibitors of important drug targets.

Hierarchical combination of VS methods

Hierarchical combination of ligand and structure-based VS approaches generally involves sequential execution of dissimilar VS methods. Mostly, computationally inexpensive ligand based approaches such as similarity search and pharmacophore screening are used during initial steps of an HLVS protocol. Methods demanding comparatively high computational resources such as molecular docking and molecular dynamics (MD) simulation are used once the number of compounds to be screened decreases to a reasonable number. The final step in a majority of HLVS campaigns incorporates the visual selection of compounds by expert researcher commonly known as “cherry picking”. In this step, ranking from VS methods is combined with expert chemical intuition and with literature-based knowledge. The HLVS can be classified in three categories based on the combination of VS methods: ligand-based HLVS, structure-based HLVS and hybrid HLVS, which will be described in detail below. A few successful cases of small molecule discovery using these three classes of HLVS are summarized in Table 1 .

Table 1

A few cases of successful hit identification utilizing HLVS protocols.

Drug target	Role	Structure of best compound	Activity of best compound	HLVS methods used	Reference
LB-HLVS
Liver X Receptor	Cholesterol metabolism		Three fold activation at 10 μM	Shape similarity, Pharmacophore modeling	Temml et al. [42]
NAADP receptor	NAADP signaling		IC50 = 2 μM	Shape similarity and electrostatic potential matching	Naylor et al. [45]
Enoyl reductase	Antibacterials		IC50 = 0.3 μM	Shape similarity, electrostatic potential matching, 2D fingerprints	Hevener et al. [44]
Melanin-concentrating hormone receptor 1 (MCHr1)	Obesity		MCHr1 binding IC50 = 2.2 nM	Shape similarity, electrostatic potential matching, 2D fingerprints	Muchmore et al. [43]
Dihydrofolate reductase and thymidylate synthase	Ovarian cancer		IC50 = 42.5 and 26.7 μM against cDDP-sensitive and IC50 = 33 and 35.7 μM against resistant ovarian cancer cell line	Physiochemical properties and pharmacophore fingerprint similarity	Carosati et al. [46]
Human bitter taste receptors (TAS2Rs)	Asthma and infection		0.5 μM effective concentration	1D molecular descriptors, fingerprint similarity, pharmacophore modeling and shape similarity	Levit et al. [48]
Cytohesin	Vesicle trafficking and insulin signaling		IC50 = 3.1 μM, Binding Kd = 7.5 μM	2D fingerprints and support vector machine	Stumpfe et al. [50]
SB-HLVS
SUMO activating enzyme 1	Sumoylation pathway		IC50 = 11.1 μM	Molecular docking, molecular dynamics and MM-PBSA binding free energy	Kumar et al. [51]
Tubulin	Cancer		IC50 = 30 μM	Molecular docking, molecular dynamics and MM-PBSA and MM-GBSA binding free energy	Cao et al. [55]
HIV-1 protease	AIDS		IC50 = 14 μM	Molecular docking and molecular dynamics simulation	Kunze et al. [56]
Cytochrome P450 aromatase	Estrogen-dependent breast cancer		IC50 = 9.4 nM	Multistep hierarchical docking	Caporuscio et al. [65]
FabI	Bacterial fatty acid biosynthesis		IC50 = 3.4 μM	Molecular docking, molecular dynamics and MM-PBSA binding free energy	Hu et al. [152]
LBSB-HLVS
Serotonin transporter	Neurotransmission		Ki = 1.5 nM	2D fingerprints, ADMET filtering, 3D pharmacophore modeling and molecular docking	Gabrielsen et al. [66]
B-RafV600E	Ras/Raf/MEK/ERK signaling pathway		IC50 = 0.3 μM	SHAFT 3D ligand similarity and molecular docking	Kong et al. [67]
Protein kinase CK2	Cancer		85% inhibition at 10 μM	Bayesian modeling, pharmacophore modeling and molecular docking	Di-wu et al. [69]
Coagulation factor VIII	Anticoagulants		IC50 = 3.5 μM	Pharmacophore modeling and molecular docking	Nicolaes et al. [70]
SUMO specific protease 2	Sumoylation pathway		IC50 = 3.7 μM	Shape similarity, electrostatic potential matching and molecular docking	Kumar et al. [108]
Insulin-like growth factor-1 receptor (IGF-1R)	Cell growth, proliferation and apoptosis		IC50 = 57 nM	Pharmacophore modeling and molecular docking	Liu et al. [72]
DNA G-quadruplex	Cellular aging and cancer		Ability to bind and stabilize telomeric G-quadruplex shown using fluorescence and CD methods.	2D fingerprints, shape similarity and molecular docking	Alcaro et al. [89]
Cruzain	Cysteine protease, Chagas disease		IC50 = 48.8 μM	Shape similarity, molecular docking, molecular hologram QSAR	Wiggers et al. [111]

Ligand based HLVS (LB-HLVS)

LB-HLVS sequentially combines methods based on similarity search or compound classification. Similarity based techniques include methods accessing the similarity of one or a few experimentally identified hits with molecules in a large library in terms of their physicochemical properties [33], structural fingerprints [34], 3D-shape [35], electrostatic potential [36] and pharmacophore features [37] etc. Compound classification techniques include clustering [38], [39] or machine learning based methods such as Bayesian methods and support vector machines [18], [40]. Although most of the ligand-based VS methods are used either standalone or in combination with structure-based VS methods, only ligand-based methods were reported to have been effectively combined in prospective applications. Yao et al. [41] reported an efficient multi-step ligand-based VS protocol that included physicochemical property filtering, pharmacophore-based screening, protein–ligand interaction fingerprint similarity analysis and 2D-fingerprint structural similarity search. Their protocol significantly improved the hit rate when compared with individual methods. Among the prospective applications, LB-HLVS that included a combination of shape-based VS and pharmacophore modeling has been used by Temml et al. [42] to identify two agonist of liver X receptor. The reported agonist activated both subtypes of liver X receptor (LXR α and β). Shape similarity has also been combined with electrostatic potential matching in the discovery of melanin-concentrating hormone receptor 1 antagonist [43], Francisella tularensis enoyl-reductase inhibitors [44] and chemical probe for nicotinic acid adenine dinucleotide phosphate (NAADP) [45]. In another application [46], ZINC database [47] subset enriched with quinoxaline scaffold was filtered based on pharmacokinetic properties. The resulting molecules were again investigated for pharmacophore fingerprint similarity with known quinoxaline based inhibitors of folate cycle proteins. Associated biological assay resulted in three compounds, which interfered with dihydrofolate reductase and thymidylate synthase and reduced their levels. Levit et al. [48] integrated 1D molecular descriptors, 2D fingerprint-based molecular similarity, ligand-based pharmacophore models and a shape-based VS method to identify activators of human bitter taste receptor TAS2R14. Bayesian analysis has been used with pharmacophore modeling to identify inhibitors of breast cancer resistance protein (BCRP) [49]. VS was carried out against 2000 FDA-approved drugs and 19 drugs were found to exhibit significant effect on BCRP transport function. Another machine learning technique, support vector machine (SVM) in combination with fingerprint similarity search was also used in the identification of molecular probes for the study of cytohesins [50]. Ligand-based methods have the advantage of being computationally inexpensive and have the ability to retrieve more potent hits than the structure-based methods. However, most of these methods suffer from the limitation of being biased towards active chemical scaffolds and thereby generate less diverse hits. Moreover it is difficult to decide the right input query as it affects the quality of output results. In this scenario, hierarchical combination of ligand-based methods such as pharmacophore modeling and similarity-based methods assure chemical diversity with potency.

Structure-based HLVS (SB-HLVS)

Structure-based methods are also occasionally combined in hierarchical style. Structure-based methods are used standalone or in combination with ligand-based methods. Hierarchical combination of structure-based methods suffered from the limitation of being computationally expensive. However, recent computational developments that enabled VS to be executed on massively parallel computing architectures allowed successful implementation of CPU intensive approaches such as flexible molecular docking and MD simulation in SB-HLVS protocols. Using structure-based methods in hierarchical manner, our group identified compounds with biaryl urea scaffold as inhibitors of small ubiquitin-like modifier (SUMO) activating protein 1 (SUMO E1) [51], [52]. We docked a small molecule library to the ATP binding site in SUMO E1 using a two-step molecular docking strategy. Initially, a fast rigid protein docking was used to dock Maybridge small molecule library. The top ranking compounds were then re-docked using a method that incorporates both ligand and protein flexibilities. Top ranking hits from docking were subsequently prioritized for biological assay using the MD simulation-based binding free energy calculation. Similar hierarchical docking and molecular-mechanics Poisson–Boltzmann surface area (MM-PBSA) scoring [53], [54] procedure has been also used in the identification of novel tubulin inhibitors [55]. Recently, Kunze et al. [56] reported low molecular weight inhibitors of human immunodeficiency virus 1 protease subtype B utilizing SB-HLVS. MD simulation was used to identify transient surface pockets, one of which was translated into a structure-based pharmacophore query and used to screen compound databases. Molecular docking was then used to prioritize compounds for biological assay that resulted in the identification of two active low molecular weight inhibitors. Apart from sequentially combining multiple structure-based approaches, molecular docking based VS was sometimes performed in multi-step hierarchical fashion [52], [57], [58]. One such approach is implemented in Glide program [59], [60], [61], [62] which starts with coarse-grained docking followed by finer levels of sampling with more strict description of scoring functions. This hierarchical docking approach has been successfully used to identify highly potent inhibitors [57], [58], [63], [64], [65]. The main advantage of hierarchical docking is that it can rapidly screen large small molecule libraries without requiring high performance computing infrastructures. Fast and simple docking methods reduce a big library to a focused set with comparatively fewer molecules that are later docked using more computationally expensive methods.

Hybrid HLVS (LBSB-HLVS)

A hybrid approach that combines the ligand and structure-based HLVS (LBSB-HLVS) is the most commonly used among the three types of HLVS. LBSB-HLVS holds significant potential as it takes advantage of available structural and ligand information in hit identification. In general, fast computational approaches are used initially followed by computationally expensive approaches. Historically, LBSB-HLVS has enjoyed more successes than either LB-HLVS or SB-HLVS. Recent identification of serotonin transporter ligands [66] is an excellent example of LBSB-HLVS where five small molecule databases were screened using a VS workflow that combines 2D fingerprint based screening, ADMET filtering, 3D pharmacophore-based screening and molecular docking. In silico screened compounds were tested using in vitro screening and full binding assays, and 74 active compounds belonging to 16 structural classes were obtained. In another study [67], a hybrid 3D molecular similarity approach, SHAFTS [68] was combined with docking-based VS to identify nine low μM inhibitors of B-RafV600E. Analog search and subsequent chemical optimization of hits resulted in compounds with nM potency. Di-wu et al. [69] also employed a LBSB-HLVS comprising of Bayesian modeling, pharmacophore modeling and molecular docking to identify casein kinase 2 inhibitors. Nicolaes et al. [70] prefiltered a Chembridge compound collection for drug-likeness and subsequently used it to screen four pharmacophore models developed from known actives. Molecular docking was then used to prioritize a list of 1000 compounds for surface plasmon resonance based biological assay. Extensive biological characterization resulted in the identification of 9 compounds targeting the C2 domain of coagulation factor VIII with activities ranging from 2.1 to 19.9 μM. We have earlier developed an LBSB-HLVS workflow that utilized protein ligand interaction information to prioritize hits from ligand and structure-based VS. This workflow was employed to identify compounds with potent anti-tubercular activities [71]. In another excellent LBSB-HLVS study, Liu et al. [72] reported potent inhibitors of insulin-like growth factor-1 receptor (IGF-1R). Structure-based pharmacophore modeling combined with molecular docking resulted in two hits with IGF-1R inhibitory activity. These initial hits were subjected to substructure search to retrieve structural analogs that were prioritized using pharmacophore constrained molecular docking. Subsequent biological assays identified 15 more hits with good IGF-1R inhibitory activity. The implementation of LBSB-HLVS is particularly advantageous in cases where both ligand and structural information are available, as it takes into account all possible information to prioritize hits for biological assay. It has been already established that no single structure or ligand-based method performs well on all targets. The hierarchical fusion of ligand and structure based methods improve the VS performance by decreasing the number of false positives while increasing true positive hits.

Molecular docking is a key component in most HLVS

Molecular docking is one of the commonly used VS method and a central component in the majority of HLVS schemes. Molecular docking algorithms predict the orientation and conformation of ligands with respect to target receptor binding site. Hundreds of binding poses are generated that are evaluated by scoring functions. Last decade has witnessed tremendous growth in the development of molecular docking methodologies, algorithms and programs [9], [10]. These methods adequately use protein structural information to prioritize hits for biological assay. Some of these methodologies were effectively utilized to identify active hits from screening libraries, though there were reports that molecular docking failed to outperform ligand-based virtual screening methods even in the presence of protein structural information [73], [74], [75]. Most molecular docking programs are quite fast and are able to dock a small library in reasonable amount of time. However, docking a library of a few million compounds still requires huge computational resources and can be effectively performed only on massively parallel computer clusters. To circumvent this issue, HLVS approaches combining molecular docking with different levels of filtering have been proposed to speed-up the time consuming procedure of structure-based VS. Several of these approaches use structural (2D or 3D) similarity, shape or electrostatic similarity with the ligand to filter a large screening library. Others utilize ligand or structure-based pharmacophore modeling to reduce a large small molecule library to a reasonable number (a few hundred to tens of thousands of compounds) for molecular docking. In the following section, we will outline these fast VS methods that were effectively combined with molecular docking.

Similarity search – molecular docking

Similarity search is a simple and computationally inexpensive ligand-based method that evaluates the similarity of a library of small molecules with known ligands in terms of their 2D or 3D properties [15], [16]. Similarity analysis can be performed using topological indices [76], substructure search [77], [78], molecular fingerprints [79], [80], [81], molecular descriptors [82], shape and electrostatic properties [73], [83], [84], [85], [86]. Among all these similarity measures, substructure search is sometimes performed to screen large ligand libraries for the molecules that contain defined molecular fragment [87], [88]. 2D and 3D molecular fingerprints are routinely used to filter compounds for molecular docking in HLVS campaigns. The molecular access system (MACCS) fingerprint [81] is one of the most popular molecular fingerprint [15], [16]. MACCS fingerprint has been successfully integrated with molecular docking in several virtual screening studies [89]. Additionally, MACCS fingerprints are routinely used to improve the VS hit rate by identifying structural analogs of hit compounds [52], [90], [91], [92]. The extended connectivity fingerprints [80] from Accelrys and 2-point (TGD) and two 3-point (TGT and GpiDAPH3) pharmacophore-based fingerprints [93] from molecular operating environment (MOE) are other fingerprints that are commonly combined with molecular docking in several VS studies [88], [94], [95]. Steric complementarity between ligand and receptor is one of the key factors determining the strength of binding. This forms the basis of a hypothesis that two molecules with similar shape should have similar biological activities. Although the shape similarity alone was sufficient to retrieve highly potent binders [44], [45], [96], [97], [98], [99], the true potential of shape similarity was realized when it was combined with molecular docking in hierarchical manner. In one study, Acharya et al. [100] screened ZINC database [47] by coupling ROCS [73], [85] shape based approach with GOLD docking [101] method to identify small molecule inhibitors of chemokine receptor 5 (CCR5)-N terminus binding to gp120 protein. Xia et al. [102] employed HLVS protocol combining shape-based approach with molecular docking to discover selective 11β-hydroxysteroid dehydrogenase 1 (11β-HSD1) inhibitors. By performing VS against the Maybridge database, they identified seven hits showing IC50 values lower than 100 nM. Salo et al. [103] performed VS against the Chembridge library using shape-based filtering and flexible docking to identify two novel SIRT3 scaffolds. Similar VS approaches have also been used to identify novel non-steroidal FXR ligands, PPARγ partial agonists, Grb7-based antitumor agents and fungal trihydroxynaphthalene reductase inhibitors [104], [105], [106], [107]. In some studies, instead of using shape-based similarity alone with molecular docking, it has been complemented with approaches such as 2D similarity search, electrostatic potential matching and pharmacophore modeling. Alcaro et al. [89] merged compounds retrieved by 2D MACCS fingerprint and ROCS similarity search using seven known actives. Drug-likeness filter and molecular docking were used to further prioritize these compounds. Biological testing of 28 selected compounds resulted in a few compounds that interact with the human telomeric G-quadruplex. Recently, our group has identified 1,2,5-oxadiazoles as novel inhibitors of SUMO specific protease 2 (SENP2) utilizing a VS protocol that combines shape and electrostatic potential based similarity searches with molecular docking [108]. This inhibitor discovery process utilizing LBSB-HLVS is outlined in Fig. 2 . A similar HLVS approach was also used to identify inhibitors of biotin carboxylase [109] and F. tularensis enoyl-reductase [44]. In another study, Hamza et al. [110] employed a VS protocol that integrated a pharmacophore model, consensus molecular shape patterns of active ligands and molecular docking to identify inhibitors with μM potency against mycobacterial mycosin protease-1. Furthermore, Wiggers et al. [111] developed a multi-step cascade strategy that integrated shape based approaches with molecular docking and molecular hologram QSAR. This protocol successfully identified non-peptidic cruzain inhibitors with trypanocidal activity.

Fig. 2

Outline of the discovery process of novel SENP2 inhibitors utilizing LBSB-HLVS. A four million compound small molecule library was filtered based on shape and electrostatic similarity with a fragment query prepared from the conjugate of SUMO1 C-terminal residues and substrate protein lysine. Molecular docking further prioritized the hits that were tested using a FRET based assay. Biological testing revealed two scaffolds that were later optimized for potency by identifying analogs.

Pharmacophore – molecular docking

As recommended by International Union of Pure and Applied Chemistry (IUPAC), the term pharmacophore represents an “ensemble of steric and electronic features that is necessary to ensure the optimal supra-molecular interactions with a specific biological target structure and to trigger (or to block) its biological response” [112]. In recent years, several pharmacophore modeling approaches and programs were developed [113], [114], [115]. Pharmacophore modeling has been extensively used in drug discovery either as a standalone VS method [116], [117], [118] or in combination with other VS approaches such as molecular docking. With a goal to identify inhibitors of matrix metalloproteinase 2 (MMP2), Pizio et al. [119] created and validated the seven-feature pharmacophore query based on the co-crystal structure of MMP8. This query was used to screen a library of ∼300,000 drug-like compounds. Molecular docking was then used to rank-order hits from pharmacophore-based screenings. Enzyme inhibition assay revealed nine active compounds. In another study, novel PKR-like endoplasmic reticulum kinase (PERK) inhibitors were discovered using hierarchical combination of pharmacophore screening and molecular docking [120]. Recently, similar VS protocols have been used to identify inhibitors of SARS-CoV 3-chymotrypsin-like protease [121], carbonic anhydrase VII [122], acetohydroxyacid synthase [123], malarial cysteine protease [124] and 5-lipoxygenase-activating protein [125]. Occasionally, studies combining structure and ligand based pharmacophore modeling have been published [126], [127], [128]. Using p53-MDM2 complex co-crystal structures, Xue et al. [126] built a structure-based pharmacophore model that was refined using a ligand-based pharmacophore model. The consensus model was used to screen a library of ∼239,735 compounds. Cascade docking further prioritized the hits and six active compounds were discovered. Ligand and receptor based pharmacophore screening in combination with molecular docking has been also used to discover CXC chemokine receptor 4 (CXCR4) antagonist [127] and p38α mitogen-activated protein kinase inhibitors [128]. Pharmacophore modeling has been also used as a post-filtering tool for molecular docking results [129], [130], [131]. In order to discover inhibitors of Met tyrosine kinase, Peach et al. [129] utilized a series of pharmacophore queries to filter out the docking poses that did not form critical interactions. Voet et al. [130] utilized pharmacophore modeling both as pre- and post-filtering tool for molecular docking. In a VS study to identify novel human androgen receptor (hAR) antagonist, they first created two pharmacophore queries, one for antagonist and the other for agonist, based on all hAR ligands co-crystal structures and homology modeling. The pharmacophore queries were then used to filter a library of commercially available small molecules to select compounds that fulfill the antagonist query and then to remove all those that match the agonist query. Molecular docking, pharmacophore post-filtering and visual selection were performed to select 31 compounds for evaluation of biological activity. The study resulted in the identification of two novel chemotypes with antagonist activities in low μM range while devoid of any agonist activity. As described above, the hierarchical combination of pharmacophore modeling and molecular docking has been extensively exploited in drug discovery. However, the performance of this HLVS scheme highly depends upon the accuracy of pharmacophore modeling. There are number of free and commercial pharmacophore modeling programs including Phase [132], [133], MOE-pharmacophore [93], LigandScout [134], Gasp [135], Catalyst [136], PharmID [137], etc. As different pharmacophore modeling programs vary in their performance, so the strength and weakness of each approach should be considered before applying them in HLVS protocols. The advantages and shortcomings of utilizing these pharmacophore modeling methods in drug discovery have been reported by Sanders et al. [138].

Molecular docking – molecular dynamics simulation – free energy

Recent advances in high performance computing has allowed the implementation of computationally demanding but more precise methods such flexible molecular docking and molecular dynamics simulation as VS tools. The success of a VS protocol depends on the accurate modeling of protein–ligand interactions, which require the consideration of protein and ligand in a dynamic environment. Several methodologies have been developed to account for protein flexibilities [139] that include soft docking [140], [141], sampling side-chain rotamer libraries and docking to an ensemble of receptor conformations [142], [143]. Molecular docking in combination with MD simulation has been especially recommended [144], [145]. MD simulation can be used to generate an ensemble of protein conformations for molecular docking [146], [147]. Moreover, MD simulation and free energy based ranking can be utilized in prioritizing molecular docking hits [148]. Both of these strategies have been effectively used in HLVS schemes for the identification of potent inhibitors. Wang et al. [149] carried out a docking-based VS using MD simulation to generate an ensemble of protein conformations. MD simulations optimized the docked complexes and led to the selection of compounds to be tested for human aldose reductase inhibition. Finally, in vitro assay identified two sub-micromolar compounds with novel chemotypes. Apart from generating receptor ensembles for molecular docking, MD simulation has also been used to identify functionally relevant receptor conformations, e.g. inactive DFG-out conformations in kinases to achieve selectivity against various kinases or the active state of GPCRs for the identification of agonists. These receptor conformations were then employed in corresponding HLVS campaigns [150], [151]. Hu et al. [152] reported F. tularensis FabI inhibitors identified from a HLVS protocol that integrated molecular docking with MD simulations. In their study, around 1000 hits from docking-based VS were subjected to a short MD simulation in implicit solvent. MM-PBSA binding free energy [53], [54] was used to rank-order compounds and 75 hits were singled out for biological assay that revealed 3 active hits. Using a similar approach, our group also successfully identified inhibitors of SUMO E1 protein [52]. In another study, Cao et al. [55] discovered novel tubulin inhibitors by taking advantage of HLVS protocol. They selected nine compounds for biological assay after screening an in-house library of around 100,000 compounds sequentially using molecular docking and MD simulation based MM-PBSA binding free energy calculation. Five out of 9 compounds demonstrated tubulin inhibitory activity. To efficiently take advantage of MM-PBSA binding free energy based ranking, computationally inexpensive VS methodologies have been developed [153], [154], [155]. Some of these studies reduce the time required by calculating MM-PBSA binding free energy from energy minimized single structure. These strategies were successfully integrated with molecular docking in HLVS schemes and resulted in the identification of potent inhibitors of plasmepsin II and Myt1 kinase [156], [157]. Although the integration of molecular docking with molecular dynamics simulation and binding free energy calculations seems to be an attractive option due to precise handling of protein–ligand-solvent interactions, it is still limited by some challenges. One obvious need is the requirement of huge computational resources to facilitate adequate sampling for accurate free energy calculations. Apart from computational requirements, the force fields used in molecular dynamics simulations are approximations and need further refinement [158]. Further, the calculation of binding free energy using relatively fast methods such as MM-PBSA rescoring [148], [153] relies on the ligand’s starting pose predicted by molecular docking and a wrong starting pose may lead to erroneous results. Rescoring of 10 docking poses has been suggested as a possible solution to achieve good accuracy [159].

HLVS for the identification of SMPPIIs

Recent advances in the discovery of small molecule protein–protein interaction inhibitors (SMPPIIs) have made the modulation of protein–protein interactions (PPIs) achievable [160], [161], [162], [163], [164]. VS approaches have played a significant role in changing the notion about the undruggability of PPIs. Several examples of VS studies have been reported in literature that targeted PPIs, including for HIV-1 IN–LEDGF/p75 [165], [166], [167], [168], XIAP–Caspase9 [169], Nrf2–Keap1 [64], [170], BCL2–BAK BH3 domain [171], uPAR–uPA [147], ΔF508–CFTR [172], RUNX1/ETO tetramerization [173], ERCC1–XPF [174], Annexin A2–S100A10 [175], SUMO–SIM [176], IFN-α–IFNAR [177], HIV Nef–SH3 [178] and p53–HDM2 [179]. Some of these VS studies took advantage of the hierarchical approach to identify initial hits that were later optimized into highly potent SMPPIIs [147], [165], [167], [168], [175], [176], [177], [178]. LBSB-HLVS was employed in some of these studies while other studies reported SB-HLVS approach. Substructure search in combination with molecular docking has been used to identify compounds disrupting the interaction of 14-3-3σ with aminopeptidase N [180]. Christ et al. [165] presented an excellent example of the employment of LBSB-HLVS for the identification of SMPPIIs for HIV-1 IN–LEDGF/p75 interaction. They have utilized pharmacophore-based screening, molecular docking and pharmacophore post-filtering to select 25 diverse molecules for biological assay. Optimization of initial hits first by testing analogs and then by medicinal chemistry led to the identification of compounds with potency in low μM range. Using a similar LBSB-HLVS approach comprising of pharmacophore screening and molecular docking, Reddy et al. [175] identified substituted 1,2,4-triazoles as inhibitors of the Annexin A2–S100A10 PPI. Voet et al. [176] translated the PPI between SUMO and SUMO interaction motif (SIM) into a pharmacophore query. This query was then used to screen a small molecule library and the resulting compounds were further prioritized using molecular docking and pharmacophore post-filtering. Multiple biological assays and HSQC-NMR revealed low μM inhibitors of SUMO-SIM PPI. Betzi et al. [178] applied drug-like filter, docking and pharmacophoric filter to identify inhibitors of HIV Nef–SH3 interaction with potency in μM range. Furthermore, novel SMPPIIs of Nrf2–Keap1 interaction were identified using HLVS scheme comprising of pharmacophore modeling and molecular docking. The most promising compound reported in this study disrupted the Nrf2–Keap1 interaction in vitro EC50 of 9.80 μM [170]. The Nrf2–Keap1 interaction was again targeted for SMPPIIs discovery, but this time utilizing an SB-HLVS approach [64]. Glide hierarchical docking protocol was used to identify nine hits with three novel scaffolds. In another study, Lavecchia et al. [57] built a theoretical model of frataxin/ubiquitin complex to reveal potential ubiquitin-binding site in frataxin. This interface was used to perform hierarchical docking-based VS that resulted in the identification of one compound that prevented frataxin ubiquitination and degradation. Subsequent scaffold searches identified a series of compounds able to disrupt frataxin/ubiquitin PPI with μM affinity. A similar hierarchical docking approach has been also used to identify inhibitors of PA–PB1 [181], p53–MDM2 [126], Cdc42–intersectin [182] and TWEAK–Fn14 [183] PPI. Using SB-HLVS, Khanna et al. [147] reported the discovery of SMPPIIs inhibiting the uPAR–uPA PPI. The compounds for biological assay were selected utilizing multi-step docking-based VS performed on a conformational ensemble of uPAR sampled from MD simulation. Extensive biochemical characterization revealed a compound inhibiting uPAR–uPA PPI with an IC50 of 10 μM. In an SB-HLVS study, Hu et al. [168] employed induced fit docking and MM-GBSA [184], [185] rescoring to identify 18 compounds targeting HIV-1 IN–LEDGF/p75 PPI. The two most potent compounds displayed IC50 values at 0.32 and 0.26 μM.

HLVS for the identification of GPCR modulators

G protein-coupled receptors (GPCRs) constitute the largest family of membrane proteins [186] and are the most important targets for the treatment of various diseases. GPCRs were the subject of many VS studies [94], [127], [187], [188], [189], [190], [191], [192], [193]. Several of these VS studies reported GPCR modulators following a hierarchical multi-step protocol where a large small molecule library was filtered to a smaller set to be tested by various assays. Due to the unavailability of crystal structures for most GPCRs, a majority of these HLVS were carried out using homology models [94], [151], [190], [192], [194]. Tikhonova et al. [94] used a structural model of free fatty acid receptor 1 (FFAR1) to perform VS employing 2D similarity search, 3D pharmacophore search and molecular docking. The biological assay of 70 compounds revealed 15 compounds acting as either antagonists or agonists of FFAR1. Engel et al. [190] created a structure-based pharmacophore that described interaction potential within the thyrotropin-releasing hormone (TRH) receptor binding pocket. ZINC small molecule library [47] was screened and the resulting compounds were prioritized for assay using flexible molecular docking and molecular descriptors. Five structurally diverse antagonists of the TRH receptors were identified. Recently using multiple crystal structure templates, Yoshikawa et al. [192] developed a homology model of CXCR7. The validated homology model was utilized for VS employing molecular docking and pharmacophore post-filtering. The study resulted in 21 CXCR7 ligands with novel scaffolds and IC50 ranging from 1.29 to 11.4 μM. In another recent report [151], dynamic pharmacophore modeling and molecular docking were used to discover agonists of serotonin receptor subtype 1A (5-HT1AR). Homology model of 5-HT1AR based on the active state of beta-2 adrenergic receptor was used for VS. Furthermore, de Graaf et al. [194] applied LBSB-HLVS to identify modulators of glucagon receptor (GLR) and glucagon-like peptide 1 receptor (GLP-1R). Both GLR and GLP-1R are challenging targets belonging to class B GPCRs. A library of 1.9 million drug-like compounds was filtered based on structural similarity with known antagonists. Subsequent hits were then docked to the trans-membrane cavity in a GLR homology model and finally selected for biological assay utilizing protein–ligand interaction information. Two of the 23 evaluated compounds displayed in vitro binding to GLR and inhibited glucagon effect. Interestingly, one compound inactive against GLR showed in vitro binding to GLP-1R. Lin et al. [195] combined homology modeling, molecular docking and MD simulations to establish kinase inhibitor Sorafenib as a nanomolar antagonist of 5-HT receptor. Other successful studies discovering GPCR modulators by performing HLVS virtual screening on homology models include reports for alpha1A adrenergic receptor [191], neurokinin-1 receptor [196] and nonpeptide CCR5 receptor [197]. Recent advances in the crystallography of GPCRs [198], [199] opened excellent opportunities for structure-based VS studies [200], [201], [202], [203], [204], [205], [206], [207]. A majority of these studies utilized only molecular docking as VS tool [201], [202], [203], [204], [206], [207]. Some of these studies applying HLVS approach on crystal structures include the discovery of histamine H1 receptor [200], β-2 adrenoreceptor [208], adenosine A2A receptor [208] and sphingosine 1-phosphate receptor [208] ligands. Among these studies, de Graaf et al. [200] developed and validated a HLVS method that combined molecular docking with interaction fingerprint (IFP) based scoring. The method was successfully applied to the identification of 19 human histamine H1 receptor ligands with potency ranging from 10 μM to 6 nM. In another study, Sanders et al. [208] combined ligand and structure-based methods to carry out prospective VS against three GPCRs (β-2 adrenoreceptor, adenosine A2A receptor and sphingosine 1-phosphate receptor). Biological assay revealed 6, 18 and 3 active ligands for β-2 adrenoreceptor, adenosine A2A receptor and sphingosine 1-phosphate receptor respectively. As described above, the majority of HLVS efforts were focused on the identification of small molecule modulators of class A or rhodopsin-like GPCRs. The reason may be the availability of plentiful ligand information and ever increasing crystal structures for this class of receptors. Barring a few studies [194], [209], [210], HLVS has not been exploited in the discovery of small molecule modulators for class B and C GPCRs, which may be due to the unavailability of sufficient structural information about receptor’s trans-membrane domain. However, the recent availability of trans-membrane domain crystal structures of class B [211], [212] and class C [213] GPCRs should augment drug discovery efforts. Specially, these developments should fuel efforts exploiting the power of LBSB-HLVS in the discovery of small molecular modulators of class B and C GPCRs.

Conclusion

We have witnessed tremendous growth in successful applications of VS in various drug discovery campaigns in recent years. VS methods relying on ligand and structure-based approaches have been developed and successfully employed in prospective applications. In order to enhance VS efficiency and improve hit rate, integration of several ligand and structure-based approaches has been suggested. A survey of recently published VS studies revealed hierarchical combination as the most common way of integrating ligand and structure-based methods. Mostly, computationally inexpensive ligand-based methods are used at the beginning of a multi-step hierarchical procedure while computationally demanding methods are used towards the end. Pharmacophore modeling and molecular docking are two extensively employed VS methods in HLVS schemes. The hierarchical combination of these two methods has been proven effective in the discovery of inhibitors of difficult drug targets such as PPIs. GPCRs are another class of important drug targets where HLVS has shown potential in retrieving promising hits, especially, considering the recent advances in GPCR crystallography.