Ab initio folding simulation of Trpcage by replica exchange with hybrid Hamiltonian

Replica-exchange molecular dynamics simulations with hybrid Hamiltonian in explicit solvent were performed to study the folding of a designed 20-residue miniprotein, Trpcage, from a fully extended structure. During the simulations several folding/unfolding events happened. In the folded states the majority of experimentally observed NMR NOE restraints are satisfied. The folded structures have root mean squared deviation of 2.0 A with respect to the NMR structures considering all heavy atoms. The free-energy surface constructed by the hybrid Hamiltonian simulations is similar to the one built by a standard replica-exchange simulation which started from the native structure. Consistent with previous experimental observation, a pre-existing hydrophobic collapse in the unfolded state is detected by investigating the desolvation behavior of Trpcage. At room temperature, an intermediate state featured by a misfolded core, a nearly formed alpha-helix segment and an absence of 3(10)-helix is found. The replica exchange with hybrid Hamiltonian method is shown here to be capable of resolving the folding picture of the miniprotein.     

Beta-hairpin folding mechanism of a nine-residue peptide revealed from molecular dynamics simulations in explicit water

The beta-hairpin fold mechanism of a nine-residue peptide, which is modified from the beta-hairpin of alpha-amylase inhibitor tendamistat (residues 15-23), is studied through direct folding simulations in explicit water at native folding conditions. Three 300-nanosecond self-guided molecular dynamics (SGMD) simulations have revealed a series of beta-hairpin folding events. During these simulations, the peptide folds repeatedly into a major cluster of beta-hairpin structures, which agree well with nuclear magnetic resonance experimental observations. This major cluster is found to have the minimum conformational free energy among all sampled conformations. This peptide also folds into many other beta-hairpin structures, which represent some local free energy minimum states. In the unfolded state, the N-terminal residues of the peptide, Tyr-1, Gln-2, and Asn-3, have a confined conformational distribution. This confinement makes beta-hairpin the only energetically favored structure to fold. The unfolded state of this peptide is populated with conformations with non-native intrapeptide interactions. This peptide goes through fully hydrated conformations to eliminate non-native interactions before folding into a beta-hairpin. The folding of a beta-hairpin starts with side-chain interactions, which bring two strands together to form interstrand hydrogen bonds. The unfolding of the beta-hairpin is not simply the reverse of the folding process. Comparing unfolding simulations using MD and SGMD methods demonstrate that SGMD simulations can qualitatively reproduce the kinetics of the peptide system.     

Direct folding simulation of alpha-helices and beta-hairpins based on a single all-atom force field with an implicit solvation model

Recently, we have shown that a modified energy model based on the param99 force field with the generalized Born (GB) solvation model produces reliable free energy landscapes of mini-proteins with a betabetaalpha motif (BBA5, 1FSD, and 1PSV), with the native structures of the mini-proteins located in their lowest free energy minimum states. One of the main features in the modified energy model is a significant improvement for more balanced treatments of alpha and beta strands in proteins. In this study, using the replica exchange molecular dynamics (REMD) simulation method with this new force field, we have carried out extensive ab initio folding studies of several well-known peptides with alpha or beta strands (C-peptide, EK-peptide, le0q, and gbl). Starting from fully extended conformations as the initial conditions, all of the native-like structures of the target peptides were successfully identified by REMD, with reasonable representations of free energy surfaces. The present simulation results with the modified energy model are consistent with experiments, demonstrating an extended applicability of the energy model to folding studies of a variety of alpha-helices, beta-strands, and alpha/beta proteins.     

Stability and amino acid preferences of type VIII reverse turn: the most common turn in peptides?

Free energies of the alpha(r)beta and betabeta conformations of 14 tetrapeptides, based on the sequence SALN and protein X-ray structures, were calculated using molecular dynamics simulations and MM-PBSA calculations. The alphaalpha conformations of five of the tetrapeptides were also studied. SALN has been earlier shown by molecular dynamics simulations and NMR spectroscopy to have a tendency to form an alpha(r)beta turn. The gas-phase energy of the molecular mechanical force field (CHARMM), the electrostatic and non-polar solvation free energies and solute entropies were used to explain the free energy differences of the alphaalpha, betabeta and alpha(r)beta conformations of the peptides. The alpha(r)beta conformation of SALN and SATN was predicted to be slightly more stable than the extended conformation (betabeta), in agreement with experimental results. The SALN mutants SAIN, SAVN, SATN, SSIN and MSHV, were also predicted to be potential alpha(r)beta turn-forming peptides. We report also revised positional potentials for the type VIII turn, based on a non-homologous set of protein structures. This protein databank analysis confirms the main results of the earlier analyses and reveals several new amino acid residues with a significant positional preference. The results of this work led us to suggest that the alpha(r)beta turn may be the most common turn type in peptides. Such turns may be readily formed in aqueous solution and thereby play important roles in the protein folding process by serving as an initiation point for structure formation.     

Topological determinants of protein domain swapping

Protein domain swapping has been repeatedly observed in a variety of proteins and is believed to result from destabilization due to mutations or changes in environment. Based on results from our studies and others, we propose that structures of the domain-swapped proteins are mainly determined by their native topologies. We performed molecular dynamics simulations of seven different proteins, known to undergo domain swapping experimentally, under mildly denaturing conditions and found in all cases that the domain-swapped structures can be recapitulated by using protein topology in a simple protein model. Our studies further indicated that, in many cases, domain swapping occurs at positions around which the protein tends to unfold prior to complete unfolding. This, in turn, enabled prediction of protein structural elements that are responsible for domain swapping. In particular, two distinct domain-swapped dimer conformations of the focal adhesion targeting domain of focal adhesion kinase were predicted computationally and were supported experimentally by data obtained from NMR analyses.     

Correlation between knowledge-based and detailed atomic potentials: application to the unfolding of the GCN4 leucine zipper.

The relationship between the unfolding pseudo free energies of reduced and detailed atomic models of the GCN4 leucine zipper is examined. Starting from the native crystal structure, a large number of conformations ranging from folded to unfolded were generated by all-atom molecular dynamics unfolding simulations in an aqueous environment at elevated temperatures. For the detailed atomic model, the pseudo free energies are obtained by combining the CHARMM all-atom potential with a solvation component from the generalized Born, surface accessibility, GB/SA, model. Reduced model energies were evaluated using a knowledge-based potential. Both energies are highly correlated. In addition, both show a good correlation with the root mean square deviation, RMSD, of the backbone from native. These results suggest that knowledge-based potentials are capable of describing at least some of the properties of the folded as well as the unfolded states of proteins, even though they are derived from a database of native protein structures. Since only conformations generated from an unfolding simulation are used, we cannot assess whether these potentials can discriminate the native conformation from the manifold of alternative, low-energy misfolded states. Nevertheless, these results also have significant implications for the development of a methodology for multiscale modeling of proteins that combines reduced and detailed atomic models.     

Loop propensity of the sequence YKGQP from staphylococcal nuclease: implications for the folding of nuclease.

Recently we performed molecular dynamics (MD) simulations on the folding of the hairpin peptide DTVKLMYKGQPMTFR from staphylococcal nuclease in explicit water. We found that the peptide folds into a hairpin conformation with native and nonnative hydrogen-bonding patterns. In all the folding events observed in the folding of the hairpin peptide, loop formation involving the region YKGQP was an important event. In order to trace the origins of the loop propensity of the sequence YKGQP, we performed MD simulations on the sequence starting from extended, polyproline II and native type I' turn conformations for a total simulation length of 300 ns, using the GROMOS96 force field under constant volume and temperature (NVT) conditions. The free-energy landscape of the peptide YKGQP shows minima corresponding to loop conformation with Tyr and Pro side-chain association, turn and extended conformational forms, with modest free-energy barriers separating the minima. To elucidate the role of Gly in facilitating loop formation, we also performed MD simulations of the mutated peptide YKAQP (Gly --> Ala mutation) under similar conditions starting from polyproline II conformation for 100 ns. Two minima corresponding to bend/turn and extended conformations were observed in the free-energy landscape for the peptide YKAQP. The free-energy barrier between the minima in the free-energy landscape of the peptide YKAQP was also modest. Loop conformation is largely sampled by the YKGQP peptide, while extended conformation is largely sampled by the YKAQP peptide. We also explain why the YKGQP sequence samples type II turn conformation in these simulations, whereas the sequence as part of the hairpin peptide DTVKLMYKGQPMTFR samples type I' turn conformation both in the X-ray crystal structure and in our earlier simulations on the folding of the hairpin peptide. We discuss the implications of our results to the folding of the staphylococcal nuclease.     

The MUMO (minimal under-restraining minimal over-restraining) method for the determination of native state ensembles of proteins

While reliable procedures for determining the conformations of proteins are available, methods for generating ensembles of structures that also reflect their flexibility are much less well established. Here we present a systematic assessment of the ability of ensemble-averaged molecular dynamics simulations with ensemble-averaged NMR restraints to simultaneously reproduce the average structure of proteins and their associated dynamics. We discuss the effects that under-restraining (overfitting) and over-restraining (underfitting) have on the structures generated in ensemble-averaged molecular simulations. We then introduce the MUMO (minimal under-restraining minimal over-restraining) method, a procedure in which different observables are averaged over a different number of molecules. As both over-restraining and under-restraining are significantly reduced in the MUMO method, it is possible to generate ensembles of conformations that accurately characterize both the structure and the dynamics of native states of proteins. The application of the MUMO method to the protein ubiquitin yields a high-resolution structural ensemble with an RDC Q-factor of 0.19.     

A fast method for predicting amino acid mutations that lead to unfolding

Amino acid mutation(s) that cause(s) partial or total unfolding of a protein can lead to disease states and failure to produce mutants. It is therefore very useful to be able to predict which mutations can retain the conformation of a wild-type protein and which mutations will lead to local or global unfolding of the protein. We have developed a fast and reasonably accurate method based on a backbone-dependent side-chain rotamer library to predict the (folded or unfolded) conformation of a protein upon mutation. This method has been tested on proteins whose wild-type 3D structures are known and whose mutant conformations have been experimentally characterized to be folded or unfolded. Furthermore, for the cases studied here, the predicted partially folded or denatured mutant conformation correlate with a decrease in the stability of the mutant relative to the wild-type protein. The key advantage of our method is that it is very fast and predicts locally or globally unfolded states fairly accurately. Hence, it may prove to be useful in designing site-directed mutagenesis, X-ray crystallography and drug design experiments as well as in free energy simulations by helping to ascertain whether a mutation will alter or retain the wild-type conformation.     

Thermodynamics and mechanism of alpha helix initiation in alanine and valine peptides

We used molecular dynamics simulations to study the folding/unfolding of one of turn of an alpha helix in Ac-(Ala)3-NHMe and Ac-(Val)3-NHMe. Using specialized sampling techniques, we computed free energy surfaces as functions of a conformational coordinate that corresponds to alpha helices at small values and to extended conformations at large values. Analysis of the peptide conformations populated during the simulations showed that alpha helices, reverse turns, and extended conformations correspond to minima on the free energy surfaces of both peptides. The free energy difference between alpha helix and extended conformations, determined from the equilibrium constants for helix unfolding, is approximately -1 kcal/mol for Ac-(Ala)3-NHMe and -5 kcal/mol for Ac-(Val)3-NHMe. The mechanism observed in our simulations, which includes reverse turns as important intermediates along the helix folding/unfolding pathway, is consistent with a mechanism proposed previously. Our results predict that both peptides (but especially the Ala peptide) have a much larger equilibrium constant for helix initiation than is predicted by the helix-coil transition theory with the host-guest parameters. We also predict a much greater difference in the equilibrium constants than the theory predicts. Insofar as helix initiation is concerned, our results suggest that the large difference between the helical propensities of Ala and Val cannot be explained by simple concepts such as side-chain rotamer restriction or unfavorable steric interactions. Rather, the origin of the difference appears to be quite complicated because it involves subtle differences in the solvation of the two peptides. The two peptides have similar turn-extended equilibria but very different helix-turn equilibria, and the difference in helical propensities reflects the fact that the helix-turn equilibrium strongly favors the turns in Ac-(Val)3-NHMe, while it favors the helices in Ac-(Ala)3-NHMe. We also computed thermodynamic decompositions of the free energy surfaces, and these revealed that the helix-turn equilibria are vastly different primarily because the changes in peptide-water interactions that accompany helix-to-turn conformational changes are qualitatively different for the two peptides.     

On prediction of folding nuclei in globular proteins

The approach described in this paper on the prediction of folding nuclei in globular proteins with known three dimensional structures is based on a search of the lowest saddle points through the barrier separating the unfolded state from the native structure on the free-energy landscape of protein chain. This search is performed by a dynamic programming method. Comparison of theoretical results with experimental data on the folding nuclei of two dozen of proteins shows that our model provides good phi value predictions for proteins whose structures have been determined by X-ray analysis, with a less limited success for proteins whose structures have been determined by NMR techniques only. Consideration of a full ensemble of transition states results in more successful prediction than consideration of only the transition states with the minimal free energy. In conclusion we have predicted the localization of folding nuclei for three dimensional protein structures for which kinetics of folding is studied now but the localization of folding nuclei is still unknown.     

Design and Structure of an Equilibrium Protein Folding Intermediate: A Hint into Dynamical Regions of Proteins

Partly unfolded protein conformations close to the native state may play important roles in protein function and in protein misfolding. Structural analyses of such conformations which are essential for their fully physicochemical understanding are complicated by their characteristic low populations at equilibrium. We stabilize here with a single mutation the equilibrium intermediate of apoflavodoxin thermal unfolding and determine its solution structure by NMR. It consists of a large native region identical with that observed in the X-ray structure of the wild-type protein plus an unfolded region. Small-angle X-ray scattering analysis indicates that the calculated ensemble of structures is consistent with the actual degree of expansion of the intermediate. The unfolded region encompasses discontinuous sequence segments that cluster in the 3D structure of the native protein forming the FMN cofactor binding loops and the binding site of a variety of partner proteins. Analysis of the apoflavodoxin inner interfaces reveals that those becoming destabilized in the intermediate are more polar than other inner interfaces of the protein. Natively folded proteins contain hydrophobic cores formed by the packing of hydrophobic surfaces, while natively unfolded proteins are rich in polar residues. The structure of the apoflavodoxin thermal intermediate suggests that the regions of natively folded proteins that are easily responsive to thermal activation may contain cores of intermediate hydrophobicity.     

A Free-Energy Approach for All-Atom Protein Simulation

All-atom free-energy methods offer a promising alternative to kinetic molecular mechanics simulations of protein folding and association. Here we report an accurate, transferable all-atom biophysical force field (PFF02) that stabilizes the native conformation of a wide range of proteins as the global optimum of the free-energy landscape. For 32 proteins of the ROSETTA decoy set and six proteins that we have previously folded with PFF01, we find near-native conformations with an average backbone RMSD of 2.14 Å to the native conformation and an average Z-score of −3.46 to the corresponding decoy set. We used nonequilibrium sampling techniques starting from completely extended conformations to exhaustively sample the energy surface of three nonhomologous hairpin-peptides, a three-stranded β-sheet, the all-helical 40 amino-acid HIV accessory protein, and a zinc-finger ββα motif, and find near-native conformations for the minimal energy for each protein. Using a massively parallel evolutionary algorithm, we also obtain a near-native low-energy conformation for the 54 amino-acid engrailed homeodomain. Our force field thus stabilized near-native conformations for a total of 20 proteins of all structure classes with an average RMSD of only 3.06 Å to their respective experimental conformations.     

Computer simulations of protein folding by targeted molecular dynamics

We have performed 128 folding and 45 unfolding molecular dynamics runs of chymotrypsin inhibitor 2 (CI2) with an implicit solvation model for a total simulation time of 0.4 microseconds. Folding requires that the three-dimensional structure of the native state is known. It was simulated at 300 K by supplementing the force field with a harmonic restraint which acts on the root-mean-square deviation and allows to decrease the distance to the target conformation. High temperature and/or the harmonic restraint were used to induce unfolding. Of the 62 folding simulations started from random conformations, 31 reached the native structure, while the success rate was 83% for the 66 trajectories which began from conformations unfolded by high-temperature dynamics. A funnel-like energy landscape is observed for unfolding at 475 K, while the unfolding runs at 300 K and 375 K as well as most of the folding trajectories have an almost flat energy landscape for conformations with less than about 50% of native contacts formed. The sequence of events, i.e., secondary and tertiary structure formation, is similar in all folding and unfolding simulations, despite the diversity of the pathways. Previous unfolding simulations of CI2 performed with different force fields showed a similar sequence of events. These results suggest that the topology of the native state plays an important role in the folding process.     

Conformation of the Ras-binding domain of Raf studied by molecular dynamics and free energy simulations

Recognition of Ras by its downstream target Raf is mediated by a Ras-recognition region in the Ras-binding domain (RBD) of Raf. Residues 78–89 in this region occupy two different conformations in the ensemble of NMR solution structures of the RBD: a fully α-helical one, and one where 87–90 form a type IV β-turn. Molecular dynamics simulations of the RBD in solution were performed to explore the stability of these and other possible conformations of both the wild-type RBD and the R89K mutant, which does not bind Ras. The simulations sample a fully helical conformation for residues 78–89 similar to the NMR helical structures, a conformation where 85–89 form a 3₁₀-helical turn, and a conformation where 87–90 form a type I |iB-turn, whose free energies are all within 0.3 kcal/mol of each other. NOE patterns and H^α chemical shifts from the simulations are in reasonable agreement with experiment. The NMR turn structure is calculated to be 3 kcal/mol higher than the three above conformations. In a simulation with the same implicit solvent model used in the NMR structure generation, the turn conformation relaxes into the fully helical conformation, illustrating possible structural artifacts introduced by the implicit solvent model. With the Raf R89K mutant, simulations sample a fully helical and a turn conformation, the turn being 0.9 kcal/mol more stable. Thus, the mutation affects the population of RBD conformations, and this is expected to affect Ras binding. For example, if the fully helical conformation of residues 78–89 is required for binding, its free energy increase in R89K will increase the binding free energy by about 0.6 kcal/mol. Proteins 31:186–200, 1998. © 1998 Wiley-Liss, Inc.     

Recent developments in structural proteomics for protein structure determination

The major challenges in structural proteomics include identifying all the proteins on the genome-wide scale, determining their structure-function relationships, and outlining the precise three-dimensional structures of the proteins. Protein structures are typically determined by experimental approaches such as X-ray crystallography or nuclear magnetic resonance (NMR) spectroscopy. However, the knowledge of three-dimensional space by these techniques is still limited. Thus, computational methods such as comparative and de novo approaches and molecular dynamic simulations are intensively used as alternative tools to predict the three-dimensional structures and dynamic behavior of proteins. This review summarizes recent developments in structural proteomics for protein structure determination; including instrumental methods such as X-ray crystallography and NMR spectroscopy, and computational methods such as comparative and de novo structure prediction and molecular dynamics simulations.     

Stability analysis for the cavity-filling mutations of the Myb DNA-binding domain utilizing free-energy calculations

We have analyzed the effect of cavity-filling mutations on protein stability by means of free-energy calculations based on molecular dynamics simulations to identify the factors contributing to stability changes caused by the mutations. We have studied the DNA-binding domain of Myb, which has a cavity in one of three homologous repeat units, and analyzed a series of mutations with nonnatural and natural amino acids at a single site, which change the size of the cavity. We found that the calculated free-energy changes caused by the mutations are in excellent agreement with experimental data (correlation coefficient 0.98). The free-energy changes in the native and denatured states were independently compared with the unfolding free-energy change (deltadeltaG) and cavity-volume changes (deltaV), and it was found that deltadeltaG and deltaV correlate with the native-state free-energy changes but not with the denatured-state free-energy changes. Further analyses in terms of enthalpy and entropy show that compensation between entropy and enthalpy occurs in the denatured state but not in the native state. The main contribution to the native-state free energy was found to be van der Waals interactions associated with the cavity. We estimate that the decrease in free energy per methylene group, which results from filling the cavity, is about 2 to 3 kcal/mol. These results suggest that the stabilization of a protein by cavity-filling mutations be determined primarily by the free energy associated with the cavity volume in the native state.     

Effective energy function for proteins in solution

A Gaussian solvent-exclusion model for the solvation free energy is developed. It is based on theoretical considerations and parametrized with experimental data. When combined with the CHARMM 19 polar hydrogen energy function, it provides an effective energy function (EEF1) for proteins in solution. The solvation model assumes that the solvation free energy of a protein molecule is a sum of group contributions, which are determined from values for small model compounds. For charged groups, the self-energy contribution is accounted for primarily by the exclusion model. Ionic side-chains are neutralized, and a distance-dependent dielectric constant is used to approximate the charge-charge interactions in solution. The resulting EEF1 is subjected to a number of tests. Molecular dynamics simulations at room temperature of several proteins in their native conformation are performed, and stable trajectories are obtained. The deviations from the experimental structures are similar to those observed in explicit water simulations. The calculated enthalpy of unfolding of a polyalanine helix is found to be in good agreement with experimental data. Results reported elsewhere show that EEF1 clearly distinguishes correctly from incorrectly folded proteins, both in static energy evaluations and in molecular dynamics simulations and that unfolding pathways obtained by high-temperature molecular dynamics simulations agree with those obtained by explicit water simulations. Thus, this energy function appears to provide a realistic first approximation to the effective energy hypersurface of proteins.     

Protein-inhibitor flexible docking by a multicanonical sampling: native complex structure with the lowest free energy and a free-energy barrier distinguishing the native complex from the others

Flexible docking between a protein (lysozyme) and an inhibitor (tri-N-acetyl-D-glucosamine, tri-NAG) was carried out by an enhanced conformational sampling method, multicanonical molecular dynamics simulation. We used a flexible all-atom model to express lysozyme, tri-NAG, and water molecules surrounding the two bio-molecules. The advantages of this sampling method are as follows: the conformation of system is widely sampled without trapping at energy minima, a thermally equilibrated conformational ensemble at an arbitrary temperature can be reconstructed from the simulation trajectory, and the thermodynamic weight can be assigned to each sampled conformation. During the simulation, exchanges between the binding and free (i.e., unbinding) states of the protein and the inhibitor were repeatedly observed. The conformational ensemble reconstructed at 300 K involved various conformational clusters. The main outcome of the current study is that the most populated conformational cluster (i.e., the cluster of the lowest free energy) was assigned to the native complex structure (i.e., the X-ray complex structure). The simulation also produced non-native complex structures, where the protein and the inhibitor bound with different modes from that of the native complex structure, as well as the unbinding structures. A free-energy barrier (i.e., activation free energy) was clearly detected between the native complex structures and the other structures. The thermal fluctuations of tri-NAG in the lowest free-energy complex correlated well with the X-ray B-factors of tri-NAG in the X-ray complex structure. The existence of the free-energy barrier ensures that the lowest free-energy structure can be discriminated naturally from the other structures. In other words, the multicanonical molecular dynamics simulation can predict the native complex structure without any empirical objective function. The current study also manifested that the flexible all-atom model and the physico-chemically defined atomic-level force field can reproduce the native complex structure. A drawback of the current method is that it requires a time consuming computation due to the exhaustive conformational sampling. We discussed a possibility for combining the current method with conventional docking methods.     

Conformations of an adenine bulge in a DNA octamer and its influence on DNA structure from molecular dynamics simulations.

Molecular dynamics simulations have been applied to the DNA octamer d(GCGCA-GAAC). d(GTTCGCGC), which has an adenine bulge at the center to determine the pathway for interconversion between the stacked and extended forms. These forms are known to be important in the molecular recognition of bulges. From a total of ~35 ns of simulation time with the most recent CHARMM27 force field a variety of distinct conformations and subconformations are found. Stacked and fully looped-out forms are in excellent agreement with experimental data from NMR and x-ray crystallography. Furthermore, in a number of conformations the bulge base associates with the minor groove to varying degrees. Transitions between many of the conformations are observed in the simulations and used to propose a complete transition pathway between the stacked and fully extended conformations. The effect on the surrounding DNA sequence is investigated and biological implications of the accessible conformational space and the suggested transition pathway are discussed, in particular for the interaction of the MS2 replicase operator RNA with its coat protein.