Free energy landscapes of two model peptides: α-helical and β-hairpin peptides explored with Brownian dynamics simulation

We applied an atomistic Brownian dynamics (BD) simulation with multiple time step method for the folding simulation of a 13-mer α-helical peptide and a 12-mer β-hairpin peptide, giving successful folding simulations. In this model, the driving energy contribution towards folding came from both electrostatic and van der Waals interactions for the α-helical peptide and from van der Waals interactions for the β-hairpin peptide. Although, many non-native structures having the same or lower energy than that of native structure were observed, the folded states formed the most populated cluster when the structures obtained by the BD simulations were subjected to the cluster analysis based on distance-based root mean square deviation of side-chains between different structures. This result indicates that we can predict the native structures from conformations sampled by BD simulation.     

Multiple time step diffusive Langevin dynamics for proteins

We present an algorithm for simulating the long time scale dynamics of proteins and other macromolecules. Our method applies the concept of multiple time step integration to the diffusive Langevin equation, in which short time scale dynamics are replaced by friction and noise. The macromolecular force field is represented at atomic resolution. Slow motions are modeled by constrained Langevin dynamics with very large time steps, while faster degrees of freedom are kept in local thermal equilibrium. In the limit of a sufficiently large molecule, our algorithm is shown to reduce the CPU time required by two orders of magnitude. We test the algorithm on two systems, alanine dipeptide and bovine pancreatic trypsin inhibitor (BPTI), and find that it accurately calculates a variety of equilibrium and dynamical properties. In the case of BPTI, the CPU time required is reduced by nearly a factor of 60 compared to a conventional, unconstrained Langevin simulation using the same force field. Proteins 30:215–227, 1998. © 1998 Wiley-Liss, Inc.     

Brownian dynamics simulation of protein folding: a study of the diffusion-collision model

The dynamic aspects of protein folding are described by a series of diffusion-collision steps involving structural units (microdomains) of various sizes that combine to form the protein in its native state. A method is introduced for obtaining the rate constants for the basic diffusion-collision step by use of Brownian dynamics. The method is applied to an investigation of the folding dynamics of two α-helices connected by a flexible (random-coil) polypeptide chain. The results of this full three-dimensional treatment are compared with simplified model calculations for the diffusion-collision step. Of particular interest are the nature of the collision dynamics and the role of the intervening peptide chain.     

Molecular dynamics simulations of peptides and proteins with amplified collective motions

We present a novel method that uses the collective modes obtained with a coarse-grained model/anisotropic network model to guide the atomic-level simulations. Based on this model, local collective modes can be calculated according to a single configuration in the conformational space of the protein. In the molecular dynamics simulations, the motions along the slowest few modes are coupled to a higher temperature by the weak coupling method to amplify the collective motions. This amplified-collective-motion (ACM) method is applied to two test systems. One is an S-peptide analog. We realized the refolding of the denatured peptide in eight simulations out of 10 using the method. The other system is bacteriophage T4 lysozyme. Much more extensive domain motions between the N-terminal and C-terminal domain of T4 lysozyme are observed in the ACM simulation compared to a conventional simulation. The ACM method allows for extensive sampling in conformational space while still restricting the sampled configurations within low energy areas. The method can be applied in both explicit and implicit solvent simulations, and may be further applied to important biological problems, such as long timescale functional motions, protein folding/unfolding, and structure prediction.     

Role of hydrophobic interactions and salt-bridges in β-hairpin folding

β-Hairpins are the simplest form of β-sheets which, due to the presence of long-range interactions, can be considered as tertiary structures. Molecular dynamics simulation is a powerful tool that can unravel whole pathways of protein folding/unfolding at atomic resolution. We have performed several molecular dynamics simulations, to a total of over 250 ns, of a β-hairpin peptide in water using GROMACS. We show that hydrophobic interactions are necessary for initiating the folding of the peptide. Once formed, the peptide is stabilized by hydrogen bonds and disruption of hydrophobic interactions in the folded peptide does not denature the structure. In the absence of hydrophobic interactions, the peptide fails to fold. However, the introduction of a salt-bridge compensates for the loss of hydrophobic interactions to a certain extent. Figure Model of b-hairpin folding: Folding is initiated by hydrophobic interactions (Brown circles). The folded structure, once formed, is stabilized by hydrogen bonds (red lines) and is unaffected by loss of hydrophobic contacts     

Outer membrane proteins: comparing X-ray and NMR structures by MD simulations in lipid bilayers

The structures of three bacterial outer membrane proteins (OmpA, OmpX and PagP) have been determined by both X-ray diffraction and NMR. We have used multiple (7 × 15 ns) MD simulations to compare the conformational dynamics resulting from the X-ray versus the NMR structures, each protein being simulated in a lipid (DMPC) bilayer. Conformational drift was assessed via calculation of the root mean square deviation as a function of time. On this basis the ‘quality’ of the starting structure seems mainly to influence the simulation stability of the transmembrane β-barrel domain. Root mean square fluctuations were used to compare simulation mobility as a function of residue number. The resultant residue mobility profiles were qualitatively similar for the corresponding X-ray and NMR structure-based simulations. However, all three proteins were generally more mobile in the NMR-based than in the X-ray simulations. Principal components analysis was used to identify the dominant motions within each simulation. The first two eigenvectors (which account for >50% of the protein motion) reveal that such motions are concentrated in the extracellular loops and, in the case of PagP, in the N-terminal α-helix. Residue profiles of the magnitude of motions corresponding to the first two eigenvectors are similar for the corresponding X-ray and NMR simulations, but the directions of these motions correlate poorly reflecting incomplete sampling on a ∼10 ns timescale.     

Backbone entropy of loops as a measure of their flexibility: Application to a Ras protein simulated by molecular dynamics

The flexibility of surface loops plays an important role in protein–protein and protein–peptide recognition; it is commonly studied by Molecular Dynamics or Monte Carlo simulations. We propose to measure the relative backbone flexibility of loops by the difference in their backbone conformational entropies, which are calculated here with the local states (LS) method of Meirovitch. Thus, one can compare the entropies of loops of the same protein or, under certain simulation conditions, of different proteins. These loops should be equal in size but can differ in their sequence of amino acids residues. This methodology is applied successfully to three segments of 10 residues of a Ras protein simulated by the stochastic boundary molecular dynamics procedure. For the first time estimates of backbone entropy differences are obtained, and their correlation with B factors is pointed out; for example, the segments which consist of residues 60–65 and 112–117 have average B factors of 67 and 18 Ų, respectively, and entropy difference T ΔS = 5.4 ± 0.1 kcal/mol at T = 300 K. In a large number of recent publications the entropy due to the fast motions (on the ps-ns time scale) of N–H and C–H vectors has been obtained from their order parameter, measured in nuclear magnetic resonance spin relaxation experiments. This enables one to estimate differences in the entropy of protein segments due to folding–unfolding transitions, for example. However, the vectors are assumed to be independent, and the effect of the neglected correlations is unknown; our method is expected to become an important tool for assessing this approximation. The present calculations, obtained with the LS method, suggest that the errors involved in experimental entropy differences might not be large; however, this should be verified in each case. Potential applications of entropy calculations to rational drug design are discussed. Proteins 29:127–140, 1997. © 1997 Wiley-Liss, Inc.     

Calcium binding and allosteric signaling mechanisms for the sarcoplasmic reticulum Ca(2+) ATPase

The sarcoplasmic reticulum Ca(2+) ATPase (SERCA) is a membrane-bound pump that utilizes ATP to drive calcium ions from the myocyte cytosol against the higher calcium concentration in the sarcoplasmic reticulum. Conformational transitions associated with Ca(2+) -binding are important to its catalytic function. We have identified collective motions that partition SERCA crystallographic structures into multiple catalytically-distinct states using principal component analysis. Using Brownian dynamics simulations, we demonstrate the important contribution of surface-exposed, polar residues in the diffusional encounter of Ca(2+) . Molecular dynamics simulations indicate the role of Glu309 gating in binding Ca(2+) , as well as subsequent changes in the dynamics of SERCA's cytosolic domains. Together these data provide structural and dynamical insights into a multistep process involving Ca(2+) binding and catalytic transitions.     

Predicting Peptide Structures in Native Proteins from Physical Simulations of Fragments

It has long been proposed that much of the information encoding how a protein folds is contained locally in the peptide chain. Here we present a large-scale simulation study designed to examine the extent to which conformations of peptide fragments in water predict native conformations in proteins. We perform replica exchange molecular dynamics (REMD) simulations of 872 8-mer, 12-mer, and 16-mer peptide fragments from 13 proteins using the AMBER 96 force field and the OBC implicit solvent model. To analyze the simulations, we compute various contact-based metrics, such as contact probability, and then apply Bayesian classifier methods to infer which metastable contacts are likely to be native vs. non-native. We find that a simple measure, the observed contact probability, is largely more predictive of a peptide''s native structure in the protein than combinations of metrics or multi-body components. Our best classification model is a logistic regression model that can achieve up to 63% correct classifications for 8-mers, 71% for 12-mers, and 76% for 16-mers. We validate these results on fragments of a protein outside our training set. We conclude that local structure provides information to solve some but not all of the conformational search problem. These results help improve our understanding of folding mechanisms, and have implications for improving physics-based conformational sampling and structure prediction using all-atom molecular simulations.     

Insights into nucleic acid conformational dynamics from massively parallel stochastic simulations

The helical hairpin is one of the most ubiquitous and elementary secondary structural motifs in nucleic acids, capable of serving functional roles and participating in long-range tertiary contacts. Yet the self-assembly of these structures has not been well-characterized at the atomic level. With this in mind, the dynamics of nucleic acid hairpin formation and disruption have been studied using a novel computational tool: large-scale, parallel, atomistic molecular dynamics simulation employing an inhomogeneous distributed computer consisting of more than 40,000 processors. Using multiple methodologies, over 500 micro s of atomistic simulation time has been collected for a large ensemble of hairpins (sequence 5'-GGGC[GCAA]GCCU-3'), allowing characterization of rare events not previously observable in simulation. From uncoupled ensemble dynamics simulations in unperturbed folding conditions, we report on 1), competing pathways between the folded and unfolded regions of the conformational space; 2), observed nonnative stacking and basepairing traps; and 3), a helix unwinding-rewinding mode that is differentiated from the unfolding and folding dynamics. A heterogeneous transition state ensemble is characterized structurally through calculations of conformer-specific folding probabilities and a multiplexed replica exchange stochastic dynamics algorithm is used to derive an approximate folding landscape. A comparison between the observed folding mechanism and that of a peptide beta-hairpin analog suggests that although native topology defines the character of the folding landscape, the statistical weighting of potential folding pathways is determined by the chemical nature of the polymer.     

Determining Intrachain Diffusion Coefficients for Biopolymer Dynamics from Single-Molecule Force Spectroscopy Measurements

The conformational diffusion coefficient for intrachain motions in biopolymers, D, sets the timescale for structural dynamics. Recently, force spectroscopy has been applied to determine D both for unfolded proteins and for the folding transitions in proteins and nucleic acids. However, interpretation of the results remains unsettled. We investigated how instrumental effects arising from the force probes used in the measurement can affect the value of D recovered via force spectroscopy. We compared estimates of D for the folding of DNA hairpins found from measurements of rates and energy landscapes made using optical tweezers with estimates obtained from the same single-molecule trajectories via the transition path time. The apparent D obtained from the rates was much lower than the result found from the same data using transition time analysis, reflecting the effects of the mechanical properties of the force probe. Deconvolution of the finite compliance effects on the measurement allowed the intrinsic value to be recovered. These results were supported by Brownian dynamics simulations of the effects of force-probe compliance and bead size.     

Determining Intrachain Diffusion Coefficients for Biopolymer Dynamics from Single-Molecule Force Spectroscopy Measurements

The conformational diffusion coefficient for intrachain motions in biopolymers, D, sets the timescale for structural dynamics. Recently, force spectroscopy has been applied to determine D both for unfolded proteins and for the folding transitions in proteins and nucleic acids. However, interpretation of the results remains unsettled. We investigated how instrumental effects arising from the force probes used in the measurement can affect the value of D recovered via force spectroscopy. We compared estimates of D for the folding of DNA hairpins found from measurements of rates and energy landscapes made using optical tweezers with estimates obtained from the same single-molecule trajectories via the transition path time. The apparent D obtained from the rates was much lower than the result found from the same data using transition time analysis, reflecting the effects of the mechanical properties of the force probe. Deconvolution of the finite compliance effects on the measurement allowed the intrinsic value to be recovered. These results were supported by Brownian dynamics simulations of the effects of force-probe compliance and bead size.     

A directed essential dynamics simulation of peptide folding

We present a directed essential dynamics (DED) method for peptide and protein folding. DED is a molecular dynamics method based on the essential dynamics sampling and the principal component analysis. The main idea of DED is to use principal component analysis to determine the direction of the most active collective motion of peptides at short intervals of time (20 fs) during the folding process and then add an additional force along it to adjust the folding direction. This method can make the peptides avoid being trapped in the local minima for a long time and enhance the sampling efficiency in conformational space during the simulation. An S-peptide with 15 amino acids is used to demonstrate the DED method. The results show that DED can lead the S-peptide to fold quickly into the native state, whereas traditional molecular dynamics needs more time to do this.     

The triplex-hairpin transition in cytosine-rich DNA

We present a theoretical study of the self-complementary single-stranded 30-mer d(TC*TTC*C*TTTTCCTTCTC*CCGAGAAGGTTTT) (PDB ID: 1b4y) that was designed to form an intramolecular triplex by folding back twice on itself. At neutral pH the molecule exists in a duplex hairpin conformation, whereas at acidic pH the cytosines labeled by an asterisk (*) are protonated, forming Hoogsteen hydrogen bonds with guanine of a GC Watson-Crick basepair to generate a triplex. As a first step in an investigation of the energetics of the triplex-hairpin transition, we applied the Bashford-Karplus multiple site model of protonation to calculate the titration curves for the two conformations. Based on these data, a two-state model is used to study the equilibrium properties of transition. Although this model properly describes the thermodynamics of the protonation-deprotonation steps that drive the folding-unfolding of the oligomer, it cannot provide insight into the time-dependent mechanism of the process. A series of molecular dynamics simulations using the ff94 force field of the AMBER 6.0 package was therefore run to explore the dynamics of the folding/unfolding pathway. The molecular dynamics method was combined with Poisson-Boltzmann calculations to determine when a change in protonation state was warranted during a trajectory. This revealed a sequence of elementary protonation steps during the folding/unfolding transition and suggests a strong coupling between ionization and folding in cytosine-rich triple-helical triplexes.     

Weighted-ensemble Brownian dynamics simulations for protein association reactions.

A new method, weighted-ensemble Brownian dynamics, is proposed for the simulation of protein-association reactions and other events whose frequencies of outcomes are constricted by free energy barriers. The method features a weighted ensemble of trajectories in configuration space with energy levels dictating the proper correspondence between "particles" and probability. Instead of waiting a very long time for an unlikely event to occur, the probability packets are split, and small packets of probability are allowed to diffuse almost immediately into regions of configuration space that are less likely to be sampled. The method has been applied to the Northrup and Erickson (1992) model of docking-type diffusion-limited reactions and yields reaction rate constants in agreement with those obtained by direct Brownian simulation, but at a fraction of the CPU time (10(-4) to 10(-3), depending on the model). Because the method is essentially a variant of standard Brownian dynamics algorithms, it is anticipated that weighted-ensemble Brownian dynamics, in conjunction with biophysical force models, can be applied to a large class of association reactions of interest to the biophysics community.     

Protein folding stabilizing time measurement: a direct folding process and three-dimensional random walk simulation

Protein particles undergo Brownian motion and collisions in solution. The diffusive collisions may lead to aggregation. For proteins to fold successfully the process has to occur quickly and before significant collision takes place. The speed of protein folding was deduced by studying the correlation time of a lysozyme refolding process from autocorrelation function analysis of the mean collision time and aggregation/soluble ratio of protein. It is a measure of time before which an aggregate can be formed and also is the time measure for a protein to fold into a stable state. We report on the protein folding stabilizing time of a lysozyme system to be 25.5-27.5 micros (<+/-4%) between 295 and 279K via direct folding experimental studies, supported by a three-dimensional random walk simulation of diffusion-limited aggregation model. Aggregation is suppressed when the protein is folded to a stable form. Spontaneous folding and diffusion-limited aggregation are antagonistic in nature. Meanwhile, the resultant aggresome, suggested by Raman and mass spectroscopy, may be formed by cross-linkages of disulfide bonds and hydrophobic interactions.     

Systematically Constructing Kinetic Transition Network in Polypeptide from Top to Down: Trajectory Mapping

Molecular dynamics (MD) simulation is an important tool for understanding bio-molecules in microscopic temporal/spatial scales. Besides the demand in improving simulation techniques to approach experimental scales, it becomes more and more crucial to develop robust methodology for precisely and objectively interpreting massive MD simulation data. In our previous work [J Phys Chem B 114, 10266 (2010)], the trajectory mapping (TM) method was presented to analyze simulation trajectories then to construct a kinetic transition network of metastable states. In this work, we further present a top-down implementation of TM to systematically detect complicate features of conformational space. We first look at longer MD trajectory pieces to get a coarse picture of transition network at larger time scale, and then we gradually cut the trajectory pieces in shorter for more details. A robust clustering algorithm is designed to more effectively identify the metastable states and transition events. We applied this TM method to detect the hierarchical structure in the conformational space of alanine-dodeca-peptide from microsecond to nanosecond time scales. The results show a downhill folding process of the peptide through multiple pathways. Even in this simple system, we found that single common-used order parameter is not sufficient either in distinguishing the metastable states or predicting the transition kinetics among these states.     

Rational design of antithrombotic peptides to target the von Willebrand Factor (vWf) - GPIb integrin interaction

Conventional antithrombotic drug discovery requires testing of large numbers of drug candidates. We used computer-aided macromolecular interaction assessment (MIAX) to select antithrombotic molecules that mimic and therefore block platelet GPIb’s binding to von Willebrand factor (vWf), an early step in thrombus formation. We screened a random array of 15-mer D-amino acid peptides for binding vWf. Structures of 4 candidate peptides were inferred by comparison to sequences in protein databases, conversion from the L to D conformations and molecular dynamics (MD) determinations of those most energetically stable. By MIAX, we deduced the amino acids and intermolecular hydrogen bonds contributing to the GPIb-vWf interaction interface. We docked the peptides onto vWf in silico to localize their binding sites and consequent potential for preventing GPIb-vWf binding. In vitro inhibition of ristocetin-initiated platelet agglutination confirmed peptide function and suitability for antithrombotic development, thereby validating this novel approach to drug discovery.