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.     

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.     

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.     

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

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

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.     

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.     

Crystallization of helical oligomers with chirality selection. I. A molecular dynamics simulation for bare helix

Helical polymers often exhibit pronounced chirality recognition during crystallization. By molecular dynamics simulation, we have already shown that the helical polymers crystallize with or without marked chirality selection depending on structural details of the polymer molecules. We have there classified the helical polymers into two categories: the bare helices made of only backbone atoms which show rather tolerant chirality selection, and the general helices with large side groups showing strict chirality recognition. Polymer crystallization is in general largely hampered and retarded by slow dynamics of the entangled chains, and therefore short helical oligomers are very suitable models for studying the chiral crystallization. We here report on molecular simulations of crystallization in the bare helical oligomer molecules by the use of Monte Carlo and molecular dynamics simulations. First we confirm the low temperature chiral crystal phase and the reversible order-disorder transition. We also observe frequent inversions of the helical sense, and the helix reversal defects propagating along the chains. Then we investigate crystallization from the melt into the chiral crystal phase. We find that the crystallization rate depends very sensitively on the degree of undercooling. The crystallization is found to be the first order transition that conforms well to the traditional picture of crystal growth in small molecules. Even when the crystallization directly into the chiral crystal phase is conducted, marked chirality selections are not observed at the early stage of crystallization; the chains adhere to the crystal surfaces selecting their helical senses rather at random resulting in racemic crystallites. The isothermal crystallization for a sufficiently long time, however, yields lamellar crystals composed of well-developed chiral domains, the growth of which seems to be accomplished through the transition back into the ordered chiral crystal phase.     

Modelling and molecular dynamics simulation studies on a hexagonal glycolipid assembly

The hexagonal columnar phase (H_I) of an aqueous formulation of octyl β-glucoside with 67 % lipid content was modelled and 15-ns molecular dynamics simulation was performed. Initial investigations on the aggregation size led to good correlation of simulation and experimental d-spacing for a 12 molecule cylinder core. The corresponding hexagonal phase was stable over the entire simulation time and provided conclusive local density profiles. Hydrogen bonding analyses showed only minor differences in the bonding profile between the hexagonal and a previously reported micellar phase. However, the glycoside interaction decreases with increasing curvature, i.e. from a lamellar assembly over the hexagonal phase to the micelle, while the opposite behaviour applies for interactions with water. A view into the water dynamics revealed an anisotropic-correlated diffusion process with higher mobility along the cylinder axes than perpendicular to them.     

Structural basis of specific binding between Aurora A and TPX2 by molecular dynamics simulations

In the present study, the impacts of G198N and W128F mutations on the recognition between Aurora A and targeting protein of Xenopus kinesin-like protein 2 (TPX2) were investigated using molecular dynamics (MD) simulations, free energy calculations, and free energy decomposition analysis. The predicted binding free energy of the wild-type complex is more favorable than those of three mutants, indicating that both single and double mutations are unfavorable for the Aurora A and TPX2 binding. It is also observed that the mutations alternate the binding pattern between Aurora A and TPX2, especially the downstream of TPX2. An intramolecular hydrogen bond between the atom OD of Asp11(TPX2) and the atom HE1 of Trp34(TPX2) disappear in three mutants and thus lead to the instability of the secondary structure of TPX2. The combination of different molecular modeling techniques is an efficient way to understand how mutation has impacts on the protein-protein binding and our work gives valuable information for the future design of specific peptide inhibitors for Aurora A.     

Studies of the mechanism of selectivity of protein tyrosine phosphatase 1B (PTP1B) bidentate inhibitors using molecular dynamics simulations and free energy calculations

Bidentate inhibitors of protein tyrosine phosphatase 1B (PTP1B) are considered as a group of ideal inhibitors with high binding potential and high selectivity in treating type II diabetes. In this paper, the binding models of five bidentate inhibitors to PTP1B, TCPTP, and SHP-2 were investigated and compared by using molecular dynamics (MD) simulations and free energy calculations. The binding free energies were computed using the Molecular Mechanics/Poisson-Boltzmann Surface Area (MM/PBSA) methodology. The calculation results show that the predicted free energies of the complexes are well consistent with the experimental data. The Molecular Mechanics/Generalized Born Surface Area (MM/GBSA) free energy decomposition analysis indicates that the residues ARG24, ARG254, and GLN262 in the second binding site of PTP1B are essential for the high selectivity of inhibitors. Furthermore, the residue PHE182 close to the active site is also important for the selectivity and the binding affinity of the inhibitors. According to our analysis, it can be concluded that in most cases the polarity of the portion of the inhibitor that binds to the second binding site of the protein is positive to the affinity of the inhibitors while negative to the selectivity of the inhibitors. We expect that the information we obtained here can help to develop potential PTP1B inhibitors with more promising specificity.     

The effect of hydrogen bonding on oligoleucine structure in water: A molecular dynamics simulation study

Molecular dynamics simulations were conducted in order to improve our understanding of the forces that determine polyleucine chains conformations and govern polyleucine self-assembly in aqueous solutions. Simulations of 10 repeat unit oligoleucine in aqueous solution were performed using the optimized potential for liquid simulations (OPLS) - all atom force field using the canonical ensemble for a minimum of 1.3 ns. These simulations provided information on conformations, chain collapse and intermolecular aggregation. Simulations indicate that single isotactic oligoleucine chains in dilute solution assume tightly packed, regular hairpin conformations while atactic oligoleucine assumes a much less regular and less compact structure. The regular, compact collapsed isotactic chain exhibited a greater degree of intramolecular hydrogen bonding and an increased level of hydrophobic t-butyl functional group aggregation compared to the atactic chain. This occurs at the expense of reduced leucine-water hydrogen bonding.     

Backbone motions of free and pheromone-bound major urinary protein I studied by molecular dynamics simulation

Molecular motions of free and pheromone-bound mouse major urinary protein I, previously investigated by NMR relaxation, were simulated in 30 ns molecular dynamics (MD) runs. The backbone flexibility was described in terms of order parameters and correlation times, commonly used in the NMR relaxation analysis. Special attention was paid to the effect of conformational changes on the nanosecond time scale. Time-dependent order parameters were determined in order to separate motions occurring on different time scales. As an alternative approach, slow conformational changes were identified from the backbone torsion angle variances, and "conformationally filtered" order parameters were calculated for well-defined conformation states. A comparison of the data obtained for the free and pheromone-bound protein showed that some residues are more rigid in the bound form, but a larger portion of the protein becomes more flexible upon the pheromone binding. This finding is in general agreement with the NMR results. The higher flexibility observed on the fast (fs-ps) time scale was typically observed for the residues exhibiting higher conformational freedom on the ns time scale. An inspection of the hydrogen bond network provided a structural explanation for the flexibility differences between the free and pheromone-bound proteins in the simulations.