Monte Carlo studies of the thermodynamics and kinetics of reduced protein models: Application to small helical, β, and α/β proteins

Artificial intelligence-based multi-objective optimization protocol for protein structure refinement

MOTIVATION

Protein structure refinement is an important step of protein structure prediction. Existing approaches have generally used a single scoring function combined with Monte Carlo method or Molecular Dynamics algorithm. The one-dimension optimization of a single energy function may take the structure too far away without a constraint. The basic motivation of our study is to reduce the bias problem caused by minimizing only a single energy function due to the very diversity of different protein structures.

RESULTS

We report a new Artificial Intelligence-based protein structure Refinement method called AIR. Its fundamental idea is to use multiple energy functions as multi-objectives in an effort to correct the potential inaccuracy from a single function. A multi-objective particle swarm optimization algorithm-based structure refinement is designed, where each structure is considered as a particle in the protocol. With the refinement iterations, the particles move around. The quality of particles in each iteration is evaluated by three energy functions, and the non-dominated particles are put into a set called Pareto set. After enough iteration times, particles from the Pareto set are screened and part of the top solutions are outputted as the final refined structures. The multi-objective energy function optimization strategy designed in the AIR protocol provides a different constraint view of the structure, by extending the one-dimension optimization to a new three-dimension space optimization driven by the multi-objective particle swarm optimization engine. Experimental results on CASP11, CASP12 refinement targets and blind tests in CASP 13 turn to be promising.

AVAILABILITY AND IMPLEMENTATION

The AIR is available online at: www.csbio.sjtu.edu.cn/bioinf/AIR/.

SUPPLEMENTARY INFORMATION

Supplementary data are available at Bioinformatics online.

Balancing simulation accuracy and efficiency with the Amber united atom force field

We have analyzed the quality of a recently proposed Amber united-atom model and its overall efficiency in ab initio folding and thermodynamic sampling of two stable beta-hairpins. It is found that the mean backbone structures are quite consistent between the simulations in the united-atom and its corresponding all-atom models in Amber. More importantly, the simulated beta turns are also consistent between the two models. Finally, the chemical shifts on H alpha are highly consistent between simulations in the two models, although the simulated chemical shifts are lower than experiment, indicating less structured peptides, probably due to the omission of the hydrophobic term in the simulations. More interestingly, the stabilities of both beta-hairpins at room temperature are similar to those derived from the NMR measurement, whether the united-atom or the all-atom model is used. Detailed analysis shows high percentages of backbone torsion angles within the beta region and high percentages of native contacts. Given the reasonable quality of the united-atom model with respect to experimental data, we have further studied the simulation efficiency of the united-atom model over the all-atom model. Our data shows that the united-atom model is a factor of 6-8 faster than the all-atom model as measured with the ab initio first pass folding time for the two tested beta-hairpins. Detailed structural analysis shows that all ab initio folded trajectories enter the native basin, whether the united-atom model or the all-atom model is used. Finally, we have also studied the simulation efficiency of the united-atom model as measured in terms of how fast thermodynamic convergence can be achieved. It is apparent that the united-atom simulations reach convergence faster than the all-atom simulations with respect to both mean potential energies and mean native contacts. These findings show that the efficiency of the united-atom model is clearly beyond the per-step dynamics simulation of about 2 over the all-atom model. Thus, reasonable reduction of a protein model can be achieved with improved sampling efficiency while still preserving a high level of accuracy for applications in both ab initio folding and thermodynamic sampling. This study motivates us to develop more simplified protein models with sufficient consistency with the all-atom models for enhanced conformational sampling.

Molecular simulations of peptides: a useful tool for the development of new drugs and for the study of molecular recognition

The study of the molecular recognition and self-organization properties of peptides has emerged in recent years as a very active and diverse field of research, ranging from biomedicine to biotechnology and even to material sciences. In the case of biomedicine, peptides can be used as ligands of biological receptors to gain insights into the structural, dynamical, and chemical determinants underlying the formation of complexes and identify new effectors of biological processes of interest. In the case of biotechnology and material science, short sequences have been used to understand the sequence determinants of the formation of ordered supra-molecular structures of nanoscale dimensions. In this work, we will describe our research activities in these two areas of modern chemical biology. In the first part, we will describe the development of a new, specific, potent, and selective anticancer peptide and its use to obtain the information needed to identify a non-peptidic small molecular lead to be used as an inhibitor of cancer growth. In the second part, we will describe the introduction of a new method for the description of the self-organization process at the basis of the growth of ordered supra-molecular structures held together by weak, non-covalent, yet specific interactions.

Overcoming entropic barrier with coupled sampling at dual resolutions

An enhanced sampling method is proposed for ab initio protein folding simulations. The new method couples a high-resolution model for accuracy and a low-resolution model for efficiency. It aims to overcome the entropic barrier found in the exponentially large protein conformational space when a high-resolution model, such as an all-atom molecular mechanics force field, is used. The proposed method is designed to satisfy the detailed balance condition so that the Boltzmann distribution can be generated in all sampling trajectories in both high and low resolutions. The method was tested on model analytical energy functions and ab initio folding simulations of a beta-hairpin peptide. It was found to be more efficient than replica-exchange method that is used as its building block. Analysis with the analytical energy functions shows that the number of energy calculations required to find global minima and to converge mean potential energies is much fewer with the new method. Ergodic measure shows that the new method explores the conformational space more rapidly. We also studied imperfect low-resolution energy models and found that the introduction of errors in low-resolution models does decrease its sampling efficiency. However, a reasonable increase in efficiency is still observed when the global minima of the low-resolution models are in the vicinity of the global minimum basin of the high-resolution model. Finally, our ab initio folding simulation of the tested peptide shows that the new method is able to fold the peptide in a very short simulation time. The structural distribution generated by the new method at the equilibrium portion of the trajectory resembles that in the equilibrium simulation starting from the crystal structure.

Conformational transition free energy profiles of an adsorbed, lattice model protein by multicanonical Monte Carlo simulation

Proteins often undergo changes in internal conformation upon interacting with a surface. We investigate the thermodynamics of surface induced conformational change in a lattice model protein using a multicanonical Monte Carlo method. The protein is a linear heteropolymer of 27 segments (of types A and B) confined to a cubic lattice. The segmental order and nearest neighbor contact energies are chosen to yield, in the absence of an adsorbing surface, a unique 3x3x3 folded structure. The surface is a plane of sites interacting either equally with A and B segments (equal affinity surface) or more strongly with the A segments (A affinity surface). We use a multicanonical Monte Carlo algorithm, with configuration bias and jump walking moves, featuring an iteratively updated sampling function that converges to the reciprocal of the density of states 1/Omega(E), E being the potential energy. We find inflection points in the configurational entropy, S(E)=k ln Omega(E), for all but a strongly adsorbing equal affinity surface, indicating the presence of free energy barriers to transition. When protein-surface interactions are weak, the free energy profiles F(E)=E-TS(E) qualitatively resemble those of a protein in the absence of a surface: a free energy barrier separates a folded, lowest energy state from globular, higher energy states. The surface acts in this case to stabilize the globular states relative to the folded state. When the protein surface interactions are stronger, the situation differs markedly: the folded state no longer occurs at the lowest energy and free energy barriers may be absent altogether.

The protein folding 'speed limit'

How fast can a protein possibly fold? This question has stimulated experimentalists to seek fast folding proteins and to engineer them to fold even faster. Proteins folding at or near the speed limit are prime candidates for all-atom molecular dynamics simulations. They may also have no free energy barrier, allowing the direct observation of intermediate structures on the pathways from the unfolded to the folded state. Both experimental and theoretical approaches predict a speed limit of approximately N/100micros for a generic N-residue single-domain protein, with alpha proteins folding faster than beta or alphabeta. The predicted limits suggest that most known ultrafast folding proteins can be engineered to fold more than ten times faster.

Inter-residue and solvent-residue interactions in proteins: a statistical study on experimental structures

A large set of protein structures resolved by X-ray or NMR techniques has been extracted from the Protein Data Bank and analyzed using statistical methods. In particular, we investigate the interactions between side chains and the interactions between solvent and side chains, pointing out on the possibility of including the solvent as part of a knowledge-based potential. The solvent-residue contacts are accounted for on the basis of the Voronoi's polyhedron analysis. Our investigation confirms the importance of hydrophobic residues in determining the protein stability. We observe that in general hydrophobic-hydrophobic interactions and, more specifically, aromatic-aromatic contacts tend to be increasingly distally separated in the primary sequence of proteins, thus connecting distinct secondary structure elements. A simple relation expressing the dependence of the protein free energy by the number of residues is proposed. Such a relation includes both the residue-residue and the solvent-residue contributions. The former is dominant for large size proteins, whereas for small sizes (number of residues less than 100) the two terms are comparable. Gapless threading experiments show that the solvent-residue knowledge-based potential yields a significant contribution with respect to discriminating the native structure of proteins. Such contribution is important especially for proteins of small size and is similar to that given by the most favorable residue-residue knowledge-based potential referring to hydrophobic-hydrophobic interactions such as isoleucine-leucine. In general, the inclusion of the solvent-residue interaction produces a relevant increase of the free energy gap between the native structures and decoys.

Sequence-based study of two related proteins with different folding behaviors

Z(SPA-1) is an engineered protein that binds to its parent, the three-helix-bundle Z domain of staphylococcal protein A. Uncomplexed Z(SPA-1) shows a reduced helix content and a melting behavior that is less cooperative, compared with the wild-type Z domain. Here we show that the difference in folding behavior between these two sequences can be partly understood in terms of an off-lattice model with 5-6 atoms per amino acid and a minimalistic potential, in which folding is driven by backbone hydrogen bonding and effective hydrophobic attraction.

Ab Initio Folding of Helix Bundle Proteins Using Molecular Dynamics Simulations

We have demonstrated that ab initio fast folding simulations at 400 K using a GB implicit solvent model with an all-atom based force field can describe the spontaneous formation of nativelike structures for the 36-residue villin headpiece and the 46-residue fragment B of Staphylococcal protein A. An implicit solvent model combined with high-temperature MD makes it possible to perform direct folding simulations of small- to medium-sized proteins by reducing the computational requirements tremendously. In the early stage of folding of the villin headpiece and protein A, initial hydrophobic collapse and rapid formation of helices were found to play important roles. For protein A, the third helix forms first in the early stage of folding and exhibits higher stability. The free energy profiles calculated from the folding simulations suggested that both of the helix-bundle proteins show a two-state thermodynamic behavior and protein A exhibits rather broad native basins.

Two-state folding over a weak free-energy barrier

We present a Monte Carlo study of a model protein with 54 amino acids that folds directly to its native three-helix-bundle state without forming any well-defined intermediate state. The free-energy barrier separating the native and unfolded states of this protein is found to be weak, even at the folding temperature. Nevertheless, we find that melting curves to a good approximation can be described in terms of a simple two-state system, and that the relaxation behavior is close to single exponential. The motion along individual reaction coordinates is roughly diffusive on timescales beyond the reconfiguration time for a single helix. A simple estimate based on diffusion in a square-well potential predicts the relaxation time within a factor of two.

A simple lattice model that exhibits a protein-like cooperative all-or-none folding transition

In a recent paper (D. Gront et al., Journal of Chemical Physics, Vol. 115, pp. 1569, 2001) we applied a simple combination of the Replica Exchange Monte Carlo and the Histogram methods in the computational studies of a simplified protein lattice model containing hydrophobic and polar units and sequence-dependent local stiffness. A well-defined, relatively complex Greek-key topology, ground (native) conformations was found; however, the cooperativity of the folding transition was very low. Here we describe a modified minimal model of the same Greek-key motif for which the folding transition is very cooperative and has all the features of the "all-or-none" transition typical of real globular proteins. It is demonstrated that the all-or-none transition arises from the interplay between local stiffness and properly defined tertiary interactions. The tertiary interactions are directional, mimicking the packing preferences seen in proteins. The model properties are compared with other minimal protein-like models, and we argue that the model presented here captures essential physics of protein folding (structurally well-defined protein-like native conformation and cooperative all-or-none folding transition).

The dual role of a loop with low loop contact distance in folding and domain swapping

Alpha helices, beta strands, and loops are the basic building blocks of protein structure. The folding kinetics of alpha helices and beta strands have been investigated extensively. However, little is known about the formation of loops. Experimental studies show that for some proteins, the formation of a single loop is the rate-determining step for folding, whereas for others, a loop (or turn) can misfold to serve as the hinge loop region for domain-swapped species. Computer simulations of an all-atom model of fragment B of Staphylococcal protein A found that the formation of a single loop initiates the dominant folding pathway. On the other hand, the stability analysis of intermediates suggests that the same loop is a likely candidate to serve as a hinge loop for domain swapping. To interpret the simulation result, we developed a simple structural parameter: the loop contact distance (LCD), or the sequence distance of contacting residues between a loop and the rest of the protein. The parameter is applied to a number of other proteins, including SH3 domains and prion protein. The results suggest that a locally interacting loop (low LCD) can either promote folding or serve as the hinge region for domain swapping. Thus, there is an intimate connection between folding and domain swapping, a possible cause of misfolding and aggregation.

Folding of a small helical protein using hydrogen bonds and hydrophobicity forces

A reduced protein model with five to six atoms per amino acid and five amino acid types is developed and tested on a three-helix-bundle protein, a 46-amino acid fragment from staphylococcal protein A. The model does not rely on the widely used Go approximation, which ignores non-native interactions. We find that the collapse transition is considerably more abrupt for the protein A sequence than for random sequences with the same composition. The chain collapse is found to be at least as fast as helix formation. Energy minimization restricted to the thermodynamically favored topology gives a structure that has a root-mean-square deviation of 1.8 A from the native structure. The sequence-dependent part of our potential is pairwise additive. Our calculations suggest that fine-tuning this potential by parameter optimization is of limited use.

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.

Three-helix-bundle protein in a Ramachandran model

We study the thermodynamic behavior of a model protein with 54 amino acids that forms a three-helix bundle in its native state. The model contains three types of amino acids and five to six atoms per amino acid and has the Ramachandran torsional angles phi(i), psi(i) as its degrees of freedom. The force field is based on hydrogen bonds and effective hydrophobicity forces. For a suitable choice of the relative strength of these interactions, we find that the three-helix-bundle protein undergoes an abrupt folding transition from an expanded state to the native state. Also shown is that the corresponding one- and two-helix segments are less stable than the three-helix sequence.

Two state behavior in a solvable model of beta-hairpin folding

Understanding the mechanism of protein secondary structure formation is an essential part of the protein-folding puzzle. Here we describe a simple model for the formation of a beta hairpin, motivated by the fact that folding of a beta hairpin captures much of the basic physics of protein folding. The modeled hairpin is composed of two interacting Gaussian chains with one pairwise (two-body) and two many-body interactions. We show that these many-body interactions, arising from side chain packing effects, are responsible for producing an "all-or-none" folding transition. We also estimate the (single exponential) folding/unfolding rate via calculating the thermodynamic weight of the "critical" droplet/bubble.

Folding of a model three-helix bundle protein: a thermodynamic and kinetic analysis

The kinetics and thermodynamics of an off-lattice model for a three-helix bundle protein are investigated as a function of a bias gap parameter that determines the energy difference between native and non-native contacts. A simple dihedral potential is used to introduce the tendency to form right-handed helices. For each value of the bias parameter, 100 trajectories of up to one microsecond are performed. Such statistically valid sampling of the kinetics is made possible by the use of the discrete molecular dynamics method with square-well interactions. This permits much faster simulations for off-lattice models than do continuous potentials. It is found that major folding pathways can be defined, although ensembles with considerable structural variation are involved. The large gap models generally fold faster than those with a smaller gap. For the large gap models, the kinetic intermediates are non-obligatory, while both obligatory and non-obligatory intermediates are present for small gap models. Certain large gap intermediates have a two-helix microdomain with one helix extended outward (as in domain-swapped dimers); the small gap intermediates have more diverse structures. The importance of studying the kinetic, as well as the thermodynamics, of folding for an understanding of the mechanism is discussed and the relation between kinetic and equilibrium intermediates is examined. It is found that the behavior of this model system has aspects that encompass both the "new" view and the "old" view of protein folding.

Exploring the origins of topological frustration: design of a minimally frustrated model of fragment B of protein A

Topological frustration in an energetically unfrustrated off-lattice model of the helical protein fragment B of protein A from Staphylococcus aureus was investigated. This G-type model exhibited thermodynamic and kinetic signatures of a well-designed two-state folder with concurrent collapse and folding transitions and single exponential kinetics at the transition temperature. Topological frustration is determined in the absence of energetic frustration by the distribution of Fersht phi values. Topologically unfrustrated systems present a unimodal distribution sharply peaked at intermediate phi, whereas highly frustrated systems display a bimodal distribution peaked at low and high phi values. The distribution of phi values in protein A was determined both thermodynamically and kinetically. Both methods yielded a unimodal distribution centered at phi = 0.3 with tails extending to low and high phi values, indicating the presence of a small amount of topological frustration. The contacts with high phi values were located in the turn regions between helices I and II and II and III, intimating that these hairpins are in large part required in the transition state. Our results are in good agreement with all-atom simulations of protein A, as well as lattice simulations of a three- letter code 27-mer (which can be compared with a 60-residue helical protein). The relatively broad unimodal distribution of phi values obtained from the all-atom simulations and that from the minimalist model for the same native fold suggest that the structure of the transition state ensemble is determined mostly by the protein topology and not energetic frustration.

Designing potential energy functions for protein folding

By following a consistent line of physical reasoning, some fundamental understanding about the foldability of proteins has been achieved. In recent years, this has led to the development of a number of successful algorithms for optimizing potential energy functions for folding protein models. The differences between the folding mechanisms of simple, contact-based lattice proteins and more traditional, realistic protein models, however, still call for further development of the potentials in addition to the optimization approaches.

Folding simulation of protein models on the structure-based cubo-octahedral lattice with the Contact Interactions algorithm

Computer simulations of protein models on lattices have been widely used as an aid in the study of protein folding process. Following the suggestion of Raghunathan and Jernigan (1997, Protein Sci. 6:2072-2083) that the cubooctahedral lattice can allow a more realistic representation of proteins than other lattices, we propose here the use of a new set of internal coordinates theta for the description of a protein model on this lattice. An easy procedure for the conversion of the theta coordinates to the Cartesian coordinates is also described. When the Contact Interaction approach, already proposed by us for simulations on square or cubic lattices, was applied to the cubo-octahedral lattice, the system obeyed the correct thermodynamics derived from the definition of energy. Thus, lattice simulations of protein models, in which secondary structure elements such as alpha-helices or beta-strands can be easily identifiable, can be performed.

Reduced protein models and their application to the protein folding problem

One of the most important unsolved problems of computational biology is prediction of the three-dimensional structure of a protein from its amino acid sequence. In practice, the solution to the protein folding problem demands that two interrelated problems be simultaneously addressed. Potentials that recognize the native state from the myriad of misfolded conformations are required, and the multiple minima conformational search problem must be solved. A means of partly surmounting both problems is to use reduced protein models and knowledge-based potentials. Such models have been employed to elucidate a number of general features of protein folding, including the nature of the energy landscape, the factors responsible for the uniqueness of the native state and the origin of the two-state thermodynamic behavior of globular proteins. Reduced models have also been used to predict protein tertiary and quaternary structure. When combined with a limited amount of experimental information about secondary and tertiary structure, molecules of substantial complexity can be assembled. If predicted secondary structure and tertiary restraints are employed, low resolution models of single domain proteins can be successfully predicted. Thus, simplified protein models have played an important role in furthering the understanding of the physical properties of proteins.