Finding the needle in the haystack: towards solving the protein-folding problem computationally
Prediction of protein tertiary structures from amino acid sequence and understanding the mechanisms of how proteins fold, collectively known as \xe2\x80\x9cthe protein folding problem,\xe2\x80\x9d has been a grand challenge in molecular biology for over half a century. Theories have been developed that provide us with an unprecedented understanding of protein folding mechanisms. However, computational simulation of protein folding is still difficult, and prediction of protein tertiary structure from amino acid sequence is an unsolved problem. Progress toward a satisfying solution has been slow due to challenges in sampling the vast conformational space and deriving sufficiently accurate energy functions. Nevertheless, several techniques and algorithms have been adopted to overcome these challenges, and the last two decades have seen exciting advances in enhanced sampling algorithms, computational power and tertiary structure prediction methodologies. This review aims at summarizing these computational techniques, specifically conformational sampling algorithms and energy approximations that have been frequently used to study protein-folding mechanisms or to de novo predict protein tertiary structures. We hope that this review can serve as an overview on how the protein-folding problem can be studied computationally and, in cases where experimental approaches are prohibitive, help the researcher choose the most relevant computational approach for the problem at hand. We conclude with a summary of current challenges faced and an outlook on potential future directions.
Introduction<!>Thermodynamics of protein folding<!>Simulation of protein folding and tertiary structure prediction are very different subproblems<!>Chemical kinetics of protein folding: mechanisms and pathways<!>Prediction of protein-folding rates<!>Conformational sampling is a bottleneck<!>Unbiased molecular dynamics simulations<!>Enhanced sampling techniques in MD<!>Monte Carlo simulation<!>Genetic algorithms<!>Energy functions are evolving objects<!>Physics-based force fields<!>Knowledge-based potentials<!>Improving sampling and scoring with restraints<!>Sparse experimental data as restraints<!>Predicted contacts as restraints<!>Examples of methods for de novo tertiary structure prediction<!>FRAGFOLD<!>Rosetta<!>I-TASSER<!>QUARK<!>BCL::Fold<!>Outlook
<p>Protein folding is a process of molecular self-assembly during which a disordered polypeptide chain collapses to form a compact and well-defined three-dimensional (3 D) tertiary structure. A grand challenge in biochemistry has been to understand the process by which proteins fold into their functional tertiary structure (folding mechanism) and to predict this tertiary structure from amino acid sequence (structure prediction), two tasks that are collectively known as "the protein-folding problem" (Chan and Dill 1993; Dill et al. 2008; Dill and MacCallum 2012). Solving this problem is of far-reaching impact as it will not only reveal the missing link between sequence and structure but also provide molecular biologists with a theoretical framework and practical tools for applications such as drug design and protein engineering. As a result, an enormous amount of effort has been contributed to study the protein-folding problem by the scientific community. This is illustrated by Figure 1, which shows the striking growth in the number of articles published each year on this problem since Anfinsen's "thermodynamic hypothesis" of protein folding, that protein native state resides in the global minimum of Gibbs free energy, was formally stated in 1973 (Anfinsen 1973). A comprehensive review of the study of this problem is deemed impossible for an article of this kind. As many excellent review articles on the theories of protein folding and their experimental validation have been published over the years (Dill et al. 1995; Onuchic et al. 1997; Dobson et al. 1998; Dobson and Karplus 1999; Radford 2000; Onuchic and Wolynes 2004; Bartlett and Radford 2009; Bowman et al. 2011; Englander and Mayne 2014; Wolynes 2015), here we focus our discussion on computational methods for studying folding mechanisms and predicting tertiary structures. Specifically, we limit our discussion to protein-folding simulations and de novo protein structure prediction at atomic detail, as methods based on coarse-grained representation of protein structures were recently comprehensively reviewed (Kmiecik et al. 2016). In addition, due to space limitations, we are not able to cover the complete literature of this topic, and we apologize to those whose contributions have not received the deserved attention.</p><p>Nevertheless, the two key components of any folding simulation or structure prediction methods are efficient sampling of conformational space and accurate evaluation of the energy of sampled conformations. Hence, the main body of this article is devoted to discussing different algorithms and their advances toward efficient sampling of conformational space followed by approaches and progress toward accurate energy functions. To put the discussion under the theoretical framework of protein folding, we first briefly summarize different views on mechanisms of protein folding. The interplay between sampling algorithms and energy functions is concretely illustrated by discussing some representative methods shown to be relatively successful in the Critical Assessment of protein Structure Prediction (CASP) experiment (Moult et al. 1995; Tai et al. 2014). Finally, we present a summary on the progress and outline specific challenges that future development in the field will likely overcome.</p><!><p>When a protein folds, it experiences constant counteractions between the effective energy, which favors the native state, and the configurational entropy, which favors unfolded states (Karplus 2011). The term "effective energy" refers to the free energy of the system (protein plus solvent) that consists of the intramolecular energy of the protein in vacuum plus the solvation free energy (the free energy of transfer of the protein from the gas phase to solution). The Gibbs free energy of the protein-solvent system is the sum of the effective energy and the configurational entropy (Lazaridis and Karplus 1999, 2000, 2003) (Figure 2). At equilibrium, both folded and unfolded states can be characterized by their Gibbs free energy. The difference in Gibbs free energies between the native state and unfolded states is termed the free energy of folding.</p><p>Both the enthalpic and entropic contributions to ΔGfolding mainly arise from intramolecular and protein-solvent nonbonded interactions and rearrangement of solvent molecules. As calculating the exact Gibbs free energy from first principles is prohibitive (Leach 2001), a simplified energy function is used in practical computer simulations of protein dynamics, folding and structure prediction. Broadly speaking, there are two different types of approaches to a simplified energy function. The first is a classical mechanical model that describes the potential energy, which is parameterized by analyzing the fundamental forces between particles; the second is a statistical model parameterized on data derived from statistical analysis of pair interactions and other properties in known protein structures (Lazaridis and Karplus 2000). In this review, we will use the term "energy" frequently when we discuss various implementations of energy functions for evaluating the "energy" of sampled conformations; however, the reader is advised to keep in mind that such energy approximations are not physically realistic Gibbs free energies.</p><p>Solvation can be accounted for by either immersing the protein into explicit solvent molecules or including in the energy function a term that implicitly models solvation free energy. The former approach is often adopted in molecular dynamics simulations and is desirable especially in cases where the purpose is to study structural details about protein–solvent interaction. Two major limitations of this approach are that the computational expense is high and the effective energy of a protein conformation is not known. The latter approach, often referred to as implicit solvation, is typically orders of magnitude faster and compatible with more sampling techniques than corresponding simulations with explicit solvent (Lazaridis and Karplus 2003).</p><!><p>While both the prediction of protein tertiary structure and simulation of folding require efficient search of conformational space and accurate evaluation of the energy of sampled conformations, it needs to be emphasized that these two subproblems are rather different, with distinct solutions and limitations. Methods for tertiary structure prediction generally create 3D models by assembling small structural fragments or motifs, quite often, with physically unrealistic trajectory of conformational search and evaluate the energy of sampled conformations using statistical potentials. While this approach has worked quite successfully in creating models that are close to native structures (Bradley et al. 2005a; Zhang 2009; Moult et al. 2016), it has very little chance of giving insight into the mechanisms of folding. It is doubtless that if one could simulate actual folding processes, both subproblems would be solved. However, as will be explained in later sections, this is only possible for relatively small proteins using molecular dynamics simulations. Thus, methods for simulating folding mechanisms, while often employ physically realistic energy functions and can reveal important thermodynamics and kinetics about folding, are generally not useful for predicting structures for all but only small proteins.</p><!><p>The conformational space accessible to a polypeptide chain is astronomically large; a systematic search for the functional structure of a polypeptide chain with 100 residues would take an amount of time even longer than the age of the universe. The fact that proteins fold on a biologically meaningful timescale, with some attaining their functional structures in just a few microseconds, led Levinthal to conclude that there must be well-defined folding mechanisms and pathways to the native state (Levinthal 1968, 1969), so that protein folding is under "kinetic control." A full characterization of the folding process requires elucidation of the mechanisms by which transition states and intermediates, if any, are formed and the determination of whether there is a single defined pathway or multiple pathways to the native state.</p><p>The "classical view" of protein folding assumes a sequential model and postulates a well-defined sequence of intermediates that follow one another to carry the protein from the unfolded random coil to a uniquely folded native state (Levinthal 1968; Kim and Baldwin 1982, 1990). In the search for such a single mechanism of protein folding, several models have been proposed about how folding gets started and native contacts and structure are subsequently formed (Baldwin 1989; Fersht 1997; Daggett and Fersht 2003a, 2003b). The "nucleation" model postulated that a folding-initiating local secondary structure, or nucleus, is formed slowly followed by the rapid propagation of native structure in a stepwise manner (Wetlaufer 1973). However, this model was dropped from favor as it predicts the absence of folding intermediates. The "framework" model and the related "diffusion-collision" model proposed that secondary structures segments are preformed independently of tertiary structure before they diffuse and collide to give stable tertiary structure (Kim and Baldwin 1990; Karplus and Weaver 1994). The "hydrophobic collapse" model hypothesized that folding starts with a rapid collapse around hydrophobic residues to the molten globule state (compact denatured state), that narrows down the conformational exploration to the native state significantly (Baldwin 1989; Ptitsyn 1996). An essential feature of these latter models is that they predict the presence of folding intermediates. However, the fact that some proteins fold by simple two-state kinetics, without the accumulation of folding intermediates, and that secondary and tertiary structure form simultaneously led to the formulation of the "nucleation-condensation" model (Fersht 1997; Daggett and Fersht 2003a, 2003b). This model assumed the concerted formation of local and nonlocal structures and was considered a "unifying" mechanism of protein folding (Daggett and Fersht 2003 b). It should be noted, however, that the "nucleation-condensation" model does not preclude the presence of folding intermediates (Daggett and Fersht 2003a, 2003b).</p><p>The observation that molten globules form asynchronously over a range of timescales fostered the concept of protein-folding funnel (Frauenfelder et al. 1991; Bryngelson et al. 1995; Dill and Chan 1997; Onuchic et al. 1997; Dobson et al. 1998; Brooks et al. 2001; Wolynes 2015) (Figure 2). In this "new view," it is inferred that proteins must fold into their unique native state through multiple unpredictable pathways that involve the progressive organization of an ensemble of partially folded intermediates on a rugged effective energy hypersurface that resembles a funnel. The funnel shape arises from the fact that the number of accessible configurations, which determine the configurational entropy, decreases as the energy decreases (Karplus (2011). A more recent formulation of the mechanism of protein folding is centered around of concept of foldons (Englander et al. 2007; Englander and Mayne (2014). In what's called the foldon-based hypothesis, a protein starts folding by forming an initial seed foldon through unguided search, and it follows a foldon-determined folding pathway as the seed foldon guides subsequent foldons in a "folding upon binding" way. While this hypothesis states that proteins fold along a definite path after formation of the initial foldon, the foldon formation at the initial stage is assumed to be accomplished through a disordered multitrack search (Englander and Mayne 2014).</p><!><p>Folding rate is an essential parameter for characterizing protein-folding kinetics. What factors determine whether a protein will be a slow or fast folder? Can we predict folding rates from amino acid sequences as well? Theoretical studies have suggested that the size, native topology and stability of a protein influence the rate and mechanisms by which it folds. In searching for a causal relationship, a key advance made in 1998 was that in a set of 12 nonhomologous single-domain proteins folding rate shows a significant correlation with a simple measure of topological complexity of the native fold, the so-called contact order, which is defined as the average sequence separation between all pairs of native contacts normalized by sequence length (Plaxco et al. 1998). In contrast, the correlations between size or native state stability and folding rate are weak to nonexistent (Plaxco et al. 1998). Based on this observation, another parameter called long-range order, which counts long-range contacts (contacts that are close in space but distant in sequence), was proposed and found to be a strong predictor of the folding rates of two-state proteins (Gromiha and Selvaraj 2001). Contact order and long-range order have also been combined to form a parameter called total contact distance that has better correlation with folding rates (Zhou and Zhou 2002b). Folding rates were also found to be inversely correlated with a parameter called multiple contact index that measures the number of residues with multiple long-range contacts (Gromiha 2009). The correlations between the various topological parameters just discussed and folding rates suggest that it is viable to predict folding rates from amino acid sequences because native topologies are determined by amino acid sequences (Baker 2000). In fact, several bioinformatics tools have been developed for this purpose (Ivankov and Finkelstein 2004; Gromiha et al. 2006; Ouyang and Liang 2008; Chou and Shen 2009; Guo and Rao 2011), and two notable web servers are FOLD-RATE (Gromiha et al. 2006), FoldRate (Chou and Shen 2009).</p><!><p>A polypeptide chain with a typical size can adopt an astronomical number of conformations. It is agreed that conformational sampling remains to be a bottleneck of de novo structure prediction (Jones 1997a; Baker and Sali 2001; Bradley et al. 2005b; Zhang 2008; Kim et al. 2009; Maximova et al. 2016). Nevertheless, there has been exciting improvement in sampling algorithms based on statistical mechanical principles or guided by experimental or predicted restraints, all of which are further accelerated by improvements in hardware speed and power (Maximova et al. 2016). For the convenience of discussion, we divide conformational search methods into the following three broad categories: molecular dynamics simulations, Monte Carlo simulations and genetic algorithms. For each category of algorithms, we give a general formulation of the algorithm and a summary of the latest studies in which the algorithm was applied to study protein folding mechanism or de novo protein structure prediction.</p><!><p>Molecular dynamics (MD) simulation is a widely used computational technique for exploring the macroscopic properties of molecular systems through explicit computation of microscopic particle motions. MD has had enormously influential applications in biomolecular systems and has been heavily used to study motion-related phenomena such as protein folding, conformational flexibility, protein structure determination from NMR, ligand–protein interaction and protein–membrane interaction (Karplus and Petsko 1990; van Gunsteren and Berendsen 1990; Karplus and McCammon 2002; Gumbart et al. 2005; Karplus and Kuriyan 2005; Lindahl and Sansom 2008; Klepeis et al. 2009; Durrant and McCammon 2011; Periole 2017). The two essential elements of an MD simulation are the interaction potential for the particles and the equations of motion governing the dynamics of the particles (Leach 2001; Rapaport 2004). Interaction potentials will be discussed in the section: Energy functions are evolving objects. Here, we describe how MD simulations explore the phase space of a molecular system.</p><p>A typical MD run involves generation of successive microstates of a molecular system by solving Newton's equations of motion for all atoms simultaneously with femtosecond time steps (Equation (2)).
(2)mid2ridt2=Fi=−∂U(r1,r2,…,rN)∂ri
where ri and U(r1,
r2,…, rN) denote position vector and potential energy of point mass i, respectively. Fi denotes the force acted upon point mass i. The result of the simulation is a trajectory of microstates that specify how the system evolves in phase space (Leach 2001). In principle, equilibrium properties can be computed by averaging over the trajectory if it is of sufficient length to give a representative ensemble of the microstates of the system. Unfortunately, the usefulness of MD in studying long timescale biological phenomena is often limited due to inadequate sampling of all relevant conformational states of a system. Even when the energy barriers separating two topologically different low-energy regions of the conformational space are of order kBT, traversing them by random thermal fluctuation cannot be achieved within a reasonable amount of time.</p><p>A wide range of biologically interesting phenomena occurs over timescales in the order of milliseconds, several orders of magnitude beyond the reach of conventional MD simulations. As a result, studying processes that involve major conformational changes, such as protein folding, activation and deactivation, by MD simulations has been traditionally challenging (Gruebele 2002). The very first protein folding simulation via MD at the microsecond timescale was notably made by Duan and Kollman (1998), who simulated the folding process of the villin headpiece (a 36 mer) in explicit solvent for two months on parallel supercomputers. The simulation showed a mechanism for the peptide to reach a marginally stable state with a main chain RMSD of 5.7 Å from the native state (Duan and Kollman 1998). This peptide was later de novo folded by Zagrovic and coworkers (2002) to an ensemble of states whose average Cα RMSD is 1.7 Å from the native state. The total simulation time was 300 μs or approximately 1000 CPU years with the help of worldwide-distributed computers (Zagrovic et al. 2002).</p><p>Substantial progress has been made during the past decade or so to extend the folding times accessible by conventional MD simulations through efficient parallelization of MD codes or MD-specialized hardware (Lane et al. 2013) (Figure 3). The MD-specialized software package Desmond and the massively parallelized machine Anton, both developed recently at D.E. Shaw Research, have allowed for conducting millisecond timescale MD simulations of systems with tens of thousands of atoms in just a few weeks (Bowers et al. 2006; Shaw et al. 2007, 2014). Desmond is a collection of codes that implement novel parallel algorithms and numerical techniques to perform high-throughput and accurate MD simulations on conventional computational clusters, general-purpose supercomputers and GPUs (Bowers et al. 2006). Anton is built on MD-specific ASICs (application-specific integrated circuits) that interact in a tightly coupled manner using a high-speed communication network. Its ability to efficiently perform simulations on the timescales over which many physiologically relevant processes take place expands the set of problems for which the use of MD is tractable (Shaw et al. 2007, 2014). Armed with this specialized set of software and hardware, researchers at D.E. Shaw Research have been able to simulate protein folding from extended random coils (Shaw et al. 2010; Lindorff-Larsen et al. 2011) and study structural origin of slow diffusion in protein folding (Chung et al. 2015), protein-ligand recognition (Dror et al. 2011; Shan et al. 2011), mechanism of nucleotide exchange in G proteins (Dror et al. 2015), and mechanisms of kinase activation and inhibition (Shan et al. 2014; Ingram et al. 2015) at realistic timescales. The de novo folding simulations conducted at D.E. Shaw Research generated computational insights in favor of the single-pathway view of protein folding (Figure 2). For example, equilibrium simulations of WW domain captured multiple folding and unfolding events that consistently follow a well-defined folding pathway (Shaw et al. 2010). However, subsequent folding simulations of 12 fast-folding proteins showed that although a majority of them fold along a single dominant route, differing "transition state classes" were observed for two proteins (Lindorff-Larsen et al. 2011).</p><p>A different approach to overcome the sampling challenge of MD is through statistical analysis of multiple independent trajectories or aggregating independent short simulations using Markov state models (MSM) to make a complete model of system dynamics (Pande et al. 2010; Prinz et al. 2011; Lane et al. 2013). The MSM effectively pieces together this complete model from independent trajectories, allowing for prediction of kinetic phenomena on timescales much longer than the individual trajectories used to construct the model (Lane et al. 2013). While the MSM-based "multitrajectory" approach has some advantages over the reaction coordinate-based single trajectory analysis, such as identifying areas of phase space for adaptive sampling (Bowman et al. 2010, Weber and Pande, 2011), insights gained from MSM analysis does not always agree with the single-pathway view of folding. For example, while it was shown via single-trajectory analysis that folding of the WW domain follows a definite pathway where the first hairpin folds first (Shaw et al. 2010), a parallel statistically significant pathway where the second hairpin of the WW domain folds first was detected using MSM to analyze the same simulation trajectories (Lane et al. 2011). Similar analysis conducted on the MD trajectories of 12 small fast-folding proteins (Beauchamp et al. 2012), while showed that two-state model is inadequate for the same set of systems as described by a previous study (Lindorff-Larsen et al. 2011), revealed a richer picture of populated states for some more complicated systems.</p><!><p>The ruggedness of energy landscapes with many local minima separated by high-energy barriers makes adequate conformational sampling a challenging task. MD trajectories often do not reach all biologically relevant conformations, a problem that can be addressed by employing enhanced sampling algorithms (Okamoto 2004; Bernardi et al. 2015). Two popular enhanced sampling techniques in the simulations of biological systems are replica-exchange molecular dynamics (REMD) and metadynamics (Bernardi et al. 2015). While we focus our discussion on the REMD along temperature, which is also known as parallel tempering, several variants of replica exchange protocols have also been reported (Fukunishi et al. 2002; Itoh et al. 2011; Wu et al. 2012).</p><p>The replica-exchange method was developed to overcome the multitude of local minima separated by high energy barriers (Sugita and Okamoto 1999). Many molecular simulation scenarios require ergodic sampling of energy landscapes that feature many minima, and barriers between minima can be difficult to overcome at ambient temperatures over accessible simulation timescales. Replica-exchange simulations seek to enhance the sampling in such scenarios by running n non-interacting copies of the system Ci (i = 1, …, n) in parallel each at a different temperature Ti in the canonical ensemble (Figure 4). The non-interacting nature of this artificial compound system (C1, C2, …, Cn) ensures that each state's weight factor is given by the product of Boltzmann factors of each copy.</p><p>Compared to a standard Monte Carlo simulation, which affects the conformation of only one copy, REMD explores the energy landscape by periodically exchanging the conformations of replicas. The probability of transition of a compound system such that the conformations between a pair of copies (Ci, Cj) are exchanged is
(4)p=min(1,eΔ)
where
(5)Δ=(βj−βi)(Uj−Ui)
</p><p>In most cases, exchange of the conformations of replicas decreases auto-correlation, thus enabling replicas to reach thermal equilibrium faster than without exchange. However, for protein folding simulation, a recent study showed that the efficiency of REMD is not much higher than that of conventional MD if the folding rate is not very temperature-dependent (Rosta and Hummer 2009). While it is not necessary to restrict the exchange to copies with neighboring temperature (e.g. j=i+ 1), doing so will be optimal, since the transition probability decreases exponentially with the difference in temperature between copies (Hansmann 1997). It is also worth noted that while exchange of conformations between copies must be conducted in a Monte Carlo way, there is no restriction on which algorithms are used for updating the conformation of an individual copy locally. In fact, several variants of REMD have been developed (Mori et al. 2016). For example, a replica-exchange Monte Carlo (REMC) technique was implemented in the threading-based structure prediction pipeline QUARK and tested in CASP11 (Zhang et al. 2016).</p><p>Metadynamics is a class of methods that eases sampling by introducing a time-dependent biasing potential that acts on a selected number of coarse-grained order parameters, often referred to as collective variables (CVs) (Laio and Parrinello 2002; Piana and Laio 2007; Barducci et al. 2011; Valsson et al. 2016). CVs are generally nonlinear functions of the atomic positions of the simulated system that should ideally distinguish between all relevant metastable states. Some simple but informative CVs used in protein-folding simulations are number of Cα contacts, number of backbone H-bonds and helicity of the backbone and the free energy surface is usually plotted as a function of these CVs (Piana and Laio 2007). The added biasing potential is introduced through successive addition of small repulsive Gaussian kernels deposited along the system trajectory in CV space (Figure 4) (Barducci et al. 2011; Valsson et al. 2016). The added Gaussian kernel is a function of the current position and the previous position of the system in the CV space, and its intended purpose is to discourage the system from revisiting configurations that have already been sampled, thus accelerating sampling. The final summation of the deposited Gaussian kernels also gives an unbiased estimate of the free energy landscape of the system. In contrast to these advantages, it is, however, far from trivial to decide when to stop a simulation and find a set of CVs proper for describing the process of interest (Barducci et al. 2011; Valsson et al. 2016).</p><p>Both REMD and metadynamics have been used to de novo fold several small peptides and proteins. The first example of using REMD to sample a folded structure starting from a completely unfolded state is probably the study of Rhee et al. (Rheeand Pande 2003) where a 23-residue BBA5 protein was folded by what's called multiplexed REMD. Using REMD simulations in implicit solvent, Pitera et al. (Pitera and Swope 2003) folded a 20-residue designed Trp-cage peptide starting from an extended coil to a state <1.0 Å Cα RMSD from conformations in the NMR ensemble. Recently,Jiang etal.(JiangandWu2014)foldedadiverseset of 14 fast folding proteins from their unfolded states using REMD with a residue-specific force field. A similar study by Nguyen et al. (Nguyen etal. 2014) included a larger set of 17 proteins;while they successfully folded most proteins, mis-folded structures are thermodynamically preferred for 3 proteins.</p><!><p>MD simulation is without a doubt a required technique if one wishes to study folding pathway or kinetics computationally. However, for tertiary structure prediction of large proteins whose energy landscapes are populated with many local minima separated by high barriers, Monte Carlo (MC) simulation can be much more efficient (Figure 5(A)). It is, in fact, the underlying search engine of some of the most successful de novo tertiary structure prediction methods (Simons et al. 1997; Bradley et al. 2005a; Xu and Zhang 2012; Zhang et al. 2016) and our method BCL::Fold (Karakas et al. 2012). Unlike MD simulations where successive conformations of the system are connected through time, in a MC simulation, each new conformation of the system depends only upon its immediate predecessor. The technique of MC simulation was introduced as the first computer simulation of a molecular system in 1952 (Metropolis et al. 1953). Nowadays, the term "Monte Carlo" is often used to describe a simulation whenever random sampling is performed.</p><p>A MC simulation explores the phase space of a system by randomly perturbing the current conformation by actions such as moving a single atom or molecule or adjusting dihedral angles. The energy of the new conformation is then evaluated using an energy function. If the new conformation is lower in energy than its predecessor, it is accepted as a starting conformation for the next iteration. If the energy is higher, the new conformation is accepted with a probability based on the famous Metropolis criterion (Metropolis et al. 1953) (Equation (4). This is often done by comparing the Boltzmann factor of the new conformation to a random number between 0 and 1, and the new conformation is accepted if its Boltzmann factor is greater than the random number and rejected otherwise. While the essential search algorithm of MC-based structure prediction methods is the same, they differ in the starting components for assembling 3 D models and in the repertoire of MC moves implemented for perturbing the model (Vitalis and Pappu 2009).</p><p>Primitive MC sampling can be computationally expensive and thus inefficient at finding global energy minimum. Typically, these methods are coupled with some optimization technique that vastly decreases computational expense by directing the progression of the MC simulation toward global energy minimum. One optimization technique is gradient-based sampling, where MC iterations are directed down local property gradients, that is, the potential next state with the lowest energy is selected. For instance, gradients can be calculated based on side-chain rotameric states (Hu et al. 2010) or, in the HP-lattice model (Dill et al. 1995), the movement of a residue in various directions (Hu et al. 2009). However, when the conformational space is continuous rather than discrete, gradient descent becomes unfeasible because the energy cannot be calculated for every step forward. The most popular optimization approach shown to effectively accelerate the convergence of a MC simulation is probably simulated annealing (Kirkpatrick et al. 1983). The essential feature of this technique is that it combines MC sampling of conformational space at an initially elevated temperature with a proper cooling scheme over the course of the simulation. The cooling scheme, if gentle enough, theoretically ensures the system will reach the global minimum. In turn, the probability of a higher energy step being accepted decreases over time, and models are directed toward the global energy minimum (Tsallis and Stariolo 1996). Many powerful de novo tertiary structure prediction methods integrate this MC-simulated annealing approach (Kmiecik et al. 2016); we include a detailed discussion on some selected examples of such methods (see Examples of methods for de novo tertiary structure prediction).</p><!><p>Genetic algorithms (GAs) are an optimization procedure based on the process of evolution that occurs in nature. GAs have been used in a variety of applications. Some prominent ones include automatic programing, machine learning and population genetics (Goldberg 1989). Generally, a GA initializes the optimization process by randomly generating an initial population of trial solutions each encoded as a string of bits, also called a chromosome (Figure 5(B)). Offspring are produced by applying nature-inspired operations, namely mutations and crossovers on bit strings. Mutations are introduced into strings by flipping one or more bits, whereas crossovers between two individuals consist of randomly selecting a crossover site and exchanging the left segment of one string with the right segment of the other (Figure 5(B)). The fittest offspring are selected for continual refinement via the iteration of multiple generations (Schulze-Kremer 2000).</p><p>A large number of studies on the use of GAs for de novo protein structure prediction and protein folding simulation have been made (Pedersen and Moult 1996; Cui et al. 1998; Schulze-Kremer 2000; Custodio et al. 2004; Unger 2004; Hoque et al. 2009; Huang et al. 2010; Zhang et al. 2010; Custodio et al. 2014; Boskovic and Brest 2016; Rashid et al. 2016) since the pioneering work of Dendekar and Argos (Dandekar and Argos 1992) on de novo folding simulation of a model protein of a four β-strand bundle and that of Unger and Moult (Unger and Moult 1993) on searching for global energy minimum on the 2D HP lattice model. The simplest protein representations used in GAs is the 2 D HP model developed by Lau and Dill (Lau and Dill 1989). In this model, amino acids are of only two types: hydrophobic (H) or polar (P). The sequence is folded on a 2 D square lattice on which bonds are orthogonal to each other. Folded structures are evaluated by a so-called "hydrophobic potential" where each pair of nonbonded direct hydrophobic contact (occupying neighboring nondiagonal lattice vertices) receives −1. Using HP lattice models avoids the computational cost needed for all-atom models while still capturing the general principles that govern protein folding, and they can be extended to account for physicochemical characteristics of individual residues such as size, hydrophobicity and charge. In more detailed models, proteins can be represented as a sequence of pairs of dihedral angles that describe the backbone degrees of freedom of each residue. Mutations can be introduced simply by changing the dihedral angle of a residue and crossovers by swapping randomly assigned sections of two sequences (Schulze-Kremer 2000; Unger 2004).</p><!><p>An essential part of almost all successful protein-folding simulations or protein tertiary structure predictions is an energy function that is a good approximation to the energy landscape of real proteins. Energy functions can be roughly divided into two classes: physics-based force fields and knowledge-based potentials (Lazaridis and Karplus 2000). Historically, physics-based forces fields are coupled with MD or MC simulations to study protein dynamics or calculate free energies (Wang et al. 2001; Ponder and Case 2003; Mackerell 2004; Lopes et al. (2015), whereas knowledge-based potentials are mostly used for fold recognition or tertiary structure prediction (Sippl 1995; Godzik 1996; Skolnick 2006). Before we give a detailed account on them, we remind the reader that both of these two types of energy functions are evolving objects. To improve accuracy, further parameter optimization for physics-based force fields is required and statistics need to be rederived for knowledge-based potentials when energy function deficiencies are identified or data sets of better qualities become available.</p><!><p>Physics-based force fields are classical mechanical models that approximate the potential energy of chemical systems. Force field models ignore the electronic motions in a system and only consider interactions among nuclei. Compared to ab initio quantum mechanical methods, force fields are much more computationally efficient while giving an acceptable level of accuracy. A force field has a functional form and a (usually very large) set of associated parameters that, taken together, model bonded and nonbonded interactions in a system. The functional form of a force field is often a compromise between accuracy and computational efficiency and depends on the level of resolution (all atom or coarse grained), chemical nature (inorganic, small organic or biomolecular) and target properties of the systems to be modeled. Nevertheless, most force fields have five components (Equation (6)). The first three of them, so-called bond stretching, angle bending, and torsion, model bonded interactions. The last two components describe electrostatic and van der Waals nonbonded interactions (Leach 2001).</p><p>This functional form looks simple, but we must keep in mind that the set of parameters associated with it is very large. For example, the term that models bond stretching (a harmonic potential) has a different force constant kb and an equilibrium bond length l0 for each bond type. These parameters must be determined by fitting the force field to a given set of data obtained from experiments or quantum mechanical calculations. Depending on the size of the data set, parameter optimization may be conducted in a number of ways: trial and error, least-squares fitting (Lifson and Warshel 1968), or, recently, machine-learning algorithms (Behler (2016).</p><p>Well-known examples of force fields intended for modeling proteins include CHARMM (Gelin and Karplus 1979; Brooks et al. 1983; MacKerell et al. 1998, 2004; Brooks et al. 2009; Best et al. 2012), AMBER (Weiner and Kollman 1981; Weiner et al. 1984; Li and Bruschweiler 2010; Lindorff-Larsen et al. 2010), OPLS (Jorgensen and Tirado-Rives, 1988, Robertson et al. 2015), GROMOS (Van Gunsteren and Berendsen 1987; van Gunsteren et al. 1998), MARTINI (Marrink et al. 2007; Monticelli et al. 2008). These force fields were previously compared in-depth (Ponder and Case 2003), we note here that while the functional forms of these force fields invariably contain the five terms of Equation (6), some of them or their different versions may differ in specifics in the treatment of non-bonded interactions and the levels of resolution covered. For example, although more recent versions of the CHARMM and AMBER force fields do not model hydrogen-bonding energetics explicitly, originally CHARMM and AMBER force fields both incorporated a 12–10 Lennard-Jones potential to model hydrogen-bonding (Gelin and Karplus 1979; Weiner et al. 1984). The need for more efficient evaluation of nonbonded interactions arises when the number of interaction sites is large. One straightforward way to improve efficiency is to absorb aliphatic hydrogens into the carbon atom to which they are bonded to form "united atoms" as was done in the united-atom version of the CHARMM and OPLS force fields, or to use a coarse-graining approach where a group of heavy atoms are combined to form a representative virtual interaction site. The MARTINI force field aims at providing a simple model that is computationally fast and easy to use, and it adopted a "four-to-one" coarse-grain-ing scheme, meaning that on average four heavy atoms are represented by one interaction site (Marrink et al. 2007; Monticelli et al. 2008; Marrink and Tieleman (2013). Although all-atom simulations are often more desirable, if special care is taken during calibration of the building blocks and parameterization, a level of accuracy comparable to all-atom simulations may be possible in reproducing some thermodynamic properties with reduced representations while achieving considerable computational savings (Baron et al. 2006a, 2006b, 2007; Marrink et al. 2007; Monticelli et al. 2008; Marrink and Tieleman 2013). Coarse-grained protein models and their applications was recently reviewed in detail (Kmiecik et al. 2016).</p><p>Physics-based force fields are traditionally coupled with MD in simulating protein dynamics and folding (McCammon et al. 1977). There have been a plethora of such studies where the utility of force fields for protein tertiary structure prediction or the accuracy of reproducing experimental data were reported (Duan and Kollman 1998; Zagrovic et al. 2002; Pande et al. 2003; Summa and Levitt 2007; Lindorff-Larsen et al. 2011; Patapati and Glykos 2011; Lindorff-Larsen et al. 2012; Huang and MacKerell 2013; Piana et al. (2013). However, no agreement has been reached regarding whether force fields are sufficiently robust for these applications (Lee et al. 2009; Piana et al. (2014). Early analysis concluded that MD simulations under physics-based force fields are not particularly successful in structure prediction (Lee et al. 2009). However, for small, fast folding proteins that are also very stable, evidence has been accumulating that demonstrates that physics-based force fields are sufficiently accurate for predicting native-state structures and folding rates (Shaw et al. 2010; Lindorff-Larsen et al. 2011; Piana et al. 2012, 2013, 2014; Chung et al. (2015). In particular, it was pointed out the prediction of tertiary structures, folding rates and melting temperatures appears to be more robust than the prediction of the enthalpy and heat capacity of folding or that of the radii of gyration of unfolded states (Piana et al. 2014). It needs to be pointed out, however, that whether these force fields hold accurate for simulating larger proteins remains to be studied.</p><!><p>Unlike physics-based force fields, which model interactions found in the most basic molecular systems using fundamental laws of physics explicitly and separately, knowledge-based potentials (KBPs) are energy functions derived from statistical analyses of known protein structures and the application of the inverse Boltzmann relation to the probability distribution of geometries (Wodak 2002; Sippl 1993, 1995). The physical meaning of KBPs has been under vigorous debate since their introduction (Finkelstein et al. 1995; Thomas and Dill 1996; Ben-Naim 1997; Moult 1997; Shortle 2003; Hamelryck et al. 2010), although justifications of KBPs as "potentials of mean force" have been provided by analogy to the reversible work theorem in statistical thermodynamics (Sippl et al. 1996) or on the basis of probabilistic arguments (Simons et al. 1997; Hamelryck et al. 2010). Nevertheless, KBPs are widely used and surprisingly effective in scenarios including but not limited to protein structure prediction (Simons et al. 1997; Lu and Skolnick 2001; Shen and Sali 2006; Xu and Zhang (2012), refinement of NMR structures (Kuszewski et al. 1996; Yang et al. 2012), fold recognition (Kocher et al. 1994; Majek and Elber 2009), protein–ligand or protein–protein interactions (Gohlke et al. 2000; Zhang et al. 2005; Huang and Zou 2006a; Huang and Zou 2006b) and protein design (Poole and Ranganathan 2006). Thus, in this article, we summarize the formalism of KBPs, specific implementations of different types of potentials, and their applications instead of concerning about the physical interpretation of KBPs.</p><p>A KBP energy function is a linear combination of individual potentials with each capturing a specific type of interaction. The most common formulation of such energy functions is as follows:
(7)E(C∣S)=∑ijwij(−kTlnp(cj∣si)p(cj))
where E(C∣S) is the energy of conformation C given that the underlying amino acid sequence is S. p(cj∣si) is the probability that a given sequence si adopts conformation cj, whereas p(cj) is an unconditional probability that any sequence fragment adopts conformation cj· p(cj∣si)p(cj) can be thought of as an "equilibrium constant" of a hypothetical chemical reaction: random sequence, unique conformation → unique sequence, unique conformation (Shortle 2003). In addition to the above inverse Boltzmann formulation, other formulations of individual KBP terms have also been widely used. For example, the KBP under the modeling package Rosetta was formulated based on the Bayes' theorem (Simons et al. 1997). This approach was also adopted by Woetzel et al. recently to derive the KBP for a secondary structure element (SSE)-based protein structure prediction algorithm (Karakas et al. 2012; Woetzel et al. 2012; Weiner et al. 2013; Fischer et al. 2016). In their discrete optimized protein energy, or DOPE, Shen and Sali computed the negative logarithm of the joint probability density function of a given protein (Shen and Sali 2006).</p><p>The types of individual potentials incorporated into a KBP energy function are essentially only limited by the type of statistical relations that can be practically extracted from known protein structures. Depending on its intended purpose, a KBP may include individual potentials that fall into one or several categories. We elaborate three such potentials in the following and refer the reader to references (Simons et al. 1997; Woetzel et al. 2012; Xu and Zhang 2012) for examples of other potentials.</p><p>(1) Pairwise distance-dependent potential that approximates residue contact energies (Wodak 2002; Sippl, 1990, 1993, 1995). Such contact potentials are based on native inter-residue contacts that play a key role in determining folding kinetics and native state stability (Gromiha and Selvaraj 2004). The concept of pairwise distance-dependent potentials was first introduced in the pioneering work of Tanaka and Scheraga (Tanaka and Scheraga 1976), who related residue contact frequencies to the free energies of formation of corresponding interactions using the simple relationship between free energy and equilibrium constant. Their work was followed by that of Miyazawa and Jernigan (Miyazawa and Jernigan 1985, 1996), who formalized the theory of residue contact potentials using quasichemical approximation. However, these early implementations of contact potentials are not, in fact, distance-dependent, except that a single cutoff distance was used to define residue contact. A real pairwise distance-dependent potential was first introduced by Sipp (Sippl 1990), and this was followed by an explosion of different statistical potentials (Hendlich et al. 1990; Kocher et al. 1994; Park and Levitt 1996; Bahar and Jernigan 1997; Melo and Feytmans 1997; Park et al. 1997; Reva et al. 1997; Rooman and Gilis 1998; Samudrala and Moult 1998; Betancourt and Thirumalai 1999; Lu and Skolnick 2001; Zhou and Zhou 2002a; Fang and Shortle 2005; Qiu and Elber 2005; Summa et al. 2005; Dehouck et al. 2006; Shen and Sali 2006; Woetzel et al. 2012). Such pair potentials are usually formulated at residue level, where inter-residue distances are measured between Cβ atoms or side-chain centroids in reduced representation of amino acid residues to promote computational efficiency. However, atomic-level formulation usually gives better discriminatory power albeit at the cost of more computational resource (Sippl 1996; Sippl et al. 1996; Melo and Feytmans 1997; Samudrala and Moult 1998; Lu and Skolnick 2001; Shen and Sali 2006).</p><p>(2) Solvent accessibility-based environment potentials that represent the interactions of individual residues with their local environment (Bowie et al. 1991; Kocher et al. 1994; DeLuca et al. 2011; Xu and Zhang 2012). Residue environment potentials are often included to account for solvation effects. Precise calculation of solvent accessibility requires full atomic structure and is time-consuming. In tertiary structure prediction scenarios where reduced representations of residues are used, good approximations to solvent accessibility, such as residue contact numbers, provide significant computational savings (Durham et al. 2009; Woetzel et al. 2012; Fischer et al. 2015; Li et al. 2016). It should be noted that in addition to transforming solvent accessibility statistics to energy-like potentials using the inverse Boltzmann relation, they have also been incorporated into KBP energy functions as a penalty term to disfavor models where residue-specific solvent accessibilities disagree with expected solvent accessibilities (Xu and Zhang 2012, Li et al. 2017).</p><p>(3) potentials of torsion angles that evaluate backbone φ, ψ torsion angles and/or the preference of side-chain rotamers (Kocher et al. 1994, Kuszewski et al. 1996, Betancourt and Skolnick 2004; Fang and Shortle 2005; Amir et al. 2008; Yang et al. 2012; Kim et al. 2013). It is well known that only certain combinations of φ, ψ torsion angles are populated in proteins (Ramakrishnan and Ramachandran 1965) and significant correlations exist between side-chain torsion angle probabilities and backbone φ, ψ angles (Dunbrack and Karplus 1993). Including such potentials has been shown to enable the energy function to exclude conformations that have unlikely combinations of torsion angles. In a study by Kocher et al. (Kocher et al. 1994) where several types of potentials were tested to recognize protein native folds, potentials representing backbone torsion angle preferences recognized as many as 68 protein chains out of a total of 74. This result was striking given the fact that backbone torsion potentials consider solely local interactions along the chain and are well known to be incapable of determining the full 3 D fold (Kocher et al. 1994). Potentials of torsion angles have also been used to refine structures generated from NMR data (Kuszewski et al. 1996, Yang et al. 2012). Kuszewski et al. (1996) incorporated a database-derived torsion angle potential into the target function for NMR structure refinement, resulting in a significant improvement in various quantitative measures of quality (Ramachandran plot, side-chain torsion angles, and overall packing. In a similar way, Yang et al. (2012) constructed a database of 2405 refined NMR structures.</p><!><p>Due to their intrinsic inaccuracies, a common issue with energy functions is that incorrect conformations may be scored comparably to (or even better than) the native state (Skolnick 2006), lending the energy function inability to recognize the native state (Figure 6(A)). This issue be remedied by incorporating sparse experimental data as restraints, which offers some structural information that by itself is insufficient to completely determine the protein's structure (Figure 6 (B,C)).</p><!><p>Restraints from sparse experimental data drastically decrease the conformational space that needs to be sampled to only those structures consistent with the data. Many software suites implement algorithms to couple their de novo prediction methods with limited experimental data, including those from nuclear magnetic resonance (NMR), electron paramagnetic resonance (EPR), cross-linking mass spectrometry (XL-MS) and electron microscopy (EM).</p><p>NMR rivals X-ray crystallography as a technique by which an entire protein structure can be unambiguously determined. Solution-state NMR can determine the structure of relatively small proteins (< ~20kDa), but intensive experimental techniques and analysis of NMR spectra are required to determine a high-quality structure of a protein. Each residue typically requires upwards of 15 constraints. Oftentimes, NMR spectroscopy can provide some degree of low-resolution information about the global conformation of a protein, even for larger proteins (Venters et al. 1995; Battiste and Wagner 2000). These sparse restraints, including chemical shifts (CSs), nuclear overhauser enhancements (NOEs) and residual dipolar couplings (RDCs) do not provide enough information to fully determine the structure of a protein, but they can be used in conjunction with computational protein structure prediction software. CSs provide information about the protein backbone conformation, while NOEs and RDCs give information about the global fold of the protein. De novo protein structure prediction software can take advantage of just CSs (Latek et al. 2007), CSs and NOEs (Bowers et al. 2000) or all three types of restraints (Weiner et al. 2014).</p><p>Site-directed spin labeling (SDSL) and EPR can be used to glean information about proteins of nearly any size in their native environments. In addition, only a small amount of sample is required for structural interrogation by EPR. The accessibility and mobility of the spin labels can be used to determine the exposure and topology of SSEs (Farahbakhsh et al. 1992; Altenbach et al. 2005). Distances between spin labels can be detected up to 60Å and can give insight into the overall fold of the protein as well as different conformational states (Rabenstein and Shin 1995; Borbat et al. 2002). However, it is not feasible to use EPR to determine the full structure of a protein. EPR is experimentally intensive, as it requires the introduction of unpaired electrons at selected sites within proteins. This is usually done by cysteine substitution mutagenesis followed by modification of the sulfhydryl group with a nitroxide reagent. However, nonsense suppressor methodology, solid-phase peptide synthesis or "click-chemistry" have also been used (Klare and Steinhoff 2009). This technique will only give a small part of structural information about the protein, so these sparse EPR data can be used in conjunction with computational protein structure prediction methods (Alexander et al. 2008; Hirst et al. 2011; Fischer et al. 2015). The selection of sites to spin label is integral to the efficacy of structure determination by EPR (Alexander et al. 2008).</p><p>Similarly, XL-MS experiments can be used to determine interatomic distances that serve as experimental restraints. XL-MS can be used with proteins in their native states, and it has proven to be compatible with relatively large proteins, flexible proteins and membrane proteins (Kalkhof et al. 2005, Jacobsen et al. 2006, Lasker et al. 2012). In addition, the samples used can be heterogeneous and dynamic, as the output of XL-MS experiments is an average. The basis of XL-MS is the ability of two functional groups of a protein to form covalent bonds if they are within a certain distance of one another. These cross links can occur both inter- and intramolecularly. The proteins are then enzymatically digested, and MS is used to identify these cross links and surface labels (Young et al. 2000; Back et al. 2003; Sinz 2003).</p><p>EM provides data similar in format to that of X-ray crystallography, that is, a density map of a protein or complex. The data are thus less sparse than many of the aforementioned experimental techniques, but EM has historically provided lower-resolution density maps, from which an atomic structure cannot be gleaned. However, even low-resolution EM density maps are integral for identifying the overall organization of large molecular complexes. In recent years, EM technologies have progressed such that density maps with resolutions in the range of 4–8 Å can regularly be attained, at which level SSEs can be visualized and even some side-chain character can be visualized (Bihnstein and Melanie 2015). Many computational modeling methods have been developed that work with EM density maps (Lindert et al. 2009b), either in fitting previously solved structures into density maps, determining the topology and location of SSEs (Jiang et al. 2001; Abeysinghe et al. 2008), performing comparative modeling and de novo protein structure prediction (Lindert et al. 2009a, 2009c, 2012a, 2012b, 2012c; Woetzel et al. 2011).</p><p>Most de novo protein structure prediction algorithms require the use of a segmented density map, which can be accomplished with the use of various segmentation algorithms (Baker et al. 2006; Pintilie et al. 2010; Burger et al. 2011). Then, SSEs can be extracted from the density map either manually or with the use of algorithms that automate the selection of helices and/or sheets from a segmented density map (Jiang et al. 2001, Kong and Ma 2003; Kong et al. 2004; Baker et al. 2007). Next, de novo modeling algorithms can use these data with the density map and primary sequence of the protein in order to create a full structural model either via optimization (Chen et al. 2016) or using Monte Carlo methods (Lindert et al. 2012b, 2009c; Wang et al. 2015).</p><!><p>If no experimental restraints are available for the protein, secondary and tertiary structural restraints can be predicted from an amino acid sequence based on existing structures. Secondary structures can be predicted using machine learning methods. Artificial neural networks (ANNs) can be used to predict secondary structures from position-specific scoring matrices (Jones 1999; Yan et al. 2013), reduced amino acid representation (Leman et al. 2013), or multiple sequence alignments (MSAs) (Rost and Sander 1993a; Rost et al. 1993b). Methods have also been developed specifically to predict membrane protein topology from amino acid sequence using ANNs (Viklund et al. 2008; Viklund and Elofsson 2008; Leman et al. 2013), support vector machines (SVMs) (Nugent and Jones 2009), or Hidden Markov Models (HMMs) (Krogh et al. 2001; Kahsay et al. 2005).</p><p>It is a long-standing observation that 3D protein folds can be predicted from sufficient information regarding the protein's inter-residue contacts (Gobel et al. 1994; Olmea and Valencia 1997; de Juan et al. (2013); the addition of even relatively sparse information about tertiary contacts into an algorithm's scoring function can help improve protein models (Kim et al. 2014). Recently, the incorporation of long range contact predictions has resulted in some of the the most effective de novo protein structure prediction algorithms (Monastyrskyy et al. 2016; Moult et al. 2016). Several algorithms have been devised to predict these contacts using the principle of correlated mutations (de Juan et al. 2013). In general, amino acid contacts that stabilize the protein fold are assumed to evolve complementarily – if one residue of a contact is mutated, the other will likely also mutate to a reasonable interaction partner.</p><p>In order to identify pairs of correlated mutations, amino acid pairs can be scored based on their physicochemical similarity using the McLachlan matrix (McLachlan 1971), which is based on the frequencies of observed mutations in homologous proteins. Correlated mutations can also be scored by mutual information between MSAs based on the equation
(8)MI=∑abf(aibj)logf(aibj)f(ai)f(bj)
</p><p>The above equation indicates that the mutual information between two protein sites i and j is computed by summing over amino acid pairs ab for every amino acid type a and b, where f(aibj) is the observed relative frequency of ab at columns ij and f(ai) is the observed relative frequency of amino acid type a at position i. The identification of these correlated mutations is used in many methods of multiple sequence alignment (Göbel et al. 1994; Neher 1994; Pollock and Taylor 1997; Ashkenazy and Kliger 2010; Hopf et al. 2014), from which tertiary contact predictions can be extrapolated.</p><p>In recent years, numerous algorithms have come out that account for covariance caused by indirect inter-residue coupling effects, which has led to improvement in prediction of correlated mutations (Burger and van Nimwegen 2010; Marks et al. 2011, 2012; Morcos et al. 2011; Jones et al. 2012b; Kamisetty et al. 2013; Skwark et al. 2013; Ekeberg et al. 2014; Hopf et al. 2014; Kaján et al. 2014; Michel et al. 2014; Ovchinnikov et al. 2014; Skwark et al. 2014; Jones et al. 2015). These methods were developed to resolve the issue that two residues aligned in multiple sequence alignments may exhibit statistical dependencies even though they are distant in physical space, which usually arises from chains of interacting pairs of residues. Also, information regarding the conservation of certain residues regardless of their tertiary contacts must be considered for correlated mutations to properly represent actual 3 D contacts. Many methods have been devised that decouple direct from indirect residue coevolution, primarily based on statistical methods. Covariation-based contact prediction has also proven successful as a scoring metric for de novo folding (Morcos et al. 2011; Kamisetty et al. 2013).</p><p>Machine learning methods, including ANNs (Fariselli and Casadio 1999; Fariselli et al. 2001; Shackelford and Karplus 2007; Tegge et al. 2009; Xue et al. 2009), genetic algorithms (MacCallum 2004; Chen and Li 2010), ran-dom forests (Li et al. 2011), HMMs (Björkholm et al. 2009; Lippi and Frasconi 2009) and SVMs (Cheng and Baldi 2007; Wu and Zhang 2008), have also arisen as successful methods to predict 3 D contacts. These methods use various features to predict contact maps. Some of the most successful of these machine learning methods for contact prediction are hybrid methods that predict contacts based on both physicochemical features and evolutionary features, using MSAs as part of their training data sets (Wallner and Elofsson 2006; Stout et al. 2008; Ma et al. 2013; Kosciolek and Jones 2015).</p><!><p>Protein structure prediction methods can be broadly grouped into template-based modeling, where construction of target models involves threading the target sequence through the structure of homologous proteins (templates) and de novo structure prediction, where target models are constructed from sequence alone, without relying on similarity at fold level between the target sequence and any of the known structures (Baker and Sali 2001; Bonneau and Baker 2001; Hardin et al. 2002; Lee et al. 2009). Template-based modeling is based on the premise that tertiary structures of proteins in the same family are more conserved than their primary sequences (Chothia and Lesk 1986; Fiser et al. 2002; Illergard et al. 2009). While it can produce accurate models for target sequences if templates with sequence identity >25% are used (Cavasotto and Phatak 2009) and can be practically useful (Xiong et al. 2011; Zhan et al. 2011; Li et al. 2012), it is nevertheless purely mechanical in that it does not provide a general understanding of the role of particular interactions in maintaining the stability of protein structure (Baker and Sali 2001; Cavasotto and Phatak 2009). Thus, one could not gain insights into the physicochemical principles underlying protein folding (Pillardy et al. 2001; Lee et al. 2009). On the contrary, de novo methods sample and energy – evaluate the folded conformations as thoroughly as computational resource permits, and they assume the native conformation is the one with the lowest energy. Logically, two of the most crucial factors that dictate whether a de novo tertiary structure prediction method will be successful are its coverage of the conformational space and how accurate its energy function is. In this section, we discuss in detail some selected examples of de novo tertiary structure prediction methods and highlight some successful cases from the history of CASP (Figure 7). Note that this selected set of methods is by no means exhaustive. The interested reader is referred to proceedings of CASP experiments (http://predictioncenter.org/ index.cgi?page= proceedings), which cover a wider spectrum of methods and in more detail.</p><!><p>FRAGFOLD was developed based on the rationale that proteins tend to have common structural motifs at the super-secondary structural level (Jones 1997b, 2001; Jones and McGuffin, 2003). In FRAGFOLD, 3 D models are built by assembling super-secondary structural fragments from a library of highly resolved protein structures with MC-simulated annealing and evaluated with a knowledge-based energy function. FRAGFOLD was initially tested in CASP2 (Jones 1997b), and later in CASP4 (Jones 2001) and CASP5 (Jones and McGuffin 2003). Its success in predicting the fold of NK-Lysin marked the first correct de novo blind prediction of a protein's fold (Jones 1997c).</p><p>The super-secondary structural fragments considered by FRAGFOLD include α-hairpin, α-corner, β-hairpin, β-corner, β-α-β unit, and split β-α-β unit. Favorable super-secondary structural fragments are selected based on the quality of threading. Threads that contradict the reliable regions of predicted secondary structure by PSIPRED (Jones 1999) are skipped. In addition to this sequence-specific fragment list, a general fragment list that consists of all tripeptide, tetrapeptide and pentapeptide fragments is also constructed from a library of highly resolved protein structures. The knowledge-based energy function in FRAGFOLD initially consists of a set of pairwise potentials, a solvation potential, a term for penalizing noncompact folds, a term for penalizing steric clashes and a term that accounts for hydrogen bonding (Jones 1997b). This energy function was recently complemented with predicted contacts as restraints (Kosciolek and Jones 2014). Kosciolek and coworkers (Kosciolek and Jones 2014) found that combining statistical potentials with contacts predicted by PSICOV (Jones et al. 2012a) is significantly better than either statistical potentials or predicted contacts alone.</p><!><p>The Rosetta algorithm for de novo protein structure prediction employs MC-simulated annealing to assemble protein-like 3 D models from fragments of unrelated protein structures with similar local sequences using an energy function based on Bayes' theorem (Simons et al. 1997; Rohl et al. 2004). The algorithm is based on the experimental observation that local sequence preferences bias, but do not uniquely determine, the local structure of a protein (Rohl et al. 2004). Rosetta has turned out to be one of the most successful methods indicated by results from CASP experiments (Bradley et al. 2003, 2005a; Jauch et al. 2007) and several other studies (Bradley et al. 2005b; Ovchinnikov et al. 2017) (Figure 7(A) for an example).</p><p>Model construction in Rosetta is performed via a sequence of fundamental conformation modification operations termed "fragment insertion". For each fragment insertion, a sequence segment of three or nine residues is selected, and the torsion angles of these residues are replaced with the torsion angles of a homologous fragment selected from a ranked list of fragments of known structure (Simons et al. 1997). Fragment insertions that decrease the energy of the resulting conformation are accepted and those that increase the energy are accepted according to the Metropolis criterion (Metropolis et al. 1953). Derivation of the Rosetta energy function was based on a Bayesian separation of the total energy into components that describe the like-lihood of a particular structure, independent of sequence and those that describe the fitness of the sequence given a particular structure (Simons et al. 1997).</p><p>The original Rosetta energy function is coarse grained: terms corresponding to solvation and electrostatic effects are based on observed residue distributions derived from known protein structure databases, and hydrogen bonding is not explicitly described. However, preferences of β-strand pairing geometries and β-sheet patterns are included. Steric clashes are penalized, while van der Waals interactions are not explicitly modeled. A more physically realistic, atomic-level energy function was developed later for applications requiring more detailed structurally information. In this "fine-grained" version of the energy function, van der Waals interactions are modeled with a 6–12 Lennard–Jones potential. Solvation effects are included, using the Lazaridis–Karplus model (Lazaridis and Karplus 1999), and hydrogen–bonding is explicitly accounted for using a secondary structure- and orientation-dependent potential derived from high-resolution protein structures (Kortemme et al. 2003). Energetics of local interactions are described using an amino acid- and secondary structure-dependent potential for backbone torsion angles. The reader is referred to reference (Rohl et al. 2004) for a more mathematically detailed description of the Rosetta energy function.</p><!><p>Recent CASP experiments have shown significant advantages of integrating various techniques such as threading, de novo modeling and atomic-level structure refinement approaches into a single pipeline of tertiary structure prediction (Battey et al. 2007; Jauch et al. 2007; Zhang 2009; Kinch et al. 2011, 2016; Tai et al. (2014). The I-TASSER method, (Wu et al. 2007; Roy et al. 2010; Yang et al. 2015) which implements TASSER (Zhang and Skolnick 2004) in an iterative mode, is one example of the composite approaches. I-TASSER has been particularly successful as shown by recent CASP experiments (Zhang 2009; Roy et al. 2010; Yang et al. 2015; Zhang et al. 2016) (Figure 7(B) for an example).</p><p>I-TASSER uses a sophisticated threading scheme, which compares the target sequence with template structures using profile–profile alignment, for the selection of the most probable structure fragments. Aligned regions of the target sequence are modeled by connecting template fragments through a random walk of Cα-Cα bond vectors of variable lengths. Unaligned regions are simulated on a cubic lattice system for computational efficiency. Initial full-length coarse-grained models are refined via REMC simulation where two kinds of moves are implemented: off-lattice rigid fragment translations and rotations of the aligned regions and on-lattice 2–6 bond movements and multibond sequence shifts of unaligned regions (Zhang and Skolnick 2004). The models of the first-round TASSER simulation are clustered and the cluster centroids are submitted to a second-round TASSER simulation to remove physically unrealistic interactions. Finally, back-bone atoms and side-chain rotamers are added to the model with the lowest energy from the second round (Wu et al. 2007). The energy function of I-TASSER includes the original TASSER knowledge-based potential and a new burial potential based on neural network-predicted accessible surface area (ASA) (Wu et al. 2007). The original TASSER potential consists of long-range pair interactions of side-chain centers of mass, local Cα correlations, hydrogen-bond, hydrophobic burial interactions, propensities for predicted secondary structures, protein-specific pair potentials of side-chain centers of mass and tertiary contact restraints extracted from the threading templates (Zhang et al. 2003).</p><!><p>QUARK is an algorithm for denovo protein structure prediction using REMC simulations guided by a consensus knowledge-based energy function. In contrast with Rosetta and I-TASSER that assemble fragments of fixed sizes, QUARK assembles 3 D models from small structure fragments of multiple sizes from 1 to 20 residues. To increase the structural flexibility and the efficiency of conformational search, QUARK also implements a set of MC moves consisting of free-chain constructions and fragment substitutions between decoy and fragment structures (Xu and Zhang 2012). The QUARK algorithm has been shown to be highly successful in recent CASP experiments (Xu and Zhang 2012; Zhang et al. 2016) (Figure 7(C) for an example).</p><p>QUARK generates structure fragments for target sequences by threading sequence segments through a library of nonhomologous experimental structures. Multiple features such as solvent accessibility, real-value φ and ψ angles, and secondary structure types as predicted from back-propagation neural networks are used to improve generation of structure fragments. Optimization of 3D models is performed via REMC simulations that start with initial models assembled by chaining randomly selected fragments with varied sizes. Conformational sampling of each replica is done through residue-level, segment-level and topology-level movements. After each running cycle, the conformations between every two adjacent replicas are exchanged according to the Metropolis criterion (Metropolis et al. 1953). Protein structure models built by QUARK are evaluated by a composite knowledge-based energy function consisting of atomic-level pair potentials, hydrogenbonding potential, SSE packing potentials, heuristic terms that account for excluded volumes, solvent accessibility, and radius of gyration (Xu and Zhang 2012).</p><!><p>The BCL::Fold algorithm developed in our group seeks to overcome the limitations of protein size and fold complexity by assembling idealized SSEs (secondary structure elements) into 3 D models. This algorithm was developed under the framework model of protein folding. As discussed previously, while the framework model is not always true, it is straightforward to implement. In addition, as shown by our benchmark study (Karakas et al. 2012; Weiner et al. 2013), BCL::Fold facilitates the sampling of nonlocal contacts. Thus, BCL::Fold may be a promising tool for the structure prediction of proteins with high contact order (Plaxco et al. 1998; Baker 2000; Bonneau et al. 2002). It is also worth mentioning that in contrast to the other four methods, which heavily reply on the availability of homologous template structural fragments (short or long), BCL::Fold is "truly" de novo in the sense that no template structure is needed at any stage of the algorithm. While BCL::Fold was not ranked among the most successful methods, we would still like to highlight the CASP11 target T0769. While this protein is in the category of template-based modeling, meaning that a suitable template can be identified that covers all or nearly all residues of the target, BCL::Fold predicted a model with a Cα-RMSD of 1.8 Å to the released solution NMR structure without relying on any homologous templates (Fischer et al. 2016) (Figure 7(D)).</p><p>In BCL::Fold, the necessary complexity reduction of the conformational space is achieved by assembling SSEs from a predetermined pool of SSEs using MC-simulated annealing and omitting more flexible loop regions. A high-quality pool of SSEs can be readily created using machine learning-based secondary structure prediction methods such as PSIPRED (Jones 1999). BCL::Fold implements a comprehensive list of SSE-based MC moves, which are categorized into six main categories: adding SSEs, removing SSEs, swapping SSEs, single SSE moves, SSE pair moves and moving domains consisting of multiple SSEs (Karakas et al. 2012). Models generated by BCL::Fold are evaluated by a knowledge-based consensus energy function called BCL::Score (Woetzel et al. 2012), which consists of potentials of residue pair interaction, residue environment, SSE packing, β-strand pairing, loop length, radius of gyration, contact order, secondary structure prediction agreement. Separate penalizing energy terms were also included to exclude conformations with clashes between amino acids or SSEs and loops that cannot be closed (Woetzel et al. 2012). BCL::Score can also be complemented with experimental or predicted restraints to improve selection of native-like models (Weiner et al. 2014; Fischer et al. 2015; Li et al. 2017).</p><!><p>In the past decade, we have seen hardware and algorithmic advances that enabled researchers to perform millisecond timescale simulations of protein folding, and we have also seen development of methodologies that predicted tertiary structure with better accuracy for proteins with larger size. Despite these achievements, there is still a long list of challenges on the way toward a solution to the protein folding problem.</p><p>On the folding mechanism side, even though long simulations have been available, unambiguous scientific results learned from such simulations have thus far been modest (Lane et al. 2013). First, it is still being debated whether proteins fold via a single definite pathway or multiple parallel pathways. Although both views have received support from simulations and experiments (Englander and Mayne 2014; Wolynes (2015), additional simulations with more robust trajectory analysis and experimental validation are required to disambiguate conflicting results. Second, realistic folding simulations have thus far been limited to small proteins (<100 residues), it is questionable whether folding mechanisms revealed by these simulations are generalizable to larger proteins. Thus, simulating the folding of larger proteins will likely be a major trend for the next decade. Finally, as far as we are aware, a theory that is quantitative and makes specific prediction about how a protein would fold is not yet available. The some-what loosely defined models of hierarchical (framework) folding, nucleation-condensation and foldons are difficult to validate or invalidate either by experiments or by simulations. Nevertheless, closer interaction between simulations and experiments such that simulations be tested by experiments and in turn aid in the interpretation of experimental results and guide the design of future experiments will have greater impact on the field.</p><p>On the structure prediction side, larger proteins, especially those with multidomains, stay a significant challenge to de novo structure prediction methodologies. These proteins are often characterized by their high contact order and long folding time (Plaxco et al. 1998; Paci et al. 2005). Conformational sampling of these proteins is usually inefficient and is complicated not only by protein size, but also by the considerable number of non-local contacts, which are formed by residues far apart in sequence but usually critical for structural stability (Moult 2005; Kim et al. 2009). Consequently, tools for de novo structure prediction are not likely to become practically useful for structure prediction for any but very small, sometimes medium-sized proteins (Jones 1997a; Baker and Sali 2001). Other challenging targets, especially for methods whose energy functions heavily rely on statistics extracted from known structures, may also include proteins with rare and unusual folds (Kinch et al. 2016). Accurate prediction of tertiary structure for these challenging targets certainly requires the joined forces of high-performance hardware, efficient algorithms for conformational sampling, accurate energy functions, and, last but not least, valuable experimental restraints.</p>