Water Exchange from the Buried Binding Sites of Cytochrome P450 Enzymes 1A2, 2D6, and 3A4 Correlates with Conformational Fluctuations.

Human cytochrome P450 enzymes (CYPs) are critical for the metabolism of small-molecule pharmaceuticals (drugs). As such, the prediction of drug metabolism by and drug inhibition of CYP activity is an important component of the drug discovery and design process. Relative to the availability of a wide range of experimental atomic-resolution CYP structures, the development of structure-based CYP activity models has been limited. To better characterize the role of CYP conformational fluctuations in CYP activity, we perform multiple microsecond-scale all-atom explicit-solvent molecular dynamics (MD) simulations on three CYP isoforms, 1A2, 2D6, and 3A4, which together account for the majority of CYP-mediated drug metabolism. The MD simulations employ a variety of positional restraints, ranging from keeping all CYP atoms close to their experimentally determined coordinates to allowing full flexibility. We find that, with full flexibility, large fluctuations in the CYP binding sites correlate with efficient water exchange from these buried binding sites. This is especially true for 1A2, which, when restrained to its crystallographic conformation, is unable to exchange water between the binding site and bulk solvent. These findings imply that, in addition to crystal structures, a representative ensemble of conformational states ought to be included when developing structure-based CYP activity models.

Site Identification by Ligand Competitive Saturation-Biologics Approach for Structure-Based Protein Charge Prediction.

Protein-based therapeutics typically require high concentrations of the active protein, which can lead to protein aggregation and high solution viscosity. Such solution behaviors can limit the stability, bioavailability, and manufacturability of protein-based therapeutics and are directly influenced by the charge of a protein. Protein charge is a system property affected by its environment, including the buffer composition, pH, and temperature. Thus, the charge calculated by summing the charges of each residue in a protein, as is commonly done in computational methods, may significantly differ from the effective charge of the protein as these calculations do not account for contributions from bound ions. Here, we present an extension of the structure-based approach termed site identification by ligand competitive saturation-biologics (SILCS-Biologics) to predict the effective charge of proteins. The SILCS-Biologics approach was applied on a range of protein targets in different salt environments for which membrane-confined electrophoresis-determined charges were previously reported. SILCS-Biologics maps the 3D distribution and predicted occupancy of ions, buffer molecules, and excipient molecules bound to the protein surface in a given salt environment. Using this information, the effective charge of the protein is predicted such that the concentrations of the ions and the presence of excipients or buffers are accounted for. Additionally, SILCS-Biologics also produces 3D structures of the binding sites of ions on the proteins, which enable further analyses such as the characterization of protein surface charge distribution and dipole moments in different environments. Notable is the capability of the method to account for competition between salts, excipients, and buffers on the calculated electrostatic properties in different protein formulations. Our study demonstrates the ability of the SILCS-Biologics approach to predict the effective charge of proteins and its applicability in uncovering protein-ion interactions and their contributions to protein solubility and function.

Atomic-Resolution Experimental Structural Biology and Molecular Dynamics Simulations of Hyaluronan and Its Complexes.

This review summarizes the atomic-resolution structural biology of hyaluronan and its complexes available in the Protein Data Bank, as well as published studies of atomic-resolution explicit-solvent molecular dynamics simulations on these and other hyaluronan and hyaluronan-containing systems. Advances in accurate molecular mechanics force fields, simulation methods and software, and computer hardware have supported a recent flourish in such simulations, such that the simulation publications now outnumber the structural biology publications by an order of magnitude. In addition to supplementing the experimental structural biology with computed dynamic and thermodynamic information, the molecular dynamics studies provide a wealth of atomic-resolution information on hyaluronan-containing systems for which there is no atomic-resolution structural biology either available or possible. Examples of these summarized in this review include hyaluronan pairing with other hyaluronan molecules and glycosaminoglycans, with ions, with proteins and peptides, with lipids, and with drugs and drug-like molecules. Despite limitations imposed by present-day computing resources on system size and simulation timescale, atomic-resolution explicit-solvent molecular dynamics simulations have been able to contribute significant insight into hyaluronan's flexibility and capacity for intra- and intermolecular non-covalent interactions.

Pyranose Ring Puckering Thermodynamics for Glycan Monosaccharides Associated with Vertebrate Proteins.

The conformational properties of carbohydrates can contribute to protein structure directly through covalent conjugation in the cases of glycoproteins and proteoglycans and indirectly in the case of transmembrane proteins embedded in glycolipid-containing bilayers. However, there continue to be significant challenges associated with experimental structural biology of such carbohydrate-containing systems. All-atom explicit-solvent molecular dynamics simulations provide a direct atomic resolution view of biomolecular dynamics and thermodynamics, but the accuracy of the results depends on the quality of the force field parametrization used in the simulations. A key determinant of the conformational properties of carbohydrates is ring puckering. Here, we applied extended system adaptive biasing force (eABF) all-atom explicit-solvent molecular dynamics simulations to characterize the ring puckering thermodynamics of the ten common pyranose monosaccharides found in vertebrate biology (as represented by the CHARMM carbohydrate force field). The results, along with those for idose, demonstrate that the CHARMM force field reliably models ring puckering across this diverse set of molecules, including accurately capturing the subtle balance between 4C1 and 1C4 chair conformations in the cases of iduronate and of idose. This suggests the broad applicability of the force field for accurate modeling of carbohydrate-containing vertebrate biomolecules such as glycoproteins, proteoglycans, and glycolipids.

Constructing 3-Dimensional Atomic-Resolution Models of Nonsulfated Glycosaminoglycans with Arbitrary Lengths Using Conformations from Molecular Dynamics.

Glycosaminoglycans (GAGs) are the linear carbohydrate components of proteoglycans (PGs) and are key mediators in the bioactivity of PGs in animal tissue. GAGs are heterogeneous, conformationally complex, and polydisperse, containing up to 200 monosaccharide units. These complexities make studying GAG conformation a challenge for existing experimental and computational methods. We previously described an algorithm we developed that applies conformational parameters (i.e., all bond lengths, bond angles, and dihedral angles) from molecular dynamics (MD) simulations of nonsulfated chondroitin GAG 20-mers to construct 3-D atomic-resolution models of nonsulfated chondroitin GAGs of arbitrary length. In the current study, we applied our algorithm to other GAGs, including hyaluronan and nonsulfated forms of dermatan, keratan, and heparan and expanded our database of MD-generated GAG conformations. Here, we show that individual glycosidic linkages and monosaccharide rings in 10- and 20-mers of hyaluronan and nonsulfated dermatan, keratan, and heparan behave randomly and independently in MD simulation and, therefore, using a database of MD-generated 20-mer conformations, that our algorithm can construct conformational ensembles of 10- and 20-mers of various GAG types that accurately represent the backbone flexibility seen in MD simulations. Furthermore, our algorithm efficiently constructs conformational ensembles of GAG 200-mers that we would reasonably expect from MD simulations.

Rigidity and flexibility in the tetrasaccharide linker of proteoglycans from atomic-resolution molecular simulation.

Proteoglycans (PGs) are covalent conjugates between protein and carbohydrate (glycosaminoglycans). Certain classes of glycosaminoglycans such as chondroitin sulfate/dermatan sulfate and heparan sulfate utilize a specific tetrasaccharide linker for attachment to the protein component: GlcAß1-3Galß1-3Galß1-4Xylß1-O-Ser. Toward understanding the conformational preferences of this linker, the present work used all-atom explicit-solvent molecular dynamics (MD) simulations combined with Adaptive Biasing Force (ABF) sampling to determine high-resolution, high-precision conformational free energy maps ΔG(φ, ψ) for each glycosidic linkage between constituent disaccharides, including the variant where GlcA is substituted with IdoA. These linkages are characterized by single, predominant (> 97% occupancy), and broad (45° × 60° for ΔG(φ, ψ) < 1 kcal/mol) free-energy minima, while the Xyl-Ser linkage has two such minima similar in free-energy, and additional flexibility from the Ser sidechain dihedral. Conformational analysis of microsecond-scale standard MD on the complete tetrasaccharide-O-Ser conjugate is consistent with ABF data, suggesting (φ, ψ) probabilities are independent of the linker context, and that the tetrasaccharide acts as a relatively rigid unit whereas significant conformational heterogeneity exists with respect to rotation about bonds connecting Xyl to Ser. © 2017 Wiley Periodicals, Inc.

Folding of Fibroblast Growth Factor 1 Is Critical for Its Nonclassical Release.

Fibroblast growth factor 1 (FGF1), a ubiquitously expressed pro-angiogenic protein that is involved in tissue repair, carcinogenesis, and maintenance of vasculature stability, is released from the cells via a stress-dependent nonclassical secretory pathway. FGF1 secretion is a result of transmembrane translocation of this protein. It correlates with the ability of FGF1 to permeabilize membranes composed of acidic phospholipids. Like several other nonclassically exported proteins, FGF1 exhibits ß-barrel folding. To assess the role of folding of FGF1 in its secretion, we applied targeted mutagenesis in combination with a complex of biophysical methods and molecular dynamics studies, followed by artificial membrane permeabilization and stress-induced release experiments. It has been demonstrated that a mutation of proline 135 located in the C-terminus of FGF1 results in (i) partial unfolding of FGF1, (ii) a decrease in FGF1's ability to permeabilize bilayers composed of phosphatidylserine, and (iii) drastic inhibition of stress-induced FGF1 export. Thus, folding of FGF1 is critical for its nonclassical secretion.

Terminal sialic acids on CD44 N-glycans can block hyaluronan binding by forming competing intramolecular contacts with arginine sidechains.

Specific sugar residues and their linkages form the basis of molecular recognition for interactions of glycoproteins with other biomolecules. Seemingly small changes, like the addition of a single monosaccharide in the covalently attached glycan component of glycoproteins, can greatly affect these interactions. For instance, the sialic acid capping of glycans affects protein-ligand binding involved in cell-cell and cell-matrix interactions. CD44 is a single-pass transmembrane glycoprotein whose binding with its carbohydrate ligand hyaluronan (HA), an extracellular matrix component, mediates processes such as leukocyte homing, cell adhesion, and tumor metastasis. This binding is highly regulated by glycosylation of the N-terminal extracellular hyaluronan-binding domain (HABD); specifically, sialic acid capped N-glycans of HABD inhibit ligand binding. However, the molecular mechanism behind this sialic acid mediated regulation has remained unknown. Two of the five N-glycosyation sites of HABD have been previously identified as having the greatest inhibitory effect on HA binding, but only if the glycans contain terminal sialic acid residues. These two sites, Asn25 and Asn120, were chosen for in silico glycosylation in this study. Here, from extensive standard molecular dynamics simulations and biased simulations, we propose a molecular mechanism for this behavior based on spontaneously-formed charge-paired hydrogen bonding interactions between the negatively-charged sialic acid residues and positively-charged Arg sidechains known to be critically important for binding to HA, which itself is negatively charged. Such intramolecular hydrogen bonds would preclude associations critical to hyaluronan binding. This observation suggests how CD44 and related glycoprotein binding is regulated by sialylation as cellular environments fluctuate.

Supramolecular assemblies of a conjugate of nucleobase, amino acids, and saccharide act as agonists for proliferation of embryonic stem cells and development of zygotes.

The synthetic challenges in glycobiology and glycochemistry hamper the development of glycobiomaterials for biomedicine. Here we report the use of molecular self-assembly to sidestep the laborious synthesis of complex glycans for promoting the proliferation of murine embryonic stem (mES) cells. Our study shows that the supramolecular assemblies of a small molecule conjugate of nucleobase, amino acids, and saccharide, as a de novo glycoconjugate, promote the proliferation of mES cells and the development of zygotes into blastocysts of mouse. Molecular engineering confirms that each motif (i.e., adenine, Arg-Gly-Asp (RGD) domain, and glucosamine) is indispensable for the observed activity of the conjugate. As the first example of using assemblies of the molecular conjugates of multiple fundamental biological building blocks to control cell behaviors, this work illustrates an unprecedented approach to use supramolecular assemblies as multifunctional mimics of glycoconjugates.

Effect of Lipid Bilayer Anchoring on the Conformational Properties of the Cytochrome P450 2D6 Binding Site.

Human cytochrome P450 (CYP) proteins metabolize 75% of small-molecule pharmaceuticals, which makes structure-based modeling of CYP metabolism and inhibition, bolstered by the current availability of X-ray crystal structures of CYP globular catalytic domains, an attractive prospect. Accounting for this broad metabolic capacity is a combination of the existence of multiple different CYP proteins and the capacity of a single CYP protein to metabolize multiple different small molecules. It is thought that structural plasticity and flexibility contribute to this latter property; therefore, incorporating diverse conformational states of a particular CYP is likely an important consideration in structure-based CYP metabolism and inhibition modeling. While all-atom explicit-solvent molecular dynamics simulations can be used to generate conformational ensembles under biologically relevant conditions, existing CYP crystal structures are of the globular domain only, whereas human CYPs contain N-terminal transmembrane and linker peptides that anchor the globular catalytic domain to the endoplasmic reticulum. To determine whether this can cause significant differences in the sampled binding site conformations, microsecond scale all-atom explicit-solvent molecular dynamics simulations of the CYP2D6 globular domain in an aqueous environment were compared with those of the full-length protein anchored in a POPC lipid bilayer. While bilayer-anchoring damped some structural fluctuations in the globular domain relative to the aqueous simulations, none of the affected residues included binding site pocket residues. Furthermore, clustering of molecular dynamics snapshots based on either pairwise binding site pocket RMSD or volume differences demonstrated a lack of separation of snapshots from the two simulation conditions into different clusters. These results suggest the substantially simpler and computationally cheaper aqueous simulation approach can be used to generate a relevant conformational ensemble of the CYP2D6 binding site for structure-based metabolism and inhibition modeling.

CD44 receptor unfolding enhances binding by freeing basic amino acids to contact carbohydrate ligand.

The extracellular carbohydrate-binding domain of the Type I transmembrane receptor CD44 is known to undergo affinity switching, where change in conformation leads to enhanced binding of its carbohydrate ligand hyaluronan. Separate x-ray crystallographic and NMR experiments have led to competing explanations, with the former supporting minor conformational changes at the binding site and the latter a major order-to-disorder unfolding transition distant from the binding site. Here, all-atom explicit-solvent molecular dynamics studies employing adaptive biasing force sampling revealed a substantial favorable free-energy change associated with contact formation between the Arg(41) side chain and hyaluronan at the binding site, independent of whether the distant site was ordered or disordered. Analogous computational experiments on Arg(41)Ala mutants showed loss of this favorable free-energy change, consistent with existing experimental data. More provocatively, the simulation data revealed the molecular mechanism by which the order-to-disorder transition enhances hyaluronan binding: in the disordered state, a number of basic residues gain sufficient conformational freedom-lacking in the ordered state-to spontaneously form side-chain contacts with hyaluronan. Mutation of these residues to Ala had been known to decrease binding affinity, but there had previously been no structural explanation, given their lack of proximity to the carbohydrate-binding site in existing structures of the complex.

Coarse grain lipid-protein molecular interactions and diffusion with MsbA flippase.

Coarse-grained (CG) modeling has proven effective for simulating lipid bilayer dynamics on scales of biological interest. Modeling the dynamics of flexible membrane proteins within the bilayer, on the other hand, poses a considerable challenge due to the complexity of the folding or conformational landscape. In the present work, the multiscale coarse-graining method is applied to atomistic peptide-lipid "soup" simulations to develop a general set of CG protein-lipid interaction potentials. The reduced model was constructed to be compatible with recent solvent-free CG models developed for protein-protein folding and lipid-lipid model bilayer interactions. The utility of the force field was demonstrated by molecular dynamics simulation of the MsbA ABC transporter in a mixed DOPC/DOPE bilayer. An elastic network was parameterized to restrain the MsbA dimer in its open, closed and hydrolysis intermediate conformations and its impact on domain flexibility was examined. Conformational stability enabled long-time dynamics simulation of MsbA freely diffusing in a 25 nm membrane patch. Three-dimensional density analysis revealed that a shell of weakly bound "annular lipids" solvate the membrane accessible surface of MsbA and its internal substrate-binding chamber. The annular lipid binding modes, along with local perturbations in head group structure, are a function of the orientation of grooves formed between transmembrane helices and may influence the alternating access mechanism of substrate entry and translocation.

Balancing target flexibility and target denaturation in computational fragment-based inhibitor discovery.

Accounting for target flexibility and selecting "hot spots" most likely to be able to bind an inhibitor continue to be challenges in the field of structure-based drug design, especially in the case of protein-protein interactions. Computational fragment-based approaches using molecular dynamics (MD) simulations are a promising emerging technology having the potential to address both of these challenges. However, the optimal MD conditions permitting sufficient target flexibility while also avoiding fragment-induced target denaturation remain ambiguous. Using one such technology (Site Identification by Ligand Competitive Saturation, SILCS), conditions were identified to either prevent denaturation or identify and exclude trajectories in which subtle but important denaturation was occurring. The target system used was the well-characterized protein cytokine IL-2, which is involved in a protein-protein interface and, in its unliganded crystallographic form, lacks surface pockets that can serve as small-molecule binding sites. Nonetheless, small-molecule inhibitors have previously been discovered that bind to two "cryptic" binding sites that emerge only in the presence of ligand binding, highlighting the important role of IL-2 flexibility. Using the above conditions, SILCS with hydrophobic fragments was able to identify both sites based on favorable fragment binding while avoiding IL-2 denaturation. An important additional finding was that acetonitrile, a water-miscible fragment, fails to identify either site yet can induce target denaturation, highlighting the importance of fragment choice.

Structural mechanism associated with domain opening in gain-of-function mutations in SHP2 phosphatase.

The SHP2 phosphatase plays a central role in a number of signaling pathways were it dephosphorylates various substrate proteins. Regulation of SHP2 activity is, in part, achieved by an intramolecular interaction between the PTP domain of the protein, which contains the catalytic site, and the N-SH2 domain leading to a "closed" protein conformation and autoinhibition. Accordingly, "opening" of the N-SH2 and PTP domains is required for the protein to become active. Binding of phosphopeptides to the N-SH2 domain is known to induce the opening event, while a number of gain-of-function (GOF) mutants, implicated in Noonan's Syndrome and childhood leukemias, are thought to facilitate opening. In the present study, a combination of computational and experimental methods are used to investigate the structural mechanism of opening of SHP2 and the impact of three GOF mutants, D61G, E76K, and N308D, on the opening mechanism. Calculated free energies of opening indicate that opening must be facilitated by effector molecules, possibly the protein substrates themselves, as the calculated free energies preclude spontaneous opening. Simulations of both wild type (WT) SHP2 and GOF mutants in the closed state indicate GOF activity to involve increased solvent exposure of selected residues, most notably Arg362, which in turn may enhance interactions of SHP2 with its substrate proteins and thereby aid opening. In addition, GOF mutations cause structural changes in the phosphopeptide-binding region of the N-SH2 domain leading to conformations that mimic the bound state. Such conformational changes are suggested to enhance binding of phosphopeptides and/or decrease interactions between the PTP and N-SH2 domains thereby facilitating opening. Experimental assays of the impact of effector molecules on SHP2 phosphatase activity against both small molecule and peptide substrates support the hypothesized mechanism of GOF mutant action. The present calculations also suggest a role for the C-SH2 domain of SHP2 in stabilizing the overall conformation of the protein in the open state, thereby aiding conformational switching between the open active and closed inactive states.

Reproducing crystal binding modes of ligand functional groups using Site-Identification by Ligand Competitive Saturation (SILCS) simulations.

The applicability of a computational method, Site Identification by Ligand Competitive Saturation (SILCS), to identify regions on a protein surface with which different types of functional groups on low-molecular weight inhibitors interact is demonstrated. The method involves molecular dynamics (MD) simulations of a protein in an aqueous solution of chemically diverse small molecules from which probability distributions of fragments types, termed FragMaps, are obtained. In the present application, SILCS simulations are performed with an aqueous solution of 1 M benzene and propane to map the affinity pattern of the protein for aromatic and aliphatic functional groups. In addition, water hydrogen and oxygen atoms serve as probes for hydrogen-bond donor and acceptor affinity, respectively. The method is tested using a set of 7 proteins for which crystal structures of complexes with several high affinity inhibitors are known. Good agreement is obtained between FragMaps and the positions of chemically similar functional groups in inhibitors as observed in the X-ray crystallographic structures. Quantitative capabilities of the SILCS approach are demonstrated by converting FragMaps to free energies, termed Grid Free Energies (GFE), and showing correlation between the GFE values and experimental binding affinities. For proteins for which ligand decoy sets are available, GFE values are shown to typically score the crystal conformation and conformations similar to it more favorable than decoys. Additionally, SILCS is tested for its ability to capture the subtle differences in ligand affinity across homologous proteins, information which may be of utility toward specificity-guided drug design. Taken together, our results show that SILCS can recapitulate the known location of functional groups of bound inhibitors for a number of proteins, suggesting that the method may be of utility for rational drug design.

Sulfation and Calcium Favor Compact Conformations of Chondroitin in Aqueous Solutions.

The effects of sulfation and calcium cations (Ca2+) on the atomic-resolution conformational properties of chondroitin carbohydrate polymers in aqueous solutions are not well studied owing to experimental challenges. Here, we compare all-atom explicit-solvent molecular dynamics simulations results for pairs of O-type (nonsulfated) and A-type (GlcNAc 4-O-sulfated) chondroitin 20-mers in 140 mM NaCl with and without Ca2+ and find that both sulfation and Ca2+ favor more compact polymer conformations. We also show that subtle differences in force-field parametrization can have dramatic effects on Ca2+ binding to chondroitin carboxylate and sulfate functional groups and thereby determine Ca2+-mediated intra- and interstrand association. In addition to providing an atomic-resolution picture of the interaction of Ca2+ with sulfated and nonsulfated chondroitin polymers, the molecular dynamics data emphasize the importance of careful force-field parametrization to balance ion-water and ion-chondroitin interactions and suggest additional parametrization efforts to tune interactions involving sulfate.

Computational fragment-based binding site identification by ligand competitive saturation.

Fragment-based drug discovery using NMR and x-ray crystallographic methods has proven utility but also non-trivial time, materials, and labor costs. Current computational fragment-based approaches circumvent these issues but suffer from limited representations of protein flexibility and solvation effects, leading to difficulties with rigorous ranking of fragment affinities. To overcome these limitations we describe an explicit solvent all-atom molecular dynamics methodology (SILCS: Site Identification by Ligand Competitive Saturation) that uses small aliphatic and aromatic molecules plus water molecules to map the affinity pattern of a protein for hydrophobic groups, aromatic groups, hydrogen bond donors, and hydrogen bond acceptors. By simultaneously incorporating ligands representative of all these functionalities, the method is an in silico free energy-based competition assay that generates three-dimensional probability maps of fragment binding (FragMaps) indicating favorable fragment:protein interactions. Applied to the two-fold symmetric oncoprotein BCL-6, the SILCS method yields two-fold symmetric FragMaps that recapitulate the crystallographic binding modes of the SMRT and BCOR peptides. These FragMaps account both for important sequence and structure differences in the C-terminal halves of the two peptides and also the high mobility of the BCL-6 His116 sidechain in the peptide-binding groove. Such SILCS FragMaps can be used to qualitatively inform the design of small-molecule inhibitors or as scoring grids for high-throughput in silico docking that incorporate both an atomic-level description of solvation and protein flexibility.

Efficient Construction of Atomic-Resolution Models of Non-Sulfated Chondroitin Glycosaminoglycan Using Molecular Dynamics Data.

Glycosaminoglycans (GAGs) are linear, structurally diverse, conformationally complex carbohydrate polymers that may contain up to 200 monosaccharides. These characteristics present a challenge for studying GAG conformational thermodynamics at atomic resolution using existing experimental methods. Molecular dynamics (MD) simulations can overcome this challenge but are only feasible for short GAG polymers. To address this problem, we developed an algorithm that applies all conformational parameters contributing to GAG backbone flexibility (i.e., bond lengths, bond angles, and dihedral angles) from unbiased all-atom explicit-solvent MD simulations of short GAG polymers to rapidly construct models of GAGs of arbitrary length. The algorithm was used to generate non-sulfated chondroitin 10- and 20-mer ensembles which were compared to MD-generated ensembles for internal validation. End-to-end distance distributions in constructed and MD-generated ensembles have minimal differences, suggesting that our algorithm produces conformational ensembles that mimic the backbone flexibility seen in simulation. Non-sulfated chondroitin 100- and 200-mer ensembles were constructed within a day, demonstrating the efficiency of the algorithm and reduction in time and computational cost compared to simulation.

CHARMM additive all-atom force field for aldopentofuranoses, methyl-aldopentofuranosides, and fructofuranose.

An additive all-atom empirical force field for aldopentofuranoses, methyl-aldopentofuranosides (Me-aldopentofuranosides), and fructofuranose carbohydrates, compatible with existing CHARMM carbohydrate parameters, is presented. Building on existing parameters transferred from cyclic ethers and hexopyranoses, parameters were further developed using target data for complete furanose carbohydrates as well as O-methyl tetrahydrofuran. The bond and angle equilibrium parameters were adjusted to reproduce target geometries from a survey of furanose crystal structures, and dihedral parameters were fit to over 1700 quantum mechanical (QM) MP2/cc-pVTZ//MP2/6-31G(d) conformational energies. The conformational energies were for a variety of complete furanose monosaccharides and included two-dimensional ring pucker energy surfaces. Bonded parameter optimization led to the correct description of the ring pucker for a large set of furanose compounds, while furanose-water interaction energies and distances reproduced QM HF/6-31G(d) results for a number of furanose monosaccharides, thereby validating the nonbonded parameters. Crystal lattice unit cell parameters and volumes, aqueous-phase densities, and aqueous NMR ring pucker and exocyclic data were used to validate the parameters in condensed-phase environments. Conformational sampling analysis of the ring pucker and exocyclic group showed excellent agreement with experimental NMR data, demonstrating that the conformational energetics in aqueous solution are accurately described by the optimized force field. Overall, the parameters reproduce available experimental data well and are anticipated to be of utility in future computational studies of carbohydrates, including in the context of proteins, nucleic acids, and/or lipids when combined with existing CHARMM biomolecular force fields.

Additive empirical force field for hexopyranose monosaccharides.

We present an all-atom additive empirical force field for the hexopyranose monosaccharide form of glucose and its diastereomers allose, altrose, galactose, gulose, idose, mannose, and talose. The model is developed to be consistent with the CHARMM all-atom biomolecular force fields, and the same parameters are used for all diastereomers, including both the alpha- and beta-anomers of each monosaccharide. The force field is developed in a hierarchical manner and reproduces the gas-phase and condensed-phase properties of small-molecule model compounds corresponding to fragments of pyranose monosaccharides. The resultant parameters are transferred to the full pyranose monosaccharides, and additional parameter development is done to achieve a complete hexopyranose monosaccharide force field. Parametrization target data include vibrational frequencies, crystal geometries, solute-water interaction energies, molecular volumes, heats of vaporization, and conformational energies, including those for over 1800 monosaccharide conformations at the MP2/cc-pVTZ//MP2/6-31G(d) level of theory. Although not targeted during parametrization, free energies of aqueous solvation for the model compounds compare favorably with experimental values. Also well-reproduced are monosaccharide crystal unit cell dimensions and ring pucker, densities of concentrated aqueous glucose systems, and the thermodynamic and dynamic properties of the exocyclic torsion in dilute aqueous systems. The new parameter set expands the CHARMM additive force field to allow for simulation of heterogeneous systems that include hexopyranose monosaccharides in addition to proteins, nucleic acids, and lipids.