Advances and Challenges in Protein-Ligand Docking

1Molecular recognitions including enzyme-substrate, drug-protein, drug-nucleic acid, protein-nucleic acid, and protein-protein interactions play important roles in many biological processes such as signal transduction, cell regulation, and other macromolecular assemblies. Therefore, determination of the binding mode and affinity between the constituent molecules in molecular recognition is crucial to understanding the interaction mechanisms and to designing therapeutic interventions. Due to the difficulties and economic cost of the experimental methods for determining the structures of complexes, computational methods such as molecular docking are desired for predicting putative binding modes and affinities. In molecular docking, based on the protein structures, thousands of possible poses of association are tried and evaluated; the pose with the lowest energy score is predicted as the “best match”, i.e., the binding mode. Since Kuntz and colleagues’ pioneering work, significant progress has been made in docking research to improve the computational speed and accuracy. Among them, protein-ligand docking is a particularly vibrant research area because of its importance to structure-based drug design and will be the subject of the present review.

2A protein-ligand docking program consists of two essential components, sampling and scoring. Sampling refers to the generation of putative ligand binding orientations/conformations near a binding site of a protein and can be further divided into two aspects, ligand sampling and protein flexibility. Scoring is the prediction of the binding tightness for individual ligand orientations/conformations with a physical or empirical energy function. The top orientation/conformation, namely the one with the lowest energy score, is predicted as the binding mode. Here, we will review the recent advances in protein-ligand docking on three important aspects: protein flexibility, ligand sampling, and scoring function, as illustrated in Figure 1. Challenges and future directions will also be discussed.

3Ligand binding commonly induces protein conformational changes (referred to as “induced fit”), which range from local rearrangements of side-chains to large domain motions. Due to the large size and many degrees of freedom of proteins, their flexibility may be the most challenging issue in molecular docking. Current methods to account for protein flexibility can be grouped into four categories: soft docking, side-chain flexibility, molecular relaxation, and protein ensemble docking.

4Soft docking is the simplest method which considers protein flexibility implicitly. It works by allowing for a small degree of overlap between the ligand and the protein through softening the interatomic van der Waals interactions in docking calculations. The advantages of soft docking are its computational efficiency and easiness for implementation. However soft docking can account for only small conformational changes.

5Many of the early attempts to incorporate certain protein conformational changes into molecular docking focused on side-chain flexibility, in which backbones are kept fixed and side-chain conformations are sampled. One of the earliest studies is the ligand docking algorithm developed by Leach, in which discrete side-chain flexibility is included by using a rotamer library. Since then, researchers have proposed many improved techniques to incorporate continuous or discrete side-chain flexibility in ligand docking.

6The third type of methods account for protein flexibility by firstly using rigid-body docking to place the ligand into the binding site and then relaxing the protein backbone and side-chain atoms nearby. Specifically, the initial rigid-body docking allows for atomic clashes between the protein and the placed ligand orientations/conformations in order to consider the protein conformational changes. Then, the formed complexes are relaxed or minimized by Monte Carlo (MC), Molecular Dynamic simulations, or other methods. The advantage of the molecular relaxation method is the inclusion of certain backbone flexibility in addition to the side-chain conformational changes. However, compared to the side-chain flexibility methods, the relaxation method is more demanding for the scoring function because it involves not only the side-chain movement but also the more challenging backbone sampling, thereby inaccuracies in the scoring function may lead to artifacts (e.g., improper backbone torsions) in the relaxed protein conformations. Moreover, the relaxation method is time-consuming.

7The most widely-used type of methods for incorporating protein flexibility utilize an ensemble of protein structures to represent different possible conformational changes. One of the earliest studies was done by Knegtel et al., in which an averaged energy grid was constructed by combining the energy grids generated from individual experimentally-determined protein structures using a weighting scheme, followed by standard ligand docking. Osterberg et al. extended the method to AutoDock with a larger ensemble consisting of 21 different conformations of the HIV-1 protease. The averaging nature of the method may miss the geometric accuracy of the protein. Claussen et al. developed a docking program FlexE to dock ligands into an ensemble of protein structures, in which the similar segments of the protein structures in the ensemble are aligned and merged while the dissimilar segments are used to combinatorially create new possible protein conformations for docking. In Wei et al.’s algorithm, a protein was decomposed into a rigid part and several flexible parts according to the crystal structures of the protein in the ensemble. For a given ligand placement, only the best-fit local conformer was kept for each flexible part in the protein, assuming the flexible regions move independently. The selected local conformers were joined with the rigid part to form the best-fit protein conformation. Compared to FlexE, this algorithm is significantly faster and scales linearly rather than exponentially with the protein flexibility. However, the algorithm may rely on the quality of ligand orientational/conformational sampling, namely whether proper initial ligand placement is included. Huang and Zou developed a fast ensemble docking algorithm by treating the protein conformational ensemble as an additional dimension to the traditional six degrees of freedom (three translational plus three rotational) for ligand energy optimization. The algorithm is almost as fast as single docking whereas keeping the accuracy of sequential docking. The ensemble docking algorithm is not used for generating new protein structures, but instead for selecting the induced-fit structure from a given protein ensemble. Following a similar procedure, Abagyan and colleagues expanded Huang and Zou’s algorithm to create ICM’s ensemble docking algorithm, referred to as four-dimensional (4D) docking. In addition to experimental structures such as NMR structures or crystal structures bound of the protein with different ligands, ensembles of protein conformations generated by molecular dynamic simulations, Monte Carlo simulations, or structure prediction have also been used to account for protein flexibility.

8Ligand sampling is the most basic element in protein-ligand docking. Given a protein target, the sampling algorithm generates putative ligand orientations/conformations (i.e., poses) around the chosen binding site of the protein. The binding site can be the experimentally determined active site, a dimer interface or other site of interest. Ligand sampling is the most successful area being developed in protein-ligand docking. Roughly, there are three types of ligand sampling algorithms: shape matching, systematic search, and stochastic algorithms.

9The shape matching method is one of the simplest sampling algorithms which is often used in the early stages of the docking process or in the first step of other more advanced ligand sampling methods. It places the ligand using the criterion that the molecular surface of the placed ligand must complement the molecular surface of the binding site on the protein. The six degrees of freedom (three translational and three rotational) of the ligand allow for many possible ligand binding orientations. Therefore, how to quickly place the ligand in the binding site with a good shape complementarity is the goal of the shape matching algorithm. Examples of docking programs that use shape matching include DOCK, FRED, FLOG, EUDOC, LigandFit, Surflex, MS-DOCK, and MDock. The major advantage of shape matching is its computational efficiency. However, the conformation of the ligand is normally fixed during the shape matching process. Therefore, flexible-ligand docking with the shape matching method is usually performed by docking an ensemble of pre-generated ligand conformations into the protein, followed by merging/reranking of the docked poses from different docking runs according to their energy scores (see below).

10Systematic search algorithms are normally used for flexible-ligand docking, which generate all possible ligand binding conformations by exploring all degrees of freedom of the ligand. There are three types of systematic search methods: exhaustive search, fragmentation and conformational ensemble.

11The most straightforward systematic algorithms are exhaustive search methods, in which flexible-ligand docking is performed by systematically rotating all possible rotatable bonds of the ligand at a given interval. Despite its sampling completeness for ligand conformations, the number of the combinations can be huge with the increase of the rotatable bonds. Therefore, to make the docking process practical, geometric/chemical constrains are normally applied to the initial screening of ligand poses, and the filtered ligand conformations are further subject to the more accurate refinement/optimization procedures. Glide and FRED are two typical examples of this type of hierarchical sampling methods.

12In fragmentation methods, the ligand is first divided into different rigid parts/fragments. Then, the ligand binding conformation is incrementally grown by placing one fragment at a time in the binding site or by docking all the fragments into the binding site and linking them covalently. DOCK, LUDI, FlexX, ADAM, and eHiTs are example docking programs that use fragmentation methods.

13In conformational ensemble methods, ligand flexibility is represented by rigidly docking an ensemble of pre-generated ligand conformations with other programs such as OMEGA (OpenEye Scientific Inc, NM). Then, ligand binding modes from different docking runs are collected and ranked according to their binding energy scores. Examples of the conformational ensemble methods for docking include FLOG, DOCK3.5, PhDOCK, MS-DOCK, MDock, and Q-Dock.

14In stochastic algorithms, ligand binding orientations and conformations are sampled by making random changes to the ligand at each step in both the conformational space and the translational/rotational space of the ligand, respectively. The random change will be accepted or rejected according to a probabilistic criterion. There are four types of stochastic algorithms: Monte Carlo (MC) methods, evolutionary algorithms (EA), Tabu search methods, and swarm optimization (SO) methods.

15In a Monte Carlo method, the probability to accept a random change is calculated by using the following Boltzmann probability function:

16where E0 and E1 stand for the energy scores of the ligand before and after the random change, respectively, kB is the Boltzmann constant, and T is the absolute temperature of the system. The docking programs that use the MC methods include DockVision, ICM, QXP, Prodock, and MCDOCK.

17Evolutionary algorithms (EAs) search for the correct ligand binding mode based on the idea from the evolutionary process in biological systems. The most popular type of EAs is the genetic algorithms (GAs). GOLD, AutoDock, DIVALI, DARWIN, MolDock, PSI-DOCK, FLIPDock, Lead finder, and EADock are the examples that have implemented evolution algorithms.

18In Tabu search methods, the probability of acceptance depends on the previously explored areas in the conformational space of the ligand. The random change will be rejected if the RMSD between the current ligand binding conformation and any of the previously recorded solutions is less than a cutoff; otherwise, the random change will be accepted. Example docking programs are PRO_LEADS and PSI-DOCK.

19Swarm optimization (SO) algorithms attempt to find an optimal solution in a search space by modeling swarm intelligence. In the method, movements of a ligand mode through the search space are guided by the information of the best positions of its neighbors. Examples of docking programs that use swarm optimization algorithms include SODOCK, Tribe-PSO, PSO@Autodock, and PLANTS.

20The scoring function is a key element of a protein-ligand docking algorithm, because it directly determines the accuracy of the algorithm. Speed and accuracy are the two important aspects of a scoring function. An ideal scoring function would be both computationally efficient and reliable. Numerous scoring functions have been developed in the past decades and can be grouped into three basic categories according to their methods of derivation: force field, empirical, and knowledge-based scoring functions.

21Force field (FF) scoring functions are based on decomposition of the ligand binding energy into individual interaction terms such as van der Waals (VDW) energies, electrostatic energies, bond stretching/bending/torsional energies, etc., using a set of derived force-field parameters such as AMBER or CHARMM force fields. One of the major challenges in FF scoring functions is how to account for the solvent effect. The simplest method is to use a distance-dependent dielectric constant ɛ (rij) such as the force field scoring function in DOCK:

22where rij stands for the distance between protein atom i and ligand atom j, Aij and Bij are the VDW parameters, and qi and qj are the atomic charges. ɛ (rij) is usually set to 4rij, reflecting the screening effect of water on electrostatic interactions.

23The most rigorous FF methods are to treat water molecules explicitly such as FEP and TI (see for review). However, these methods, together with their simplified approaches such as LIEPROFEC, and OWFEG are computationally expensive. To reduce the computational expense, accelerated methods have been developed while preserving the reasonable accuracy by treating water as a continuum dielectric medium. The Poisson-Boltzmann/surface area (PB/SA) models and the generalized-Born/surface area (GB/SA) models are typical examples of such implicit solvent models.

24In addition to the challenge on solvent effect, how to accurately account for entropic effect is an even more severe challenge for FF scoring functions. Moreover, whether the free energy of ligand binding can be decomposed into a linear combination of individual interaction terms without calculating the partition function (“ensemble average”) also remains in question, referred to as the nonadditive problem.

25In empirical scoring functions, the binding energy score of a complex is calculated by summing up a set of weighted empirical energy terms such as VDW energy, electrostatic energy, hydrogen bonding energy, desolvation term, entropy term, hydrophobicity term, etc.

26where {ΔGi} represent individual empirical energy terms, and the corresponding coefficients {Wi} are determined by reproducing the binding affinity data of a training set of protein-ligand complexes with known three-dimensional structures, using least squares fitting. Compared to the force field scoring functions, the empirical scoring functions are normally much more computationally efficient due to their simple energy terms. However, the general applicability of an empirical scoring function depends on the training set due to the nature of its fitting to known binding affinities of its training set. With the rapid increase in the number of crystal structures of diverse protein-ligand complexes with known binding affinities, a general empirical scoring function could be developed by training on the binding constants of thousands of protein-ligand complexes. GlideScore, PLP, SYBYL/F-Score, LigScore, LUDI, SCORE, X-Score, ChemScore, MedusaScore, AIScore, and SFCscore are examples of empirical scoring functions.

27The potential parameters of knowledge-based scoring functions are directly derived from the structural information in experimentally determined protein-ligand complexes. The principle behind knowledge-based scoring functions is the potential of mean force, which is defined by the inverse Boltzmann relation

28where kB is the Boltzmann constant, T is the absolute temperature of the system, ρ(r) is the number density of the protein-ligand atom pair at distance r in the training set, and ρ*(r) is the pair density in a reference state where the interatomic interactions are zero. After the potential parameters w(r) are derived, the energy of ligand binding for a given complex is simply the sum of the interaction terms for all the protein-ligand atom pairs in the complex. Based on the early idea of Tanaka and Scheraga, a number of knowledge-based scoring functions have been developed for protein-ligand interactions. Compared to the force field and empirical scoring functions, the knowledge-based scoring functions offer a good balance between accuracy and speed. Namely, because the potentials in Equation. (4) are extracted from a large number of structures rather than attempting to reproduce the known affinities by fitting, the knowledge-based scoring functions are relatively robust and general. Their pairwise characteristic also enables the scoring process to be as fast as empirical scoring functions.

29As the ideal reference state is inaccessible for complicated systems like proteins, one major challenge for knowledge-based scoring functions is the calculation for the afore-mentioned reference state; based on the methods knowledge-based scoring functions can be classified into three categories: traditional atom-randomized reference state, corrected reference state, and circumventing the reference state.

30Traditional methods to approximate the reference state are randomization of the atoms in the training set. Examples include DrugScore, SMoG, BLEEP, GOLD/ASP, MScore, and KScore. The disadvantage of the atom-randomization approximation is the neglection of the effects of excluded volume, interatomic connectivity, etc.. Methods like PMF and DFIRE have introduced correction terms for the reference state. Yet, binding mode prediction and virtual database screening are main problems for most knowledge-based scoring functions as a result of the reference state problem. To circumvent the long-standing reference state problem, Huang and Zou have developed a physics-based iterative method and derived the ITScore scoring function, which has been extensively tested with multiple benchmarks for binding mode prediction, affinity prediction and virtual screening.

31Other challenges for knowledge-based scoring functions include extension of the pairwise interactions to many-body interactions to account for hydrogen bonding and other directional interactions, development of an accurate method for entropic calculations, etc.

32Consensus scoring is not really a specific type of scoring function but a technique in protein-ligand docking. It improves the probability of finding a correct solution by combining the scoring information from multiple scoring functions in hopes of balancing out the errors of the individual scoring functions. Therefore, the main issue in consensus scoring is how to make the combination rule for individual scores so that the true binders can be discriminated from others according to the consensus rule. MultiScore and X-Cscore are two examples of consensus scoring methods.

33In addition to consensus scoring, another technique to improve the performances of scoring functions is clustering-based scoring methods, which incorporate the entropic effects by dividing generated ligand binding modes into different clusters. The entropic contribution in each cluster is measured by the configurational space covered by the ligand poses or the number of the ligand poses in the cluster. One restriction in clustering-based scoring methods is that its performance depends on the ligand sampling protocol that is used, i.e., it is docking program-dependent. These methods in combination with ligand conformational sampling using AutoDock have significantly improved binding mode prediction.

34We have reviewed three important aspects of protein-ligand docking: protein flexibility, ligand sampling, and scoring functions. Rapid advances in the last two decades have almost solved the ligand sampling issue. Although equal or even more efforts have been paid to scoring function development, entropy and desolvation effects remain the two major challenging issues for current scoring functions, particularly for the force field scoring functions. Speed and accuracy are the two important characteristics of a scoring function. Because of the rapid increase in computing power, how to improve the accuracy is the future direction for scoring function development. In contrast to ligand sampling and scoring functions which have been extensively studied for more than two decades, protein flexibility has only been addressed recently because of the difficulty resulting from the enormous degrees of freedom and the limitation of the computing power. The development of computational methods for protein flexibility is still in its infancy and thereby remains one of the major future directions in protein-ligand docking. Finally, how to evaluate different docking programs and scoring functions is another active area. Although many comparison studies for docking and scoring have been published, publicly available docking benchmarks such as DUD (http://dud.docking.org/) and CSAR (http://www.csardock.org/) are extremely valuable for systematic and consistent evaluation and improvement of new and existing docking algorithms.