Scoring a diverse set of high-quality docked conformations: a metascore based on electrostatic and desolvation interactions

Predicting protein-protein interactions involves sampling and scoring docked conformations. Barring some large structural rearrangement, rapidly sampling the space of docked conformations is now a real possibility, and the limiting step for the successful prediction of protein interactions is the scoring function used to reduce the space of conformations from billions to a few, and eventually one high affinity complex. An atomic level free-energy scoring function that estimates in units of kcal/mol both electrostatic and desolvation interactions (plus van der Waals if appropriate) of protein-protein docked conformations is used to rerank the blind predictions (860 in total) submitted for six targets to the community-wide Critical Assessment of PRediction of Interactions (CAPRI; http://capri.ebi.ac.uk). We found that native-like models often have varying intermolecular contacts and atom clashes, making unlikely that one can construct a universal function that would rank all these models as native-like. Nevertheless, our scoring function is able to consistently identify the native-like complexes as those with the lowest free energy for the individual models of 16 (out of 17) human predictors for five of the targets, while at the same time the modelers failed to do so in more than half of the cases. The scoring of high-quality models developed by a wide variety of methods and force fields confirms that electrostatic and desolvation forces are the dominant interactions determining the bound structure. The CAPRI experiment has shown that modelers can predict valuable models of protein-protein complexes, and improvements in scoring functions should soon solve the docking problem for complexes whose backbones do not change much upon binding. A scoring server and programs are available at http://structure.pitt.edu.     

Protein docking using continuum electrostatics and geometric fit

The computer program DOT quickly finds low-energy docked structures for two proteins by performing a systematic search over six degrees of freedom. A novel feature of DOT is its energy function, which is the sum of both a Poisson-Boltzmann electrostatic energy and a van der Waals energy, each represented as a grid-based correlation function. DOT evaluates the energy of interaction for many orientations of the moving molecule and maintains separate lists scored by either the electrostatic energy, the van der Waals energy or the composite sum of both. The free energy is obtained by summing the Boltzmann factor over all rotations at each grid point. Three important findings are presented. First, for a wide variety of protein-protein interactions, the composite-energy function is shown to produce larger clusters of correct answers than found by scoring with either van der Waals energy (geometric fit) or electrostatic energy alone. Second, free-energy clusters are demonstrated to be indicators of binding sites. Third, the contributions of electrostatic and attractive van der Waals energies to the total energy term appropriately reflect the nature of the various types of protein-protein interactions studied.     

Optimal clustering for detecting near-native conformations in protein docking

Clustering is one of the most powerful tools in computational biology. The conventional wisdom is that events that occur in clusters are probably not random. In protein docking, the underlying principle is that clustering occurs because long-range electrostatic and/or desolvation forces steer the proteins to a low free-energy attractor at the binding region. Something similar occurs in the docking of small molecules, although in this case shorter-range van der Waals forces play a more critical role. Based on the above, we have developed two different clustering strategies to predict docked conformations based on the clustering properties of a uniform sampling of low free-energy protein-protein and protein-small molecule complexes. We report on significant improvements in the automated prediction and discrimination of docked conformations by using the cluster size and consensus as a ranking criterion. We show that the success of clustering depends on identifying the appropriate clustering radius of the system. The clustering radius for protein-protein complexes is consistent with the range of the electrostatics and desolvation free energies (i.e., between 4 and 9 Angstroms); for protein-small molecule docking, the radius is set by van der Waals interactions (i.e., at approximately 2 Angstroms). Without any a priori information, a simple analysis of the histogram of distance separations between the set of docked conformations can evaluate the clustering properties of the data set. Clustering is observed when the histogram is bimodal. Data clustering is optimal if one chooses the clustering radius to be the minimum after the first peak of the bimodal distribution. We show that using this optimal radius further improves the discrimination of near-native complex structures.     

Hydrophobicity at the surface of proteins

A new method is presented to quantitatively estimate and graphically display the propensity of nonpolar groups to bind at the surface of proteins. It is based on the calculation of the binding energy, i.e., van der Waals interaction plus protein electrostatic desolvation, of a nonpolar probe sphere rolled over the protein surface, and on the color coding of this quantity on a smooth molecular surface (hydrophobicity map). The method is validated on ten protein-ligand complexes and is shown to distinguish precisely where polar and nonpolar groups preferentially bind. Comparisons with existing approaches, like the display of the electrostatic potential or the curvature, illustrate the advantages and the better predictive power of the present method. Hydrophobicity maps will play an important role in the characterization of binding sites for the large number of proteins emerging from the genome projects and structure modeling approaches.     

Investigation of Atomic Level Patterns in Protein—Small Ligand Interactions

Background

Shape complementarity and non-covalent interactions are believed to drive protein-ligand interaction. To date protein-protein, protein-DNA, and protein-RNA interactions were systematically investigated, which is in contrast to interactions with small ligands. We investigate the role of covalent and non-covalent bonds in protein-small ligand interactions using a comprehensive dataset of 2,320 complexes.

Methodology and Principal Findings

We show that protein-ligand interactions are governed by different forces for different ligand types, i.e., protein-organic compound interactions are governed by hydrogen bonds, van der Waals contacts, and covalent bonds; protein-metal ion interactions are dominated by electrostatic force and coordination bonds; protein-anion interactions are established with electrostatic force, hydrogen bonds, and van der Waals contacts; and protein-inorganic cluster interactions are driven by coordination bonds. We extracted several frequently occurring atomic-level patterns concerning these interactions. For instance, 73% of investigated covalent bonds were summarized with just three patterns in which bonds are formed between thiol of Cys and carbon or sulfur atoms of ligands, and nitrogen of Lys and carbon of ligands. Similar patterns were found for the coordination bonds. Hydrogen bonds occur in 67% of protein-organic compound complexes and 66% of them are formed between NH- group of protein residues and oxygen atom of ligands. We quantify relative abundance of specific interaction types and discuss their characteristic features. The extracted protein-organic compound patterns are shown to complement and improve a geometric approach for prediction of binding sites.

Conclusions and Significance

We show that for a given type (group) of ligands and type of the interaction force, majority of protein-ligand interactions are repetitive and could be summarized with several simple atomic-level patterns. We summarize and analyze 10 frequently occurring interaction patterns that cover 56% of all considered complexes and we show a practical application for the patterns that concerns interactions with organic compounds.     

Continuum electrostatic model for the binding of cytochrome c2 to the photosynthetic reaction center from Rhodobacter sphaeroides

Electrostatic interactions are important for protein-protein association. In this study, we examined the electrostatic interactions between two proteins, cytochrome c(2) (cyt c(2)) and the reaction center (RC) from the photosynthetic bacterium Rhodobacter sphaeroides, that function in intermolecular electron transfer in photosynthesis. Electrostatic contributions to the binding energy for the cyt c(2)-RC complex were calculated using continuum electrostatic methods based on the recent cocrystal structure [Axelrod, H. L., et al. (2002) J. Mol. Biol. 319, 501-515]. Calculated changes in binding energy due to mutations of charged interface residues agreed with experimental results for a protein dielectric constant epsilon(in) of 10. However, the electrostatic contribution to the binding energy for the complex was close to zero due to unfavorable desolvation energies that compensate for the favorable Coulomb attraction. The electrostatic energy calculated as a function of displacement of the cyt c(2) from the bound position showed a shallow minimum at a position near but displaced from the cocrystal configuration. These results show that although electrostatic steering is present, other short-range interactions must be present to contribute to the binding energy and to determine the structure of the complex. Calculations made to model the experimental data on association rates indicate a solvent-separated transition state for binding in which the cyt c(2) is displaced approximately 8 A above its position in the bound complex. These results are consistent with a two-step model for protein association: electrostatic docking of the cyt c(2) followed by desolvation to form short-range van der Waals contacts for rapid electron transfer.     

On the role of the crystal environment in determining protein side-chain conformations

The role of crystal packing in determining the observed conformations of amino acid side-chains in protein crystals is investigated by (1) analysis of a database of proteins that have been crystallized in different unit cells (space group or unit cell dimensions) and (2) theoretical predictions of side-chain conformations with the crystal environment explicitly represented. Both of these approaches indicate that the crystal environment plays an important role in determining the conformations of polar side-chains on the surfaces of proteins. Inclusion of the crystal environment permits a more sensitive measurement of the achievable accuracy of side-chain prediction programs, when validating against structures obtained by X-ray crystallography. Our side-chain prediction program uses an all-atom force field and a Generalized Born model of solvation and is thus capable of modeling simple packing effects (i.e. van der Waals interactions), electrostatic effects, and desolvation, which are all important mechanisms by which the crystal environment impacts observed side-chain conformations. Our results are also relevant to the understanding of changes in side-chain conformation that may result from ligand docking and protein-protein association, insofar as the results reveal how side-chain conformations change in response to their local environment.     

Modeling side-chains using molecular dynamics improve recognition of binding region in CAPRI targets

The CAPRI-II experiment added an extra level of complexity to the problem of predicting protein-protein interactions by including 5 targets for which participants had to build or complete the 3-dimensional (3D) structure of either the receptor or ligand based on the structure of a close homolog. In this article, we describe how modeling key side-chains using molecular dynamics (MD) in explicit solvent improved the recognition of the binding region of a free energy- based computational docking method. In particular, we show that MD is able to predict with relatively high accuracy the rotamer conformation of the anchor side-chains important for molecular recognition as suggested by Rajamani et al. (Proc Natl Acad Sci USA 2004;101:11287-11292). As expected, the conformations are some of the most common rotamers for the given residue, while latch side-chains that undergo induced fit upon binding are forced into less common conformations. Using these models as starting conformations in conjunction with the rigid-body docking server ClusPro and the flexible docking algorithm SmoothDock, we produced valuable predictions for 6 of the 9 targets in CAPRI-II, missing only the 3 targets that underwent significant structural rearrangements upon binding. We also show that our free energy- based scoring function, consisting of the sum of van der Waals, Coulombic electrostatic with a distance-dependent dielectric, and desolvation free energy successfully discriminates the nativelike conformation of our submitted predictions. The latter emphasizes the critical role that thermodynamics plays on our methodology, and validates the generality of the algorithm to predict protein interactions.     

Electrostatic contribution to the binding stability of protein-protein complexes

To investigate roles of electrostatic interactions in protein binding stability, electrostatic calculations were carried out on a set of 64 mutations over six protein-protein complexes. These mutations alter polar interactions across the interface and were selected for putative dominance of electrostatic contributions to the binding stability. Three protocols of implementing the Poisson-Boltzmann model were tested. In vdW4 the dielectric boundary between the protein low dielectric and the solvent high dielectric is defined as the protein van der Waals surface and the protein dielectric constant is set to 4. In SE4 and SE20, the dielectric boundary is defined as the surface of the protein interior inaccessible to a 1.4-A solvent probe, and the protein dielectric constant is set to 4 and 20, respectively. In line with earlier studies on the barnase-barstar complex, the vdW4 results on the large set of mutations showed the closest agreement with experimental data. The agreement between vdW4 and experiment supports the contention of dominant electrostatic contributions for the mutations, but their differences also suggest van der Waals and hydrophobic contributions. The results presented here will serve as a guide for future refinement in electrostatic calculation and inclusion of nonelectrostatic effects.     

Comparative binding energy analysis of the substrate specificity of haloalkane dehalogenase from Xanthobacter autotrophicus GJ10

Comparative binding energy (COMBINE) analysis was conducted for 18 substrates of the haloalkane dehalogenase from Xanthobacter autotrophicus GJ10 (DhlA): 1-chlorobutane, 1-chlorohexane, dichloromethane, 1,2-dichloroethane, 1,2-dichloropropane, 2-chloroethanol, epichlorohydrine, 2-chloroacetonitrile, 2-chloroacetamide, and their brominated analogues. The purpose of the COMBINE analysis was to identify the amino acid residues determining the substrate specificity of the haloalkane dehalogenase. This knowledge is essential for the tailoring of this enzyme for biotechnological applications. Complexes of the enzyme with these substrates were modeled and then refined by molecular mechanics energy minimization. The intermolecular enzyme-substrate energy was decomposed into residue-wise van der Waals and electrostatic contributions and complemented by surface area dependent and electrostatic desolvation terms. Partial least-squares projection to latent structures analysis was then used to establish relationships between the energy contributions and the experimental apparent dissociation constants. A model containing van der Waals and electrostatic intermolecular interaction energy contributions calculated using the AMBER force field explained 91% (73% cross-validated) of the quantitative variance in the apparent dissociation constants. A model based on van der Waals intermolecular contributions from AMBER and electrostatic interactions derived from the Poisson-Boltzmann equation explained 93% (74% cross-validated) of the quantitative variance. COMBINE models predicted correctly the change in apparent dissociation constants upon single-point mutation of DhlA for six enzyme-substrate complexes. The amino acid residues contributing most significantly to the substrate specificity of DhlA were identified; they include Asp124, Trp125, Phe164, Phe172, Trp175, Phe222, Pro223, and Leu263. These residues are suitable targets for modification by site-directed mutagenesis.     

What determines the van der Waals coefficient beta in the LIE (linear interaction energy) method to estimate binding free energies using molecular dynamics simulations?

Recently a semiempirical method has been proposed by Aqvist et al. to calculate absolute and relative binding free energies. In this method, the absolute binding free energy of a ligand is estimated as deltaGbind = alpha + beta, where Vel(bound) and Vvdw(bound) are the electrostatic and van der Waals interaction energies between the ligand and the solvated protein from an molecular dynamics (MD) trajectory with ligand bound to protein and Vel(free) and Vel(free) and Vvdw(free) are the electrostatic and van der Waals interaction energies between the ligand and the water from an MD trajectory with the ligand in water. A set of values, alpha = 0.5 and beta = 0.16, was found to give results in good agreement with experimental data. Later, however, different optimal values of beta were found in studies of compounds binding to P450cam and avidin. The present work investigates how the optimal value of beta depends on the nature of binding sites for different protein-ligand interactions. By examining seven ligands interacting with five proteins, we have discovered a linear correlation between the value of beta and the weighted non-polar desolvation ratio (WNDR), with a correlation coefficient of 0.96. We have also examined the ability of this correlation to predict optimal values of beta for different ligands binding to a single protein. We studied twelve neutral compounds bound to avidin. In this case, the WNDR approach gave a better estimate of the absolute binding free energies than results obtained using the fixed value of beta found for biotin-avidin. In terms of reproducing the relative binding free energy to biotin, the fixed-beta value gave better results for compounds similar to biotin, but for compounds less similar to biotin, the WNDR approach led to better relative binding free energies.     

A computational study of the protein-ligand interactions in CDK2 inhibitors: using quantum mechanics/molecular mechanics interaction energy as a predictor of the biological activity

We report a combined quantum mechanics/molecular mechanics (QM/MM) study to determine the protein-ligand interaction energy between CDK2 (cyclin-dependent kinase 2) and five inhibitors with the N(2)-substituted 6-cyclohexyl-methoxy-purine scaffold. The computational results in this work show that the QM/MM interaction energy is strongly correlated to the biological activity and can be used as a predictor, at least within a family of substrates. A detailed analysis of the protein-ligand structures obtained from molecular dynamics simulations shows specific interactions within the active site that, in some cases, have not been reported before to our knowledge. The computed interaction energy gauges the strength of protein-ligand interactions. Finally, energy decomposition and multiple regression analyses were performed to check the contribution of the electrostatic and van der Waals energies to the total interaction energy and to show the capabilities of the computational model to identify new potent inhibitors.     

Predicting protein-protein interactions on a proteome scale by matching evolutionary and structural similarities at interfaces using PRISM

Prediction of protein-protein interactions at the structural level on the proteome scale is important because it allows prediction of protein function, helps drug discovery and takes steps toward genome-wide structural systems biology. We provide a protocol (termed PRISM, protein interactions by structural matching) for large-scale prediction of protein-protein interactions and assembly of protein complex structures. The method consists of two components: rigid-body structural comparisons of target proteins to known template protein-protein interfaces and flexible refinement using a docking energy function. The PRISM rationale follows our observation that globally different protein structures can interact via similar architectural motifs. PRISM predicts binding residues by using structural similarity and evolutionary conservation of putative binding residue 'hot spots'. Ultimately, PRISM could help to construct cellular pathways and functional, proteome-scale annotation. PRISM is implemented in Python and runs in a UNIX environment. The program accepts Protein Data Bank-formatted protein structures and is available at http://prism.ccbb.ku.edu.tr/prism_protocol/.     

Analytically defined surfaces to analyze molecular interaction properties

Molecular surfaces are widely used for characterizing molecules and displaying and quantifying their interaction properties. Here we consider molecular surfaces defined as isocontours of a function (a sum of exponential functions centered on each atom) that approximately represents electron density. The smoothness is advantageous for surface mapping of molecular properties (e.g., electrostatic potential). By varying parameters, these surfaces can be constructed to represent the van der Waals or solvent-accessible surface of a molecular with any accuracy. We describe numerical algorithms to operate on the analytically defined surfaces. Two applications are considered: (1) We define and locate extremal points of molecular properties on the surfaces. The extremal points provide a compact representation of a property on a surface, obviating the necessity to compute values of the property on an array of surface points as is usually done; (2) a molecular surface patch or interface is projected onto a flat surface (by introducing curvilinear coordinates) with approximate conservation of area for analysis purposes. Applications to studies of protein-protein interactions are described.     

Study of protein-protein interaction using conformational space annealing

We apply conformational space annealing (CSA), an efficient global optimization method, to the study of protein-protein interaction. The CSA is incorporated into the Tinker molecular modeling package along with a B-spline method for CAPRI Round 5 experiments. We have used an energy function for the protein-protein interaction that consists of electrostatic interaction, van der Waals interaction, and solvation energy terms represented by the occupancy desolvation method. The parameters of the AMBER94 all-atom empirical force field are used. Each energy term is calculated by precalculated grid potentials and B-spline method approximation. The ligand protein is placed inside a sphere of 50 A radius centered at an appropriate location, and the CSA rigid docking studies are carried out to find stable complexes. Up to 10 complexes are selected using the K-mean clustering method and biological information when available. These complexes are energy-minimized for further refinement by considering the flexibility of interacting proteins. The results show that the CSA method has a potential for the study of protein-protein interaction.     

An all atom energy based computational protocol for predicting binding affinities of protein-ligand complexes

We report here a computationally fast protocol for predicting binding affinities of non-metallo protein-ligand complexes. The protocol builds in an all atom energy based empirical scoring function comprising electrostatics, van der Waals, hydrophobicity and loss of conformational entropy of protein side chains upon ligand binding. The method is designed to ensure transferability across diverse systems and has been validated on a heterogenous dataset of 161 complexes consisting of 55 unique protein targets. The scoring function trained on a dataset of 61 complexes yielded a correlation of r=0.92 for the predicted binding free energies against the experimental binding affinities. Model validation and parameter analysis studies ensure the predictive ability of the scoring function. When tested on the remaining 100 protein-ligand complexes a correlation of r=0.92 was recovered. The high correlation obtained underscores the potential applicability of the methodology in drug design endeavors. The scoring function has been web enabled at as binding affinity prediction of protein-ligand (BAPPL) server.     

Characterization of domain-peptide interaction interface: a case study on the amphiphysin-1 SH3 domain

Many important protein-protein interactions are mediated by peptide recognition modular domains, such as the Src homology 3 (SH3), SH2, PDZ, and WW domains. Characterizing the interaction interface of domain-peptide complexes and predicting binding specificity for modular domains are critical for deciphering protein-protein interaction networks. Here, we propose the use of an energetic decomposition analysis to characterize domain-peptide interactions and the molecular interaction energy components (MIECs), including van der Waals, electrostatic, and desolvation energy between residue pairs on the binding interface. We show a proof-of-concept study on the amphiphysin-1 SH3 domain interacting with its peptide ligands. The structures of the human amphiphysin-1 SH3 domain complexed with 884 peptides were first modeled using virtual mutagenesis and optimized by molecular mechanics (MM) minimization. Next, the MIECs between domain and peptide residues were computed using the MM/generalized Born decomposition analysis. We conducted two types of statistical analyses on the MIECs to demonstrate their usefulness for predicting binding affinities of peptides and for classifying peptides into binder and non-binder categories. First, combining partial least squares analysis and genetic algorithm, we fitted linear regression models between the MIECs and the peptide binding affinities on the training data set. These models were then used to predict binding affinities for peptides in the test data set; the predicted values have a correlation coefficient of 0.81 and an unsigned mean error of 0.39 compared with the experimentally measured ones. The partial least squares-genetic algorithm analysis on the MIECs revealed the critical interactions for the binding specificity of the amphiphysin-1 SH3 domain. Next, a support vector machine (SVM) was employed to build classification models based on the MIECs of peptides in the training set. A rigorous training-validation procedure was used to assess the performances of different kernel functions in SVM and different combinations of the MIECs. The best SVM classifier gave satisfactory predictions for the test set, indicated by average prediction accuracy rates of 78% and 91% for the binding and non-binding peptides, respectively. We also showed that the performance of our approach on both binding affinity prediction and binder/non-binder classification was superior to the performances of the conventional MM/Poisson-Boltzmann solvent-accessible surface area and MM/generalized Born solvent-accessible surface area calculations. Our study demonstrates that the analysis of the MIECs between peptides and the SH3 domain can successfully characterize the binding interface, and it provides a framework to derive integrated prediction models for different domain-peptide systems.     

Understanding beta-hairpin formation by molecular dynamics simulations of unfolding

We have studied the mechanism of formation of a 16-residue beta-hairpin from the protein GB1 using molecular dynamics simulations in an aqueous environment. The analysis of unfolding trajectories at high temperatures suggests a refolding pathway consisting of several transient intermediates. The changes in the interaction energies of residues are related with the structural changes during the unfolding of the hairpin. The electrostatic energies of the residues in the turn region are found to be responsible for the transition between the folded state and the hydrophobic core state. The van der Waals interaction energies of the residues in the hydrophobic core reflect the behavior of the radius of gyration of the core region. We have examined the opposing influences of the protein-protein (PP) energy, which favors the native state, and the protein-solvent (PS) energy, which favors unfolding, in the formation of the beta-hairpin structure. It is found that the behavior of the electrostatic components of PP and PS energies reflects the structural changes associated with the loss of backbone hydrogen bonding. Relative changes in the PP and PS van der Waals interactions are related with the disruption of the hydrophobic core of a protein. The results of the simulations support the hydrophobic collapse mechanism of beta-hairpin folding.     

Alignment of non-covalent interactions at protein-protein interfaces

Background

The study and comparison of protein-protein interfaces is essential for the understanding of the mechanisms of interaction between proteins. While there are many methods for comparing protein structures and protein binding sites, so far no methods have been reported for comparing the geometry of non-covalent interactions occurring at protein-protein interfaces.

Methodology/Principal Findings

Here we present a method for aligning non-covalent interactions between different protein-protein interfaces. The method aligns the vector representations of van der Waals interactions and hydrogen bonds based on their geometry. The method has been applied to a dataset which comprises a variety of protein-protein interfaces. The alignments are consistent to a large extent with the results obtained using two other complementary approaches. In addition, we apply the method to three examples of protein mimicry. The method successfully aligns respective interfaces and allows for recognizing conserved interface regions.

Conclusions/Significance

The Galinter method has been validated in the comparison of interfaces in which homologous subunits are involved, including cases of mimicry. The method is also applicable to comparing interfaces involving non-peptidic compounds. Galinter assists users in identifying local interface regions with similar patterns of non-covalent interactions. This is particularly relevant to the investigation of the molecular basis of interaction mimicry.     

Do electrostatic interactions destabilize protein–nucleic acid binding?

The negatively charged phosphates of nucleic acids are often paired with positively charged residues upon binding proteins. It was thus counter-intuitive when previous Poisson-Boltzmann (PB) calculations gave positive energies from electrostatic interactions, meaning that they destabilize protein-nucleic acid binding. Our own PB calculations on protein-protein binding have shown that the sign and the magnitude of the electrostatic component are sensitive to the specification of the dielectric boundary in PB calculations. A popular choice for the boundary between the solute low dielectric and the solvent high dielectric is the molecular surface; an alternative is the van der Waals (vdW) surface. In line with results for protein-protein binding, in this article, we found that PB calculations with the molecular surface gave positive electrostatic interaction energies for two protein-RNA complexes, but the signs are reversed when the vdW surface was used. Therefore, whether destabilizing or stabilizing effects are predicted depends on the choice of the dielectric boundary. The two calculation protocols, however, yielded similar salt effects on the binding affinity. Effects of charge mutations differentiated the two calculation protocols; PB calculations with the vdW surface had smaller deviations overall from experimental data.