Electrostatic contributions to protein-protein binding affinities: application to Rap/Raf interaction

Analysis of Protein-Protein Interactions for Intermolecular Bond Prediction

Protein-protein interactions often involve a complex system of intermolecular interactions between residues and atoms at the binding site. A comprehensive exploration of these interactions can help reveal key residues involved in protein-protein recognition that are not obvious using other protein analysis techniques. This paper presents and extends DiffBond, a novel method for identifying and classifying intermolecular bonds while applying standard definitions of bonds in chemical literature to explain protein interactions. DiffBond predicted intermolecular bonds from four protein complexes: Barnase-Barstar, Rap1a-raf, SMAD2-SMAD4, and a subset of complexes formed from three-finger toxins and nAChRs. Based on validation through manual literature search and through comparison of two protein complexes from the SKEMPI dataset, DiffBond was able to identify intermolecular ionic bonds and hydrogen bonds with high precision and recall, and identify salt bridges with high precision. DiffBond predictions on bond existence were also strongly correlated with observations of Gibbs free energy change and electrostatic complementarity in mutational experiments. DiffBond can be a powerful tool for predicting and characterizing influential residues in protein-protein interactions, and its predictions can support research in mutational experiments and drug design.

Electrostatic influence on IL-1 transport through the GSDMD pore

A variety of signals, including inflammasome activation, trigger the formation of large transmembrane pores by gasdermin D (GSDMD). There are primarily two functions of the GSDMD pore, to drive lytic cell death, known as pyroptosis, and to permit the release of leaderless interleukin-1 (IL-1) family cytokines, a process that does not require pyroptosis. We are interested in the mechanism by which the GSDMD pore channels IL-1 release from living cells. Recent studies revealed that electrostatic interaction, in addition to cargo size, plays a critical role in GSDMD-dependent protein release. Here, we determined computationally that to enable electrostatic filtering against pro-IL-1β, acidic lipids in the membrane need to effectively neutralize positive charges in the membrane-facing patches of the GSDMD pore. In addition, we predicted that salt has an attenuating effect on electrostatic filtering and then validated this prediction using a liposome leakage assay. A calibrated electrostatic screening factor is necessary to account for the experimental observations, suggesting that ion distribution within the pore may be different from the bulk solution. Our findings corroborate the electrostatic influence of IL-1 transport exerted by the GSDMD pore and reveal extrinsic factors, including lipid and salt, that affect the electrostatic environment.

DiffBond: A Method for Predicting Intermolecular Bond Formation

Many tools that explore models of protein complexes are also able to analyze interactions between specific residues and atoms. A comprehensive exploration of these interactions can often uncover aspects of protein-protein recognition that are not obvious using other protein analysis techniques. This paper describes DiffBond, a novel method for searching for intermolecular interactions between protein complexes while differentiating between three different types of interaction: hydrogen bonds, ionic bonds, and salt bridges. DiffBond incorporates textbook definitions of these three interactions while contending with uncertainties that are inherent in computational models of interacting proteins. We used it to examine the barnase-barstar, Rap1a-raf, and Smad2-Smad4 complexes, as well as a subset of protein complexes formed between three-finger toxins and nAChRs. Based on electrostatic interactions established by previous experimental studies, DiffBond was able to identify ionic and hydrogen bonds with high precision and recall, and identify salt bridges with high precision. In combination with other electrostatic analysis methods, DiffBond can be a useful tool in helping predict influential amino acids in protein-protein interactions and characterizing the type of interaction.

Predicting Mutational Effects on Receptor Binding of the Spike Protein of SARS-CoV-2 Variants

The pandemic caused by SARS-CoV-2 has cost millions of lives and tremendous social/financial loss. The virus continues to evolve and mutate. In particular, the recently emerged "UK", "South Africa", and Delta variants show higher infectivity and spreading speed. Thus, the relationship between the mutations of certain amino acids and the spreading speed of the virus is a problem of great importance. In this respect, understanding the mutational mechanism is crucial for surveillance and prediction of future mutations as well as antibody/vaccine development. In this work, we used a coarse-grained model (that was used previously in predicting the importance of mutations of N501) to calculate the free energy change of various types of single-site or combined-site mutations. This was done for the UK, South Africa, and Delta mutants. We investigated the underlying mechanisms of the binding affinity changes for mutations at different spike protein domains of SARS-CoV-2 and provided the energy basis for the resistance of the E484 mutant to the antibody m396. Other potential mutation sites were also predicted. Furthermore, the in silico predictions were assessed by functional experiments. The results establish that the faster spreading of recently observed mutants is strongly correlated with the binding-affinity enhancement between virus and human receptor as well as with the reduction of the binding to the m396 antibody. Significantly, the current approach offers a way to predict new variants and to assess the effectiveness of different antibodies toward such variants.

Fragment-based quantum mechanical calculation of protein-protein binding affinities

The electrostatically embedded generalized molecular fractionation with conjugate caps (EE-GMFCC) method has been successfully utilized for efficient linear-scaling quantum mechanical (QM) calculation of protein energies. In this work, we applied the EE-GMFCC method for calculation of binding affinity of Endonuclease colicin-immunity protein complex. The binding free energy changes between the wild-type and mutants of the complex calculated by EE-GMFCC are in good agreement with experimental results. The correlation coefficient (R) between the predicted binding energy changes and experimental values is 0.906 at the B3LYP/6-31G*-D level, based on the snapshot whose binding affinity is closest to the average result from the molecular mechanics/Poisson-Boltzmann surface area (MM/PBSA) calculation. The inclusion of the QM effects is important for accurate prediction of protein-protein binding affinities. Moreover, the self-consistent calculation of PB solvation energy is required for accurate calculations of protein-protein binding free energies. This study demonstrates that the EE-GMFCC method is capable of providing reliable prediction of relative binding affinities for protein-protein complexes. © 2018 Wiley Periodicals, Inc.

Trypsin-Ligand binding affinities calculated using an effective interaction entropy method under polarized force field

Molecular dynamics (MD) simulation in the explicit water is performed to study the interaction mechanism of trypsin-ligand binding under the AMBER force field and polarized protein-specific charge (PPC) force field combined the new developed highly efficient interaction entropy (IE) method for calculation of entropy change. And the detailed analysis and comparison of the results of MD simulation for two trypsin-ligand systems show that the root-mean-square deviation (RMSD) of backbone atoms, B-factor, intra-protein and protein-ligand hydrogen bonds are more stable under PPC force field than AMBER force field. Our results demonstrate that the IE method is superior than the traditional normal mode (Nmode) method in the calculation of entropy change and the calculated binding free energy under the PPC force field combined with the IE method is more close to the experimental value than other three combinations (AMBER-Nmode, AMBER-IE and PPC-Nmode). And three critical hydrogen bonds between trypsin and ligand are broken under AMBER force field. However, they are well preserved under PPC force field. Detailed binding interactions of ligands with trypsin are further analyzed. The present work demonstrates that the polarized force field combined the highly efficient IE method is critical in MD simulation and free energy calculation.

Simulating the catalytic effect of a designed mononuclear zinc metalloenzyme that catalyzes the hydrolysis of phosphate triesters

One of the greatest challenges in biotechnology and in biochemistry is the ability to design efficient enzymes. In fact, such an ability would be one of the most convincing manifestations of a full understanding of the origin of enzyme catalysis. Despite some progress on this front, most of the advances have been made by placing the reacting fragments in the proper places rather than by optimizing the preorganization of the environment, which is the key factor in enzyme catalysis. A rational improvement of the preorganization and a consistent assessment of the effectiveness of different design options require approaches capable of evaluating reliably the actual catalytic effect. In this work we examine the ability of the empirical valence bond (EVB) to reproduce the results of directed evolution improvements of the catalysis of diethyl 7-hydroxycoumarinyl by a designed mononuclear zinc metalloenzyme. Encouragingly, our study reproduced the catalytic effect obtained by directed evolution and offers a good start for further studies of this system.

Multiscale modeling of biological functions: from enzymes to molecular machines (Nobel Lecture)

A detailed understanding of the action of biological molecules is a pre-requisite for rational advances in health sciences and related fields. Here, the challenge is to move from available structural information to a clear understanding of the underlying function of the system. In light of the complexity of macromolecular complexes, it is essential to use computer simulations to describe how the molecular forces are related to a given function. However, using a full and reliable quantum mechanical representation of large molecular systems has been practically impossible. The solution to this (and related) problems has emerged from the realization that large systems can be spatially divided into a region where the quantum mechanical description is essential (e.g. a region where bonds are being broken), with the remainder of the system being represented on a simpler level by empirical force fields. This idea has been particularly effective in the development of the combined quantum mechanics/molecular mechanics (QM/MM) models. Here, the coupling between the electrostatic effects of the quantum and classical subsystems has been a key to the advances in describing the functions of enzymes and other biological molecules. The same idea of representing complex systems in different resolutions in both time and length scales has been found to be very useful in modeling the action of complex systems. In such cases, starting with coarse grained (CG) representations that were originally found to be very useful in simulating protein folding, and augmenting them with a focus on electrostatic energies, has led to models that are particularly effective in probing the action of molecular machines. The same multiscale idea is likely to play a major role in modeling of even more complex systems, including cells and collections of cells.

Predicting the pKa and stability of organic acids and bases at an oil-water interface

We have used density functional theory and the implicit solvent model, COSMO-RS, to investigate how the acidity constant, pKa, of organic acids and bases adsorbed at the organic compound-aqueous solution interface changes, compared to its value in the aqueous phase. The pKa determine the surface charge density of the molecules that accumulate at the fluid-fluid interface. We have estimated the pKa by comparing the stability of the protonated and unprotonated forms of a series of molecules in the bulk aqueous solution and at an interface where parts of each molecule reside in the hydrophobic phase and the rest remains in the hydrophilic phase. We found that the pKa for acids is shifted by ∼1 pH unit to higher values compared to the bulk water pKa, whereas they are shifted to lower values by a similar amount for bases. Because this pKa shift is similar in magnitude for each of the molecules studied, we propose that the pKa for molecules at a water-organic compound interface can easily be predicted by adding a small shift to the aqueous pKa. This shift is general and correlates with the functional group. We also found that the relative composition of molecules at the fluid-fluid interface is not the same as in the bulk. For example, species such as carboxylic acids are enriched at the interface, where they can dominate surface properties, even when they are a modest component in the bulk fluid. For high surface concentrations of carboxylic acid groups at an interface, such as a self-assembled monolayer, we have demonstrated that the pKa depends on the degree of deprotonation through direct hydrogen bonding between protonated and deprotonated acidic headgroups.

Predicting pKa for proteins using COSMO-RS

We have used the COSMO-RS implicit solvation method to calculate the equilibrium constants, pK_a, for deprotonation of the acidic residues of the ovomucoid inhibitor protein, OMTKY3. The root mean square error for comparison with experimental data is only 0.5 pH units and the maximum error 0.8 pH units. The results show that the accuracy of pK_a prediction using COSMO-RS is as good for large biomolecules as it is for smaller inorganic and organic acids and that the method compares very well to previous pK_a predictions of the OMTKY3 protein using Quantum Mechanics/Molecular Mechanics. Our approach works well for systems of about 1000 atoms or less, which makes it useful for small proteins as well as for investigating portions of larger proteins such as active sites in enzymes.

Modeling catalytic promiscuity in the alkaline phosphatase superfamily

In recent years, it has become increasingly clear that promiscuity plays a key role in the evolution of new enzyme function. This finding has helped to elucidate fundamental aspects of molecular evolution. While there has been extensive experimental work on enzyme promiscuity, computational modeling of the chemical details of such promiscuity has traditionally fallen behind the advances in experimental studies, not least due to the nearly prohibitive computational cost involved in examining multiple substrates with multiple potential mechanisms and binding modes in atomic detail with a reasonable degree of accuracy. However, recent advances in both computational methodologies and power have allowed us to reach a stage in the field where we can start to overcome this problem, and molecular simulations can now provide accurate and efficient descriptions of complex biological systems with substantially less computational cost. This has led to significant advances in our understanding of enzyme function and evolution in a broader sense. Here, we will discuss currently available computational approaches that can allow us to probe the underlying molecular basis for enzyme specificity and selectivity, discussing the inherent strengths and weaknesses of each approach. As a case study, we will discuss recent computational work on different members of the alkaline phosphatase superfamily (AP) using a range of different approaches, showing the complementary insights they have provided. We have selected this particular superfamily, as it poses a number of significant challenges for theory, ranging from the complexity of the actual reaction mechanisms involved to the reliable modeling of the catalytic metal centers, as well as the very large system sizes. We will demonstrate that, through current advances in methodologies, computational tools can provide significant insight into the molecular basis for catalytic promiscuity, and, therefore, in turn, the mechanisms of protein functional evolution.

A multiscale coarse-grained polarizable solvent model for handling long tail bulk electrostatics

A multiscale coarse-grained approach able to handle efficiently the solvation of microscopic solutes in extended chemical environment is described. That approach is able to compute readily and efficiently very long-range solute/solvent electrostatic microscopic interactions, up to the 1-μm scale, by considering a reduced amount of computational resources. All the required parameters are assigned to reproduce available data concerning the solvation of single ions. Such a strategy makes it possible to reproduce with good accuracy the solvation properties concerning simple ion pairs in solution (in particular, the asymptotic behavior of the ion pair potentials of mean force). This new method represents an extension of the polarizable pseudoparticle solvent model, which has been recently improved to account for the main features of hydrophobic effects in liquid water (Masella et al., J. Comput. Chem. 2011, 32, 2664). This multiscale approach is well suited to be used for computing the impact of charge changes in free energy computations, in terms of both accuracy and efficiency.

Why nature really chose phosphate

Phosphoryl transfer plays key roles in signaling, energy transduction, protein synthesis, and maintaining the integrity of the genetic material. On the surface, it would appear to be a simple nucleophile displacement reaction. However, this simplicity is deceptive, as, even in aqueous solution, the low-lying d-orbitals on the phosphorus atom allow for eight distinct mechanistic possibilities, before even introducing the complexities of the enzyme catalyzed reactions. To further complicate matters, while powerful, traditional experimental techniques such as the use of linear free-energy relationships (LFER) or measuring isotope effects cannot make unique distinctions between different potential mechanisms. A quarter of a century has passed since Westheimer wrote his seminal review, 'Why Nature Chose Phosphate' (Science 235 (1987), 1173), and a lot has changed in the field since then. The present review revisits this biologically crucial issue, exploring both relevant enzymatic systems as well as the corresponding chemistry in aqueous solution, and demonstrating that the only way key questions in this field are likely to be resolved is through careful theoretical studies (which of course should be able to reproduce all relevant experimental data). Finally, we demonstrate that the reason that nature really chose phosphate is due to interplay between two counteracting effects: on the one hand, phosphates are negatively charged and the resulting charge-charge repulsion with the attacking nucleophile contributes to the very high barrier for hydrolysis, making phosphate esters among the most inert compounds known. However, biology is not only about reducing the barrier to unfavorable chemical reactions. That is, the same charge-charge repulsion that makes phosphate ester hydrolysis so unfavorable also makes it possible to regulate, by exploiting the electrostatics. This means that phosphate ester hydrolysis can not only be turned on, but also be turned off, by fine tuning the electrostatic environment and the present review demonstrates numerous examples where this is the case. Without this capacity for regulation, it would be impossible to have for instance a signaling or metabolic cascade, where the action of each participant is determined by the fine-tuned activity of the previous piece in the production line. This makes phosphate esters the ideal compounds to facilitate life as we know it.

Effect of interprotein polarization on protein-protein binding energy

Molecular dynamics simulation in explicit water for the binding of the benchmark barnase-barstar complex was carried out to investigate the effect polarization of interprotein hydrogen bonds on its binding free energy. Our study is based on the AMBER force field but with polarized atomic charges derived from fragment quantum mechanical calculation for the protein complex. The quantum-derived atomic charges include the effect of polarization of interprotein hydrogen bonds, which was absent in the standard force fields that were used in previous theoretical calculations of barnase-barstar binding energy. This study shows that this polarization effect impacts both the static (electronic) and dynamic interprotein electrostatic interactions and significantly lowers the free energy of the barnase-barstar complex.

Validating the vitality strategy for fighting drug resistance

The current challenge in designing effective drugs against HIV-1 is to find novel candidates with high potency, but with a lower susceptibility to mutations associated with drug resistance. Trying to address this challenge, we developed in our previous study (Ishikita and Warshel, Angew Chem Int Ed Engl 2008; 47:697-700) a novel computational strategy for fighting drug resistance by predicting the likely moves of the virus through constraints on binding and catalysis. This has been based on calculating the ratio between the vitality values ((K(i) k(cat)/K(M))(mutant)/(K(i) k(cat)/K(M))(wild-type)) and using it as a guide for predicting the moves of the virus. The corresponding calculations of the binding affinity, K(i), were carried out using the semi-macroscopic version of the protein dipole Langevin dipole (PDLD/S) in its linear response approximation (LRA) in its β version (PDLD/S-LRA/β). We also calculate the proteolytic efficiency, k(cat)/K(M), by evaluating the transition state (TS) binding free energies using the PDLD/S-LRA/β method. Here we provide an extensive validation of our strategy by calculating the vitality of six existing clinical and experimental drug candidates. It is found that the computationally determined vitalities correlate reasonably well with those derived from the corresponding experimental data. This indicates that the calculated vitality may be used to identify mutations that would be most effective for the survival of the virus. Thus, it should be possible to use our approach in screening for mutations that would provide the most effective resistance to any proposed antiviral drug. This ability should be very useful in guiding the design of drug molecules that will lead to the slowest resistance.

Towards quantitative computer-aided studies of enzymatic enantioselectivity: the case of Candida antarctica lipase A

The prospect for consistent computer-aided refinement of stereoselective enzymes is explored by simulating the hydrolysis of enantiomers of an α-substituted ester by wild-type and mutant Candida antarctica lipase A, using several strategies. In particular, we focused on the use of the empirical valence bond (EVB) method in a quantitative screening for enantioselectivity, and evaluate both k(cat) and k(cat)/K(M) of the R and S stereoisomers. We found that an extensive sampling is essential for obtaining converging results. This requirement points towards possible problems with approaches that use a limited conformational sampling. However, performing the proper sampling appears to give encouraging results and to offer a powerful tool for the computer-aided design of enantioselective enzymes. We also explore faster strategies for identifying mutations that will help in augmenting directed-evolution experiments, but these approaches require further refinement.

Probing the structural and energetic basis of kinesin-microtubule binding using computational alanine-scanning mutagenesis

Kinesin-microtubule (MT) binding plays a critical role in facilitating and regulating the motor function of kinesins. To obtain a detailed structural and energetic picture of kinesin-MT binding, we performed large-scale computational alanine-scanning mutagenesis based on long-time molecular dynamics (MD) simulations of the kinesin-MT complex in both ADP and ATP states. First, we built three all-atom kinesin-MT models: human conventional kinesin bound to ADP and mouse KIF1A bound to ADP and ATP. Then, we performed 30 ns MD simulations followed by kinesin-MT binding free energy calculations for both the wild type and mutants obtained after substitution of each charged residue of kinesin with alanine. We found that the kinesin-MT binding free energy is dominated by van der Waals interactions and further enhanced by electrostatic interactions. The calculated mutational changes in kinesin-MT binding free energy are in excellent agreement with results of an experimental alanine-scanning study with a root-mean-square error of ~0.32 kcal/mol [Woehlke, G., et al. (1997) Cell 90, 207-216]. We identified a set of important charged residues involved in the tuning of kinesin-MT binding, which are clustered on several secondary structural elements of kinesin (including well-studied loops L7, L8, L11, and L12, helices α4, α5, and α6, and less-explored loop L2). In particular, we found several key residues that make different contributions to kinesin-MT binding in ADP and ATP states. The mutations of these residues are predicted to fine-tune the motility of kinesin by modulating the conformational transition between the ADP state and the ATP state of kinesin.

Exploring challenges in rational enzyme design by simulating the catalysis in artificial kemp eliminase

One of the fundamental challenges in biotechnology and in biochemistry is the ability to design effective enzymes. Doing so would be a convincing manifestation of a full understanding of the origin of enzyme catalysis. Despite an impressive progress, most of the advances on this front have been made by placing the reacting fragments in the proper places, rather than by optimizing the environment preorganization, which is the key factor in enzyme catalysis. Rational improvement of the preorganization would require approaches capable of evaluating reliably the actual catalytic effect. This work takes previously designed kemp eliminases as a benchmark for a computer aided enzyme design, using the empirical valence bond as the main screening tool. The observed absolute catalytic effect and the effect of directed evolution are reproduced and analyzed (assuming that the substrate is in the designed site). It is found that, in the case of kemp eliminases, the transition state charge distribution makes it hard to exploit the active site polarity, even with the ability to quantify the effect of different mutations. Unexpectedly, it is found that the directed evolution mutants lead to the reduction of solvation of the reactant state by water molecules rather that to the more common mode of transition state stabilization used by naturally evolved enzymes. Finally it is pointed out that our difficulties in improving Kemp eliminase are not due to overlooking exotic effect, but to the challenge in designing a preorganized environment that would exploit the small change it charge distribution during the formation of the transition state.

Absolute binding free energy calculations: on the accuracy of computational scoring of protein-ligand interactions

Calculating the absolute binding free energies is a challenging task. Reliable estimates of binding free energies should provide a guide for rational drug design. It should also provide us with deeper understanding of the correlation between protein structure and its function. Further applications may include identifying novel molecular scaffolds and optimizing lead compounds in computer-aided drug design. Available options to evaluate the absolute binding free energies range from the rigorous but expensive free energy perturbation to the microscopic linear response approximation (LRA/beta version) and related approaches including the linear interaction energy (LIE) to the more approximated and considerably faster scaled protein dipoles Langevin dipoles (PDLD/S-LRA version) as well as the less rigorous molecular mechanics Poisson-Boltzmann/surface area (MM/PBSA) and generalized born/surface area (MM/GBSA) to the less accurate scoring functions. There is a need for an assessment of the performance of different approaches in terms of computer time and reliability. We present a comparative study of the LRA/beta, the LIE, the PDLD/S-LRA/beta, and the more widely used MM/PBSA and assess their abilities to estimate the absolute binding energies. The LRA and LIE methods perform reasonably well but require specialized parameterization for the nonelectrostatic term. The PDLD/S-LRA/beta performs effectively without the need of reparameterization. Our assessment of the MM/PBSA is less optimistic. This approach appears to provide erroneous estimates of the absolute binding energies because of its incorrect entropies and the problematic treatment of electrostatic energies. Overall, the PDLD/S-LRA/beta appears to offer an appealing option for the final stages of massive screening approaches.

Toward accurate screening in computer-aided enzyme design

The ability to design effective enzymes is one of the most fundamental challenges in biotechnology and in some respects in biochemistry. In fact, such ability would be one of the most convincing manifestations of a full understanding of the origin of enzyme catalysis. In this work, we explore the reliability of different simulation approaches, in terms of their ability to rank different possible active site constructs. This validation is done by comparing the ability of different approaches to evaluate the catalytic contributions of various residues in chorismate mutase. It is demonstrated that the empirical valence bond (EVB) model can serve as a practical yet accurate tool in the final stages of computer-aided enzyme design (CAED). Other approaches for fast screening are also examined and found to be less accurate and mainly useful for qualitative screening of ionized residues. It is pointed out that accurate ranking of different options for enzyme design cannot be accomplished by approaches that cannot capture the electrostatic preorganization effect. This is in particular true with regard to current design approaches that use gas phase or small cluster calculations and then estimate the interaction between the enzyme and the transition state (TS) model rather than the TS binding free energy or the relevant activation free energy. The ability of the EVB model to provide a tool for quantitative ranking in the final stage of CAED may help in progressing toward the design of enzymes whose catalytic power is closer to that of native enzymes than to that of the current generation of designer enzymes.

Insights into affinity and specificity in the complexes of alpha-lytic protease and its inhibitor proteins: binding free energy from molecular dynamics simulation

We report the binding free energy calculation and its decomposition for the complexes of alpha-lytic protease and its protein inhibitors using molecular dynamics simulation. Standard mechanism serine protease inhibitors eglin C and OMTKY3 are known to have strong binding affinity for many serine proteases. Their binding loops have significant similarities, including a common P1 Leu as the main anchor in the binding interface. However, recent experiments demonstrate that the two inhibitors have vastly different affinity towards alpha-lytic protease (ALP), a bacterial serine protease. OMTKY3 inhibits the enzyme much more weakly (by approximately 10(6) times) than eglin C. Moreover, a variant of OMTKY3 with five mutations, OMTKY3M, has been shown to inhibit 10(4) times more strongly than the wild-type inhibitor. The underlying mechanisms for the unusually large difference in binding affinities and the effect of mutation are not well understood. Here we use molecular dynamics simulation with molecular mechanics-Poisson Boltzmann/surface area method (MM-PB/SA) to investigate quantitatively the binding specificity. The calculated absolute binding free energies correctly differentiate the thermodynamic stabilities of these protein complexes, but the magnitudes of the binding affinities are systematically overestimated. Analysis of the binding free energy components provides insights into the molecular mechanism of binding specificity. The large DeltaDeltaG(bind) between eglin C and wild type OMTKY3 towards ALP is mainly attributable to the stronger nonpolar interactions in the ALP-eglin C complex, arising from a higher degree of structural complementarity. Here the electrostatic interaction contributes to a lesser extent. The enhanced inhibition in the penta-mutant OMTKY3M over its wild type is entirely due to an overall improvement in the solvent-mediated electrostatic interactions in the ALP-OMTKY3M complex. The results suggest that for these protein-complexes and similar enzyme-inhibitor systems (1) the binding is driven by nonpolar interactions, opposed by overall electrostatic and solute entropy contributions; (2) binding specificity can be tuned by improving the complementarity in electrostatics between two associating proteins. Binding free energy decomposition into contributions from individual protein residues provides additional detailed information on the structural determinants and subtle conformational changes responsible for the binding specificity.

On the electrostatic component of protein-protein binding free energy

Calculations of electrostatic properties of protein-protein complexes are usually done within framework of a model with a certain set of parameters. In this paper we present a comprehensive statistical analysis of the sensitivity of the electrostatic component of binding free energy (DeltaDeltaGel) with respect with different force fields (Charmm, Amber, and OPLS), different values of the internal dielectric constant, and different presentations of molecular surface (different values of the probe radius). The study was done using the largest so far set of entries comprising 260 hetero and 2148 homo protein-protein complexes extracted from a previously developed database of protein complexes (ProtCom). To test the sensitivity of the energy calculations with respect to the structural details, all structures were energy minimized with corresponding force field, and the energies were recalculated. The results indicate that the absolute value of the electrostatic component of the binding free energy (DeltaDeltaGel) is very sensitive to the force field parameters, the minimization procedure, the values of the internal dielectric constant, and the probe radius. Nevertheless our results indicate that certain trends in DeltaDeltaGel behavior are much less sensitive to the calculation parameters. For instance, the fraction of the homo-complexes, for which the electrostatics was found to oppose binding, is 80% regardless of the force fields and parameters used. For the hetero-complexes, however, the percentage of the cases for which electrostatics opposed binding varied from 43% to 85%, depending on the protocol and parameters employed. A significant correlation was found between the effects caused by raising the internal dielectric constant and decreasing the probe radius. Correlations were also found among the results obtained with different force fields. However, despite of the correlations found, the absolute DeltaDeltaGel calculated with different force field parameters could differ more than tens of kcal/mol in some cases. Set of rules of obtaining confident predictions of absolute DeltaDeltaGel and DeltaDeltaGel sign are provided in the conclusion section.PACS codes: 87.15.A-, 87.15. km.

Combining a polarizable force-field and a coarse-grained polarizable solvent model: application to long dynamics simulations of bovine pancreatic trypsin inhibitor

The dynamic coupling between a polarizable protein force field and a particle-based implicit solvent model is described. The polarizable force field, TCPEp, developed recently to simulate protein systems, is characterized by a reduced number of polarizable sites, with a substantial gain in efficiency for an equal chemical accuracy. The Polarizable Pseudo-Particle (PPP) solvent model represents the macroscopic solvent polarization by induced dipoles placed on mobile Lennard-Jones pseudo-particles. The solvent-induced dipoles are sensitive to the solute electric field, but not to each other, so that the computational cost of solvent-solvent interactions is basically negligible. The solute and solvent induced dipoles are determined self-consistently and the equations of motion are solved using an efficient iterative multiple time step procedure. The solvation cost with respect to vacuum simulations is shown to decrease with solute size: the estimated multiplicative factor is 2.5 for a protein containing about 1000 atoms, and as low as 1.15 for 8000 atoms. The model is tested for six 20 ns molecular dynamics trajectories of a traditional benchmark system: the hydrated Bovine Pancreatic Trypsin Inhibitor (BPTI). Even though the TCPEp parameters have not been refined to be used with the solvent PPP model, we observe a good conservation of the BPTI structure along the trajectories. Moreover, our approach is able to provide a description of the protein solvation thermodynamic at the same accuracy as the standard Poisson-Boltzman continuum methods. It provides in addition a good description of the microscopic structural aspects concerning the solute/solvent interaction.

Molecular dynamics-solvated interaction energy studies of protein-protein interactions: the MP1-p14 scaffolding complex

Using the MP1-p14 scaffolding complex from the mitogen-activated protein kinase signaling pathway as model system, we explored a structure-based computational protocol to probe and characterize binding affinity hot spots at protein-protein interfaces. Hot spots are located by virtual alanine-scanning consensus predictions over three different energy functions and two different single-structure representations of the complex. Refined binding affinity predictions for select hot-spot mutations are carried out by applying first-principle methods such as the molecular mechanics generalized Born surface area (MM-GBSA) and solvated interaction energy (SIE) to the molecular dynamics (MD) trajectories for mutated and wild-type complexes. Here, predicted hot-spot residues were actually mutated to alanine, and crystal structures of the mutated complexes were determined. Two mutated MP1-p14 complexes were investigated, the p14(Y56A)-mutated complex and the MP1(L63A,L65A)-mutated complex. Alternative ways to generate MD ensembles for mutant complexes, not relying on crystal structures for mutated complexes, were also investigated. The SIE function, fitted on protein-ligand binding affinities, gave absolute binding affinity predictions in excellent agreement with experiment and outperformed standard MM-GBSA predictions when tested on the MD ensembles of Ras-Raf and Ras-RalGDS protein-protein complexes. For wild-type and mutant MP1-p14 complexes, SIE predictions of relative binding affinities were supported by a yeast two-hybrid assay that provided semiquantitative relative interaction strengths. Results on the MP1-mutated complex suggested that SIE predictions deteriorate if mutant MD ensembles are approximated by just mutating the wild-type MD trajectory. The SIE data on the p14-mutated complex indicated feasibility for generating mutant MD ensembles from mutated wild-type crystal structure, despite local structural differences observed upon mutation. For energetic considerations, this would circumvent costly needs to produce and crystallize mutated complexes. The sensitized protein-protein interface afforded by the p14(Y56A) mutation identified here has practical applications in screening-based discovery of first-generation small-molecule hits for further development into specific modulators of the mitogen-activated protein kinase signaling pathway.

Protein-protein binding is often associated with changes in protonation state

pK(a) values of ionizable residues have been calculated using the PROPKA method and structures of 75 protein-protein complexes and their corresponding free forms. These pK(a) values were used to compute changes in protonation state of individual residues, net changes in protonation state of the complex relative to the uncomplexed proteins, and the correction to a binding energy calculated assuming standard protonation states at pH 7. For each complex, two different structures for the uncomplexed form of the proteins were used: the X-ray structures determined for the proteins in the absence of the other protein and the individual protein structures taken from the structure of the complex (referred to as unbound and bound structures, respectively). In 28 and 77% of the cases considered here, protein-protein binding is accompanied by a complete (>95%) or significant (>50%) change in protonation state of at least one residue using unbound structures. Furthermore, in 36 and 61% of the cases, protein-protein binding is accompanied by a complete or significant net change in protonation state of the complex relative to the separated monomers. Using bound structures, the corresponding values are 12, 51, 20, and 48%. Comparison to experimental data suggest that using unbound and bound structures lead to over- and underestimation of binding-induced protonation state changes, respectively. Thus, we conclude that protein-protein binding is often associated with changes in protonation state of amino acid residues and with changes in the net protonation state of the proteins. The pH-dependent correction to the binding energy contributes at least one order of magnitude to the binding constant in 45 and 23%, using unbound and bound structures, respectively.

Exploring the role of large conformational changes in the fidelity of DNA polymerase beta

The relationships between the conformational landscape, nucleotide insertion catalysis and fidelity of DNA polymerase beta are explored by means of computational simulations. The simulations indicate that the transition states for incorporation of right (R) and wrong (W) nucleotides reside in substantially different protein conformations. The protein conformational changes that reproduce the experimentally observed fidelity are significantly larger than the small rearrangements that usually accompany motions from the reactant state to the transition state in common enzymatic reactions. Once substrate binding has occurred, different constraints imposed on the transition states for insertion of R and W nucleotides render it highly unlikely that both transition states can occur in the same closed structure, because the predicted fidelity would then be many orders of magnitude too large. Since the conformational changes reduce the transition state energy of W incorporation drastically they decrease fidelity rather than increase it. Overall, a better agreement with experimental data is attained when the R is incorporated through a transition state in a closed conformation and W is incorporated through a transition state in one or perhaps several partially open conformations. The generation of free energy surfaces for R and W also allow us to analyze proposals about the relationship between induced fit and fidelity.

Optimization of electrostatic interactions in protein-protein complexes

In this article, we present a statistical analysis of the electrostatic properties of 298 protein-protein complexes and 356 domain-domain structures extracted from the previously developed database of protein complexes (ProtCom, http://www.ces.clemson.edu/compbio/protcom). For each structure in the dataset we calculated the total electrostatic energy of the binding and its two components, Coulombic and reaction field energy. It was found that in a vast majority of the cases (>90%), the total electrostatic component of the binding energy was unfavorable. At the same time, the Coulombic component of the binding energy was found to favor the complex formation while the reaction field component of the binding energy opposed the binding. It was also demonstrated that the components in a wild-type (WT) structure are optimized/anti-optimized with respect to the corresponding distributions, arising from random shuffling of the charged side chains. The degree of this optimization was assessed through the Z-score of WT energy in respect to the random distribution. It was found that the Z-scores of Coulombic interactions peak at a considerably negative value for all 654 cases considered while the Z-score of the reaction field energy varied among different types of complexes. All these findings indicate that the Coulombic interactions within WT protein-protein complexes are optimized to favor the complex formation while the total electrostatic energy predominantly opposes the binding. This observation was used to discriminate WT structures among sets of structural decoys and showed that the electrostatic component of the binding energy is not a good discriminator of the WT; while, Coulombic or reaction field energies perform better depending upon the decoy set used.

COMBINE analysis of the specificity of binding of Ras proteins to their effectors

The small guanosine triphosphate (GTP)-binding proteins of the Ras family are involved in many cellular pathways leading to cell growth, differentiation, and apoptosis. Understanding the interaction of Ras with other proteins is of importance not only for studying signalling mechanisms but also, because of their medical relevance as targets, for anticancer therapy. To study their selectivity and specificity, which are essential to their signal transfer function, we performed COMparative BINding Energy (COMBINE) analysis for 122 different wild-type and mutant complexes between the Ras proteins, Ras and Rap, and their effectors, Raf and RalGDS. The COMBINE models highlighted the amino acid residues responsible for subtle differences in binding of the same effector to the two different Ras proteins, as well as more significant differences in the binding of the two different effectors (RalGDS and Raf) to Ras. The study revealed that E37, D38, and D57 in Ras are nonspecific hot spots at its effector interface, important for stabilization of both the RalGDS-Ras and Raf-Ras complexes. The electrostatic interaction between a GTP analogue and the effector, either Raf or RalGDS, also stabilizes these complexes. The Raf-Ras complexes are specifically stabilized by S39, Y40, and D54, and RalGDS-Ras complexes by E31 and D33. Binding of a small molecule in the vicinity of one of these groups of amino acid residues could increase discrimination between the Raf-Ras and RalGDS-Ras complexes. Despite the different size of the RalGDS-Ras and Raf-Ras complexes, we succeeded in building COMBINE models for one type of complex that were also predictive for the other type of protein complex. Further, using system-specific models trained with only five complexes selected according to the results of principal component analysis, we were able to predict binding affinities for the other mutants of the particular Ras-effector complex. As the COMBINE analysis method is able to explicitly reveal the amino acid residues that have most influence on binding affinity, it is a valuable aid for protein design.

Quantification and rationalization of the higher affinity of sodium over potassium to protein surfaces

For a series of different proteins, including a structural protein, enzyme, inhibitor, protein marker, and a charge-transfer system, we have quantified the higher affinity of Na+ over K+ to the protein surface by means of molecular dynamics simulations and conductivity measurements. Both approaches show that sodium binds at least twice as strongly to the protein surface than potassium does with this effect being present in all proteins under study. Different parts of the protein exterior are responsible to a varying degree for the higher surface affinity of sodium, with the charged carboxylic groups of aspartate and glutamate playing the most important role. Therefore, local ion pairing is the key to the surface preference of sodium over potassium, which is further demonstrated and quantified by simulations of glutamate and aspartate in the form of isolated amino acids as well as short oligopeptides. As a matter of fact, the effect is already present at the level of preferential pairing of the smallest carboxylate anions, formate or acetate, with Na+ versus K+, as shown by molecular dynamics and ab initio quantum chemical calculations. By quantifying and rationalizing the higher preference of sodium over potassium to protein surfaces, the present study opens a way to molecular understanding of many ion-specific (Hofmeister) phenomena involving protein interactions in salt solutions.

Interaction energy decomposition in protein-protein association: a quantum mechanical study of barnase-barstar complex

Protein-protein interactions are very important in the function of a cell. Computational studies of these interactions have been of interest, but often they have utilized classical modelling techniques. In recent years, quantum mechanical (QM) treatment of entire proteins has emerged as a powerful approach to study biomolecular systems. Herein, we apply a semi-empirical divide and conquer (DC) methodology coupled with a dielectric continuum model for the solvent, to explore the contribution of electrostatics, polarization and charge transfer to the interaction energy between barnase and barstar in their complex form. Molecular dynamic (MD) simulation was performed to account for the dynamic behavior of the complex. The results show that electrostatics, charge transfer and polarization favor the formation of the complex. Our study shows that electrostatics dominates the interaction between barnase and barstar ( approximately 73%), while charge transfer and polarization are approximately 21% and approximately 6%, respectively. Close inspection of the polarization and charge-transfer effects on the charge distribution of the complex reveals the existence of two, well localized, regions in barstar. The first region includes the residues between P27 and Y47 and the second region is between N65 and D83. Since no such regions could be detected in barnase clearly suggests that barstar is well optimized for efficiently binding barnase. Furthermore, using our interaction energy decomposition scheme, we were able to identify all residues that have been experimentally determined to be important for the complex formation and to suggest other residues never have been investigated. This suggests that our approach will be useful as an aid in further understanding protein-protein contacts for the ultimate goal to produce successful inhibitors for protein complexes.

The role of residue stability in transient protein-protein interactions involved in enzymatic phosphate hydrolysis. A computational study

Finding why protein-protein interactions (PPIs) are so specific can provide a valuable tool in a variety of fields. Statistical surveys of so-called transient complexes (like those relevant for signal transduction mechanisms) have shown a tendency of polar residues to participate in the interaction region. Following this scheme, residues in the unbound partners have to compete between interacting with water or interacting with other residues of the protein. On the other hand, several works have shown that the notion of active site electrostatic preorganization can be used to interpret the high efficiency in enzyme reactions. This preorganization can be related to the instability of the residues important for catalysis. In some enzymes, in addition, conformational changes upon binding to other proteins lead to an increase in the activity of the enzymatic partner. In this article the linear response approximation version of the semimacroscopic protein dipoles Langevin dipoles (PDLD/S-LRA) model is used to evaluate the stability of several residues in two phosphate hydrolysis enzymes upon complexation with their activating partners. In particular, the residues relevant for PPI and for phosphate hydrolysis in the CDK2/Cyclin A and Ras/GAP complexes are analyzed. We find that the evaluation of the stability of residues in these systems can be used to identify not only active site regions but it can also be used as a guide to locate "hot spots" for PPIs. We also show that conformational changes play a major role in positioning interfacing residues in a proper "energetic" orientation, ready to interact with the residues in the partner protein surface. Thus, we extend the preorganization theory to PPIs, extrapolating the results we obtained from the above-mentioned complexes to a more general case. We conclude that the correlation between stability of a residue in the surface and the likelihood that it participates in the interaction can be a general fact for transient PPIs.

Unraveling a hotspot for TCR recognition on HLA-A2: evidence against the existence of peptide-independent TCR binding determinants

T cell receptor (TCR) recognition of peptide takes place in the context of the major histocompatibility complex (MHC) molecule, which accounts for approximately two-thirds of the peptide/MHC buried surface. Using the class I MHC HLA-A2 and a large panel of mutants, we have previously shown that surface mutations that disrupt TCR recognition vary with the identity of the peptide. The single exception is Lys66 on the HLA-A2 alpha1 helix, which when mutated to alanine disrupts recognition for 93% of over 250 different T cell clones or lines, independent of which peptide is bound. Thus, Lys66 could serve as a peptide-independent TCR binding determinant. Here, we have examined the role of Lys66 in TCR recognition of HLA-A2 in detail. The structure of a peptide/HLA-A2 molecule with the K66A mutation indicates that although the mutation induces no major structural changes, it results in the exposure of a negatively charged glutamate (Glu63) underneath Lys66. Concurrent replacement of Glu63 with glutamine restores TCR binding and function for T cells specific for five different peptides presented by HLA-A2. Thus, the positive charge on Lys66 does not serve to guide all TCRs onto the HLA-A2 molecule in a manner required for productive signaling. Furthermore, electrostatic calculations indicate that Lys66 does not contribute to the stability of two TCR-peptide/HLA-A2 complexes. Our findings are consistent with the notion that each TCR arrives at a unique solution of how to bind a peptide/MHC, most strongly influenced by the chemical and structural features of the bound peptide. This would not rule out an intrinsic affinity of TCRs for MHC molecules achieved through multiple weak interactions, but for HLA-A2 the collective mutational data place limits on the role of any single MHC amino acid side-chain in driving TCR binding in a peptide-independent fashion.

Conformational energy landscape of the acyl pocket loop in acetylcholinesterase: a Monte Carlo-generalized Born model study

The X-ray crystal structure of the reaction product of acetylcholinesterase (AChE) with the inhibitor diisopropylphosphorofluoridate (DFP) showed significant structural displacement in a loop segment of residues 287-290. To understand this conformational selection, a Monte Carlo (MC) simulation study was performed of the energy landscape for the loop segment. A computational strategy was applied by using a combined simulated annealing and room temperature Metropolis sampling approach with solvent polarization modeled by a generalized Born (GB) approximation. Results from thermal annealing reveal a landscape topology of broader basin opening and greater distribution of energies for the displaced loop conformation, while the ensemble average of conformations at 298 K favored a shift in populations toward the native by a free-energy difference in good agreement with the estimated experimental value. Residue motions along a reaction profile of loop conformational reorganization are proposed where Arg-289 is critical in determining electrostatic effects of solvent interaction versus Coulombic charging.

Simulating redox coupled proton transfer in cytochromecoxidase: Looking for the proton bottleneck

Gaining a detailed understanding of the molecular nature of the redox coupled proton transfer in cytochrome c oxidase (COX) is one of the challenges of modern biophysics. The present work addresses this by integrating approaches for simulations of proton transport (PTR) and electron transfer (ET). The resulting method converts the electrostatic energies of different charge configurations and reorganization energies to free-energy profiles for different PTR and ET pathways. This approach provides for the first time a tool to study the actual activation barriers and kinetics of different feasible PTR processes in the cycle of COX. Using this tool, we explore the PTR through the bottleneck water molecules. It is found that a stepwise PTR along this commonly assumed path leads to far too high barriers and is, thus, inconsistent with the observed kinetics. Furthermore, the simulated free-energy profile does not provide a simple gating mechanism. Fortunately, we obtain reasonable kinetics when we consider a PTR that involves a concerted transfer of protons to and from E286. Finally, semi-qualitative considerations of the forward and backward barriers point toward open questions about the actual gating process and offer a feasible pumping mechanism. Although further studies are clearly needed, we believe that our approach offers a general and effective tool for correlating the structure of COX with its function.

The role of electrostatic interactions in calmodulin-peptide complex formation

The complex between calmodulin and the calmodulin-binding portion of smMLCKp has been studied. Electrostatic interactions have been anticipated to be important in this system where a strongly negative protein binds a peptide with high positive charge. Electrostatic interactions were probed by varying the pH in the range from 4 to 11 and by charge deletions in CaM and smMLCKp. The change in net charge of CaM from approximately -5 at pH 4.5 to -15 at pH 7.5 leaves the binding constant virtually unchanged. The affinity was also unaffected by mutations in CaM and charge substitutions in the peptide. The insensitivity of the binding constant to pH may seem surprising, but it is a consequence of the high charge on both protein and peptide. At low pH it is further attenuated by a charge regulation mechanism. That is, the protein releases a number of protons when binding the positively charged peptide. We speculate that the role of electrostatic interactions is to discriminate against unbound proteins rather than to increase the affinity for any particular target protein.