CHARMM: the biomolecular simulation program

Modeling anesthetic binding sites within the glycine alpha one receptor based on prokaryotic ion channel templates: the problem with TM4

Ligand-gated ion channels (LGICs) significantly modulate anesthetic effects. Their exact molecular structure remains unknown. This has led to ambiguity regarding the proper amino acid alignment within their 3D structure and, in turn, the location of any anesthetic binding sites. Current controversies suggest that such a site could be located in either an intra- or intersubunit locale within the transmembrane domain of the protein. Here, we built a model of the glycine alpha one receptor (GlyRa1) based on the open-state structures of two new high-resolution ion channel templates from the prokaryote, Gloebacter violaceus (GLIC). Sequence scoring suggests reasonable homology between GlyRa1 and GLIC. Three of the residues notable for modulating anesthetic action are on transmembrane segments 1-3 (TM1-3): (ILE229, SER 267, and ALA 288). They line an intersubunit interface, in contrast to previous models. However, residues from the fourth transmembrane domain (TM4) that are known to modulate a variety of anesthetic effects are quite distant from this putative anesthetic binding site. While this model can account for a large proportion of the physicochemical data regarding such proteins, it cannot readily account for the alterations on anesthetic effects that are due to mutations within TM4.

Probing the structural determinants for the function of intracellular loop 2 in structurally cognate G-protein-coupled receptors

Intracellular loop 2 (IL2) in G-protein-coupled receptors (GPCRs) is functionally important, e.g., in binding to G-protein and β-arrestin. Differences in secondary structure of IL2 in the crystal structures of the very similar β(1)- and β(2)-adrenergic receptors (β(1)AR and β(2)AR, respectively), i.e., an α-helix and an L-shaped strand, respectively, emphasize the need to understand the structural basis for IL2 functionality. We studied the properties of IL2 in the context of experimental data using a Monte Carlo-based ab initio method. The procedure was validated first by verifying that the IL2 structures in β(1)AR and β(2)AR crystals were correctly reproduced, even after conformational ensemble searches at >1200 K where most secondary structure had been lost. We found that IL2 in β(1)AR and β(2)AR sampled each other's conformation but adopted different energetically preferred conformations, consistent with the crystal structures. The results indicate a persistent contextual preference for the structure of IL2, which was conserved when the IL2 sequences were interchanged between the receptors. We conclude that the protein environment, more than the IL2 sequence, regulates the IL2 structures. We extended the approach to the molecular model of 5-HT(2A)R for which no crystal structure is available and found that IL2 is predominantly helical, similar to IL2 in β(1)AR. Because the P3.57A mutation in IL2 had been shown to decrease β-arrestin binding and internalization, we predicted the effects of the mutation and found that it decreased the propensity of IL2 to form helix, identifying the helical IL2 as a component of the GPCR active form.

Applications of computational science for understanding enzymatic deconstruction of cellulose

Understanding the molecular-level mechanisms that enzymes employ to deconstruct plant cell walls is a fundamental scientific challenge with significant ramifications for renewable fuel production from biomass. In nature, bacteria and fungi use enzyme cocktails that include processive and non-processive cellulases and hemicellulases to convert cellulose and hemicellulose to soluble sugars. Catalyzed by an accelerated biofuels R&D portfolio, there is now a wealth of new structural and experimental insights related to cellulases and the structure of plant cell walls. From this background, computational approaches commonly used in other fields are now poised to offer insights complementary to experiments designed to probe mechanisms of plant cell wall deconstruction. Here we outline the current status of computational approaches for a collection of critical problems in cellulose deconstruction. We discuss path sampling methods to measure rates of elementary steps of enzyme action, coarse-grained modeling for understanding macromolecular, cellulosomal complexes, methods to screen for enzyme improvements, and studies of cellulose at the molecular level. Overall, simulation is a complementary tool to understand carbohydrate-active enzymes and plant cell walls, which will enable industrial processes for the production of advanced, renewable fuels.

Intramolecular hydrogen bonding in articaine can be related to superior bone tissue penetration: a molecular dynamics study

Local anesthetics (LAs) are drugs that cause reversible loss of nociception during surgical procedures. Articaine is a commonly used LA in dentistry that has proven to be exceptionally effective in penetrating bone tissue and induce anesthesia on posterior teeth in maxilla and mandibula. In the present study, our aim was to gain a deeper understanding of the penetration of articaine through biological membranes by studying the interactions of articaine with a phospholipid membrane. Our approach involves Langmuir monolayer experiments combined with molecular dynamics simulations. Membrane permeability of LAs can be modulated by pH due to a titratable amine group with a pKa value close to physiological pH. A change in protonation state is thus known to act as a lipophilicity switch in LAs. Our study shows that articaine has an additional unique lipophilicity switch in its ability to form an intramolecular hydrogen bond. We suggest this intramolecular hydrogen bond as a novel and additional solvent-dependent mechanism for modulation of lipophilicity of articaine which may enhance its diffusion through membranes and connective tissue.

Conformational dynamics in human purine nucleoside phosphorylase with reactants and transition-state analogues

Dynamic motions of human purine nucleoside phosphorylase (hPNP) in complex with transition-state analogues and reactants were studied using 10 ns explicit solvent molecular dynamics simulations. hPNP is a homotrimer that catalyzes the phosphorolysis of purine 6-oxynucleosides. The ternary complex of hPNP includes the binding of a ligand and phosphate to the active site. Molecular dynamics simulations were performed on the ternary complex of six ligands including the picomolar transition-state analogues, Immucillin-H (K(d) = 56 pM), DADMe-Immucillin-H (K(d) = 8.5 pM), DATMe-Immucillin-H (K(d) = 8.6 pM), SerMe-Immucillin-H (K(d) = 5.2 pM), the substrate inosine, and a complex containing only phosphate. Protein-inhibitor complexes of the late transition-state inhibitors, DADMe-Imm-H and DATMe-Imm-H, are inflexible. Despite the structural similarity of SerMe-Imm-H and DATMe-Imm-H, the protein complex of SerMe-Imm-H is flexible, and the inhibitor is highly mobile within the active sites. All inhibitors exhibit an increased number of nonbonding interactions in the active site relative to the substrate inosine. Water density within the catalytic site is lower for DADMe-ImmH, DATMe-Imm-H, and SerMe-Imm-H than that for the substrate inosine. Tight binding of the picomolar inhibitors results from increased interactions within the active site and a reduction in the number of water molecules organized within the catalytic site relative to the substrate inosine.

Perturbing the folding energy landscape of the bacterial immunity protein Im7 by site-specific N-linked glycosylation

N-linked glycosylation modulates protein folding and stability through a variety of mechanisms. As such there is considerable interest in the development of general rules to predict the structural consequences of site-specific glycosylation and to understand how these effects can be exploited in the design and development of modified proteins with advantageous properties. In this study, expressed protein ligation is used to create site-specifically glycosylated variants of the bacterial immunity protein Im7 modified with the chitobiose disaccharide (GlcNAc-GlcNAc). Glycans were introduced at seven solvent exposed sites within the Im7 sequence and the kinetic and thermodynamic consequences of N-linked glycosylation analyzed. The ΔΔG° values for glycan incorporation were found to range from +5.2 to -3.8 kJ·mol(-1). In several cases, glycosylation influences folding by modulating the local conformational preferences of the glycosylated sequence. These locally mediated effects are most prominent in the center of α-helices where glycosylation negatively effects folding and in compact turn motifs between segments of ordered secondary structure where glycosylation promotes folding and enhances the overall stability of the native protein. The studies also provide insight into why glycosylation is commonly identified at the transition between different types of secondary structure and when glycosylation may be used to elaborate protein structure to protect disordered sequences from proteolysis or immune system recognition.

Ring puckering: a metric for evaluating the accuracy of AM1, PM3, PM3CARB-1, and SCC-DFTB carbohydrate QM/MM simulations

The puckered conformations of furanose and pyranose carbohydrate rings are central to analyzing the action of enzymes on carbohydrates. Enzyme reaction mechanisms are generally inaccessible to experiments and so have become the focus of QM(semiempirical)/MM simulations. We show that the complete free energy of puckering is required to evaluate the accuracy of semiempirical methods used to study reactions involving carbohydrates. Interestingly, we find that reducing the free energy space to lower dimensions results in near meaningless minimum energy pathways. We analyze the furanose and pyranose free energy pucker surfaces and volumes using AM1, PM3, PM3CARB-1, and SCC-DFTB. A comparison with DFT optimized structures and a HF free energy surface reveals that SCC-DFTB provides the best semiempirical description of five- and six-membered carbohydrate ring deformation.

Structural determinants for transport across the intestinal bile acid transporter using C-24 bile acid conjugates

The human apical sodium dependent bile acid transporter (hASBT) reabsorbs gram quantities of bile acid daily and is a potential prodrug target to increase oral drug absorption. In the absence of a high resolution hASBT crystal structure, 3D-QSAR modeling may prove beneficial in designing prodrug targets to hASBT. The objective was to derive a conformationally sampled pharmacophore 3D-QSAR (CSP-SAR) model for the uptake of bile acid conjugates by hASBT. A series of bile acid conjugates of glutamyl chenodeoxycholate were evaluated in terms of K(m) and normalized V(max) (normV(max)) using hASBT-MDCK cells. All monoanionic conjugates were potent substrates. Dianions, cations and zwitterions, which bound with a high affinity, were not substrates. CSP-SAR models were derived using structural and physicochemical descriptors, and evaluated via cross validation. The best CSP-SAR model for K(m) included two structural and two physiochemical descriptors, where substrate hydrophobicity enhanced affinity. A best CSP-SAR model for K(m)/normV(max) employed one structural and three physicochemical descriptors, also indicating hydrophobicity enhanced efficiency. Overall, the bile acid C-24 region accommodated a range of substituted anilines, provided a single negative charge was present near C-24. In comparing uptake findings to prior inhibition results, increased hydrophobicity enhanced activity, with dianions and zwitterions hindering activity.

Probing the periplasmic-open state of lactose permease in response to sugar binding and proton translocation

Based on the crystal structure of lactose permease (LacY) open to the cytoplasm, a hybrid molecular simulation approach with self-guided Langevin dynamics is used to describe conformational changes that lead to a periplasmic-open state. This hybrid approach consists of implicit (IM) and explicit (EX) membrane simulations and requires self-guided Langevin dynamics to enhance protein motions during the IM simulations. The pore radius of the lumen increases by 3.5 Å on the periplasmic side and decreases by 2.5 Å on the cytoplasmic side (relative to the crystal structure), suggesting a lumen that is fully open to the periplasm to allow for extracellular sugar transport and closed to the cytoplasm. Based on our simulations, the mechanism that triggers this conformational change to the periplasmic-open state is the protonation of Glu269 and binding of the disaccharide. Then, helix packing is destabilized by breaking of several side chains involved in hydrogen bonding (Asn245, Ser41, Glu374, Lys42, and Gln242). For the periplasmic-open conformations obtained from our simulations, helix-helix distances agree well with experimental measurements using double electron-electron resonance, fluorescence resonance energy transfer, and varying sized cross-linkers. The periplasmic-open conformations are also in compliance with various substrate accessibility/reactivity measurements that indicate an opening of the protein lumen on the periplasmic side on sugar binding. The comparison with these measurements suggests a possible incomplete closure of the cytoplasmic half in our simulations. However, the closure is sufficient to prevent the disaccharide from transporting to the cytoplasm, which is in accordance with the well-established alternating access model. Ser53, Gln60, and Phe354 are determined to be important in sugar transport during the periplasmic-open stage of the sugar transport cycle and the sugar is found to undergo an orientational change in order to escape the protein lumen.

Allosteric regulation of glycogen synthase kinase 3β: a theoretical study

Glycogen synthase kinase 3β (GSK-3β) is a serine-threonine kinase belonging to the CMGC family that plays a key role in many biological processes, such as glucose metabolism, cell cycle regulation, and proliferation. Like most protein kinases, GSK-3β is regulated via multiple pathways and sites. We performed all-atom molecular dynamics simulations on the unphosphorylated and phosphorylated unbound GSK-3β and the phosphorylated GSK-3β bound to a peptide substrate, its product, and a derived inhibitor. We found that GSK-3β autophosphorylation at residue Tyr(216) results in widening of the catalytic groove, thereby facilitating substrate access. In addition, we studied the interactions of the phosphorylated GSK-3β with a substrate and peptide inhibitor located at the active site and observed higher affinity of the inhibitor to the kinase. Furthermore, we detected a potential remote binding site which was previously identified in other kinases. In agreement with experiments we observed that binding of specific peptides at this remote site leads to stabilization of the activation loop located in the active site. We speculate that this stabilization could enhance the catalytic activity of the kinase. We point to this remote site as being structurally conserved and suggest that the allosteric phenomenon observed here may occur in the protein kinase superfamily.

A gain-of-function mutation of Arabidopsis cryptochrome1 promotes flowering

Plants use different classes of photoreceptors to collect information about their light environment. Cryptochromes are blue light photoreceptors that control deetiolation, entrain the circadian clock, and are involved in flowering time control. Here, we describe the cry1-L407F allele of Arabidopsis (Arabidopsis thaliana), which encodes a hypersensitive cryptochrome1 (cry1) protein. Plants carrying the cry1-L407F point mutation have elevated expression of CONSTANS and FLOWERING LOCUS T under short-day conditions, leading to very early flowering. These results demonstrate that not only the well-studied cry2, with an unequivocal role in flowering promotion, but also cry1 can function as an activator of the floral transition. The cry1-L407F mutants are also hypersensitive toward blue, red, and far-red light in hypocotyl growth inhibition. In addition, cry1-L407F seeds are hypersensitive to germination-inducing red light pulses, but the far-red reversibility of this response is not compromised. This demonstrates that the cry1-L407F photoreceptor can increase the sensitivity of phytochrome signaling cascades. Molecular dynamics simulation of wild-type and mutant cry1 proteins indicated that the L407F mutation considerably reduces the structural flexibility of two solvent-exposed regions of the protein, suggesting that the hypersensitivity might result from a reduced entropic penalty of binding events during downstream signal transduction. Other nonmutually exclusive potential reasons for the cry1-L407F gain of function are the location of phenylalanine-407 close to three conserved tryptophans, which could change cry1's photochemical properties, and stabilization of ATP binding, which could extend the lifetime of the signaling state of cry1.

Molecular dynamics simulations of protein dynamics and their relevance to drug discovery

Molecular dynamics simulations have become increasingly useful in studying biological systems of biomedical interest, and not just in the study of model or toy systems. In this article, the methods and principles of all-atom molecular dynamics will be elucidated with several examples provided of their utility to investigators interested on drug discovery.

Models of membrane-bound Alzheimer's Abeta peptide assemblies

Although it is clear that amyloid beta (Aβ) peptides play a pivotal role in the development of Alzheimer's disease, the precise molecular model of action remains unclear. Aβ peptide forms assemble both in aqueous solution and in lipid membranes. It has been proposed that deleterious effects occur when the peptides interact with membranes, possibly by forming Ca(2+) permeant ion channels. In the accompanying manuscript, we propose models in which the C-terminus third of six Aβ42 peptides forms a six-stranded β-barrel in highly toxic soluble oligomers. Here we extend this hypothesis to membrane-bound assemblies. In these Aβ models, the hydrophobic β-barrel of a hexamer may either reside on the surface of the bilayer, or span the bilayer. Transmembrane pores are proposed to form between several hexamers. Once the β-barrels of six hexamers have spanned the bilayer, they may merge to form a more stable 36-stranded β-barrel. We favor models in which parallel β-barrels formed by N-terminus segments comprise the lining of the pores. These types of models explain why the channels are selective for cations and how metal ions, such as Zn(2+) , synthetic peptides that contain histidines, and some small organic cations may block channels or inhibit formation of channels. Our models were developed to be consistent with microscopy studies of Aβ assemblies in membranes, one of which is presented here for the first time.

In silico experiments of single-chain antibody fragment against drugs of abuse

Three sets of in silico experiments have been conducted to elucidate the binding mechanics of two drugs, (+)-methamphetamine (METH) and amphetamine (AMP) to the single-chain variable fragment (scFv) recently engineered from anti-METH monoclonal antibody mAb6H4 (IgG, κlight chain, K(d)=11nM). The first set of in silico experiments are long time equilibration runs of scFv:drug complexes and of drug-free scFv both in the solution. They demonstrate how the solution structures of scFv deviate from its crystallographic form with or without drug molecules bound to it. They lead to the prediction that the Arrhenius activation barrier is nearly zero for transitions from the dissociated state to the bound state. The second set of in silico experiments are nonequilibrium dynamics of pulling the drug molecules out of the binding pocket of scFv and the equilibration runs for drugs to fall back into the binding pocket. They demonstrate that extra water molecules (in addition to the two crystallographic waters) exist inside the binding pocket, underneath the drug molecules. These extra waters must have been evaporated from the binding pockets during the crystallization process of the in vitro experiments of structural determination. The third set of in silico experiments are nonequilibrium steered molecular dynamics simulations to determine the absolute binding free energies of METH and AMP to scFv. The center of mass of a drug molecule (METH or AMP) is steered (pulled) towards (forward) and away from (reverse) the binding site, sampling forward and reverse pulling paths. Mechanic work is measured along the pulling paths. The work measurements are averaged through the Brownian dynamics fluctuation dissipation theorem to produce the free-energy profiles of the scFv:drug complexes as a function of the drug-scFv separation. These experiments lead to the theoretical prediction of absolute binding energies of METH and AMP that are in agreement with the in vitro experimental results.

Beta-barrel models of soluble amyloid beta oligomers and annular protofibrils

Both soluble and membrane-bound prefibrillar assemblies of Abeta (Aβ) peptides have been associated with Alzheimer's disease (AD). The size and nature of these assemblies vary greatly and are affected by many factors. Here, we present models of soluble hexameric assemblies of Aβ42 and suggest how they can lead to larger assemblies and eventually to fibrils. The common element in most of these assemblies is a six-stranded β-barrel formed by the last third of Aβ42, which is composed of hydrophobic residues and glycines. The hydrophobic core β-barrels of the hexameric models are shielded from water by the N-terminus and central segments. These more hydrophilic segments were modeled to have either predominantly β or predominantly α secondary structure. Molecular dynamics simulations were performed to analyze stabilities of the models. The hexameric models were used as starting points from which larger soluble assemblies of 12 and 36 subunits were modeled. These models were developed to be consistent with numerous experimental results.

Absolute binding free energy calculations of sparsomycin analogs to the bacterial ribosome

The interactions of the 50S subunit of bacterial ribosome with antibiotic sparsomycin (SPS) and five analogs (AN) are investigated through the calculation of the standard (absolute) binding free energy and the characterization of conformational dynamics. The standard binding free energies of the complexes are computed using free energy perturbation molecular dynamics (FEP/MD) simulations with explicit solvent. Restraining potentials are applied and then released during the simulation to efficiently sample the changes in translational, orientational, and conformational freedom of the ligand and receptor upon binding. The biasing effects of the restraining potentials are rigorously removed. The loss of conformational freedom of the ligand upon binding is determined by introducing a potential of mean force (PMF) as a function of the root-mean-square deviation (rmsd) of the ligand relative to its conformation in the bound state. To reduce the size of the simulated system, the binding pocket of the ribosome is simulated in the framework of the generalized solvent boundary potential (GSBP). The number of solvent molecules in the buried binding site is treated via grand canonical Monte Carlo (GCMC) during the FEP/MD simulations. The correlation coefficient between the calculated and measured binding free energies is 0.96, and the experimentally observed ranking order for the binding affinities of the six ligands is reproduced. However, while the calculated affinities of the strong binders agree well with the experimental values, those for the weak binders are underestimated.

Molecular simulations on the thermal stabilization of DNA by hyperthermophilic chromatin protein Sac7d, and associated conformational transitions

Sac7d belongs to a family of chromosomal proteins, which are crucial for thermal stabilization of DNA at higher growth temperatures. It is capable of binding DNA nonspecifically, and is responsible for the increase in the melting temperature of DNA in the bound form up to 85 °C. Molecular dynamics (MD) simulations were performed at different temperatures on two protein-DNA complexes of Sac7d. Various structural and energetic parameters were calculated to examine the DNA stability and to investigate the conformational changes in DNA and the protein-DNA interactions. Room temperature simulations indicated very good agreement with the experimental structures. The protein structure is nearly unchanged at both 300 and 360 K, and only up to five base pairs of the DNA are stabilized by Sac7d at 360 K. However, the MD simulations on DNA alone systems show that they lose their helical structures at 360 K further supporting the role of Sac7d in stabilizing the oligomers. At higher temperatures (420 and 480 K), DNA undergoes denaturation in the presence and the absence of the protein. The DNA molecules were found to undergo B- to A-form transitions consistent with experimental studies, and the extent of these transitions are examined in detail. The extent of sampling B- and A-form regions was found to show temperature and sequence dependence. Multiple MD simulations yielded similar results validating the proposed model. Interaction energy calculations corresponding to protein-DNA binding indicates major contribution due to DNA backbone, explaining the nonspecific interactions of Sac7d.

In vivo patterns of resistance to the HIV attachment inhibitor BMS-488043

Attachment inhibitors (AI) are a novel class of HIV-1 antivirals, with little information available on clinical resistance. BMS-488043 is an orally bioavailable AI that binds to gp120 of HIV-1 and abrogates its binding to CD4(+) lymphocytes. A clinical proof-of-concept study of the AI BMS-488043, administered as monotherapy for 8 days, demonstrated significant viral load reductions. In order to examine the effects of AI monotherapy on HIV-1 sensitivity, phenotypic sensitivity assessment of baseline and postdosing (day 8) samples was performed. These analyses revealed that four subjects had emergent phenotypic resistance (a 50% effective concentration [EC(50)] >10-fold greater than the baseline value) and four had high baseline EC(50)s (>200 nM). Population sequencing and sequence determination of cloned envelope genes uncovered five gp120 mutations at four loci (V68A, L116I, S375I/N, and M426L) associated with BMS-488043 resistance. Substitution at the 375 locus, located near the CD4 binding pocket, was the most common (maintained in 5/8 subjects at day 8). The five substitutions were evaluated for their effects on AI sensitivity through reverse genetics in functional envelopes, confirming their role in decreasing sensitivity to the drug. Additional analyses revealed that these substitutions did not alter sensitivity to other HIV-1 entry inhibitors. Thus, our studies demonstrate that although the majority of the subjects' viruses maintained sensitivity to BMS-488043, substitutions can be selected that decrease HIV-1 susceptibility to the AI. Most importantly, the substitutions described here are not associated with resistance to other approved antiretrovirals, and therefore, attachment inhibitors could complement the current arsenal of anti-HIV agents.

Performance of ZDOCK and ZRANK in CAPRI rounds 13-19

We report the performance of the ZDOCK and ZRANK algorithms in CAPRI rounds 13-19 and introduce a novel measure atom contact frequency (ACF). To compute ACF, we identify the residues that most often make contact with the binding partner in the complete set of ZDOCK predictions for each target. We used ACF to predict the interface of the proteins, which, in combination with the biological data available in the literature, is a valuable addition to our docking pipeline. Furthermore, we incorporated a straightforward and efficient clustering algorithm with two purposes: (1) to determine clusters of similar docking poses (corresponding to energy funnels) and (2) to remove redundancies from the final set of predictions. With these new developments, we achieved at least one acceptable prediction for targets 29 and 36, at least one medium-quality prediction for targets 41 and 42, and at least one high-quality prediction for targets 37 and 40; thus, we succeeded for six out of a total of 12 targets.

A reoptimized GROMOS force field for hexopyranose-based carbohydrates accounting for the relative free energies of ring conformers, anomers, epimers, hydroxymethyl rotamers, and glycosidic linkage conformers

This article presents a reoptimization of the GROMOS 53A6 force field for hexopyranose-based carbohydrates (nearly equivalent to 45A4 for pure carbohydrate systems) into a new version 56A(CARBO) (nearly equivalent to 53A6 for non-carbohydrate systems). This reoptimization was found necessary to repair a number of shortcomings of the 53A6 (45A4) parameter set and to extend the scope of the force field to properties that had not been included previously into the parameterization procedure. The new 56A(CARBO) force field is characterized by: (i) the formulation of systematic build-up rules for the automatic generation of force-field topologies over a large class of compounds including (but not restricted to) unfunctionalized polyhexopyranoses with arbritrary connectivities; (ii) the systematic use of enhanced sampling methods for inclusion of experimental thermodynamic data concerning slow or unphysical processes into the parameterization procedure; and (iii) an extensive validation against available experimental data in solution and, to a limited extent, theoretical (quantum-mechanical) data in the gas phase. At present, the 56A(CARBO) force field is restricted to compounds of the elements C, O, and H presenting single bonds only, no oxygen functions other than alcohol, ether, hemiacetal, or acetal, and no cyclic segments other than six-membered rings (separated by at least one intermediate atom). After calibration, this force field is shown to reproduce well the relative free energies of ring conformers, anomers, epimers, hydroxymethyl rotamers, and glycosidic linkage conformers. As a result, the 56A(CARBO) force field should be suitable for: (i) the characterization of the dynamics of pyranose ring conformational transitions (in simulations on the microsecond timescale); (ii) the investigation of systems where alternative ring conformations become significantly populated; (iii) the investigation of anomerization or epimerization in terms of free-energy differences; and (iv) the design of simulation approaches accelerating the anomerization process along an unphysical pathway.

On the role of anionic lipids in charged protein interactions with membranes

We investigate the role of anionic lipids in the binding to, and subsequent movement of charged protein groups in lipid membranes, to help understand the role of membrane composition in all membrane-active protein sequences. We demonstrate a small effect of phosphatidylglycerol (PG) lipids on the ability of an arginine (Arg) side chain to bind to, and cross a lipid membrane, despite possessing a neutralizing charge. We observe similar membrane deformations in lipid bilayers composed of phosphatidylcholine (PC) and PC/PG mixtures, with comparable numbers of water and lipid head groups pulled into the bilayer hydrocarbon core, and prohibitively large ~20 kcal/mol barriers for Arg transfer across each bilayer, dropping by just 2-3 kcal/mol due to the binding of PG lipids. We explore the causes of this small effect of introducing PG lipids and offer an explanation in terms of the limited membrane interaction for the choline groups of PC lipids bound to the translocating ion. Our calculations reveal a surprising lack of preference for Arg binding to PG lipids themselves, but a small increase in interfacial binding affinity for lipid bilayers containing PG lipids. These results help to explain the nature of competitive lipid binding to charged protein sequences, with implications for a wide range of membrane binding domains and cell perturbing peptides.

Electron microscopy and persistence length analysis of semi-rigid smooth muscle tropomyosin strands

The structural mechanics of tropomyosin are essential determinants of its affinity and positioning on F-actin. Thus, tissue-specific differences among tropomyosin isoforms may influence both access of actin-binding proteins along the actin filaments and the cooperativity of actin-myosin interactions. Here, 40 nm long smooth and striated muscle tropomyosin molecules were rotary-shadowed and compared by means of electron microscopy. Electron microscopy shows that striated muscle tropomyosin primarily consists of single molecules or paired molecules linked end-to-end. In contrast, smooth muscle tropomyosin is more a mixture of varying-length chains of end-to-end polymers. Both isoforms are characterized by gradually bending molecular contours that lack obvious signs of kinking. The flexural stiffness of the tropomyosins was quantified and evaluated. The persistence lengths along the shaft of rotary-shadowed smooth and striated muscle tropomyosin molecules are equivalent to each other (approximately 100 nm) and to values obtained from molecular-dynamics simulations of the tropomyosins; however, the persistence length surrounding the end-to-end linkage is almost twofold higher for smooth compared to cardiac muscle tropomyosin. The tendency of smooth muscle tropomyosin to form semi-rigid polymers with continuous and undampened rigidity may compensate for the lack of troponin-based structural support in smooth muscles and ensure positional fidelity on smooth muscle thin filaments.

A QM/MM study of the phosphoryl transfer to the Kemptide substrate catalyzed by protein kinase A. The effect of the phosphorylation state of the protein on the mechanism

We present here a theoretical study of the phosphoryl transfer catalytic mechanism of protein kinase A, which is the best known member of the large protein kinase family. We have built different theoretical models of the complete PKA-Mg(2)-ATP-substrate system to explore the two most accepted reaction pathways, using for the first time in a reaction mechanism theoretical study, the heptapeptide substrate Kemptide, which is relevant for its high efficiency and small size. The effect of the protein configuration, as modeled by two different X-ray structures with different phosphorylation states and degrees of flexibility, has been analyzed. The results indicate that the environmental conditions can influence the availability of the pathways and thus the choice of the mechanism to be followed. In addition, the roles of the two active site conserved residues, Asp166 and Lys168, have been analyzed for each reaction mechanism.

Two mechanisms of ion selectivity in protein binding sites

A theoretical framework is presented to clarify the molecular determinants of ion selectivity in protein binding sites. The relative free energy of a bound ion is expressed in terms of the main coordinating ligands coupled to an effective potential of mean force representing the influence of the rest of the protein. The latter is separated into two main contributions. The first includes all the forces keeping the ion and the coordinating ligands confined to a microscopic subvolume but does not prevent the ligands from adapting to a smaller or larger ion. The second regroups all the remaining forces that control the precise geometry of the coordinating ligands best adapted to a given ion. The theoretical framework makes it possible to delineate two important limiting cases. In the limit where the geometric forces are dominant (rigid binding site), ion selectivity is controlled by the ion-ligand interactions within the matching cavity size according to the familiar "snug-fit" mechanism of host-guest chemistry. In the limit where the geometric forces are negligible, the ion and ligands behave as a "confined microdroplet" that is free to fluctuate and adapt to ions of different sizes. In this case, ion selectivity is set by the interplay between ion-ligand and ligand-ligand interactions and is controlled by the number and the chemical type of ion-coordinating ligands. The framework is illustrated by considering the ion-selective binding sites in the KcsA channel and the LeuT transporter.

Structure of apolipoprotein A-I N terminus on nascent high density lipoproteins

Apolipoprotein A-I (apoA-I) is the major protein component of high density lipoproteins (HDL) and a critical element of cholesterol metabolism. To better elucidate the role of the apoA-I structure-function in cholesterol metabolism, the conformation of the apoA-I N terminus (residues 6-98) on nascent HDL was examined by electron paramagnetic resonance (EPR) spectroscopic analysis. A series of 93 apoA-I variants bearing single nitroxide spin label at positions 6-98 was reconstituted onto 9.6-nm HDL particles (rHDL). These particles were subjected to EPR spectral analysis, measuring regional flexibility and side chain solvent accessibility. Secondary structure was elucidated from side-chain mobility and molecular accessibility, wherein two major α-helical domains were localized to residues 6-34 and 50-98. We identified an unstructured segment (residues 35-39) and a β-strand (residues 40-49) between the two helices. Residues 14, 19, 34, 37, 41, and 58 were examined by EPR on 7.8, 8.4, and 9.6 nm rHDL to assess the effect of particle size on the N-terminal structure. Residues 14, 19, and 58 showed no significant rHDL size-dependent spectral or accessibility differences, whereas residues 34, 37, and 41 displayed moderate spectral changes along with substantial rHDL size-dependent differences in molecular accessibility. We have elucidated the secondary structure of the N-terminal domain of apoA-I on 9.6 nm rHDL (residues 6-98) and identified residues in this region that are affected by particle size. We conclude that the inter-helical segment (residues 35-49) plays a role in the adaptation of apoA-I to the particle size of HDL.

T-cell epitope prediction and immune complex simulation using molecular dynamics: state of the art and persisting challenges

Atomistic Molecular Dynamics provides powerful and flexible tools for the prediction and analysis of molecular and macromolecular systems. Specifically, it provides a means by which we can measure theoretically that which cannot be measured experimentally: the dynamic time-evolution of complex systems comprising atoms and molecules. It is particularly suitable for the simulation and analysis of the otherwise inaccessible details of MHC-peptide interaction and, on a larger scale, the simulation of the immune synapse. Progress has been relatively tentative yet the emergence of truly high-performance computing and the development of coarse-grained simulation now offers us the hope of accurately predicting thermodynamic parameters and of simulating not merely a handful of proteins but larger, longer simulations comprising thousands of protein molecules and the cellular scale structures they form. We exemplify this within the context of immunoinformatics.

Revisiting hydrophobic mismatch with free energy simulation studies of transmembrane helix tilt and rotation

Protein-lipid interaction and bilayer regulation of membrane protein functions are largely controlled by the hydrophobic match between the transmembrane (TM) domain of membrane proteins and the surrounding lipid bilayer. To systematically characterize responses of a TM helix and lipid adaptations to a hydrophobic mismatch, we have performed a total of 5.8-mus umbrella sampling simulations and calculated the potentials of mean force (PMFs) as a function of TM helix tilt angle under various mismatch conditions. Single-pass TM peptides called WALPn (n = 16, 19, 23, and 27) were used in two lipid bilayers with different hydrophobic thicknesses to consider hydrophobic mismatch caused by either the TM length or the bilayer thickness. In addition, different flanking residues, such as alanine, lysine, and arginine, instead of tryptophan in WALP23 were used to examine their influence. The PMFs, their decomposition, and trajectory analysis demonstrate that 1), tilting of a single-pass TM helix is the major response to a hydrophobic mismatch; 2), TM helix tilting up to approximately 10 degrees is inherent due to the intrinsic entropic contribution arising from helix precession around the membrane normal even under a negative mismatch; 3), the favorable helix-lipid interaction provides additional driving forces for TM helix tilting under a positive mismatch; 4), the minimum-PMF tilt angle is generally located where there is the hydrophobic match and little lipid perturbation; 5), TM helix rotation is dependent on the specific helix-lipid interaction; and 6), anchoring residues at the hydrophilic/hydrophobic interface can be an important determinant of TM helix orientation.

Extending the PRIME model for protein aggregation to all 20 amino acids

We extend PRIME, an intermediate-resolution protein model previously used in simulations of the aggregation of polyalanine and polyglutamine, to the description of the geometry and energetics of peptides containing all 20 amino acid residues. The 20 amino acid side chains are classified into 14 groups according to their hydrophobicity, polarity, size, charge, and potential for side chain hydrogen bonding. The parameters for extended PRIME, called PRIME 20, include hydrogen-bonding energies, side chain interaction range and energy, and excluded volume. The parameters are obtained by applying a perceptron-learning algorithm and a modified stochastic learning algorithm that optimizes the energy gap between 711 known native states from the PDB and decoy structures generated by gapless threading. The number of independent pair interaction parameters is chosen to be small enough to be physically meaningful yet large enough to give reasonably accurate results in discriminating decoys from native structures. The most physically meaningful results are obtained with 19 energy parameters.

Coarse-grained model for simulation of RNA three-dimensional structures

The accurate prediction of an RNA's three-dimensional structure from its "primary structure" will have a tremendous influence on the experimental design and its interpretation and ultimately our understanding of the many functions of RNA. This paper presents a general coarse-grained (CG) potential for modeling RNA 3-D structures. Each nucleotide is represented by five pseudo atoms, two for the backbone (one for the phosphate and another for the sugar) and three for the base to represent base-stacking interactions. The CG potential has been parametrized from statistical analysis of 688 RNA experimental structures. Molecular dynamic simulations of 15 RNA molecules with the length of 12-27 nucleotides have been performed using the CG potential, with performance comparable to that from all-atom simulations. For ~75% of systems tested, simulated annealing led to native-like structures at least once out of multiple repeated runs. Furthermore, with weak distance restraints based on the knowledge of three to five canonical Watson-Crick pairs, all 15 RNAs tested are successfully folded to within 6.5 Å of native structures using the CG potential and simulated annealing. The results reveal that with a limited secondary structure model the current CG potential can reliably predict the 3-D structures for small RNA molecules. We also explored an all-atom force field to construct atomic structures from the CG simulations.

A novel side-chain orientation dependent potential derived from random-walk reference state for protein fold selection and structure prediction

BACKGROUND

An accurate potential function is essential to attack protein folding and structure prediction problems. The key to developing efficient knowledge-based potential functions is to design reference states that can appropriately counteract generic interactions. The reference states of many knowledge-based distance-dependent atomic potential functions were derived from non-interacting particles such as ideal gas, however, which ignored the inherent sequence connectivity and entropic elasticity of proteins.

METHODOLOGY

We developed a new pair-wise distance-dependent, atomic statistical potential function (RW), using an ideal random-walk chain as reference state, which was optimized on CASP models and then benchmarked on nine structural decoy sets. Second, we incorporated a new side-chain orientation-dependent energy term into RW (RWplus) and found that the side-chain packing orientation specificity can further improve the decoy recognition ability of the statistical potential.

SIGNIFICANCE

RW and RWplus demonstrate a significantly better ability than the best performing pair-wise distance-dependent atomic potential functions in both native and near-native model selections. It has higher energy-RMSD and energy-TM-score correlations compared with other potentials of the same type in real-life structure assembly decoys. When benchmarked with a comprehensive list of publicly available potentials, RW and RWplus shows comparable performance to the state-of-the-art scoring functions, including those combining terms from multiple resources. These data demonstrate the usefulness of random-walk chain as reference states which correctly account for sequence connectivity and entropic elasticity of proteins. It shows potential usefulness in structure recognition and protein folding simulations. The RW and RWplus potentials, as well as the newly generated I-TASSER decoys, are freely available in http://zhanglab.ccmb.med.umich.edu/RW.

Enhanced thermostability of methyl parathion hydrolase from Ochrobactrum sp. M231 by rational engineering of a glycine to proline mutation

Protein thermostability can be increased by some glycine to proline mutations in a target protein. However, not all glycine to proline mutations can improve protein thermostability, and this method is suitable only at carefully selected mutation sites that can accommodate structural stabilization. In this study, homology modeling and molecular dynamics simulations were used to select appropriate glycine to proline mutations to improve protein thermostability, and the effect of the selected mutations was proved by the experiments. The structure of methyl parathion hydrolase (MPH) from Ochrobactrum sp. M231 (Ochr-MPH) was constructed by homology modeling, and molecular dynamics simulations were performed on the modeled structure. A profile of the root mean square fluctuations of Ochr-MPH was calculated at the nanosecond timescale, and an eight-amino acid loop region (residues 186-193) was identified as having high conformational fluctuation. The two glycines nearest to this region were selected as mutation targets that might affect protein flexibility in the vicinity. The structures and conformational fluctuations of two single mutants (G194P and G198P) and one double mutant (G194P/G198P) were modeled and analyzed using molecular dynamics simulations. The results predicted that the mutant G194P had the decreased conformational fluctuation in the loop region and might increase the thermostability of Ochr-MPH. The thermostability and kinetic behavior of the wild-type and three mutant enzymes were measured. The results were consistent with the computational predictions, and the mutant G194P was found to have higher thermostability than the wild-type enzyme.

A molecular mechanics model for imatinib and imatinib:kinase binding

Imatinib is an important anticancer drug, which binds specifically to the Abl kinase and blocks its signalling activity. To model imatinib:protein interactions, we have developed a molecular mechanics force field for imatinib and four close analogues, which is consistent with the CHARMM force field for proteins and nucleic acids. Atomic charges and Lennard-Jones parameters were derived from a supermolecule ab initio approach. We considered the ab initio energies and geometries of a probe water molecule interacting with imatinib fragments at 32 different positions. We considered both a neutral and a protonated imatinib. The final RMS deviation between the ab initio and force field energies, averaged over both forms, was 0.2 kcal/mol. The model also reproduces the ab initio geometry and flexibility of imatinib. To apply the force field to imatinib:Abl simulations, it is also necessary to determine the most likely imatinib protonation state when it binds to Abl. This was done using molecular dynamics free energy simulations, where imatinib is reversibly protonated during a series of MD simulations, both in solution and in complex with Abl. The simulations indicate that imatinib binds to Abl in its protonated, positively-charged form. To help test the force field and the protonation prediction, we did MD free energy simulations that compare the Abl binding affinities of two imatinib analogs, obtaining good agreement with experiment. Finally, two new imatinib variants were considered, one of which is predicted to have improved Abl binding. This variant could be of interest as a potential drug.

Assessment of the validity of the double superhelix model for reconstituted high density lipoproteins: a combined computational-experimental approach

For several decades, the standard model for high density lipoprotein (HDL) particles reconstituted from apolipoprotein A-I (apoA-I) and phospholipid (apoA-I/HDL) has been a discoidal particle ∼100 Å in diameter and the thickness of a phospholipid bilayer. Recently, Wu et al. (Wu, Z., Gogonea, V., Lee, X., Wagner, M. A., Li, X. M., Huang, Y., Undurti, A., May, R. P., Haertlein, M., Moulin, M., Gutsche, I., Zaccai, G., Didonato, J. A., and Hazen, S. L. (2009) J. Biol. Chem. 284, 36605-36619) used small angle neutron scattering to develop a new model they termed double superhelix (DSH) apoA-I that is dramatically different from the standard model. Their model possesses an open helical shape that wraps around a prolate ellipsoidal type I hexagonal lyotropic liquid crystalline phase. Here, we used three independent approaches, molecular dynamics, EM tomography, and fluorescence resonance energy transfer spectroscopy (FRET) to assess the validity of the DSH model. (i) By using molecular dynamics, two different approaches, all-atom simulated annealing and coarse-grained simulation, show that initial ellipsoidal DSH particles rapidly collapse to discoidal bilayer structures. These results suggest that, compatible with current knowledge of lipid phase diagrams, apoA-I cannot stabilize hexagonal I phase particles of phospholipid. (ii) By using EM, two different approaches, negative stain and cryo-EM tomography, show that reconstituted apoA-I/HDL particles are discoidal in shape. (iii) By using FRET, reconstituted apoA-I/HDL particles show a 28-34-Å intermolecular separation between terminal domain residues 40 and 240, a distance that is incompatible with the dimensions of the DSH model. Therefore, we suggest that, although novel, the DSH model is energetically unfavorable and not likely to be correct. Rather, we conclude that all evidence supports the likelihood that reconstituted apoA-I/HDL particles, in general, are discoidal in shape.

Free energy surfaces for the interaction of D-glucose with planar aromatic groups in aqueous solution

Multidimensional potentials of mean force for the interactions in aqueous solution of both anomers of D-glucopyranose with two planar aromatic molecules, indole and para-methyl-phenol, have been calculated using molecular dynamics simulations with umbrella sampling and were subsequently used to estimate binding free energies. Indole and para-methyl-phenol serve as models for the side chains of the amino acids tryptophan and tyrosine, respectively. In all cases, a weak affinity between the glucose molecules and the flat aromatic surfaces was found. The global minimum for these interactions was found to be for the case when the pseudoplanar face of β-D-glucopyranose is stacked against the planar surfaces of the aromatic residues. The calculated binding free energies are in good agreement with both experiment and previous simulations. The multidimensional free energy maps suggest a mechanism that could lend kinetic stability to the complexes formed by sugars bound to sugar-binding proteins.

Simulating polarizable molecular ionic liquids with Drude oscillators

The Drude oscillator model is applied to the molecular ionic liquid 1-ethyl-3-methyl-imidazolium triflate. The range of manageable Drude charges is tested. The strength of the polarizability is systematically varied from 0% to 100%. The influence on the structure, single particle dynamics, and collective dielectric properties is investigated. The generalized dielectric constant can be decomposed into a dielectric permittivity, a dielectric conductivity, and an optical dielectric constant ɛ(∞). The major part of the static generalized dielectric constant comes from the collective rotation of the ions, i.e., the dielectric permittivity. The translational contribution from the dielectric conductivity is about 58% of the dielectric permittivity. For the evaluation of the optical dielectric contribution, the computational dielectric theory was adapted to the case of heterogeneous polarizabilities. In case of 100% polarizability, it reaches a value of approximately 2.

Water evaporation and conformational changes from partially solvated ubiquitin

Using molecular dynamics simulation, we study the evaporation of water molecules off partially solvated ubiquitin. The evaporation and cooling rates are determined for a molecule at the initial temperature of 300 K. The cooling rate is found to be around 3 K/ns, and decreases with water temperature in the course of the evaporation. The conformation changes are monitored by studying a variety of intermediate partially solvated ubiquitin structures. We find that ubiquitin shrinks with decreasing hydration shell and exposes more of its hydrophilic surface area to the surrounding.

Comparing simulated and experimental translation and rotation constants: range of validity for viscosity scaling

Proper simulation of dynamic properties, including molecular diffusion, is an important goal of empirical force fields. However, the widely used TIP3P water model does not reproduce the experimental viscosity of water. Consequently, scaling of simulated diffusion constants of solutes in aqueous solutions is required to effectively compare them with experiment. It is proposed that scaling by the ratio of viscosities of model and real water is appropriate in the regime where the concentration dependence of simulated and experimental solution viscosities is parallel. With this ansatz, viscosity scaling can be carried out for glucose and trehalose up to 20 wt % for simulations carried out with the CHARMM additive carbohydrate force field C35 and TIP3P water; above this value, the concentration dependence of simulated viscosities lags that of experiment, and scaling is not advised. Scaled translational diffusion constants for glucose and the disaccharides trehalose, maltose, and melibiose at low concentration agree nearly quantitatively with experiment, as do NMR (13)C T(1)'s for glucose, trehalose, and maltose; these results support the use of C35 for simulations of sugar transport properties at low concentration. At high concentrations the scaled diffusion constants for glucose and trehalose underestimate and overestimate experiment, respectively. Hydrodynamic bead model calculations indicate a hydration level of approximately 1 water/hydroxyl for glucose. Patterns for the disaccharides are more complicated, though trehalose binds 0.5 to 1 more water than does maltose depending on the analysis.