Learning To Fold Proteins Using Energy Landscape Theory

4.1 All-atom models versus template-free coarse-grained models

When starting to build a model to perform protein structure prediction, it might initially seem appealing to include as many details of the structure and of the interactions as possible. Once the form of interatomic interactions has been determined, the only task left to perform would seem to be a tuning of parameters which dictate the relative strength of the interactions. Classical mechanical, pairwise additive, explicit solvent varieties of these “all-atom” models exist [9, 13], and have recently found some success in folding small proteins [67]. One of the main disadvantages of the full blown atomistic approach is, however, that because of the very many degrees of freedom and the roughness of fully atomistic energy landscapes, especially before they were properly tuned, the search for the global free energy minimum is computationally expensive. This has made the development of all-atom force fields an arduous task, which has however met with success [67, 7].

The use of coarse-grained models [63, 38, 69, 54, 20] can simultaneously ameliorate both of the problems that lead to difficult computations. By describing fewer degrees of freedom, the forces involved in a coarse-grained model are much faster to compute and by keeping only the important degrees of freedom many of the local minima are also eliminated. For example, the solvent minima corresponding to amorphous ice and ice clathrates are avoided. These are involved in many barriers in real folding. Fortunately, so long as enough structural detail is kept, and an appropriate optimization procedure is performed, predictive, transferable coarse-grained models can be derived. Strategies for how to optimize and refine these types of models will be discussed in subsequent sections.

Most of the speedup achieved by coarse-grained models actually comes from integrating over the solvent degrees of freedom, not the missing protein atoms. This is reasonable because solvent motions are typically fast compared to protein backbone motions that involve crossing dihedral angle barriers, and those motions that are not relatively fast can be aliased onto the Hamiltonian that only explicitly depends on the coordinates of the backbone atoms. The Hamiltonian is therefore a solvent-averaged free energy function, not literally a potential energy function, and the amount of time computing forces is dramatically reduced by going to the coarse-grained description.

It is important to remember that when building coarse-grained models for the purpose of structure prediction alone you do not always have to follow all the rules that nature follows. For example, it may be useful when designing a Hamiltonian that performs structure prediction via simulated annealing from an extended conformation to intentionally lower the barriers for dihedral flips in order to allow the configurations to be explored more rapidly than actually would occur in nature. Likewise, drying effects that come from expelling solvent when preformed subunits approach each other can give large barriers that slow kinetics for some protein folding processes [15]. We can also take advantage of the robustness of the funnel in evolution: given that we know many sequences that fold into the same structure, imperfections in the optimized parameters trained on a finite set of proteins can be overcome by averaging the basic force field for a given sequence over the Hamiltonians of many sequences within a family of related proteins [52, 53, 46]. By working through the development of coarse-grained models, and seeing what succeeds and what fails, it has also been possible to get a feel for what the most important aspects of protein physics are and what amount of information and detail is actually necessary to not only predict structure but also understand folding kinetic mechanisms where getting the barriers right is often more important.

As the number of experimentally resolved structures increases, the fraction of newly resolved structures with genuinely novel folds that have not been seen before has been decreasing. A higher and higher fraction of sequences are therefore becoming good candidates for homology modeling, and homology modeling will remain a most reliable way of performing structure prediction for most sequences for some time. Nevertheless, the kinds of physically motivated coarse-grained models discussed in this review are useful for more than just tertiary structure prediction. The ability of a coarse-grained model to perform structure prediction is an important benchmark of its adequacy, and can be taken as a sign that the model is realistic enough to be used for other purposes; several examples of such applications to mechanistic questions are given in Section 9. Many interesting molecular biological phenomena, especially at the cellular level, are still well beyond the reach of all-atom simulations in terms of their time and length scales; optimized coarse-grained models are thus likely to be useful for many years to come. The combination of coarse-grained and all-atom models has already proved useful in many structure prediction applications [92].

4.2 Tertiary interactions, non-additivity and cooperativity

As more structural details of a protein model are integrated over, the appropriate form of the model energy function becomes increasingly less obvious. Building a coarse-grained model still retains elements of an art. Nonetheless, certain statistical mechanical principles are useful to keep in mind. For example, as a model becomes more and more coarse-grained, the pairwise additive approximation between degrees of freedom becomes less and less safe. It is therefore frequently useful explicitly to include contextual information about the local sequence and its environment to modulate otherwise pairwise additive interactions. An example from our own work is the water mediated interaction introduced in the AMW model. The motivation for implementing the water mediated interaction’s particular functional form came from the desire to test whether the binding landscapes of hydrophilic and hydrophobic protein-protein interfaces were funneled or perhaps were rugged, leading to difficult binding search problems. Without the use of the water mediated interaction, even landscape optimized model energy functions that correctly predicted many monomer protein structures could only correctly predict the structure of hydrophobic interfaces [82]. It was discovered that having water mediated interactions lead to better funneled folding landscapes [81], outside of the binding context in which it was originally motivated. The water mediated interaction is a sequence dependent pairwise contact interaction that switches smoothly between two different interaction weights depending upon the degree of burial of the interacting residues. It is therefore a non-additive potential. The switching function is illustrated in Figure 3. If either of the two residues participating in the interaction are buried, the residues are assumed to be interacting indirectly through protein and that interaction is assigned a particular weight. If, on the other hand, both residues are exposed, the residues are assumed to be interacting indirectly through a water molecule, and are assigned a different interaction weight.

The water mediated interaction story illustrates the necessary interplay between observations, prediction quality analysis and implementation when developing coarse-grained models. The story also illustrates how the particular chosen functional form of a model determines the ultimate success of structure prediction. Part of the water mediated interaction that was introduced was a plausible one-residue burial propensity potential, another example of a non-pairwise additive potential, which sorts residues into their preferred burial environment (buried, partially buried or exposed). Residue-residue contact potentials and a single-residue burial propensity are two of the most commonly used styles of coarse-grained interactions. The motivation of this form is clear simply from looking at protein structures, where one finds a la Kauzmann, that certain residues prefer to be buried while others prefer to be exposed. At the next level of description, the pair level, it is clear that some residues prefer to be close to each other while others do not. Yet success requires the pair interactions to be modulated by burial, then leading to a very non-additive functional form.

Protein folding of small globular domains is empirically a cooperative process, but sub-folding events are also known to be cooperative. One example of a cooperative subfolding event is the formation of hydrogen bonds between two β strands. This cooperativity can be explicitly introduced into coarse-grained models [56, 46]. Explicitly introducing cooperativity into structure based models has been shown to be useful in achieving realistic barriers to folding [26] and in predicting and understanding hydrogen-exchange experiments [19]. Having realistically large barriers may actually hinder structure prediction schemes by impeding the search for low energy states. Nevertheless, achieving a realistic degree of cooperativity in predictive coarse-grained models remains an important challenge. Recent work indicates that the spatial range of the interaction potentials in coarse-grained models play a dominant role in determining the cooperativity of the model and that realistic cooperativity is obtained by making the interaction ranges consistent with desolvation physics [51].

4.3 Backbone, steric and short range in sequence interactions

The strongest constraints on protein structure are given by its very polymeric nature and by the specific chemical nature of its backbone which leads to stereo-typical possibilities of hydrogen bonding. Due to the Pauli exclusion principle, atoms essentially never overlap at biologically relevant energy and temperature scales, so it is important in coarse-grained models also to try to minimize the overlap that would occur if a higher resolution (all-atom) model were being used. Of course overlap at higher resolution cannot be entirely avoided when only a subset of the protein’s atoms are considered, however. Thankfully, the local structures preferred by the protein backbone are surprisingly few. Most residues in natural protein structures can be unambiguously classified as being in either a helical or sheet conformation. In a coarse-grained model, these dominant local patterns can be satisfied by imposing a potential that acts on the dihedral angles of the backbone. Even when these constraints are taken into consideration, however, the specificity of the local conformation and how it relates to the local sequence is hard to capture in coarse-grained models due to the importance of local steric effects from overlap of sidechains with internal conformational freedom. These excluded volume effects are not explicitly included when the sidechain is treated with a unified representation. As a result, it is common to supplement the Hamiltonian for tertiary interactions with more information at the local-in-sequence level. Simply using fragments of possibly related experimentally determined structures is a popular method. Another method, used in AMW, is to perform bioinformatic alignments of global sequences to experimentally determined complete structures so as to enhance local compatibility in a mean field sense. This local information can then be used to impose gentle constraints on the relative distances of the backbone atoms nearby in sequence based on the analogous distances in the input candidate structures. Also one can use all-atom simulations of peptides to get structures that can determine the local biases [61]. Similar methods can be used to impose constraints on sidechain rotamer conformations.

The overall fold of a protein is determined by the trace of its backbone atoms. The trace of the backbone atoms is in turn determined by the secondary structural elements that form and the packing of these secondary structural elements. The representation of the backbone, therefore, is crucial to a good coarse-grained model of proteins. Ideally a backbone model restricts the backbone atoms to arranging themselves in realistic conformations while using a minimum number of degrees of freedom in order to keep computational efficiency. To this end, assumptions about bond lengths and planarity of the peptide bond are often employed, and dihedral angle potentials are used to mimic local steric effects. Steric effects in general can be an important factor in determining a coarse-grained model’s ability to discriminate between what would be allowed and what would be disallowed configurations in a more detailed model. Usually coarse-grained models allow too many configurations that might lead to conflicts at an all-atom level, and screening the results of coarse-grained models with more detailed models is a useful practice. Nevertheless, the steric constraints of all-atom models can give an appearance of a much more rugged landscape than is correct simply because small flexible adjustments in all-atom conformation [71] can usually avoid the worst clashes.

5.1 Decoy structures

The epigram at the beginning of this article is deliberately ambiguous. It can be understood in multiple ways in the context of developing energy functions for protein structure prediction. First, through the continual development and refinement of energy functions for various structure prediction tasks, guided by the principles of energy landscape theory as well as by new experiments, misconceptions about protein physics and errors in the functional form of the potentials have turned into an increasingly coherent understanding of protein folding at the same time as more accurate and useful energy functions are developed. A more important lesson from the motto, however, is that the mistakes that have to be corrected when performing protein structure prediction correspond to the many possible misfolded configurations that must be simultaneously destabilized as the native state is preferentially stabilized. Generating mistakes is thus an important part of the landscape learning process. Generating mistakes has not only been a problem for people making predictions but also for evolution which has solved the problem, through the trial and error of natural selection. These misfolded configurations are sometimes called “decoy structures”. This term often brings to mind a small, fixed set while in fact the decoys span a cosmologically large space. So the problem is how to use a small set to say something about the whole space: thus, statistical energy landscape theory.

Decoy structures can be generated in a number of ways, each of which has advantages and disadvantages. The simplest way of generating decoy structures is merely to shuffle the sequence of a protein while keeping its structure fixed. This would imitate the pairings in a highly mixed set of molten globule structures. This is very computationally inexpensive but is unrealistic, particularly since there are, in reality, strong constraints on where particular types of amino acids can be placed. One such example is that of membrane proteins. In membrane proteins, the residues that reside in the hydrocarbon layer are almost completely hydrophobic, while the residues in the phosphate layer are enriched in polar and charged amino acids. Therefore, shuffling the sequence completely randomly creates unrealistic decoys that exaggerate the contribution of certain types of interactions to the stabilization of the native state (such as the polar-polar or charge-charge interactions, in this case). Another related way of generating decoy structures is to simply thread the native sequence, without shuffling it, but possibly allowing gaps, over a series of unrelated, known tertiary structures. Care must of course be taken to only thread over structures with at least as many amino acids as are in the native sequence, and this method is much better than simple shuffling. Both of these methods of finding a set of representative decoys still have the assumption that the set of misfolded structures is independent of the parameters of the Hamiltonian. This is, of course, only an approximation, since energy landscapes based on pairs end up being correlated landscapes. The “right decoys” must be found, leading to an iterative procedure.

5.2 Self-consistent optimization

The preferred method for generating decoy structures would be to explicitly generate those decoys that the Hamiltonian you are trying to optimize would actually have as kinetic traps in folding. While computationally expensive, requiring iteration to self-consistency, this repeated process of trial and error only needs to be done a few times for a given form of the Hamiltonian. Explicit decoy generation and self-consistent iteration is the best way to take into account the correlations that are present in the sequences and landscapes of natural proteins. Thus, it is the best way of optimizing a set of parameters that will discriminate against the realistic decoy structures that would appear in a prediction attempt. Carrying out self-consistent optimization also requires an initial guess for the parameters in order to explicitly generate the decoys. Using shuffling or threading is a good way of generating an initial guess for the set of parameters in the Hamiltonian.

5.3 Optimization functionals

Given an (unparameterized) energy function, a set of decoy structures and a set of native structures, you might think that all that is needed is to be certain that the native fold is a bit more stable than the decoys. This would be linear programming problem [70]. The problem is, however, you never have a complete set of decoys and to be transferable, the potential needs to yield a sizable energy gap. One must be able to compute an average quantity that guarantees success with all possible decoys as competitors requiring using knowledge about the whole configuration space while knowing only a sample of the permissible decoy space. Energy landscape theory fortunately tells us the right quantity that should be optimized in order to best discriminate between native and misfolded structures that will make the native state kinetically accessible. This quantity is the ratio T_f/T_g, or, equivalently, the ratio of the energy gap between the folded and misfolded states divided by the standard deviation of the misfolded energies δE/ΔE. If this parameter is large enough, the landscape is funneled and will provide a good thermodynamic discrimination between the misfolded and folded states at temperatures where the dynamics of search are still fast enough to access kinetically the native state (far from T_g). Simultaneously optimizing an appropriate average of the ratio of T_f/T_g for a set of training proteins helps to ensure that the optimized parameters are as transferable as possible to proteins outside of the training set. This process requires some type of averaging over example proteins, which is to some extent arbitrary, but tricks such as weighting the contribution of each protein to the average in a way inversely proportional to its T_f/T_g help to prevent the average from being dominated by only a few proteins with large T_f/T_g and thus favors being able to fold the worst cases [113]. Even when a large training set is used, there is a large number of interaction parameters to be determined. Care, therefore, must also be taken to avoid noisy assignment of parameters. To prevent noisy interactions from dominating the contribution to the energy, a filtering scheme based on eigenvalue decomposition of the interaction matrices is used so that for examples of an interaction that is sampled poorly, their influence in the learning is minimized. This is much like the way one tries to avoid learning superstitions by asking for robustness to coincidence.

5.4 Constraints on optimization

The generic polymeric nature of proteins which allows them to be collapsed or random coil, or microphase separated, for example, as well as specific structural peculiarities of proteins, make it useful to constrain other properties characterizing the landscape while optimizing the parameters to give funneled landscapes. This leads to a constrained optimization problem. The additional constraints are necessary to control additional phase transitions such as the collapse transition. These constraints enter into the optimization through Lagrange multipliers. Any physical consideration which can be expressed via summation over energy terms in the Hamiltonian, such as measures of local rigidity, collapse and contributions from various distance-range interactions, can be constrained in this manner. The degree of collapse is particularly relevant to control because the kinetics within a collapsed and misfolded ensemble is considerably slower from even just the steric constraints than it is within an expanded and also unfolded state, and of course, the larger the number of strong interactions that form the deeper the trap that results. When generating decoy structures, a generic collapse bias can therefore be applied to ensure sampling of the “worst case” - i.e., folding from the collapsed state. The preformation of secondary structure can lower the barrier to folding but also has the effect of slowing down rearrangements. There are simple physical arguments that suggest that the relative contributions of the local-in-sequence versus long range interactions are comparable [98]. This idea can be used to constrain relative contributions of these terms in the potential parameters. The variance of the local-in-sequence or long range interaction energy distributions also needs to be constrained so as to minimize the probability of there being a glass transition on short length scales before the global one occurs. For example, a large variance in the local-in-sequence energy distributions may lead to dynamical freezing of those local-in-sequence interactions at T > T_g [88]. Besides T_f/T_g optimization, other relatedx optimization functionals have been proposed and tested. Most of these include maximizing the free energy gap or native state occupancy in some way. In general, statistical landscape theory shows these measures of landscape funneledness are monotonically related to the T_f/T_g criterion which we have generally employed.

5.5 The mathematics of the optimization

The mathematics of the optimization is simplest when the parameters that enter the energy function do so in a linear fashion, E = Σ_i γ_iϕ_i. The γ_i’s are the strengths of the interaction terms whereas the ϕ_i’s are monomials encoding the basic forms of the interaction potential. The stability gap can be written as δE_s = Aγ whereas the energetic variance can be written as a general quadratic function ΔE² = γBγ. A and γ are vectors of dimensionality equal to the number of interaction types while B is a matrix. A and B are defined as

where 〈ϕ_i〉_mg and ϕ_n are, for a particular interaction type, the average ϕ_i of the decoy states and native state, respectively. The optimization of $\delta E_s/\mathrm{\Delta }E=\mathit{A}\mathit{\gamma }/\sqrt{\mathit{\gamma }\mathit{B}\mathit{\gamma }}$ with respect to the set of γ_i’s is equivalent to the maximization of Aγ under the linear constraint that $\sqrt{\mathit{\gamma }\mathit{B}\mathit{\gamma }}$ is constant. That is, one optimizes with respect to the vector γ the functional, $R=\mathit{A}\mathit{\gamma }-{\lambda }_1\sqrt{\mathit{\gamma }\mathit{B}\mathit{\gamma }}$ where the Lagrange multiplier, λ₁, sets the energy scale. The solution of this geometric problem amounts to solving a system of linear equations γ = B⁻¹A up to a scalar multiple. This is worked out in the Appendix. We may also control the collapse temperature T_c = A′γ where $A_i^{\prime }={\langle {\varphi }_i\rangle }_{mg}/N$, the average ϕ_i of the decoy states divided by the number of residues in the protein, by imposing an additional constraint on our optimization functional. This ensures the decoys to be generated by the Hamiltonian will have a similar degree of collapse. Optimizing the new functional, $R=[\mathit{A}-{\lambda }_2{\mathit{A}}^{\prime }]\gamma -{\lambda }_1\sqrt{\mathit{\gamma }\mathit{B}\mathit{\gamma }}$ again yields a solution, γ = B⁻¹[A − λ₂A′] up to a scalar multiple. The Lagrange multiplier, λ₂ can be chosen to maintain the ratio of T_f/T_c = 1. The variance of the energy of the molten globules coming from different length scales can also be controlled by imposing an additional constraint on a new fluctuation matrix B′ which is constructed using only the contributions from a particular length scale in the potential, such as the local-in-sequence interactions. The Lagrange multiplier constraining this new fluctuation matrix may be chosen to hold γB′γ/γBγ constant. The mean energy of the molten globules coming from different length scales can also be controlled by separating out the contributions from the different length scales in the term, A′. For the constraint, ${\sum }_{n=1}^k{\lambda }_k{\mathit{A}}_{\mathit{k}}^{\prime }\mathit{\gamma }$, where the index n denotes the contributions from the kth length scale, the Lagrange multipliers λ_k can be chosen to make the contributions equal.

For the simplest case, where γ = B⁻¹A up to a scalar multiple, how is the optimization carried out in practice? One begins with constructing a training set of experimentally determined native protein structures, taking into consideration factors such as sequence homology, sequence length, whether the proteins are globular or membrane proteins, and whether the proteins require cofactors. Typically A and B⁻¹ are computed using decoys generated via shuffling and averaged over the entire training set. The initial γ and resulting Hamiltonian is then used for explicit generation of decoy structures by molecular dynamics. The decoys often include biasing to low Q or constraints on the radius of gyration to ensure the configurations sampled have a similar degree of collapse to the native structure. These explicitly generated decoys are then used to compute a new γ, and the process is iterated until γ converges, thus the optimization is self consistent. During each round of optimization, filtering of the B matrix may be required in order to minimize noise arising from poor statistics of certain types of interaction. One recomputes B⁻¹ by first computing an eigenvalue decomposition, B⁻¹ = PΛ⁻¹P⁻¹ where the columns of the P matrix are eigenvectors and Λ⁻¹ is the the inverse diagonal matrix of eigenvalues. The contributions coming from unreliable eigenmodes are discarded by zeroing out the corresponding eigenvalues in B⁻¹ below a cut-off. One may also apply a damping condition when combining γ’s between rounds of optimization in order to accelerate convergence, as is shown in Equation 17.

5.6 Physical interpretation of optimized parameters

Although no experimental information except the native structures of a limited training set has been used in parameterizing the Hamiltonian using a self-consistent T_f/T_g optimization, the resulting parameters can sometimes be usefully compared to simple physical measurements, such as hydrophobicity and secondary structure propensity [39]. An example of one such comparison is given in Figure 4. One can see in Figure 4 that the burial energy in the prediction Hamiltonian correlates with the hydrophobocity scale determined by transfer free energies from water while the secondary structure energy term correlates with the Chou-Fasman secondary structure propensity.

The physico-chemical interpretation of the the AMW/AWSEM contact potential interaction parameters has been discussed previously in detail [82, 81]. The three interaction matrices utilized by the AWSEM contact potential are shown in Figure 5, where the more positive the interaction weight, the greater the energetic stabilization. Polar interactions differ considerably among the different interaction types. Differences in polar interactions between water mediated and direct contact strengths suggest a large desolvation penalty upon formation of direct contacts. While polar residues are thus biased to water mediated contacts, the desolvation penalty, as expected, decreases the less polar the interaction. Hydrophobic interactions are strongest for direct contacts, also as expected.

We know the interaction matrices are not maximally complex from the general themes described about the physical chemistry of proteins already, but how complex is the protein folding code that comes from these optimized parameters?. To assess the effective number of amino acid flavors the energy function encodes, one can inspect the eigenvalue decomposition of each individual interaction matrix, as suggested by Wingreen and coworkers [64]. They found that the Miyazawa-Jernigan matrix could be accurately reconstructed using only the two largest eigenmodes corresponding to a hydrophobic-polar code, explaining why so much of a general nature about stability follows from hydrophobicity. On the other hand, we find approximately 10 eigenmodes are necessary to reconstruct the interactions matrices employed by AWSEM as summarized in Figure 6. Clearly the forces involved in structure selection and specificity go far beyond the hydrophobic scale.

5.7 Systematic refinement of parameters

The optimization scheme described above is an efficient way of parameterizing linear parameters in a Hamiltonian. It is sometimes advantageous to also perform further refinements on parameters that enter in the Hamiltonian in a nonlinear way such as the ranges of the interactions once a reasonable set of initial parameters has been determined. One way of performing such optimizations is to simply scan through a set of parameters and perform simulations with each Hamiltonian. This, however, can be very computationally intensive. A faster way of reliably extrapolating the results from a simulation with one Hamiltonian to one using another energy function uses the idea of computing the properties of a perturbed Hamiltonian on a set of structures that were generated using an unperturbed Hamiltonian. The Free Energy Perturbation scheme of Zwanzig [130], and more general but related cumulant expansion methods [24], provide systematic ways of extrapolating thermodynamic quantities.

In some cases, the small perturbation idea is not valid - e.g., when adjusting the excluded volume radii of the particles in a simulation. For such cases, more sophisticated statistical mechanical methods can be used. The Mayer cluster expansion method is one such example [25]. Consistent with the optimization method described above, the Mayer cluster expansion method can be used to estimate how the folding and glass transition temperatures will change when the excluded volume parameters are changed in a residue specific way.

5.8 Optimization and design

It should not escape the notice of the readers that optimizing a Hamiltonian to fold a set of amino acid sequences into a given set of structures is very similar in spirit to the problem of designing a sequence to fold into a particular structure given a fixed Hamiltonian. The problems of design and optimization are dual to each other and much of the machinery described above can be and has been adapted to the protein design problem [95].

6.1 Simulated annealing

Simulated annealing [91], the gradual cooling of a system from above to below its ordering transition temperature, has proven to be a general method for solving high dimensional optimization problems when the landscapes are funneled. Simulated annealing is typically carried out either via Monte Carlo moves or using molecular dynamics to explore configuration space. Molecular dynamics is straightforward to implement especially in parallel computer architectures but requires a Hamiltonian with continuous derivatives. Even when the the folding model is globally funneled, frustration in the native basin still can cause simulations to fall out of equilibrium at low temperature and thereby prevent the simulation from reaching the absolute lowest energy state [45]. If the landscape is sufficiently funneled, one nevertheless will still find a structure which is closely related to the true global minimum.

6.2 Biased sampling techniques and equilibration

Molecular dynamics combined with biased sampling techniques, such as umbrella sampling, can be used to obtain a global picture of the landscape and can be very helpful for assessing the structure prediction Hamiltonian itself and provides a way to assess the likely quality of the scheme’s predictions, as will be discussed in Section 7. Obtaining reliable equilibrium properties with finite sampling remains a problem, but techniques exist which can be used to test for equilibration. In one such method, samples from independent simulations are compared using the Kolmogorov-Smirnov test [25].

6.3 Flexible backbones versus fragment-assembly

Coarse-grained molecular dynamics models with flexible backbones allow configurations that would be disallowed at higher resolution due to steric clashes. Fragment-based Monte Carlo methods benefit from having realistic backbone structures, and this often translates into high quality prediction results, particularly for smaller proteins. These methods have elements of de novo prediction at the global level mixed with homology modeling at the local level. The process of folding a protein is inherently coupled to collapse. Monte Carlo search, which requires large fragments to remain rigid during moves, suffers from an inability of simple moves to rearrange in the collapsed state resulting in an inability to avoid steric clashes [47]. Then sampling even a well-funneled landscape in this way can be problematic for large structures if a simple move set is used. Direct fragment assembly also has issues with protein geometries that require small rearrangements such as β-strand alignment. The use of human insight in game playing has led to promising ways of avoiding these ultimately local dynamical problems [55].

6.4 Advanced sampling techniques

Beyond molecular dynamics and Monte Carlo methods, many advanced sampling techniques exist. Methods such as basin hopping flatten the landscape making it easier to explore and find global energy minima, but this is done at the cost of breaking detailed balance [90]. For this reason, such methods cannot be naively applied if the goal of the study is to obtain equilibrium properties, but these are a good choice if the goal is limited to finding the lowest energy states. Advanced sampling methods for high-dimensional macromolecular systems were recently reviewed in [93].

7.1 Free energy and energy landscape analysis

In order to obtain free energy profiles that can be used to objectively evaluate structure prediction Hamiltonians, it is necessary first to perform simulations biased by forces that constrain quantitative measures of the proximity of states to the known native structure to fixed values at multiple temperatures. These can then be combined and unbiased to obtain multidimensional free energy profiles using methods such as WHAM [60] or MBAR [104]. Biased simulations can be performed in many ways, and one popular method is umbrella sampling along a predefined reaction coordinate [111, 107].

Free energy landscape analysis is useful because it gives an overview of which parts of the landscape will be sampled and which will not during prediction runs, or more specifically how long it will take to sample each part of the landscape, during an unbiased equilibrium simulation. Regions that are high in free energy will not be sampled very often at equilibrium; the system will spend most of its time in ensembles corresponding to low free energy regions. Given equilibrium free energy curves, it becomes possible to make quantitative estimates about how long it would take a model to find the native basin. It is also possible to make estimates of how many distinct structures exist at various points along the reaction coordinate [23].

In an ideally funneled landscape, such as those used in structure-based modeling, the native basin will be low in free energy below the folding temperature. The folding transition temperature itself is controlled by the interplay between entropy loss and energetic stabilization when going from the unfolded to the folded state. It is therefore useful to compare the free energy profiles obtained using structure-based perfect funnel models to those obtained using predictive transferable energy functions [23].

Methods that allow for the calculation of free energy profiles can also be extended to the calculation of expectation values to characterize partially folded ensembles. One of the most interesting expectation values is that of energy versus degree of foldedness as measured by similarity of the contact map. If the expectation value of the energy versus similarity to the native state shows decreasing energy as the configurations become more native, and it follows that the gap between the native basin and unfolded basin is large in units of the variance of energies in the unfolded basin, the landscape is said to be funneled. Care must be taken in the interpretation of the results, however. Plots of E(Q) are often funneled at high temperatures where the occupation of the native state is entropically disfavored. Therefore, the degree of funneling also must be checked at or below the folding transition temperature of the model, where the energetic stabilization in the native state does indeed overcome the entropic stabilization of the unfolded states. An example of energy landscape analysis in the context of a structure prediction application is shown in Figure 7. This example was taken from [112], a structure prediction study comparing designed and natural proteins, which will be discussed further in Section 9.

7.2 Choice of order parameters

The choice of order parameters must also be done with care. Root-mean-square deviation (RMSD) is often used as a measure of distance to the native state, owing to its importance in x-ray crystallographic refinement. RMSD can be quite useful as a way comparing structures within the native basin. However, when looking at structures synoptically across the landscape, reaction coordinates like Q are more useful because they are better correlated with the (largely contact-based) energy functions that are used. They also are more indicative of folding mechanism since contacts can form while the global topology is not yet set. The RMSD of structures with quite reasonable but partial contact maps can be very large. This may give the false impression that the folding landscape resembles a golf-course and is not funneled. The fact that Q is unitless may make some people uneasy, but with a little practice it is easy to obtain an intuitive understanding of what different Q values mean. Structures around Q = 0.25 tend to be almost random looking, whereas Q = 0.4 structures have reasonably well formed secondary structure and have significant topological similarity to the native state. At Q = 0.55 the structure becomes easily recognizable and at Q = 0.7 the structure is typically less than a few Å RMSD from the native structure. This rough picture of the meaning of Q is transferable across a range of protein sizes, which aids in the comparison of different systems.

When looking at global folding reactions, Q has been found to be a useful reaction coordinate for kinetic analysis when the landscape is well funneled [17]. Likewise, when investigating specific phenomena it is often useful to create new reaction coordinates that are specific to the problem that you are investigating. For example, if you are interested in knowing how important short range in sequence interactions are in the folding, it makes sense to design additional reaction coordinates that single out the contribution to the folding from short range in sequence interactions alone. Also if a particular trap is important because of a specific frustrated part of the landscape, other coordinates may prove useful [109].

A generic formula for Q is given in Equation 20.

In Equation 20, N^pairs is a normalization factor equal to the number of terms in the sum. r_ij is the instantaneous distance between atoms or groups i and j, and $r_{ij}^N$ is this same distance but in the reference structure. σ_ij is a (potentially sequence separation dependent) width of the Gaussian which sets the degree of tolerance to deviation from the reference structure and is typically on the order of an Å. The list of pairs {ij ∈ pairs} can be chosen to be either all possible pair distances, only those pairs in contact in the native state, or in other ways - such as those in specific foldon units [100] - depending upon the application.

8.1 AMH

The AMH/AMC/AMW/AWSEM family of models is a series of coarse-grained protein folding models that have been continually developed, mostly in the Wolynes group, more recently with Clementi and in the Papoian group, over the last 24 years. The original version of the model, the Associative Memory Hamiltonian (AMH), was motivated by the neural network models of Hopfield and Little [48, 68]. In the AMH model [37], different residues have “charges”, and interactions between residues depend on the value of these charges. In early papers, these charges were empirically defined using concepts such as hydrophobic vs. hydrophilic tendencies of the residues. Much as the way spin models were setup to “recall” a particular configuration given a different configuration that was nearby in configuration space, these early models were able to “recall” the native structure from a database of input structures given the sequence or a closely similar one as input. An example prediction using this model is shown in Figure 8. It was found using analytical theory and confirmed in simulations that beyond a certain number of candidate structures (the “capacity” [38]), the simplest associative memory model becomes unable to faithfully recall the native structure given an input sequence. Such multiple memory models are useful for describing allostery [76, 77, 49, 66], however. It was later found that, in some cases, these types of models could “generalize”, i.e., could produce predicted structures that were closer to the native structure of the input sequence than any of the homologs in the database of input structures [36]. The ability to generalize was found to depend critically upon the correct choice of charges. In this way, generalization is closely related to the problem of finding the symmetries between different amino acids that are allowed by evolution. Grouping sequences by assigning similar values of the charge to “biologically symmetric” residue types effectively increases the number of sequences for which the structure can be reliably predicted. Successful predictions are shown in Figures 9 and 10. However, simple intuition about what the correct choice of charges proved to be insufficient, and so a systematic way of optimizing parameters in these types of coarse-grained models of proteins was implemented.

8.2 AMC

As discussed in Sections 3 and 5, energy landscape theory can be used to optimize parameters [42]. This was first done in the context of the AMH model to assign weights to different interactions coming from memories. In particular, self-consistent optimization of the Associative Memory Hamiltonian using T_f /T_g as the optimization function was shown to yield quantitatively correct structures. Later, the functional form of the potential was modified so as to use the Associative Memory interactions only for the short-range in sequence interactions while using a contact interaction for long-range in sequence interactions (Associative Memory with Contact, AMC [43]). This allows generalization to include arbitrary length insertions or deletions. The contact energy based model was reoptimized using a similar scheme (also based on an optimized local energy function for the initial threading to find memories [57, 58]) and was successfully applied to the problem of α-helical protein structure prediction without the use of any homologs in the Associative Memory database [44] documenting then “de novo” structure prediction capability. Two AMC structure predictions are shown in Figures 11 and 13. The original version of the AMC model had three contact wells extending out as far as 15Å and used a global alignment scheme to obtain lists of “memories” for the Associative Memory interactions. Further refinements were added later including the addition of an explicit β-hydrogen bonding potential with cooperativity between nearby hydrogen bonds and a refined Ramachandran potential that could be further biased based on secondary structure predictions to obtain more realistic distributions of the backbone dihedral angles [45].

8.3 AMW

The next major innovation came in new physics the form of the introduction of explicit water mediated interactions (Associative Memory with Water Mediated Interactions, AMW). The necessity of using the water mediated interaction emerged from a study of dimeric interfaces by Papoian and Wolynes [82]. They showed that a simple pairwise additive contact potential proved insufficient for recognizing dimeric interfaces when those interfaces contained a significant amount of water. Knowing the problem - in a sense learning from errors - the functional form of the potential could then be modified to allow for the possibility of two residues interacting with different interaction weights depending on the local density of residues around each of them, i.e., whether or not the interacting residues are buried or exposed. If both residues are exposed, and they are separated by a distance that is sufficient to allow a water molecule to fit between them, then they are said to be participating in a water mediated interaction. Once the model was reoptimized in this form, it was able to recognize both dry and wet interfaces. There were some surprises in the resulting potential. For example, it turns out that oppositely charged residues sometimes have an effectively favorable interaction because they are interacting with the same water molecule or perhaps an ion from the solvent. The AMW model was later used to predict the structures of monomeric proteins and the same water mediated interactions which allowed wet interfaces in dimers to be recognized were shown to also be important in the way larger proteins fold [81]. Figure 12 shows a prediction on an early CASP competition target.

To see whether a completely physically motivated model lacking any local structural biases as input could be used to determine local structures, a model having only water mediated interactions along with an α-helical hydrogen bonding terms was also used to predict the structure of α-helical proteins [78]. While some results from this pure physics based model are fine, in general the prediction results obtained were significantly worse than predictions that also used bioinformatic alignments to determine associative memory forces to guide the formation of local structures. The difference between using secondary structure prediction information and not using any secondary structure input can be seen clearly in Figure 14.

8.4 AWSEM

There are many ways to choose the short range guiding forces: all-atom physical simulations of fragments [61], global sequence threading onto templates [57], local sequence similarity of peptide fragments [47, 20]. The associative memory, water mediated, structure and energy model, or AWSEM, is a concrete open-source instantiation of the AMW model. It uses rather simple direct bioinformatic searches to identify local fragments to use as memories and uses these memories to dictate the short range in sequence interactions. Otherwise, the interactions are made up of a combination of contact terms including the water mediated interactions and α and β hydrogen bonding potentials. The primary advantage of AWSEM over previous instantiations of the model is that AWSEM is integrated into the LAMMPS molecular dynamics package [86] and is therefore fully open source. A summary of the Hamiltonian is given in Equation 21.

In Equation 21, V_backbone consists of several terms which are responsible for maintaining the connectivity of the polymeric peptide backbone, ensuring correct chirality of the amino acids, restricting the conformations to reasonable dihedral angles and preventing overlap of the chain with itself. V_contact consists of a direct, pairwise-additive contact term as well as the nonpairwise-additive water mediated interaction discussed previously. V_burial sorts amino acid types into their preferred burial environments. V_HB is made up of both α-helical and β-strand hydrogen bonding terms. V_FM uses information from bioinformatic alignments to bias local-in-sequence configurations. A complete description of the model is given in the Supplementary Information of a recent paper by Davtyan et al. [20].

9.1 Prediction of monomeric protein structures using simulated annealing of AWSEM with and without the use of homology

How well can the structures of monomeric proteins be predicted today? We tested the ability of AWSEM to predict the structures of single domain proteins when varying degrees of homology information was assumed to be known. The associative memory interactions come from the structures in the database of “memories” to bias local in sequence structures and thus exercises assuming no homology (ab initio) or acknowledging homology can be set up [20]. The actual native structure of target sequence was never used to inform any of the predictions we discuss below. Figure 15 shows a summary of the results.

The results in Figure 15 were obtained by performing multiple independent simulated annealing runs starting from an unfolded and extended conformation. The “homologs excluded” results correspond to de novo predictions in the sense that no structural information from any homologous sequence was used to inform the prediction. Homologs are rigorously excluded from the memory list in this exercise. A few example structures predicted without any homology input are shown aligned to their native structures in Figure 16 for the smallest (1r69) and largest (2fha) proteins studied. These structures are close to what a simple matching with a homolog would give if it were to be made available. As expected, including more homology information (“homologs allowed”) improves the structure prediction quality but only modestly if these are not recognized as homologs to start with. If homologs are already recognized, we can then employ “homologs only” as memories. For sequences where homologous sequences had experimentally solved structures, using exclusively structures from homologous sequences in the database of associative memories to bias the local structure formation significantly improves the best sampled structures. The best sampled structures when homologs only are used as short range input were always within a few Å of the best structures produced by MODELLER [32], a popular all-atom homology modeling tool.

9.2 Natural versus Designed Protein Landscapes, Full versus Simplified Amino Acid Alphabets

The AWSEM potential takes amino acid identity seriously but we have seen the interaction matrix is not maximally complex (i.e., only 10 principal components can be used to largely reconstitute the interaction matrix). Can a simplified folding alphabet then be used, as in the earlier concept of “biological symmetries”? Likewise, is evolution the main story (through the fragments) or can proteins designed by humans also be predicted? To address these issues, the energy landscapes of evolved and some designed proteins were studied using AWSEM. To look at specific kinetic issues we also employed a non-additive structure based model that has a higher (and more realistic) degree of cooperativity than AWSEM now has [112]. The designed sequences chosen for the study were Top7, from the Baker group, and two sequences designed and synthesized by the Takada group [59, 50]. Top7 was designed to fold to a novel topology starting from a “sketch” of the topology and its initial sequence was generated by using fragments with consistent secondary structure [59]. The design procedure was then iterated by the Baker group using Monte Carlo based sequence design and gradient based backbone optimization for multiple rounds. The two sequences designed by the Takada group began with a target scaffold of a relaxed structure of protein G-related albumin binding domain and then sequences that were expected to fold to this structure were found by a search in sequence space motivated by two criteria inspired by landscape theory [50]. One sequence, which we call TakadaE, was designed based on a scheme that merely minimized the target structure energy over sequences, while the other sequence, TakadaZ, was designed using a T_f /T_g criterion. Takada’s designs were based on an energy function for structure prediction that used many of the landscape tools and ideas we have already sketched. Analysis of Top7 in the study was carried out alongside S6, a natural protein having similar secondary structure elements but a different wiring. Top7 has been the focus of other theoretical investigations that reach similar conclusions to our own study [125, 126]. We also then compared the designs TakadaE and TakadaZ to the behavior of the natural sequence of protein G-related albumin binding domain, which we refer to as TakadaN. We also can study with AWSEM how robust the of folding of these sequences is to simplification of their amino acid code.

AWSEM with memories derived from bioinformatic alignments (while excluding fragments from homologous sequences, “homologs excluded”) was applied to predict the structure of Top7 and S6 via simulated annealing. For Top7, we found several predicted structures which were of excellent quality, the best structure having 2.1Å C_α RMSD (Figure 17). To study kinetics, a flavor of AWSEM which uses the native structure as the only “memory” for the short range in sequence associative memory interaction was used to compute free energy profiles as functions of radius of gyration (R_g) and fraction of native pairwise distances, Q_w of Top7 and S6. The free energy profiles looked remarkably similar for Top7 and S6 with the only distinguishing feature being a modestly larger range of R_g over which Top7 is low in free energy. We also found, however, several topologically distinct structures which were comparable in energy to the best predicted structure. These structures, characterized by mispairing of β strands, are consistent with the suggestion of Baker’s group from their kinetic studies [115] and the discussion of their possible role by Chan [125, 126]. The quality of S6 structure prediction is low in comparison to that observed in the survey discussed in the previous section. We concluded that this lower quality was primarily attributed to a poor bioinformatic prediction of secondary structure, an input to the model which biases AWSEM’s Ramachandran potential and β hydrogen bonding terms. AWSEM was able to successfully predict the structures of TakadaE, TakadaZ, and TakadaN. The quality of the predictions for natural sequence was the highest. The better structure prediction quality of TakadaE when compared to TakadaZ was largely due to greater funneling of the contact energy for TakadaE. Apparently the AWSEM energy function is not completely correlated with the one used in the original design by Takada.

Modern natural proteins use twenty types of amino acids, but many evolved proteins must, in addition to being foldable, be functional and bind to other proteins, so perhaps that is why the folding palette is so complex. The evolutionary constraints on function are possible sources of energetic frustration in folding [29]. Simplification schemes to two letter and to five letter codes were adapted from Wang et al. [114, 65] and the structures of the simplified sequences were predicted with AWSEM. These structure prediction studies suggest that a two-letter code would be insufficient to fold either the natural or the designed proteins, but a five-letter code may indeed be sufficient for designed sequences. The structure prediction quality for the five letter simplified variant of TakadaN was considerably poorer than the native sequence, whereas the predictions for five letter simplified variants of Top7, TakadaE, and TakadaZ do not degrade in quality with respect to that found for their respective full sequences. TakadaN’s poorer structure prediction quality can be explained by looking at the total energy as a function of Q_w (Figure 18). While the full sequence is funneled to high Q, the five letter sequence displays an energetic trap, which arises from competing secondary structures.

9.3 Prediction of Protein Binding Sites and Structure

The water mediated interaction was introduced by the landscape analysis of protein binding interfaces by Papoian et al. [82]. Recently we have investigated the dynamics of protein-protein association using the Associative-memory, Water mediated, Structure and Energy Model (AWSEM), a coarse-grained protein folding model which tests those ideas inspired by the Principle of Minimal Frustration. The parameters used in the AWSEM code were optimized by maximizing the ratio of the folding temperature to the glass transition temperature in order to create funneled folding landscapes for individual monomeric proteins, even though that optimization started with information about shuffled interfaces in dimers. Can dimer interfaces still be predicted? As shown in Figs. 19, simulated annealing using the AWSEM code is indeed able to predict successfully the native interfaces of 8 homodimers and 4 heterodimers; thus, AWSEM amounts to a flexible docking algorithm [127]. In these examples, the memory terms use local structural information about the monomers, much like the “homologs only” example discussed above. The success of the model in predicting binding sites and complete binding structures while the training set contains only monomers buttresses the idea that the same energy landscape principles that are applicable monomeric folding also apply to binding processes.

In addition to allowing interface and binding site prediction, the potential also allows us to study the role of non-native intermonomeric contacts in the process of dimer formation. Homodimers are often categorized as being either obligatory or nonobligatory dimers, meaning that the monomers must associate in order to complete their folding (obligatory) or are stably folded even in isolation at physiological temperature (nonobligatory). Non-native interactions play different roles for obligatory and non-obligatory dimers as seen in Fig. 20. An example of the free energy profile of an obligatory dimer, Arc repressor, is shown on the top of Fig. 20. States stabilized by non-native interactions correspond to on-pathway intermediates that catalyze the association process through a fly-casting mechanism [105]; the individual monomers, which are both in extended conformations before the association, have significantly larger capture radii than those of the folded monomers. The large capture radius increases the rate of binding. In the case of non-obligatory dimers as in the lower panel of Fig. 20, however, the states with non-native contacts generally appear to be off-pathway and impede binding by acting as kinetic traps.

9.4 Misfolding and frustration

As we have seen in optimizing potentials for structure prediction, protein folding and misfolding bear a yin-yang relationship in the energy landscape theory approach [121]. Can we say anything using predictive landscapes about misfolding in the laboratory? From a purely physical viewpoint, the driving forces must essentially be of the same type for misfolding as those for proper folding. For evolved proteins, which satisfy the Principle of Minimal Frustration, domain-swapped interactions are therefore the most obvious candidate for specific interactions that drive misfolding [124]. Like their counterparts in the monomer, the native contacts, domain-swapped contacts are in general stronger than other contacts. These same strong interactions can also drive oligomerization via domain-swapping with nearby domains, as suggested in [4]. Misfolding can occur, however, in a different way. AWSEM simulations show that the formation of self-recognition contacts, which are strong contacts formed between the same sequence segments from different polypeptides [99], is also a possibility. These self-recognition contacts between two segments in different molecules can be extremely strong since the segments are pretty rigid locally and these strong interactions can act like Velcro to hold two different molecules together. Unlike domain swapping interactions, self-recognition contacts have no exact counterpart in the native structure and therefore can only be involved in misfolding. A survey shows they have been avoided but not completed eliminated by evolution [101, 40].

These ideas were explored using the AWSEM code by studying fused dimers consisting of the 27th Ig domain of human cardiac titin (I27; PDB ID 1TIU). Fig. 21 shows the energy and free energy profiles for the I27–I27 fused dimer using two order parameters: the fraction of native contacts Q and the sum of the number of self-recognition contacts N_self and the number of swapped contacts N_swap. The misfolded state I is energetically less stable than the native state N but is entropically more favored since it is disordered in the other non-Velcro parts of the molecule. This type of metastable ensemble acts as a kinetic trap even below the folding temperature. When a fused dimer construct becomes trapped in this metastable state, it can act to initiate aggregation when other fused dimers are present nearby.

9.5 Initiation and branching in aggregation

Mature fibrils are a very striking feature found in patients suffering from protein aggregation related diseases, but how these structures relate to disease pathology remains open [94, 16, 3, 5]. Recent evidence supports the notion that in some cases fibers are byproducts of a process that starts with oligomers that are themselves toxic and pathogenic. It is then vitally important to understand the early stages of aggregation that are invisible to many experimental techniques. AWSEM simulations are proving useful for proposing candidate misfolded structures and other oligomeric species. Recently we have studied this problem for higher oligomers of I27. In the sequence of I27, there are two segments that can recognize themselves. According to the AWSEM energy function [129], one of these is weaker than the other. AWSEM simulations show that misfolded multidomain structures can arise from multiple independent polymer units that are cross-linked to other chains via these self-recognition interactions. Because β-strands can be linked to two other β-strands via backbone hydrogen bonding, a single domain of I27 can be thought of as a polymer unit having two reactive groups [33, 108] corresponding to the self-recognizing fragment. In the case where there is only a single self-recognition segment, each chain can only form two crosslinks. Therefore the polymerization would proceed linearly through a combination of elongation and breakup events, and later through the association of protofibril structures into mature fibrils. The simulations show [128], however, that when the number of reactive groups per chain exceeds 2, as is the case here, branched structures become possible and completely ordered fibrils must compete with these other structures. In AWSEM simulations, simulated annealing of the tetramer yielded protofibril-like structures 30% of the time and branched structures 33% of the time. The fibrillar and branched structures are illustrated in Fig. 22. The growth kinetics and molecular weight distribution of branched and linear aggregates should be very different, and this difference should be discernible in experiment. Structure prediction algorithms thus can give inspiration to new studies of the misfolding kinetics. The specificity of the Velcro self-recognizing segments and the local structural uniqueness of the oligomers also give us hope for druggability of these misfolded structures with other peptides [106, 123] or small molecules.