Molecular Docking in Modern Drug Discovery: Principles and Recent Applications

Abstract

The process of hunt of a lead molecule is a long and a tedious process and one is often demoralized by the endless possibilities one has to search through. Fortunately, computational tools have come to the rescue and have undoubtedly played a pivotal role in rationalizing the path to drug discovery. Of all techniques, molecular docking has played a crucial role in computer aided drug design and has swiftly gained ranks to secure a valuable position in the modern scenario of structure-based drug design. In this chapter, the principle, sampling algorithms, scoring functions and diverse available software’s for molecular docking have been summarized. We demonstrate the interplay of docking, classical techniques of structure-based design and X-ray crystallography in the process of drug discovery. In addition, we dwell upon some of the limitations faced in docking studies. Finally, several success stories of molecular docking approaches in drug discovery have been highlighted, concluding with remarks on molecular docking for the future.

1. Introduction

The path to drug discovery is a long, challenging & arduous task not to mention the overburdening finances it demands. As of 2014, the average cost of developing a new drug from scratch was found to be an estimated $2.5 billion, an increase of 145% from the previous study done by the same organization in 2003. The major reasons for this drastic increase in the cost is mainly attributed to high failure rate of drugs among others [1]. Understanding of the drug discovery process is important to handle the challenges faced by the pharma companies in terms of cost and innovation.

The process of identifying a target, synthesizing an active compound with suitable characteristics like minimal toxicity, high bioavailability, cost-effective synthesis, etc., and finally developing it to introduce in the market is a time-consuming, extremely complex and risky endeavor [2]. Initially, a target is identified which plays a key role in progress of the disease. Once a link between the target and the disease has been established, the next step is to identify potential candidates which can stop or reverse the progress of the disease [3]. This process starts with the discovery of molecules that show efficacy in a simple screen, called “hits.” Screening is a process in which normally a large number of compounds from natural products and online databases are examined for biological activity in high-throughput assays. This step in the drug discovery process is very crucial and demands maintaining huge molecular libraries and carrying out thousands or millions of assays, which leaves the academicians and small pharmaceutical companies at a disadvantage and also shoots up the cost for larger industries. Next, the “hits” found are chemically modified to give improved pharmaceutical properties, such compounds are often called “leads.” But, it is quite apparent that the method stated above for discovery of a drug has a number of pitfalls. From an academic point of view, carrying out high throughput screens (HTS) is costly, time-consuming and not feasible; while, from an industrial perspective, it does nothing to improve the eminent danger of market saturation.

Truly innovative and blockbuster drugs are what drive the pharmaceutical industry forward but, over the past few years introduction of new molecular entities (NMEs) has vastly reduced. For example, in 2007 only 19 NMEs were approved by the US Food and Drug Administration (FDA), the least since 1983 [4]. Currently, and even in the future it is expected that only slight modifications of the existing blockbuster drugs would be carried out which would further aggravate this problem [5]. HTS would not help in either curbing the rising costs of discovering hits or the problem of finding truly innovative and blockbuster NMEs, the two major hurdles that the pharmaceutical industry faces now-a-days.

To overcome these challenges, molecular docking is an exemplary tool. During the first step to find hits from existing chemicals for a drug discovery and development project, virtual screening (VS) is a perfectly viable and an alternative or complementary approach to HTS for fulfilling the screening of thousands or millions of compounds within a few days. In addition, the speed of VS helps in kick-starting projects for newer targets for which no leads are available [6]. Molecular docking is one of the most applied virtual screening methods and has become increasingly useful overtime on account of immense growth in 3D X-ray and NMR structures and their improved resolution (physics and knowledge based docking algorithms depend on it) reported in the Protein Data Bank (PDB). As an example, in total 46,541 X-ray structures were reported at the end of 2008 in PDB, but by the end of 2018 it had grown to a staggering figure of 131,993 [7]. In addition, it is a resource saving technique which provides accessibility of screening to academia and small industries which were earlier limited to large pharmaceutical giants.

In this chapter, we will discuss a particular class of molecular design, i.e., “Docking” along with the various algorithms, techniques, success stories and limitations related to it. In the end, we will conclude with its scope in the near future.

2. Molecular docking

Two molecules can interact in a number of ways let alone the interaction of a protein and protein or a protein and small molecule. Molecular docking helps us in predicting the intermolecular framework formed between a protein and a small molecule or a protein and protein and suggest the binding modes responsible for inhibition of the protein. To accurately carry out docking studies one requires the high-resolution X-ray, NMR or homology-modeled structure with known/predicted binding site in the biomolecule. To date, 148,827 are available in the database (PDB) [3]. Docking methods fit a ligand into a binding site by combining and optimizing variables like steric, hydrophobic and electrostatic complementarity and also estimating the free energy of binding (scoring) [8].

There are two basic components which distinguish the variety of docking softwares available to choose from—One is sampling algorithm and the other is scoring function, these are discussed in detail.

2.1 Sampling algorithm

As pointed out earlier, there are a huge number of modes of binding between two molecules and even with advances in parallel computing and higher clock speeds of modern computers it would be expensive and time-consuming to generate all the possible modes. Therefore, algorithms were needed which could fish out valuable conformations from the fruitless ones.

Various algorithms were developed in this regard and can be classified by the number of degrees of freedom they ignore. The simplest of the algorithms introduced treated the molecules as two rigid bodies thereby reducing the degree of freedom to just six (three translational and three rotational). A very well cited example of a program using this algorithm is DOCK [9]. This program was designed to find molecules which had a huge extent of shape similarity to the pockets/grooves or binding sites. It derives an image of suspected binding sites present on the surface of the protein. This image consists of several overlapping spheres of varying radii which touch the molecular surface of the macromolecule at just two points. The ligand molecule is also considered as a set of spheres which approximately fill the space occupied by the ligand. Once the respective representations of the protein surface and the ligand in terms of sphere are complete, the pairing rule is applied. The pairing rule is based on the principle that ligand sphere can be paired with a protein sphere if the internal distances of all the spheres in the ligand set match all the internal distances within the protein set, allowing some user specified tolerance. Thus, it allows the program to identify geometrically similar cluster of spheres on the protein site and the ligand. Many other programs were developed later which make use of such matching algorithm (MA) which include LibDock [8], LIDAEUS [10], PhDOCK [11], Ph4DOCK [12], Q-fit [13], SANDOCK [14], etc. All these programs based on MA have the advantage of speed but have several limitations such as prior need for detailed receptor geometry and lack of molecular flexibility which does not accurately define many aspects of ligand-protein interactions.

The second algorithm is that of incremental construction (IC), wherein the ligand is fragmented from rotatable bonds into various segments. One of the segments is anchored to the receptor surface. The anchor is generally considered to be the fragment which shows maximum interactions with the receptor surface, has minimum number of alternate conformations and fairly rigid such as the ring system. Once the base/anchor has been established, the next step is to add each of the fragments step by step. Ideally, those fragments are added first which have a greater chance of showing interactions like hydrogen bonding since they are directional in nature and are responsible for specificity of the ligand. In addition, hydrogen bonds lead to more accurate prediction of geometry. Once a particular fragment is added, the poses with the least energies are considered for the next iteration, making the algorithm extremely fast and robust [15]. IC has been used in programs like DOCK 4.0 [16], FlexX [15], Hammerhead [17], SLIDE [18] and eHiTS [19], SKELGEN [20], ProPose [21], PatchDock [22], MacDock [23], FLOG [24], etc. One major limitation of this program is that it is restricted to medium sized ligands and is not feasible for large size ligands where the number of fragments generated pose a problem.

Another useful algorithm is the Monte Carlo (MC) technique. In this approach, a ligand is modified gradually using bond rotation and translation or rotation of the entire ligand. More than one parameter can also be changed at a time to get a particular conformation. That conformation is then evaluated at the binding site based on energy calculation using molecular mechanics and is then rejected or accepted for the next iteration based on Boltzmann’s probability constant. Acceptance or rejection of the conformation is a function of the change in energy with respect to a parameter T, which can be physically interpreted as temperature (simulated annealing). This criterion of acceptance or rejection makes this method strikingly different than the others. Whereas the other algorithm favor decrease in energy, in MC method increases are also possible. For higher values of T increases are likely. If one starts at a high value of T, then small energy barriers can be jumped and the configuration can move beyond local minima and is therefore particularly useful in situations where a global minimum is sought among many local minima [25]. An interesting spin-off of the MC approach is the Tabu search, which maintains a record of the search space of the binding site which has already been visited and thus ensures that the binding site is explored to the maximum [26]. MC approach has been made use of in programs like DockVision 1.0.3 [25], FDS [27], GlamDock [28], ICM [29], MCDOCK [30], PRODOCK [31], QXP [32], ROSETTALIGAND [33], RiboDock [34], Yucca [35], AutoDock [36], etc. One of the major concerns with MC approach is the uncertainty of convergence, which can be improved by performing multiple independent runs.

Genetic algorithm (GA) is quite similar to MC method and is basically used to find the global minima [37]. These are much inspired by the Darwin’s Theory of Evolution [38]. GA maintains a population of ligands with an associated fitness determined by the scoring function. Each ligand represents a potential hit. The GA alters the ligands of the population by mutation or crossover. In the first stage, a new population is created by accessing and then selecting the more fit ligands from the previous step. The members of the populations are then transformed in the alteration step. The mutation operator creates new ligands from a single ligand by randomly changing a fragment in its representation while the crossover operator exchanges information between two (occasionally more) members of the population [39, 40, 41]. GA has been incorporated in programs like Autodock 4.0 [42], DARWIN [43], DIVALI [39], FITTED [44], FLIPDock [45], GAMBLER [46], GAsDock [47], GOLD 3.1 [48], PSI-DOCK [49]. GA has a similar limitation like that of MC method wherein the uncertainty of convergence is a major drawback.

Another approach is the hierarchical method. In this approach, the low energy conformations of the ligand are pre-computed and aligned. The populations of the pre-generated ligand conformations are merged into a hierarchy such that similar conformations are positioned adjacent to each other within the hierarchy. Afterwards, on carrying out rotation or translation of the ligand, the docking program will make use of this hierarchical data structure and thus minimize the outcomes. Let us understand with a simple example—if an atom near the rigid center of the ligand is found to clash with the protein in a given rotation/translation, then this approach can reject all of the conformations lying below in the hierarchy to that of the conformation under scrutiny, because the descendants must contain the same clash as well [50]. GLIDE software makes use of the hierarchical method [51, 52].

2.2 Scoring functions

Sampling changes among varying degrees of freedom must be performed with sufficient accuracy to identify a conformation that best matches the receptor structure, and also must be fast enough to permit the evaluation of millions of compounds in a set computational time. This is taken care by the variety of algorithms discussed above. Algorithms are further complemented by scoring functions.

The evaluation and ranking of predicted ligand conformations is a crucial aspect of VS. When we are interested in only how a single ligand binds to a biomolecule, then the scoring function needs to predict the docked orientation which most accurately represents the “true” structure of the intermolecular complex. On the other hand, if we are interested to evaluate multiple ligands, in that scenario the scoring function should not only identify the “true” docking pose but also be able to rank one ligand relative to another. Therefore, the design of reliable scoring functions and schemes which can rank different poses is of fundamental importance [53].

The scoring functions usually estimate binding energy of complex using many assumptions and simplifications to arrive as close as possible to actual binding energy in minimum time. Popular scoring functions have an adequate balance between accurate estimation of binding energy and computational cost in terms of time. There have been a number of scoring functions developed over the past many years and can be classified into three main categories—force field, empirical and knowledge based [54].

Force field functions: force field (FF) scoring functions are developed based on physical atomic interactions like van der Waals interactions, electrostatic interactions and bond lengths, bond angles and torsions [55]. Force field functions and parameters are usually derived from both experimental data and ab initio quantum mechanical calculations according to the principles of physics.

$$E=\sum _i\sum _j\left(\frac{A_{\mathit{ij}}}{r_{\mathit{ij}}^{12}}-\frac{B_{\mathit{ij}}}{r_{\mathit{ij}}^6}+\frac{q_i{q}_j}{\epsilon \left(r_{\mathit{ij}}\right)r_{\mathit{ij}}}\right)$$E1

Here, r_ij stands for the distance between protein atom i and ligand atom j, A_ij and B_ij are the van der Waal parameters, q_i and q_j are the atomic charges and ε(r_ij) is the distance dependent dielectric constant.

One common example of a FF scoring function is that of the program DOCK [56] represented in Eq. (1), where, the effect of solvent is indirectly considered by the distance dependent dielectric constant e(r_ij) seen in the Coulombic potential. One major drawback of this function is that it does not consider an important solvent effect that charged groups favor aqueous environments whereas non-polar groups tend to stay in non-aqueous environments, commonly referred to as the desolvation effect [57]. Ignorance could lead to biased results as the function would now be totally dependent on the coulombic interactions and would thus favor highly charged ligands. In other words, it only takes into account the interaction of protein and ligand, which is inadequate. To build a more robust function one needs to also evaluate how both interact with water before the formation of the complex and how water mediates this process.

Later the Shoichet group [58] improved upon the existing function by adding the effects of the solvent on protein-ligand interactions using implicit solvent models. They employed the Poisson-Boltzmann approach to model the electrostatic potential of the protein. The van der Waals interactions were calculated using the Lennard-Jones potential; the electrostatic interaction between the ligand and the protein was estimated using a precomputed receptor potential determined with DelPhi [59]. Ligand desolvation penalties were calculated with HYDREN [60]. The solvent-corrected scores were found to be closer to experimental binding free energies than the DOCK program scores, but were still too favorable. The overestimation of complex stability could be due to the neglect of solute entropic terms [58].

There a few scoring functions which be classified in this category such as DockScore [56], GoldScore [61], HADDOCK Score [62], ICM SF [29], QXP SF [32], etc.

Empirical scoring functions: the basis of this scoring function is that the binding energies of a complex can be approximated by a sum of individual uncorrelated terms. The coefficients of the various terms involved in calculation of binding energy are obtained from regression analysis using experimentally determined binding energies or potentially from X-ray structural information. Empirical functions have simpler energy terms to evaluate when compared to force field scoring functions and thus are much faster in binding score calculations.

The first empirical scoring function developed to predict binding free energies was implemented in LUDI, credited to the pioneering work of Bohm [63]. The energy was derived using experimental binding free energies and protein-ligand crystal structures for 45 complexes.

$$\begin{array}{cc}∆G_{\mathit{\text{bind}}}& =∆G_O+∆G_{\mathit{hb}}\sum _{h-\mathit{\text{bonds}}}f\left(∆R,∆\alpha \right)+∆G_{\mathit{\text{ionic}}}\sum _{\mathit{\text{ionic}}\mathit{int}.}f\left(∆R,∆\alpha \right)+∆G_{\mathit{\text{lipo}}}\left|A_{\mathit{\text{lipo}}}\right|\\ & +∆G_{\mathit{rot}}\mathit{\text{NROT}}.\end{array}$$E2

Here, ∆G_o is the binding energy independent of protein interactions, ∆G_hb describes contribution to binding energy from hydrogen bonds, ∆G_ionic denotes contribution to binding energy from unperturbed ionic interactions, ∆G_lipo considers contribution to binding energy through lipophilic interactions while A_lipo is the lipophilic contact surface between the protein and the ligand, ∆G_rot describes the loss of binding energy due to freezing of internal degrees of freedom in the ligand while NROT represents number of rotatable bonds and f(∆R, ∆α) is a penalty function that accounts for large deviations from ideal hydrogen bond and salt bridge geometry.

As shown in Eq. (2), the binding free energy is modeled using hydrogen bonds, salt bridges, the hydrophobic effect, and solute entropy terms. The hydrogen bond and salt bridge terms are modified by a penalty function which accounts for deviation from ideal geometry. Entropy loss of the ligand upon complex formation is based on the Number of ROTatable bonds (NROT) in the ligand [64, 65]. Eldridge et al. presented an empirical scoring function referred to as ChemScore by taking into account different energetic parameters like hydrogen bonds, the lipophilic effects of atoms, the effective number of rotatable bonds in the ligand among others. The scoring function was calibrated using 82 ligand-receptor complexes with known binding affinities [66].

By including different empirical energy terms, many different empirical scoring functions have been developed such as SCORE2 [67], ChemScore [66], RankScore [68], LigScore [69], GlideScore [51], HINT [70], etc. The empirical scoring functions take into account many different energy terms and thus the problem of unknowingly double-counting of certain energy terms difficult issue to tackle.

Knowledge based scoring functions: these are derived from the structural information embedded in experimentally determined atomic structures. The functions use statistical analysis on crystal structures of complexes to obtain the interatomic contact frequencies between the protein and the ligand based on the presumption that stronger an interaction is, the greater the frequency of its occurrence will be. The overall score is calculated with the help of Eq. (3) by accounting for favorable contacts and repulsive interactions between each atom in the ligand and protein lying within a sphere with a specified cutoff [71, 72, 73, 74, 75, 76, 77, 78].

$$w\left(r\right)=-k_{B}Tln\left[g\left(r\right)\right],g\left(r\right)=\left(r\right)\rho \left(r\right)/{\rho }^{\ast }\left(r\right)$$E3

Here, k_B is the Boltzmann constant, T is the absolute temperature of the system, ρ(r) is the number density of the protein-ligand atom at distance r, ρ*(r) is the pair density in the reference state where interatomic interactions are zero and g(r) is pair distribution function.

Popular knowledge based functions include DrugScore [79], PMF [72, 80], MScore [81], SMoG [71], BLEEP [74], ITScore/SE [75], etc. The computational simplicity of such functions is a major advantage especially when one has large databases at hand however, the accuracy of predicting the reference state and underrepresentation of interactions with halogens and metals are the major hurdles.

Each of the above classified have their inherent drawbacks, in this regard, combination of more than one scoring functions has given improved results. This approach has been widely regarded as “Consensus Scoring” [46].

Another set of scoring functions which have recently started to attract attention are based on machine learning. One of the programs based on functions incorporating machine learning was able to achieve an astounding hit rate of 88.6% [82]. The nexus of machine learning and scoring functions is promising but the development of such a tool is slow owing to its complexity.

In order to compare the variety of scoring functions that have been developed up until now, comparative assessment of scoring functions (CASF) is an incredible platform to begin with [83].

There is another set of classification proposed for the scoring functions namely physics-based methods, empirical scoring functions, knowledge-based potentials, and descriptor-based scoring functions but there is still no clear consensus on which classification of scoring functions would be appropriate [84].

3. Applications

Molecular docking has been developed and improving for many years, but its ability to generate a viable drug is still generally questioned. In the section below, you will find examples where docking approach lead to recognition of active hits for a variety of different receptors/targets.

HIV 1 Integrase—a new binding site for drugs treating AIDS was discovered by Schames et al. using docking while considering the flexibility of the receptor through molecular dynamics. The group used AutoDock in conjunction with the relaxed-complex method to discover novel mode of inhibition of HIV integrase [85].

α1A Adrenergic receptor—Evers et al. generated a model of the receptor using homology modeling based on the X-ray crystallographic structure of bovine rhodopsin. Hierarchical virtual screening method was performed by them on the Aventis in-house compound repository in a stepwise manner. 22,950 filtered compounds were then docked into the α1A receptor homology model with the program GOLD and scored with PMF. The top scoring compounds were finally clustered according to their unity fingerprint similarity, and a diverse set of 80 compounds was tested in a radio ligand displacement assay. Thirty-seven compounds displayed a Ki < 10 μM with the most active having Ki = 1.4 nM [86].

Type I TGF-beta receptor kinase—A striking example and a proof of the benefit of in silico approach over classical high-throughput screening involves the discovery of novel Type I TGF-beta receptor kinase inhibitor. The same molecule (HTS-466284); Figure 1, a 27 nM inhibitor, was discovered independently using virtual screening [87] and also by traditional enzyme and cell-based high-throughput screening in the same year [88]. The compound discovered experimentally required in vitro screening of a large library of compounds in a TGF-β-dependent cell-based assay which required more time, proved to be costlier and required usage of a variety of chemicals when compared to its computational counterpart.

Figure 1.

Structure of HTS-46628, type I TGF-beta receptor kinase inhibitor.

Aurora Kinase A—A major improvement was seen in the inhibitory activity of Aurora Kinase A inhibitors which were designed using in silico techniques by Park et al. [89]. This research group made use of a genetic algorithm to carry out the sampling while the scoring function involved the energy terms from the AutoDock program with a slight modification of the dehydration energy term. The design strategy and tools used to carry out the study proved to be immensely successful with some inhibitors revealing exceptionally high potency at low picomolar levels; Figure 2 [89].

Figure 2.

Structures for Aurora Kinase A inhibitor with IC₅₀ 12 and 43 pM respectively.

Dopamine D3 receptor—The 3D structure of the Dopamine 3 (D3) subtype receptor was modeled by Varady et al. from the X-ray crystallographic structure of rhodopsin and validated using experimental data. A D3 pharmacophore model was devised by them from 10 selective and potent known D3 receptor ligands. Using their model, 250,251 compound were screened from the National Cancer Institute (NCI) 3D database. The hit list of 2478 potential ligands was then filtered for known chemotypes. After removal of all compounds that were structurally similar to known D3 receptor ligands, 1314 candidates remained. At the end, 20 compounds supplied by NCI to the group were tested, out of which eight had Ki values below 500 nM, among which one of the compounds had Ki = 11 nM; Figure 3 [90].

Figure 3.

Structure of dopamine D3 receptor inhibitor with Ki = 11 nM.

Serotonin receptor (5HT1A)—Due to lack of structural information available for the receptor, Becker et al. made use of PREDICT, to develop a unique non-homology model for building a virtual 3D structure of the receptor. Using the model, 40,000 compounds from Predix’s compound library were screened for molecular docking and 78 virtual hits were discovered and then purchased by them from respective vendors. The in vitro 5-HT1A binding assays elucidated that 16 of the 78 compounds tested by the group were found to be hits with Ki < 5 μM, reflecting a 21% hit rate, 9 of which had a Ki < 1 μM. The most potent molecule had Ki = 1 nM (Figure 4) and was selected as a lead molecule for further optimization. One significant feature of the study which highlights the utility of docking was that the complete discovery process, i.e., from in silico screening through lead optimization, preclinical, and into clinical studies, was very rapid, requiring less than a couple of years from program initiation to Phase I clinical trial [91].

Figure 4.

Structure of serotonin receptor inhibitor with Ki = 1 nM.

Crystal structure prediction challenge—The International Blind Test is a challenge organized by the Cambridge Crystallographic Data Center wherein a previously determined crystal structure is only revealed once all the participants submit their respective structures. In the Fifth International Blind Test, the challenge was toughened by including flexible molecules with 50–60 atoms. The successful prediction by two participants of the crystal structure of molecule XX in the blind test indicated that search methods and models for lattice energy are capable of providing worthwhile results, both in terms of the range of structures considered in the search and relative energies of the structures and thus can act as efficient ranking systems [92].

Muscarinic M3 receptor—A pharmacophore model was constructed by Marriot et al. from the known molecules showing significant M3 potency [93]. The research group utilized the program DISCO, which generated five models. Three models were rejected based on structural overlay. 3D screening was performed by Unity 3D of the Astra compound database. The first model developed by them gave 176 hits while the second model gave 173 hits; 172 compounds were common to the two sets and were tested for their M3-antagonistic potency. Several compounds with micromolar and even submicromolar activities resulted, for example, compound below had A50 M3 antagonism ≈ 0.2 μM; pA₂ = 6.67; Figure 5 [93].

Figure 5.

Structure of muscarinic M3 receptor antagonist.

Checkpoint Kinase 1—Lyne et al. utilized virtual screening to discover Checkpoint Kinase 1 (Chk-1) inhibitors [94]. Compounds with molecular weight > 600 or with more than 10 rotatable bonds were excluded from the database. Then 3D structures of the ligands were generated using Corina and a maximum of 8 stereoisomers were generated for each molecule. A 3D pharmacophore search was performed with their in-house program Plurality to eliminate compounds that do not have the typical binding motif for the kinase region. The remainder of the compounds were docked into the ATP binding site of Chk-1, using the program FlexX-Pharm, which considers full flexibility of the ligand but treats the protein as a rigid structure. The research group then utilized consensus scoring to identify molecules which were consistently giving good score with different scoring functions. Finally, visual inspection by the group of the 250 highest scoring hits for unfavorable interactions with the binding site or compounds with unrealistic conformations resulted in a list of 103 compounds for biological testing. Thirty-six hits were identified with IC₅₀ ranging from 110 nM to 68 μM; Figure 6 [94].

Figure 6.

Structures of checkpoint kinase 1 inhibitor with IC₅₀ 450 nM and 4 μM respectively.

Human Cathepsin K—Schröder et al. presented the implementation of a docking-based virtual screening workflow for the retrieval of covalent binders, human cathepsin K was utilized as a test case [95]. By using the filter of electrophilic war heads, a database with two million structurally diverse compounds with a variety of functional groups was reduced to a data set of just 343 test compounds. Molecular docking was performed by them and the top scoring poses of the GoldScore ranking list were taken into account for the manual selection of the virtual hits based on visual inspection of the appropriate fit of the molecule in the active site. A data set of 44 compounds including the five low scoring compounds were finally selected for experimental evaluation. The activity of 21 out of the selected 39 in silico hits was experimentally confirmed and four out of the five structures predicted as inactive showed no activity on cathepsin K. This study demonstrated to a huge extent the ability of docking to generate positive outcomes (Figure 7) [95].

Figure 7.

Respective structures for active and inactive covalent binders of human cathepsin K.

Human aldose reductase (ALR2)—ALR2 catalyzes a key reaction in the polyol pathway of glucose metabolism, a process implicated in the long-term complications of diabetes. Its inhibitors were designed by Wang et al. using molecular dynamic (MD) simulations and virtual screening [96]. A major challenge encountered by them in the in silico studies was that the binding site of the enzyme underwent large conformational changes and adopted distinct configurations upon binding different classes of ligands. To address this issue, the group sampled potentially accessible binding site conformations by MD simulations based on the available crystallographic structures of ALR2. After this procedure, three average conformations were selected for the docking. FlexX was utilized to carry out docking of 7200 compounds of which 128 compounds were selected by them for further screening. Out of these 72 molecules were selected which had RMSD < 3.00 A for experimental assay, of which 15 novel ALR2 inhibitors hits were discovered. The most potent inhibitor had an IC₅₀ = 0.24 μM; Figure 8 [96].

Figure 8.

Structure of human aldose reductase inhibitor with IC₅₀ = 0.24 μM.

Cyclooxygenase-2 (COX-2) and β-amyloid aggregation inhibitors—Dadashpour et al. made use of AutoDock4.2 to carry out docking studies of designed molecules based on diaryltriazine as lead. To validate the enzyme-inhibitor complex, the key molecular interactions and calculated binding energy were considered by them. Among the designed molecules, one of the compounds (Figure 9) showed an IC50 of 10.1 μM in experimental COX-2 assay. In addition, it showed potent anti-aggregation activity on β peptides [97].

Figure 9.

Structure of cyclooxygenase-2 inhibitor with IC₅₀ = 10.1 μM.

4. Limitations

The major limitation of molecular docking is due to the lack of confidence on the ability of scoring functions to give accurate binding energies. This stems from the fact that some intermolecular interaction terms are hardly predicted accurately, such as solvation effect and entropy change [98]. In addition, some intermolecular interactions are rarely considered in scoring functions which have been proven to be of significance. For instance, halogen bonding is verified to make a contribution to protein-ligand binding affinity [99] and so do guanidine-arginine interactions [100], but are not considered.

Transthyretin-thyroxine complex—One critical example wherein energy functions failed is that of transthyretin-thyroxine complex. The docking simulations with energy functions resulted in generation of two binding modes, one similar to the native binding mode of thyroxine and the other belonging to an alternate binding domain with a root mean square deviation (RMSD) of 8.97 Å from native binding state. The energy simulation was carried out and the lower energy solution picked by the docking program was the one with higher RMSD. Thus, in this case molecular docking failed to make the correct prediction of binding mode. Thereby, it would be fair to conclude that we might get many false negatives during the process of VS. [101].

It is still an unsolved problem to accurately deal with the water molecules in binding pocket during docking process, which is tough task and needs a lot of attention in the near future due to two reasons. Firstly, the x-ray crystal structures lack the coordinate information of hydrogen, due to inefficient scattering by smaller atoms. Not knowing the exact position of hydrogen leads to inaccuracies in identifying water molecules which might be acting as a bridging molecule between the ligand and the receptor. Secondly, no reliable theoretical approach is available to accurately predict how water molecules are affected by ligands and how strong the effect is. On top of that, it impossible with our current knowledge to predict how many water molecules in the binding pocket would be replaced by potential ligands and how the hydrogen bonding network would be disturbed by ligand binding [102].

One of the major challenges faced in the field of docking is that of rigid receptor. A protein can adopt many different conformations depending upon the ligand to which it binds. As a result, docking performed using a rigid receptor will correspond to a single receptor conformation, which leads to false negatives in many cases where later the ligand was found to be active. This happens because a protein can exist in constant motion between different conformational states having similar energies, which is usually neglected in docking [58].

Finally, the spectrum of activity against off-target proteins is something rarely seen even in computational screens and is only dealt by animal and human trials.

5. Conclusion

Thus, it is quite evident from the case studies highlighted above and many more success stories that one can find in literature related to computer aided drug design, that in silico approaches in combination with biophysical data, experimental high throughput screening and biology/toxicology/clinical studies are an indispensable tool in the process of drug discovery. It assists in decision making, conceptualizing new ideas and exploring them in a rapid manner to test a hypothesis, bringing solutions to problems that cannot be assessed experimentally either because the experiments is too difficult to design or because it would cost too much.

Undoubtedly, many challenges still remain to be addressed such as role of water molecules, solvent effects, entropic effects, and receptor flexibility.

There is more than sufficient information now that proves the utility of computational tools in drug design and there is no scope for any debate regarding the effectiveness and advantage of computational tools in the process of drug discovery.