On achieving experimental accuracy from molecular dynamics simulations of flexible molecules: aqueous glycerol

Structure and dynamics in aqueous mixtures of glycerol: insights from molecular dynamics simulations

Aqueous glycerol mixtures are investigated over the whole concentration range of glycerol xGLY = 0.1-0.9 via molecular dynamics (MD) simulations at ambient pressure and temperature. Two glycerol force fields are used: an all-atom (AA) and a united-atom (UA) model. Structural changes upon different mixing ratios are discussed through the site-site radial distribution functions (RDFs), coordination numbers and cluster analysis. As both species are hydrogen bonded, they form an almost perfect H-bonded network, with no observed clusters. There are, however, noticeable changes in the RDFs. Glycerol correlations grow stronger with increasing glycerol content, as do water correlations. There is significant transformation in dynamics as well, as evidenced by the self-diffusion coefficients, the velocity autocorrelation functions and the rotational autocorrelation functions. Diffusion of both species slows down with increasing glycerol content. Rotational relaxation is also altered depending on the mixture composition and there is a slow-down at the lower end of glycerol content.

Free Energy Differences from Molecular Simulations: Exact Confidence Intervals from Transition Counts

Here, we demonstrate a method to estimate the uncertainty (confidence intervals and standard errors) of free energy differences calculated by molecular simulations. The widths of confidence intervals and standard errors can be calculated solely from temperature and the number of transitions between states. Uncertainty (95% confidence interval) lower than ±1 kcal/mol can be achieved by a simulation with four forward and four reverse transitions. For a two-state Markovian system, the confidence interval is exact, regardless the number of transitions.

1H-NMR Karplus Analysis of Molecular Conformations of Glycerol under Different Solvent Conditions: A Consistent Rotational Isomerism in the Backbone Governed by Glycerol/Water Interactions

Glycerol is a symmetrical, small biomolecule with high flexibility in molecular conformations. Using a 1H-NMR spectroscopic Karplus analysis in our way, we analyzed a rotational isomerism in the glycero backbone which generates three kinds of staggered conformers, namely gt (gauche-trans), gg (gauche-gauche), and tg (trans-gauche), at each of sn-1,2 and sn-2,3 positions. The Karplus analysis has disclosed that the three rotamers are consistently equilibrated in water keeping the relation of 'gt:gg:tg = 50:30:20 (%)' at a wide range of concentrations (5 mM~540 mM). The observed relation means that glycerol in water favors those symmetric conformers placing 1,2,3-triol groups in a gauche/gauche geometry. We have found also that the rotational isomerism is remarkably changed when the solvent is replaced with DMSO-d6 or dimethylformamide (DMF-d7). In these solvents, glycerol gives a relation of 'gt:gg:tg = 40:30:30 (%)', which means that a remarkable shift occurs in the equilibrium between gt and tg conformers. By this shift, glycerol turns to also take non-symmetric conformers orienting one of the two vicinal diols in an antiperiplanar geometry.

Mechanisms of Interaction of Small Hydroxylated Cryosolvents with Dehydrated Model Cell Membranes: Stabilization vs Destruction

The mechanism by which cryosolvents such as alcohols modify and penetrate cell membranes as a function of their concentration and hydration state remains poorly understood. We conducted molecular dynamics simulations of 1,2-dioleoyl-sn-glycero-3-phosphocholine bilayers in the presence of aqueous solutions of four common penetrating hydroxylated cryosolvents (methanol, ethylene glycol, propylene glycol, and glycerol) at varying concentration ranges and across three different hydration states. All cryosolvents were found to preferentially replace water at the bilayer interface, and a reduction in hydration state correlates with a higher proportion of cryosolvent at the interface for relative concentrations. Minor differences in chemical structure had a profound effect on cryosolvent-membrane interactions, as the lone methyl groups of methanol and propylene glycol enhanced their membrane localization and penetration, but with increasing concentrations acted to destabilize the membrane structure in a process heightened at higher hydration states. By contrast, ethylene glycol and glycerol promoted and retained membrane structural integrity by forming hydrogen-bonded lipid bridges via distally located hydroxyl groups. Glycerol exhibited the highest capacity to cross-link lipids at relative concentrations, as well as promoted a bilayer structure consistent with a fully hydrated bilayer in the absence of cryosolvent for all hydration states investigated.

Hydrogen Bonding and the Structural Properties of Glycerol-Water Mixtures with a Microwave Field: a Molecular Dynamics Study

In a microwave field, the dielectric properties, molecular structures, and hydrogen bonding dynamics of glycerol in its mixtures with water were determined by the molecular dynamics simulation method. The dipole-dipole correlation of glycerol is linked to the field intensity of microwaves. The results show that as the field intensity is increased, even glycerol in the second coordination shell can become correlated with each other. The structures of up to 35 glycerol molecules are observed. More than that, it was observed that lifetimes of glycerol-glycerol hydrogen bonds were prolonged, while the average hydrogen bond number was also increased. Besides, the structures in a strong microwave field mimic the weak C-H⋯O hydrogen bonds seen in high-glycerol concentration mixtures, yet the concentration is lower. These results indicate that with the assistance of the microwave field, glycerol molecules become concentrated and are more likely to establish stable interactions with others. As a consequence, the spherical clusters composed by glycerol molecules in our nanosheet synthesis experiment are easier to form.

Predicting the Structures of Glycans, Glycoproteins, and Their Complexes

Complex carbohydrates are ubiquitous in nature, and together with proteins and nucleic acids they comprise the building blocks of life. But unlike proteins and nucleic acids, carbohydrates form nonlinear polymers, and they are not characterized by robust secondary or tertiary structures but rather by distributions of well-defined conformational states. Their molecular flexibility means that oligosaccharides are often refractory to crystallization, and nuclear magnetic resonance (NMR) spectroscopy augmented by molecular dynamics (MD) simulation is the leading method for their characterization in solution. The biological importance of carbohydrate-protein interactions, in organismal development as well as in disease, places urgency on the creation of innovative experimental and theoretical methods that can predict the specificity of such interactions and quantify their strengths. Additionally, the emerging realization that protein glycosylation impacts protein function and immunogenicity places the ability to define the mechanisms by which glycosylation impacts these features at the forefront of carbohydrate modeling. This review will discuss the relevant theoretical approaches to studying the three-dimensional structures of this fascinating class of molecules and interactions, with reference to the relevant experimental data and techniques that are key for validation of the theoretical predictions.

Extension and validation of the GLYCAM force field parameters for modeling glycosaminoglycans

Glycosaminoglycans (GAGs) are an important class of carbohydrates that serve critical roles in blood clotting, tissue repair, cell migration and adhesion, and lubrication. The variable sulfation pattern and iduronate ring conformations in GAGs influence their polymeric structure and nature of interaction. This study characterizes several heparin-like GAG disaccharides and tetrasaccharides using NMR and molecular dynamics simulations to assist in the development of parameters for GAGs within the GLYCAM06 force field. The force field additions include parameters and charges for a transferable sulfate group for O- and N-sulfation, neutral (COOH) forms of iduronic and glucuronic acid, and Δ4,5-unsaturated uronate (ΔUA) residues. ΔUA residues frequently arise from the enzymatic digestion of heparin and heparin sulfate. Simulations of disaccharides containing ΔUA reveal that the presence of sulfation on this residue alters the relative populations of 1H2 and 2H1 ring conformations. Simulations of heparin tetrasaccharides containing N-sulfation in place of N-acetylation on glucosamine residues influence the ring conformations of adjacent iduronate residues.

Glass polymorphism in glycerol-water mixtures: I. A computer simulation study

We perform out-of-equilibrium molecular dynamics (MD) simulations of water-glycerol mixtures in the glass state. Specifically, we study the transformations between low-density (LDA) and high-density amorphous (HDA) forms of these mixtures induced by compression/decompression at constant temperature. Our MD simulations reproduce qualitatively the density changes observed in experiments. Specifically, the LDA-HDA transformation becomes (i) smoother and (ii) the hysteresis in a compression/decompression cycle decreases as T and/or glycerol content increase. This is surprising given the fast compression/decompression rates (relative to experiments) accessible in MD simulations. We study mixtures with glycerol molar concentration χ(g) = 0-13% and find that, for the present mixture models and rates, the LDA-HDA transformation is detectable up to χ(g) ≈ 5%. As the concentration increases, the density of the starting glass (i.e., LDA at approximately χ(g) ≤ 5%) rapidly increases while, instead, the density of HDA remains practically constant. Accordingly, the LDA state and hence glass polymorphism become inaccessible for glassy mixtures with approximately χ(g) > 5%. We present an analysis of the molecular-level changes underlying the LDA-HDA transformation. As observed in pure glassy water, during the LDA-to-HDA transformation, water molecules within the mixture approach each other, moving from the second to the first hydration shell and filling the first interstitial shell of water molecules. Interestingly, similar changes also occur around glycerol OH groups. It follows that glycerol OH groups contribute to the density increase during the LDA-HDA transformation. An analysis of the hydrogen bond (HB)-network of the mixtures shows that the LDA-HDA transformation is accompanied by minor changes in the number of HBs of water and glycerol. Instead, large changes in glycerol and water coordination numbers occur. We also perform a detailed analysis of the effects that the glycerol force field (FF) has on our results. By comparing MD simulations using two different glycerol models, we find that glycerol conformations indeed depend on the FF employed. Yet, the thermodynamic and microscopic mechanisms accompanying the LDA-HDA transformation and hence, our main results, do not. This work is accompanied by an experimental report where we study the glass polymorphism in glycerol-water mixtures prepared by isobaric cooling at 1 bar.

Computational NMR coupling constants: shifting and scaling factors for evaluating 1JCH

Optimized shifting and/or scaling factors for calculating one-bond carbon-hydrogen spin-spin coupling constants have been determined for 35 combinations of representative functionals (PBE, B3LYP, B3P86, B97-2 and M06-L) and basis sets (TZVP, HIII-su3, EPR-III, aug-cc-pVTZ-J, ccJ-pVDZ, ccJ-pVTZ, ccJ-pVQZ, pcJ-2 and pcJ-3) using 68 organic molecular systems with 88 (1)JCH couplings including different types of hybridized carbon atoms. Density functional theory assessment for the determination of (1)JCH coupling constants is examined, comparing the computed and experimental values. The use of shifting constants for obtaining the calculated coupling improves substantially the results, and most models become qualitatively similar. Thus, for the whole set of couplings and for all approaches excluding those using the M06 functional, the root-mean-square deviations lie between 4.7 and 16.4 Hz and are reduced to 4-6.5 Hz when shifting constants are considered. Alternatively, when a specific rovibrational contribution of 5 Hz is subtracted from the experimental values, good results are obtained with PBE, B3P86 and B97-2 functionals in combination with HIII-su3, aug-cc-pVTZ-J and pcJ-2 basis sets.

Polarizable empirical force field for acyclic polyalcohols based on the classical Drude oscillator

A polarizable empirical force field for acyclic polyalcohols based on the classical Drude oscillator is presented. The model is optimized with an emphasis on the transferability of the developed parameters among molecules of different sizes in this series and on the condensed-phase properties validated against experimental data. The importance of the explicit treatment of electronic polarizability in empirical force fields is demonstrated in the cases of this series of molecules with vicinal hydroxyl groups that can form cooperative intra- and intermolecular hydrogen bonds. Compared to the CHARMM additive force field, improved treatment of the electrostatic interactions avoids overestimation of the gas-phase dipole moments resulting in significant improvement in the treatment of the conformational energies and leads to the correct balance of intra- and intermolecular hydrogen bonding of glycerol as evidenced by calculated heat of vaporization being in excellent agreement with experiment. Computed condensed phase data, including crystal lattice parameters and volumes and densities of aqueous solutions are in better agreement with experimental data as compared to the corresponding additive model. Such improvements are anticipated to significantly improve the treatment of polymers in general, including biological macromolecules.

Motional timescale predictions by molecular dynamics simulations: case study using proline and hydroxyproline sidechain dynamics

We propose a new approach for force field optimizations which aims at reproducing dynamics characteristics using biomolecular MD simulations, in addition to improved prediction of motionally averaged structural properties available from experiment. As the source of experimental data for dynamics fittings, we use ¹³C NMR spin-lattice relaxation times T₁ of backbone and sidechain carbons, which allow to determine correlation times of both overall molecular and intramolecular motions. For structural fittings, we use motionally averaged experimental values of NMR J couplings. The proline residue and its derivative 4-hydroxyproline with relatively simple cyclic structure and sidechain dynamics were chosen for the assessment of the new approach in this work. Initially, grid search and simplexed MD simulations identified large number of parameter sets which fit equally well experimental J couplings. Using the Arrhenius-type relationship between the force constant and the correlation time, the available MD data for a series of parameter sets were analyzed to predict the value of the force constant that best reproduces experimental timescale of the sidechain dynamics. Verification of the new force-field (termed as AMBER99SB-ILDNP) against NMR J couplings and correlation times showed consistent and significant improvements compared to the original force field in reproducing both structural and dynamics properties. The results suggest that matching experimental timescales of motions together with motionally averaged characteristics is the valid approach for force field parameter optimization. Such a comprehensive approach is not restricted to cyclic residues and can be extended to other amino acid residues, as well as to the backbone. Proteins 2014; 82:195–215. © 2013 Wiley Periodicals, Inc.

Contribution of high-energy conformations to NMR chemical shifts, a DFT-BOMD study

This paper highlights the relevance of including the high-energy conformational states sampled by Born-Oppenheimer molecular dynamics (BOMD) in the calculation of time-averaged NMR chemical shifts. Our case study is the very flexible glycerol molecule that undergoes interconversion between conformers in a nonrandom way. Along the sequence of structures from one backbone conformer to another, transition states have been identified. The three (13)C NMR chemical shifts of the molecule were estimated by averaging their calculated values over a large set of BOMD snapshots. The simulation time needed to obtain a good agreement with the two signals present in the experimental spectrum is shown to be dependent on the atomic orbital basis set used for the dynamics, with a necessary longer trajectory for the most extended basis sets. The large structural deformations with respect to the optimized conformer geometries that occur along the dynamics are related to a kinetically driven conformer distribution. Calculated conformer type populations are in good agreement with experimental gas phase microwave results.

Conformational flexibility of N-glycans in solution studied by REMD simulations

Protein-glycan recognition regulates a wide range of biological and pathogenic processes. Conformational diversity of glycans in solution is apparently incompatible with specific binding to their receptor proteins. One possibility is that among the different conformational states of a glycan, only one conformer is utilized for specific binding to a protein. However, the labile nature of glycans makes characterizing their conformational states a challenging issue. All-atom molecular dynamics (MD) simulations provide the atomic details of glycan structures in solution, but fairly extensive sampling is required for simulating the transitions between rotameric states. This difficulty limits application of conventional MD simulations to small fragments like di- and tri-saccharides. Replica-exchange molecular dynamics (REMD) simulation, with extensive sampling of structures in solution, provides a valuable way to identify a family of glycan conformers. This article reviews recent REMD simulations of glycans carried out by us or other research groups and provides new insights into the conformational equilibria of N-glycans and their alteration by chemical modification. We also emphasize the importance of statistical averaging over the multiple conformers of glycans for comparing simulation results with experimental observables. The results support the concept of "conformer selection" in protein-glycan recognition.

Molecular dynamics simulation study of glycerol-water liquid mixtures

To study the effects of water on conformational dynamics of polyalcohols, Molecular Dynamics simulations of glycerol-water liquid mixtures have been carried out at different concentrations: 42.9 and 60.0 wt % of glycerol, respectively. On the basis of the analysis of backbone conformer distributions, it is found that the surrounding water molecules have a large impact on the populations of the glycerol conformers. While the local structure of water in the liquid mixture is surprisingly close to that in pure liquid water, the behavior of glycerols can be divided into three different categories where roughly 25% of them occur in a structure similar to that in pure liquid of glycerol, ca. 25% of them exist as monomers, solvated by water, and the remaining 50% of glycerols in the mixture form H-bonded strings as remains of the glycerol H-bond network. The typical glycerol H-bond network still exists even at the lower concentration of 40 wt % of glycerol. The microheterogeneity of water-glycerol mixtures is analyzed using time-averaged distributions of the sizes of the water aggregates. At 40 wt % of glycerol, the cluster sizes from 3 to 10 water molecules are observed. The increase of glycerol content causes a depletion of clusters leading to smaller 3-5 molecule clusters domination. Translational diffusion coefficients have been calculated to study the dynamical behavior of both glycerol and water molecules. Rotational-reorientational motion is studied both in overall and in selected substructures on the basis of time correlation functions. Characteristic time scales for different motional modes are deduced on the basis of the calculated correlation times. The general conclusion is that the presence of water increases the overall mobility of glycerol, while glycerol slows the mobility of water.

Experimental verification of force fields for molecular dynamics simulations using Gly-Pro-Gly-Gly

Experimental NMR verification of MD simulations using 12 different force fields (AMBER, CHARMM, GROMOS, and OPLS-AA) and 5 different water models has been undertaken to identify reliable MD protocols for structure and dynamics elucidations of small open chain peptides containing Gly and Pro. A conformationally flexible tetrapeptide Gly-Pro-Gly-Gly was selected for NMR (3)J-coupling, chemical shift, and internuclear distance measurements, followed by their calculations using 2 μs long MD simulations in water. In addition, Ramachandran population maps for Pro-2 and Gly-3 residues of GPGG obtained from MD simulations were used for detailed comparisons with similar maps from the protein data bank (PDB) for large number of Gly and Pro residues in proteins. The MD simulations revealed strong dependence of the populations and geometries of preferred backbone and side chain conformations, as well as the time scales of the peptide torsional transitions on the force field used. On the basis of the analysis of the measured and calculated data, AMBER99SB is identified as the most reliable force field for reproducing NMR measured parameters, which are dependent on the peptide backbone and the Pro side chain geometries and dynamics. Ramachandran maps showing the dependence of conformational populations as a function of backbone ϕ/ψ angles for Pro-2 and Gly-3 residues of GPGG from MD simulations using AMBER99SB, AMBER03, and CHARMM were found to resemble similar maps for Gly and Pro residues from the PDB survey. Three force fields (AMBER99, AMBER99ϕ, and AMBER94) showed the least satisfactory agreement with both the solution NMR and the PDB survey data. The poor performance of these force fields is attributed to their propensity to overstabilize helical peptide backbone conformations at the Pro-2 and Gly-3 residues. On the basis of the similarity of the MD and PDB Ramachandran plots, the following sequence of transitions is suggested for the Gly backbone conformation: α(L) ⇆ β(PR) ⇆ β(S) ⇆ β(P) ⇆ α, where backbone secondary structures α(L) and α are associated with helices and turns, β(P) and β(PR) correspond to the left- and right-handed polyproline II structures and β(S) denotes the fully stretched backbone conformation. Compared to the force field dependence, less significant, but noteworthy, variations in the populations of the peptide backbone conformations were observed. For different solvent models considered, a correlation was noted between the number of torsional transitions in GPGG and the water self-diffusion coefficient on using TIP3P, TIP4P, and TIP5P models. In addition to MD results, we also report DFT derived Karplus relationships for Gly and Pro residues using B972 and B3LYP functionals.

Structure, dynamics, and energetics of lysobisphosphatidic acid (LBPA) isomers

Lysobisphosphatidic acid (LBPA), or bis(monoacylglycerol)phosphate, is a very interesting lipid, that is mainly found in late endosomes. It has several intriguing characteristics, which differ from those of other animal glycerophospholipids, that may be related to its specific functions, particularly in the metabolism of cholesterol. Its phosphodiester group is bonded at the sn-1 (sn-1') positions of the glycerols rather than at sn-3 (sn-3'); the position of the two fatty acid chains is still under debate but, increasingly, arguments favor the sn-2, sn-2' position in the native molecule, whereas isolation procedures or acidic conditions lead to the thermodynamically more stable sn-3, sn-3' structure. Because of these peculiar features, it can be expected that LBPA shape and interactions with membrane lipids and proteins are related to its structure at the molecular level. We applied quantum mechanical methods to study the structures and stabilities of the 2,2' and 3,3' LBPA isomers, using a step-by-step procedure from glycerol to precursors (in vitro syntheses) and to the final isoforms. The structures of the two positional LBPA isomers are substantially different, showing that the binding positions of the fatty acid chains on the glycerol backbone determine the shape of the LBPA molecule and thus, possibly, its functions. The 3,3' LBPA structures obtained are more stable with respect to the 2,2' form, as expected from experiment. If one argues that the in vivo synthesis starts from the present glycerol conformers and considering the most stable bis(glycero)phosphate structures, the 2,2' isoform should be the most probable isomer.

Computational glycoscience: characterizing the spatial and temporal properties of glycans and glycan-protein complexes

Modern computational methods offer the tools to provide insight into the structural and dynamic properties of carbohydrate-protein complexes, beyond that provided by experimental structural biology. Dynamic properties such as the fluctuation of inter-molecular hydrogen bonds, the residency times of bound water molecules, side chain motions and ligand flexibility may be readily determined computationally. When taken with respect to the unliganded states, these calculations can also provide insight into the entropic and enthalpic changes in free energy associated with glycan binding. In addition, virtual ligand screening may be employed to predict the three dimensional (3D) structures of carbohydrate-protein complexes, given 3D structures for the components. In principle, the 3D structure of the protein may itself be derived by modeling, leading to the exciting--albeit high risk--realm of virtual structure prediction. This latter approach is appealing, given the difficulties associated with generating experimental 3D structures for some classes of glycan binding proteins; however, it is also the least robust. An unexpected outcome of the development of algorithms for modeling carbohydrate-protein interactions has been the discovery of errors in reported experimental 3D structures and a heightened awareness of the need for carbohydrate-specific computational tools for assisting in the refinement and curation of carbohydrate-containing crystal structures. Here we present a summary of the basic strategies associated with employing classical force field based modeling approaches to problems in glycoscience, with a focus on identifying typical pitfalls and limitations. This is not an exhaustive review of the current literature, but hopefully will provide a guide for the glycoscientist interested in modeling carbohydrates and carbohydrate-protein complexes, as well as the computational chemist contemplating such tasks.

Molecular simulations of carbohydrates and protein-carbohydrate interactions: motivation, issues and prospects

The characterization of the 3D structure of oligosaccharides, their conjugates and analogs is particularly challenging for traditional experimental methods. Molecular simulation methods provide a basis for interpreting sparse experimental data and for independently predicting conformational and dynamic properties of glycans. Here, we summarize and analyze the issues associated with modeling carbohydrates, with a detailed discussion of four of the most recently developed carbohydrate force fields, reviewed in terms of applicability to natural glycans, carbohydrate-protein complexes and the emerging area of glycomimetic drugs. In addition, we discuss prospectives and new applications of carbohydrate modeling in drug discovery.

The conformational properties of methyl alpha-(2,8)-di/trisialosides and their N-acyl analogues: implications for anti-Neisseria meningitidis B vaccine design

The conformational properties of di- and trisaccharide fragments of the polysialic acid O-antigen capsular polysaccharide (CPS) of Neisseria meningitidis B (NmB) have been investigated by a combination of solution phase NMR spectroscopy and explicit-solvent molecular dynamics (MD) simulations. Simulations employing 100 ns of conventional MD, as well as 160 ns of replica exchange MD (REMD), with the GLYCAM06 force field were shown to be in agreement with experimental NMR scalar J-coupling and NOE values. The presence of conformational families has been determined by monitoring interglycosidic torsion angles, by comparing structural superimpositions, as well as via a Bayesian statistical analysis of the torsional data. Attempts to augment the immunogenicity of NmB CPS often involve chemical modifications of the N-acetyl moiety. Here the effects of these chemical group modifications on the conformational properties of the trisialoside have been probed via REMD simulations of the N-glycolyl, N-propionyl, N-propyl and N-butanoyl analogues. Although there were conformational families unique to each non-native analogue, the chemical modifications resulted in largely equivalent overall conformational phase-spaces compared to the native trisialoside. On the basis of the conformational distributions, these shared conformational properties suggest that a recurrent global conformational epitope may be present in both the native and chemically modified CPS fragments. Explanations are therefore provided for monoclonal antibody cross-reactivity, in terms of recognition of a shared global CPS conformation, as well as for lack of cross-reactivity, in terms of fine structural differences associated with the N-acyl groups, which may be dominant in highly matured antibody responses.