Analysis of the activation mechanism of the guinea-pig Histamine H<Subscript>1</Subscript>-receptor

The Histamine H(1)-receptor (H1R), belonging to the amine receptor-class of family A of the G-protein coupled receptors (GPCRs) gets activated by agonists. The consequence is a conformational change of the receptor, which may involve the binding-pocket. So, for a good prediction of the binding-mode of an agonist, it is necessary to have knowledge about these conformational changes. Meanwhile some experimental data about the structural changes of GPCRs during activation exist. Based on homology modeling of the guinea-pig H1R (gpH1R), using the crystal structure of bovine rhodopsin as template, we performed several MD simulations with distance restraints in order to get an inactive and an active structure of the gpH1R. The calculations led to a Phe6.44/Trp6.48/Phe6.52-switch and linearization of the proline kinked transmembrane helix VI during receptor activation. Our calculations showed that the Trp6.48/Phe6.52-switch induces a conformational change in Phe6.44, which slides between transmembrane helices III and VI. Additionally we observed a hydrogen bond interaction of Ser3.39 with Asn7.45 in the inactive gpH1R, but because of a counterclockwise rotation of transmembrane helix III Ser3.39 establishes a water-mediated hydrogen bond to Asp2.50 in the active gpH1R. Additionally we simulated a possible mechanism for receptor activation with a modified LigPath-algorithm.     

Role of the N-Terminal Transmembrane Helix Contacts in the Activation of FGFR3

Fibroblast growth factor receptor 3 (FGFR3) is a member of receptor tyrosine kinases, which is involved in skeletal cell growth, differentiation, and migration. FGFR3 transduces biochemical signals from the extracellular ligand-binding domain to the intracellular kinase domain through the conformational changes of the transmembrane (TM) helix dimer. Here, we apply generalized replica exchange with solute tempering method to wild type (WT) and G380R mutant (G380R) of FGFR3. The dimer interface in G380R is different from WT and the simulation results are in good agreement with the solid-state nuclear magnetic resonance (NMR) spectroscopy. TM helices in G380R are extended more than WT, and thereby, G375 in G380R contacts near the N-termini of the TM helix dimer. Considering that both G380R and G375C show the constitutive activation, the formation of the N-terminal contacts of the TM helices can be generally important for the activation mechanism. © 2019 Wiley Periodicals, Inc.     

Interaction specific binding hotspots in Endonuclease colicin-immunity protein complex from MD simulations

The binding of Endonuclease colicin 9 (E9) by Immunity protein 9 (Im9) was found to involve some hotspots from helix III of Im9 on protein-protein interface that contribute the dominant binding energy to the complex.In the current work,MD simulations of the WT and three hotspot mutants (D51A,Y54A and Y55A of Im9) of the E9-Im9 complexes were carried out to investigate specific interaction mechanisms of these three hotspot residues.The changes of binding energy between the WT and mutants of the complex were computed by the MM/PBSA method using a polarized force field and were in excellent agreement with experiment values,verifying that these three residues were indeed hotspots of the binding complex.Energy decomposition analysis revealed that binding by D51 to E9 was dominated by electrostatic interaction due to the presence of the carboxyl group of Asp51 which hydrogen bonds to K89.For binding by hotspots Y54 and Y55,van der Waals interaction from the aromatic side chain of tyrosine provided the dominant interaction.For comparison,calculation by using the standard (nonpolarizable) AMBER99SB force field produced binding energy changes from these mutations in opposite direction to the experimental observation.Dynamic hydrogen bond analysis showed that conformations sampled from MD simulation in the standard AMBER force field were distorted from the native state and they disrupted the inter-protein hydrogen bond network of the protein-protein complex.The current work further demonstrated that electrostatic polarization plays a critical role in modulating protein-protein binding.     

Solvophobic and steric effects of side groups on polymer folding: molecular modeling studies of amine-functionalized m-poly(phenyleneethynylene) foldamers in aqueous solution

The folding behavior of five different amine-functionalized m-poly(phenyleneethynylene) (m-PPE) oligomers containing 24 phenyl rings (12 residues, where a residue includes 2 phenyl rings) in water was examined by using a combination of molecular dynamics (MD) and replica exchange molecular dynamics (REMD) simulation techniques. The REMD method employed the highly parallelized GROMACS MD software and a modified OPLS-AA force field to simulate 44 replicas of each solvated system in parallel, with temperatures ranging from 300 to 577 K. Our results showed that the REMD method was more effective in predicting the helical conformation of the m-PPE in water, from an extended structure, than canonical MD methods in the same simulation time. Furthermore, we observed from canonical MD simulations of the explicitly solvated helical m-PPEs at 300 K that the radius of gyration, average helix inner diameter, and average helix pitch of the helical structure all pass through a minima when the side group is R = OC(2)H(5) as R is changed from R = H through OC(4)H(9).     

Generation of a homology model of the human histamine H<Subscript>3</Subscript> receptor for ligand docking and pharmacophore-based screening

The human histamine H(3) receptor (hH(3)R) is a G-protein coupled receptor (GPCR), which modulates the release of various neurotransmitters in the central and peripheral nervous system and therefore is a potential target in the therapy of numerous diseases. Although ligands addressing this receptor are already known, the discovery of alternative lead structures represents an important goal in drug design. The goal of this work was to study the hH(3)R and its antagonists by means of molecular modelling tools. For this purpose, a strategy was pursued in which a homology model of the hH(3)R based on the crystal structure of bovine rhodopsin was generated and refined by molecular dynamics simulations in a dipalmitoylphosphatidylcholine (DPPC)/water membrane mimic before the resulting binding pocket was used for high-throughput docking using the program GOLD. Alternatively, a pharmacophore-based procedure was carried out where the alleged bioactive conformations of three different potent hH(3)R antagonists were used as templates for the generation of pharmacophore models. A pharmacophore-based screening was then carried out using the program Catalyst. Based upon a database of 418 validated hH(3)R antagonists both strategies could be validated in respect of their performance. Seven hits obtained during this screening procedure were commercially purchased, and experimentally tested in a [(3)H]N(alpha)-methylhistamine binding assay. The compounds tested showed affinities at hH(3)R with K ( i ) values ranging from 0.079 to 6.3 muM.     

The Molecular Mechanism of P2Y1 Receptor Activation

Human purinergic G protein‐coupled receptor P2Y₁ (P2Y₁R) is activated by adenosine 5′‐diphosphate (ADP) to induce platelet activation and thereby serves as an important antithrombotic drug target. Crystal structures of P2Y₁R revealed that one ligand (MRS2500) binds to the extracellular vestibule of this GPCR, whereas another (BPTU) occupies the surface between transmembrane (TM) helices TM2 and TM3. We introduced a total of 20 μs all‐atom long‐timescale molecular dynamic (MD) simulations to inquire why two molecules in completely different locations both serve as antagonists while ADP activates the receptor. Our results indicate that BPTU acts as an antagonist by stabilizing extracellular helix bundles leading to an increase of the lipid order, whereas MRS2500 blocks signaling by occupying the ligand binding site. Both antagonists stabilize an ionic lock within the receptor. However, binding of ADP breaks this ionic lock, forming a continuous water channel that leads to P2Y₁R activation.     

Structure-based prediction of subtype selectivity of histamine H3 receptor selective antagonists in clinical trials

Histamine receptors (HRs) are excellent drug targets for the treatment of diseases, such as schizophrenia, psychosis, depression, migraine, allergies, asthma, ulcers, and hypertension. Among them, the human H(3) histamine receptor (hH(3)HR) antagonists have been proposed for specific therapeutic applications, including treatment of Alzheimer's disease, attention deficit hyperactivity disorder (ADHD), epilepsy, and obesity. However, many of these drug candidates cause undesired side effects through the cross-reactivity with other histamine receptor subtypes. In order to develop improved selectivity and activity for such treatments, it would be useful to have the three-dimensional structures for all four HRs. We report here the predicted structures of four HR subtypes (H(1), H(2), H(3), and H(4)) using the GEnSeMBLE (GPCR ensemble of structures in membrane bilayer environment) Monte Carlo protocol, sampling ～35 million combinations of helix packings to predict the 10 most stable packings for each of the four subtypes. Then we used these 10 best protein structures with the DarwinDock Monte Carlo protocol to sample ～50?000 × 10(20) poses to predict the optimum ligand-protein structures for various agonists and antagonists. We find that E206(5.46) contributes most in binding H(3) selective agonists (5, 6, 7) in agreement with experimental mutation studies. We also find that conserved E5.46/S5.43 in both of hH(3)HR and hH(4)HR are involved in H(3)/ H(4) subtype selectivity. In addition, we find that M378(6.55) in hH(3)HR provides additional hydrophobic interactions different from hH(4)HR (the corresponding amino acid of T323(6.55) in hH(4)HR) to provide additional subtype bias. From these studies, we developed a pharmacophore model based on our predictions for known hH(3)HR selective antagonists in clinical study [ABT-239 1, GSK-189,254 2, PF-3654746 3, and BF2.649 (tiprolisant) 4] that suggests critical selectivity directing elements are: the basic proton interacting with D114(3.32), the spacer, the aromatic ring substituted with the hydrophilic or lipophilic groups interacting with lipophilic pockets in transmembranes (TMs) 3-5-6 and the aliphatic ring located in TMs 2-3-7. These 3D structures for all four HRs should help guide the rational design of novel drugs for the subtype selective antagonists and agonists with reduced side effects.     

Molecular dynamics and thermodynamics of protein-RNA interactions: mutation of a conserved aromatic residue modifies stacking interactions and structural adaptation in the U1A-stem loop 2 RNA complex.

Molecular dynamics (MD) simulations and free energy component analysis have been performed to evaluate the molecular origins of the 5.5 kcal/mol destabilization of the complex formed between the N-terminal RNP domain of U1A and stem loop 2 of U1 snRNA upon mutation of a conserved aromatic residue, Phe56, to Ala. MD simulations, including counterions and water, have been carried out on the wild type and Phe56Ala peptide-stem loop 2 RNA complexes, the free wild type and Phe56Ala peptides, and the free stem loop 2 RNA. The MD structure of the Phe56Ala-stem loop 2 complex is similar to that of the wild type complex except the stacking interaction between Phe56 and A6 of stem loop 2 is absent and loop 3 of the peptide is more dynamic. However, the MD simulations predict large changes in the structure and dynamics of helix C and increased dynamic range of loop 3 for the free Phe56Ala peptide compared to the wild type peptide. Since helix C and loop 3 are highly variable regions of RNP domains, this indicates that a significant contribution to the reduced affinity of the Phe56Ala peptide for RNA results from cooperation between highly conserved and highly variable regions of the RNP domain of U1A. Surprisingly, these structural effects, which are manifested as cooperative free energy changes, occur in the free peptide, rather than in the complex, and are revealed only by study of both the initial and final states of the complexation process. Free energy component analysis correctly accounts for the destabilization of the Phe56Ala-stem loop 2 complex, and indicates that approximately 80% of the destabilization is due to the loss of the stacking interaction and approximately 20% is due to differences in U1A adaptation.     

Conformational transition and energy landscape of ErbB4 activated by neuregulin1β: one microsecond molecular dynamics simulations

ErbB4, a receptor tyrosine kinase of the ErbB family, plays crucial roles in cell growth and differentiation, especially in the development of the heart and nervous system. Ligand binding to its extracellular region could modulate the activation process. To understand the mechanism of ErbB4 activation induced by ligand binding, we performed one microsecond molecular dynamics (MD) simulations on the ErbB4 extracellular region (ECR) with and without its endogenous ligand neuregulin1β (NRG1β). The conformational transition of the ECR-ErbB4/NRG1β complex from a tethered inactive conformation to an extended active-like form has been observed, while such large and function-related conformational change has not been seen in the simulation on the ECR-ErbB4, suggesting that ligand binding is indeed the active inducing force for the conformational transition and further dimerization. On the basis of MD simulations and principal component analysis, we constructed a rough energy landscape for the conformational transition of ECR-ErbB4/NRG1β complex, suggesting that the conformational change from the inactive state to active-like state involves a stable conformation. The energy barrier for the tether opening was estimated as ~2.7 kcal/mol, which is very close to the experimental value (1-2 kcal/mol) reported for ErbB1. On the basis of the simulation results, an atomic mechanism for the ligand-induced activation of ErbB4 was postulated. The present MD simulations provide a new insight into the conformational changes underlying the activation of ErbB4.     

In silico models for the human alpha4beta2 nicotinic acetylcholine receptor

The neuronal alpha4beta2 nicotinic acetylcholine receptor (nAChR) is one of the most widely expressed nAChR subtypes in the brain. Its subunits have high sequence identity (54 and 46% for alpha4 and beta2, respectively) with alpha and beta subunits in Torpedo nAChR. Using the known structure of the Torpedo nAChR as a template, the closed-channel structure of the alpha4beta2 nAChR was constructed through homology modeling. Normal-mode analysis was performed on this closed structure and the resulting lowest frequency mode was applied to it for a "twist-to-open" motion, which increased the minimum pore radius from 2.7 to 3.4 A and generated an open-channel model. Nicotine could bind to the predicted agonist binding sites in the open-channel model but not in the closed one. Both models were subsequently equilibrated in a ternary lipid mixture via extensive molecular dynamics (MD) simulations. Over the course of 11 ns MD simulations, the open channel remained open with filled water, but the closed channel showed a much lower water density at its hydrophobic gate comprised of residues alpha4-V259 and alpha4-L263 and their homologous residues in the beta2 subunits. Brownian dynamics simulations of Na+ permeation through the open channel demonstrated a current-voltage relationship that was consistent with experimental data on the conducting state of alpha4beta2 nAChR. Besides establishment of the well-equilibrated closed- and open-channel alpha4beta2 structural models, the MD simulations on these models provided valuable insights into critical factors that potentially modulate channel gating. Rotation and tilting of TM2 helices led to changes in orientations of pore-lining residue side chains. Without concerted movement, the reorientation of one or two hydrophobic side chains could be enough for channel opening. The closed- and open-channel structures exhibited distinct patterns of electrostatic interactions at the interface of extracellular and transmembrane domains that might regulate the signal propagation of agonist binding to channel opening. A potential prominent role of the beta2 subunit in channel gating was also elucidated in the study.     

芬太尼类化合物与阿片μ受体相互作用的分子对接与分子动力学模拟

采用分子对接和分子动力学(MD)模拟方法研究了芬太尼类化合物与阿片μ受体的相互作用机制.先用AutoDock4.0程序将芬太尼类化合物对接到同源模建的阿片μ受体结构中,再用GROMACS程序包在水溶液体系中分别对12个芬太尼激动剂和阿片μ受体蛋白复合物进行了MD模拟研究,优化对接复合物的结构,最后利用MM-PBSA方法,在APBS程序中计算芬太尼类衍生物与阿片μ受体的结合自由能,计算出的受体配合物结合常数(Ki)与其实验值吻合较好,并预测了化合物的活性排序.结果表明,复合物蛋白结构与空载受体蛋白结构有较大差异,特别是胞内区IL2、IL3和跨膜区段TM4骨架构象变化较大,不同的化合物对受体结构影响也有差异,活性较好的化合物会增加蛋白特定区域结构的柔性.芬太尼类化合物可能是通过和受体结合后诱导阿片μ受体构象转变为活性构象,引起一系列的信号传导激活G蛋白,从而引发生理效应.     

Protein engineering of nitrile hydratase activity of papain: molecular dynamics study of a mutant and wild-type enzyme

The mechanism of hydrolysis of the nitrile (N-acetyl-phenylalanyl-2-amino-propionitrile, I) catalyzed by Gln19Glu mutant of papain has been studied by nanosecond molecular dynamics (MD) simulations. MD simulations of the complex of mutant enzyme with I and of mutant enzyme covalently attached to both neutral (II) and protonated (III) thioimidate intermediates were performed. An MD simulation with the wild-type enzyme.I complex was undertaken as a reference. The ion pair between protonated His159 and thiolate of Cys25 is coplanar, and the hydrogen bonding interaction S(-)(25).HD1-ND1(159) is observed throughout MD simulation of the mutant enzyme.I complex. Such a sustained hydrogen bond is absent in nitrile-bound wild-type papain due to the flexibility of the imidazole ring of His159. The nature of the residue at position 19 plays a critical role in the hydrolysis of the covalent thioimidate intermediate. When position 19 represents Glu, the imidazolium ion of His159-ND1(+).Cys25-S(-) ion pair is distant, on average, from the nitrile nitrogen of substrate I. Near attack conformers (NACs) have been identified in which His159-ImH(+) is positioned to initiate a general acid-catalyzed addition of Cys-S(-) to nitrile. Though Glu19-CO(2)H is distant from nitrile nitrogen in the mutant.I structure, MD simulations of the mutant.II covalent adduct finds Glu19-CO(2)H hydrogen bonded to the thioimide nitrogen of II. This hydrogen bonded species is much less stable than the hydrogen bonded Glu19-CO(2)(-) with mutant-bound protonated thioimidate (III). This observation supports Glu19-CO(2)H general acid catalysis of the formation of mutant.III. This is the commitment step in the Gln19Glu mutant catalysis of nitrile hydrolysis.     

Conserved Tyr223(5.58) plays different roles in the activation and G-protein interaction of rhodopsin

Rhodopsin, a seven transmembrane helix (TM) receptor, binds its ligand 11-cis-retinal via a protonated Schiff base. Coupling to the G-protein transducin (G(t)) occurs after light-induced cis/trans-retinal isomerization, which leads through photoproducts into a sequence of metarhodopsin (Meta) states: Meta I ? Meta IIa ? Meta IIb ? Meta IIbH(+). The structural changes behind this three-step activation scheme are mediated by microswitch domains consisting of conserved amino acids. Here we focus on Tyr223(5.58) as part of the Y(5.58)X(7)K(R)(5.66) motif. Mutation to Ala, Phe, or Glu results in specific impairments of G(t)-activation measured by intrinsic G(t) fluorescence. UV-vis/FTIR spectroscopy of rhodopsin and its complex with a C-terminal G(t)α peptide allows the assignment of these deficiencies to specific steps in the activation path. Effects of mutation occur already in Meta I but do not directly influence deprotonation of the Schiff base during formation of Meta IIa. Absence of the whole phenol ring (Y223A) allows the activating motion of TM6 in Meta IIb but impairs the coupling to G(t). When only the hydroxyl group is lacking (Y223F), Meta IIb does not accumulate, but the activity toward G(t) remains substantial. From the FTIR features of Meta IIbH(+) we conclude that proton uptake to Glu134(3.49) is mandatory for Tyr223(5.58) to engage in the interaction with the key player Arg135(3.50) predicted by X-ray analysis. This polar interaction is partially recovered in Y223E, explaining its relatively high activity. Only the phenol side chain of tyrosine provides all characteristics for accumulation of the active state and G-protein activation.     

基于分子模拟方法预测PIP2在双孔钾通道TREK-1上结合位点的研究

磷脂酰肌醇4,5-二磷酸酯(PIP₂)是一类分布在质膜内层的信号磷脂分子, 对钾、 钠和氯等离子通道和转运蛋白等多种跨膜蛋白具有调节作用. TREK-1是一类重要的背景钾通道, 受温度、 机械拉伸及胞内pH等多种因素调节, PIP₂在特定浓度范围内可激活TREK-1通道, 在内面向外膜片钳记录TREK-1通道电流中使用PIP₂抗结剂(如多聚赖氨酸)可导致TREK-1通道关闭. 利用分子对接和全原子分子动力学模拟探索了PIP₂与双孔钾通道TREK-1的相互作用. 分子对接计算结果表明, PIP₂在TREK-1通道上有两个可能的结合位点. 进一步的分子动力学模拟和均力势(PMF)计算结果表明, 其中位于螺旋M4和螺旋M1的位点可能是PIP₂激活TREK-1的优先结合位点. 模拟展示了PIP₂与TREK-1结合的可能构象. PIP₂的肌醇头部磷酸根与位于M1和M4上的碱性残基K45, K304和R311形成稳定盐桥; M1螺旋上的一系列疏水残基对稳定PIP₂的脂肪长链具有关键作用.     

Molecular dynamics simulations of helix-forming, amine-functionalized m-poly(phenyleneethynylene)s

We present here the results of all-atom and united-atom molecular dynamics (MD) simulations that were used to examine the folding behavior of an amine-functionalized m-poly(phenyleneethynylene) (m-PPE) oligomer in aqueous environment. The parallelized GROMACS MD simulation code and OPLS force field were used for multiple MD simulations of m-PPE oligomers containing 24 phenyl rings in extended, coiled and helix conformations separately in water to determine the minimum energy conformation of the oligomer in aqueous solvent and what interactions are most important in determining this structure. Simulation results showed that the helix is the preferred minimum energy conformation of a single oligomer in water and that Lennard-Jones interactions are the dominant forces for the stabilization of the helix. In addition, these solvophobic interactions are strong enough to maintain the helix conformation at temperatures up to 523 K.     

Conserved patterns in backbone torsional changes allow for single base flipping from duplex DNA with minimal distortion of the double helix

Base flipping is a structural mechanism common to many DNA processing and repair enzymes. Changes in the local backbone torsions that occur during base flipping and the effect of environment on their behavior are of particular interest in understanding different base flipping mechanisms. In the present study, structures sampled during umbrella sampling molecular dynamics (MD) simulations of base flipping in aqueous and protein-bound environments, carried out with two different MD simulation strategies, are analyzed to find the most significant phosphodiester backbone distortions in the vicinity of the flipping base. Torsional sampling on the 5' side of the flipping base during flipping through the major groove shows similarities to the torsional sampling on the 3' side during flipping through the minor groove and vice versa. In differing environments, this behavior varies only marginally. These compensating torsional changes in the DNA backbone on 5' and 3' sides of the flipping base limit overall distortion of the DNA double helix during single base flipping. Rotameric intermediate states observed during base flipping are identified and postulated to be metastable states implicated in both large-scale structural changes and functional effects of chemical modifications in DNA.     

力场/溶剂水模型对EGFRvⅢ抗原-MR1(scFv)抗体复合物的MM-PBSA计算精度的影响与实验验证

以表皮生长因子Ⅲ型突变体(EGFRvⅢ)抗原多肽与其抗体(MR1)及其人源化突变体的复合物结构为出发点,采用分子动力学中的6种常用力场及3种常用溶剂水模型,分别对上述抗原-抗体复合物进行100ns的分子动力学模拟与分子力学和连续介质模型计算自由能(MM-PBSA),并在实验上利用等温滴定量热(ITC)仪测定了抗原和抗体相互作用的热力学参数.通过在结构变化、能量变化及野生型与突变体比较等几个方面进行综合分析,给出了最佳的计算模型.对不同力场及水模型计算精度等相关问题进行了探讨.     

The Dynamics of the Neuropeptide Y Receptor Type 1 Investigated by Solid-State NMR and Molecular Dynamics Simulation

We report data on the structural dynamics of the neuropeptide Y (NPY) G-protein-coupled receptor (GPCR) type 1 (Y1R), a typical representative of class A peptide ligand GPCRs, using a combination of solid-state NMR and molecular dynamics (MD) simulation. First, the equilibrium dynamics of Y1R were studied using ¹⁵N-NMR and quantitative determination of ¹H-¹³C order parameters through the measurement of dipolar couplings in separated-local-field NMR experiments. Order parameters reporting the amplitudes of the molecular motions of the C-H bond vectors of Y1R in DMPC membranes are 0.57 for the Cα sites and lower in the side chains (0.37 for the CH₂ and 0.18 for the CH₃ groups). Different NMR excitation schemes identify relatively rigid and also dynamic segments of the molecule. In monounsaturated membranes composed of longer lipid chains, Y1R is more rigid, attributed to a higher hydrophobic thickness of the lipid membrane. The presence of an antagonist or NPY has little influence on the amplitude of motions, whereas the addition of agonist and arrestin led to a pronounced rigidization. To investigate Y1R dynamics with site resolution, we conducted extensive all-atom MD simulations of the apo and antagonist-bound state. In each state, three replicas with a length of 20 μs (with one exception, where the trajectory length was 10 μs) were conducted. In these simulations, order parameters of each residue were determined and showed high values in the transmembrane helices, whereas the loops and termini exhibit much lower order. The extracellular helix segments undergo larger amplitude motions than their intracellular counterparts, whereas the opposite is observed for the loops, Helix 8, and termini. Only minor differences in order were observed between the apo and antagonist-bound state, whereas the time scale of the motions is shorter for the apo state. Although these relatively fast motions occurring with correlation times of ns up to a few µs have no direct relevance for receptor activation, it is believed that they represent the prerequisite for larger conformational transitions in proteins.     

Classical and quantum mechanical/molecular mechanical molecular dynamics simulations of alanine dipeptide in water: comparisons with IR and vibrational circular dichroism spectra

We have implemented the combined quantum mechanical (QM)/molecular mechanical (MM) molecular dynamics (MD) simulations of alanine dipeptide in water along with the polarizable and nonpolarizable classical MD simulations with different models of water. For the QM/MM MD simulation, the alanine dipeptide is treated with the AM1 or PM3 approximations and the fluctuating solute dipole moment is calculated by the Mulliken population analysis. For the classical MD simulations, the solute is treated with the polarizable or nonpolarizable AMBER and polarizable CHARMM force fields and water is treated with the TIP3P, TIP4P, or TIP5P model. It is found that the relative populations of right-handed alpha-helix and extended beta and P(II) conformations in the simulation trajectory strongly depend on the simulation method. For the QM/MM MD simulations, the PM3/MM shows that the P(II) conformation is dominant, whereas the AM1/MM predicts that the dominant conformation is alpha(R). Polarizable CHARMM force field gives almost exclusively P(II) conformation and other force fields predict that both alpha-helical and extended (beta and P(II)) conformations are populated with varying extents. Solvation environment around the dipeptide is investigated by examining the radial distribution functions and numbers and lifetimes of hydrogen bonds. Comparing the simulated IR and vibrational circular dichroism spectra with experimental results, we concluded that the dipeptide adopts the P(II) conformation and PM3/MM, AMBER03 with TIP4P water, and AMBER polarizable force fields are acceptable for structure determination of the dipeptide considered in this paper.     

An efficient method for computing excess free energy of liquid

We present a new theoretical method for efficient calculation of free energy of liquid. This interaction entropy method allows one to compute entropy and free energy of liquid from standard single step MD (molecular dynamics) simulation directly in liquid state without the need to perform MD simulations at many intermediate states as required in thermodynamic integration or free energy perturbation methods. In this new approach, one only needs to evaluate the interaction energy of a single (fixed) liquid molecule with the rest of liquid molecules as a function of time from a standard MD simulation of liquid and the fluctuation of distribution of this interaction energy is then used to calculate the interaction entropy of the liquid. Explicit theoretical derivation of this interaction entropy approach is provided and numerical calculations for the benchmark liquid water system were carried out using three different water models. Numerical analysis of the result was performed and comparison of the computational result with experimental data and other theoretical results were provided. Excellent agreement of calculated free energies with the experimental data using TIP4P model is obtained for liquid water.