New approach to free energy of solvation applying continuum models to molecular dynamics simulation

A new approach to the calculation of the free energy of solvation from trajectories obtained by molecular dynamics simulation is presented. The free energy of solvation is computed as the sum of three contributions originated at the cavitation of the solute by the solvent, the solute-solvent nonpolar (repulsion and dispersion) interactions, and the electrostatic solvation of the solute. The electrostatic term is calculated based on ideas developed for the broadly used continuum models, the cavitational contribution from the excluded volume by the Claverie-Pierotti model, and the Van der Waals term directly from the molecular dynamics simulation. The proposed model is tested for diluted aqueous solutions of simple molecules containing a variety of chemically important functions: methanol, methylamine, water, methanethiol, and dichloromethane. These solutions were treated by molecular dynamics simulations using SPC/E water and the OPLS force field for the organic molecules. Obtained free energies of solvation are in very good agreement with experimental data.     

Coupled semiempirical quantum mechanics and molecular mechanics (QM/MM) calculations on the aqueous solvation free energies of ionized molecules

The aqueous solvation free energies of ionized molecules were computed using a coupled quantum mechanical and molecular mechanical (QM/MM) model based on the AM1, MNDO, and PM3 semiempirical molecular orbital methods for the solute molecule and the TIP3P molecular mechanics model for liquid water. The present work is an extension of our model for neutral solutes where we assumed that the total free energy is the sum of components derived from the electrostatic/polarization terms in the Hamiltonian plus an empirical “nonpolar” term. The electrostatic/polarization contributions to the solvation free energies were computed using molecular dynamics (MD) simulation and thermodynamic integration techniques, while the nonpolar contributions were taken from the literature. The contribution to the electrostatic/polarization component of the free energy due to nonbonded interactions outside the cutoff radii used in the MD simulations was approximated by a Born solvation term. The experimental free energies were reproduced satisfactorily using variational parameters from the vdW terms as in the original model, in addition to a parameter from the one-electron integral terms. The new one-electron parameter was required to account for the short-range effects of overlapping atomic charge densities. The radial distribution functions obtained from the MD simulations showed the expected H-bonded structures between the ionized solute molecule and solvent molecules. We also obtained satisfactory results by neglecting both the empirical nonpolar term and the electronic polarization of the solute, i.e., by implementing a nonpolarization model. ©1999 John Wiley & Sons, Inc. J Comput Chem 20: 1028–1038, 1999     

Base stacking in cytosine dimer. A comparison of correlated ab initio calculations with three empirical potential models and density functional theory calculations

Ab initio MP2/6-31G* interaction energies were calculated for more than 80 geometries of stacked cytosine dimer. Diffuse polarization functions were used to properly cover the dispersion energy. The results of ab initio calculations were compared with those obtained from three electrostatic empirical potential models, constructed as the sum of a Lennard-Jones potential (covering dispersion and repulsion contributions) and the electrostatic term. Point charges and point multipoles of the electrostatic term were also obtained at the MP2/6-31G* level of theory. The point charge MEP model (atomic charges derived from molecular electrostatic potential) satisfactorily reproduced the ab initio data. Addition of π-charges localized below and above the cytosine plane did not affect the calculated energies. The model employing the distributed multipole analysis gave worse agreement with the ab initio data than the MEP approach. The MP2 MEP charges were also derived using larger sets of atomic orbitals: cc-pVDZ, 6-311 + G(2d, p), and aug-cc-pVDZ. Differences between interaction energies calculated using these three sets of point charges and the MP2/6-31G* charges were smaller than 0.8 kcal/mol. The correlated ab initio calculations were also compared with the density functional theory (DFT) method. DFT calculations well reproduced the electrostatic part of interaction energy. They also covered some nonelectrostatic short-range effects which were not reproduced by the empirical potentials. The DFT method does not include the dispersion energy. This energy, approximated by an empirical term, was therefore added to the DFT interaction energy. The resulting interaction energy exhibited an artifact secondary minimum for a 3.9-4.0 vertical separation of bases. This defect is inherent in the DFT functionals, because it is not observed for the Hartree-Fock + dispersion interaction energy.© 1996 John Wiley & Sons, Inc.     

AGBNP: an analytic implicit solvent model suitable for molecular dynamics simulations and high-resolution modeling

We have developed an implicit solvent effective potential (AGBNP) that is suitable for molecular dynamics simulations and high-resolution modeling. It is based on a novel implementation of the pairwise descreening Generalized Born model for the electrostatic component and a new nonpolar hydration free energy estimator. The nonpolar term consists of an estimator for the solute-solvent van der Waals dispersion energy designed to mimic the continuum solvent solute-solvent van der Waals interaction energy, in addition to a surface area term corresponding to the work of cavity formation. AGBNP makes use of a new parameter-free algorithm to calculate the scaling coefficients used in the pairwise descreening scheme to take into account atomic overlaps. The same algorithm is also used to calculate atomic surface areas. We show that excellent agreement is achieved for the GB self-energies and surface areas in comparison to accurate, but much more expensive, numerical evaluations. The parameter-free approach used in AGBNP and the sensitivity of the AGBNP model with respect to large and small conformational changes makes the model suitable for high-resolution modeling of protein loops and receptor sites as well as high-resolution prediction of the structure and thermodynamics of protein-ligand complexes. We present illustrative results for these kinds of benchmarks. The model is fully analytical with first derivatives and is computationally efficient. It has been incorporated into the IMPACT molecular simulation program.     

Binding of cationic and neutral phenanthridine intercalators to a DNA oligomer is controlled by dispersion energy: quantum chemical calculations and molecular mechanics simulations

Correlated ab initio as well as semiempirical quantum chemical calculations and molecular dynamics simulations were used to study the intercalation of cationic ethidium, cationic 5-ethyl-6-phenylphenanthridinium and uncharged 3,8-diamino-6-phenylphenanthridine to DNA. The stabilization energy of the cationic intercalators is considerably larger than that of the uncharged one. The dominant energy contribution with all intercalators is represented by dispersion energy. In the case of the cationic intercalators, the electrostatic and charge-transfer terms are also important. The DeltaG of ethidium intercalation to DNA was estimated at -4.5 kcal mol(-1) and this value agrees well with the experimental result. Of six contributions to the final free energy, the interaction energy value is crucial. The intercalation process is governed by the non-covalent stacking (including charge-transfer) interaction while the hydrogen bonding between the ethidium amino groups and the DNA backbone is less important. This is confirmed by the evaluation of the interaction energy as well as by the calculation of the free energy change. The intercalation affects the macroscopic properties of DNA in terms of its flexibility. This explains the easier entry of another intercalator molecule in the vicinity of an existing intercalation site.     

Molecular dynamics studies of ion distributions for DNA duplexes and DNA clusters: salt effects and connection to DNA melting

We present extensive molecular dynamics simulations of the ion distributions for DNA duplexes and DNA clusters using the Amber force field with implicit water. The distribution of ions and the electrostatic energy of ions around an isolated DNA duplex and clusters of DNA duplexes in different salt (NaCl) concentrations over the range 0.2-1.0 mol/L are determined on the basis of the simulation results. Using the electrostatic energy profile, we determine a local net charge fraction phi, which is found to increase with increasing of salt concentration. For DNA clusters containing two DNA duplexes (DNA pair) or four DNA duplexes, phi increases as the distance between the duplexes decreases. Combining this result with experimental results for the dependence of the DNA melting temperature on bulk salt concentration, we conclude that for a pair of DNA duplexes the melting temperature increases by 5-10 K for interaxis separations of 25-40 A. For a cluster of four DNA duplexes, an even larger melting temperature increase should occur. We argue that this melting temperature increase in dense DNA clusters is responsible for the cooperative melting mechanism in DNA-linked nanoparticle aggregates and DNA-linked polymer aggregates.     

Theoretical calculation of the coiled-coil stability in water in the context of its possible use as a molecular rack

The coiled-coil stability and rigidity may be of importance for molecular electronics (electronically bistable molecules). The coiled-coil binding free energy has been calculated using molecular dynamics (MD). The energy has been computed as a difference of the appropriate free energies; derived for the coiled-coil and isolated alpha-helices separately. All MD simulations have been performed using an explicit model of the solvent, whereas the continuum solvent approach has been applied to analyze the MD trajectories. The computed stability of the coiled-coil is of the order of -87 kcal/mol, i.e., about -1.2 kcal/mol per amino acid residue, and arises mainly from the electrostatic interactions and hydrophobic effect. The entropy term has been roughly estimated to be of the order of -22 kcal/mol. This assures that coiled-coil polypeptides may be used as a stable molecular scaffolding.     

平面五元水分子簇的量子化学研究

选用Gaussian03的B3LYP/6-31G(d,p)、DMol3的BLYP/DNP和deMon的BLYP/TZVP等方法计算了甲烷水合物(结构-1)中平面五元水分子簇的结合能和氢键能,作了基组重叠误差(BSSE)和色散能(dispersion)的修正,估算了次级相互作用的贡献.在DMol3程序中使用了大型数值基组DNP,将基组重叠误差降至最低.在Gaussi-an03的B3LYP/6-31G(d,p)计算中,采用平衡法(Counterpoise)校正基组重叠误差.两种计算方法给出了一致的结果,证实了在使用6-31G(d,p)基组时,一对水分子在平衡距离的基组重叠误差高达8 kJ/mol.为估算色散能的贡献,使用了新近发展的包含色散能的密度泛函的DFT程序deMon计算了五元水分子簇.用多种方法计算出了经基组重叠误差和色散能修正的五元水分子簇的分子间结合能和氢键能的较为精确的势能超曲面,为甲烷和其他气体水合物的分子动力学模拟提供了依据.     

Recognition of RNA by amide modified backbone nucleic acids: molecular dynamics simulations of DNA-RNA hybrids in aqueous solution

Thermodynamic and structural properties of a chemically modified DNA-RNA hybrid in which a phosphodiester linkage is replaced by a neutral amide-3 linkage (3'-CH(2)-CONH-5') were investigated using UV melting experiments, molecular dynamics simulations in explicit water, and continuum solvent models. van't Hoff analysis of the experimental UV melting curves suggests that the significant increase of the thermodynamic stability of a 15-mer DNA-RNA with seven alternated amide-3 modifications (+11 degrees C) is mainly due to an increased binding enthalpy. To further evaluate the origin in the observed affinities differences, the electrostatic contribution to the binding free energy was calculated by solving the Poisson-Boltzmann equation numerically. The nonelectrostatic contribution was estimated as the product of a hydrophobic surface tension coefficient and the surface area that is buried upon double strand formation. Structures were taken from 10 ns molecular dynamics simulations computed in a consistent fashion using explicit solvent, counterions, and the particle-mesh Ewald procedure. The present preliminary thermodynamic study suggests that the favorable binding free energy of the amide-3 DNA single strand to the complementary RNA is equally driven by electrostatic and nonpolar contributions to the binding compared to their natural analogues. In addition, molecular dynamics simulations in explicit water were performed on an amide-3 DNA single strand and the corresponding natural DNA. Results from the conformations cluster analysis of the simulated amide-3 DNA single strand ensembles suggest that the 25% of the population sampled within 10 ns has a pre-organized conformation where the sugar C3' endo pucker is favored at the 3'-flanking nucleotides. These structural and thermodynamic features contribute to the understanding of the observed increased affinities of the amide-3 DNA-RNA hybrids at the microscopic level.     

Binding of Cationic and Neutral Phenanthridine Intercalators to a DNA Oligomer Is Controlled by Dispersion Energy: Quantum Chemical Calculations and Molecular Mechanics Simulations

Correlated ab initio as well as semiempirical quantum chemical calculations and molecular dynamics simulations were used to study the intercalation of cationic ethidium, cationic 5‐ethyl‐6‐phenylphenanthridinium and uncharged 3,8‐diamino‐6‐phenylphenanthridine to DNA. The stabilization energy of the cationic intercalators is considerably larger than that of the uncharged one. The dominant energy contribution with all intercalators is represented by dispersion energy. In the case of the cationic intercalators, the electrostatic and charge‐transfer terms are also important. The ΔG of ethidium intercalation to DNA was estimated at ?4.5 kcal mol^?1 and this value agrees well with the experimental result. Of six contributions to the final free energy, the interaction energy value is crucial. The intercalation process is governed by the non‐covalent stacking (including charge‐transfer) interaction while the hydrogen bonding between the ethidium amino groups and the DNA backbone is less important. This is confirmed by the evaluation of the interaction energy as well as by the calculation of the free energy change. The intercalation affects the macroscopic properties of DNA in terms of its flexibility. This explains the easier entry of another intercalator molecule in the vicinity of an existing intercalation site.     

Optimization of partial multicanonical molecular dynamics simulations applied to an alanine dipeptide in explicit water solvent

The partial multicanonical algorithm for molecular dynamics and Monte Carlo simulations samples a wide range of an important part of the potential energy. Although it is a strong technique for structure prediction of biomolecules, the choice of the partial potential energy has not been optimized. In order to find the best choice, partial multicanonical molecular dynamics simulations of an alanine dipeptide in explicit water solvent were performed with 15 trial choices for the partial potential energy. The best choice was found to be the sum of the electrostatic, Lennard-Jones, and torsion-angle potential energies between solute atoms. In this case, the partial multicanonical simulation sampled all of the local-minimum free-energy states of the P(II), C(5), α(R), α(P), α(L), and C states and visited these states most frequently. Furthermore, backbone dihedral angles ? and ψ rotated very well. It is also found that the most important term among these three terms is the electrostatic potential energy and that the Lennard-Jones term also helps the simulation to overcome the steric restrictions. On the other hand, multicanonical simulation sampled all of the six states, but visited these states fewer times. Conventional canonical simulation sampled only four of the six states: The P(II), C(5), α(R), and α(P) states.     

Density functional theory for efficient ab initio molecular dynamics simulations in solution

We present a density functional for first-principles molecular dynamics simulations that includes the electrostatic effects of a continuous dielectric medium. It allows for numerical simulations of molecules in solution in a model polar solvent. We propose a smooth dielectric model function to model solvation into water and demonstrate its good numerical properties for total energy calculations and constant energy molecular dynamics.     

Re-examining the properties of the aqueous vapor-liquid interface using dispersion corrected density functional theory

First-principles molecular dynamics simulations, in which the forces are computed from electronic structure calculations, have great potential to provide unique insight into structure, dynamics, electronic properties, and chemistry of interfacial systems that is not available from empirical force fields. The majority of current first-principles simulations are driven by forces derived from density functional theory with generalized gradient approximations to the exchange-correlation energy, which do not capture dispersion interactions. We have carried out first-principles molecular dynamics simulations of air-water interfaces employing a particular generalized gradient approximation to the exchange-correlation functional (BLYP), with and without empirical dispersion corrections. We assess the utility of the dispersion corrections by comparison of a variety of structural, dynamic, and thermodynamic properties of bulk and interfacial water with experimental data, as well as other first-principles and force field-based simulations.     

Nature of adhesion of condensed organic films on platinum by first-principles simulations

Understanding the nature of the adhesion of an organic liquid on a metal surface is of paramount importance for elucidating the stability and chemical reactivity at these complex interfaces. However, to date, the morphology, layering and chemical properties at organic liquid metal interfaces have been rarely known. Using semi-empirical dispersion corrected density functional theory calculations and ab initio molecular dynamics simulations, we show that carbon tetrachloride and ethanol films confined to a platinum surface alter their intrinsic properties and exhibit interfacial reactivity. A few interface carbon tetrachloride (ethanol) molecules adsorb dissociatively (molecularly) on platinum thanks to the surrounding medium. The adsorption strength of the interfacial molecules is consequently increased in the condensed phase as compared to the gas phase. This remarkable effect is rationalized by an interaction energy decomposition model and an electrostatic potential analysis.     

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.     

Tryptophan side chain electrostatic interactions determine edge-to-face vs parallel-displaced tryptophan side chain geometries in the designed beta-hairpin "trpzip2"

The interaction geometries of the four tryptophan (Trp) side chains in the 12-residue designed beta-hairpin trpzip2 are investigated using all-atom explicit-solvent molecular dynamics simulations. The experimentally observed edge-to-face (EtF) pairwise interaction geometries are stable on a time scale of 10 ns. However, removing the electrostatic multipoles of the Trp side chains while retaining the dipoles of the side chains' NH moieties induces a conformational change to a geometry in which three of the four side chains interact in a parallel-displaced (PD) manner. Free energy simulations of the Etf to PD conformational change reveal that, with the side chain multipole moments intact (+MP), the EtF conformation is preferred by 5.79 kcal/mol. Conversely, with only the dipole moments of the side chain NH moieties intact (-MP), the PD conformation's free energy is more favorable by 1.71 kcal/mol. In contrast to energetic similarities for Trp side chain-water electrostatic and Trp side chain-Trp side chain and Trp side chain-water van der Waals, +MP Trp side chain-Trp side chain electrostatic interactions are more favorable by 4.21 kcal/mol in the EtF conformation, while in the -MP case the EtF and PD conformations' Trp side chain-Trp side chain electrostatic energies are nearly identical. The results highlight the importance of electrostatic multipole moments in determining aromatic-aromatic interaction geometries in aqueous biomolecular systems and argue for the inclusion of this physics in simplified models used for protein-ligand docking and protein structure prediction, possibly through a truncated Coulomb term between aromatic moieties.     

A hybrid all-atom/coarse grain model for multiscale simulations of DNA

Hybrid simulations of molecular systems, which combine all-atom (AA) with simplified (or coarse grain, CG) representations, propose an advantageous alternative to gain atomistic details on relevant regions while getting profit from the speedup of treating a bigger part of the system at the CG level. Here we present a reduced set of parameters derived to treat a hybrid interface in DNA simulations. Our method allows us to forthrightly link a state-of-the-art force field for AA simulations of DNA with a CG representation developed by our group. We show that no modification is needed for any of the existing residues (neither AA nor CG). Only the bonding parameters at the hybrid interface are enough to produce a smooth transition of electrostatic, mechanic and dynamic features in different AA/CG systems, which are studied by molecular dynamics simulations using an implicit solvent. The simplicity of the approach potentially permits us to study the effect of mutations/modifications as well as DNA binding molecules at the atomistic level within a significantly larger DNA scaffold considered at the CG level. Since all the interactions are computed within the same classical Hamiltonian, the extension to a quantum/classical/coarse-grain multilayer approach using QM/MM modules implemented in widely used simulation packages is straightforward.