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).     

Drug-DNA recognition: energetics and implications for design

In this article we review thermodynamic studies designed to examine the interaction of low molecular weight ligands or drugs with DNA. Over the past 10 years there has been an increase in the number of rigorous biophysical studies of DNA-drug interactions and considerable insight has been gained into the energetics of these binding reactions. The advent of high-sensitivity calorimetric techniques has meant that the energetics of DNA-drug association reactions can be probed directly and enthalpic and entropic contributions to the binding free energy established. There are two principal consequences arising from this type of work, firstly three-dimensional structures of DNA-drug complexes from X-ray and NMR studies can be put into a thermodynamic context and the energetics responsible for stabilizing the observed structures can be more fully understood. Secondly, any rational approach to structure-based drug design requires a fundamental base of knowledge where structural detail and thermodynamic data on complex formation are intimately linked. Therefore these types of studies allow a set of general guidelines to be established, which can then be used to develop drug design algorithms. In this review we describe recent breakthroughs in duplex DNA-directed drug design and also discuss how similar principles are now being used to target higher-order DNA molecules, for example, triplex (three-stranded) and tetraplex (four-stranded) structures.     

Molecular dynamics simulations and MM–PBSA calculations of the lectin from snowdrop (Galanthus nivalis)

Galanthus nivalis agglutinin (GNA), a mannose-specific lectin from snowdrop bulbs, is a member of the monocot mannose-specific lectin family and exhibits antiviral activity toward HIV. In the present study, molecular dynamics (MD) simulations were performed to study the interaction between GNA and its carbohydrate ligand over a specific time span. By analysis of the secondary structures, it was observed that the GNA conformation maintains rather stable along the trajectories and the high fluctuations were only centered on the carbohydrate recognition domains. Our MD simulations also reproduced most of the hydrogen bonds observed in the x-ray crystal structure. Furthermore, the obtained MD trajectories were used to estimate the binding free energy of the complex using the molecular mechanics/Poisson Boltzmann surface area (MM-PBSA) method. It was revealed by the inspection of the binding free energy components that the major contributions to the complex stability arose from electrostatic interactions.     

Probing the dynamics of restriction endonuclease NgoMIV‐DNA interaction by single‐molecule FRET

Many type II restriction endonucleases require two copies of their recognition sequence for optimal activity. Concomitant binding of two DNA sites by such an enzyme produces a DNA loop. Here we exploit single‐molecule Förster resonance energy transfer (smFRET) of surface‐immobilized DNA fragments to study the dynamics of DNA looping induced by tetrameric endonuclease NgoMIV. We have employed a DNA fragment with two NgoMIV recognition sites and a FRET dye pair such that upon protein‐induced DNA looping the dyes are brought to close proximity resulting in a FRET signal. The dynamics of DNA ‐ NgoMIV interactions proved to be heterogeneous, with individual smFRET trajectories exhibiting broadly different average looped state durations. Distinct types of the dynamics were attributed to different types of DNA ‐ protein complexes, mediated either by one NgoMIV tetramer simultaneously bound to two specific sites (“slow” trajectories) or by semi‐specific interactions of two DNA‐bound NgoMIV tetramers (“fast” trajectories), as well as to conformational heterogeneity of individual NgoMIV molecules.     

Ranking ligand binding affinities with avidin: a molecular dynamics-based interaction energy study.

The binding of 14 biotin analogues to avidin is examined to evaluate the viability of calculating binding free energy based on molecular dynamics (MD) trajectories. Two approaches were investigated in this work. The first one uses the linear interaction energy approximation, while the other approach utilizes the interaction free energy. The results obtained from these two methods were found to correlate well with the experimental binding free energy data for 10 out of 14 ligands. For the other four ligands, both methods overestimate their binding strength by more than 7 kcal/mol. Free energy calculations using the thermodynamic integration method are employed to understand this overestimation. The effect of protein flexibility on binding free energy calculation and the effect of charged or neutral ligands on the calculated results are discussed. MD simulations are shown to be able to provide insight into the interactions occurring in the active site and the origins of variations in binding free energy.     

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.     

Do water molecules mediate protein-DNA recognition?

A comprehensive analysis of interfacial water molecules in the structures of 109 unique protein-DNA complexes is presented together with a new view on their role in protein-DNA recognition. Location of interfacial water molecules as reported in the crystal structures and as emerging from a series of molecular dynamics studies on protein-DNA complexes with explicit solvent and counterions, was analyzed based on their acceptor, donor hydrogen bond relationships with the atoms and residues of the macromolecules, electrostatic field calculations and packing density considerations. Water molecules for the purpose of this study have been categorized into four classes: viz. (I) those that contact both the protein and the DNA simultaneously and thus mediate recognition directly; (II) those that contact either the protein or the DNA exclusively via hydrogen bonds solvating each solute separately; (III) those that contact the hydrophobic groups in either the protein or the DNA; and, lastly (IV) those that contact another water molecule. Of the 17,963 crystallographic water molecules under examination, about 6% belong to class I and 76% belong to class II. About three-fourths of class I and class II water molecules are exclusively associated with hydrogen bond acceptor atoms of both protein and DNA. Noting that DNA is polyanionic, it is significant that a majority of the crystallographically observed water molecules as well as those from molecular dynamics simulations should be involved in facilitating binding by screening unfavorable electrostatics. Less than 2% of the reported water molecules occur between hydrogen bond donor atoms of protein and acceptor atoms of DNA. These represent cases where protein atoms cannot reach out to DNA to make favorable hydrogen bond interactions due to packing/structural restrictions and interfacial water molecules provide an extension to side-chains to accomplish hydrogen bonding.     

Mechanism of the difference in the binding affinity of E. coli tRNAGln to glutaminyl-tRNA synthetase caused by noninterface nucleotides in variable loop

Aminoacyl-tRNA synthetases (ARSs) distinguish their cognate tRNAs from many other kinds of tRNAs, despite the very similar tertiary structures of tRNAs. Many researchers have supported the view that this recognition is achieved by intermolecular interactions between tRNA and ARS. However, one of the aptamers of Escherichia coli glutamine specific tRNA, var-AGGU, has a higher affinity to ARS than the wild-type, although the sequence difference only lies in the variable loop located on the opposite side of the binding interface with ARS. To understand the reason for the difference in affinity, we did molecular dynamics simulations on tRNAs and their complexes with ARS. We calculated the enthalpic and entropic contributions to the binding free energy with the molecular mechanics-Poisson-Boltzmann/surface area method and found that the entropic difference plays an important role in the difference in binding free energies. During the molecular dynamics simulations, dynamic rearrangements of hydrogen bonds occurred in the tertiary core region of the wild-type tRNA, whereas they were not observed in the free var-AGGU simulation. Since the internal mobility was suppressed upon complex formation with ARS, the entropy loss in the wild-type was larger than that of the aptamer. We therefore concluded that the sequence difference in the variable loop caused the difference in the internal mobility of the tertiary core region tRNAs and led to the difference in the affinity to ARS through the entropy term.     

Molecular dynamics simulations and free energy calculations of netropsin and distamycin binding to an AAAAA DNA binding site

Molecular dynamics simulations have been performed on netropsin in two different charge states and on distamycin binding to the minor groove of the DNA duplex d(CGCGAAAAACGCG)·d(CGCGTTTTTCGCG). The relative free energy of binding of the two non-covalently interacting ligands was calculated using the thermodynamic integration method and reflects the experimental result. From 2 ns simulations of the ligands free in solution and when bound to DNA, the mobility and the hydrogen-bonding patterns of the ligands were studied, as well as their hydration. It is shown that even though distamycin is less hydrated than netropsin, the loss of ligand–solvent interactions is very similar for both ligands. The relative mobilities of the ligands in their bound and free forms indicate a larger entropic penalty for distamycin when binding to the minor groove compared with netropsin, partially explaining the lower binding affinity of the distamycin molecule. The detailed structural and energetic insights obtained from the molecular dynamics simulations allow for a better understanding of the factors determining ligand–DNA binding.     

Computer Simulation of Macromolecules

Computer simulation techniques are now an essential part of modern structural molecular biology. They are used in many different ways in order to study the conformation, dynamics and interactions of proteins and nucleic acids. In this paper, I shall review several of these applications and then focus on three specific areas, namely the conformation and dynamics of proteins including the use of free energy perturbation methods to study mutant proteins, the conformation and dynamics of DNA and DNA-drug complexes, and the use of computers with parallel architectures. Although simulation of molecules as large and complex as proteins and nucleic acids may be considered a grand challenge in itself, there are even greater challenges for the future.     

Molecular dynamic simulations of cisplatin- and oxaliplatin-d(GG) intrastand cross-links reveal differences in their conformational dynamics

Mismatch repair proteins, DNA damage-recognition proteins and translesion DNA polymerases discriminate between Pt-GG adducts containing cis-diammine ligands (formed by cisplatin (CP) and carboplatin) and trans-RR-diaminocyclohexane ligands (formed by oxaliplatin (OX)) and this discrimination is thought to be important in determining differences in the efficacy, toxicity and mutagenicity of these platinum anticancer agents. We have postulated that these proteins recognize differences in conformation and/or conformational dynamics of the DNA containing the adducts. We have previously determined the NMR solution structure of OX-DNA, CP-DNA and undamaged duplex DNA in the 5'-d(CCTCAGGCCTCC)-3' sequence context and have shown the existence of several conformational differences in the vicinity of the Pt-GG adduct. Here we have used molecular dynamics simulations to explore differences in the conformational dynamics between OX-DNA, CP-DNA and undamaged DNA in the same sequence context. Twenty-five 10 ns unrestrained fully solvated molecular dynamics simulations were performed starting from two different DNA conformations using AMBER v8.0. All 25 simulations reached equilibrium within 4 ns, were independent of the starting structure and were in close agreement with previous crystal and NMR structures. Our data show that the cis-diammine (CP) ligand preferentially forms hydrogen bonds on the 5' side of the Pt-GG adduct, while the trans-RR-diaminocyclohexane (OX) ligand preferentially forms hydrogen bonds on the 3' side of the adduct. In addition, our data show that these differences in hydrogen bond formation are strongly correlated with differences in conformational dynamics, specifically the fraction of time spent in different DNA conformations in the vicinity of the adduct, for CP- and OX-DNA adducts. We postulate that differential recognition of CP- and OX-GG adducts by mismatch repair proteins, DNA damage-recognition proteins and DNA polymerases may be due, in part, to differences in the fraction of time that the adducts spend in a conformation favorable for protein binding.     

Minor groove deformability of DNA: a molecular dynamics free energy simulation study

The conformational deformability of nucleic acids can influence their function and recognition by proteins. A class of DNA binding proteins including the TATA box binding protein binds to the DNA minor groove, resulting in an opening of the minor groove and DNA bending toward the major groove. Explicit solvent molecular dynamics simulations in combination with the umbrella sampling approach have been performed to investigate the molecular mechanism of DNA minor groove deformations and the indirect energetic contribution to protein binding. As a reaction coordinate, the distance between backbone segments on opposite strands was used. The resulting deformed structures showed close agreement with experimental DNA structures in complex with minor groove-binding proteins. The calculated free energy of minor groove deformation was approximately 4-6 kcal mol(-1) in the case of a central TATATA sequence. A smaller equilibrium minor groove width and more restricted minor groove mobility was found for the central AAATTT and also a significantly ( approximately 2 times) larger free energy change for opening the minor groove. The helical parameter analysis of trajectories indicates that an easier partial unstacking of a central TA versus AT basepair step is a likely reason for the larger groove flexibility of the central TATATA case.     

Thermodynamic computational approach to capture molecular recognition in the binding of different inhibitors to the DNA gyrase B subunit from Escherichia coli

DNA gyrase subunit B, that catalyzes the hydrolysis of ATP, is an attractive target for the development of antibacterial drugs. This work is intended to rationalize molecular recognition at DNA gyrase B enzyme – inhibitor binding interface through the evaluation of different scoring functions in finding the correct pose and scoring properly 50 Escherichia coli DNA Gyrase B inhibitors belonging to five different classes. Improving the binding free energy calculation accuracy is further attempted by using rescoring schemes after short molecular dynamic simulations of the obtained docked complexes. These data are then compared with the corresponding experimental enzyme activity data. The results are analyzed from a structural point of view emphasizing the strengths and limitations of the techniques applied in the study.     

Computational analyses of interactions between ALK-5 and bioactive ligands: insights for the design of potential anticancer agents

Activin Receptor-Like Kinase 5 (ALK-5) is related to some types of cancer, such as breast, lung, and pancreas. In this study, we have used molecular docking, molecular dynamics simulations, and free energy calculations in order to explore key interactions between ALK-5 and six bioactive ligands with different ranges of biological activity. The motivation of this work is the lack of crystal structure for inhibitor–protein complexes for this set of ligands. The understanding of the molecular structure and the protein–ligand interaction could give support for the development of new drugs against cancer. The results show that the calculated binding free energy using MM-GBSA, MM-PBSA, and SIE is correlated with experimental data with r² = 0.88, 0.80, and 0.94, respectively, which indicates that the calculated binding free energy is in excellent agreement with experimental data. In addition, the results demonstrate that H bonds with Lys232, Glu245, Tyr249, His283, Asp351, and one structural water molecule play an important role for the inhibition of ALK-5. Overall, we discussed the main interactions between ALK-5 and six inhibitors that may be used as starting points for designing new molecules to the treatment of cancer.     

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.     

Determinative factors in inhibition of aquaporin by different pharmaceuticals: Atomic scale overview by molecular dynamics simulation

The inhibition of water permeation through aquaporins by ligands of pharmaceutical compounds is considered as a method to control the cell lifetime. The inhibition of aquaporin 1 (AQP1) by bacopaside-I and torsemide, was explored and its atomistic nature was elucidated by molecular docking and molecular dynamics (MD) simulation collectively along with Poisson-Boltzmann surface area (PBSA) method. Docking results revealed that torsemide has a lower level of docking energy in comparison with bacopaside-I at the cytoplasmic side. Furthermore, the effect of steric constraints on water permeation was accentuated. Bacopaside-I inhibits the channel properly due to the strong interaction with the channel and larger spatial volume, whereas torsemide blocks the cytoplasmic side of the channel imperfectly. The most probable active sites of AQP1 for the formation of hydrogen bonds between the inhibitor and the channel were identified by numerical analysis of the bonds. Eventually, free energy assessments indicate that binding of both inhibitors is favorable in complex with AQP1, and van der Waals interaction has an important contribution in stabilizing the complexes.     

Shear-driven rolling of DNA-adhesive microspheres

Many biologically important cell binding processes, such as the rolling of leukocytes in the vasculature, are multivalent, being mediated by large numbers of weak binding ligands. Quantitative agreement between experiments and models of rolling has been elusive and often limited by the poor understanding of the binding and unbinding kinetics of the ligands involved. Here, we present a cell-free experimental model for such rolling, consisting of polymer microspheres whose adhesion to a glass surface is mediated by ligands with well-understood force-dependent binding free energy—short complementary DNA strands. We observe robust rolling activity for certain values of the shear rate and the grafted DNA strands’ binding free energy and force sensitivity. The simulation framework developed to model leukocyte rolling, adhesive dynamics, quantitatively captures the mean rolling velocity and lateral diffusivity of the experimental particles using known values of the experimental parameters. Moreover, our model captures the velocity variations seen within the trajectories of single particles. Particle-to-particle variations can be attributed to small, plausible differences in particle characteristics. Overall, our findings confirm that state-of-the-art adhesive dynamics simulations are able to capture the complex physics of particle rolling, boding well for their extension to modeling more complex systems of rolling cells.     

Molecular recognition of H3/H4 histone tails by the tudor domains of JMJD2A: a comparative molecular dynamics simulations study

Background

Histone demethylase, JMJD2A, specifically recognizes and binds to methylated lysine residues at histone H3 and H4 tails (especially trimethylated H3K4 (H3K4me3), trimethylated H3K9 (H3K9me3) and di,trimethylated H4K20 (H4K20me2, H4K20me3)) via its tandem tudor domains. Crystal structures of JMJD2A-tudor binding to H3K4me3 and H4K20me3 peptides are available whereas the others are not. Complete picture of the recognition of the four histone peptides by the tandem tudor domains yet remains to be clarified.

Methodology/Principal Findings

We report a detailed molecular dynamics simulation and binding energy analysis of the recognition of JMJD2A-tudor with four different histone tails. 25 ns fully unrestrained molecular dynamics simulations are carried out for each of the bound and free structures. We investigate the important hydrogen bonds and electrostatic interactions between the tudor domains and the peptide molecules and identify the critical residues that stabilize the complexes. Our binding free energy calculations show that H4K20me2 and H3K9me3 peptides have the highest and lowest affinity to JMJD2A-tudor, respectively. We also show that H4K20me2 peptide adopts the same binding mode with H4K20me3 peptide, and H3K9me3 peptide adopts the same binding mode with H3K4me3 peptide. Decomposition of the enthalpic and the entropic contributions to the binding free energies indicate that the recognition of the histone peptides is mainly driven by favourable van der Waals interactions. Residue decomposition of the binding free energies with backbone and side chain contributions as well as their energetic constituents identify the hotspots in the binding interface of the structures.

Conclusion

Energetic investigations of the four complexes suggest that many of the residues involved in the interactions are common. However, we found two receptor residues that were related to selective binding of the H3 and H4 ligands. Modifications or mutations on one of these residues can selectively alter the recognition of the H3 tails or the H4 tails.