A knowledge-based move set for protein folding

The free energy landscape of protein folding is rugged, occasionally characterized by compact, intermediate states of low free energy. In computational folding, this landscape leads to trapped, compact states with incorrect secondary structure. We devised a residue-specific, protein backbone move set for efficient sampling of protein-like conformations in computational folding simulations. The move set is based on the selection of a small set of backbone dihedral angles, derived from clustering dihedral angles sampled from experimental structures. We show in both simulated annealing and replica exchange Monte Carlo (REMC) simulations that the knowledge-based move set, when compared with a conventional move set, shows statistically significant improved ability at overcoming kinetic barriers, reaching deeper energy minima, and achieving correspondingly lower RMSDs to native structures. The new move set is also more efficient, being able to reach low energy states considerably faster. Use of this move set in determining the energy minimum state and for calculating thermodynamic quantities is discussed.     

Generalized-ensemble algorithms for molecular simulations of biopolymers.

In complex systems with many degrees of freedom such as peptides and proteins, there exists a huge number of local-minimum-energy states. Conventional simulations in the canonical ensemble are of little use, because they tend to get trapped in states of these energy local minima. A simulation in generalized ensemble performs a random walk in potential energy space and can overcome this difficulty. From only one simulation run, one can obtain canonical-ensemble averages of physical quantities as functions of temperature by the single-histogram and/or multiple-histogram reweighting techniques. In this article we review uses of the generalized-ensemble algorithms in biomolecular systems. Three well-known methods, namely, multicanonical algorithm, simulated tempering, and replica-exchange method, are described first. Both Monte Carlo and molecular dynamics versions of the algorithms are given. We then present three new generalized-ensemble algorithms that combine the merits of the above methods. The effectiveness of the methods for molecular simulations in the protein folding problem is tested with short peptide systems.     

Using robotics to fold proteins and dock ligands

The problems of protein folding and ligand docking have been explored largely using molecular dynamics or Monte Carlo methods. These methods are very compute intensive because they often explore a much wider range of energies, conformations and time than necessary. In addition, Monte Carlo methods often get trapped in local minima. We initially showed that robotic motion planning permitted one to determine the energy of binding and dissociation of ligands from protein binding sites (Singh et al., 1999). The robotic motion planning method maps complicated three-dimensional conformational states into a much simpler, but higher dimensional space in which conformational rearrangements can be represented as linear paths. The dimensionality of the conformation space is of the same order as the number of degrees of conformational freedom in three-dimensional space. We were able to determine the relative energy of association and dissociation of a ligand to a protein by calculating the energetics of interaction for a few thousand conformational states in the vicinity of the protein and choosing the best path from the roadmap. More recently, we have applied roadmap planning to the problem of protein folding (Apaydin et al., 2002a). We represented multiple conformations of a protein as nodes in a compact graph with the edges representing the probability of moving between neighboring states. Instead of using Monte Carlo simulation to simulate thousands of possible paths through various conformational states, we were able to use Markov methods to calculate the steady state occupancy of each conformation, needing to calculate the energy of each conformation only once. We referred to this Markov method of representing multiple conformations and transitions as stochastic roadmap simulation or SRS. We demonstrated that the distribution of conformational states calculated with exhaustive Monte Carlo simulations asymptotically approached the Markov steady state if the same Boltzman energy distribution was used in both methods. SRS permits one to calculate contributions from all possible paths simultaneously with far fewer energy calculations than Monte Carlo or molecular dynamics methods. The SRS method also permits one to represent multiple unfolded starting states and multiple, near-native, folded states and all possible paths between them simultaneously. The SRS method is also independent of the function used to calculate the energy of the various conformational states. In a paper to be presented at this conference (Apaydin et al., 2002b) we have also applied SRS to ligand docking in which we calculate the dynamics of ligand-protein association and dissociation in the region of various binding sites on a number of proteins. SRS permits us to determine the relative times of association to and dissociation from various catalytic and non-catalytic binding sites on protein surfaces. Instead of just following the best path in a roadmap, we can calculate the contribution of all the possible binding or dissociation paths and their relative probabilities and energies simultaneously.     

Folding simulations of small proteins

Understanding how a protein folds is a long-standing challenge in modern science. We have used an optimized atomistic model (united-residue force field) to simulate folding of small proteins of various structures: HP-36 (alpha protein), protein A (beta), 1fsd (alpha+beta), and betanova (beta). Extensive Monte Carlo folding simulations (ten independent runs with 10(9) Monte Carlo steps at a temperature) starting from non-native conformations are carried out for each protein. In all cases, proteins fold into their native-like conformations at appropriate temperatures, and glassy transitions occur at low temperatures. To investigate early folding trajectories, 200 independent runs with 10(6) Monte Carlo steps are also performed at a fixed temperature for a protein. There are a variety of possible pathways during non-equilibrium early processes (fast process, approximately 10(4) Monte Carlo steps). Finally, these pathways converge to the point unique for each protein. The convergence point of the early folding pathways can be determined only by direct folding simulations. The free energy surface, an equilibrium thermodynamic property, dictates the rest of the folding (slow process, approximately 10(8) Monte Carlo steps).     

A heteropolymer model study for the mechanism of protein folding

A heteropolymer model of randomly self-interacting chains in two dimensions is studied with numerical simulations in order to elucidate the folding mechanism of protein. We find that the model occasionally shows folding propensity depending on the sequence of random numbers given to the chain. We study the thermodynamic and kinematic roles in the folding mechanism by grouping the local energy minima found in the simulations into clusters according to the similarity of their conformations. It is suggested that the local minima to which some heteropolymers show a folding tendency are always the lowest energy states of the energy spectrum within a cluster, though which cluster is selected depends on the sequence. For the eight random sequences we study, we find that the energy gap between the ground state and excited states is little correlated with folding or nonfolding. We rather find that folding propensities are correlated with the global structure of the average energy surface, implying a dominant kinetic role in the folding mechanism, although thermal factors cannot be ignored as the mechanism of choosing the ground state within a cluster of states connected by small deformations. We suggest that a hierarchical cluster structure plays an important role in selecting a unique folded state out of the huge number of local minima of heteropolymers. © 1997 John Wiley & Sons, Inc.     

Determination of the conformation of folding initiation sites in proteins by computer simulation

Experimental evidence and theoretical models both suggest that protein folding begins by specific short regions of the polypeptide chain intermittently assuming conformations close to their final ones. The independent folding properties and small size of these folding initiation sites make them suitable subjects for computational methods aimed at deriving structure from sequence. We have used a torsion space Monte Carlo procedure together with an all-atom free energy function to investigate the folding of a set of such sites. The free energy function is derived by a potential of mean force analysis of experimental protein structures. The most important contributions to the total free energy are the local main chain electrostatics, main chain hydrogen bonds, and the burial of nonpolar area. Six proposed independent folding units and four control peptides 11–14 residues long have been investigated. Thirty Monte Carlo simulations were performed on each peptide, starting from different random conformations. Five of the six folding units adopted conformations close to the experimental ones in some of the runs. None of the controls did so, as expected. The generated conformations which are close to the experimental ones have among the lowest free energies encountered, although some less native like low free energy conformations were also found. The effectiveness of the method on these peptides, which have a wide variety of experimental conformations, is encouraging in two ways: First, it provides independent evidence that these regions of the sequences are able to adopt native like conformations early in folding, and therefore are most probably key components of the folding pathways. Second, it demonstrates that available simulation methods and free energy functions are able to produce reasonably accurate structures. Extensions of the methods to the folding of larger portions of proteins are suggested. © 1995 Wiley-Liss, Inc.     

Replica-exchange methods and predictions of helix configurations of membrane proteins

A simulation in generalized ensemble is based on a non-Boltzmann weight factor and performs a random walk in potential energy space, which allows the simulation to avoid getting trapped in states of local-minimum energy states. In this article, we review uses of the generalized-ensemble algorithms. Three well-known methods, namely, multicanonical algorithm (MUCA), simulated tempering (ST) and replica-exchange method (REM), are described first. Both Monte Carlo (MC) and molecular dynamics (MD) versions of the algorithms are given. We then present the results of the application of replica-exchange MC method to the predictions of membrane protein structures.     

Computer simulations of de novo designed helical proteins.

In the context of reduced protein models, Monte Carlo simulations of three de novo designed helical proteins (four-member helical bundle) were performed. At low temperatures, for all proteins under consideration, protein-like folds having different topologies were obtained from random starting conformations. These simulations are consistent with experimental evidence indicating that these de novo designed proteins have the features of a molten globule state. The results of Monte Carlo simulations suggest that these molecules adopt four-helix bundle topologies. They also give insight into the possible mechanism of folding and association, which occurs in these simulations by on-site assembly of the helices. The low-temperature conformations of all three sequences have the features of a molten globule state.     

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.     

Structural basis of hierarchical multiple substates of a protein. I: Introduction

A computer experiment of protein dynamics is carried out, which consists of two steps: (1) A Monte Carlo simulation of thermal fluctuations in the native state of a globular protein, bovine pancreatic trypsin inhibitor; and (2) a simulation of the quick freezing of fluctuating conformations into energy minima by minimization of the energy of a number of conformations sampled in the Monte Carlo simulation. From the analysis of results of the computer experiment is obtained the following picture of protein dynamics: multiple energy minima exist in the native state, and they are distributed in clusters in the conformational space. The dynamics has a hierarchical structure which has at least two levels. In the first level, dynamics is restricted within one of the clusters of minima. In the second, transitions occur among the clusters. Local parts of a protein molecule, side chains and local main chain segments, can take multiple locally stable conformations in the native state. Many minima result from combinations of these multiple local conformations. The hierarchical structure in the dynamics comes from interactions among the local parts. Protein molecules have two types of flexibility, each associated with elastic and plastic deformations, respectively.     

Membrane adsorption, folding, insertion and translocation of synthetic trans-membrane peptides

Spontaneous membrane adsorption, folding and insertion of the synthetic WALP16 and KALP16 peptides was studied by computer simulations starting from completely extended conformations. The peptides were simulated using an unmodified all-atom force field in combination with an efficient Monte Carlo sampling algorithm. The membrane is represented implicitly as a hydrophobic zone inside a continuum solvent modelled using the generalized Born theory of solvation. The method was previously parameterized to match insertion energies of hydrophobic side chain analogs into cyclohexane and no parameters were optimized for the present simulations. Both peptides rapidly precipitate out of bulk solution and adsorb to the membrane surface. Interfacial folding into a helical conformation is followed by membrane insertion. Both the peptide conformations and their location in the membrane are strongly temperature dependent. The temperature dependent behaviour can be summarized by fitting to a four-state model, separating the system into folded and unfolded conformers, which are either inserted into the membrane or located at the interfaces. As the temperature is lowered the dominant peptide conformation of the system changes from unfolded surface bound configurations to folded surface bound states. Folded trans-membrane conformers represent the dominant configuration at low temperatures. The analysis allows direct estimates of the free energies of peptide folding and membrane insertion. In the case of WALP the quality of the fit is excellent and the thermodynamic behaviour is in good agreement with expected theoretical consideration. For KALP the fit is more problematic due to the large solvation energies of the charged lysine residues.     

Membrane adsorption,folding, insertion and translocation of synthetic trans-membrane peptides

Spontaneous membrane adsorption, folding and insertion of the synthetic WALP16 and KALP16 peptides was studied by computer simulations starting from completely extended conformations. The peptides were simulated using an unmodified all-atom force field in combination with an efficient Monte Carlo sampling algorithm. The membrane is represented implicitly as a hydrophobic zone inside a continuum solvent modelled using the generalized Born theory of solvation. The method was previously parameterized to match insertion energies of hydrophobic side chain analogs into cyclohexane and no parameters were optimized for the present simulations. Both peptides rapidly precipitate out of bulk solution and adsorb to the membrane surface. Interfacial folding into a helical conformation is followed by membrane insertion. Both the peptide conformations and their location in the membrane are strongly temperature dependent. The temperature dependent behaviour can be summarized by fitting to a four-state model, separating the system into folded and unfolded conformers, which are either inserted into the membrane or located at the interfaces. As the temperature is lowered the dominant peptide conformation of the system changes from unfolded surface bound configurations to folded surface bound states. Folded trans-membrane conformers represent the dominant configuration at low temperatures. The analysis allows direct estimates of the free energies of peptide folding and membrane insertion. In the case of WALP the quality of the fit is excellent and the thermodynamic behaviour is in good agreement with expected theoretical consideration. For KALP the fit is more problematic due to the large solvation energies of the charged lysine residues.     

Phylogenetic tree construction using sequential stochastic approximation Monte Carlo

Monte Carlo methods have received much attention recently in the literature of phylogenetic tree construction. However, they often suffer from two difficulties, the curse of dimensionality and the local-trap problem. The former one is due to that the number of possible phylogenetic trees increases at a super-exponential rate as the number of taxa increases. The latter one is due to that the phylogenetic tree has often a rugged energy landscape. In this paper, we propose a new phylogenetic tree construction method, which attempts to alleviate these two difficulties simultaneously by making use of the sequential structure of phylogenetic trees in conjunction with stochastic approximation Monte Carlo (SAMC) simulations. The use of the sequential structure of the problem provides substantial help to reduce the curse of dimensionality in simulations, and SAMC effectively prevents the system from getting trapped in local energy minima. The new method is compared with a variety of existing Bayesian and non-Bayesian methods on simulated and real datasets. Numerical results are in favor of the new method in terms of quality of the resulting phylogenetic trees.     

ARTIST: an activated method in internal coordinate space for sampling protein energy landscapes

We present the first applications of an activated method in internal coordinate space for sampling all-atom protein conformations, the activation-relaxation technique for internal coordinate space trajectories (ARTIST). This method differs from all previous internal coordinate-based studies aimed at folding or refining protein structures in that conformational changes result from identifying and crossing well-defined saddle points connecting energy minima. Our simulations of four model proteins containing between 4 and 47 amino acids indicate that this method is efficient for exploring conformational space in both sparsely and densely packed environments, and offers new perspectives for applications ranging from computer-aided drug design to supramolecular assembly.     

Folding cooperativity in a three-stranded beta-sheet model

The thermodynamic behavior of a previously designed three-stranded beta-sheet was studied via several microseconds of standard and replica exchange molecular dynamics simulations. The system is shown to populate at least four thermodynamic minima, including two partially folded states in which only a single hairpin is formed. Simulated melting curves show different profiles for the C and N-terminal hairpins, consistent with differences in secondary structure content in published NMR and CD/FTIR measurements, which probed different regions of the chain. Individual beta-hairpins that comprise the three-stranded beta-sheet are observed to form cooperatively. Partial folding cooperativity between the component hairpins is observed, and good agreement between calculated and experimental values quantifying this cooperativity is obtained when similar analysis techniques are used. However, the structural detail in the ensemble of conformations sampled in the simulations permits a more direct analysis of this cooperativity than has been performed on the basis of experimental data. The results indicate the actual folding cooperativity perpendicular to strand direction is significantly larger than the lower bound obtained previously.     

Structural basis of hierarchical multiple substates of a protein. III: Side chain and main chain local conformations

An analysis is carried out of differences in the minimum energy conformations obtained in the previous paper by energy minimization starting from conformations sampled by a Monte Carlo simulation of conformational fluctuations in the native state of a globular protein, bovine pancreatic trypsin inhibitor. Main conformational differences in each pair of energy minima are found usually localized in several side chains and in a few local main chain segments. Such side chains and local main chain segments are found to take a few distinct local conformations in the minimum energy conformations. Energy minimum conformations can thus be described in terms of combinations of these multiple local conformations.     

Non-Boltzmann thermodynamic integration (NBTI) for macromolecular systems: relative free energy of binding of trypsin to benzamidine and benzylamine

The relative free energies of binding of trypsin to two amine inhibitors, benzamidine (BZD) and benzylamine (BZA), were calculated using non-Boltzmann thermodynamic integration (NBTI). Comparison of the simulations with the crystal structures of both complexes, trypsin-BZD and trypsin-BZA, shows that NBTI simulations better sample conformational space relative to thermodynamic integration (TI) simulations. The relative binding free energy calculated using NBTI was much closer to the experimentally determined value than that obtained using TI. The error in the TI simulation was found to be primarily due to incorrect sampling of BZA's conformation in the binding pocket. In contrast, NBTI produces a smooth mutation from BZD to BZA using a surrogate potential, resulting in a much closer agreement between the inhibitors' conformations and the omit electron density maps. This superior agreement between experiment and simulation, of both relative binding free energy differences and conformational sampling, demonstrates NBTI's usefulness for free energy calculations in macromolecular simulations.     

Folding of Top7 in unbiased all‐atom Monte Carlo simulations

For computational studies of protein folding, proteins with both helical and β‐sheet secondary structure elements are very challenging, as they expose subtle biases of the physical models. Here, we present reproducible folding of a 92 residue α/β protein (residues 3–94 of Top7, PDB ID: 1QYS) in computer simulations starting from random initial conformations using a transferable physical model which has been previously shown to describe the folding and thermodynamic properties of about 20 other smaller proteins of different folds. Top7 is a de novo designed protein with two α‐helices and a five stranded β‐sheet. Experimentally, it is known to be unusually stable for its size, and its folding transition distinctly deviates from the two‐state behavior commonly seen in natural single domain proteins. In our all‐atom implicit solvent parallel tempering Monte Carlo simulations, Top7 shows a rapid transition to a group of states with high native‐like secondary structure, and a much slower subsequent transition to the native state with a root mean square deviation of about 3.5 Å from the experimentally determined structure. Consistent with experiments, we find Top7 to be thermally extremely stable, although the simulations also find a large number of very stable non‐native states with high native‐like secondary structure. Proteins 2013; 81:1446–1456. © 2013 Wiley Periodicals, Inc.     

Structural basis of hierarchical multiple substates of a protein. II: Monte Carlo simulation of native thermal fluctuations and energy minimization

Conformational fluctuations in a globular protein, bovine pancreatic trypsin inhibitor, in the time range between picoseconds and nanoseconds are studied by a Monte Carlo simulation method. Multiple energy minima are derived from sampled conformations by minimizing their energy. They are distributed in clusters in the conformational space. A hierarchical structure is observed in the simulated dynamics. In the time range between 10(-14) and 10(-10) seconds dynamics is well represented by a superposition of vibrational motions within an energy well with transitions among minima within each cluster. Transitions among clusters take place in the time range of nanoseconds or longer.     

A Modified Valley Restrained Monte Carlo Method to Efficiently Search the Low Energy Structures of Peptides

Abstract

A new Monte Carlo sampling scheme, namely the Modified Valley Restrained Monte Carlo procedure, is used to obtain the global energy minimum conformations for polypeptides, such as Met-enkephalin and Melittin. For each peptide, we found close agreement with previous results from both theoretical and experimental studies. The simple idea for controlling the step size according to the Valley Function, provides useful suggestions in searching the global energy minimum structures, and furthermore helps solve the multiple minima problem.