Generating stereochemically acceptable protein pathways

We describe a new method for rapidly generating stereochemically acceptable pathways in proteins. The method, called geometric targeting, is publicly available at the webserver http://pathways.asu.edu , and includes tools for visualization of the pathway and creating movie files for use in presentations. The user submits an initial structure and a target structure, and a pathway between the two input states is generated automatically. Besides visualization, the structural quality of the pathways makes them useful as input pathways into pathway refinement techniques and further computations. The approach in geometric targeting is to gradually change the system's RMSD relative to the target structure while enforcing a set of geometric constraints. The generated pathways are not minimum free energy pathways, but they are geometrically plausible pathways that maintain good covalent bond distances and angles, keep backbone dihedral angles in allowed Ramachandran regions, avoid eclipsed side‐chain torsion angles, avoid non‐bonded overlap, and maintain a set of hydrogen bonds and hydrophobic contacts. Resulting pathways for over 20 proteins featuring a wide variety of conformational changes are reported here, including the very large GroEL complex. Proteins 2010. © 2010 Wiley‐Liss, Inc.     

Exploration of the activation pathway of Δα-Chymotrypsin with molecular dynamics simulations and correlation with kinetic experiments

Correlating the experimentally observed kinetics of protein conformational changes with theoretical predictions is a formidable and challenging task, due to the multitude of degrees of freedom (>5,000) in a protein and the huge gap between the timescale of the kinetic event of interest (ms) and the typical timescale of computer simulations (ns). In this study we show that using the targeted molecular dynamics (TMD) method it is possible to simulate conformational changes of the ms time range and to correlate multiple simulations of single pathways with ensemble experiments on both the structural and energetic basis. As a model system we chose to study the conformational change of rat-Δα-chymotrypsin from its inactive to its active conformation. This activation process has been analyzed previously by experimental and theoretical methods, i.e. fluorescence stopped-flow spectroscopy (FSF), molecular dynamics (MD) and TMD. Inspired by the results of these studies on the wild type (WT) enzyme, several mutants were constructed to alter the conformational pathway and studied by FSF measurements. In the present work WT and mutant N18G were subjected to multiple MD and subsequent TMD simulations. We report the existence of two main activation pathways, a feature of chymotrypsin activation that has been hitherto unknown. A method to correlate the energetics of the different pathways calculated by TMD and the kinetic parameters observed by experimental methods such as FSF is presented. Our work is relevant for experimental single molecule studies of enzymes in general. Electronic supplementary material  The online version of this article (doi:) contains supplementary material, which is available to authorized users.     

Conformational pathway for the kissing complex-->extended dimer transition of the SL1 stem-loop from genomic HIV-1 RNA as monitored by targeted molecular dynamics techniques

HIV-1 retroviral genomic RNA dimerization is initiated by loop-loop interactions between the SL1 stem-loops of two identical RNA molecules. The SL1-SL1 unstable resulting kissing complex (KC) then refolds irreversibly into a more stable complex called extended dimer (ED). Although the structures of both types of complex have been determined, very little is known about the conformational pathway corresponding to the transition, owing to the difficulty of observing experimentally intermediate conformations. In this study, we applied targeted molecular dynamics simulation techniques (TMD) to the phosphorus atoms for monitoring this pathway for the backbone, and a two-step strategy was adopted. In a first step, called TMD(-1), the dimer structure was constrained to progressively move away from KC without indicating the direction, until the RMSD from KC reaches 36A. A total of 20 TMD(-1) simulations were performed under different initial conditions and different simulation parameters. For RMSD ranging between 0 and 22A, the whole set of TMD(-1) simulations follows a similar pathway, then divergences are observed. None of the simulations leads to the ED structure. At RMSD=22A, the dimers look like two parallel Us, still linked by the initial loop-loop interaction, but the strands of the stems (the arms of the Us) are positioned in such a manner that they can form intramolecular as well as intermolecular Watson-Crick base-pairs. This family of structure is called UU. In a second step (TMD simulations), 18 structures were picked up along the pathways generated with TMD(-1) and were constrained to move toward ED by decreasing progressively their RMSD from ED. We found that only structures from the UU family are able to easily reach ED-like conformations of the backbones without exhibiting a large constraint energy.     

The conformational transition pathways of ATP-binding cassette transporter BtuCD revealed by targeted molecular dynamics simulation

BtuCD is a member of the ATP-binding cassette transporters in Escherichia coli that imports vitamin B(12) into the cell by utilizing the energy of ATP hydrolysis. Crystal structures of BtuCD and its homologous protein HI1470/1 in various conformational states support the "alternating access" mechanism which proposes the conformational transitions of the substrate translocation pathway at transmembrane domain (TMD) between the outward-facing and inward-facing states. The conformational transition at TMD is assumed to couple with the movement of the cytoplasmic nucleotide-binding domains (NBDs) driven by ATP hydrolysis/binding. In this study, we performed targeted molecular dynamics (MD) simulations to explore the atomic details of the conformational transitions of BtuCD importer. The outward-facing to inward-facing (O→I) transition was found to be initiated by the conformational movement of NBDs. The subsequent reorientation of the substrate translocation pathway at TMD began with the closing of the periplasmic gate, followed by the opening of the cytoplamic gate in the last stage of the conformational transition due to the extensive hydrophobic interactions at this region, consistent with the functional requirement of unidirectional transport of the substrates. The reverse inward-facing to outward-facing (I→O) transition was found to exhibit intrinsic diversity of the conformational transition pathways and significant structural asymmetry, suggesting that the asymmetric crystal structure of BtuCD-F is an intermediate state in this process.     

Discrete restraint-based protein modeling and the Calpha-trace problem

We present a novel de novo method to generate protein models from sparse, discretized restraints on the conformation of the main chain and side chain atoms. We focus on Calpha-trace generation, the problem of constructing an accurate and complete model from approximate knowledge of the positions of the Calpha atoms and, in some cases, the side chain centroids. Spatial restraints on the Calpha atoms and side chain centroids are supplemented by constraints on main chain geometry, phi/xi angles, rotameric side chain conformations, and inter-atomic separations derived from analyses of known protein structures. A novel conformational search algorithm, combining features of tree-search and genetic algorithms, generates models consistent with these restraints by propensity-weighted dihedral angle sampling. Models with ideal geometry, good phi/xi angles, and no inter-atomic overlaps are produced with 0.8 A main chain and, with side chain centroid restraints, 1.0 A all-atom root-mean-square deviation (RMSD) from the crystal structure over a diverse set of target proteins. The mean model derived from 50 independently generated models is closer to the crystal structure than any individual model, with 0.5 A main chain RMSD under only Calpha restraints and 0.7 A all-atom RMSD under both Calpha and centroid restraints. The method is insensitive to randomly distributed errors of up to 4 A in the Calpha restraints. The conformational search algorithm is efficient, with computational cost increasing linearly with protein size. Issues relating to decoy set generation, experimental structure determination, efficiency of conformational sampling, and homology modeling are discussed.     

On the roles of substrate binding and hinge unfolding in conformational changes of adenylate kinase

We characterized the conformational change of adenylate kinase (AK) between open and closed forms by conducting five all-atom molecular-dynamics simulations, each of 100 ns duration. Different initial structures and substrate binding configurations were used to probe the pathways of AK conformational change in explicit solvent, and no bias potential was applied. A complete closed-to-open and a partial open-to-closed transition were observed, demonstrating the direct impact of substrate-mediated interactions on shifting protein conformation. The sampled configurations suggest two possible pathways for connecting the open and closed structures of AK, affirming the prediction made based on available x-ray structures and earlier works of coarse-grained modeling. The trajectories of the all-atom molecular-dynamics simulations revealed the complexity of protein dynamics and the coupling between different domains during conformational change. Calculations of solvent density and density fluctuations surrounding AK did not show prominent variation during the transition between closed and open forms. Finally, we characterized the effects of local unfolding of an important hinge near Pro¹⁷⁷ on the closed-to-open transition of AK and identified a novel mechanism by which hinge unfolding modulates protein conformational change. The local unfolding of Pro¹⁷⁷ hinge induces alternative tertiary contacts that stabilize the closed structure and prevent the opening transition.     

Building proteins from C alpha coordinates using the dihedral probability grid Monte Carlo method.

Dihedral probability grid Monte Carlo (DPG-MC) is a general-purpose method of conformational sampling that can be applied to many problems in peptide and protein modeling. Here we present the DPG-MC method and apply it to predicting complete protein structures from C alpha coordinates. This is useful in such endeavors as homology modeling, protein structure prediction from lattice simulations, or fitting protein structures to X-ray crystallographic data. It also serves as an example of how DPG-MC can be applied to systems with geometric constraints. The conformational propensities for individual residues are used to guide conformational searches as the protein is built from the amino-terminus to the carboxyl-terminus. Results for a number of proteins show that both the backbone and side chain can be accurately modeled using DPG-MC. Backbone atoms are generally predicted with RMS errors of about 0.5 A (compared to X-ray crystal structure coordinates) and all atoms are predicted to an RMS error of 1.7 A or better.     

Conformational analysis of GT1B ganglioside and its interaction with botulinum neurotoxin type B: a study by molecular modeling and molecular dynamics

The conformational property of oligosaccharide GT1B in aqueous environment was studied by molecular dynamics (MD) simulation using all-atom model. Based on the trajectory analysis, three prominent conformational models were proposed for GT1B. Direct and water-mediated hydrogen bonding interactions stabilize these structures. The molecular modeling and 15 ns MD simulation of the Botulinum Neuro Toxin/B (BoNT/B) - GT1B complex revealed that BoNT/B can accommodate the GT1B in the single binding mode. Least mobility was seen for oligo-GT1B in the binding pocket. The bound conformation of GT1B obtained from the MD simulation of the BoNT/B-GT1B complex bear a close conformational similarity with the crystal structure of BoNT/A-GT1B complex. The mobility noticed for Arg 1268 in the dynamics was accounted for its favorable interaction with terminal NeuNAc. The internal NeuNAc1 tends to form 10 hydrogen bonds with BoNT/B, hence specifying this particular site as a crucial space for the therapeutic design that can restrict the pathogenic activity of BoNT/B.     

Targeted molecular dynamics: A new approach for searching pathways of conformational transitions

Molecular dynamics simulations have proven to be a valuable tool to investigate the dynamic behavior of stable macromolecules at finite temperatures. However, considerable conformational transitions take place during a simulation only accidentally or at exceptionally high temperatures far from the range of experimental conditions. Targeted molecular dynamics (TMD) is a method to induce a conformational change to a known target structure at ordinary temperature by applying a time-dependent, purely geometrical constraint. The transition is enforced independently of the height of energy barriers, while the dynamics of the molecule is only minimally influenced by the constraint. Simulations of decaalanine and insulin show the ability of the method to explore the configurational space for pathways accessible at a given temperature. The transitions studied at insulin comprise unfolding of an α-helical portion and, in the reverse direction, refolding from an extended conformation. A possible application of TMD is the search for energy barriers and stable intermediates from rather local changes up to protein denaturation.     

A probabilistic fragment-based protein structure prediction algorithm

Conformational sampling is one of the bottlenecks in fragment-based protein structure prediction approaches. They generally start with a coarse-grained optimization where mainchain atoms and centroids of side chains are considered, followed by a fine-grained optimization with an all-atom representation of proteins. It is during this coarse-grained phase that fragment-based methods sample intensely the conformational space. If the native-like region is sampled more, the accuracy of the final all-atom predictions may be improved accordingly. In this work we present EdaFold, a new method for fragment-based protein structure prediction based on an Estimation of Distribution Algorithm. Fragment-based approaches build protein models by assembling short fragments from known protein structures. Whereas the probability mass functions over the fragment libraries are uniform in the usual case, we propose an algorithm that learns from previously generated decoys and steers the search toward native-like regions. A comparison with Rosetta AbInitio protocol shows that EdaFold is able to generate models with lower energies and to enhance the percentage of near-native coarse-grained decoys on a benchmark of [Formula: see text] proteins. The best coarse-grained models produced by both methods were refined into all-atom models and used in molecular replacement. All atom decoys produced out of EdaFold's decoy set reach high enough accuracy to solve the crystallographic phase problem by molecular replacement for some test proteins. EdaFold showed a higher success rate in molecular replacement when compared to Rosetta. Our study suggests that improving low resolution coarse-grained decoys allows computational methods to avoid subsequent sampling issues during all-atom refinement and to produce better all-atom models. EdaFold can be downloaded from http://www.riken.jp/zhangiru/software/.     

Solution structure and backbone dynamics of the holo form of the frenolicin acyl carrier protein

During polyketide biosynthesis, acyl carrier proteins (ACPs) perform the central role of transferring polyketide intermediates between active sites of polyketide synthase. The 4'-phosphopantetheine prosthetic group of a holo-ACP is a long and flexible arm that can reach into different active sites and provide a terminal sulfhydryl group for the attachment of acyl groups through a thioester linkage. We have determined the solution structure and characterized backbone dynamics of the holo form of the frenolicin acyl carrier protein (fren holo-ACP) by nuclear magnetic resonance (NMR). Unambiguous assignments were made for 433 hydrogen atoms, 333 carbon atoms, and 84 nitrogen atoms, representing a total of 94.6% of the assignable atoms in this protein. From 879 meaningful NOEs and 45 angle constraints, a family of 24 structures has been calculated. The solution structure is composed of three major alpha-helices packed in a bundle with three additional short helices in intervening loops; one of the short helices slowly exchanges between two conformations. Superposition of the major helical regions on the mean structure yields average atomic rmsd values of 0.49 +/- 0.09 and 0.91 +/- 0.08 A for backbone and non-hydrogen atoms, respectively. Although the three-helix bundle fold is conserved among acyl carrier proteins involved in fatty acid synthases and polyketide synthases, a detailed comparison revealed that ACPs from polyketide biosynthetic pathways are more related to each other in tertiary fold than to their homologues from fatty acid biosynthetic pathways. Comparison of the free form of ACPs (NMR structures of fren ACP and the Bacillus subtilis ACP) with the substrate-bound form of ACP (crystal structure of butyryl-ACP from Escherichia coli) suggests that conformational exchange plays a role in substrate binding.     

Extending the capabilities of targeted molecular dynamics: simulation of a large conformational transition in plasminogen activator inhibitor 1

Plasminogen activator inhibitor type 1 (PAI-1) is an inhibitor of plasminogen activators such as tissue-type plasminogen activator or urokinase-type plasminogen activator. For this molecule, different conformations are known. The inhibiting form that interacts with the proteinases is called the active form. The noninhibitory, noncleavable form is called the latent form. X-ray and modeling studies have revealed a large change in position of the reactive center loop (RCL), responsible for the interaction with the proteinases, that is inserted into a beta-sheet (s4A) in the latent form. The mechanism underlying this spontaneous conformational change (half-life = 2 h at 37 degrees C) is not known in detail. This investigation attempts to predict a transition path from the active to the latent structure at the atomic level, by using simulation techniques. Together with targeted molecular dynamics (TMD), a plausible assumption on a rigid body movement of the RCL was applied to define an initial guess for an intermediate. Different pathways were simulated, from the active to the intermediate, from the intermediate to the latent structure and vice versa under different conditions. Equilibrium simulations at different steps of the path also were performed. The results show that a continuous pathway from the active to the latent structure can be modeled. This study also shows that this approach may be applied in general to model large conformational changes in any kind of protein for which the initial and final three-dimensional structure is known.     

From helix-coil transitions to protein folding

An evolution of procedures to simulate protein structure and folding pathways is described. From an initial focus on the helix-coil transition and on hydrogen-bonding and hydrophobic interactions, our original attempts to determine protein structure and folding pathways were based on an experimental approach. Experiments on the oxidative folding of reduced bovine pancreatic ribonuclease A (RNase A) led to a mechanism by which the molecule folded to the native structure by a minimum of four different pathways. The experiments with RNase A were followed by development of a molecular mechanics approach, first, making use of global optimization procedures and then with molecular dynamics (MD), evolving from an all-atom to a united-residue model. This hierarchical MD approach facilitated probing of the folding trajectory to longer time scales than with all-atom MD, and hence led to the determination of complete folding trajectories, thus far for a protein containing as many as 75 amino acid residues. With increasing refinement of the computational procedures, the computed results are coming closer to experimental observations, providing an understanding as to how physics directs the folding process.     

A molecular model for proflavine-DNA intercalation.

A molecular model has been derived for the intercalation of proflavine into the CpG site of the decamer duplex of d(GATACGATAC). The starting geometry of the intercalation site was taken from previous crystallographic studies on the d(CpG)-proflavine complex, and molecular mechanics used to obtain a stereochemically acceptable structure. This has widened grooves compared to standard A- or B- double helices, as well as distinct conformational, roll, twist and tilt features.     

Flexible fitting of high-resolution x-ray structures into cryoelectron microscopy maps using biased molecular dynamics simulations

A methodology for flexible fitting of all-atom high-resolution structures into low-resolution cryoelectron microscopy (cryo-EM) maps is presented. Flexibility of the modeled structure is simulated by classical molecular dynamics and an additional effective potential is introduced to enhance the fitting process. The additional potential is proportional to the correlation coefficient between the experimental cryo-EM map and a synthetic map generated for an all-atom structure being fitted to the map. The additional forces are calculated as a gradient of the correlation coefficient. During the molecular dynamics simulations under the additional forces, the molecule undergoes a conformational transition that maximizes the correlation coefficient, which results in a high-accuracy fit of all-atom structure into a cryo-EM map. Using five test proteins that exhibit structural rearrangement during their biological activity, we demonstrate performance of our method. We also test our method on the experimental cryo-EM of elongation factor G and show that the model obtained is comparable to previous studies. In addition, we show that overfitting can be avoided by assessing the quality of the fitted model in terms of correlation coefficient and secondary structure preservation.     

Exploring the Conformational Transitions of Biomolecular Systems Using a Simple Two-State Anisotropic Network Model

Biomolecular conformational transitions are essential to biological functions. Most experimental methods report on the long-lived functional states of biomolecules, but information about the transition pathways between these stable states is generally scarce. Such transitions involve short-lived conformational states that are difficult to detect experimentally. For this reason, computational methods are needed to produce plausible hypothetical transition pathways that can then be probed experimentally. Here we propose a simple and computationally efficient method, called ANMPathway, for constructing a physically reasonable pathway between two endpoints of a conformational transition. We adopt a coarse-grained representation of the protein and construct a two-state potential by combining two elastic network models (ENMs) representative of the experimental structures resolved for the endpoints. The two-state potential has a cusp hypersurface in the configuration space where the energies from both the ENMs are equal. We first search for the minimum energy structure on the cusp hypersurface and then treat it as the transition state. The continuous pathway is subsequently constructed by following the steepest descent energy minimization trajectories starting from the transition state on each side of the cusp hypersurface. Application to several systems of broad biological interest such as adenylate kinase, ATP-driven calcium pump SERCA, leucine transporter and glutamate transporter shows that ANMPathway yields results in good agreement with those from other similar methods and with data obtained from all-atom molecular dynamics simulations, in support of the utility of this simple and efficient approach. Notably the method provides experimentally testable predictions, including the formation of non-native contacts during the transition which we were able to detect in two of the systems we studied. An open-access web server has been created to deliver ANMPathway results.     

Dislocation of a type I membrane protein requires interactions between membrane-spanning segments within the lipid bilayer

The human cytomegalovirus gene product US11 causes rapid degradation of class I major histocompatibility complex (MHCI) heavy chains by inducing their dislocation from the endoplasmic reticulum (ER) and subsequent degradation by the proteasome. This set of reactions resembles the endogenous cellular quality control pathway that removes misfolded or unassembled proteins from the ER. We show that the transmembrane domain (TMD) of US11 is essential for MHCI heavy chain dislocation, but dispensable for MHCI binding. A Gln residue at position 192 in the US11 TMD is crucial for the ubiquitination and degradation of MHCI heavy chains. Cells that express US11 TMD mutants allow formation of MHCI-beta2m complexes, but their rate of egress from the ER is significantly impaired. Further mutagenesis data are consistent with the presence of an alpha-helical structure in the US11 TMD essential for MHCI heavy chain dislocation. The failure of US11 TMD mutants to catalyze dislocation is a unique instance in which a polar residue in the TMD of a type I membrane protein is required for that protein's function. Targeting of MHCI heavy chains for dislocation by US11 thus requires the formation of interhelical hydrogen bonds within the ER membrane.     

Molecular dynamics study of the M412 intermediate of bacteriorhodopsin.

Molecular dynamics simulations have been carried out to study the M412 intermediate of bacteriorhodopsin's (bR) photocycle. The simulations start from two simulated structures for the L550 intermediate of the photocycle, one involving a 13-cis retinal with strong torsions, the other a 13,14-dicis retinal, from which the M412 intermediate is initiated through proton transfer to Asp-85. The simulations are based on a refined structure of bR568 obtained through all-atom molecular dynamics simulations and placement of 16 waters inside the protein. The structures of the L550 intermediates were obtained through simulated photoisomerization and subsequent molecular dynamics, and simulated annealing. Our simulations reveal that the M412 intermediate actually comprises a series of conformations involving 1) a motion of retinal; 2) protein conformational changes; and 3) diffusion and reconfiguration of water in the space between the retinal Schiff base nitrogen and the Asp-96 side group. (1) turns the retinal Schiff base nitrogen from an early orientation toward Asp-85 to a late orientation toward Asp-96; (2) disconnects the hydrogen bond network between retinal and Asp-85 and tilts the helix F of bR, enlarging bR's cytoplasmic channel; (3) adds two water molecules to the three water molecules existing in the cytoplasmic channel at the bR568 stage and forms a proton conduction pathway. The conformational change (2) of the protein involves a 60 degrees bent of the cytoplasmic side of helix F and is induced through a break of a hydrogen bond between Tyr-185 and a water-side group complex in the counterion region.     

Comparison of Generalised Born/Surface Area with Periodic Boundary Simulations to Study Protein Unfolding

Simulations comparing the rapid unfolding behaviour of the model protein barnase under explicit and implicit solvent systems have been undertaken in order to validate a faster implicit solvent method for studying proteins which are kinetically stable in silico. A comparison is made between all-atom explicitly solvated simulations of barnase undertaken using Particle Mesh Ewald electrostatic interactions with all-atom implicit solvent simulations undertaken using the generalised born/surface area (GBSA) method with a long non-bonded cut-off. The two explicit solvent unfolding trajectories appear to explore slightly different pathways showing the importance of having statistically valid ensembles which are not accessible from a single trajectory. The 500?K GBSA trajectory is unsuitable for exploring intermediate structures on the unfolding pathway of barnase, as the protein almost immediately jumps to a predominately random coil conformation. However, dropping the temperature to 400?K gives rise to trajectories where the protein is unable to climb out of the energy well containing the first intermediate state, in a reasonable timescale. A similar pattern to the explicit solvent unfolding trajectories is seen in 450?K GBSA runs, with the intermediate states differing between trajectories. The development of computer simulation methods suitable for application to more kinetically stable proteins will offer insight into the atomic detail of the conformational changes associated with protein unfolding diseases.     

Transition state of a SH3 domain detected with principle component analysis and a charge-neutralized all-atom protein model

The src SH3 domain has been known to be a two-state folder near room temperature. However, in a previous study with an all-atom model simulation near room temperature, the transition state of this protein was not successfully detected on a free-energy profile using two axes: the radius of gyration (R(g)) and native contact reproduction ratio (Q value). In this study, we focused on an atom packing effect to characterize the transition state and tried another analysis to detect it. To explore the atom packing effect more efficiently, we introduced a charge-neutralized all-atom model, where all of the atoms in the protein and water molecules were treated explicitly, but their partial atomic charges were set to zero. Ten molecular dynamics simulations were performed starting from the native structure at 300 K, where the simulation length of each run was 90 ns, and the protein unfolded in all runs. The integrated trajectories (10 x 90 = 900 ns) were analyzed by a principal component analysis (PCA) and showed a clear free-energy barrier between folded- and unfolded-state conformational clusters in a conformational space generated by PCA. There were segments that largely deformed when the conformation passed through the free-energy barrier. These segments correlated well with the structural core regions characterized by large phi-values, and the atom-packing changes correlated with the conformational deformations. Interestingly, using the same simulation data, no significant barrier was found in a free-energy profile using the R(g) and Q values for the coordinate axes. These results suggest that the atom packing effect may be one of the most important determinants of the transition state.