Molecular targets are dynamic in their physiological environment, which are often crucial for various biological functions. The target binding pocket often adapts upon ligand binding to fit the ligands through various conformational changes ranging from small side-chain flip to large loop shift. Nevertheless, the experimentally solved target structures or ligand–target complex structures are basically static snapshots. Though previous works have shown that proper consideration of target flexibility can improve DBVS results (³³), it still represents one of the greatest challenges for current docking programs (³⁴) and becomes a hot issue in recent DBVS studies (³⁵–³⁹).

Ensemble docking that takes advantage of multiple target conformers has emerged as a partial solution to account for target flexibility in docking. The MultiCopyMD method developed by Okamoto et al. can generate a target ensemble through molecular dynamics (MD) with multiple ligands in the target binding site simultaneously (⁴⁰). Applying this target ensemble in their SBVS for novel inhibitors of death-associated protein kinase (DAPK), they discovered a highly potent (IC₅₀ = 69 nM) and selective inhibitor for DAPK1. To select appropriate target conformers, Rueda et al. have suggested a simple recipe by choosing the target conformers co-crystallized with the largest ligands (⁴¹), providing higher selectivity and better results than randomly picked ones when combined in ensemble. Using cyclin-dependent kinase 2 (CDK2) as a test example, Sperandio et al. have demonstrated normal mode analysis as an effective tool to select relevant target conformations with diverse binding sites (⁴²). Generally in ensemble docking, an individual docking run is required for each target conformation, which is thus computationally inefficient. To address this issue, Bottegoni et al. have proposed a 4D docking approach that allows fast and accurate account of target conformational ensembles in a single docking simulation (⁴³). This is achieved by merging 3D grids from optimally superimposed multiple target conformers into a single 4D object.

Some targets, such as metalloproteins, contain transition metal ions in their binding sites. The binding of ligand to these targets can be substantially distinct from other target types since such metal ions often coordinate ligand polar atoms, which may help to place and orient the ligand correctly in the binding sites. However, it is nontrivial to take metal ions into account accurately in current docking/scoring algorithms. The neglection of them would inevitably lead to underestimation of the metal–ligand interaction or even incorrectly docked ligands. Therefore, increasing attentions are being paid to metal ions in recent DBVS.

Röhrig et al. have studied the irons in heme proteins and demonstrated their importance for DBVS (⁴⁴). Two docking runs were performed in parallel by using a test set of 50 heme-containing complexes with iron–ligand contact. In one standard docking using EADock, a success rate of only 28% was achieved, clearly indicating the underestimation of the role of iron–ligand interactions. They then introduced the Morse-like metal binding potentials into EADock, which were fitted to reproduce density functional theory calculations. As a result, the success rate was doubled to 62%. To evaluate the reliability of the chosen docking protocol for screening potent cytochrome P450 aromatase inhibitors (AIs), Caporuscio et al. investigated a set of known imidazole and triazole AIs and found that the Glide docking program failed to predict a correct binding mode in all cases where the azole nitrogen coordinates the heme iron (⁴⁵). This observation inspired them to set up a metal constraint in Glide, which requires that a ligand atom lies within a certain region of the binding site in order to interact with specific target functionalities. Their structure-based design efforts eventually resulted in several novel AIs with IC₅₀ activity in the range of 21.7 μM to 9.4 nM.

Missing parameters of zinc ions is another common barrier for docking many metalloenzymes including histone deacetylases (HDACs). In seek of novel HDAC inhibitors, Park et al. derived potential parameters for zinc ions following a standard procedure (⁴⁶), in which geometry optimization of a simplified structural model was conducted for the active-site zinc ion cluster in complex with a hydroxamate-based inhibitor at the B3LYP/6–31 G** theory level. With these zinc parameters, they discovered six novel HDAC inhibitors with IC₅₀ value ranging from 1 to 100 μM.

There is a recognition that active-site water molecules play an important role in ligand-target binding (⁴⁷). Such water molecules can significantly contribute enthalpically and entropically to ligand–target binding. The most known role of water molecules is to mediate the ligand–target interaction by forming hydrogen bonds at the interface between the ligand and the target. On the other hand, the presence (or absence) and the location of water molecules may vary largely among ligands (⁴⁸). Despite their critical role, accounting for water molecules accurately in docking is a long-standing challenge. Several very recent studies directly targeted this issue.

Abel and coworkers have developed a unique approach WaterMap (⁴⁹) to account for the contribution of the displacement of water molecules by ligand to binding free energy. It first identifies “hydration sites” in the active site by clustering the trajectories from MD simulation of a solvated target with explicit water molecules. Inhomogeneous solvation theory is then applied to compute the thermodynamic properties of these active-site solvents including enthalpic and entropic changes. A displaced solvent functional is derived to estimate the relative binding free energies of a series of congeneric ligands based on their measured free energies by displacing active-site water molecules. This feature has made WaterMap particularly suitable for (and thus also limited to) lead optimization by providing insightful guidance to medicinal chemistry. More recently, WaterMap has been augmented by the introduction of an additional term attributable to the occupation of the dry regions in the target active site by ligand atoms (⁵⁰).

Lie et al. have proposed a very interesting approach that attached water molecules to ligand during docking (⁵¹). In their method, ligand polar atoms are solvated with maximum number of water molecules, which are then retained or displaced depending on energy contributions during docking simulation. The novelty of their method is that each water molecule is treated as a flexible on/off part of the ligand, instead of being a static part of the target. In such a manner, water molecules are sampled with the same flexibility as the ligand itself. Their method has been evaluated with considerable improvement by using 12 structurally diverse complexes, where several water molecules bridge the ligand and the target.

Rossato et al. have introduced a directional approach, AcquaAlta, to consider the solvation of ligand–target complexes (⁵²). Through an extensive analysis of the Cambridge Structural Database, they derived a geometric criteria defining interactions of water molecules with ligand and target. They also evaluated the propensity of ligand hydration through ab initio calculations. AcquaAlta has been validated with 20 crystal structures and reproduced 76% of the positions of water molecules that were experimentally observed.

Understanding of the interactions essential for ligand–target binding is critical to the success of lead discovery and optimization. For example, in a recent attempt to identify novel inhibitors of trihydroxynaphthalene reductase (3HNR) (⁵³), the authors first overlaid the known 3HNR inhibitors and then constructed a pharmacophore model that consists of several key interaction points within the active site: H-bonds with Ser149, Tyr163, Met200, and Tyr201 and π-stacking with Tyr208. In accordance to these interactions, the docking experiment was conducted in such a way that it only considered docking solutions that predicted π-stacking with Tyr208 and an optional H-bond with Ser149. The most potent hit compound they found exhibited a K_i of 5.3 ± 0.3 μM against 3HNR.

As revealed by the crystal structures of kinases in complex with ATP-competing inhibitors, such inhibitors typically form at least one hydrogen bond with backbone amide or carbonyl groups in the hinge region. Therefore, introducing relevant constraints with the hinge region for the molecules docked into the ATP sites of kinases would improve the chance of finding active compounds. This has been practiced by Ravindranathan et al. in the hit discovery of fibroblast growth factor receptor 1 (FGFR1) (⁵⁴). Among the 23 purchasable compounds suggested by a virtual screening experiment against 2.2 million compounds, two were identified to inhibit FGFR1 kinase with medium potency (IC₅₀ = 23 and 50 μM, respectively).

For certain target or ligand system, specifically designed methods may be more efficient. For example, Lang et al. recently have optimized DOCK 6 for docking small molecules to RNA targets (⁵⁵) and obtained a success rate of 70% for the ligands with less than seven rotatable bonds at the 2-Å heavy-atom root-mean-squared deviation threshold. The BALLDock/SLICK developed by Kerzmann is a ligand-specific docking approach for docking carbohydrate or carbohydrate-like compounds, which are often problematic for standard docking programs (⁵⁶).

The most straightforward way to obtain a target-biased scoring function is, probably, to re-calibrate an existing all-purpose scoring function directly on certain target classes. For example, DrugScore-RNA (⁸²) adopts the same framework as DrugScore (⁶⁹) but is derived from 670 crystal structures of nucleic acid–ligand and nucleic acid–protein complexes. Similar idea has been implemented in the kinase family-specific potential of mean force (kinase-PMF) (⁶⁸), a kinase-targeted scoring function adjusted from the original PMF04 (⁶⁷).

Tweaking the parameters in original scoring functions toward specific targets is also a prevalent strategy to derive target-biased scoring functions. For example, Teramoto and Fukunishi have applied a supervised scoring model to tailor the FlexX scoring function (F-score), which outperformed its former version on three of the five tested targets (⁸³). The TOP approach suggested by Seifert (⁸⁴) have employed iterative taboo search to optimize the scoring function in ProPose and the original Böhm scoring function against three targets, including CDK2, estrogen receptor, and COX2. By adding negative data of ligands that are known not to bind particular target, Pham and Jain have tuned the scoring function in Surflex-Dock and observed substantially enhanced screening enrichment for HIV protease and poly(ADP-ribose) polymerase (⁸⁵). An augmented Flo+ scoring function has been developed by Catana and Stouten using N-way partial least squares (PLS) (⁸⁶), which significantly improved the correlation between observed and calculated pK_i values from R² = 0.5 to 0.8 on a relatively diverse set of ligand–target complexes spanning seven protein families. Therefore, it would be attractive if scoring functions offer extendable or customizable features.

The above-mentioned target-biased scoring functions typically require re-parameterization or special treatment of established scoring functions. Too often, existing scoring functions are available to end-users as black boxes, hence it is not readily possible to adjust their parameters by any optimization algorithm. Several approaches have been proposed to address this issue. One of the earliest examples is the MultiScore that employs the raw scores from eight scoring functions to characterize the observed pK_i (⁸⁷), which has been found to work better for matrix metalloproteinases. The implied idea is slightly different from that of consensus scoring (⁷⁹,⁸⁰) in that it assumes uneven contributions from individual scoring functions. In a similar way, the AutoShim method has incorporated the original Flo+ score as well as additional target-specific pharmacophore points (shims) as descriptors in PLS analysis (⁸⁸). More recently, Cheng et al. have proposed a knowledge-guided strategy (KGS) based on the similarity principle aiming to improve the accuracy of binding affinity prediction of current scoring functions (⁸⁹). The KGS strategy computes the binding affinity of a query ligand–target complex based on the known binding affinity of an appropriate reference complex, which is required to share a similar pattern of key ligand–target interactions to that of the query complex of interest. The KGS strategy has been validated with both observed and docked ligand–target complex structures. Moreover, it can in principle work in concert with any scoring method, and its application is not limited to specific classes of ligand–target complexes.