Ensemble docking of multiple protein structures: considering protein structural variations in molecular docking

One approach to incorporate protein flexibility in molecular docking is the use of an ensemble consisting of multiple protein structures. Sequentially docking each ligand into a large number of protein structures is computationally too expensive to allow large-scale database screening. It is challenging to achieve a good balance between docking accuracy and computational efficiency. In this work, we have developed a fast, novel docking algorithm utilizing multiple protein structures, referred to as ensemble docking, to account for protein structural variations. The algorithm can simultaneously dock a ligand into an ensemble of protein structures and automatically select an optimal protein structure that best fits the ligand by optimizing both ligand coordinates and the conformational variable m, where m represents the m-th structure in the protein ensemble. The docking algorithm was validated on 10 protein ensembles containing 105 crystal structures and 87 ligands in terms of binding mode and energy score predictions. A success rate of 93% was obtained with the criterion of root-mean-square deviation <2.5 A if the top five orientations for each ligand were considered, comparable to that of sequential docking in which scores for individual docking are merged into one list by re-ranking, and significantly better than that of single rigid-receptor docking (75% on average). Similar trends were also observed in binding score predictions and enrichment tests of virtual database screening. The ensemble docking algorithm is computationally efficient, with a computational time comparable to that for docking a ligand into a single protein structure. In contrast, the computational time for the sequential docking method increases linearly with the number of protein structures in the ensemble. The algorithm was further evaluated using a more realistic ensemble in which the corresponding bound protein structures of inhibitors were excluded. The results show that ensemble docking successfully predicts the binding modes of the inhibitors, and discriminates the inhibitors from a set of noninhibitors with similar chemical properties. Although multiple experimental structures were used in the present work, our algorithm can be easily applied to multiple protein conformations generated by computational methods, and helps improve the efficiency of other existing multiple protein structure(MPS)-based methods to accommodate protein flexibility.     

Flexible ligand docking using conformational ensembles.

Molecular docking algorithms suggest possible structures for molecular complexes. They are used to model biological function and to discover potential ligands. A present challenge for docking algorithms is the treatment of molecular flexibility. Here, the rigid body program, DOCK, is modified to allow it to rapidly fit multiple conformations of ligands. Conformations of a given molecule are pre-calculated in the same frame of reference, so that each conformer shares a common rigid fragment with all other conformations. The ligand conformers are then docked together, as an ensemble, into a receptor binding site. This takes advantage of the redundancy present in differing conformers of the same molecule. The algorithm was tested using three organic ligand protein systems and two protein-protein systems. Both the bound and unbound conformations of the receptors were used. The ligand ensemble method found conformations that resembled those determined in X-ray crystal structures (RMS values typically less than 1.5 A). To test the method's usefulness for inhibitor discovery, multi-compound and multi-conformer databases were screened for compounds known to bind to dihydrofolate reductase and compounds known to bind to thymidylate synthase. In both cases, known inhibitors and substrates were identified in conformations resembling those observed experimentally. The ligand ensemble method was 100-fold faster than docking a single conformation at a time and was able to screen a database of over 34 million conformations from 117,000 molecules in one to four CPU days on a workstation.     

Conformer selection and induced fit in flexible backbone protein-protein docking using computational and NMR ensembles

Accommodating backbone flexibility continues to be the most difficult challenge in computational docking of protein-protein complexes. Towards that end, we simulate four distinct biophysical models of protein binding in RosettaDock, a multiscale Monte-Carlo-based algorithm that uses a quasi-kinetic search process to emulate the diffusional encounter of two proteins and to identify low-energy complexes. The four binding models are as follows: (1) key-lock (KL) model, using rigid-backbone docking; (2) conformer selection (CS) model, using a novel ensemble docking algorithm; (3) induced fit (IF) model, using energy-gradient-based backbone minimization; and (4) combined conformer selection/induced fit (CS/IF) model. Backbone flexibility was limited to the smaller partner of the complex, structural ensembles were generated using Rosetta refinement methods, and docking consisted of local perturbations around the complexed conformation using unbound component crystal structures for a set of 21 target complexes. The lowest-energy structure contained > 30% of the native residue-residue contacts for 9, 13, 13, and 14 targets for KL, CS, IF, and CS/IF docking, respectively. When applied to 15 targets using nuclear magnetic resonance ensembles of the smaller protein, the lowest-energy structure recovered at least 30% native residue contacts in 3, 8, 4, and 8 targets for KL, CS, IF, and CS/IF docking, respectively. CS/IF docking of the nuclear magnetic resonance ensemble performed equally well or better than KL docking with the unbound crystal structure in 10 of 15 cases. The marked success of CS and CS/IF docking shows that ensemble docking can be a versatile and effective method for accommodating conformational plasticity in docking and serves as a demonstration for the CS theory—that binding-competent conformers exist in the unbound ensemble and can be selected based on their favorable binding energies.     

In silico assessment of potential druggable pockets on the surface of α1-antitrypsin conformers

The search for druggable pockets on the surface of a protein is often performed on a single conformer, treated as a rigid body. Transient druggable pockets may be missed in this approach. Here, we describe a methodology for systematic in silico analysis of surface clefts across multiple conformers of the metastable protein α(1)-antitrypsin (A1AT). Pathological mutations disturb the conformational landscape of A1AT, triggering polymerisation that leads to emphysema and hepatic cirrhosis. Computational screens for small molecule inhibitors of polymerisation have generally focused on one major druggable site visible in all crystal structures of native A1AT. In an alternative approach, we scan all surface clefts observed in crystal structures of A1AT and in 100 computationally produced conformers, mimicking the native solution ensemble. We assess the persistence, variability and druggability of these pockets. Finally, we employ molecular docking using publicly available libraries of small molecules to explore scaffold preferences for each site. Our approach identifies a number of novel target sites for drug design. In particular one transient site shows favourable characteristics for druggability due to high enclosure and hydrophobicity. Hits against this and other druggable sites achieve docking scores corresponding to a K(d) in the μM-nM range, comparing favourably with a recently identified promising lead. Preliminary ThermoFluor studies support the docking predictions. In conclusion, our strategy shows considerable promise compared with the conventional single pocket/single conformer approach to in silico screening. Our best-scoring ligands warrant further experimental investigation.     

Simultaneous single-structure and bundle representation of protein NMR structures in torsion angle space

A method is introduced to represent an ensemble of conformers of a protein by a single structure in torsion angle space that lies closest to the averaged Cartesian coordinates while maintaining perfect covalent geometry and on average equal steric quality and an equally good fit to the experimental (e.g. NMR) data as the individual conformers of the ensemble. The single representative ‘regmean structure’ is obtained by simulated annealing in torsion angle space with the program CYANA using as input data the experimental restraints, restraints for the atom positions relative to the average Cartesian coordinates, and restraints for the torsion angles relative to the corresponding principal cluster average values of the ensemble. The method was applied to 11 proteins for which NMR structure ensembles are available, and compared to alternative, commonly used simple approaches for selecting a single representative structure, e.g. the structure from the ensemble that best fulfills the experimental and steric restraints, or the structure from the ensemble that has the lowest RMSD value to the average Cartesian coordinates. In all cases our method found a structure in torsion angle space that is significantly closer to the mean coordinates than the alternatives while maintaining the same quality as individual conformers. The method is thus suitable to generate representative single structure representations of protein structure ensembles in torsion angle space. Since in the case of NMR structure calculations with CYANA the single structure is calculated in the same way as the individual conformers except that weak positional and torsion angle restraints are added, we propose to represent new NMR structures by a ‘regmean bundle’ consisting of the single representative structure as the first conformer and all but one original individual conformers (the original conformer with the highest target function value is discarded in order to keep the number of conformers in the bundle constant). In this way, analyses that require a single structure can be carried out in the most meaningful way using the first model, while at the same time the additional information contained in the ensemble remains available.     

Protein flexibility in ligand docking and virtual screening to protein kinases

The main complicating factor in structure-based drug design is receptor rearrangement upon ligand binding (induced fit). It is the induced fit that complicates cross-docking of ligands from different ligand-receptor complexes. Previous studies have shown the necessity to include protein flexibility in ligand docking and virtual screening. Very few docking methods have been developed to predict the induced fit reliably and, at the same time, to improve on discriminating between binders and non-binders in the virtual screening process.We present an algorithm called the ICM-flexible receptor docking algorithm (IFREDA) to account for protein flexibility in virtual screening. By docking flexible ligands to a flexible receptor, IFREDA generates a discrete set of receptor conformations, which are then used to perform flexible ligand-rigid receptor docking and scoring. This is followed by a merging and shrinking step, where the results of the multiple virtual screenings are condensed to improve the enrichment factor. In the IFREDA approach, both side-chain rearrangements and essential backbone movements are taken into consideration, thus sampling adequately the conformational space of the receptor, even in cases of large loop movements.As a preliminary step, to show the importance of incorporating protein flexibility in ligand docking and virtual screening, and to validate the merging and shrinking procedure, we compiled an extensive small-scale virtual screening benchmark of 33 crystal structures of four different protein kinases sub-families (cAPK, CDK-2, P38 and LCK), where we obtained an enrichment factor fold-increase of 1.85±0.65 using two or three multiple experimental conformations. IFREDA was used in eight protein kinase complexes and was able to find the correct ligand conformation and discriminate the correct conformations from the “misdocked” conformations solely on the basis of energy calculation. Five of the generated structures were used in the small-scale virtual screening stage and, by merging and shrinking the results with those of the original structure, we show an enrichment factor fold increase of 1.89±0.60, comparable to that obtained using multiple experimental conformations.Our cross-docking tests on the protein kinase benchmark underscore the necessity of incorporating protein flexibility in both ligand docking and virtual screening. The methodology presented here will be extremely useful in cases where few or no experimental structures of complexes are available, while some binders are known.     

A fast protein-ligand docking algorithm based on hydrogen bond matching and surface shape complementarity

With the rapid development of structural determination of target proteins for human diseases, high throughout virtual screening based drug discovery is gaining popularity gradually. In this paper, a fast docking algorithm (H-DOCK) based on hydrogen bond matching and surface shape complementarity was developed. In H-DOCK, firstly a divide-and-conquer strategy based enumeration approach is applied to rank the intermolecular modes between protein and ligand by maximizing their hydrogen bonds matching, then each docked conformation of the ligand is calculated according to the matched hydrogen bonding geometry, finally a simple but effective scoring function reflecting mainly the van der Waals interaction is used to evaluate the docked conformations of the ligand. H-DOCK is tested for rigid ligand docking and flexible one, the latter is implemented by repeating rigid docking for multiple conformations of a small molecule and ranking all together. For rigid ligands, H-DOCK was tested on a set of 271 complexes where there is at least one intermolecular hydrogen bond, and H-DOCK achieved success rate (RMSD<2.0?Å) of 91.1%. For flexible ligands, H-DOCK was tested on another set of 93 complexes, where each case was a conformation ensemble containing native ligand conformation as well as 100 decoy ones generated by AutoDock [1], and the success rate reached 81.7%. The high success rate of H-DOCK indicates that the hydrogen bonding and steric hindrance can grasp the key interaction between protein and ligand. H-DOCK is quite efficient compared with the conventional docking algorithms, and it takes only about 0.14 seconds for a rigid ligand docking and about 8.25 seconds for a flexible one on average. According to the preliminary docking results, it implies that H-DOCK can be potentially used for large scale virtual screening as a pre-filter for a more accurate but less efficient docking algorithm.     

Minimal ensembles of side chain conformers for modeling protein–protein interactions

The goal of this article is to reduce the complexity of the side chain search within docking problems. We apply six methods of generating side chain conformers to unbound protein structures and determine their ability of obtaining the bound conformation in small ensembles of conformers. Methods are evaluated in terms of the positions of side chain end groups. Results for 68 protein complexes yield two important observations. First, the end‐group positions change less than 1 Å on association for over 60% of interface side chains. Thus, the unbound protein structure carries substantial information about the side chains in the bound state, and the inclusion of the unbound conformation into the ensemble of conformers is very beneficial. Second, considering each surface side chain separately in its protein environment, small ensembles of low‐energy states include the bound conformation for a large fraction of side chains. In particular, the ensemble consisting of the unbound conformation and the two highest probability predicted conformers includes the bound conformer with an accuracy of 1 Å for 78% of interface side chains. As more than 60% of the interface side chains have only one conformer and many others only a few, these ensembles of low‐energy states substantially reduce the complexity of side chain search in docking problems. This approach was already used for finding pockets in protein–protein interfaces that can bind small molecules to potentially disrupt protein–protein interactions. Side‐chain search with the reduced search space will also be incorporated into protein docking algorithms. Proteins 2012. © 2011 Wiley Periodicals, Inc.     

Pre-docking filter for protein and ligand 3D structures

Virtual drug screening using protein-ligand docking techniques is a time-consuming process, which requires high computational power for binding affinity calculation. There are millions of chemical compounds available for docking. Eliminating compounds that are unlikely to exhibit high binding affinity from the screening set should speed-up the virtual drug screening procedure. We performed docking of 6353 ligands against twenty-one protein X-ray crystal structures. The docked ligands were ranked according to their calculated binding affinities, from which the top five hundred and the bottom five hundred were selected. We found that the volume and number of rotatable bonds of the top five hundred docked ligands are similar to those found in the crystal structures and corresponded with the volume of the binding sites. In contrast, the bottom five hundred set contains ligands that are either too large to enter the binding site, or too small to bind with high specificity and affinity to the binding site. A pre-docking filter that takes into account shapes and volumes of the binding sites as well as ligand volumes and flexibilities can filter out low binding affinity ligands from the screening sets. Thus, the virtual drug screening procedure speed is increased.     

Monte Carlo replica‐exchange based ensemble docking of protein conformations

A replica‐exchange Monte Carlo (REMC) ensemble docking approach has been developed that allows efficient exploration of protein–protein docking geometries. In addition to Monte Carlo steps in translation and orientation of binding partners, possible conformational changes upon binding are included based on Monte Carlo selection of protein conformations stored as ordered pregenerated conformational ensembles. The conformational ensembles of each binding partner protein were generated by three different approaches starting from the unbound partner protein structure with a range spanning a root mean square deviation of 1–2.5 Å with respect to the unbound structure. Because MC sampling is performed to select appropriate partner conformations on the fly the approach is not limited by the number of conformations in the ensemble compared to ensemble docking of each conformer pair in ensemble cross docking. Although only a fraction of generated conformers was in closer agreement with the bound structure the REMC ensemble docking approach achieved improved docking results compared to REMC docking with only the unbound partner structures or using docking energy minimization methods. The approach has significant potential for further improvement in combination with more realistic structural ensembles and better docking scoring functions. Proteins 2017; 85:924–937. © 2016 Wiley Periodicals, Inc.     

Effects of initial settings on computational protein–ligand docking accuracies for several docking programs

In this study, the influences of initial settings, i.e. initial conformations, configurations and docking parameters, on docking results were investigated. The conformations used in the study were generated by the CAMDAS program. After the conformational search calculations, five structures were selected from the conformer groups according to their conformation energies and root mean square deviations against crystal structures; for example, the lowest energy conformer, as well as the closest and farthest conformers to the crystal structure, was retrieved. Several docking parameter settings were used (default, high speed, generating 50 poses). In this study, docking calculations were conducted using the GOLD, eHiTS, AutoDock, AutoDock vina, FRED and DOCK programs. The success rates of GOLD, eHiTS and FRED were better than those of AutoDock, AutoDock vina and DOCK. The docking results using the farthest conformations were worse than those obtained using other conformations, indicating that some conformation search for the ligand molecule should be performed before the docking calculations.     

Sampling of the native conformational ensemble of myoglobin via structures in different crystalline environments

Proteins sample multiple conformational substates in their native environment, but the process of crystallization selects the conformers that allow for close packing. The population of conformers can be shifted by varying the environment through a range of crystallization conditions, often resulting in different space groups and changes in the packing arrangements. Three high resolution structures of myoglobin (Mb) in different crystal space groups are presented, including one in a new space group P6(1)22 and two structures in space groups P2(1)2(1)2(1) and P6. We compare coordinates and anisotropic displacement parameters (ADPs) from these three structures plus an existing structure in space group P2(1). While the overall changes are small, there is substantial variation in several external regions with varying patterns of crystal contacts across the space group packing arrangements. The structural ensemble containing four different crystal forms displays greater conformational variance (Calpha rmsd of 0.54-0.79 A) in comparison to a collection of four Mb structures with different ligands and mutations in the same crystal form (Calpha rmsd values of 0.28-0.37 A). The high resolution of the data enables comparison of both the magnitudes and directions of ADPs, which are found to be suppressed by crystal contacts. A composite dynamic profile of Mb structural variation from the four structures was compared with an independent structural ensemble developed from NMR refinement. Despite the limitations and biases of each method, the ADPs of the crystallographic ensemble closely match the positional variance from the solution NMR ensemble with linear correlation of 0.8. This suggests that crystal packing selects conformers representative of the solution ensemble, and several different crystal forms give a more complete view of the plasticity of a protein structure.     

Applying conformational selection theory to improve crossdocking efficiency in 3‐phosphoinositide dependent protein kinase‐1

The emerging picture of biomolecular recognition is that of conformational selection followed by induced‐fit. Conformational selection theory states that binding partners exist in various conformations in solution, with binding involving a “selection” between complementary conformers. In this study, we devise a docking protocol that mimics conformational selection in protein–ligand binding and demonstrate that it significantly enhances crossdocking accuracy over Glide's flexible docking protocol, which is widely used in the pharmaceutical industry. Our protocol uses a pregenerated conformational ensemble to simulate ligand flexibility. The ensemble was generated by thorough conformational sampling coupled with conformer minimization. The generated conformers were then rigidly docked in the active site of the protein along with a postdocking minimization step that allows limited induced fit effects to be modeled for the ligand. We illustrate the improved performance of our protocol through crossdocking of 31 ligands to cocomplexed proteins of the kinase 3‐phosphoinositide dependent protein kinase‐1 extracted from the crystal structures 1H1W (ATP bound), 1OKY (staurosporine bound) and 3QD0 (bound to a potent inhibitor). Consistent with conformational selection theory, the performance of our protocol was the best for crossdocking to the cognate protein bound to the natural ligand, ATP. Proteins 2014; 82:436–451. © 2013 Wiley Periodicals, Inc.     

Fluctuations in ion pairs and their stabilities in proteins

This report investigates the effect of systemic protein conformational flexibility on the contribution of ion pairs to protein stability. Toward this goal, we use all NMR conformer ensembles in the Protein Data Bank (1) that contain at least 40 conformers, (2) whose functional form is monomeric, (3) that are nonredundant, and (4) that are large enough. We find 11 proteins adhering to these criteria. Within these proteins, we identify 22 ion pairs that are close enough to be classified as salt bridges. These are identified in the high-resolution crystal structures of the respective proteins or in the minimized average structures (if the crystal structures are unavailable) or, if both are unavailable, in the "most representative" conformer of each of the ensembles. We next calculate the electrostatic contribution of each such ion pair in each of the conformers in the ensembles. This results in a comprehensive study of 1,201 ion pairs, which allows us to look for consistent trends in their electrostatic contributions to protein stability in large sets of conformers. We find that the contributions of ion pairs vary considerably among the conformers of each protein. The vast majority of the ion pairs interconvert between being stabilizing and destabilizing to the structure at least once in the ensembles. These fluctuations reflect the variabilities in the location of the ion pairing residues and in the geometric orientation of these residues, both with respect to each other, and with respect to other charged groups in the remainder of the protein. The higher crystallographic B-factors for the respective side-chains are consistent with these fluctuations. The major conclusion from this study is that salt bridges observed in crystal structure may break, and new salt bridges may be formed. Hence, the overall stabilizing (or, destabilizing) contribution of an ion pair is conformer population dependent.     

Towards in silico lead optimization: scores from ensembles of protein/ligand conformations reliably correlate with biological activity

Accurately ranking protein/ligand interactions and distinguishing subtle differences between homologous compounds in a virtual focused library in silico is essential in a structure-based drug discovery program. In order to establish a predictive model to design novel inhibitors of dihydrofolate reductase (DHFR) from the parasitic protozoa, Cryptosporidium hominis, we docked a series of 30 DHFR inhibitors with measured inhibition constants against the crystal structure of the protein. By including protein flexibility and averaging the energies of the 25 lowest protein/ligand conformers we obtained more accurate total nonbonded energies from which we calculated a predicted biological activity. The calculated and measured biological activities showed reliable correlations of 72.9%. Additionally, visual analysis of the ensemble of protein/ligand conformations revealed alternative ligand binding pockets in the active site. Using the same principles we then created a homology model of DHFR from Toxoplasma gondii and docked 11 inhibitors. A correlation of 50.2% between docking score and activity validates both the method and the model. The correlations presented here are particularly compelling considering the high structural similarity of the ligands and the fact that we have used structures derived from crystallographic data and homology modeling. These docking principles may be useful in any lead optimization study where accurate ranking of similar compounds is desired.     

In pursuit of virtual lead optimization: the role of the receptor structure and ensembles in accurate docking

Accurate ranking during in silico lead optimization is critical to drive the generation of new ligands with higher affinity, yet it is especially difficult because of the subtle changes between analogs. In order to assess the role of the structure of the receptor in delivering accurate lead ranking results, we docked a set of forty related inhibitors to structures of one species of dihydrofolate reductase (DHFR) derived from crystallographic, NMR solution data, and homology models. In this study, the crystal structures yielded the superior results: the compounds were placed in the active site in the conserved orientation and the docking scores for 80% percent of the compounds clustered into the same bins as the measured affinity. Single receptor structures derived from NMR data or homology models did not serve as accurate docking receptors. To our knowledge, these are the first experiments that assess ranking of homologous lead compounds using a variety of receptor structures. We then extended the study to investigate whether ensembles, either computationally or experimentally derived, of all of the single starting structures aid, hinder or have no effect on the performance of the starting template. Impressively, when ensembles of receptor structures derived from NMR data or homology models were employed, docking accuracy improved to a level equal to that of the high resolution crystal structures. The same experiments using a second species of DHFR and set of ligands confirm the results. A comparison of the structures of the individual ensemble members to the starting structures shows that the effect of the ensembles can be ascribed to protein flexibility in addition to absorption of computational error.     

Yucca: an efficient algorithm for small-molecule docking

In this paper, we present a new algorithm, which is based on an efficient heuristic for local search, for rigid protein-small-molecule docking. We tested our algorithm, called Yucca, on the recent 100-complex benchmark, using the conformer generator OMEGA to generate a set of low-energy conformers. The results showed that Yucca is competitive both in terms of algorithm efficiency and docking accuracy.