Some practical aspects of free energy calculations from molecular dynamics simulation

In this work, we address two critical aspects of calculation of the free energy differences in molecular systems from molecular simulations. The first aspect involves checking whether the calculated free energy difference depends significantly on the extent of perturbation used for accomplishment of a given transformation. The second aspect of interest is to verify if the sampling errors in calculating the free energy differences between the wild-type molecule and a mutated one in its free state and in a complex are similar, or not, for a finite-length dynamic simulation. The reliability of the free energy estimates obtained from molecular simulations using thermodynamic cycles depends in part on this fact. For investigating these aspects, we use a self-transformation scheme in which a transformation of a part of a molecular system into itself is considered. We perform MD simulations of DNA fragments in which a part of a specific base is subjected to such a self-transformation. Results indicate that the estimated free energy differences do not depend significantly on the extent of perturbation used to achieve the transformation. Interestingly, the variation in the cumulative free energy difference, ΔA, with the coupling parameter, λ, depends significantly on the extent of perturbation. We examine the physical basis of the observed nature of the variation of the accumulated free energy difference, ΔA, against the λ value in the case of a self-transformation. In a thermodynamic cycle, the sampling errors due to the finite-length simulation for the molecular system are found to be similar to each other for the two perturbations (free and in a complex) justifying the use of such approach in calculating ΔΔA in molecular complexes. ©1999 John Wiley & Sons, Inc. J Comput Chem 20: 877–885, 1999     

Adaptively biased molecular dynamics for free energy calculations

We present an adaptively biased molecular dynamics (ABMD) method for the computation of the free energy surface of a reaction coordinate using nonequilibrium dynamics. The ABMD method belongs to the general category of umbrella sampling methods with an evolving biasing potential and is inspired by the metadynamics method. The ABMD method has several useful features, including a small number of control parameters and an O(t) numerical cost with molecular dynamics time t. The ABMD method naturally allows for extensions based on multiple walkers and replica exchange, where different replicas can have different temperatures and/or collective variables. This is beneficial not only in terms of the speed and accuracy of a calculation, but also in terms of the amount of useful information that may be obtained from a given simulation. The workings of the ABMD method are illustrated via a study of the folding of the Ace-GGPGGG-Nme peptide in a gaseous and solvated environment.     

Computational studies on pseudorotaxanes by molecular dynamics and free energy perturbation simulations

A computational scheme that comprises the utilization of the AMBER force field with RESP charges and an explicit solvent model for acetonitrile proved to be useful for studying the structures and energetics of pseudorotaxanes of benzidine and 4,4'-biphenol with cyclobis(paraquat-p-phenylene). The scheme can be further utilized for modeling [2]rotaxanes.     

A negative cooperativity mechanism of human CYP2E1 inferred from molecular dynamics simulations and free energy calculations

Human cytochrome P450 2E1 (CYP2E1) participates in the metabolism of over 2% of all the oral drugs. A hallmark peculiar feature of this enzyme is that it exhibits a pronounced negative cooperativity in substrate binding. However the mechanism by which the negative cooperativity occurs is unclear. Here, we performed molecular dynamics simulations and free energy calculations on human CYP2E1 to examine the structural differences between the substrate-free and the enzymes with one and two aniline molecules bound. Our results indicate that although the effector substrate does not bind in the active site cavity, it still can directly interact with the active site residues of human CYP2E1. The interaction of the effector substrate with the active site leads to a reorientation of active site residues, which thereby weakens the interactions of the active substrate with this site. We also identify a conserved residue T303 that plays a crucial role in the negative cooperative binding on the short-range effects. This residue is a key factor in the positioning of substrates and in proton delivery to the active site. Additionally, a long-range effect of the effector substrate is identified in which F478 is proposed to play a key role. As located in the interface between the active and effector sites, this residue structurally links the active and effector sites and is found to play a significant role in affecting substrate access and ligand positioning within the active site. In the negative cooperative binding, this residue can decrease the interactions of the active substrate with the active site by π-π stacking which then lowers the hydroxylation activity for the active substrate. These findings are in agreement with previous experimental observations and thus provide detailed atomistic insight into the poorly understood mechanism of the negative cooperativity in human CYP2E1.     

Aggregation mechanism of polyglutamine diseases revealed using quantum chemical calculations, fragment molecular orbital calculations, molecular dynamics simulations, and binding free energy calculations

Polyglutamine (polyQ) diseases, including Huntington’s disease (HD), are caused by expansion of polyQ-encoding repeats within otherwise unrelated gene products. The aggregation mechanism of polyQ diseases, the inhibition mechanism of Congo red, and the alleviation mechanism of trehalose were proposed here based on quantum chemical calculations and molecular dynamics simulations. The calculations and simulations revealed the following. The effective molecular bonding is between glutamine (Gln) and Gln (Gln + Gln), between Gln and Congo red (Gln + Congo red), and between Gln and trehalose (Gln + trehalose). The bonding strength is −13.1 kcal/mol for Gln + Gln, −24.4 kcal/mol for Gln + Congo red, and −12.0 kcal/mol for Gln + trehalose. In the polyQ region, both the number of intermolecular Gln + Gln formations and the total calories generated by the Gln + Gln formation are proportional to the number of repetitions of Gln. We propose an aggregation mechanism whose heat generated by the intermolecular Gln + Gln formation causes the pathogeny of polyQ disease. In our aggregation mechanism, this generated heat collapses the host protein and promotes fibrillogenesis. Without contradiction, our mechanism can explain all the experimental results reported to date. Our mechanism can also explain the inhibition mechanism by Congo red as an inhibitor of polyglutamine-induced protein aggregation and the alleviation mechanism by trehalose as an alleviator of that aggregation. The inhibition mechanism by Congo red is explained by the strong interaction with Gln and by the characteristic structure of Congo red.     

Studies of chirality effect of 4-(phenylamino)-pyrrolo[2,1-f][1,2,4]triazine on p38α by molecular dynamics simulations and free energy calculations

4-(Phenylamino)-pyrrolo[2,1-f][1,2,4]triazines have been discovered as inhibitors of p38α. Experimental assays have proven that the configuration of α-Me-benzyl connected with amide at C6 is essential for the binding affinity. The S-configured inhibitor (11j) displays 80 times more potency than the R-configured one (11k). Here we investigated the mechanism how different configurations influence the binding affinity using molecular dynamics simulations, free energy calculations and free energy decomposition analysis. We found that the van der Waals interactions play the most important role in differentiating the activities between 11j and 11k with p38α. The difference of the van der Waals interactions is primarily determined by two residues, LEU108 and LEU167. Consequently stabilization of pyrrolo[2,1-f][1,2,4]triazine ring is important for the activities of inhibitors. Meanwhile we observed that the different configuration of the α-Me-benzyl group leads to the difference of binding between 11j and 11k. In conclusion, our work shows that it is feasible to analyze the chirality effect of inhibitors with different configurations by molecular dynamics simulations and free energy calculations, and provides useful information for drug design.     

Molecular dynamics simulation, free energy calculation and structure-based 3D-QSAR studies of B-RAF kinase inhibitors

(V600E)B-RAF kinase is the most frequent onco-genic protein kinase mutation in melanoma and is a promising target to treat malignant melanoma. In this work, a molecular modeling study combining QM-polarized ligand docking, molecular dynamics, free energy calculation, and three-dimensional quantitative structure-activity relationships (3D-QSAR) was performed on a series of pyridoimidazolone compounds as the inhibitors of (V600E)B-RAF kinase to understand the binding mode between the inhibitors and (V600E)B-RAF kinase and the structural requirement for the inhibiting activity. 3D-QSAR models, including CoMFA and CoMSIA, were developed from the conformations obtained by QM-polarized ligand docking strategy. The obtained models have a good predictive ability in both internal and external validation. Furthermore, molecular dynamics simulation and free energy calculations were employed to determine the detailed binding process and to compare the binding mode of the inhibitors with different activities. The binding free energies calculated by MM/PBSA gave a good correlation with the experimental biological activity. The decomposition of free energies by MM/GBSA indicates the van der Waals interaction is the major driving force for the interaction between the inhibitors and (V600E)B-RAF kinase. The hydrogen bond interactions between the inhibitors with Glu501 and Asp594 of the (V600E)B-RAF kinase help to stabilize the DFG-out conformation. The results from this study can provide some insights into the development of novel potent (V600E)B-RAF kinase inhibitors.     

Fast and flexible gpu accelerated binding free energy calculations within the amber molecular dynamics package

Alchemical free energy (AFE) calculations based on molecular dynamics (MD) simulations are key tools in both improving our understanding of a wide variety of biological processes and accelerating the design and optimization of therapeutics for numerous diseases. Computing power and theory have, however, long been insufficient to enable AFE calculations to be routinely applied in early stage drug discovery. One of the major difficulties in performing AFE calculations is the length of time required for calculations to converge to an ensemble average. CPU implementations of MD‐based free energy algorithms can effectively only reach tens of nanoseconds per day for systems on the order of 50,000 atoms, even running on massively parallel supercomputers. Therefore, converged free energy calculations on large numbers of potential lead compounds are often untenable, preventing researchers from gaining crucial insight into molecular recognition, potential druggability and other crucial areas of interest. Graphics Processing Units (GPUs) can help address this. We present here a seamless GPU implementation, within the PMEMD module of the AMBER molecular dynamics package, of thermodynamic integration (TI) capable of reaching speeds of >140 ns/day for a 44,907‐atom system, with accuracy equivalent to the existing CPU implementation in AMBER. The implementation described here is currently part of the AMBER 18 beta code and will be an integral part of the upcoming version 18 release of AMBER. © 2018 Wiley Periodicals, Inc.     

Direct calculation of the crystal-melt interfacial free energy via molecular dynamics computer simulation

We review our recent work on the direct calculation of the interfacial free energy, gamma, of the crystal-melt interface via molecular dynamics computer simulation for a number of model systems. The value of gamma as a function of crystal orientation is determined using a thermodynamic integration technique employing moving cleaving walls [Phys. Rev. Lett. 2000, 85, 4751]. The calculation is sufficiently precise to resolve the small anisotropy in gamma, which is crucial in determining the kinetics and morphology of dendritic growth. We report values of gamma for the hard-sphere and Lennard-Jones systems, as well as recent results on the series of inverse-power potentials. For the inverse sixth-, seventh-, and eighth-power systems, we determine gamma for both fcc and bcc crystal structures. For these systems, the bcc-melt gamma is lower than that for fcc crystals by about 25%, consistent with recent experiments and computer simulations on fcc-forming systems that show preferential formation of bcc nuclei in the initial stages of crystallization. In addition, we show that our results give a molecular interpretation of Turnbull's rule, which is the empirical relationship between gamma and the enthalpy of fusion.     

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.     

麋鹿朊蛋白结构稳定性的分子动力学研究

朊蛋白病是一种致命且具有高度传染性的神经退行性疾病.糜鹿是目前已经报道的哺乳类动物中较易发生朊蛋白病的物种之一.作者使用分子动力学和操控式分子动力学模拟相结合的方法对糜鹿正常朊蛋白的结构稳定性进行了研究.发现了麋鹿朊蛋白结构中的不稳定结构域分布以及热动力学性质,揭示了糜鹿朊蛋白稳定性的分子结构基础以及力学特征.     

Insights into the phosphoryl-transfer mechanism of cAMP-dependent protein kinase from quantum chemical calculations and molecular dynamics simulations

To investigate the molecular details of the phosphoryl-transfer mechanism catalyzed by cAMP-dependent protein kinase, we performed quantum mechanical (QM) calculations on a cluster model of the active site and molecular dynamics (MD) simulations of a ternary complex of the protein with Mg(2)ATP and a 20-residue peptide substrate. Overall, our theoretical results confirm the participation of the conserved aspartic acid, Asp(166), as an acid/base catalyst in the reaction mechanism catalyzed by protein kinases. The MD simulation shows that the contact between the nucleophilic serine side chain and the carboxylate group of Asp(166) is short and dynamically stable, whereas the QM study indicates that an Asp(166)-assisted pathway is structurally and energetically feasible and is in agreement with previous experimental results.     

Minimum free energy pathways and free energy profiles for conformational transitions based on atomistic molecular dynamics simulations

An efficient method for the calculation of minimum free energy pathways and free energy profiles for conformational transitions is presented. Short restricted perturbation-targeted molecular dynamics trajectories are used to generate an approximate free energy surface. Approximate reaction pathways for the conformational change are constructed from one-dimensional line segments on this surface using a Monte Carlo optimization. Accurate free energy profiles are then determined along the pathways by means of one-dimensional adaptive umbrella sampling simulations. The method is illustrated by its application to the alanine "dipeptide." Due to the low computational cost and memory demands, the method is expected to be useful for the treatment of large biomolecular systems.     

3D-QSAR,docking, molecular dynamics simulation and free energy calculation studies of some pyrimidine derivatives as novel JAK3 inhibitors

Janus kinase 3 (JAK3) is a promising drug target for the treatment of inflammatory diseases, autoimmune disorders, organ transplant rejection and various cancers. In the present study, 3D-QSAR, docking, MD simulation and MM/PBSA studies were performed on a series of pyrimidine-based JAK3 inhibitors. A reliable COMSIA (q² = 0.717 and r² = 0.986) model was developed and validated using external validation test set, bootstrapping, progressive scrambling and r_m² metrics analyses. Structural requirements identified through contour maps of the model were strategically utilized to computationally design 170 novel JAK3 inhibitors with improved potency. Docking studies were performed on the selected data set and newly designed compounds to show their binding mode and to identify important interacting residues inside the active site of JAK3. In addition, docking results of the selected designed compounds inside the active sites of JAK1, JAK2 and TYK2 indicated their JAK3 selectivity. MD simulation (100 ns) on the docked complex of compound 28 (one of highly active compounds of the data set) assisted in the further exploration of the binding interactions. Some crucial residues like Lys830 (glycine-rich loop), Val836, Ala853, Leu905 (hinge region), Cys909, Asn954, Leu956 and Ala966 were identified. Hydrogen bond interactions with hinge residue Leu905 were critical for the binding of JAK3 inhibitors. Additionally, MM/PBSA calculation provided the binding free energy of the compound 28. Newly designed molecules showed promising results in the preliminary in silico ADMET evaluations. Outcomes of the study can further be exploited to develop potent JAK3 inhibitors.     

Steric and allosteric effects of fatty acids on the binding of warfarin to human serum albumin revealed by molecular dynamics and free energy calculations

Human serum albumin (HSA) binds with drugs and fatty acids (FAs). This study was initiated to elucidate the relationship between the warfarin binding affinity of HSA and the positions of bound FA molecules. Molecular dynamics simulations of 11 HSA-warfarin-myristate complexes were performed. HSA-warfarin binding free energy was then calculated for each of the complexes by the molecular mechanics-Poisson-Boltzmann surface area (MM-PBSA) method. The results indicated that the magnitude of the binding free energy was smaller in HSA-warfarin complexes that had 4 or more myristate molecules than in complexes with no myristate molecules. The unfavorable effect on the HSA-warfarin binding affinity was caused sterically by the binding of a myristate molecule to the FA binding site closest to the warfarin binding site. On the other hand, the magnitude of HSA-warfarin binding free energy was largest when 3 myristate molecules were bound to the high-affinity sites. The strongest HSA-warfarin binding was attributable to favorable entropic contribution related to larger atomic fluctuations of the amino acid residues at the warfarin binding site. In the binding of 2 myristate molecules to the sites with the highest and second-highest affinities, allosteric modulation that enhanced electrostatic interactions between warfarin and some of the amino acid residues around the warfarin binding site was observed. This study clarified the structural and energetic properties of steric/allosteric effects of FAs on the HSA-warfarin binding affinity and illustrated the approach to analyze protein-ligand interactions in situations such that multiple ligands bind to the other sites of the protein.     

The effect of molecular chirality on the incidence of twisted smectic A phases

Four chiral chloro-esters were synthesized in order to examine systematically the effect that molecular chirality has on the temperature range and incidence of the twisted smectic A phase or the twist grain boundary phase as it has become known. We find that when the motion of the chiral centre is restricted by rotational hindrance, thereby increasing the chirality of the system, the twisted S_A phase range is increased. At low chirality the twisted S_A phase disappears altogether. Furthermore, it is shown that the transition temperatures of the chiral compounds are lower than their racemic analogues.     

Calculation of the entropy and free energy of peptides by molecular dynamics simulations using the hypothetical scanning molecular dynamics method

Hypothetical scanning (HS) is a method for calculating the absolute entropy S and free energy F from a sample generated by any simulation technique. With this approach each sample configuration is reconstructed with the help of transition probabilities (TPs) and their product leads to the configuration's probability, hence to the entropy. Recently a new way for calculating the TPs by Monte Carlo (MC) simulations has been suggested, where all system interactions are taken into account. Therefore, this method--called HSMC--is in principle exact where the only approximation is due to insufficient sampling. HSMC has been applied very successfully to liquid argon, TIP3P water, self-avoiding walks on a lattice, and peptides. Because molecular dynamics (MD) is considered to be significantly more efficient than MC for a compact polymer chain, in this paper HSMC is extended to MD simulations as applied to peptides. Like before, we study decaglycine in vacuum but for the first time also a peptide with side chains, (Val)(2)(Gly)(6)(Val)(2). The transition from MC to MD requires implementing essential changes in the reconstruction process of HSMD. Results are calculated for three microstates, helix, extended, and hairpin. HSMD leads to very stable differences in entropy TDeltaS between these microstates with small errors of 0.1-0.2 kcal/mol (T=100 K) for a wide range of calculation parameters with extremely high efficiency. Various aspects of HSMD and plans for future work are discussed.     

Free energy perturbation and molecular dynamics calculations of copper binding to azurin

Free energy perturbation/molecular dynamics simulations have been carried out on copper/azurin systems calculating the binding affinities of copper (II) ion to azurin either in the native or in the unfolded state. In order to test the validity of the strategy adopted for the calculations and to establish what force field is suitable for these kinds of calculations, three different force fields, AMBER, CVFF, and CFF, have been alternatively used for the calculations and the results have been compared with experimental data obtained by spectroscopic titrations of copper (II)/azurin solutions and denaturation experiments. Our findings have pointed out that only CFF gives satisfactory results, thus providing a reliable tool for copper binding simulations in copper protein.     

Extra precision docking,free energy calculation and molecular dynamics studies on glutamic acid derivatives as MurD inhibitors

The binding modes of well known MurD inhibitors have been studied using molecular docking and molecular dynamics (MD) simulations. The docking results of inhibitors 1-30 revealed similar mode of interaction with Escherichia coli-MurD. Further, residues Thr36, Arg37, His183, Lys319, Lys348, Thr321, Ser415 and Phe422 are found to be important for inhibitors and E. coli-MurD interactions. Our docking procedure precisely predicted crystallographic bound inhibitor 7 as evident from root mean square deviation (0.96 Å). In addition inhibitors 2 and 3 have been successfully cross-docked within the MurD active site, which was pre-organized for the inhibitor 7. Induced fit best docked poses of 2, 3, 7 and 15/2Y1O complexes were subjected to 10 ns MD simulations to determine the stability of the predicted binding conformations. Induce fit derived docked complexes were found to be in a state of near equilibrium as evident by the low root mean square deviations between the starting complex structure and the energy minimized final average MD complex structures. The results of molecular docking and MD simulations described in this study will be useful for the development of new MurD inhibitors with high potency.