Testing a new Monte Carlo algorithm for protein folding

We demonstrate that the recently proposed pruned-enriched Rosenbluth method (PERM) (Grassberger, Phys. Rev. E 56:3682, 1997) leads to extremely efficient algorithms for the folding of simple model proteins. We test it on several models for lattice heteropolymers, and compare it to published Monte Carlo studies of the properties of particular sequences. In all cases our method is faster than the previous ones, and in several cases we find new minimal energy states. In addition to producing more reliable candidates for ground states, our method gives detailed information about the thermal spectrum and thus allows one to analyze thermodynamic aspects of the folding behavior of arbitrary sequences. Proteins 32:52–66, 1998. © 1998 Wiley-Liss, Inc.     

Monte Carlo studies on equilibrium globular protein folding. II. Beta-barrel globular protein models

In the context of dynamic Monte Carlo simulations on a model protein confined to a tetrahedral lattice, the interplay of protein size and tertiary structure, and the requirements for an all-or-none transition to a unique native state, are investigated. Small model proteins having a primary sequence consisting of a central bend neutral region flanked by two tails having an alternating hydrophobic/hydrophilic pattern of residues are seen to undergo a continuous transition to a beta-hairpin collapsed state. On increasing the length of the tails, the beta-hairpin structural motif is found to be in equilibrium with a four-member beta-barrel. Further increase of the tail length results in the shift of the structural equilibrium to the four-member beta-barrel. The random coil to beta-barrel transition is of an all-or-none character, but while the central turn is always the desired native bend, the location of the turns involving the two external strands is variable. That is, beta-barrels having the external stands that are two residues out of register are also observed in the transition region. Introduction into the primary sequence of two additional regions that are at the very least neutral toward turn formation produces an all-or-none transition to the unique, native, four-member beta-barrel. Various factors that can augment the stability of the native conformation are explored. Overall, these folding simulations strongly indicate that the general rules of globular protein folding are rather robust--namely, one requires a general pattern of hydrophobic/hydrophilic residues that allow the protein to have a well-defined interior and exterior and the presence of regions in the amino acid sequence that at the very least are locally indifferent to turn formation. Since no site-specific interactions between hydrophobic and hydrophilic residues are required to produce a unique four-member beta-barrel, these simulations strongly suggest that site specificity is involved in structural fine-tuning.     

A replica exchange Monte Carlo algorithm for protein folding in the HP model

Background  

The ab initio protein folding problem consists of predicting protein tertiary structure from a given amino acid sequence by minimizing an energy function; it is one of the most important and challenging problems in biochemistry, molecular biology and biophysics. The ab initio protein folding problem is computationally challenging and has been shown to be -hard even when conformations are restricted to a lattice. In this work, we implement and evaluate the replica exchange Monte Carlo (REMC) method, which has already been applied very successfully to more complex protein models and other optimization problems with complex energy landscapes, in combination with the highly effective pull move neighbourhood in two widely studied Hydrophobic Polar (HP) lattice models.     

Monte Carlo studies on equilibrium globular protein folding. III. The four helix bundle

The nature of the equilibrium conformational transition from the denatured state to a four-member alpha-helical bundle was studied employing a dynamic Monte Carlo algorithm in which the model protein chain was confined to a tetrahedral lattice. The model chain was allowed to hunt over all phase space, the target native state was not assumed a priori, and no site-specific interactions were introduced. The exterior vs the interior part of the protein is distinguished by the pattern of hydrophilic and hydrophobic interactions encoded into the primary sequence. The importance of a statistical preference for forming bends, as a function of bend location in the primary sequence, and helical wheel type cooperative interactions were examined, and the necessary conditions for collapse of the chain to the unique native structure were investigated. It was found that an amphipathic pattern of hydrophobic/hydrophilic interactions along with a statistical preference of the central residues for bend formation are sufficient to obtain the four-helix bundle. The transition to the native state has an all-or-none character.     

Confinement effects on the thermodynamics of protein folding: Monte Carlo simulations

The effects of chaperonin-like cage-induced confinement on protein stability have been studied for molecules of varying sizes and topologies. Minimalist models based on Gō-like interactions are employed for the proteins, and density-of-states-based Monte Carlo simulations are performed to accurately characterize the thermodynamic transitions. This method permits efficient sampling of conformational space and yields precise estimates of free energy and entropic changes associated with protein folding. We find that confinement-driven stabilization is not only dependent on protein size and cage radius, but also on the specific topology. The choice of the confining potential is also shown to have an effect on the observed stabilization and the scaling behavior of the stabilization with respect to the cage size.     

Dynamic Monte Carlo simulations of globular protein folding/unfolding pathways. II. Alpha-helical motifs

Dynamic Monte Carlo simulations of the folding pathways of alpha-helical protein motifs have been undertaken in the context of a diamond lattice model of globular proteins. The first question addressed in the nature of the assembly process of an alpha-helical hairpin. While the hairpin could, in principle, be formed via the diffusion-collision-adhesion of isolated performed helices, this is not the dominant mechanism of assembly found in the simulations. Rather, the helices that form native hairpins are constructed on-site, with folding initiating at or near the turn in almost all cases. Next, the folding/unfolding pathways of four-helix bundles having tight bends and one and two long loops in the native state are explored. Once again, an on-site construction mechanism of folding obtains, with a hairpin forming first, followed by the formation of a three-helix bundle, and finally the fourth helix of the native bundle assembles. Unfolding is essentially the reverse of folding. A simplified analytic theory is developed that reproduces the equilibrium folding transitions obtained from the simulations remarkably well and, for the dominant folding pathway, correctly identifies the intermediates seen in the simulations. The analytic theory provides the free energy along the reaction co-ordinate and identifies the transition state for all three motifs as being quite close to the native state, with three of the four helices assembled, and approximately one turn of the fourth helix in place. The transition state is separated from the native conformation by a free-energy barrier of mainly energetic origin and from the denatured state by a barrier of mainly entropic origin. The general features of the folding pathway seen in all variants of the model four-helix bundles are similar to those observed in the folding of beta-barrel, Greek key proteins; this suggests that many of the qualitative aspects of folding are invariant to the particular native state topology and secondary structure.     

All-atom Monte Carlo simulation of GCAA RNA folding

We report a detailed all-atom simulation of the folding of the GCAA RNA tetraloop. The GCAA tetraloop motif is a very common and thermodynamically stable secondary structure in natural RNAs. We use our simulation methods to study the folding behavior of a 12-base GCAA tetraloop structure with a four-base helix adjacent to the tetraloop proper. We implement an all-atom Monte Carlo (MC) simulation of RNA structural dynamics using a Go potential. Molecular dynamics (MD) simulation of RNA and protein has realistic energetics and sterics, but is extremely expensive in terms of computational time. By coarsely treating non-covalent energetics, but retaining all-atom sterics and entropic effects, all-atom MC techniques are a useful method for the study of protein and now RNA. We observe a sharp folding transition for this structure, and in simulations at room temperature the state histogram shows three distinct minima: an unfolded state (U), a more narrow intermediated state (I), and a narrow folded state (F). The intermediate consists primarily of structures with the GCAA loop and some helix hydrogen bonds formed. Repeated kinetic folding simulations reveal that the number of helix base-pairs forms a simple 1D reaction coordinate for the I-->N transition.     

Monte Carlo studies on equilibrium globular protein folding. I. Homopolymeric lattice models of beta-barrel proteins

Dynamic Monte Carlo studies have been performed on various diamond lattice models of β-proteins. Unlike previous work, no bias toward the native state is introduced; instead, the protein is allowed to freely hunt through all of phase space to find the equilibrium conformation. Thus, these systems may aid in the elucidation of the rules governing protein folding from a given primary sequence; in particular, the interplay of short- vs long-range interaction can be explored. Three distinct models (A? C) were examined. In model A, in addition to the preference for trans (t) over gauche states (g⁺ and g^?) (thereby perhaps favoring β-sheet formation), attractive interactions are allowed between all nonbonded, nearest neighbor pairs of segments. If the molecules possess a relatively large fraction of t states in the denatured form, on cooling spontaneous collapse to a well-defined β-barrel is observed. Unfortunately, in model A the denatured state exhibits too much secondary structure to correctly model the globular protein collapse transition. Thus in models B and C, the local stiffness is reduced. In model B, in the absence of long-range interactions, t and g states are equally weighted, and cooperativity is introduced by favoring formation of adjacent pairs of nonbonded (but not necessarily parallel) t states. While the denatured state of these systems behaves like a random coil, their native globular structure is poorly defined. Model C retains the cooperativity of model B but allows for a slight preference of t over g states in the short-range interactions. Here, the denatured state is indistinguishable from a random coil, and the globular state is a well-defined β-barrel. Over a range of chain lengths, the collapse is well represented by an all-or-none model. Hence, model C possesses the essential qualitative features observed in real globular proteins. These studies strongly suggest that the uniqueness of the globular conformation requires some residual secondary structure to be present in the denatured state.     

Monte Carlo simulation studies on the prediction of protein folding types from amino acid composition.

In the methodology development for statistical prediction of protein structures, the founders of different methods usually selected different sets of proteins to test their predicted results. Therefore, it is hard to make a fair comparison according to the results they reported. Even if the predictions by different methods are performed for the same set of proteins, there is still such a problem: a method better that the other for one set of proteins would not necessarily remain so when applied to another set of proteins. To tackle this problem, a Monte Carlo simulation method is proposed to establish an objective criterion to measure the accuracy of prediction for the protein folding type. Such an objective accuracy is actually corresponding to the asymptotical limit genereated during the Monte Carlo simulation process. Based on that, it has been found that the average objective accuracy for predicting the all-alpha, all-beta, alpha + beta, and alpha/beta proteins by the least Euclid's distance method (Nakashima, H., K. Nishikawa, and T. Ooi. 1986. J. Biochem. 99:152-162) is 73.0% and that by the least Minkowski's distance method (Chou, P.Y. 1989. Prediction in Protein Structure and the Principles of Protein Conformation. Plenum Press. New York. 549-586) is 70.9%, indicating that the former is better than the latter. However, according to the original reports, the latter claimed a rate of correct prediction with 79.7% but the former with only 70.2%, leading to a completely opposite conclusion. This indicates the necessity of establishing an objective criterion, and a comparison is meaningful only when it is based on the objective criterion. The simulation method and the idea developed here also can be applied to examine any other statistical prediction methods.     

Alternative implementations of Monte Carlo EM algorithms for likelihood inferences

Two methods of computing Monte Carlo estimators of variance components using restricted maximum likelihood via the expectation-maximisation algorithm are reviewed. A third approach is suggested and the performance of the methods is compared using simulated data.     

Monte Carlo simulation studies on the prediction of protein folding types from amino acid composition

A number of methods to predicting the folding type of a protein based on its amino acid composition have been developed during the past few years. In order to perform an objective and fair comparison of different prediction methods, a Monte Carlo simulation method was proposed to calculate the asymptotic limit of the prediction accuracy [Zhang and Chou (1992),Biophys. J.63, 1523–1529, referred to as simulation method I]. However, simulation method I was based on an oversimplified assumption, i.e., there are no correlations between the compositions of different amino acids. By taking into account such correlations, a new method, referred to as simulation method II, has been proposed to recalculate the objective accuracy of prediction for the least Euclidean distance method [Nakashimaet al. (1986),J. Biochem.99, 152–162] and the least Minkowski distance method [Chou (1989),Prediction in Protein Structure and the Principles of Protein Conformation, Plenum Press, New York, pp. 549–586], respectively. The results show that the prediction accuracy of the former is still better than that of the latter, as found by simulation method I; however, after incorporating the correlative effect, the objective prediction accuracies become lower for both methods. The reason for this phenomenon is discussed in detail. The simulation method and the idea developed in this paper can be applied to examine any other statistical prediction method, including the computersimulated neural network method.     

Polypeptide folding with off-lattice Monte Carlo dynamics: the method

The MC dynamics of an off-lattice all-atom protein backbone model with rigid amide planes are studied. The only degrees of freedom are the dihedral angle pairs of the C-atoms. Conformational changes are generated by Monte Carlo (MC) moves. The MC moves considered are single rotations (simple moves, SM's) giving rise to global conformational changes or, alternatively, cooperative rotations in a window of amide planes (window moves, WM's) generating local conformational changes in the window. Outside the window the protein conformation is kept invariant by constraints. These constraints produce a bias in the distribution of dihedral angles. The WM's are corrected for this bias by suitable Jacobians. The energy function used is derived from the CHARMM force field. In a first application to polyalanine it is demonstrated that WM's sample the conformational space more efficiently than SM's.Abbreviations CPU Central Processing Unit - MC Monte Carlo - MCD Monte Carlo Dynamics - MD Molecular Dynamics - RMS Root-Mean-Square - RMSD Root-Mean-Square-Deviation - SM Simple Move - WM Window Move     

Dynamic Monte Carlo simulations of globular protein folding/unfolding pathways. I. Six-member, Greek key beta-barrel proteins

In the context of a simplified diamond lattice model of a six-member, Greek key beta-barrel protein that is closely related in topology to plastocyanin, the nature of the folding and unfolding pathways have been investigated using dynamic Monte Carlo techniques. The mechanism of Greek key assembly is best described as punctuated "on site construction". Folding typically starts at or near a beta-turn, and then the beta-strands sequentially form by using existing folded structure as a scaffold onto which subsequent tertiary structure assembles. On average, beta-strands tend to zip up from one tight bend to the next. After the four-member, beta-barrel assembles, there is a long pause as the random coil portion of the chain containing the long loop thrahes about trying to find the native state. Thus, there is an entropic barrier that must be surmounted. However, while a given piece of the protein may be folding, another section may be unfolding. A competition therefore exists to assemble a fairly stable intermediate before it dissolves. Folding may initiate at any of the tight turns, but the turn closer to the N terminus seems to be preferred due to well-known excluded volume effects. When the protein first starts to fold, there are a multiplicity of folding pathways, but the number of options is reduced as the system gets closer to the native state. In the early stages, the excluded volume effect exerted by the already assembled protein helps subsequent assembly. Then, near the native conformation, the folded parts reduce the accessible conformational space available to the remaining unfolded sections. Unfolding essentially occurs in reverse. Employing a simple statistical mechanical theory, the configurational free energy along the reaction co-ordinate for this model has been constructed. The free energy surface, in agreement with the simulations, provides the following predictions. The transition state is quite near the native state, and consists of five of the six beta-strands being fully assembled, with the remaining long loop plus sixth beta-strand in place, but only partially assembled. It is separated from the beta-barrel intermediate by a free energy barrier of mainly entropic origin and from the native state by a barrier that is primarily energetic in origin. The latter feature is in agreement with the "Cardboard Box" model described by Goldenberg and Creighton but, unlike their model, the transition state is not a high-energy distorted form of the native state.(ABSTRACT TRUNCATED AT 250 WORDS)     

Monte Carlo study of substrate-induced folding and refolding of lattice proteins

Many proteins can switch from one conformation to another under the influence of an external driving force, such as the binding to a specific substrate. Using a simple lattice model we show that it is feasible to design protein-like lattice proteins that can have two different conformations, depending on whether or not they are bound to a substrate. We give three different examples of such substrate-induced refolding. In addition, we have explored substrate-induced folding of lattice proteins that do not fold when free in solution. We show that such proteins can bind with the same high specificity as prefolded protein, but have a considerably lower binding free energy. In this way proteins can bind to a substrate in a way that is highly specific, yet reversible.     

Potential of genetic algorithms in protein folding and protein engineering simulations.

Genetic algorithms are very efficient search mechanisms which mutate, recombine and select amongst tentative solutions to a problem until a near optimal one is achieved. We introduce them as a new tool to study proteins. The identification and motivation for different fitness functions is discussed. The evolution of the zinc finger sequence motif from a random start is modelled. User specified changes of the lambda repressor structure were simulated and critical sites and exchanges for mutagenesis identified. Vast conformational spaces are efficiently searched as illustrated by the ab initio folding of a model protein of a four beta strand bundle. The genetic algorithm simulation which mimicked important folding constraints as overall hydrophobic packaging and a propensity of the betaphilic residues for trans positions achieved a unique fold. Cooperativity in the beta strand regions and a length of 3-5 for the interconnecting loops was critical. Specific interaction sites were considerably less effective in driving the fold.     

Dynamic Monte Carlo simulations of a new lattice model of globular protein folding, structure and dynamics

A long-standing problem of molecular biology is the prediction of globular protein tertiary structure from the primary sequence. In the context of a new, 24-nearest-neighbor lattice model of proteins that includes both alpha and beta-carbon atoms, the requirements for folding to a unique four-member beta-barrel, four-helix bundles and a model alpha/beta-bundle have been explored. A number of distinct situations are examined, but the common requirements for the formation of a unique native conformation are tertiary interactions plus the presence of relatively small (but not irrelevant) intrinsic turn preferences that select out the native conformer from a manifold of compact states. When side-chains are explicitly included, there are many conformations having the same or a slightly greater number of side-chain contacts as in the native conformation, and it is the local intrinsic turn preferences that produce the conformational selectivity on collapse. The local preference for helix or beta-sheet secondary structure may be at odds with the secondary structure ultimately found in the native conformation. The requisite intrinsic turn populations are about 0.3% for beta-proteins, 2% for mixed alpha/beta-proteins and 6% for helix bundles. In addition, an idealized model of an allosteric conformational transition has been examined. Folding occurs predominantly by a sequential on-site assembly mechanism with folding initiating either at a turn or from an isolated helix or beta-strand (where appropriate). For helical and beta-protein models, similar folding pathways were obtained in diamond lattice simulations, using an entirely different set of local Monte Carlo moves. This argues strongly that the results are universal; that is, they are independent of lattice, protein model or the particular realization of Monte Carlo dynamics. Overall, these simulations demonstrate that the folding of all known protein motifs can be achieved in the context of a single class of lattice models that includes realistic backbone structures and idealized side-chains.     

Improving genetic algorithms for protein folding simulations by systematic crossover

To improve protein folding simulations, we investigated a new search strategy in combination with the simple genetic algorithm on a two-dimensional lattice model. This search strategy, we called systematic crossover, couples the best individuals, tests every possible crossover point, and takes the two best individuals for the next generation. We compared the standard genetic algorithm with and without this new implementation for various chain lengths and showed that this strategy finds local minima with better energy values and is significantly faster in identifying the global minimum than the standard genetic algorithm.     

Monte Carlo studies of folding, dynamics, and stability in alpha-helices

Folding simulations of polyalanine peptides were carried out using an off-lattice Monte Carlo simulation technique. The peptide was represented as a chain of residues, each of which contains two interaction sites: one corresponding to the C(alpha) atom and the other to the side chain. A statistical potential was used to describe the interaction between these sites. The preferred conformations of the peptide chain on the energy surface, starting from several initial conditions, were searched by perturbations on its generalized coordinates with the Metropolis criterion. We observed that, at low temperatures, the effective energy was low and the helix content high. The calculated helix propagation (s) and nucleation (sigma) parameters of the Zimm-Bragg model were in reasonable agreement with the empirical data. Exploration of the energy surface of the alanine-based peptides (AAQAA)(3) and AAAAA(AAARA)(3)A demonstrated that their behavior is similar to that of polyalanine, in regard to their effective energy, helix content, and the temperature-dependence of their helicity. In contrast, stable secondary structures were not observed for (Gly)(20) at similar temperatures, which is consistent with the nonfolder nature of this peptide. The fluctuations in the slowest dynamics mode, which describe the elastic behavior of the chain, showed that as the temperature decreases, the polyalanine peptides become stiffer and retain conformations with higher helix content. Clustering of conformations during the folding phase implied that polyalanine folds into a helix through fewer numbers of intermediate conformations as the temperature decreases.     

3D-Lattice Monte Carlo simulations of model proteins. Size effects on folding thermodynamics and kinetics

Recently, we devised an energy scale to vary systematically amino-acid residue-solvent interactions for Monte Carlo simulations of lattice-model proteins in water. For 27-mer proteins, the folding behavior varies appreciably with the choice of interaction parameters. We now perform similar simulations with 64-mers to study the size dependence of the optimal energy parameter set for representing realistic behavior typical of many real proteins (i.e. fast folding and high cooperativity for single chains). We find that 64-mers are considerably more stable and more cooperative compared to 27-mers. The optimal interfacial-interaction-energy parameter set, however, is relatively size independent.