2.1 Conformational variability

Our understanding of receptor–ligand binding has significantly advanced from the original lock-and-key model proposed by Fischer (Reference Fischer1890). Early experimental studies showed that the act of ligand binding influences the protein conformation, referred to as conformational induction or induced fit (Koshland, Reference Koshland1958). Another model of ligand binding is conformational selection, wherein the ligand chooses a binding partner from among available states in the conformational ensemble, thereby shifting the population distribution (Berger et al. Reference Berger, Weber-Bornhauser, Eggenberger, Hanes, Pluckthun and Bosshard1999; Foote & Milstein, Reference Foote and Milstein1994; Kumar et al. Reference Kumar, Ma, Tsai, Sinha and Nussinov2000; Monod et al. Reference Monod, Wyman and Changeux1965; Tsai et al. Reference Tsai, Kumar, Ma and Nussinov1999a, Reference Tsai, Ma and Nussinovb, Reference Tsai, Ma, Sham, Kumar and Nussinov2001). Recently, several papers have been published examining the evidence for whether receptor binding occurs because of induced fit or conformational selection. Sullivan & Holyoak (Reference Sullivan and Holyoak2008) studied the kinetics of phosphoenolpyruvate carboxykinase and concluded that induced fit and conformational selection were not mutually exclusive, rather they were complementary avenues for binding. Weikl & von Deuster (Reference Weikl and von Deuster2009) developed equations of binding kinetics using four states, the dominant apo and holo protein conformations without a bound ligand and those same conformations with a ligand bound, to describe the protein–ligand complex. They used the equilibrium constants from this model to distinguish between the induced and selected fit. Hammes et al. (Reference Hammes, Chang and Oas2009) examined the binding pathways for dihydrofolate reductase (DHFR) and flavodoxin and determined that a mixed binding mechanism was most likely. They noted that the relative importance of induced fit or conformational selection for a particular case could be analyzed by comparing the flux through the reaction path of ligand-binding and conformational change. Both mechanisms of binding will produce the same result; it is important only that some mechanism of receptor conformational change be incorporated in docking simulations. This is particularly useful in SBDD because it implies that computationally inexpensive methods that include protein flexibility should correctly predict binding modes.

2.2 Allostery

Protein–ligand flexibility upon binding is crucial for proper understanding of allosteric regulation. Nussinov and co-workers postulated that most proteins exist in an ensemble of states, and thus most proteins have the potential for allostery (Gunasekaran et al. Reference Gunasekaran, Ma and Nussinov2004). Based on the experimental literature, they demonstrated that the binding of an allosteric ligand shifted the population of conformational substates, thereby influencing the ability of other ligands to bind an alternate site.

Drug resistance is a growing problem that calls for new approaches to drug therapy. Exploring allosteric control in protein targets provides exciting opportunities for finding new modes of action and hence, overcoming emergent resistance, and developing cocktails of drugs to improve treatment.

3.1 The cross-docking problem

Research into protein flexibility and allostery has lent support to the importance of representing multiple states in binding studies (Fig. 1). Mobley & Dill (Reference Mobley and Dill2009) noted that binding free energy (ΔG _bind) and entropy are influenced by the shape and width of the entire conformational landscape, rather than a single rigid pose. Murray et al. (Reference Murray, Baxter and Frenkel1999) examined approaches for using rigid receptors in docking studies. When the authors attempted to dock a known ligand into a protein structure solved in the presence of a different ligand (referred to as cross-docking), they found that the active site was biased toward to the native ligand (Fig. 2). A variety of differences in the surface of the binding site were identified for the same protein solved with different ligands or in the absence of a ligand. Movement was observed in the backbone, side chain (both dependent and independent of the backbone Cα), and active site metals. As a consequence, active sites were biased toward a particular ligand type, negatively impacting docking efforts. The resultant misdocking could not be overcome without accounting for critical conformational shifts. Najmanovich et al. (Reference Najmanovich, Kuttner, Sobolev and Edelman2000) examined the percentage of residues that actually undergo a change upon ligand binding, based on two non-redundant datasets containing paired holo- and apo-protein structures from the PDB. No significant correlation was found between backbone movement and side-chain flexibility, but Lys, Arg, Gln, and Met were identified as the residues most likely to demonstrate side-chain rearrangement upon ligand binding. Najmanovich et al. (Reference Najmanovich, Kuttner, Sobolev and Edelman2000) showed that 60–70% of the binding site undergoes some changes in the side-chain orientation. Taken together, these studies clearly illustrated that the lack of protein flexibility in typical SBDD efforts severely limits the identification of true ligands and accurate docked poses.

Fig. 1.

The conformational variation in holo- and apo structures of beta-secretase 1 illustrates structural differences frequently seen across multiple crystal structures of the same protein.

Fig. 2.

The same protein crystallized with different ligands. This simple two ligand-one protein set demonstrates the influence of structure on the shape of the binding pocket.

An underlying assumption in rigid docking efforts is that the complexed state is the lowest free energy state. However, the ligand-bound state is not always the lowest energy state for either the protein or the ligand. Smith et al. (Reference Smith, Redfield, Smith, Dobson, Clore, Gronenborn, Walter, Naganbushan and Wlodawer1994) were first to examine conformational similarity, they compared the four available structures of interleukin-4, two from NMR and two from crystallography. They found that differences in the experimental methods clearly affected the structure, and that the two NMR structures were more similar to the crystal structures than to each other. Further, they noted that the exposed regions of the surface, particularly the highly flexible loops, showed the largest conformational difference among the structures. Eyal et al. (Reference Eyal, Gerzon, Potapov, Edelman and Sobolev2005) evaluated a set of 659 structure pairs from the PDB and showed that experimental variation in protein structures certainly exists, as a result of refinement, crystal packing, and/or crystallization conditions. Chopra et al. (Reference Chopra, Summa and Levitt2008) employed Molecular Dynamics (MD) on a set of 75 proteins to explore the role of solvent in structural determination and concluded that solvent effects have an enormous influence on the final refined protein structure.

3.2 Sources of structural information

Use of a set of similar conformations only generates a finite amount of information on potential binding partners. Additional receptor space should be explored, especially in scenarios where alternate binding modes are of interest. X-ray crystallography, NMR structures, homology models, normal mode analysis (NMA), and molecular mechanics simulations are all potential sources of structural diversity. Advances in technology have allowed for significant increases in available structural information for a vast number of receptor targets. The PDB now contains over 76 000 structures (as of 9/27/11), up from 13 605 structures in 2000 (Berman et al. Reference Berman, Bhat, Bourne, Feng, Gilliland, Weissig and Westbrook2000). Several research groups, as well as thousands of computer users, have dedicated their computer power to distributed structure prediction of proteins with Folding@HOME, POEM@HOME, and Rosetta@HOME. Martin-Renom et al. (Reference Marti-Renom, Stuart, Fiser, Sanchez, Melo and Sali2000) estimated that almost one-third of all known protein sequences can have their structure predicted via homology modeling efforts. Schafferhans & Klebe (Reference Schafferhans and Klebe2001) found that accuracy increases when the results of several different homology modeling programs are combined. While predictive models are not as accurate as structures solved by x-ray crystallography or NMR, they enable the study of targets that are difficult to determine experimentally. No consensus has been reached on the best source of data for including structural flexibility in SBDD, and a number of studies clearly demonstrate the difficulty and variability in assessing appropriate conformations (Bolstad & Anderson, Reference Bolstad and Anderson2008; Damm & Carlson, Reference Damm and Carlson2007; Laughton et al. Reference Laughton, Orozco and Vranken2009; Philippopoulos & Lim, Reference Philippopoulos and Lim1999; Yang et al. Reference Yang, Eyal, Chennubhotla, Jee, Gronenborn and Bahar2007). Bolstad and Anderson (Reference Bolstad and Anderson2008) found that their docking protocol correctly placed all forty tested LcDHFR ligands into the crystal structure, while the scoring function correctly ranked only 80% of the ligands. Including an NMR structure, homology model, or an ensemble of structures increased the number of incorrectly placed ligands while decreasing or maintaining the correct ranking of affinities. The authors found a similar trend when docking 19 ligands to hDHFR. This suggested that both sampling and scoring failures might lead to inaccuracies within a study. Sampling failures may occur without careful selection of flexibility information, but scoring may also negatively impact the results if the function cannot correctly classify binding poses.

The following sections are organized based on different techniques for including protein flexibility in docking and de novo design (Table 1). Soft docking allows limited discovery of new conformation space while docking with flexible side-chain samples observed conformational space. More recent techniques incorporate ensembles of structures to allow for a greater measure of native flexibility or to reveal new conformations. By relying on several protein conformations, new chemical space can be explored, allosteric sites can be discovered, and more accurate binding conformations can be generated. The accepted metric for a well-docked pose is ⩽2.0 Å RMSD from the experimentally determined structure. Protein flexibility can be incorporated into a binding-site model before docking is initiated or it can be allowed post-docking through refinement of the bound complex. Docking can be performed against an average structure or ensemble, a set of receptor conformations (multiple receptor conformations (MRC)), or with conformations generated ‘on the fly’ (induced-fit docking (IFD)).

Table 1.

Techniques for incorporating full protein flexibility into docking approaches

Experimental methods for SBDD inherently include protein–ligand flexibility. However, lead development and design that is based solely on the experimental data for ligand binding can be time-consuming and prohibitively expensive since every lead optimization step must be tested. Fragment-based drug design (FBDD) has enabled the application of experimental methods to searching the interaction potential along the protein surface with small organic probes. The use of fragment-based studies in SBDD enables the identification of ligands with high affinity for targets that have been traditionally difficult to characterize. These insights have greatly influenced computational approaches to protein flexibility, where fragment mapping can be applied to flexible receptors with functional group probes to find binding hot spots and identify specific interaction potentials. The use of molecular dynamics and/or probe minimizations to search for important binding sites with small probes has been frequently applied in SBDD projects because the results enable the development of a pharmacophore model for use in virtual screening studies.

4.1 Soft docking

It was not until the mid-1990s that researchers began to tackle the problem of receptor degrees of freedom in SBDD. The original technique involved soft docking, which accounted for minimal backbone movement, and side-chain flexibility through attenuation of the Lennard–Jones repulsion term between the rigid protein and ligand (Jiang & Kim, Reference Jiang and Kim1991). This allowed for some penetration between the ligand and protein, and was followed by a rigid-body minimization. Soft docking increased the number of high scoring hits compared to the use of a single structure for docking. A study of T4 lysozyme and aldose reductase showed that soft docking was superior to docking to a static structure (Ferrari et al. Reference Ferrari, Wei, Costantino and Shoichet2004). However, the authors also found that serial docking to a structural ensemble of rigid receptors identified six novel compounds targeting aldose reductase that soft docking had ranked poorly. Of these six compounds, four were experimentally verified and one had a low micromolar IC₅₀ value. Ferrari et al. conjectured that the differences in ranking were due to the differences in the scoring function applied for soft versus serial docking.

Soft docking is a fast and easy method for including some protein flexibility into docking studies, and as a result is it included as a step in many of the methods presented below.

4.2 Relaxation methods

Refinement of the docked complex is another simple approach that adjusts for protein flexibility by modeling induced-fit effects. The incorporation of protein flexibility on the ‘back end’ can only be done when the docking technique was based on an all-atom structure; it cannot be performed on a grid, or in any other situation where the protein is not explicitly represented. Monte Carlo (MC) or MD simulations are a popular choice for refinement because they enable the optimization of docked poses, investigation of solvent effects, examination of the kinetic stability, and prediction of ΔG _bind from physics-based scoring functions. Refinement is frequently performed as a final step in many of the docking approaches discussed here; Table 2 lists protein–ligand docking methods that limit the inclusion of full flexibility to the refinement stage (Mizutani et al. Reference Mizutani, Takamatsu, Ichinose, Nakamura and Itai2006; Pencheva et al. Reference Pencheva, Lagorce, Pajeva, Villoutreix and Miteva2008; Taylor et al. Reference Taylor, Jewsbury and Essex2003).

Table 2.

Studies including full receptor flexibility after docking (post-docking refinement)

4.3 Docking with flexible side chains

Early efforts at on-the-fly docking focused on the side chains, with the use of a rotamer library based on backbone dihedral angles to describe protein flexibility. In 1993, rotamer libraries were first used to predict side-chain conformations while studying the protein folding problem (Dunbrack & Karplus, Reference Dunbrack and Karplus1993). In 1994, the use of side-chain rotamer libraries was extended to protein–ligand docking with the examination of trypsin-benzamidine and antibody McPC 603-phosphocholine binding through a modified version of AMBER 4.0 (Leach). The study restricted the protein backbone completely while the side chains within 10 Å of the ligand were allowed for sampling a set of discrete rotameric states. Leach found that the presence of the ligand revealed additional accessible conformational states by modulating the protein's potential energy surface. Similar discrete sampling methods for modeling side-chain flexibility in docking include SOFTSPOTS/PLASTIC (Anderson et al. Reference Anderson, O'Neil, Surti and Stroud2001), as presented in Table 3a.

Table 3a.

Docking studies that included side-chain flexibility through rotamer libraries. Caveats for each method are similar in that flexibility of the side chains alone has limited success in describing receptor motion for most binding events and the implementation of a rotamer library limits the conformational space that can be sampled

Side-chain motion upon ligand binding has also been approached through continuum sampling methods. Alberts et al. (Reference Alberts, Todorov and Dean2005a) found that random movement of side chains was preferable to a rotamer library for successful docking results. Techniques that apply some form of continuum sampling to side-chain motion are listed in Table 3b, including AutoDock 4.0 (Morris et al. Reference Morris, Huey, Lindstrom, Sanner, Belew, Goodsell and Olson2009), Dynasite/GOLD (Kallblad & Dean, Reference Kallblad and Dean2003), FlexX (Rarey et al. Reference Rarey, Kramer, Lengauer and Klebe1996), ICM/MC (Abagyan et al. Reference Abagyan, Totrov and Kuznetsov1994), Mining Minima (Kairys & Gilson, Reference Kairys and Gilson2002), Skelgen (Alberts et al. Reference Alberts, Todorov and Dean2005a, Reference Alberts, Todorov, Kallblad and Deanb), and SLIDE (Schnecke & Kuhn, Reference Schnecke and Kuhn2000; Zavodszky et al. Reference Zavodszky, Sanschagrin, Korde and Kuhn2002; Zavodszky & Kuhn, Reference Zavodszky and Kuhn2005). A rotamer library can be limited by the information it contains, while continuum sampling may locate unrealistic conformational states.

Table 3b.

Docking studies where selected side-chain flexibility was accounted for through continuum sampling. This approach to flexibility at the side-chain positions may explore new space but cannot describe large-scale motion

Allowing side-chain flexibility is less resource intensive than full flexibility methods and it enabled some conformational variability through the exploration of low-energy orientations of side chains. However, incorporating side-chain flexibility has failed in cases of proteins with large-scale hinge or loop rearrangements or even backbone-dependent movement of side chains, neither of which can be taken into account by discrete or continuum sampling.

5.1 Conformational generation on the fly

Significant progress has been made since the initial implementation of flexible side chains in docking. Sampling motion of the side chains on the fly increases the potential energy space but can still overlook global conformational shifts. Allowing for conformational changes during docking can be a highly accurate technique for modeling bound poses of protein–ligand complexes. IFD is important because it can allow the docking simulation to search new conformational space; however, this sampling of receptor and ligand degrees of freedom is also quite computationally intensive, which limits its application in large-scale virtual screening studies.

The first on-the-fly method to expand beyond conformational sampling of side chains without relying upon a composite representation was based on DOCK 4.0 and ITScore: the ensemble docking algorithm (Huang & Zou, Reference Huang and Zou2007). The authors included the conformational ensemble as an adjustable variable in the IFD algorithm as opposed to performing serial docking across an entire ensemble. In addition to the usual six degrees of freedom for ligand flexibility, the authors added a protein conformation term, which allowed the protocol to simultaneously optimize the ligand conformation and choose the best protein structure for binding. In order to initially orient the ligand, a reference conformation comprising the largest possible binding site was generated from the combination of existing conformations. Ten separate protein ensembles were used for development and were obtained from the set by Cavasotto & Abagyan (Reference Cavasotto and Abagyan2004), the set by Claussen et al. (Reference Claussen, Buning, Rarey and Lengauer2001), and two HIV-1 protease (HIVp) ensembles. Huang and Zou achieved an average success rate of 92% based on both pose and energy score, where the energy score was required to be within 1.0 kcal/mol of the native score or better and the five highest-ranked poses were included using a threshold for success of 2.5 Å RMSD from the experimental pose. They noted that docking accuracy based on binding pose alone was insufficient, as demonstrated by the influence of the binding energy upon results with rigid docking. Their ensemble method ran for a similar amount of time as one rigid docking calculation and was 12.9 times faster than serial docking (data averaged over ten systems, with an average of 10.5 structures).

RosettaLigand was introduced as an extension of the protein–protein docking program, RosettaDock (Meiler & Baker, Reference Meiler and Baker2006). In this paper, side chains were repacked during docking through the implementation of a rotamer library, while the vdW repulsion term was softened and ligand conformations were perturbed with MC. Their benchmark set included 100 complexes for native docking and 20 for cross-docking. They found a success rate (RMSD <2.0 Å) of 80% for native docking and 70% for cross-docking; the authors noted that their RMSD metric accounted for hydrogen atoms and side chains in the binding site, thus capturing small changes in conformation overlooked by other methods. RosettaLigand was updated to include backbone sampling through a gradient-based minimization of the torsion angles for both protein and ligand as the final step in their docking procedure (Davis & Baker, Reference Davis and Baker2009; Davis et al. Reference Davis, Raha, Head and Baker2009). Although RosettaLigand has generated good results, with protein flexibility included, it requires 40–80 processor hours to generate a bound pose for the receptor–ligand complex and therefore is too time-consuming for use with large compound libraries.

Other IFD methods for including protein flexibility include modeling flexibility for a limited number of receptor residues, Table 4a: GLIDE/PRIME (Sherman et al. Reference Sherman, Beard and Farid2006a, Reference Sherman, Day, Jacobson, Friesner and Faridb), SCaRE (Bottegoni et al. Reference Bottegoni, Kufareva, Totrov and Abagyan2008); discrete sampling, Table 4b: FlexX-Ensemble (Claussen et al. Reference Claussen, Buning, Rarey and Lengauer2001), a modified DOCK protocol (Wei et al. Reference Wei, Weaver, Ferrari, Matthews and Shoichet2004), Fleksy (Nabuurs et al. Reference Nabuurs, Wagener and De Vlieg2007), FLIPDock (Zhao & Sanner, Reference Zhao and Sanner2007, Reference Zhao and Sanner2008), FITTED (Corbeil et al. Reference Corbeil, Englebienne and Moitessier2007, Reference Corbeil, Englebienne, Yannopoulos, Chan, Das, Bilimoria, L'Heureux and Moitessier2008; Corbeil & Moitessier, Reference Corbeil and Moitessier2009), 4D Docking (Bottegoni et al. Reference Bottegoni, Kufareva, Totrov and Abagyan2009); and continuum sampling, Table 4c: F-DycoBlock (Zhu et al. Reference Zhu, Fan, Liu and Shi2001), PC-RELAX (Zacharias, Reference Zacharias2004, Reference Zacharias2008), RosettaLigand (Davis & Baker, Reference Davis and Baker2009), and implicit MD (Huang & Wong, Reference Huang and Wong2009).

Table 4a.

Docking studies that modeled 2–3 residues to represent receptor flexibility during the performance of IFD

Table 4b.

Docking studies that applied discrete sampling of flexibility through panels of receptor structures during the performance of IFD

Table 4c.

Docking studies that applied molecular mechanics or minimization to incorporate receptor flexibility during IFD

6.1 Grid-based ensembles

The first docking method to use a composite grid was performed through DOCK3.5 in order to evaluate the impact of representing conformational variability as a structural ensemble on binding pose results (Knegtel et al. Reference Knegtel, Kuntz and Oshiro1997). With HIVp, ras p21, retinol binding protein, and uteroglobin as their test cases, the authors compared the capacity of standard grids to energy- and geometry-weighted average grids to reproduce known binding poses and affinities. They found that docking scores were sensitive to the grid spacing and threshold parameters used to define the composite grid. The performance of the geometry-weighted grid was less dependent on the protein conformations that were available than the energy-weighted, to the extent that the geometry-weighted grid placed known ligands for HIVp within the top 21%, while the energy-weighted grid placed them in the top 33%. This occurred because the geometry-weighted grid ignored the repulsive potentials between flexible atoms, while the energy-weighted grid simply smoothed the repulsive potential, which can result in the retrieval of paradoxical inhibitors.

The incorporation of protein flexibility through docking to an interaction grid that represents the receptor ensemble is a common approach for SBDD, and additional methods are presented in Table 5, including those based on AutoDock 3.0 (Osterberg et al. Reference Osterberg, Morris, Sanner, Olson and Goodsell2002), Flog (Broughton, Reference Broughton2000), GRID/CPCA (Kastenholz et al. Reference Kastenholz, Pastor, Cruciani, Haaksma and Fox2000), IFREDA (Cavasotto & Abagyan, Reference Cavasotto and Abagyan2004), ISCD (Zentgraf et al. Reference Zentgraf, Fokkens and Sotriffer2006), and sets of consensus structures (van Westen et al. Reference van Westen, Wegner, Bender, Ijzerman and van Vlijmen2010).

Table 5.

Studies including receptor flexibility through grid-based ensemble docking

6.2 Structure-based ensembles

The multiple copy simultaneous search (MCSS) methodology was the earliest computational approach for mapping binding sites with functional group probes (Miranker & Karplus, Reference Miranker and Karplus1991). The authors simultaneously minimized or quenched the probes by MD to the binding site of the hemagglutinin protein to identify favorable minima and found some correspondences between the positions of the minima and the functional groups on the ligand, sialic acid. While the authors noted that their method could account for a limited amount of side-chain flexibility by including multiple copies of the side chains, the published work was performed against a rigid structure.

Stultz and Karplus (Reference Stultz and Karplus1999) explored the influence of protein flexibility on the results from MCSS where five different protocols for placing two functional group probes, methanol and methyl ammonium, were used to search the binding surface of HIVp. MCSS calculations were performed for 1: a minimized crystal structure, 2: a conformation generated from quenched MD of the initial structure, 3: an unrestrained structure, 4: the output of 1 subjected to quenched MD with functional groups restrained and multiple side-chain copies, and 5: a different rigid crystal structure of HIVp. The authors found that their results with protocols 4 and 5 yielded more favorable interaction energies than the reference protocol 1. Protocol 4 gave valuable information for the ligand design, while the different low-energy minima from protocol 5 supported the idea of performing MCSS against an ensemble of structures and comparing the results; therefore, the authors suggested a combination of 4 and 5. However, one of the limitations of MCSS and other ensemble methods is that they do not account for entropy or solvation during surface mapping, which hampers discrimination between druggable and irrelevant minima. Furthermore, as Schubert & Stultz (Reference Schubert and Stultz2009) point out in their review of MCSS, functionality maps are difficult to combine across different structures of a protein.

An original technique for structure-based ensemble docking was introduced by Carlson et al. (Reference Carlson, Masukawa, Rubins, Bushman, Jorgensen, Lins, Briggs and McCammon2000) as the dynamic pharmacophore model (Carlson et al. Reference Carlson, Masukawa and McCammon1999), later termed the multiple protein structure (MPS) method. Two pharmacophores for HIV integrase were generated, a static model based on a rigid crystal structure and a dynamic model based on MD snapshots initiated from the same crystal structure. The MPS method mapped the dynamic binding surface of the protein with common functional groups by running a series of multi-unit search for interacting conformers (MUSIC) simulations (Fig. 3), wherein hundreds of probes were used to flood the protein surface and then minimized by MC sampling. During MUSIC simulations, all of the probe–probe interactions were ignored, which allowed probes to interact with the rugged binding potential of the protein and cluster into hot-spot positions along the protein surface. These cluster sites were used to develop a pharmacophore for virtual screening, and where the rigid model was unable to locate any of the test case inhibitors, the dynamic pharmacophore model not only identified known inhibitors, but it also predicted new inhibitors that were confirmed by experiment. The MPS method is one of the few ensemble-based methods that has been successfully used to find novel leads with demonstrated activity (Bowman et al. Reference Bowman, Nikolovska-Coleska, Zhong, Wang and Carlson2007; Damm et al. Reference Damm, Ung, Quintero, Gestwicki and Carlson2008).

Fig. 3.

MUSIC simulations of benzene probes (gray) to the surface of an HIVp monomer (purple). The initial stages are universal to mapping methods, but the combination of the structures in the fourth and final frames are unique to the MPS method. Five hundred probes were flooded onto the protein surface, minimized, and clustered together. Parent probes were identified for each monomer conformation, the HIVp monomers were aligned through a weighted-RMSD function, and then the parent probes were clustered together and retained when at least half of the protein conformations contained a probe at that site. This resulted in low-energy probe clusters (colored red through purple in the final frame) of benzene that could be combined with ethane and methanol clusters to develop a pharmacophore model.

McCammon and co-workers developed the relaxed complex scheme (RCS) as a technique for incorporating receptor flexibility prior to docking (Lin et al. Reference Lin, Perryman, Schames and McCammon2002, Reference Lin, Perryman, Schames and McCammon2003; McCammon, Reference McCammon2005). Initially the method was developed to generate conformational ensembles of the apo state through MD, which were then used in AutoDock to screen compound mini-libraries and score the receptor–ligand complexes. The 2003 study by Lin et al. examined FKBP-12 by RCS and implemented a final refinement step with MM/PBSA (Kollman et al. Reference Kollman, Massova, Reyes, Kuhn, Huo, Chong, Lee, Lee, Duan, Wang, Donini, Cieplak, Srinivasan, Case and Cheatham2000) to yield final results. Further updates to this method have included extension to virtual screening and reduction of the conformational ensemble to a representative configurational set (Amaro et al. Reference Amaro, Baron and McCammon2008). RCS was successfully used to identify cryptic binding pockets and/or lead inhibitors of HIVp (Perryman et al. Reference Perryman, Lin and McCammon2006), avian influenza neuraminidase (Cheng et al. Reference Cheng, Amaro, Xu, Li, Arzberger and McCammon2008), acetylcholine binding protein (Babakhani et al. Reference Babakhani, Talley, Taylor and McCammon2009), DNA polymerase β (Barakat & Tuszynski, Reference Barakat and Tuszynski2010), MDM2/MDMX (Barakat et al. Reference Barakat, Mane, Friesen and Tuszynski2010), and cruzain (Durrant et al. Reference Durrant, Keranen, Wilson and McCammon2010). Cheng et al. noted that conformational sampling markedly influenced the re-ranking of compounds based on interactions gained or lost compared to the crystal structure. The authors noted that the most dominant structures did not have the lowest binding energy, but given the limitations of MD, they believed that the dominant clusters were good representatives of the ensemble. Babakhani et al. suggested that both the scoring function and the conformational variability contributed to several failures in reproducing known binding modes.

To determine the best process for developing a structural ensemble, Rueda et al. (Reference Rueda, Bottegoni and Abagyan2010) used ICM to examine the same set of 1068 conformations from 99 pharmaceutically relevant proteins that was used in the 4D docking study. The authors judged performance in virtual screening based on the area under the curve (AUC) of a ROC plot and found that holo conformations yielded improved AUC values compared to apo conformations. It was interesting that the authors found that ensembles of holo plus apo conformations did not significantly improve results compared to holo-only ensembles. In cases without a holo structure, the authors recommended the use of an apo structure with the largest pocket volume. The authors proposed the use of a ligand-guided approach to find the optimal protein conformation, but stated that docking to an ensemble would be an acceptable substitute in the absence of ligand activity data. Several groups have found that there tends to be a single receptor conformation that grants the best performance in docking studies (Barril & Morley, Reference Barril and Morley2005; Rueda et al. Reference Rueda, Bottegoni and Abagyan2010).

NMA has been frequently used as a tool for examining protein flexibility relevant to ligand binding (Cavasotto et al. Reference Cavasotto, Kovacs and Abagyan2005; Keseru & Kolossvary, Reference Keseru and Kolossvary2001; May & Zacharias, Reference May and Zacharias2005, Reference May and Zacharias2008; May et al. Reference May, Sieker and Zacharias2008; Zacharias, Reference Zacharias2004; Zacharias & Sklenar, Reference Zacharias and Sklenar1999; Zavodszky & Kuhn, Reference Zavodszky and Kuhn2005). It is important to note that low-frequency modes identify large-scale motions, while high-frequency modes identify small-scale motion. Rueda et al. (Reference Rueda, Bottegoni and Abagyan2009) showed that they could use NMA on elastic network models for all heavy atoms of the receptor in order to derive an ensemble of protein conformations and improve results from cross-docking in ICM. Fourteen proteins, each with two structures and their cognate ligands, were used as the benchmark test set. Optimal results involved the generation of 100 different conformations from NMA on all heavy atoms within 10 Å of the ligand-binding site. Larger conformational ensembles, such as those with 200 members, negatively affected false positive rates. Another NMA-based method for conformational selection was based on cyclin dependent kinase-2 (CDK-2), because it had a sufficient amount of available experimental data on ligand binding for validation of the NMA protocol (Sperandio et al. Reference Sperandio, Mouawad, Pinto, Villoutreix, Perahia and Miteva2010). Unlike most other NMA-based methods, Sperandio et al. included all of the protein atoms in the calculation in order to better represent conformational changes and found that too many protein conformations negatively impacted docking results by recovering additional false positives. They suggested addressing this by developing scoring protocols specifically tailored to the use of a conformational ensemble.

Of the variety of techniques that have been developed to improve binding predictions for protein–ligand complexes, the use of a conformational ensemble is perhaps the most common approach. Representative cases and novel methods based on structure-based ensemble docking are listed in Tables 6a–6d. Crystal structure and NMR-based ensemble methods include docking to a large set of conformations (Barril & Morley, Reference Barril and Morley2005; Craig et al. Reference Craig, Essex and Spiegel2010), and ICM (Rueda et al. Reference Rueda, Bottegoni and Abagyan2010), as listed in Table 6a. Procedures that rely upon molecular dynamics to provide receptor flexibility include MCSS (Stultz & Karplus, Reference Stultz and Karplus1999), MPS (Carlson et al. Reference Carlson, Masukawa and McCammon1999, Reference Carlson, Masukawa, Rubins, Bushman, Jorgensen, Lins, Briggs and McCammon2000), MD+LigBuilder (Mustata & Briggs, Reference Mustata and Briggs2002; Mustata et al. Reference Mustata, Soares and Briggs2003), RCS (Lin et al. Reference Lin, Perryman, Schames and McCammon2003), MD+AutoDock (Frembgen-Kesner & Elcock, Reference Frembgen-Kesner and Elcock2006), Flo98 (Popov et al. Reference Popov, Yee and Anderson2007), a reduced receptor ensemble (Landon et al. Reference Landon, Amaro, Baron, Ngan, Ozonoff, McCammon and Vajda2008), REMD/PLOP (Wong & Jacobson, Reference Wong and Jacobson2008), the ensemble reduction method (Bolstad & Anderson, Reference Bolstad and Anderson2009), enhanced molecular docking (Kranjc et al. Reference Kranjc, Bongarzone, Rossetti, Biarnes, Cavalli, Bolognesi, Roberti, Legname and Carloni2009) and MD+Glide (Nichols et al. Reference Nichols, Baron, Ivetac and McCammon2011), as listed in Table 6b. Ensembles have also been generated for docking through NMA, as demonstrated by FIRST/ROCK/SLIDE (Zavodszky et al. Reference Zavodszky, Lei, Thorpe, Day and Kuhn2004), low-frequency NMA (Cavasotto et al. Reference Cavasotto, Kovacs and Abagyan2005), ENM (Lindahl & Delarue, Reference Lindahl and Delarue2005), EN–NMA (Rueda et al. Reference Rueda, Bottegoni and Abagyan2009), and all atom NMA (Sperandio et al. Reference Sperandio, Mouawad, Pinto, Villoutreix, Perahia and Miteva2010), as listed in Table 6c. Finally, ligand-induced flexibility has also been used as another source for conformers, as shown in Table 6d and applied by MCSA-PCR (Ota & Agard, Reference Ota and Agard2001), Skelgen (Firth-Clark et al. Reference Firth-Clark, Willems, Williams and Harris2006), QSiCR (Subramanian et al. Reference Subramanian, Sharma and B-Rao2006, Reference Subramanian, Sharma and B-Rao2008), PhE/SVM (Leong, Reference Leong2007; Leong & Chen, Reference Leong and Chen2008; Leong et al. Reference Leong, Chen, Chen and Chen2009), structure prediction (Bisson et al. Reference Bisson, Cheltsov, Bruey-Sedano, Lin, Chen, Goldberger, May, Christopoulos, Dalton, Sexton, Zhang and Abagyan2007), and Surflex (Jain, Reference Jain2003, Reference Jain2007, Reference Jain2009). Although listed in Table 5 due to its incorporation of protein flexibility within an interaction grid, IFREDA also relies upon ligand-induced conformations for contributing flexibility.

Table 6a.

Studies including full receptor flexibility through collection of a crystal or NMR structural ensemble

Table 6b.

Studies including full receptor flexibility through a structural ensemble previously generated by MD

Table 6c.

Studies including full receptor flexibility through a normal mode analysis, or other non-MD-based conformational generation, prior to docking

Table 6d.

Studies including full receptor flexibility through ligand-induced flexibility