Insights into protein-protein binding by binding free energy calculation and free energy decomposition for the Ras-Raf and Ras-RalGDS complexes

Absolute binding free energy calculations and free energy decompositions are presented for the protein-protein complexes H-Ras/C-Raf1 and H-Ras/RalGDS. Ras is a central switch in the regulation of cell proliferation and differentiation. In our study, we investigate the capability of the molecular mechanics (MM)-generalized Born surface area (GBSA) approach to estimate absolute binding free energies for the protein-protein complexes. Averaging gas-phase energies, solvation free energies, and entropic contributions over snapshots extracted from trajectories of the unbound proteins and the complexes, calculated binding free energies (Ras-Raf: -15.0(+/-6.3)kcal mol(-1); Ras-RalGDS: -19.5(+/-5.9)kcal mol(-1)) are in fair agreement with experimentally determined values (-9.6 kcal mol(-1); -8.4 kcal mol(-1)), if appropriate ionic strength is taken into account. Structural determinants of the binding affinity of Ras-Raf and Ras-RalGDS are identified by means of free energy decomposition. For the first time, computationally inexpensive generalized Born (GB) calculations are applied in this context to partition solvation free energies along with gas-phase energies between residues of both binding partners. For selected residues, in addition, entropic contributions are estimated by classical statistical mechanics. Comparison of the decomposition results with experimentally determined binding free energy differences for alanine mutants of interface residues yielded correlations with r(2)=0.55 and 0.46 for Ras-Raf and Ras-RalGDS, respectively. Extension of the decomposition reveals residues as far apart as 25A from the binding epitope that can contribute significantly to binding free energy. These "hotspots" are found to show large atomic fluctuations in the unbound proteins, indicating that they reside in structurally less stable regions. Furthermore, hotspot residues experience a significantly larger-than-average decrease in local fluctuations upon complex formation. Finally, by calculating a pair-wise decomposition of interactions, interaction pathways originating in the binding epitope of Raf are found that protrude through the protein structure towards the loop L1. This explains the finding of a conformational change in this region upon complex formation with Ras, and it may trigger a larger structural change in Raf, which is considered to be necessary for activation of the effector by Ras.     

Ferric uptake regulator protein: binding free energy calculations and per-residue free energy decomposition

Iron homeostasis is, in many bacterial species, mediated by the ferric uptake regulator (Fur). A regulatory site able to bind iron to activate Fur for DNA binding has been described, and a structural zinc site essential for the dimerization has also been proposed. They have been localized and named site 1 and site 2, respectively, from the crystal structure of a zinc-substituted Pseudomonas aeruginosa Fur (PA-Fur). Notwithstanding the studies on Fur proteins from various species, both the precise site of iron binding and the effect on DNA binding affinity are still controversial. These issues were investigated here by molecular dynamics simulations and free energy calculations. Simulations were performed for eight molecular systems represented by the three forms of Fur, that is, apo Fur, metal-substituted Fur, and Fur complexed with DNA. Because of the lack of a Fur-DNA complex crystal structure, the recently published model based on mass spectrometry experiments on Escherichia coli Fur (EC-Fur), and the crystal structure of PA-Fur, was used, after adjustment to adopt a symmetric conformation. The simulation results suggest that the formerly proposed site 2 is, in fact, the regulatory iron-sensing site. The calculations also predict that Fe(2+) at site 2 is hexacoordinated having an octahedral environment with only nitrogen and oxygen atoms, which is in accordance with previous spectroscopic characterizations. Energy decomposition pinpoints H87 as an additional amino acid that defines the regulatory metal site. Finally, free energy decomposition analysis reveals a number of amino acids potentially important in dimerization and in DNA binding.     

Study of the insulin dimerization: binding free energy calculations and per-residue free energy decomposition

A calculation of the binding free energy for the dimerization of insulin has been performed using the molecular mechanics-generalized Born surface area approach. The calculated absolute binding free energy is -11.9 kcal/mol, in approximate agreement with the experimental value of -7.2 kcal/mol. The results show that the dimerization is mainly due to nonpolar interactions. The role of the hydrogen bonds between the 2 monomers appears to give the direction of the interactions. A per-atom decomposition of the binding free energy has been performed to identify the residues contributing most to the self association free energy. Residues B24-B26 are found to make the largest favorable contributions to the dimerization. Other residues situated at the interface between the 2 monomers were found to make favorable but smaller contributions to the dimerization: Tyr B16, Val B12, and Pro B28, and to an even lesser extent, Gly B23. The energy decomposition on a per-residue basis is in agreement with experimental alanine scanning data. The results obtained from a single trajectory (i.e., the dimer trajectory is also used for the monomer analysis) and 2 trajectories (i.e., separate trajectories are used for the monomer and dimer) are similar.     

Molecular mechanisms of resistance: free energy calculations of mutation effects on inhibitor binding to HIV-1 protease

The changes in the inhibitor binding constants due to the mutation of isoleucine to valine at position 84 of HIV-1 protease are calculated using molecular dynamics simulations. The calculations are done for three potent inhibitors--KNI-272, L-735,524 (indinavir or MK-639), and Ro 31-8959 (saquinavir). The calculations agree with the experimental data both in terms of an overall trend and in the magnitude of the resulting free energy change. HIV-1 protease is a homodimer, so each mutation causes two changes in the enzyme. The decrease in the binding free energy from each mutated side chain differs among the three inhibitors and correlates well with the size of the cavities induced in the protein interior near the mutated residue. The cavities are created as a result of a mutation to a smaller side chain, but the cavities are less than would be predicted from the wild-type structures, indicating that there is significant relaxation to partially fill the cavities.     

Investigation of inhibitory properties of some hydrazone compounds on hCA I,hCA II and AChE enzymes

Recently, inhibition of carbonic anhydrase (hCA) and acetylcholinesterase (AChE) have appeared as a promising approach for pharmacological intervention in a variety of disorders such as glaucoma, epilepsy, obesity, cancer, and Alzheimer’s disease. Keeping this in mind, N,N′-bis[(1-aryl-3-heteroaryl)propylidene]hydrazine dihydrochlorides, N1-N11, P1, P4-P8, and R1-R6, were synthesized to investigate their inhibitory activity against hCA I, hCA II, and AChE enzymes. All compounds in N, P, and R-series inhibited hCAs (I and II) and AChE more efficiently than the reference compounds acetazolamide (AZA), and tacrine. According to the activity results, the most effective inhibitory compounds were in R-series with the K_i values of 203 ± 55–473 ± 67 nM and 200 ± 34–419 ± 94 nM on hCA I, and hCA II, respectively. N,N′-Bis[1-(4-fluorophenyl)-3-(morpholine-4-yl)propylidene]hydrazine dihydrochlorides, N8, in N-series, N,N′-Bis[1-(4-hydroxyphenyl)-3-(piperidine-1-yl)propylidene]hydrazine dihydrochlorides, P4, in P-series, and N,N′-bis[1-(4-chlorophenyl)-3-(pyrrolidine-1-yl)propylidene]hydrazine dihydrochlorides, R5, in R-series were the most powerful compounds against hCA I with the K_i values of 438 ± 65 nM, 344 ± 64 nM, and 203 ± 55 nM, respectively. Similarly, N8, P4, and R5 efficiently inhibited hCA II isoenzyme with the K_i values of 405 ± 60 nM, 327 ± 80 nM, and 200 ± 34 nM, respectively. On the other hand, P-series compounds had notable inhibitory effect against AChE than the reference compound tacrine and the K_i values were between 66 ± 20 nM and 128 ± 36 nM. N,N′-Bis[1-(4-fluorophenyl)-3-(piperidine-1-yl)propylidene]hydrazine dihydrochlorides, P7, was the most potent compound on AChE with the K_i value of 66 ± 20 nM. The other most promising compounds, N,N′-bis[1-(4-hydroxyphenyl)-3-(morpholine-4-yl)propylidene]hydrazine dihydrochlorides, N4 in N-series and N,N′-bis[1-(4-hydroxyphenyl)-3-(pyrrolidine-1-yl)propylidene]hydrazine dihydrochlorides, R4 in R-series were againts AChE with the K_i values of 119 ± 20 nM, 88 ± 14 nM, respectively.     

Substrate binding mode and reaction mechanism of undecaprenyl pyrophosphate synthase deduced from crystallographic studies

Undecaprenyl pyrophosphate synthase (UPPs) catalyzes eight consecutive condensation reactions of farnesyl pyrophosphate (FPP) with isopentenyl pyrophosphate (IPP) to form a 55-carbon long-chain product. We previously reported the crystal structure of the apo-enzyme from Escherichia coli and the structure of UPPs in complex with sulfate ions (resembling pyrophosphate of substrate), Mg(2+), and two Triton molecules (product-like). In the present study, FPP substrate was soaked into the UPPs crystals, and the complex structure was solved. Based on the crystal structure, the pyrophosphate head group of FPP is bound to the backbone NHs of Gly29 and Arg30 as well as the side chains of Asn28, Arg30, and Arg39 through hydrogen bonds. His43 is close to the C2 carbon of FPP and may stabilize the farnesyl cation intermediate during catalysis. The hydrocarbon moiety of FPP is bound with hydrophobic amino acids including Leu85, Leu88, and Phe89, located on the alpha3 helix. The binding mode of FPP in cis-type UPPs is apparently different from that of trans-type and many other prenyltransferases which utilize Asprich motifs for substrate binding via Mg(2+). The new structure provides a plausible mechanism for the catalysis of UPPs.     

RosettaDDGPrediction for high-throughput mutational scans: From stability to binding

Reliable prediction of free energy changes upon amino acid substitutions (ΔΔGs) is crucial to investigate their impact on protein stability and protein–protein interaction. Advances in experimental mutational scans allow high-throughput studies thanks to multiplex techniques. On the other hand, genomics initiatives provide a large amount of data on disease-related variants that can benefit from analyses with structure-based methods. Therefore, the computational field should keep the same pace and provide new tools for fast and accurate high-throughput ΔΔG calculations. In this context, the Rosetta modeling suite implements effective approaches to predict folding/unfolding ΔΔGs in a protein monomer upon amino acid substitutions and calculate the changes in binding free energy in protein complexes. However, their application can be challenging to users without extensive experience with Rosetta. Furthermore, Rosetta protocols for ΔΔG prediction are designed considering one variant at a time, making the setup of high-throughput screenings cumbersome. For these reasons, we devised RosettaDDGPrediction, a customizable Python wrapper designed to run free energy calculations on a set of amino acid substitutions using Rosetta protocols with little intervention from the user. Moreover, RosettaDDGPrediction assists with checking completed runs and aggregates raw data for multiple variants, as well as generates publication-ready graphics. We showed the potential of the tool in four case studies, including variants of uncertain significance in childhood cancer, proteins with known experimental unfolding ΔΔGs values, interactions between target proteins and disordered motifs, and phosphomimetics. RosettaDDGPrediction is available, free of charge and under GNU General Public License v3.0, at https://github.com/ELELAB/RosettaDDGPrediction .     

Kinetic and crystallographic studies of the role of tyrosine 7 in the active site of human carbonic anhydrase II

The rate limiting step in catalysis of bicarbonate dehydration by human carbonic anhydrase II (HCA II) is an intramolecular proton transfer from His64 to the zinc-bound hydroxide. We have examined the role of Tyr7 using site-specific mutagenesis and measuring catalysis by the ¹⁸O exchange method using membrane inlet mass spectrometry. The side chain of Tyr7 in HCA II extends into the active-site cavity about 7 Å from the catalytic zinc atom. Replacement of Tyr7 with eight other amino acids had no effect on the interconversion of bicarbonate and CO₂, but in some cases caused enhancements in the rate constant of proton transfer by nearly 10-fold. The variant Y7I HCA II enhanced intramolecular proton transfer approximately twofold; its structure was determined by X-ray crystallography at 1.5 Å resolution. No changes were observed in the ordered solvent structure in the active-site cavity or in the conformation of the side chain of the proton shuttle His64. However, the first 11 residues of the amino-terminal chain in Y7I HCA II assumed an alternate conformation compared with the wild type. Differential scanning calorimetry showed variants at position 7 had a melting temperature approximately 8 °C lower than that of the wild type.     

Computation of entropy contribution to protein-ligand binding free energy

The entropy contribution ΔS to protein-ligand binding free energy is studied for nine protein-lipid complexes. The entropy effect from the loss of the translational/rotational degrees of freedom (ΔS _tr) is calculated using the ideal gas approach. The change in the vibrational entropy (ΔS _vib) is calculated using the effective quantum oscillator approach with frequencies derived from the coordinate covariance matrix, so the inharmonic effects are taken into account. The change in the entropy of solvation (ΔS _solv) is considered using the binomial cell model (developed by the authors) for the hydrophobic effect. The entropy contribution from loss of conformations that are available for the free ligand (ΔS _conf) is also estimated. It is revealed that the negative in view of binding term ΔS _tr is only partly compensated by increasing of ΔS _vib, so T(ΔS _tr + ΔS _vib + ΔS _conf) < 0 for all complexes under investigation, but taking into account ΔS _solv leads to significantly increased ΔS. For all complexes except biotin-streptavidin, the results are found to be in reasonable agreement with experimental data. Published in Russian in Biokhimiya, 2007, Vol. 72, No. 7, pp. 963–973.     

The systematic modeling studies and free energy calculations of the phenazine compounds as anti-tuberculosis agents

Abstract

Phenazine compounds have good activity against Mycobacterium tuberculosis (MTB). Based on the reported activities that were obtained in MTB H37Rv, a three-dimensional quantitative structure–activity relationship (3D-QSAR) model was built to design novel compounds against MTB. A fivefold cross-validation method and external validation were used to analyze the accuracy of forecasting. The model has a cross-validation coefficient q²=0.7 and a non-cross-validation coefficient r² = 0.903, indicating that the model has good predictive possibility. The design of anti-pneumococcus MTB compounds was guided by the obtained 3D-QSAR model, and several compounds with better activity were obtained. To test the activity of these compounds, molecular docking, molecular dynamics simulation, and post-simulation analysis of the already reported drug targets in MTB were carried out. Among the total 15 drug targets, only three targets (Rv2361c, Rv2965c, and Rv3048c) were selected based on the docking results. Initial results reported that these compounds possessed good inhibition activity for Rv2361c. The top nine complexes of Rv2361 ligands were only subjected to MD simulation which resulted in a stable dynamics of the structures and showed a residual fluctuation in inhibitors binding pocket. Free energy reported that overall, the derivatives hold strong energy against the protein target. Energetic contribution results showed that residues, Asp76, Arg80, Asn124, Arg127, Arg244, and Arg250, play a major role in total energy. Systems biology approach validates shortlisted drug effect on the entire system which might be useful to predict potential drug in wet lab as well.

Communicated by Ramaswamy H. Sarma     

Deciphering the mechanism behind the varied binding activities of COXIBs through Molecular Dynamic Simulations,MM-PBSA binding energy calculations and per-residue energy decomposition studies

COX-2 is a well-known drug target in inflammatory disorders. COX-1/COX-2 selectivity of NSAIDs is crucial in assessing the gastrointestinal side effects associated with COX-1 inhibition. Celecoxib, rofecoxib, and valdecoxib are well-known specific COX-2 inhibiting drugs. Recently, polmacoxib, a COX-2/CA-II dual inhibitor has been approved by the Korean FDA. These COXIBs have similar structure with diverse activity range. Present study focuses on unraveling the mechanism behind the 10-fold difference in the activities of these sulfonamide-containing COXIBs. In order to obtain insights into their binding with COX-2 at molecular level, molecular dynamics simulations studies, and MM-PBSA approaches were employed. Further, per-residue decomposition of these energies led to the identification of crucial amino acids and interactions contributing to the differential binding of COXIBs. The results clearly indicated that Leu338, Ser339, Arg499, Ile503, Phe504, Val509, and Ser516 (Leu352, Ser353, Arg513, Ile517, Phe518, Val523, and Ser530 in PGHS-1 numbering) were imperative in determining the activity of these COXIBs. The binding energies and energy contribution of various residues were similar in all the three simulations. The results suggest that hydrogen bond interaction between the hydroxyl group of Ser516 and five-membered ring of diarylheterocycles augments the affinity in COXIBs. The SAR of the inhibitors studied and the per-residue energy decomposition values suggested the importance of Ser516. Additionally, the positive binding energy obtained with Arg106 explains the binding of COXIBs in hydrophobic channel deep in the COX-2 active site. The findings of the present work would aid in the development of potent COX-2 inhibitors.     

Reaction path and free energy calculations of the transition between alternate conformations of HIV-1 protease

Two different structures of ligand-free HIV protease have been determined by X-ray crystallography. These structures differ in the position of two 12 residue, β-hairpin regions (or “flaps”) which cap the active site. The movements of the flaps must be involved in the binding of substrates since, in either conformation, the flaps block the binding site. One of these structures is similar to structures of the ligand-bound enzyme; however, the importance of both structures to enzyme function is unclear. This transformation takes place on a time scale too long for conventional molecular dynamics simulations, so the process was studied by first identifying a reaction path between the two structures and then calculating the free energy along this path using umbrella sampling. For the ligand-free enzyme, it is found that the two structures are nearly equally stable, with the ligand-bound-type structure being less stable, consistent with X-ray crystallography data. The more stable open structure does not have a lower potential energy, but is stabilized by entropy. The transition occurs through a collapse and reformation of the β-sheet structure of the conformationally flexible, glycine-rich flap ends. Additionally, some problems in studying conformational changes in proteins through the use of a single reaction path are addressed. Proteins 32:7–16, 1998. © 1998 Wiley-Liss, Inc.     

Insight into the key active sites of afChiA1 based on molecular dynamics simulations and free energy calculations

In this study, the binding of the enzyme chitinase A1 (afChiA1) from the plant-type Aspergillus fumigatus with four potent inhibitors, allosamidin (ASM), acetazolamide (AZM), 8-chloro-theophylline (CTP) and kinetin (KIT) is investigated by molecular docking, molecular dynamics simulation and binding free energy calculation. The results reveal that the electrostatic interactions play an important role in the stabilisation of the binding of afChiA1 with inhibitors. Based on the binding energy of afChiA1-ligands, the key residues (Gln37 and Trp312) in the active binding pocket of the complex systems are confirmed by molecular mechanics/Poisson–Boltzmann surface area method, and the active inhibitors, ASM and AZM, both could form strong interaction with Gln37 and Trp312, and the non-active ligands, CTP and KIT, could not interact with these two residues, which is consistent with the result of experimental report. Then, it is identified that Gln37 and Trp312 should be one of the important active site residues of afChiA1.     

RosettaDDGPrediction for high‐throughput mutational scans: From stability to binding

Reliable prediction of free energy changes upon amino acid substitutions (ΔΔGs) is crucial to investigate their impact on protein stability and protein–protein interaction. Advances in experimental mutational scans allow high‐throughput studies thanks to multiplex techniques. On the other hand, genomics initiatives provide a large amount of data on disease‐related variants that can benefit from analyses with structure‐based methods. Therefore, the computational field should keep the same pace and provide new tools for fast and accurate high‐throughput ΔΔG calculations. In this context, the Rosetta modeling suite implements effective approaches to predict folding/unfolding ΔΔGs in a protein monomer upon amino acid substitutions and calculate the changes in binding free energy in protein complexes. However, their application can be challenging to users without extensive experience with Rosetta. Furthermore, Rosetta protocols for ΔΔG prediction are designed considering one variant at a time, making the setup of high‐throughput screenings cumbersome. For these reasons, we devised RosettaDDGPrediction, a customizable Python wrapper designed to run free energy calculations on a set of amino acid substitutions using Rosetta protocols with little intervention from the user. Moreover, RosettaDDGPrediction assists with checking completed runs and aggregates raw data for multiple variants, as well as generates publication‐ready graphics. We showed the potential of the tool in four case studies, including variants of uncertain significance in childhood cancer, proteins with known experimental unfolding ΔΔGs values, interactions between target proteins and disordered motifs, and phosphomimetics. RosettaDDGPrediction is available, free of charge and under GNU General Public License v3.0, at https://github.com/ELELAB/RosettaDDGPrediction.     

Empirical free energy calculation: comparison to calorimetric data.

An effective free energy potential, developed originally for binding free energy calculation, is compared to calorimetric data on protein unfolding, described by a linear combination of changes in polar and nonpolar surface areas. The potential consists of a molecular mechanics energy term calculated for a reference medium (vapor or nonpolar liquid), and empirical terms representing solvation and entropic effects. It is shown that, under suitable conditions, the free energy function agrees well with the calorimetric expression. An additional result of the comparison is an independent estimate of the side-chain entropy loss, which is shown to agree with a structure-based entropy scale. These findings confirm that simple functions can be used to estimate the free energy change in complex systems, and that a binding free energy evaluation model can describe the thermodynamics of protein unfolding correctly. Furthermore, it is shown that folding and binding leave the sum of solute-solute and solute-solvent van der Waals interactions nearly invariant and, due to this invariance, it may be advantageous to use a nonpolar liquid rather than vacuum as the reference medium.     

Insights into the functional role of protonation states in the HIV-1 protease-BEA369 complex: molecular dynamics simulations and free energy calculations

The molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) method combined with molecular dynamics (MD) simulations were used to investigate the functional role of protonation in human immunodeficiency virus type 1 (HIV-1) protease complexed with the inhibitor BEA369. Our results demonstrate that protonation of two aspartic acids (Asp25/Asp25′) has a strong influence on the dynamics behavior of the complex, the binding free energy of BEA369, and inhibitor–residue interactions. Relative binding free energies calculated using the MM-PBSA method show that protonation of Asp25 results in the strongest binding of BEA369 to HIV-1 protease. Inhibitor–residue interactions computed by the theory of free energy decomposition also indicate that protonation of Asp25 has the most favorable effect on binding of BEA369. In addition, hydrogen-bond analysis based on the trajectories of the MD simulations shows that protonation of Asp25 strongly influences the water-mediated link of a conserved water molecule, Wat301. We expect that the results of this study will contribute significantly to binding calculations for BEA369, and to the design of high affinity inhibitors.     

Investigation on the binding mechanism of loratinib with the c-ros oncogene 1 (ROS1) receptor tyrosine kinase via molecular dynamics simulation and binding free energy calculations

The c-ros oncogene 1 (ROS1) has proven to be an important cancer target for the treatment of various human cancers. The anaplastic lymphoma kinase inhibitor crizotinib has been granted approval for the treatment of patients with ROS1 positive metastatic non-small-cell lung cancer by the Food and Drug Administration on 2016. However, serious resistance due to the secondary mutation of glycine 2032 to arginine (G2032R) was developed in clinical studies. Loratinib (PF-06463922), a macrocyclic analog of crizotinib, showed significantly improved inhibitory activity against wild–type (WT) ROS1 and ROS1^G2032R mutant. To provide insights into the inhibition mechanism, molecular dynamics simulations and free energy calculations were carried out for the complexes of loratinib with WT and G2032R mutated ROS1. The apo-ROS1^WT and apo-ROS1^G2032R systems showed similar RMSF distributions, while ROS1^G2032R-loratinib showed significantly higher than that of WT ROS1-loratinib, which revealed that the binding of loratinib to ROS1^G2032R significantly interfered the ?uctuation of protein. Calculations of binding free energies indicate that G2032R mutation significantly reduces the binding affinity of loratinib for ROS1, which arose mostly from the increase of conformation entropy and the decrease of solvation energy. Furthermore, detailed per-residue binding free energies highlighted the increased and decreased contributions of some residues in the G2032R mutated systems. The present study revealed the detailed inhibitory mechanism of loratinib as potent WT and G2032R mutated ROS1 inhibitor, which was expected to provide a basis for rational drug design.