Calculation of ligand binding free energies from molecular dynamics simulations

A recently developed method for predicting binding affinities in ligand–receptor complexes, based on interaction energy averaging and conformational sampling by molecular dynamics simulation, is presented. Polar and nonpolar contributions to the binding free energy are approximated by a linear scaling of the corresponding terms in the average intermolecular interaction energy for the bound and free states of the ligand. While the method originally assumed the validity of electrostatic linear response, we show that incorporation of systematic deviations from linear response derived from free energy perturbation calculations enhances the accuracy of the approach. The method is applied to complexes of wild-type and mutant human dihydrofolate reductases with 2,4-diaminopteridine and 2,4-diaminoquinazoline inhibitors. It is shown that a binding energy accuracy of about 1 kcal/mol is attainable even for multiply ionized compounds, such as methotrexate, for which electrostatic interactions energies are very large. © 1998 John Wiley & Sons, Inc. Int J Quant Chem 69: 77–88, 1998     

Estimating the temperature dependence of peptide folding entropies and free enthalpies from total energies in molecular dynamics simulations

The temperature dependence of thermodynamic quantities, such as heat capacity, entropy and free enthalpy, may be obtained by using well-known equations that relate these quantities to the enthalpy of the molecular system of interest at a range of temperatures. In turn, the enthalpy of a molecular system can be estimated from molecular dynamics simulations of an appropriate model. To demonstrate this, we have investigated the temperature dependence of the enthalpy, heat capacity, entropy and free enthalpy of a system that consists of a beta-heptapeptide in methanol and have used the statistical mechanics relationships to describe the thermodynamics of the folding/unfolding equilibrium of the peptide. The results illustrate the power of current molecular simulation force fields and techniques in establishing the link between thermodynamic quantities and conformational distributions.     

Denaturation of Drew-Dickerson DNA in a high salt concentration medium: molecular dynamics simulations

We have performed molecular dynamics simulation on B-DNA duplex (CGCGAATTGCGC) at different temperatures. The DNA was immerged in a salt-water medium with 1 M NaCl concentration to investigate salt effect on the denaturation process. At each temperature, configurational entropy is estimated using the covariance matrix of atom-positional fluctuations, from which the melting temperature (T(m)) was found to be 349 K. The calculated configuration entropy for different bases shows that the melting process involves more peeling (including fraying from the ends) conformations, and therefore the untwisting of the duplex and peeling states form the transition state of the denaturation process. There is a narrow minor groove in the AATT sequence that becomes wider by increasing temperature which disappears at high temperatures, especially above the melting temperature. We have also calculated the fraction of denatured base pairs, f-curve, from which T(m) was found to be 340 K, close to experimental value of 341 K. We found that DNA at high salt concentrations has few hydrogen bonds even at temperatures higher than the T(m). Our calculations show the fact that adding salt leads to increase of T(m) and stabilization of DNA.     

Solution Structure,Mechanism of Replication,and Optimization of an Unnatural Base Pair

As part of an ongoing effort to expand the genetic alphabet for in vitro and eventual in vivo applications, we have synthesized a wide variety of predominantly hydrophobic unnatural base pairs and evaluated their replication in DNA. Collectively, the results have led us to propose that these base pairs, which lack stabilizing edge‐on interactions, are replicated by means of a unique intercalative mechanism. Here, we report the synthesis and characterization of three novel derivatives of the nucleotide analogue d MMO2 , which forms an unnatural base pair with the nucleotide analogue d 5SICS . Replacing the para‐methyl substituent of d MMO2 with an annulated furan ring (yielding d FMO ) has a dramatically negative effect on replication, while replacing it with a methoxy (d DMO ) or with a thiomethyl group (d TMO ) improves replication in both steady‐state assays and during PCR amplification. Thus, d TMO –d 5SICS , and especially d DMO –d 5SICS , represent significant progress toward the expansion of the genetic alphabet. To elucidate the structure–activity relationships governing unnatural base pair replication, we determined the solution structure of duplex DNA containing the parental d MMO2 –d 5SICS pair, and also used this structure to generate models of the derivative base pairs. The results strongly support the intercalative mechanism of replication, reveal a surprisingly high level of specificity that may be achieved by optimizing packing interactions, and should prove invaluable for the further optimization of the unnatural base pair.     

Net charge changes in the calculation of relative ligand‐binding free energies via classical atomistic molecular dynamics simulation

The calculation of binding free energies of charged species to a target molecule is a frequently encountered problem in molecular dynamics studies of (bio‐)chemical thermodynamics. Many important endogenous receptor‐binding molecules, enzyme substrates, or drug molecules have a nonzero net charge. Absolute binding free energies, as well as binding free energies relative to another molecule with a different net charge will be affected by artifacts due to the used effective electrostatic interaction function and associated parameters (e.g., size of the computational box). In the present study, charging contributions to binding free energies of small oligoatomic ions to a series of model host cavities functionalized with different chemical groups are calculated with classical atomistic molecular dynamics simulation. Electrostatic interactions are treated using a lattice‐summation scheme or a cutoff‐truncation scheme with Barker–Watts reaction‐field correction, and the simulations are conducted in boxes of different edge lengths. It is illustrated that the charging free energies of the guest molecules in water and in the host strongly depend on the applied methodology and that neglect of correction terms for the artifacts introduced by the finite size of the simulated system and the use of an effective electrostatic interaction function considerably impairs the thermodynamic interpretation of guest‐host interactions. Application of correction terms for the various artifacts yields consistent results for the charging contribution to binding free energies and is thus a prerequisite for the valid interpretation or prediction of experimental data via molecular dynamics simulation. Analysis and correction of electrostatic artifacts according to the scheme proposed in the present study should therefore be considered an integral part of careful free‐energy calculation studies if changes in the net charge are involved. © 2013 The Authors Journal of Computational Chemistry Published by Wiley Periodicals, Inc.     

Molecular dynamics simulations and free energy calculations on the enzyme 4‐hydroxyphenylpyruvate dioxygenase

4‐Hydroxyphenylpyruvate dioxygenase is a relevant target in both pharmaceutical and agricultural research. We report on molecular dynamics simulations and free energy calculations on this enzyme, in complex with 12 inhibitors for which experimental affinities were determined. We applied the thermodynamic integration approach and the more efficient one‐step perturbation. Even though simulations seem well converged and both methods show excellent agreement between them, the correlation with the experimental values remains poor. We investigate the effect of slight modifications on the charge distribution of these highly conjugated systems and find that accurate models can be obtained when using improved force field parameters. This study gives insight into the applicability of free energy methods and current limitations in force field parameterization. © 2011 Wiley Periodicals, Inc. J Comput Chem 2011     

Molecular Switching Behavior in Isosteric DNA Base Pairs

The structures and proton‐coupled behavior of adenine–thymine (A‐T) and a modified base pair containing a thymine isostere, adenine–difluorotoluene (A‐F), are studied in different solvents by dispersion‐corrected density functional theory. The stability of the canonical Watson–Crick base pair and the mismatched pair in various solvents with low and high dielectric constants is analyzed. It is demonstrated that A‐F base pairing is favored in solvents with low dielectric constant. The stabilization and conformational changes induced by protonation are also analyzed for the natural as well as the mismatched base pair. DNA sequences capable of changing their sequence conformation on protonation are used in the construction of pH‐based molecular switches. An acidic medium has a profound influence in stabilizing the isostere base pair. Such a large gain in stability on protonation leads to an interesting pH‐controlled molecular switch, which can be incorporated in a natural DNA tract.     

Absolute free energies of binding of peptide analogs to the HIV-1 protease from molecular dynamics simulations

The constants of binding of five peptide analogs to the active site of the HIV-1 aspartic-protease are calculated based on a novel sampling scheme that is efficient and does not introduce any approximations in addition to the energy function used to describe the system. The results agree with experiments. The squared correlation coefficient of the calculated vs. the measured values is 0.79. The sampling scheme consists of a series of molecular dynamics integrations with biases. The biases are selected based on an estimate of the probability density function of the system in a way to explore the conformational space and to reduce the statistical error in the calculated binding constants. The molecular dynamics integrations are done with a vacuum potential using a short cutoff scheme. To estimate the probability density of the simulated system, the results of the molecular dynamics integrations are combined using an extension of the weighted histogram analysis method (C. Bartels, Chem. Phys. Letters 331 (2000) 446-454). The probability density of the solvated ligand-protein system is obtained by applying a correction for the use of the short cutoffs in the simulations and by taking into account solvation with an electrostatic term and a hydrophobic term. The electrostatic part of the solvation is determined by finite difference Poisson-Boltzmann calculations; the hydrophobic part of the solvation is set proportional to the solvent accessible surface area. Setting the hydrophobic surface tension parameter equal to 8 mol(-1) K(-1) A(-2), absolute binding constants are in the muM to nM range. This is in agreement with experiments. The standard errors determined from eight repeated binding constant determinations are a factor of 14 to 411. A single determination of a binding constant is done with 499700 steps of molecular dynamics integration and 4500 finite difference Poisson-Boltzmann calculations. The simulations can be analyzed with respect to conformational changes of the active site of the HIV-1 protease or the ligands upon binding and provide information that complements experiments and can be used in the drug development process.     

Efficient calculation of relative binding free energies by umbrella sampling perturbation

An important task of biomolecular simulation is the calculation of relative binding free energies upon chemical modification of partner molecules in a biomolecular complex. The potential of mean force (PMF) along a reaction coordinate for association or dissociation of the complex can be used to estimate binding affinities. A free energy perturbation approach, termed umbrella sampling (US) perturbation, has been designed that allows an efficient calculation of the change of the PMF upon modification of a binding partner based on the trajectories obtained for the wild type reference complex. The approach was tested on the interaction of modified water molecules in aqueous solution and applied to in silico alanine scanning of a peptide‐protein complex. For the water interaction test case, excellent agreement with an explicit PMF calculation for each modification was obtained as long as no long range electrostatic perturbations were considered. For the alanine scanning, the experimentally determined ranking and binding affinity changes upon alanine substitutions could be reproduced within 0.1–2.0 kcal/mol. In addition, good agreement with explicitly calculated PMFs was obtained mostly within the sampling uncertainty. The combined US and perturbation approach yields, under the condition of sufficiently small system modifications, rigorously derived changes in free energy and is applicable to any PMF calculation. © 2014 Wiley Periodicals, Inc.     

Efficient estimation of binding free energies between peptides and an MHC class II molecule using coarse‐grained molecular dynamics simulations with a weighted histogram analysis method

We estimate the binding free energy between peptides and an MHC class II molecule using molecular dynamics (MD) simulations with the weighted histogram analysis method (WHAM). We show that, owing to its more thorough sampling in the available computational time, the binding free energy obtained by pulling the whole peptide using a coarse‐grained (CG) force field (MARTINI) is less prone to significant error induced by inadequate‐sampling than using an atomistic force field (AMBER). We further demonstrate that using CG MD to pull 3–4 residue peptide segments while leaving the remaining peptide segments in the binding groove and adding up the binding free energies of all peptide segments gives robust binding free energy estimations, which are in good agreement with the experimentally measured binding affinities for the peptide sequences studied. Our approach thus provides a promising and computationally efficient way to rapidly and reliably estimate the binding free energy between an arbitrary peptide and an MHC class II molecule. © 2017 Wiley Periodicals, Inc.     

Visualisation of inter loop tertiary base pairing and stacking interactions in an unusual slipped loop DNA structure

The conformation of an unusual slipped loop DNA structure exhibited by the sequence d(GAATTCCCGAATTC)₂ is determined using a combination of geometrical and molecular mechanics methods. This sequence is known to form a B-DNA-like duplex with the central non-complementary cytosines extruded into single stranded loop regions. The unusual feature is that the interior guanine does not pair with the cytosine across, instead, it pairs with the cytosine upstream by skipping two cytosines, leading to a slipped loop DNA structure with the loops staggered by two base pairs. The two loops, despite being very small, can fold across minor or major groove symmetrically or asymmetrically disposed, with one of the loop bases partially blocking the major or minor groove. Most interestingly, for certain conformations, the loop bases approach one another at close proximity so as to engage even in base pairing as well as base stacking interactions across the major groove. While such pairing and stacking are common in the tertiary folds of RNA, this is the first time that such an interaction is visualized in a DNA. This observation demonstrates that a W-C pair can readily be accomplished in a typical slipped loop structure postulated for DNA. Such tertiary loop interaction may prevent access to regulatory proteins across the major groove of the duplex DNA, thus providing a structure-function relation for the occurrence of slipped loop structure in DNA. Contribution no. 839 from this department     

Ab initio molecular dynamics simulation of a water-hydrogen fluoride equimolar mixture.

Hydrogen fluoride and water can be mixed in any proportion. The resulting solutions have unique acidic properties. In particular, hydrogen fluoride undergoes a weak-to-strong acidity transition with increasing concentration of HF To supplement the knowledge already obtained on dilute or moderately concentrated solutions and gas-phase aggregates, an equimolar mixture is studied here by Car-Parrinello molecular dynamics. The natures of the ions and of the complexes formed in the equimolar liquid were determined. Specifically, H3O+, H5O2+, FH-OH2, and HF2- were spontaneously obtained while only hydronium and fluoride ions pre-exist in the equimolar crystal. The behaviour of the proton in the equimolar liquid was compared with mixtures of other proportions simulated previously in an attempt to relate proton dynamics to acidity. In the same way, the behaviour of HF2- was also examined. In this case, proton localization and transfer appeared to be driven by the fluctuating environment of the solvated ion.     

Comparison of thermodynamic integration and Bennett acceptance ratio for calculating relative protein‐ligand binding free energies

The performances of Bennett's acceptance ratio method and thermodynamic integration (TI) for the calculation of free energy differences in protein simulations are compared. For the latter, the standard trapezoidal rule, Simpson's rule, and Clenshaw‐Curtis integration are used as numerical integration methods. We evaluate the influence of the number and definition of intermediate states on the precision, accuracy, and efficiency of the free energy calculations. Our results show that non‐equidistantly spaced intermediate states are in some cases beneficial for the TI methods. Using several combinations of softness parameters and the λ power dependence, it is shown that these benefits are strongly dependent on the shape of the integrand. Although TI is more user‐friendly due to its simplicity, it was found that Bennett's acceptance ratio method is the more efficient method. It is also the least dependent on the choice of the intermediate states, making it more robust than TI. © 2013 Wiley Periodicals, Inc.     

The performance of ANI-ML potentials for ligand-n(H2O) interaction energies and estimation of hydration free energies from end-point MD simulations

Here, we investigate the performance of “Accurate NeurAl networK engINe for Molecular Energies” (ANI), trained on small organic compounds, on bulk systems including non-covalent interactions and applicability to estimate solvation (hydration) free energies using the interaction between the ligand and explicit solvent (water) from single-step MD simulations. The method is adopted from ANI using the Atomic Simulation Environment (ASE) and predicts the non-covalent interaction energies at the accuracy of wb97x/6-31G(d) level by a simple linear scaling for the conformations sampled by molecular dynamics (MD) simulations of ligand-n(H₂O) systems. For the first time, we test ANI potentials' abilities to reproduce solvation free energies using linear interaction energy (LIE) formulism by modifying the original LIE equation. Our results on ~250 different complexes show that the method can be accurate and have a correlation of R² = 0.88–0.89 (MAE <1.0 kcal/mol) to the experimental solvation free energies, outperforming current end-state methods. Moreover, it is competitive to other conventional free energy methods such as FEP and BAR with 15-20 × fold reduced computational cost.     

Solvation of tetraalkylammonium chlorides in acetonitrile-water mixtures: mass spectrometry and molecular dynamics simulations.

The solvation of tetramethylammonium chloride (Me4NCl) and tetra-n-butylammonium chloride (Bu4NCl) in water-acetonitrile mixtures was investigated by mass spectrometry of clusters isolated from the solution. As far as the positive ions are concerned, clusters composed of alkylammonium ions and acetonitrile molecules only were observed, even for mixtures with high water content. In contrast, for the negative ions, clusters composed of chloride with both water and/or acetonitrile molecules were observed. For the smaller system (Me4NCl) we performed quantum chemical calculations and molecular dynamics simulations. It was found that even though water is present in the solvation shell of Me4N+, only acetonitrile has a strong electrostatic interaction with the cation. Water molecules around Me4N+ form hydrogen bonds with other water molecules, and they interact with Me4N+ mainly via dispersive interactions. These results indicate that Me4N+ behaves like a hydrophobic solute. On the other hand, the interaction of Cl- with water and acetonitrile is of comparable strength and, in both cases, the electrostatic interaction dominates. Herein we demonstrate experimentally and theoretically that positive and negative ions give rise to characteristic solvation structures in mixed solvents: even a relatively small organic cation, such as Me4N+, exhibits a hydrophobic-like solvation shell.     

Electrostatic energy analysis of 8-oxoguanine DNA lesion—molecular dynamics study

One nanosecond molecular dynamics (MD) simulation was performed for two DNA segments each composed of 30 base pairs. In one DNA segment the native guanines at nucleotides positions 17 and 19 were replaced with two 8-oxoguanines (8-oxoG) (8-oxoG is mutagenic DNA oxo-lesion). The analysis of results was focused on the electrostatic energy that is supposed to be significant factor causing the disruption of DNA base stacking in DNA duplex and may also serve as a signal toward the repair enzyme informing the presence of the lesion. The repulsive interaction between 8-oxoG and the entire DNA molecule was observed, which caused the extrahelical position of 8-oxoG (position 19). The repulsive electrostatic interaction between both 8-oxoG lesions contributed to the flipping out of one 8-oxoG and to the local instability of the lesioned DNA region. The electrostatic potential at the surface of DNA close to the lesions has more negative value than the same region on the native DNA. This electrostatic potential may signal presence of the lesion to the repair enzyme. In the simulation of native DNA segment, no significant structural changes were observed and B-DNA structure was well preserved throughout the MD simulation.     

Force field validation for nucleic acid simulations: comparing energies and dynamics of a DNA dodecamer

Important questions exist regarding the quality of force fields used in molecular dynamics (MD) simulations and their interoperable use with other available MD implementations. NAMD is one of the most efficient and scalable parallel molecular dynamics codes for large-scale biomolecular simulations in the open source domain. It is the aim of this article to analyze and compare the dynamics of a benchmark DNA dodecamer d(CTTTTGCAAAAG)2 system, including its binding to a specific drug molecule arising from the use of various simulation protocols in NAMD using Amber98, with the dynamics arising from simulations of the same dodecamer using Amber98 in the AMBER package, one of the most well-established simulation codes for nucleic acids. Based upon a set of validation benchmarks, the details of which are discussed, we find that nucleic acid simulations using NAMD give meaningful results and that the essential features of the resulting dynamics are similar to those arising from the AMBER package. This sets the stage for reliable large-scale simulations of nucleic acids using NAMD.     

Hydrogen bonding in mimics of Watson-Crick base pairs involving C-H proton donor and F proton acceptor groups: a theoretical study.

We have theoretically analyzed mimics of Watson-Crick adenine-thymine (AT) and guanine-cytosine (GC) base pairs in which N-H...O and N...H-N hydrogen bonds are replaced by N-H...F and N...H-C, respectively, by using the generalized gradient approximation of density functional theory at BP86/TZ2P. The general effect of the above substitutions is an elongation and weakening of the hydrogen bonds that hold together the base pairs. However, the precise effects depend on how many and, in particular, on which hydrogen bonds are substituted in AT and GC. Another purpose of this work is to clarify the relative importance of electrostatic attraction versus orbital interaction in the weak hydrogen bonds involved in the mimics, by using a quantitative bond-energy decomposition scheme. At variance with widespread believe, the orbital interaction component in these weak hydrogen bonds is found to contribute 34-42% of the attractive interactions and is thus of the same order of magnitude as the electrostatic component, which provides the remaining attraction component.