Cutoff Errors in the Ewald Summation Formulae for Point Charge Systems

Abstract

Closed formulae for both real and reciprocal space parts of cutoff errors in the Ewald summation method in cubic periodic boundary conditions are derived. Such estimates are useful in tuning parameters in molecular simulations. Errors in both the electrostatic energy and forces are considered. The estimates apply to a disordered configuration of point charges and, with some limitations, also to point-charge molecular models. The accuracy of our estimates is tested and confirmed using simulated configurations of two systems (molten salt and diethylether) under a variety of conditions.     

Hydration effects on the electrostatic potential around tuftsin

The electrostatic potential and component dielectric constants from molecular dynamics (MD) trajectories of tuftsin, a tetrapeptide with the amino acid sequence Thr–Lys–Pro–Arg in water and in saline solution are presented. The results obtained from the analysis of the MD trajectories for the total electrostatic potential at points on a grid using the Ewald technique are compared with the solution to the Poisson–Boltzmann (PB) equation. The latter was solved using several sets of dielectric constant parameters. The effects of structural averaging on the PB results were also considered. Solute conformational mobility in simulations gives rise to an electrostatic potential map around the solute dominated by the solute monopole (or lowest order multipole). The detailed spatial variation of the electrostatic potential on the molecular surface brought about by the compounded effects of the distribution of water and ions close to the peptide, solvent mobility, and solute conformational mobility are not qualitatively reproducible from a reparametrization of the input solute and solvent dielectric constants to the PB equation for a single structure or for structurally averaged PB calculations. Nevertheless, by fitting the PB to the MD electrostatic potential surfaces with the dielectric constants as fitting parameters, we found that the values that give the best fit are the values calculated from the MD trajectories. Implications of using such field calculations on the design of tuftsin peptide analogues are discussed. © 1999 John Wiley & Sons, Inc. Biopoly 50: 133–143, 1999     

Decomposition analysis of free energy profile for Hsp90-ADP association

The functional cycle of heat shock protein 90 (Hsp90) is driven and inhibited by the association/dissociation of ligand molecules. In order to understand the molecular mechanism of the association of N-terminal domain of Hsp90 (N-Hsp90) and its ligand molecule, it is necessary to investigate which part in the target system promotes or inhibits the association of N-Hsp90 and its ligand molecule. We apply the decomposition analysis for the association free energy of N-Hsp90 and ADP. The mean force calculated by thermodynamic integration method combined with molecular dynamic simulations is divided into the contributions from molecules in the target system. Van der Waals interaction of the solvent water molecules strongly stabilises the association. Three lysine residues on the surface of the N-Hsp90 pull ADP toward the binding pocket of N-Hsp90. This study elucidates the association process of ADP from the bulk region to the binding pocket of the N-terminal domain Hsp90. This approach is applicable to elucidate the association process of biomolecules.     

Modeling loop reorganization free energies of acetylcholinesterase: a comparison of explicit and implicit solvent models

The treatment of hydration effects in protein dynamics simulations varies in model complexity and spans the range from the computationally intensive microscopic evaluation to simple dielectric screening of charge-charge interactions. This paper compares different solvent models applied to the problem of estimating the free-energy difference between two loop conformations in acetylcholinesterase. Molecular dynamics (MD) simulations were used to sample potential energy surfaces of the two basins with solvent treated by means of explicit and implicit methods. Implicit solvent methods studied include the generalized Born (GB) model, atomic solvation potential (ASP), and the distance-dependent dieletric constant. By using the linear response approximation (LRA), the explicit solvent calculations determined a free-energy difference that is in excellent agreement with the experimental estimate, while rescoring the protein conformations with GB or the Poisson equation showed inconsistent and inferior results. While the approach of rescoring conformations from explicit water simulations with implicit solvent models is popular among many applications, it perturbs the energy landscape by changing the solvent contribution to microstates without conformational relaxation, thus leading to non-optimal solvation free energies. Calculations applying MD with a GB solvent model produced results of comparable accuracy as observed with LRA, yet the electrostatic free-energy terms were significantly different due to optimization on a potential energy surface favored by an implicit solvent reaction field. The simpler methods of ASP and the distance-dependent scaling of the dielectric constant both produced considerable distortions in the protein internal free-energy terms and are consequently unreliable.     

Calculation of conformational free energies by confinement simulations in explicit water with implicit desolvation

Abstract

The confinement method is a robust and conceptually simple free energy simulation method that allows the calculation of conformational free energy differences between highly dissimilar states. Application of the method to explicitly solvated systems requires a multi-stage simulation protocol for the calculation of desolvation free energies. Here we show that these desolvation free energies can be readily obtained from an implicit treatment, which is simpler and less costly. The accuracy and robustness of this protocol was shown by the calculation of conformational free energy differences of a series of explicitly solvated test systems. Given the accuracy and ease by which these free energy differences were obtained, the confinement method is promising for the treatment of conformational changes in large and complex systems.     

Molecular Dynamics Simulations of a Calmodulin-peptide Complex in Solution

Abstract

We have performed an 4-ns MD simulation of calmodulin complexed with a target peptide in explicit water, under realistic conditions of constant temperature and pressure, in the presence of a physiological concentration of counterions and using Ewald summation to avoid truncation of long-range electrostatic forces. During the simulation the system tended to perform small fluctuations around a structure similar to, but somewhat looser than the starting crystal structure. The calmodulin-peptide complex was quite rigid and did not exhibit any large amplitude domain motions such as previously seen in apo- and calcium-bound calmodulin. We analyzed the calmodulin-peptide interactions by calculating buried surface areas, CHARMM interaction energies and continuum model interaction free energies. In the trajectory, the protein surface area buried by contact with the peptide is 1373 Ų, approximately evenly divided between the calmodulin N-terminal, C-terminal and central linker regions. A majority of this buried surface, 803 ·A², comes from nonpolar residues, in contrast to the protein as a whole, for which the surface is made up of mostly polar and charged groups. Our continuum calculations indicate that the largest favorable contribution to pep- tide binding comes from burial of molecular surface upon complex formation. Electrostatic contributions are favorable but smaller in the trajectory structures, and actually unfavorable for binding in the crystal structure. Since nonpolar groups make up most of buried surface of the protein, our calculations suggest that the hydrophobic effect is the main driving force for binding the helical peptide to calmodulin, consistent with thermodynamic analysis of experimental data. Besides the burial of nonpolar surface area, secondary contributions to peptide binding come from burial of polar surface and electrostatic interactions. In the nonpolar interactions a crucial role is played by the nine methionines of calmodulin. In the electrostatic interactions the negatively charged protein residues and positively charged peptide residues play a dominant role.     

On the Ewald Artifacts in Computer Simulations. The Test-Case of the Octaalanine Peptide With Charged Termini

Abstract

The treatment of electrostatic interactions in molecular simulations is of fundamental importance. Ewald and related methods are being increasingly used to the detriment of cutoff schemes, which are known to produce several artifacts. A potential drawback of the Ewald method is the spatial periodicity that is imposed to the system, which could produce artifacts when applied in the simulation of liquids. In this work we analyze the octaalanine peptide with charged termini in explicit solvent, for which severe effects due to the use of Ewald sums were predicted using continuum electrostatics. Molecular Dynamics simulations for a total of 158 nanoseconds were performed in cells of different sizes. From the comparison of the results of different system sizes, no significant periodicity-induced artifacts were observed. It is argued that in current biomolecular simulations, the incomplete sampling is likely to affect the results to a larger extent than the artifacts induced by the use of Ewald sums.     

Computing free energies using nested sampling-based approaches

Abstract

Nested sampling (NS) has emerged as a powerful statistical mechanical sampling technique to compute the partition function of atomic and molecular systems. From the partition function all thermodynamic quantities can be calculated in absolute terms, including absolute free energies and entropies. In this article, we provide a brief overview of NS within a Bayesian context, as well as overviews of how NS is used to compute the partition functions and thermodynamic quantities in the canonical and isothermal-isobaric ensembles. Then we introduce a new scheme, Coupling Parameter Path Nested Sampling, to estimate the free energy difference between two systems with different potential energy functions. The method uses a NS simulation to traverse the same path through phase space as would be covered in traditional coupling parameter-based methods such as thermodynamic integration and perturbation approaches. We demonstrate the new method with two case studies and confirm its accuracy by comparison to conventional methods, including Widom test particle insertion and thermodynamic integration. The proposed method provides a powerful alternative to traditional coupling parameter-based free energy simulation methods.     

Solvation free energies of nucleic acid bases in ionic liquids

The solvation free energies of five nucleic acid bases in [C_nbim]Br (where n = 2, 4, 6) ionic liquids (ILs) were computed using the Bennett acceptance ratio (BAR) method employing molecular dynamics simulations. The computed free energies using BAR were in agreement with other methods. The large and negative predicted free energies of the bases in ILs indicated that the bases were better solvated in the ILs rather than in water. Hydrogen bonding interactions between polar sites of the bases and ILs’ ions significantly contributed to the solvation mechanism.     

Standards for molecular dynamics modelling and simulation of relaxation

An attempt is made to formulate a set of requirements for simulation and modelling of relaxation in dense media. Each requirement is illustrated by examples of numerical simulation of particles with different types of interaction given by soft-sphere, Lennard–Jones, embedded atom method or Coulomb potential. The approaches developed are expected to be universal for some classes of relaxation processes in liquids, fluids, crystals and plasmas.     

Exploring the interaction between D-xylose isomerase and D-xylose by free energy calculation

The numerical quadrature thermodynamic integration method is used to investigate enzyme-substrate interaction of D-xylose isomerase. A screening function for the coulombic interaction is introduced into the simulation to correct the effect of finite cutoff radius for the non-bonded interaction. The binding free energy difference for D-xylose with D-xylose isomerase and its N184D mutant has been calculated, and the result 3.9 ± 1.2 kJ/mol agrees well with experimental data of 4.38 kJ/mol. In addition, the structure and dynamics of enzyme-substrate complex were simulated for mutant and wild-type enzyme, respectively. Analysis of the structures and intramolecular interactions of the complexes were found to be valuable for understanding the reaction mechanism of the enzyme D-xylose isomerase. © Wiley-Liss, Inc.     

Dielectric characterization of neutralized and nonneutralized chitosan upon drying

Isothermal dielectric loss spectra of neutralized and nonneutralized chitosan were acquired in successive runs from -130 degrees C up to increasing final temperatures, in a frequency range between 20 Hz and 1 MHz. Essentially, three relaxation processes were detected in the temperature range covered: (i) a beta-wet process, detected when the sample has a higher water content that vanishes after heating to 150 degrees C; (ii) a beta process, which is located at temperatures below 0 degrees C, becoming better defined and maintaining its location after annealing at 150 degrees C independently of the protonation state of the amino side group; and (iii) a sigma process that deviates to higher temperatures with drying, being more mobile in the nonneutralized form. Moreover, in dried neutralized chitosan, a fourth process was detected in the low frequency side of the secondary beta process that diminishes after annealing. Whether this process is a distinct relaxation of the dried polymer or a deviated beta-wet process due to the loss of water residues achieved by annealing is not straightforward. Only beta and sigma processes persist after annealing at 150 degrees C. The changes in molecular mobility upon drying of these two relaxation processes were evaluated.     

Docking study and free energy simulation of the complex between p53 DNA-binding domain and azurin

Molecular interaction between p53 tumor suppressor and the copper protein azurin (AZ) has been demonstrated to enhance p53 stability and hence antitumoral function, opening new perspectives in cancer treatment. While some experimental work has provided evidence for AZ binding to p53, no crystal structure for the p53-AZ complex was solved thus far. In this work the association between AZ and the p53 DNA-binding domain (DBD) was investigated by computational methods. Using a combination of rigid-body protein docking, experimental mutagenesis information, and cluster analysis 10 main p53 DBD-AZ binding modes were generated. The resulting structures were further characterized by molecular dynamics (MD) simulations and free energy calculations. We found that the highest scored docking conformation for the p53 DBD-AZ complex also yielded the most favorable free energy value. This best three-dimensional model for the complex was validated by using a computational mutagenesis strategy. In this structure AZ binds to the flexible L(1) and s(7)-s(8) loops of the p53 DBD and stabilizes them through protein-protein tight packing interactions, resulting in high degree of both surface matching and electrostatic complementarity.     

Molecular dynamics simulation studies of cytochrome b5 from outer mitochondrial and microsomal membrane

Two forms of cytochrome b(5) have been identified, associated with the outer membrane of liver mitochondria (OM cyt b(5)) and with the membrane of the endoplasmic reticulum (microsomal, Mc cyt b(5)). These proteins have very similar structures, but differ significantly in physical properties, with the OM cyt b(5) exhibiting a more negative reduction potential, higher stability, and stronger interactions with the heme. We perform molecular dynamics simulations to probe the structures and fluctuations of the two proteins in solution, to help explain the observed physical differences. We find that the structures of the two proteins, highly similar in the crystal, differ in position of a surface loop involving residues 49-51 in solution. Hydrophobic residues Ala-18, Ile-32, Leu-36, and Leu-47 tend to cluster together on the surface of rat OM cyt b(5), blocking water access to the protein interior. In bovine Mc cyt b(5), two of these positions, Ser-18 and Arg-47, are occupied by hydrophilic residues. This leads to breaking the hydrophobic cluster and allowing the protein to occupy a more open conformation. A measure of this structural transition is the opening of a cleft on the protein surface, which is 5 A wider in the OM cyt b(5) simulation compared to the Mc form. The OM protein also appears to have a more compact hydrophobic core in its beta-sheet region. These effects may be used to explain observed stability differences between the two proteins.     

Evaluating the suitability of RNA intervention mechanism exerted by some flavonoid molecules against dengue virus MTase RNA capping site: a molecular docking,molecular dynamics simulation,and binding free energy study

Abstract

Dengue virus (DENV) is one of the most dangerous mosquito-borne human pathogens known to the mankind. Currently, no vaccines or standard therapy is avaliable to treate DENV infection. This makes the drug development against DENV more significant and challenging. The MTase domain of DENV RNA RdRp NS5 is a promising drug target, because this domain hosts the RNA capping process of DENV RNA to escape from human immune system. In the present study, we have analysed the RNA intervention mechanism exerted by flavoniod molecules against NS5 MTase RNA capping site by using molecular docking, molecular dynamics simulation and the binding free energy calculations. The results from the docking analysis confirmed that the RNA intervention mecanism is exerted by the quercetagetin (QGN) molecule with all necessary intermolecular interactions and high binding affinity. Notably, QGN forms strong hydrogen bonding interactions with Asn18, Leu20 and Ser150 residues and π???π stacking interaction with Phe25 residue. The apo and QGN bound NS5 MTase and QGN-NS5 MTase complex were used for MD simulation. The results of MD simulation reveal that the RMSD and RMSF values of QGN-MTase complex have increased on comparing the apo protein due to the effect of ligand binding. The binding free energy calulation includes prediction of total binding free energy of ligand-protein complex and per-residue free energy decomposition. The QGN binding to NS5 MTase affects it’s native motion, this result is found from Principal component analysis.

Communicated by Ramaswamy H. Sarma     

Continuous metadynamics in essential coordinates as a tool for free energy modelling of conformational changes

Modelling of conformational changes in biopolymers is one of the greatest challenges of molecular biophysics. Metadynamics is a recently introduced free energy modelling technique that enhances sampling of configurational (e.g. conformational) space within a molecular dynamics simulation. This enhancement is achieved by the addition of a history-dependent bias potential, which drives the system from previously visited regions. Discontinuous metadynamics in the space of essential dynamics eigenvectors (collective motions) has been proposed and tested in conformational change modelling. Here, we present an implementation of two continuous formulations of metadynamics in the essential subspace. The method was performed in a modified version of the molecular dynamics package GROMACS. These implementations were tested on conformational changes in cyclohexane, alanine dipeptide (terminally blocked alanine, Ace-Ala-Nme) and SH3 domain. The results illustrate that metadynamics in the space of essential coordinates can accurately model free energy surfaces associated with conformational changes. Figure The conformational free energy surface of cyclohexane in the space of the two most intensive collective motions.

---

     

Molecular modeling study on the dynamical structural features of human smoothened receptor and binding mechanism of antagonist LY2940680 by metadynamics simulation and free energy calculation

Background

The smoothened (SMO) receptor, one of the Class F G protein coupled receptors (GPCRs), is an essential component of the canonical hedgehog signaling pathway which plays a key role in the regulation of embryonic development in animals. The function of the SMO receptor can be modulated by small-molecule agonists and antagonists, some of which are potential antitumour agents. Understanding the binding mode of an antagonist in the SMO receptor is crucial for the rational design of new antitumour agents.

Methods

Molecular dynamics (MD) simulation and dynamical network analysis are used to study the dynamical structural features of SMO receptor. Metadynamics simulation and free energy calculation are employed to explore the binding mechanism between the antagonist and SMO receptor.

Results

The MD simulation results and dynamical network analysis show that the conserved KTXXXW motif in helix VIII has strong interaction with helix I. The α-helical extension of transmembrane 6 (TM6) is detected as part of the ligand-binding pocket and dissociation pathway of the antagonist. The metadynamics simulation results illustrate the binding mechanism of the antagonist in the pocket of SMO receptor, and free energy calculation shows the antagonist needs to overcome about 38 kcal/mol of energy barrier to leave the binding pocket of SMO receptor.

Conclusions

The unusually long TM6 plays an important role on the binding behavior of the antagonist in the pocket of SMO receptor.

General significance

The results can not only profile the binding mechanism between the antagonist and Class F GPCRs, but also supply the useful information for the rational design of a more potential small molecule antagonist bound to SMO receptor.     

Dynamic ejection behaviour of water molecules passing through a nano-aperture nozzle

This study used molecular dynamics (MD) simulation to investigate the passage of water molecules through a composite graphene/Au nano-nozzle. Our focus was on the degree to which system temperature, extrusion speed, and nozzle diameter affect jet dynamics and the associated transient phenomena. Our findings show that high pressure and spatial confinement cause the nanojet from a small nozzle diameter (1.0?nm) to bend and twist, whereas the jets from a nozzle with a diameter of 1.5?nm present columns of greater stability. At 100?K, the H₂O nanojet froze at the outlet of the nozzle in the form of condensed icicles. At 500?K, the H₂O nanojet formed a loose spray and gaseous clusters. High extrusion speed of 55.824?m/s produced recirculating flow downstream from the nanojet with the appearance of an erupting volcano, which further prompted the jet column to thicken. Lower extrusion speeds produced jets with flow velocity insufficient to overcome the capillary force at the outlet of the nozzle, which subsequently manifests as unstable fluctuations in the flow rate.

• HIGHLIGHTS

• Water molecules through a composite graphene/Au nano-nozzle forming a nanojet is investigated.

• High pressure and spatial confinement cause the nanojet from a small nozzle diameter (≤1.0?nm) to bend and twist.

• High extrusion speed (≧55.824?m/s) produced recirculating flow downstream from the nanojet.

• Figure abstract: Schematic of the H₂O nano-jet through a nano-nozzle of graphene/Au

     

Accuracy and precision of NMR relaxation experiments and MD simulations for characterizing protein dynamics

Model-free parameters obtained from nuclear magnetic resonance (NMR) relaxation experiments and molecular dynamics (MD) simulations commonly are used to describe the intramolecular dynamical properties of proteins. To assess the relative accuracy and precision of experimental and simulated model-free parameters, three independent data sets derived from backbone ¹⁵N NMR relaxation experiments and two independent data sets derived from MD simulations of Escherichia coli ribonuclease HI are compared. The widths of the distributions of the differences between the order parameters for pairs of NMR data sets are congruent with the uncertainties derived from statistical analyses of individual data sets; thus, current protocols for analyzing NMR data encapsulate random uncertainties appropriately. Large differences in order parameters for certain residues are attributed to systematic differences between samples for intralaboratory comparisons and unknown, possibly magnetic field-dependent, experimental effects for interlaboratory comparisons. The widths of distributions of the differences between the order parameters for two NMR sets are similar to widths of distributions for an NMR and an MD set or for two MD sets. The linear correlations between the order parameters for an MD set and an NMR set are within the range of correlations observed between pairs of NMR sets. These comparisons suggest that the NMR and MD generalized order parameters for the backbone amide N—H bond vectors are of comparable accuracy for residues exhibiting motions on a fast time scale (<100 ps). Large discrepancies between NMR and MD order parameters for certain residues are attributed to the occurrence of “rare” motional events over the simulation trajectories, the disruption of an element of secondary structure in one of the simulations, and lack of consensus among the experimental data sets. Consequently, (easily detectable) severe distortions of local protein structure and infrequent motional events in MD simulations appear to be the most serious artifacts affecting the accuracy and precision, respectively, of MD order parameters relative to NMR values. In addition, MD order parameters for motions on a fast (<100 ps) timescale are more precisely determined than their NMR counterparts, thereby permitting more detailed dynamic characterization of biologically important residues by MD simulation than is sometimes possible by experimental methods. Proteins 28:481–493, 1997. © 1997 Wiley-Liss, Inc.     

Quantum-chemical and empirical calculations on Phospholipids VIII. The electrostatic potential from isolated molecules up to layer systems

Using 1-6-12 empirical functions with a solvent-averaged electrostatic contribution $\text{q}{}_{\mathrm{I}}\text{q}{}_{\mathrm{j}}\text{ε(r}{}_{\mathrm{Ij}}\text{)}\text{× r}{}_{\mathrm{Ij}}$ and electrostatic potentials from CNDO-type wavefunctions, the development of specific interactions of ions visualized by the molecular electrostatic potential of PO₄-group containing molecules was studied. Going from single molecules to monolayers made up of 37 head groups of phosphatidylcholine (PC) or phosphatidylethanolamine (PE) for quantum-chemical calculations, or of 23 head groups for empirical calculations we found decreasing potential minima. Only the inclusion of the screening effect of water, simulated by a distance dependent dielectric constant, ε(r), gives an explanation of stereospecific interactions of model membranes with ions. This finding can be compared with results of simulation calculations on water structure above a PE head group layer.