CHARMM: the biomolecular simulation program

Predicting fragment binding poses using a combined MCSS MM-GBSA approach

Improved methods are required to predict the position and orientation (pose) of binding to the target protein of low molecular weight compounds identified in fragment screening campaigns. This is particularly important to guide initial chemistry to generate structure-activity relationships for the cases where a high resolution structure cannot be obtained. We have assessed the benefit of an implicit solvent method for assessment of fragment binding poses generated by the Multiple Copy Simultaneous Search (MCSS) method in CHARMm. Additionally, the effect of using multiple receptor structures for a flexible receptor is investigated. The original MCSS performance -50% of fragment positions accurately predicted and scored - was increased up to 67% by scoring MCSS energy minima with a Molecular Mechanics Generalized Born approach with molecular volume integration and Surface Area model (MM-GBSA). The same increase in performance (but occasionally for different targets) was observed when using the docking program GOLD followed by MM-GBSA rescoring. The combined results from both methods resulted in a higher success rate emphasizing that a comparison of different docking methods can increase the correct identification of binding poses. For a receptor where multiple structures are available, Hsp90, the average performance on randomly adding receptor structures was also investigated. The results suggest that predictions using these docking methods can be used with some confidence to guide chemical optimization, if the structure of the target either remains relatively fixed on ligand binding, or if a number of crystal structures are available with diverse ligands bound and there is information on the positions of key water molecules in the binding site.

Conformational relaxation and water penetration coupled to ionization of internal groups in proteins

Molecular dynamics simulations were used to examine the effects of ionization of internal groups on the structures of eighteen variants of staphylococcal nuclease (SNase) with internal Lys, Asp, or Glu. In most cases the RMSD values of internal ionizable side chains were larger when the ionizable moieties were charged than when they were neutral. Calculations of solvent-accessible surface area showed that the internal ionizable side chains were buried in the protein interior when they were neutral and moved toward crevices and toward the protein-water interface when they were charged. The only exceptions are Lys-36, Lys-62, and Lys-103, which remained buried even after charging. With the exception of Lys-38, the number of internal water molecules surrounding the ionizable group increased upon charging: the average number of water oxygen atoms within the first hydration shell increased by 1.7 for Lys residues, by 5.2 for Asp residues, and by 3.2 for Glu residues. The polarity of the microenvironment of the ionizable group also increased when the groups were charged: the average number of polar atoms of any kind within the first hydration shell increased by 2.7 for Lys residues, by 4.8 for Asp residues, and by 4.0 for Glu residues. An unexpected correlation was observed between the absolute value of the shifts in pK(a) values measured experimentally, and several parameters of structural relaxation: the net difference in the polarity of the microenvironment of the charged and neutral forms of the ionizable groups, the net difference in hydration of the charged and neutral forms of the ionizable groups, and the difference in RMSD values of the charged and neutral forms of the ionizable groups. The effects of ionization of internal groups on the conformation of the backbone were noticeable but mostly small and localized to the area immediately next to the internal ionizable moiety. Some variants did exhibit local unfolding.

Molecular modeling and molecular dynamics simulation studies of the GSK3β/ATP/substrate complex: understanding the unique P+4 primed phosphorylation specificity for GSK3β substrates

Substrate specificity of protein kinases is of fundamental importance for the integrity and fidelity of signaling pathways. Glycogen synthase kinase 3β (GSK3β) has a unique substrate specificity that prefers phosphorylation of its substrates at the P+4 serine before it can further phosphorylate the substrate at the P0 serine in the canonical motif SXXXS(p), where S(p) is the primed phosphorylation site. The detailed phosphorylation mechanism, however, is not clearly understood. In this study, a three-dimensional (3D) model of the ternary complex of GSK3β, ATP, and the phosphorylated glycogen synthase (pGS), termed GSK3β/ATP/pGS, is constructed using a hierarchical approach and by integrating molecular modeling and molecular dynamics (MD) simulations. Based on the 3D model, the substrate primed phosphorylation mechanism is investigated via two 12 ns comparative MD simulations of the GSK3β/ATP/pGS and GSK3β/ATP/GS systems, which differ in the phosphate group bound to the P+4 serine of GS. In agreement with structural analysis, computed binding free energies reveal that the binding of pGS to GSK3β is favored in the prephosphorylated state compared with the GS native state. More importantly, comparison with the system simulated without primed phosphorylation in the GSK3β/ATP/GS complex shows that for an optimal phosphorylation reaction to occur, the pGS priming phosphate in the GSK3β/ATP/pGS system optimizes the proper orientation of the GSK3β N- and C-terminal domains and clamps the P0 serine of pGS in the appropriate configuration for interaction with the ATP γ-phosphate within the catalytic groove.

Crystal structure of the HIV-2 neutralizing Fab fragment 7C8 with high specificity to the V3 region of gp125

7C8 is a mouse monoclonal antibody specific for the third hypervariable region (V3) of the human immunodeficiency virus type 2 (HIV-2)-associated protein gp125. The three-dimensional crystal structure of the Fab fragment of 7C8, determined to 2.7 Å resolution, reveals a deep and narrow antigen-binding cleft with architecture appropriate for an elongated epitope. The highly hydrophobic cleft is bordered on one side by the negatively charged second complementarity determining region (CDR2) and the unusually long positively charged CDR3 of the heavy chain and, on the other side, by the CDR1 of the light chain. Analysis of 7C8 in complex with molecular models of monomeric and trimeric gp125 highlights the importance of a conserved stretch of residues FHSQ that is localized centrally on the V3 region of gp125. Furthermore, modeling also indicates that the Fab fragment neutralizes the virus by sterically impairing subsequent engagement of the gp125 trimer with the co-receptor on the target cell.

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.

Solvation studies of a zinc finger protein in hydrated ionic liquids

The solvation of the zinc finger protein with the PDB-ID “5ZNF” in hydrated ionic liquids was studied at varying water content. 1-Ethyl-3-methylimidazolium and trifluoromethanesulfonate were the cation and anion, respectively. The protein stability as well as the solvation structure, the shell dynamics and the shell resolved dielectric properties were investigated by means of molecular dynamics simulations. The lengths of the respective trajectories extended up to 200 nanoseconds in order to cover the complete solvent dynamics. Considering the above mentioned properties as a function of the water content they all exhibit a maximum or minimum at the very same mole fraction. While the exact value x(H(2)O) = 0.927 depends on the underlying force field, its origin may be traced back to the competition between the van der Waals and the electrostatic energy of the protein as well as to the transition from aqueous dielectric screening to ionic charge screening with decreasing water content. The parameter-free Voronoi decomposition of space served as a basis for the analysis of most results. In particular, solvation shells were naturally inferred from this concept. In addition to the molecular analysis a mesoscopic view is given in terms of dielectric properties. Thereby, the net dielectric constant is decomposed into contributions from the protein, the first and second solvation shells as well as the bulk. Cross-terms between these components are given, too.

A full-length 3D structure for MAPK/ERK kinase 2 (MEK2)

As a pivotal signal pathway, the Ras/Raf/MEK/ERK cascade can be activated by multiple extracellular stimuli and can transmit signals to diverse substrates. It remains to be elucidated how so many different signals can be variously transferred by only two MEK molecules (MEK1 and MEK2). Because of technological limitations the complete structures of the MEKs are still unavailable. Here, we report the full-length structure of MEK2 obtained by homology modeling and molecular dynamics simulations. The simulations show that the N-terminal part of MEK2 is highly flexible and this flexibility may enable MEK2 to interact with ERKs and other ligands in diverse manners that correspond to various upstream signals and downstream consequences.

Nucleotide promiscuity of 3-phosphoglycerate kinase is in focus: implications for the design of better anti-HIV analogues

The wide specificity of 3-phosphoglycerate kinase (PGK) towards its nucleotide substrate is a property that allows contribution of this enzyme to the effective phosphorylation (i.e. activation) of nucleotide-based pro-drugs against HIV. Here, the structural basis of the nucleotide-PGK interaction is characterised in comparison to other kinases, namely pyruvate kinase (PK) and creatine kinase (CK), by enzyme kinetic analysis and structural modelling (docking) studies. The results provided evidence for favouring the purine vs. pyrimidine base containing nucleotides for PGK rather than for PK or CK. This is due to the exceptional ability of PGK in forming the hydrophobic contacts of the nucleotide rings that assures the appropriate positioning of the connected phosphate-chain for catalysis. As for the D-/L-configurations of the nucleotides, the L-forms (both purine and pyrimidine) are well accepted by PGK rather than either by PK or CK. Here again the dominance of the hydrophobic interactions of the L-form of pyrimidines with PGK is underlined in comparison with those of PK or CK. Furthermore, for the l-forms, the absence of the ribose OH-groups with PGK is better tolerated for the purine than for the pyrimidine containing compounds. On the other hand, the positioning of the phosphate-chain is an even more important term for PGK in the case of both purines and pyrimidines with an L-configuration, as deduced from the present kinetic studies with various nucleotide-site mutants of PGK. These characteristics of the kinase-nucleotide interactions can provide a guideline for designing new drugs.

Crystallographic and molecular dynamics analysis of loop motions unmasking the peptidoglycan-binding site in stator protein MotB of flagellar motor

Background

The C-terminal domain of MotB (MotB-C) shows high sequence similarity to outer membrane protein A and related peptidoglycan (PG)-binding proteins. It is believed to anchor the power-generating MotA/MotB stator unit of the bacterial flagellar motor to the peptidoglycan layer of the cell wall. We previously reported the first crystal structure of this domain and made a puzzling observation that all conserved residues that are thought to be essential for PG recognition are buried and inaccessible in the crystal structure. In this study, we tested a hypothesis that peptidoglycan binding is preceded by, or accompanied by, some structural reorganization that exposes the key conserved residues.

Methodology/Principal Findings

We determined the structure of a new crystalline form (Form B) of Helicobacter pylori MotB-C. Comparisons with the existing Form A revealed conformational variations in the petal-like loops around the carbohydrate binding site near one end of the β-sheet. These variations are thought to reflect natural flexibility at this site required for insertion into the peptidoglycan mesh. In order to understand the nature of this flexibility we have performed molecular dynamics simulations of the MotB-C dimer. The results are consistent with the crystallographic data and provide evidence that the three loops move in a concerted fashion, exposing conserved MotB residues that have previously been implicated in binding of the peptide moiety of peptidoglycan.

Conclusion/Significance

Our structural analysis provides a new insight into the mechanism by which MotB inserts into the peptidoglycan mesh, thus anchoring the power-generating complex to the cell wall.

MSCALE: A General Utility for Multiscale Modeling

The combination of theoretical models of macromolecules that exist at different spatial and temporal scales has become increasingly important for addressing complex biochemical problems. This work describes the extension of concurrent multiscale approaches, introduces a general framework for carrying out calculations, and describes its implementation into the CHARMM macromolecular modeling package. This functionality, termed MSCALE, generalizes both the additive and subtractive multiscale scheme (e.g. QM/MM ONIOM-type), and extends its support to classical force fields, coarse grained modeling (e.g. ENM, GNM, etc.), and a mixture of them all. The MSCALE scheme is completely parallelized with each subsystem running as an independent, but connected calculation. One of the most attractive features of MSCALE is the relative ease of implementation using the standard MPI communication protocol. This allows external access to the framework and facilitates the combination of functionality previously isolated in separate programs. This new facility is fully integrated with free energy perturbation methods, Hessian based methods, and the use of periodicity and symmetry, which allows the calculation of accurate pressures. We demonstrate the utility of this new technique with four examples; (1) subtractive QM/MM and QM/QM calculations; (2) multi-force field alchemical free energy perturbation; (3) integration with the SANDER module of AMBER and the TINKER package to gain access to potentials not available in CHARMM; and (4) mixed resolution (i.e. coarse grain / all-atom) normal mode analysis. The potential of this new tool is clearly established and in conclusion an interesting mathematical problem is highlighted and future improvements are proposed.

GRIFFIN: A versatile methodology for optimization of protein-lipid interfaces for membrane protein simulations

As new atomic structures of membrane proteins are resolved, they reveal increasingly complex transmembrane topologies, and highly irregular surfaces with crevices and pores. In many cases, specific interactions formed with the lipid membrane are functionally crucial, as is the overall lipid composition. Compounded with increasing protein size, these characteristics pose a challenge for the construction of simulation models of membrane proteins in lipid environments; clearly, that these models are sufficiently realistic bears upon the reliability of simulation-based studies of these systems. Here, we introduce GRIFFIN, which uses a versatile framework to automate and improve a widely-used membrane-embedding protocol. Initially, GRIFFIN carves out lipid and water molecules from a volume equivalent to that of the protein, so as to conserve the system density. In the subsequent optimization phase GRIFFIN adds an implicit grid-based protein force-field to a molecular dynamics simulation of the pre-carved membrane. In this force-field, atoms inside the implicit protein volume experience an outward force that will expel them from that volume, whereas those outside are subject to electrostatic and van-der-Waals interactions with the implicit protein. At each step of the simulation, these forces are updated by GRIFFIN and combined with the intermolecular forces of the explicit lipid-water system. This procedure enables the construction of realistic and reproducible starting configurations of the protein-membrane interface within a reasonable timeframe and with minimal intervention. GRIFFIN is a standalone tool designed to work alongside any existing molecular dynamics package, such as NAMD or GROMACS.

Cy3-DNA stacking interactions strongly depend on the identity of the terminal basepair

We characterized the effect of the first basepair on the conformational dynamics of the fluorescent dye Cy3 attached to the 5' end of double-stranded DNA using gaussian-mixture adaptive umbrella sampling simulations. In the simulations, the sampling of all five dihedral angles along the linker was enhanced, so that both stacked and unstacked states were sampled. The affinity of Cy3 for a T·A basepair (with the dye attached to T) was found to be significantly less than for the other basepairs. This was verified experimentally by measuring the activation energies for cis-trans isomerization of the dye. The simulation and experimental results indicate the existence of partially unstacked conformations amenable to photoisomerization. The simulations also showed that stacking of Cy3 straightens the DNA while stabilizing the first basepair. Our findings indicate that fluorescence is modulated by Cy3-DNA interactions in a sequence-dependent manner.

Tropomyosin position on F-actin revealed by EM reconstruction and computational chemistry

Electron microscopy and fiber diffraction studies of reconstituted F-actin-tropomyosin filaments reveal the azimuthal position of end-to-end linked tropomyosin molecules on the surface of actin. However, the longitudinal z-position of tropomyosin along F-actin is still uncertain. Without this information, atomic models of F-actin-tropomyosin filaments, free of constraints imposed by troponin or other actin-binding proteins, cannot be formulated, and thus optimal interfacial contacts between actin and tropomyosin remain unknown. Here, a computational search assessing electrostatic interactions for multiple azimuthal locations, z-positions, and pseudo-rotations of tropomyosin on F-actin was performed. The information gleaned was used to localize tropomyosin on F-actin, yielding an atomic model characterized by protein-protein contacts that primarily involve clusters of basic amino acids on actin subdomains 1 and 3 juxtaposed against acidic residues on the successive quasi-repeating units of tropomyosin. A virtually identical model generated by docking F-actin and tropomyosin atomic structures into electron microscopy reconstructions of F-actin-tropomyosin validated the above solution. Here, the z-position of tropomyosin alongside F-actin was defined by matching the seven broad and narrow motifs that typify tropomyosin's twisting superhelical coiled-coil to the wide and tapering tropomyosin densities seen in surface views of F-actin-tropomyosin reconstructions. The functional implications of the F-actin-tropomyosin models determined in this work are discussed.

Design of a switchable eliminase

The active sites of enzymes are lined with side chains whose dynamic, geometric, and chemical properties have been finely tuned relative to the corresponding residues in water. For example, the carboxylates of glutamate and aspartate are weakly basic in water but become strongly basic when dehydrated in enzymatic sites. The dehydration of the carboxylate, although intrinsically thermodynamically unfavorable, is achieved by harnessing the free energy of folding and substrate binding to reach the required basicity. Allosterically regulated enzymes additionally rely on the free energy of ligand binding to stabilize the protein in a catalytically competent state. We demonstrate the interplay of protein folding energetics and functional group tuning to convert calmodulin (CaM), a regulatory binding protein, into AlleyCat, an allosterically controlled eliminase. Upon binding Ca(II), native CaM opens a hydrophobic pocket on each of its domains. We computationally identified a mutant that (i) accommodates carboxylate as a general base within these pockets, (ii) interacts productively in the Michaelis complex with the substrate, and (iii) stabilizes the transition state for the reaction. Remarkably, a single mutation of an apolar residue at the bottom of an otherwise hydrophobic cavity confers catalytic activity on calmodulin. AlleyCat showed the expected pH-rate profile, and it was inactivated by mutation of its active site Glu to Gln. A variety of control mutants demonstrated the specificity of the design. The activity of this minimal 75-residue allosterically regulated catalyst is similar to that obtained using more elaborate computational approaches to redesign complex enzymes to catalyze the Kemp elimination reaction.

Enantiospecificity of chloroperoxidase-catalyzed epoxidation: biased molecular dynamics study of a cis-β-methylstyrene/chloroperoxidase-compound I complex

Molecular dynamics simulations of an explicitly solvated cis-β-methylstyrene/chloroperoxidase-Compound I complex are performed to determine the cause of the high enantiospecificity of epoxidation. From the simulations, a two-dimensional free energy potential is calculated to distinguish binding potential wells from which reaction to 1S2R and 1R2S epoxide products may occur. Convergence of the free energy potential is accelerated with an adaptive biasing potential. Analysis of binding is followed by analysis of 1S2R and 1R2S reaction precursor structures in which the substrate, having left the binding wells, places its reactive double bond in steric proximity to the oxyferryl heme center. Structural analysis of binding and reaction precursor conformations is presented. We find that 1), a distortion of Glu(183) is important for CPO-catalyzed epoxidation as was postulated previously based on experimental results; 2), the free energy of binding does not provide significant differentiation between structures leading to the respective epoxide enantiomers; and 3), CPO's enantiospecificity toward cis-β-methylstyrene is likely to be caused by a specific group of residues which form a hydrophobic core surrounding the oxyferryl heme center.

Symmetry-based self-assembled nanotubes constructed using native protein structures: the key role of flexible linkers

We construct nanotubes using native protein structures and their native associations from structural databases. The construction is based on a shape-guided symmetric self-assembly concept. Our strategy involves fusing judiciously-selected oligomerization domains via peptide linkers. Linkers are inherently flexible, hence their choice is critical: they should position the domains in three-dimensional space in the desired orientation while retaining their own natural conformational tendencies; however, at the same time, retain the construct stability. Here we outline a design scheme which accounts for linker flexibility considerations, and present two examples. The first is HIV-1 capsid protein, which in vitro self-assembles into nanotubes and conical capsids, and its linker exists as a short flexible loop. The second involves novel nanotubes construction based on antimicrobial homodimer Magainin 2, employing linkers of distinct lengths and flexibility levels. Our strategy utilizes the abundance of unique shapes and sizes of proteins and their building blocks which can assemble into a vast number of combinations, and consequently, nanotubes of distinct morphologies and diameters. Computational design and assessment methodologies can help reduce the number of candidates for experimental validation. This is an invited paper for a special issue on protein dynamics, here focusing on flexibility in nanotube design based on protein building blocks.

Effect of flanking residues on the conformational sampling of the internal fusion peptide from Ebola virus

Fusion peptides mediate viral and host-cell membrane fusion during viral entry. The monomeric form of the internal fusion peptide from Ebola virus was studied in membrane bilayer and water environments with computer simulations using replica exchange sampling and an implicit solvent description of the environment. Wild-type Ebola fusion peptide (EFP), the W8A mutant form, and an extended construct with flanking residues were examined. It was found that the monomeric form of wild-type EFP adopts coil-helix-coil structure with a short helix from residues 8 to 11 mostly sampling orientations parallel to the membrane surface. W8A mutation disrupts the helicity in the N-terminal region of the peptide and leads to a preference for slightly oblique orientation relative to the membrane surface. The addition of flanking residues also alters the fusion peptide conformation with either a helix-break-helix structure or extended N and C-termini and reduced membrane insertion. In water, the fusion peptide is found to adopt structures with low helicity.

Examination of the α-chitin structure and decrystallization thermodynamics at the nanoscale

Chitin is the primary structural material of insect and crustacean exoskeletons and fungal and algal cell walls, and as such it is the one of the most abundant biological materials on Earth. Chitin forms linear polymers of β1,4-linked-N-acetyl-D-glucosamine (GlcNAc), and in Nature, enzyme cocktails deconstruct chitin to GlcNAc. The mechanism of chitin deconstruction, like that of cellulose deconstruction, has been under investigation due to its importance in the global carbon cycle and in production of renewable and sustainable products from biological matter. To further understand the nanoscale properties of chitin, here we simulate crystals of α-chitin, which is the most prevalent form in Nature. We find excellent agreement with the recently reported crystal structure and we report the salient features of the simulations related to crystalline stability. We also compute the thermodynamic work required to peel individual chains from α-chitin surfaces, which a chitinase enzyme must conduct to deconstruct chitin. Compared with previous simulations of native plant cellulose Iβ, α-chitin exhibits higher decrystallization work for chains in the middle of surfaces and similar work for chains on the edges of crystals. Unlike cellulose, the free energy profile is dominated by a single bifurcated hydrogen bond between chains formed by the GlcNAc side chains and the O6 atoms on the primary alcohol group. This study highlights the molecular features of chitin that make it such a tough, recalcitrant material, and provides a key thermodynamic parameter in our quantitative understanding of how enzymes contribute to the turnover of carbohydrates in the biosphere.

Observation of poly(ethylene glycol) clusters with the chlorine anion in the gas phase under electrospray conditions

It is demonstrated herein that poly(ethylene glycol) (PEG) oligomers can form stable complexes with the chlorine anion in the gas phase as evidenced by results from electrospray ionization mass spectrometry (ESI-MS) and molecular dynamics simulation. While the formation of crown-ether-like structures by acyclic polyethers in their complexes with alkali metal cations coordinated by the ether oxygen atoms has been extensively studied, the possibility of forming 'inversed' quasi-cyclic structures able to bind a monoatomic anion has not been proved till now. We have observed the formation of stable gas-phase complexes of oligomers of PEG-400 with the Cl(-) anion experimentally by ESI-MS for the first time. It is suggested that a necessary precondition for obtaining the polyether-chlorine anion clusters is the prevention of the formation of neutral ion pairs. Molecular dynamics simulation has demonstrated the wrapping of the Cl(-) anion by the PEG chain, to stabilize the PEG(n)•Cl(-) clusters in the gas phase. The conformation of the polyether chain in such quasi-cyclic or quasi-helical complexes is 'inversed' compared with that in the complexes with cations: that is its hydrogen atoms are turned towards the central anion. Awareness of the possibility of the Cl(-) anion being trapped in quasi-cyclic PEG structures may be of practical importance when considering the intermolecular interactions of PEGs.

Understanding the determinants of selectivity in drug metabolism through modeling of dextromethorphan oxidation by cytochrome P450

Cytochrome P450 enzymes play key roles in the metabolism of the majority of drugs. Improved models for prediction of likely metabolites will contribute to drug development. In this work, two possible metabolic routes (aromatic carbon oxidation and O-demethylation) of dextromethorphan are compared using molecular dynamics (MD) simulations and density functional theory (DFT). The DFT results on a small active site model suggest that both reactions might occur competitively. Docking and MD studies of dextromethorphan in the active site of P450 2D6 show that the dextromethorphan is located close to heme oxygen in a geometry apparently consistent with competitive metabolism. In contrast, calculations of the reaction path in a large protein model [using a hybrid quantum mechanical-molecular mechanics (QM/MM) method] show a very strong preference for O-demethylation, in accordance with experimental results. The aromatic carbon oxidation reaction is predicted to have a high activation energy, due to the active site preventing formation of a favorable transition-state structure. Hence, the QM/MM calculations demonstrate a crucial role of many active site residues in determining reactivity of dextromethorphan in P450 2D6. Beyond substrate binding orientation and reactivity of Compound I, successful metabolite predictions must take into account the detailed mechanism of oxidation in the protein. These results demonstrate the potential of QM/MM methods to investigate specificity in drug metabolism.

Development of molecular simulation methods to accurately represent protein-surface interactions: The effect of pressure and its determination for a system with constrained atoms

When performing molecular dynamics simulations for a system with constrained (fixed) atoms, traditional isobaric algorithms (e.g., NPT simulation) often cannot be used. In addition, the calculation of the internal pressure of a system with fixed atoms may be highly inaccurate due to the nonphysical nature of the atomic constraints and difficulties in accurately defining the volume occupied by the unconstrained atoms in the system. The inability to properly set and control pressure can result in substantial problems for the accurate simulation of condensed-phase systems if the behavior of the system (e.g., peptide/protein adsorption) is sensitive to pressure. To address this issue, the authors have developed an approach to accurately determine the internal pressure for a system with constrained atoms. As the first step in this method, a periodically extendable portion of the mobile phase of the constrained system (e.g., the solvent atoms) is used to create a separate unconstrained system for which the pressure can be accurately calculated. This model system is then used to create a pressure calibration plot for an intensive local effective virial parameter for a small volume cross section or "slab" of the system. Using this calibration plot, the pressure of the constrained system can then be determined by calculating the virial parameter for a similarly sized slab of mobile atoms. In this article, the authors present the development of this method and demonstrate its application using the CHARMM molecular simulation program to characterize the adsorption behavior of a peptide in explicit water on a hydrophobic surface whose lattice spacing is maintained with atomic constraints. The free energy of adsorption for this system is shown to be dramatically influenced by pressure, thus emphasizing the importance of properly maintaining the pressure of the system for the accurate simulation of protein-surface interactions.

"Fluctuograms" reveal the intermittent intra-protein communication in subtilisin Carlsberg and correlate mechanical coupling with co-evolution

The mechanism of intra-protein communication and allosteric coupling is key to understanding the structure-property relationship of protein function. For subtilisin Carlsberg, the Ca²⁺-binding loop is distal to substrate-binding and active sites, yet the serine protease function depends on Ca²⁺ binding. The atomic molecular dynamics (MD) simulations of apo and Ca²⁺-bound subtilisin show similar structures and there is no direct evidence that subtilisin has alternative conformations. To model the intra-protein communication due to Ca²⁺ binding, we transform the sequential segments of an atomic MD trajectory into separate elastic network models to represent anharmonicity and nonlinearity effectively as the temporal and spatial variation of the mechanical coupling network. In analogy to the spectrogram of sound waves, this transformation is termed the “fluctuogram” of protein dynamics. We illustrate that the Ca²⁺-bound and apo states of subtilisin have different fluctuograms and that intra-protein communication proceeds intermittently both in space and in time. We found that residues with large mechanical coupling variation due to Ca²⁺ binding correlate with the reported mutation sites selected by directed evolution for improving the stability of subtilisin and its activity in a non-aqueous environment. Furthermore, we utilize the fluctuograms calculated from MD to capture the highly correlated residues in a multiple sequence alignment. We show that in addition to the magnitude, the variance of coupling strength is also an indicative property for the sequence correlation observed in a statistical coupling analysis. The results of this work illustrate that the mechanical coupling networks calculated from atomic details can be used to correlate with functionally important mutation sites and co-evolution.

A hallmark of protein molecules is their machine-like behaviors while carrying out biological functions. At the molecular level, molecular signals such as binding a metal ion at an action site can cause long-range effects and alter protein function. Such phenomena are often referred to as intra-protein communication or allosteric coupling. Elucidating the underlying mechanisms could lead to novel discovery of molecular modulators to regulate protein function in a more specific and effective manner. A long-standing puzzle is the roles of the anharmonicity and nonlinearity in protein dynamics. To incorporate these characters in modeling intra-protein communication, we devise a “fluctuogram” analysis to record the choreography of allosteric coupling in an atomic molecular dynamics simulation. We show that fluctuogram analysis can bridge the results of physics-based simulation and sequence alignment in bioinformatics by capturing the residues that exhibit high correlation in a multiple sequence alignment. We also show that the fluctuograms calculated from atomic details have the potential to be applied as a tool to select mutation sites for modulating protein function.

UV resonance Raman study of TTR(105-115) structural evolution as a function of temperature

UV resonance Raman spectroscopy was used to probe the temperature dependence of the conformation of TTR(105-115) in solution. Resonance Raman spectra with excitation at 239.5 nm, show an increase in the absolute resonance Raman cross section of Tyr with an increase in temperature. This trend is associated with an increase in the hydrophobicity of the Tyr local environment, suggesting a conformational change at 28 °C. Excitation at ~200 nm is known to enhance scattering due to amide vibrations and provides insights as to the secondary structure of a peptide or protein. UVRR spectra at this excitation suggest that in solution the peptide assumes a disordered conformation with frequent formation of β-turns. Explicit-solvent replica-exchange MD simulations of the isolated peptide in the region 15 to 37 °C suggest that the dominant conformation assumed by the peptide corresponds to a coil with β-turns in the central and C-terminal region. In line with the experiments, an increase in temperature induces structural order in the peptide, reflected by an increase in the probability for the formation of β-turns and hydrophobic side-chain contacts, mainly in the 8-11 moiety, and to a lesser extent in the 4-7 moiety.

Evaluation of spectrophotometric and HPLC methods for shikimic acid determination in plants: models in glyphosate-resistant and -susceptible crops

Endogenous shikimic acid determinations are routinely used to assess the efficacy of glyphosate in plants. Numerous analytical methods exist in the public domain for the detection of shikimic acid, yet the most commonly cited comprise spectrophotometric and high-pressure liquid chromatography (HPLC) methods. This paper compares an HPLC and two spectrophotometric methods (Spec 1 and Spec 2) and assesses the effectiveness in the detection of shikimic acid in the tissues of glyphosate-treated plants. Furthermore, the study evaluates the versatility of two acid-based shikimic acid extraction methods and assesses the longevity of plant extract samples under different storage conditions. Finally, Spec 1 and Spec 2 are further characterized with respect to (1) the capacity to discern between shikimic acid and chemically related alicyclic hydroxy acids, (2) the stability of the chromophore (t1/2), (3) the detection limits, and (4) the cost and simplicity of undertaking the analytical procedure. Overall, spectrophotometric methods were more cost-effective and simpler to execute yet provided a narrower detection limit compared to HPLC. All three methods were specific to shikimic acid and detected the compound in the tissues of glyphosate-susceptible crops, increasing exponentially in concentration within 24 h of glyphosate application and plateauing at approximately 72 h. Spec 1 estimated more shikimic acid in identical plant extract samples compared to Spec 2 and, likewise, HPLC detection was more effective than spectrophotometric determinations. Given the unprecedented global adoption of glyphosate-resistant crops and concomitant use of glyphosate, an effective and accurate assessment of glyphosate efficacy is important. Endogenous shikimic acid determinations are instrumental in corroborating the efficacy of glyphosate and therefore have numerous applications in herbicide research and related areas of science as well as resolving many commercial issues as a consequence of glyphosate utilization.

Molecular-level origins of biomass recalcitrance: decrystallization free energies for four common cellulose polymorphs

Cellulose is a crystalline polymer of β1,4-D-glucose that is difficult to deconstruct to sugars by enzymes. The recalcitrance of cellulose microfibrils is a function of both the shape of cellulose microfibrils and the intrinsic work required to decrystallize individual chains, the latter of which is calculated here from the surfaces of four crystalline cellulose polymorphs: cellulose Iβ, cellulose Iα, cellulose II, and cellulose III(I). For edge chains, the order of decrystallization work is as follows (from highest to lowest): Iβ, Iα, ΙΙΙ(Ι), and II. For cellulose Iβ, we compare chains from three different locations on the surface and find that an increasing number of intralayer hydrogen bonds (from 0 to 2) increases the intrinsic decrystallization work. From these results, we propose a microkinetic model for the deconstruction of cellulose (and chitin) by processive enzymes, which when taken with a previous study [Horn et al. Proc. Natl. Acad. Sci. U.S.A. 2006, 103, 18089] identifies the thermodynamic and kinetic attributes of enzyme and substrate engineering for enhanced cellulose (or chitin) conversion. Overall, this study provides new insights into the molecular interactions that form the structural basis of cellulose, which is the primary building block of plant cell walls, and highlights the need for experimentally determining microfibril shape at the nanometer length scale when comparing conversion rates of cellulose polymorphs by enzymes.

Enhancing QM/MM molecular dynamics sampling in explicit environments via an orthogonal-space-random-walk-based strategy

Accurate prediction of molecular conformations in explicit environments, such as aqueous solution and protein interiors, can facilitate our understanding of various molecular recognition processes. Most computational approaches are limited as a result of their compromised choices between the underlying energy model and the sampling length. Taking advantage of a recent second-order generalized ensemble scheme [e.g., the orthogonal space random walk (OSRW) strategy], which can synergistically accelerate the motion of a focused region and its coupled environmental response, we are presenting a QM/MM (combined quantum mechanical/molecular mechanical)-based molecular dynamics sampling technique to explore molecular conformational landscapes in explicit environments. The present QM/MM potential scaling-based OSRW sampling scheme is employed to study the binding of DMSO to the FKBP12 protein, the conformation distribution of a novel mercaptosulfonamide inhibitor in aqueous solution, and its binding poses in zinc-containing matrix metalloproteinase-9 (MMP-9). As demonstrated, the present QM/MM second-order generalized ensemble sampling technique enables feasible usage of the QM/MM model to sample molecular conformations in condensed environments.

Unsuspected pathway of the allosteric transition in hemoglobin

Large conformational transitions play an essential role in the function of many proteins, but experiments do not provide the atomic details of the path followed in going from one end structure to the other. For the hemoglobin tetramer, the transition path between the unliganded (T) and tetraoxygenated (R) structures is not known, which limits our understanding of the cooperative mechanism in this classic allosteric system, where both tertiary and quaternary changes are involved. The conjugate peak refinement algorithm is used to compute an unbiased minimum energy path at atomic detail between the two end states. Although the results confirm some of the proposals of Perutz [Perutz MF (1970) Stereochemistry of cooperative effects in haemoglobin. Nature 228:726-734], the subunit motions do not follow the textbook description of a simple rotation of one αβ-dimer relative to the other. Instead, the path consists of two sequential quaternary rotations, each involving different subdomains and axes. The quaternary transitions are preceded and followed by phases of tertiary structural changes. The results explain the recent photodissociation measurements, which suggest that the quaternary transition has a fast (2 μs) as well as a slow (20 μs) component and provide a testable model for single molecule FRET experiments.

Understanding the physical basis of the salt dependence of the electrostatic binding free energy of mutated charged ligand-nucleic acid complexes

The predictions of the derivative of the electrostatic binding free energy of a biomolecular complex, ΔG(el), with respect to the logarithm of the 1:1 salt concentration, d(ΔG(el))/d(ln[NaCl]), SK, by the Poisson-Boltzmann equation, PBE, are very similar to those of the simpler Debye-Hückel equation, DHE, because the terms in the PBE's predictions of SK that depend on the details of the dielectric interface are small compared to the contributions from long-range electrostatic interactions. These facts allow one to obtain predictions of SK using a simplified charge model along with the DHE that are highly correlated with both the PBE and experimental binding data. The DHE-based model developed here, which was derived from the generalized Born model, explains the lack of correlation between SK and ΔG(el) in the presence of a dielectric discontinuity, which conflicts with the popular use of this supposed correlation to parse experimental binding free energies into electrostatic and nonelectrostatic components. Moreover, the DHE model also provides a clear justification for the correlations between SK and various empirical quantities, like the number of ion pairs, the ligand charge on the interface, the Coulomb binding free energy, and the product of the charges on the complex's components, but these correlations are weak, questioning their usefulness.

Insertion domain within mammalian mitochondrial translation initiation factor 2 serves the role of eubacterial initiation factor 1

Mitochondria have their own translational machineries for the synthesis of thirteen polypeptide chains that are components of the complexes that participate in the process of oxidative phosphorylation (or ATP generation). Translation initiation in mammalian mitochondria requires two initiation factors, IF2(mt) and IF3(mt), instead of the three that are present in eubacteria. The mammalian IF2(mt) possesses a unique 37 amino acid insertion domain, which is known to be important for the formation of the translation initiation complex. We have obtained a three-dimensional cryoelectron microscopic map of the mammalian IF2(mt) in complex with initiator fMet-tRNA(iMet) and the eubacterial ribosome. We find that the 37 amino acid insertion domain interacts with the same binding site on the ribosome that would be occupied by the eubacterial initiation factor IF1, which is absent in mitochondria. Our finding suggests that the insertion domain of IF2(mt) mimics the function of eubacterial IF1, by blocking the ribosomal aminoacyl-tRNA binding site (A site) at the initiation step.

The O-glycosylated linker from the Trichoderma reesei Family 7 cellulase is a flexible, disordered protein

Fungi and bacteria secrete glycoprotein cocktails to deconstruct cellulose. Cellulose-degrading enzymes (cellulases) are often modular, with catalytic domains for cellulose hydrolysis and carbohydrate-binding modules connected by linkers rich in serine and threonine with O-glycosylation. Few studies have probed the role that the linker and O-glycans play in catalysis. Since different expression and growth conditions produce different glycosylation patterns that affect enzyme activity, the structure-function relationships that glycosylation imparts to linkers are relevant for understanding cellulase mechanisms. Here, the linker of the Trichoderma reesei Family 7 cellobiohydrolase (Cel7A) is examined by simulation. Our results suggest that the Cel7A linker is an intrinsically disordered protein with and without glycosylation. Contrary to the predominant view, the O-glycosylation does not change the stiffness of the linker, as measured by the relative fluctuations in the end-to-end distance; rather, it provides a 16 Å extension, thus expanding the operating range of Cel7A. We explain observations from previous biochemical experiments in the light of results obtained here, and compare the Cel7A linker with linkers from other cellulases with sequence-based tools to predict disorder. This preliminary screen indicates that linkers from Family 7 enzymes from other genera and other cellulases within T. reesei may not be as disordered, warranting further study.

Unfolding of the amyloid β-peptide central helix: mechanistic insights from molecular dynamics simulations

Alzheimer's disease (AD) pathogenesis is associated with formation of amyloid fibrils caused by polymerization of the amyloid β-peptide (Aβ), which is a process that requires unfolding of the native helical structure of Aβ. According to recent experimental studies, stabilization of the Aβ central helix is effective in preventing Aβ polymerization into toxic assemblies. To uncover the fundamental mechanism of unfolding of the Aβ central helix, we performed molecular dynamics simulations for wild-type (WT), V18A/F19A/F20A mutant (MA), and V18L/F19L/F20L mutant (ML) models of the Aβ central helix. It was quantitatively demonstrated that the stability of the α-helical conformation of both MA and ML is higher than that of WT, indicating that the α-helical propensity of the three nonpolar residues (18, 19, and 20) is the main factor for the stability of the whole Aβ central helix and that their hydrophobicity plays a secondary role. WT was found to completely unfold by a three-step mechanism: 1) loss of α-helical backbone hydrogen bonds, 2) strong interactions between nonpolar sidechains, and 3) strong interactions between polar sidechains. WT did not completely unfold in cases when any of the three steps was omitted. MA and ML did not completely unfold mainly due to the lack of the first step. This suggests that disturbances in any of the three steps would be effective in inhibiting the unfolding of the Aβ central helix. Our findings would pave the way for design of new drugs to prevent or retard AD.

The adipophilin C terminus is a self-folding membrane-binding domain that is important for milk lipid secretion

Cytoplasmic lipid droplets (CLD) in mammary epithelial cells undergo secretion by a unique membrane envelopment process to produce milk lipids. Adipophilin (ADPH/Plin2), a member of the perilipin/PAT family of lipid droplet-associated proteins, is hypothesized to mediate CLD secretion through interactions with apical plasma membrane elements. We found that the secretion of CLD coated by truncated ADPH lacking the C-terminal region encoding a putative four-helix bundle structure was impaired relative to that of CLD coated by full-length ADPH. We used homology modeling and analyses of the solution and membrane binding properties of purified recombinant ADPH C terminus to understand how this region possibly mediates CLD secretion. Homology modeling supports the concept that the ADPH C terminus forms a four-helix bundle motif and suggests that this structure can form stable membrane bilayer interactions. Circular dichroism and protease mapping studies confirmed that the ADPH C terminus is an independently folding α-helical structure that is relatively resistant to urea denaturation. Liposome binding studies showed that the purified C terminus binds to phospholipid membranes through electrostatic dependent interactions, and cell culture studies documented that it localizes to the plasma membrane. Collectively, these data provide direct evidence that the ADPH C terminus forms a stable membrane binding helical structure that is important for CLD secretion. We speculate that interactions between the four-helix bundle of ADPH and membrane phospholipids may be an initial step in milk lipid secretion.

Modeling of ligand binding to G protein coupled receptors: cannabinoid CB1, CB2 and adrenergic β 2 AR

Cannabinoid and adrenergic receptors belong to the class A (similar to rhodopsin) G protein coupled receptors. Docking of agonists and antagonists to CB(1) and CB(2) cannabinoid receptors revealed the importance of a centrally located rotamer toggle switch and its possible participation in the mechanism of agonist/antagonist recognition. The switch is composed of two residues, F3.36 and W6.48, located on opposite transmembrane helices TM3 and TM6 in the central part of the membranous domain of cannabinoid receptors. The CB(1) and CB(2) receptor models were constructed based on the adenosine A(2A) receptor template. The two best scored conformations of each receptor were used for the docking procedure. In all poses (ligand-receptor conformations) characterized by the lowest ligand-receptor intermolecular energy and free energy of binding the ligand type matched the state of the rotamer toggle switch: antagonists maintained an inactive state of the switch, whereas agonists changed it. In case of agonists of β(2)AR, the (R,R) and (S,S) stereoisomers of fenoterol, the molecular dynamics simulations provided evidence of different binding modes while preserving the same average position of ligands in the binding site. The (S,S) isomer was much more labile in the binding site and only one stable hydrogen bond was created. Such dynamical binding modes may also be valid for ligands of cannabinoid receptors because of the hydrophobic nature of their ligand-receptor interactions. However, only very long molecular dynamics simulations could verify the validity of such binding modes and how they affect the process of activation.

The hidden energetics of ligand binding and activation in a glutamate receptor

Ionotropic glutamate receptors (iGluRs) are ligand-gated ion channels that mediate most excitatory synaptic transmission in the central nervous system. The free energy of neurotransmitter binding to the ligand-binding domains (LBDs) of iGluRs is converted into useful work to drive receptor activation. We have computed the principal thermodynamic contributions from ligand docking and ligand-induced closure of LBDs for nine ligands of GluA2 using all-atom molecular dynamics free energy simulations. We have validated the results by comparison with experimentally measured apparent affinities to the isolated LBD. Features in the free energy landscapes that govern closure of LBDs are key determinants of binding free energies. An analysis of accessible LBD conformations transposed into the context of an intact GluA2 receptor revealed that the relative displacement of specific diagonal subunits in the tetrameric structure may be key to the action of partial agonists.

Recent advances in ligand-based drug design: relevance and utility of the conformationally sampled pharmacophore approach

In the absence of three-dimensional (3D) structures of potential drug targets, ligand-based drug design is one of the popular approaches for drug discovery and lead optimization. 3D structure-activity relationships (3D QSAR) and pharmacophore modeling are the most important and widely used tools in ligand-based drug design that can provide crucial insights into the nature of the interactions between drug target and ligand molecule and provide predictive models suitable for lead compound optimization. This review article will briefly discuss the features and potential application of recent advances in ligand-based drug design, with emphasis on a detailed description of a novel 3D QSAR method based on the conformationally sample pharmacophore (CSP) approach (denoted CSP-SAR). In addition, data from a published study are used to compare the CSP-SAR approach to the Catalyst method, emphasizing the utility of the CSP approach for ligand-based model development.

A computational study of the mechanism of the selective crystallization of α- and β-glycine from water and methanol-water mixture

Understanding the control of polymorphism in organic crystals is of paramount importance to the pharmaceutical, chemical, and food industries. In this work, we investigated two mechanisms described in the literature about the selective crystallization of α- and β-glycine from water and from mixtures of water and methanol using molecular simulations. The link hypothesis (J. Phys. Chem. B 2008, 112, 7794; Cryst. Growth Des. 2006, 6, 1788; J. Inclusion Phenom. Mol. Recognit. Chem. 1990, 8, 395; J. Am. Chem. Soc. 1986, 108, 5871.), which tries to relate the structure of the polymorph obtained from crystallization to the structure of the prenucleation aggregates in the solutions, says the abundance of glycine cyclic dimers in aqueous solutions leads to the crystallization of α-glycine, the polymorph using cyclic dimers as the packing units. This hypothesis was studied first. We revisited the self-assembly of glycine molecules in solution using molecular dynamics to address the debate (Phys. Rev. Lett. 2007, 99, 115702; J. Phys. Chem. B 2008, 112, 7280; J. Am. Chem. Soc. 2008, 130, 13973.) about which is the dominating species in the glycine aqueous solutions and whether there is a link between the solution chemistry and the polymorphic outcome of crystallization. The structures of the glycine clusters were characterized using a structural parameter called cyclic dimer fraction. The glycine clusters in methanol-water mixtures have higher cyclic dimer compositions than those in the pure aqueous solutions. Moreover, the glycine open-chain dimer is more stable than the cyclic dimer regardless of the presence of methanol. All these suggest that the link hypothesis does not work for the polymorphic system of glycine, and the selective crystallization of α- and β-glycine from water and methanol-water mixture, respectively, is not due to the abundance of glycine aggregates in the solution phase with a similar structure to the crystallizing solid form. The hypothesis of the methanol inhibition on the growth of α-glycine {010} and {010} faces, proposed by Weissbuch (Angew. Chem., Int. Ed. 2005, 44, 3226.), was also studied. The interfaces between the {010} and {010} faces of both crystal forms (α and β) and both solvents (water and methanol-water 3:7 mixture) were studied using molecular simulation. No strong binding of methanol onto the {010} and {010} faces of both crystal forms was observed, and the addition of methanol dilutes the crystal-solvent interactions on all faces. Therefore, the selective crystallization of β and α-glycine with and without methanol does not follow either of the two mechanisms in the literature.

Calculation of the Vibrational Stark Effect Using a First-Principles QM/MM Approach

The proper description of the electric environment of condensed phases is a critical challenge for force field methods. To test and validate the ability of the CHARMM additive force field to describe the electric environment in aqueous solution combined QM/MM calculations have been used to calculate the vibrational Stark effect (VSE). We utilized a first principles methodology using correlated electronic structure techniques to compute the Stark shift between the gas phase and solvent environments and between two different solvent environments of three VSE probes containing acetonitrile or fluorine functionalities which have been well-characterized experimentally. Reasonable agreement with the experimentally determined Stark shifts is obtained when the MM atoms are described by the CHARMM additive force field, though it is essential to employ an anharmonic correction in the frequency calculation. In addition, the electric field created by the solvent is computed along the CN bond and a theoretical Stark tuning rate is determined for acetonitrile and shown to be in satisfactory agreement with experiment.