Application of the Free Energy Calculations to Study Drug-enzyme and Drug-dna Complexes

Applications of two free energy calculation approaches are presented to study drug-biomolecule complexes. The first method, the free energy perturbation (FEP) method and molecular dynamics simulations has been applied to study the JG-365 inhibitor bound to the HIV-aspartic protease. The FEP method has been applied to predict the consequence of replacing each of the seven peptide bonds in the JG-365 by trans-ethylene or fluoroethylene units. The necessary initial conformations of the inhibitor for "in water" perturbations have been found using neural network clustering approach applied to the long molecular dynamics trajectory of the inhibitor in water solution. The second method is applied to study binding free energies of some DNA-drug complexes and is based on analysis of long molecular dynamics trajectories by continuum solvent approach (MM/PBSA).     

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.     

Molecular dynamic simulation of wild type and mutants of the polymorphic amyloid NNQNTF segments of elk prion: structural stability and thermodynamic of association

A hexapeptide with amino acid sequence NNQNTF from the elk prion protein forms amyloid fibrils. Here we use molecular dynamic simulations of the oligomers and their single point glycine mutants to study their stability. In an effort to probe the structural stability and association thermodynamic in a realistic environment, all wildtype of NNQNTF polymorphic forms with different size and their corresponding double layer 5 strands single point glycine mutants were subjected to a total of 500 ns of explicit-solvent molecular dynamics (MD) simulation. Our results show that the structural stability of the NNQNTF oligomers increases with increasing the number of β-strands for double layers. Our results also demonstrated that hydrophobic interaction is the principle driving force to stabilize the adjacent β-strands while the steric zipper is responsible for holding the neighboring β-sheet layers together. We used MM-PBSA approach free energy calculations to determine the role of nonpolar effects, electrostatics and entropy in binding. Nonpolar effects remained consistently more favorable in wild type and mutants reinforcing the importance of hydrophobic effects in protein-protein binding. While entropy systematically opposed binding in all cases, there was no observed trend in the entropy difference between wildtype and glycine mutant. Free energy decomposition shows residues situated at the interface were found to make favorable contributions to the peptide-peptide association. The study of the wild type and mutants in an explicit solvent may provide valuable insight for amyloid aggregation inhibitor design efforts.     

Binding free energies and free energy components from molecular dynamics and Poisson-Boltzmann calculations. Application to amino acid recognition by aspartyl-tRNA synthetase

Specific amino acid binding by aminoacyl-tRNA synthetases (aaRS) is necessary for correct translation of the genetic code. Engineering a modified specificity into aminoacyl-tRNA synthetases has been proposed as a means to incorporate artificial amino acid residues into proteins in vivo. In a previous paper, the binding to aspartyl-tRNA synthetase of the substrate Asp and the analogue Asn were compared by molecular dynamics free energy simulations. Molecular dynamics combined with Poisson-Boltzmann free energy calculations represent a less expensive approach, suitable for examining multiple active site mutations in an engineering effort. Here, Poisson-Boltzmann free energy calculations for aspartyl-tRNA synthetase are first validated by their ability to reproduce selected molecular dynamics binding free energy differences, then used to examine the possibility of Asn binding to native and mutant aspartyl-tRNA synthetase. A component analysis of the Poisson-Boltzmann free energies is employed to identify specific interactions that determine the binding affinities. The combined use of molecular dynamics free energy simulations to study one binding process thoroughly, followed by molecular dynamics and Poisson-Boltzmann free energy calculations to study a series of related ligands or mutations is proposed as a paradigm for protein or ligand design.The binding of Asn in an alternate, "head-to-tail" orientation observed in the homologous asparagine synthetase is analyzed, and found to be more stable than the "Asp-like" orientation studied earlier. The new orientation is probably unsuitable for catalysis. A conserved active site lysine (Lys198 in Escherichia coli) that recognizes the Asp side-chain is changed to a leucine residue, found at the corresponding position in asparaginyl-tRNA synthetase. It is interesting that the binding of Asp is calculated to increase slightly (rather than to decrease), while that of Asn is calculated, as expected, to increase strongly, to the same level as Asp binding. Insight into the origin of these changes is provided by the component analyses. The double mutation (K198L,D233E) has a similar effect, while the triple mutation (K198L,Q199E,D233E) reduces Asp binding strongly. No binding measurements are available, but the three mutants are known to have no ability to adenylate Asn, despite the "Asp-like" binding affinities calculated here. In molecular dynamics simulations of all three mutants, the Asn ligand backbone shifts by 1-2 A compared to the experimental Asp:AspRS complex, and significant side-chain rearrangements occur around the pocket. These could reduce the ATP binding constant and/or the adenylation reaction rate, explaining the lack of catalytic activity in these complexes. Finally, Asn binding to AspRS with neutral K198 or charged H449 is considered, and shown to be less favorable than with the charged K198 and neutral H449 used in the analysis.     

Comparative analysis of binding sites of human meprins with hydroxamic acid derivative by molecular dynamics simulation study

Meprins are complex and highly glycosylated multi-domain enzymes that require post-translational modifications to reach full activity. Meprins are metalloproteases of the astacin family characterized by a conserved zinc-binding motif (HExxHxxGFxHExxRxDR). Human meprin-α and -β protease subunits are 55% identical at the amino acid level, however the substrate and peptide bond specificities vary markedly. Current work focuses on the critical amino acid residues in the non-primed subsites of human meprins-α and -β involved in inhibitor/ligand binding. To compare the molecular events underlying ligand affinity, homology modeling of the protease domain of humep-α and -β based on the astacin crystal structure followed by energy minimization and molecular dynamics simulation of fully solvated proteases with inhibitor Pro-Leu-Gly-hydroxamate in S subsites were performed. The solvent accessible surface area curve shows a decrease in solvent accessibility values at specific residues upon inhibitor binding. The potential energy, total energy, H-bond interactions, root mean square deviation and root mean square fluctuation plot reflect the subtle differences in the S subsite of the enzymes which interact with different residues at P site of the inhibitor.     

Insight into the Oseltamivir Resistance R292K Mutation in H5N1 Influenza Virus: A Molecular Docking and Molecular Dynamics Approach

H5N1 is a subtype of the influenza A virus that can cause disease in humans and many other animal species. Oseltamivir (Tamiflu) is a potent and selective antiviral drug employed to fight the flu virus in infected individuals by inhibiting neuraminidase (NA), a flu protein responsible for the release and spread of the progeny virions. However, oseltamivir resistance has become a critical problem. In particular, influenza strains with a R292K NA mutation are highly resistant to the oseltamivir. Though the biological functions of the mutations have previously been characterized, the structural basis behind the reduced catalytic activity and reduced protein level is not clear. In this study, molecular docking and molecular dynamics (MD) approach were employed to investigate the structural and dynamical effects throughout the protein structure and specifically, at the drug-binding pocket. Furthermore, potential of mean force was analyzed using explicit solvent MD simulations with the umbrella sampling method to explore the free energy of binding. It is believed that this study provides valuable guidance for the resistance management of oseltamivir and designing of more potent antiviral inhibitor.     

A computational analysis of the binding affinities of FKBP12 inhibitors using the MM-PB/SA method

The FK506-binding proteins have been targets of pharmaceutical interests over years. We have studied the binding of a set of 12 nonimmunosuppressive small-molecule inhibitors to FKBP12 through molecular dynamics simulations. Each complex was subjected to 1-ns MD simulation conducted in an explicit solvent environment under constant temperature and pressure. The binding free energy of each complex was then computed by the MM-PB/SA method in the AMBER program. Our MM-PB/SA computation produced a good correlation between the experimentally determined and the computed binding free energies with a correlation coefficient (R(2)) of 0.93 and a standard deviation as low as 0.30 kcal/mol. The vibrational entropy term given by the normal mode analysis was found to be helpful for achieving this correlation. Moreover, an adjustment to one weight factor in the PB/SA model was essential to correct the absolute values of the final binding free energies to a reasonable range. A head-to-head comparison of our MM-PB/SA model with a previously reported Linear Response Approximation (LRA) model suggested that the MM-PB/SA method is more robust in binding affinity prediction for this class of compounds.     

Conformational flexibility,binding energy,role of salt bridge and alanine-mutagenesis for c-Abl kinase complex

Abl kinase plays a decisive role in the mechanism of the most fatal human pathogen chronic mylogenous leukemia (CML). Here, we have carried out a comprehensive study about the conformational flexibility, role of salt bridge and the protein- ligand interaction for this kinase with its well-known inhibitor, Imatinib. We have performed molecular dynamics simulations for conformational behavior, investigated the salt bridges and calculated the binding free energy of Imatinib with MM-PB/SA method for Abl kinase complex. We also explored the role of salt-bridge in the kinase complex and its effect on binding activity of inhibitors. Furthermore, to investigate the importance of those residues which form salt bridges, we mutated them by Alanine with the help of Alanine scanning program. We noticed significant variations in total free energy of Imatinib in all possible mutations. The binding free energy of ligand for kinase receptor was analyzed by molecular mechanics Poission Boltzmann surface area (MM-PB/SA) method. These results suggest that conserved glutamic acid and lysine are necessary for stability of complex.     

Impact of calcium on N1 influenza neuraminidase dynamics and binding free energy

The highly pathogenic influenza strains H5N1 and H1N1 are currently treated with inhibitors of the viral surface protein neuraminidase (N1). Crystal structures of N1 indicate a conserved, high affinity calcium binding site located near the active site. The specific role of this calcium in the enzyme mechanism is unknown, though it has been shown to be important for enzymatic activity and thermostability. We report molecular dynamics (MD) simulations of calcium‐bound and calcium‐free N1 complexes with the inhibitor oseltamivir (marketed as the drug Tamiflu), independently using both the AMBER FF99SB and GROMOS96 force fields, to give structural insight into calcium stabilization of key framework residues. Y347, which demonstrates similar sampling patterns in the simulations of both force fields, is implicated as an important N1 residue that can “clamp” the ligand into a favorable binding pose. Free energy perturbation and thermodynamic integration calculations, using two different force fields, support the importance of Y347 and indicate a +3 to +5 kcal/mol change in the binding free energy of oseltamivir in the absence of calcium. With the important role of structure‐based drug design for neuraminidase inhibitors and the growing literature on emerging strains and subtypes, inclusion of this calcium for active site stability is particularly crucial for computational efforts such as homology modeling, virtual screening, and free energy methods. Proteins 2010. © 2010 Wiley‐Liss, Inc.     

Energetic contributions of amino acid residues and its cross-talk to delineate ligand-binding mechanism

Receptor-based QSAR approaches can enumerate the energetic contributions of amino acid residues toward ligand binding only when experimental binding affinity is associated. The structural data of protein-ligand complexes are witnessing a tremendous growth in the Protein Data Bank deposited with a few entries on binding affinity. We present here a new approach to compute the E nergetic CONT ributions of A mino acid residues and its possible C ross-T alk (ECONTACT) to study ligand binding using per-residue energy decomposition, molecular dynamics simulations and rescoring method without the need for experimental binding affinity. This approach recognizes potential cross-talks among amino acid residues imparting a nonadditive effect to the binding affinity with evidence of correlative motions in the dynamics simulations. The protein-ligand interaction energies deduced from multiple structures are decomposed into per-residue energy terms, which are employed as variables to principal component analysis and generated cross-terms. Out of 16 cross-talks derived from eight datasets of protein-ligand systems, the ECONTACT approach is able to associate 10 potential cross-talks with site-directed mutagenesis, free energy, and dynamics simulations data strongly. We modeled these key determinants of ligand binding using joint probability density function (jPDF) to identify cross-talks in protein structures. The top two cross-talks identified by ECONTACT approach corroborated with the experimental findings. Furthermore, virtual screening exercise using ECONTACT models better discriminated known inhibitors from decoy molecules. This approach proposes the jPDF metric to estimate the probability of observing cross-talks in any protein-ligand complex. The source code and related resources to perform ECONTACT modeling is available freely at https://www.gujaratuniversity.ac.in/econtact /.     

Molecular docking study and development of an empirical binding free energy model for phosphodiesterase 4 inhibitors

In the present work, several computational methodologies were combined to develop a model for the prediction of PDE4B inhibitors' activity. The adequacy of applying the ligand docking approach, keeping the enzyme rigid, to the study of a series of PDE4 inhibitors was confirmed by a previous molecular dynamics analysis of the complete enzyme. An exhaustive docking procedure was performed to identify the most probable binding modes of the ligands to the enzyme, including the active site metal ions and the surrounding structural water molecules. The enzyme-inhibitor interaction enthalpies, refined by using the semiempirical molecular orbital approach, were combined with calculated solvation free energies and entropy considerations in an empirical free energy model that enabled the calculation of binding free energies that correlated very well with experimentally derived binding free energies. Our results indicate that both the inclusion of the structural water molecules close to the ions in the binding site and the use of a free energy model with a quadratic dependency on the ligand free energy of solvation are important aspects to be considered for molecular docking investigations involving the PDE4 enzyme family.     

Molecular dynamics,free energy,and SPR analyses of the interactions between the SH2 domain of Grb2 and ErbB phosphotyrosyl peptides

We studied the interactions between the SH2 domain of growth factor receptor binding protein 2 (Grb2) and ErbB receptor-derived phosphotyrosyl peptides using molecular dynamics, free energy calculations, and surface plasmon resonance (SPR) analysis. Binding free energies for nine phosphotyrosyl peptides were calculated using the MM-PBSA continuum solvent method, and excellent qualitative agreement with the SPR experimental data, with a correlation coefficient of 0.92, was obtained. Consistent with previous experimental findings, phosphotyrosyl peptides with the consensus sequence pYXNX showed favorable binding affinity for the Grb2. Unexpectedly, phosphotyrosyl peptides with the consensus sequence pYQQD, which had not shown any specific binding affinity for the Grb2 in earlier studies, also showed favorable binding affinity for the Grb2 in our experimental and computational analyses. Component analysis of the calculated binding free energies revealed that van der Waals interaction between the Grb2 and the phosphotyrosyl peptide was the dominant factor for specificity and binding affinity. These results indicate that current methods of estimating binding free energies are efficient for obtaining important information about protein-protein interactions, which are essential for the transmission of signals in cellular signaling pathways.     

Probing the interaction mechanism of small molecule inhibitors with matriptase based on molecular dynamics simulation and free energy calculations

Matriptase is a serine protease associated with a wide variety of human tumors and carcinoma progression. Up to now, many promising anti-cancer drugs have been developed. However, the detailed structure–function relationship between inhibitors and matriptase remains elusive. In this work, molecular dynamics simulation and binding free energy studies were performed to investigate the biochemistry behaviors of two class inhibitors binding to matriptase. The binding free energies predicted by MM/GBSA methods are in good agreement with the experimental bioactivities, and the analysis of the individual energy terms suggests that the van der Waals interaction is the major driving force for ligand binding. The key residues responsible for achieving strong binding have been identified by the MM/GBSA free energy decomposition analysis. Especially, Trp215 and Phe99 had an important impact on active site architecture and ligand binding. The results clearly identify the two class inhibitors exist different binding modes. Through summarizing the two different modes, we have mastered some important and favorable interaction patterns between matriptase and inhibitors. Our findings would be helpful for understanding the interaction mechanism between the inhibitor and matriptase and afford important guidance for the rational design of potent matriptase inhibitors.     

Probing the binding mechanism of mercaptoguanine derivatives as inhibitors of HPPK by docking and molecular dynamics simulations

6-Hydroxymethyl-7,8-dihydropterin pyrophosphokinase (HPPK) is a promising antimicrobial target involved in the folate biosynthesis pathway. Although, the results from crystallographic studies of HPPK have attracted a great interest in the design of novel HPPK inhibitors, the mechanism of action of HPPK due to inhibitor binding remains questionable. Recently, mercaptoguanine derivatives were reported to inhibit the pyrophosphoryl transfer mechanism of Staphylococcus aureus HPPK (SaHPPK). The present study is an attempt to understand the SaHPPK-inhibitors binding mechanism and to highlight the key residues that possibly involve in the complex formation. To decipher these questions, we used the state-of-the-art advanced insilico approach such as molecular docking, molecular dynamics (MD), molecular mechanics-generalized Born surface area approach. Domain cross correlation and principle component analysis were applied to the snapshots obtained from MD revealed that the compounds with high binding affinity stabilize the conformational dynamics of SaHPPK. The binding free energy estimation showed that the van der Waals and electrostatic interactions played a vital role for the binding mechanism. Additionally, the predicted binding free energy was in good agreement with the experimental values (R² = .78). Moreover, the free energy decomposition on per-residue confirms the key residues that significantly contribute to the complex formation. These results are expected to be useful for rational design of novel SaHPPK inhibitors.     

A critical analysis of continuum electrostatics: the screened Coulomb potential--implicit solvent model and the study of the alanine dipeptide and discrimination of misfolded structures of proteins

An analysis of the screened Coulomb potential--implicit solvent model (SCP--ISM) is presented showing that general equations for both the electrostatic and solvation free energy can be derived in a continuum approach, using statistical averaging of the polarization field created by the solvent around the molecule. The derivation clearly shows how the concept of boundary, usually found in macroscopic approaches, is eliminated when the continuum model is obtained from a microscopic treatment using appropriate averaging techniques. The model is used to study the alanine dipeptide in aqueous solution, as well as the discrimination of native protein structures from misfolded conformations. For the alanine dipeptide the free energy surface in the phi--psi space is calculated and compared with recently reported results of a detailed molecular dynamics simulation using an explicit representation of the solvent, and with other available data. The study showed that the results obtained using the SCP--ISM are comparable to those of the explicit water calculation and compares favorably to the FDPB approach. Both transition states and energy minima show a high correlation (r > 0.98) with the results obtained in the explicit water analysis. The study of the misfolded structures of proteins comprised the analysis of three standard decoy sets, namely, the EMBL, Park and Levitt, and Baker's CASP3 sets. In all cases the SCP--ISM discriminated well the native structures of the proteins, and the best-predicted structures were always near-native (cRMSD approximately 2 A).     

蛋白质与极性配体复合物的绝对结合自由能计算

介绍了用分子动力学模拟与热力学积分法相结合,模拟蛋白质与配体的绝对结合自由能的方法.通过分子转换法,使蛋白质分子（包括水分子）与配体小分子之间的相互作用逐渐减弱 （或增强）至完全消失（或完全出现）. 运用体约束方法,计算了配体与受体结合后平动、转动自由度的丧失即熵效应所引起的自由能变化.以胰蛋白酶双突变体（D189G/G226D）与极性配体苯甲脒为例,研究了蛋白质活性部位与极性配体的相互作用对结合自由能的影响,该复合物绝对结合自由能的模拟结果(－15.5 kJ/mol)与实验值(－10.5 kJ/mol)相近.     

A computational analysis of pyrazole-based inhibitors binding to Hsp90 using molecular dynamics simulation and the MM-GBSA method

Owing to the key role of heat-shock protein 90 (Hsp90) in the evolution, development and disease pathogenesis of cancer, it has been an important target for anti-cancer chemotherapy over the years. A five-nanosecond molecular dynamics simulation combined with the calculation of the binding free energy was carried out to investigate the binding mechanisms of three Hsp90 inhibitors 4BH, 2E1 and 2D9 to Hsp90. The binding free energy of each complex was computed using the molecular mechanics–generalised Born surface area method. Detailed binding free energies between each inhibitor and residues of Hsp90 were calculated using a per-residue basis decomposition method. The detailed inhibitor–residue interaction provides insights into binding mechanisms and in-depth understanding of the structure–affinity relationship. This study suggests that van der Waals energy is primarily responsible for driving the binding of the inhibitors to Hsp90, and the three inhibitors bind to Hsp90 in a similar binding mode. However, a substituent in 2D9 leads to higher binding free energy than the other two inhibitors. These data may assist in designing new potent drugs to combat cancer.     

Propagation of dynamic changes in barnase upon binding of barstar: an NMR and computational study

NMR spectroscopy and computer simulations were used to examine changes in chemical shifts and in dynamics of the ribonuclease barnase that result upon binding to its natural inhibitor barstar. Although the spatial structures of free and bound barnase are very similar, binding results in changes of the dynamics of both fast side-chains, as revealed by (2)H relaxation measurements, and NMR chemical shifts in an extended beta-sheet that is located far from the binding interface. Both side-chain dynamics and chemical shifts are sensitive to variations in the ensemble populations of the inter-converting molecular states, which can escape direct structural observation. Molecular dynamics simulations of free barnase and barnase in complex with barstar, as well as a normal mode analysis of barnase using a Gaussian network model, reveal relatively rigid domains that are separated by the extended beta-sheet mentioned above. The observed changes in NMR parameters upon ligation can thus be rationalized in terms of changes in inter-domain dynamics and in populations of exchanging states, without measurable structural changes. This provides an alternative model for the propagation of a molecular response to ligand binding across a protein that is based exclusively on changes in dynamics.     

Effect of hydrophobic and hydrogen bonding interactions on the potency of ß-alanine analogs of G-protein coupled glucagon receptor inhibitors

G-protein coupled glucagon receptors (GCGRs) play an important role in glucose homeostasis and pathophysiology of Type-II Diabetes Mellitus (T2DM). The allosteric pocket located at the trans-membrane domain of GCGR consists of hydrophobic (TM5) and hydrophilic (TM7) units. Hydrophobic interactions with the amino acid residues present at TM5, found to facilitate the favorable orientation of antagonist at GCGR allosteric pocket. A statistically robust and highly predictive 3D-QSAR model was developed using 58 β-alanine based GCGR antagonists with significant variation in structure and potency profile. The correlation coefficient (R²) and cross-validation coefficient (Q²) of the developed model were found to be 0.9981 and 0.8253, respectively at the PLS factor of 8. The analysis of the favorable and unfavorable contribution of different structural features on the glucagon receptor antagonists was done by 3D-QSAR contour plots. Hydrophobic and hydrogen bonding interactions are found to be main dominating non-bonding interactions in docking studies. Presence of highest occupied molecular orbital (HOMO) in the polar part and lowest unoccupied molecular orbital (LUMO) in the hydrophobic part of antagonists leads to favorable protein-ligand interactions. Molecular mechanics/generalized born surface area (MM/GBSA) calculations showed that van der Waals and nonpolar solvation energy terms are crucial components for thermodynamically stable binding of the inhibitors. The binding free energy of highly potent compound was found to be −63.475 kcal/mol; whereas the least active compound exhibited binding energy of −41.097 kcal/mol. Further, five 100 ns molecular dynamics simulation (MD) simulations were done to confirm the stability of the inhibitor-receptor complex. Outcomes of the present study can serve as the basis for designing improved GCGR antagonists.