As second half of the genetic code, the prediction of tertiary structure of proteins from their primary amino acid sequence is one of the most important and challenging problems in contemporary structural biology. There are three classes of theoretical approaches to the problem in the recent literature (Murzin, 2001): homology modeling (Guex and Peitsch, 1997; Sanchez and Sali, 1997); threading (Bowie et al., 1991; Panchenko et al., 2000; Skolnick and Kihara, 2001); and ab initio folding (Pillardy et al., 2001; Simons et al., 2001; Kolinski and Skolnick, 1998). Although homology modeling aims to find a template protein whose sequence is clearly evolutionarily related to the query sequence, the aim of threading is to detect both evolutionary-related sequences and analogous folds, which adopt very similar structures to the query protein. Both threading and homology modeling, in principle, are capable of producing high-resolution folds based on the identified template proteins, but they suffer from the fundamental limitation that the native topology for the sequence of interest must have already been solved; and new folds cannot be predicted by these approaches. To address this issue, the most difficult and general approach is ab initio folding, where one attempts to fold a protein from a random conformation.

In principle, ab initio approaches are based on the thermodynamic hypothesis formulated by Anfinsen (Anfinsen, 1973), according to which the native structure corresponds to the global free energy minimum under the given set of conditions. The success of the approach therefore relies on the effectiveness of the following factors: 1), An implicit representation of the protein with sufficient structural fidelity and computational tractability. 2), A force field having near-native structures as its global minimum. 3), A protocol to effectively search the important regions of conformational phase space in a reasonable amount of CPU time. 4), A methodology to correctly identify near-native structure from the decoys produced by the simulation.

In this article, we will present our efforts to address all of these four issues that extend and improve our previous TOUCHSTONE approach (Kihara et al., 2001). Here, we exploit a new lattice representation for the protein structure, in which three united atom groups of C_α, C_β, and the side-group (SG) center of mass of the remaining (non-C_β) heavy atoms (CABS) are specified. Compared with our previous side-chain-only (SICHO) model (Kolinski and Skolnick, 1998) where only the side-chain centers of mass are treated, the CABS model has higher geometric fidelity.

Similar to the SICHO model (Kolinski and Skolnick, 1998), the basic force field includes energy terms describing short-range structural correlations, hydrogen-bond interactions, long-range pairwise potentials, one-body burial interactions, and a residue-contact-based environmental profile. All interactions are reconstructed in a more specific way in the new lattice model. For example, the H-bond and proteinlike conformational stiffness are more precisely constructed because of the inclusion of explicit C_α atoms in the model. The combination of C_α-C_α and SG-SG correlations provides for short-range interactions of higher amino acid specificity than the SICHO model. We also incorporate electrostatic interactions for the charged residues and a global propensity to the predicted contact order and contact number. Because these energy terms are not independent, some interactions are overcounted. To combine all of these energies, we create 60,000 decoys for each of 30 training proteins of diverse lengths (47 ~ 146 residues) and topologies. We obtain their weight factors by maximizing the correlation between the total energy and the structural similarity of decoys to the native structure, and by maximizing the energy gap between native structure and decoy ensembles. Decoy-based optimization of force fields has been exploited in previous studies that either maximize the correlation of the energy (scoring) function and RMSD to native (Simons et al., 1999) or require a lower energy of the native structure than the ensemble of decoys (Vendruscolo et al., 1999; Tobi and Elber, 2000). Here, we find that their combination provides for a better folding yield than when using either one alone. The optimized force field has a significantly improved energy versus RMSD correlation in favor of the native structure, compared to the naïve uniformly weighted combination of all the energy terms.

To effectively search the resultant energy landscape, we exploit the recently developed parallel hyperbolic sampling algorithm in our Monte Carlo (MC) simulations (Zhang et al., 2002). Previously, this protocol was shown to be more effective than general replica sampling in searching for low-energy structures, especially for proteins of large size where the energy landscape is significantly more rugged than the energy landscape of small proteins. To identify near-native structures from decoys generated in the MC simulations, we exploit the structure-clustering algorithm (SCAR) (Betancourt and Skolnick, 2001) to cluster the low-energy trajectories. We introduce two quantities dependent on the cluster density and the combination Y of energy and free energy, which are more discriminative than the generally used average energy and cluster size for the identification of near-native structures.

We apply our approach to a test set of 125 proteins (65 proteins that are the same as used in the original TOUCHSTONE paper (Kihara et al., 2001) plus an additional, harder 60-protein test set that covers a larger range of protein sizes). Using only protein sequence information, we can fold 41 cases that have structures with root mean square deviation (RMSD) from native of 1.79 ~ 6.5 Å in the top five clusters. All these foldable cases are restricted to small proteins (36 ~ 120 residues). To fold proteins of larger size and to improve the folding yield of the small proteins as well, we take the threading-based predictions of side-chain contacts as loose restraints in our force field to guide the folding simulations. These restraints are collected from consensus contacts hit by PROSPECTOR (Skolnick and Kihara, 2001). Their inclusion results in a significant improvement in the overall folding performance. There are 83 cases (70 cases with length less than 120 residues plus 13 cases with 120 ~ 174 residues) in the restraint-guided simulations that have at least one structure with a RMSD to native below 6.5 Å in the top five clusters. Especially, for the 60 harder representative proteins, the fraction of foldable cases (defined as having one of the top five clusters with RMSD from native below 6.5 Å) by the SICHO and CABS models are 1/3 and 1/2, respectively (i.e., 20/60 and 32/60), indicating a qualitative improvement of the new CABS model over the SICHO model. This improvement may, however, be partly due to the force-field optimization procedure used for the CABS model.

This article is organized as follows: we first describe the lattice representation of protein structure. Second, we give a detailed discussion of the interaction scheme and the procedure used to optimize the force field. This is followed by a description of the conformational search engine and the secondary structure prediction scheme. Then, we present the results of our approach applied to representative proteins, and our method for the evaluation of the simulation results. Finally, we summarize the key results.

FIGURE 1 Schematic representation of a three-residue fragment of polypeptide chain in the CABS model. The Cα trace is confined to the underlying cubic lattice system, whereas the Cβ atom and side-group rotamers are off-latticed and specified by (more ...)

TABLE 1 The lattice vectors employed in the representation of the main chain Cα trace

The positions of three consecutive C_αs define the local coordinate system used for the determination of the remaining two interaction units: the β-carbon (C_β) (except glycine), and the center of mass of remaining side-group heavy atoms (except glycine and alanine). A two-rotamer approximation has been assumed, depending on whether the configuration of the main chain is expanded (for instance in a β-sheet) or compact (for instance in an α-helix). The secondary structure-dependent numerical parameters for the determination of the C_β and SG positions are extracted from the protein data bank (PDB) (Berman et al., 2000).

The excluded volume of the envelope of the C_α and C_β atoms are represented as identical size hard spheres (infinite energy of overlap) of diameter 3.25 Å plus a 1/r type of soft-core potential in the range [3.25 Å, 5.0 Å]. This mimics the minimal observed cutoff distance of 4.0 Å in real proteins, and allows a few atoms to approach closer than the reality at a penalty, thereby partly remedying the coarseness of the discrete lattice model. The excluded volume of the SG units is approximated by a strong energy penalty when the distance of a side group from other units is below cutoff values specific to the interacting pair of amino acids. With the above geometric restrictions, all PDB structures can be represented with an average RMSD of 0.4 Å from native, better than that of 0.8 Å via the SICHO model (Kolinski and Skolnick, 1998). This geometric fidelity does not show any systematic dependence on protein length.

The different energy terms achieve various effects on the generation of nativelike states. The most important factors for overall folding in our force field are the secondary structure prediction propensities, hydrogen bonding, and tertiary contact restraints derived from threading. The first two terms provide a basic folding framework; the contact restraints are of critical importance in modifying the energy landscape to guide the simulations to near native states, especially for large proteins where the general proteinlike potential cannot distinguish the native state among a huge number of possible topologies.

The other terms describing short-range correlations, environment profiles, burial, and long-range pairwise interactions are helpful for refining the packing of side chains and local fragments; the bias toward predicted contact order and contact number also helps somewhat to speed up the folding processes.

Some of the above energy terms are similar to that in the SICHO model (Kolinski and Skolnick, 1998); however, the implementation is different in the new lattice model. For the sake of completeness, we present all the energy terms used in the CABS model. Below, we first describe the short-range interactions and then the long-range terms. Next, we will determine the relative weights of the energy terms in the combined force field, based on the correlation between the energy and structure quality of the decoys.

Short-range interactions

Multiple short-range correlations The potential contains both C_α-C_α and SG-SG local structure correlations derived from the PDB as the negative logarithm of the relative frequency histogram:

(1)

Here, A_i denotes the amino acid identity of the ith residue; r_i,j (s_i,j) is the C_α (SG) distance between the ith residue and the jth residue; _i denotes the local chain chirality of three consecutive C_α-C_α vectors from i to i + 3. E_i,j represents the C_α-C_α correlation of the ith and jth residues extracted from a statistical analysis of a structural data base of nonhomologous proteins (Kolinski and Skolnick, 1994, 1998). E₁₃ includes only two bins depending on the distance of r_i,i+2, which correspond to local extended and compact structures, respectively. E₁₄ and E₁₅ include more bins because more distant interactions are involved. When the predicted secondary structure is assigned in the fragments, both E₁₄ and E₁₅ are half/half combinations of the general and secondary structure specific parts extracted from generic and secondary structure specific fragments of the PDB, respectively.

E′_ij in Eq. 1 represents the local side-group correlation from the ith residue to jth residue and is derived from a set of PDB structures that have certain levels of sequence similarity to the query proteins and where homologous proteins of more than 25% sequence identity to the query proteins are excluded from the structural data base (Kolinski et al., 1998). When predicted secondary structure is assigned, E′_ij is also combined with a SG-SG correlation potential derived from secondary structure specific fragments of PDB data base. The w_i values in Eq. 1 and the equations below are the relative scale factors of these interactions that will be determined in the next section.

Local conformational stiffness This potential describes the characteristic local stiffness of global proteins and the general tendency toward regular arrangements of (predicted and nonpredicted) secondary structure:

(2)

Here, the unit tangent vector , and bisector vector , where r_i,i+1 is the C_α-C_α bond vector from vertex i to vertex i + 1. The first two terms represent the general propensities to common bond-vector orientations of α-helical and β-sheet structures as shown in Fig. 2. The third term is designed to impose further structure biases to the individual regularities of α-helical and β-sheet structures and is written as

(3)

FIGURE 2 Schematic illustration of the virtual Cα-Cα vectors for regular helical and sheet structures. , , , where ri,i+1 is the Cα-Cα bond vector from vertex i to vertex i + 1. As demonstrated in the first two terms (more ...)

Θ₂ denotes the strong tendency to form predicted secondary structures, which are taken from the combined PSIPRED (Jones, 1999) and SAM-T99 (Karplus et al., 1998) secondary structure prediction algorithms (discussed in “Secondary structure prediction”).

(4)

Θ₃ imposes a penalty to the irregular crumpled structures, i.e.

(5)

Hydrogen bonds Hydrogen bond interactions can be short range or long range depending on the secondary structures of the involved residues, although we list it here in the short-range category of interactions. Only main chain hydrogen bonds are considered. Due to the lack of the explicit positions for the peptide bond atoms, the effect of hydrogen bonds is translated into C_α packing preferences:

(6)

Here and impose a bias to the specific vertex orientation of regular H-bonds. defines the conditions when the ith residue is hydrogen bonded to the jth residue, i.e.,

(7)

Secondary structure assignments (when predicted) modify the formation of H-bonds: H-bonds between extended-assigned and helical-assigned residues and long-range H-bonds between helical-assigned residues are prohibited. Moreover, to enhance the H-bond in the better assigned secondary structure regions, we set the stiffness modulation factor λ′ to 1.5 or 1, respectively, depending on whether or not regular helix and sheet structures are predicted.

Local distant restraints The consensus local distance predictions for pairs of C_αs less than six residues along the sequence are collected from the templates and short fragments hit by our threading program PROSPECTOR (Skolnick and Kihara, 2001). These protein-specific predictions are incorporated in the force field as loose restraints on the local structure:

(8)

(9)

The accumulated normalized deviations to the predicted distant map enter into the force field as a penalty when they exceed the number of predictions, N_dp. This penalty term allows for the significant violation of a small fraction of unreasonable predictions.

Long-range interactions

Pairwise interactions The long-range pairwise interactions of C_α(_β)-C_α(_β) and that of SG-C_α(_β) are essentially the general excluded volume interactions, which as mentioned above are represented by a smaller hard-sphere potential plus a 1/r type of soft-core potential with a slightly larger range. The SG-SG interaction is written as

(10)

(11)

Here, γ_i,j denotes the relative orientations of the bisector vectors of the backbone vertices, i.e., parallel (), antiparallel (), and perpendicular (). and are the cutoff values for the hard-core excluded volume interactions and the soft-core square-well interactions, respectively. The pairwise potential is derived from a structural data base (Skolnick et al., 1997). aims to enhance the contact interactions of the short fragments that have favorable short-range correlations and therefore stabilize their local structures.

Burial interactions This potential represents the general propensity of amino acids to be buried or exposed to solvent and is only applicable to single-domain proteins. It includes contributions from both the C_αs and the SGs:

(12)

Here r_i is the radial distance of the ith C_α related to the protein center, r₀ is the average radius of gyration, which has an approximate relationship with the length of protein N, i.e., . is a statistical potential derived from PDB data base, where r_i was divided into five bins for each amino acid A_i. is a half/half combination of the two most commonly used hydrophilicity scales of Kyte-Doolittle (Kyte and Doolittle, 1982) and Hopp-Woods (Hopp and Woods, 1981). μ_i is the burial factor relative to hydrophobic core defined as

(13)

Electrostatic interactions We consider electrostatic interactions among four charged residues, i.e., Asp (−), Glu (−), Lys (+), and Arg (+), in a Debye-Huckel form:

(14)

Here, κ is the inverse Debye length that is sensitive to solvent conditions (Zhang et al., 2001). Through examination of the potential on structure decoys, we found that a value of 1/κ ~ 15 Å produces the best correlation between the RMSD and energy.

Environment profile The potential describing the contact environment of individual residues is written as

(15)

Contact order and contact number We also include biases to the expected contact order and contact number:

(16)

Here, is the contact order of the given structure, defined as average sequence separation of residues in contact (Baker, 2000). The expected contact order has an approximately linear dependence on protein length N, i.e., , where α is a protein secondary structure specific parameter that is derived from the PDB data base, which was divided into three categories of helix, sheet, and helix/sheet proteins. is the number of contacting residues, and is an approximate estimate of the contact number according to the PDB.

Contact restraints Consensus tertiary contact predictions are collected from templates hit by the threading program PROSPECTOR (Skolnick and Kihara, 2001), where sequence homologs have been excluded from the data base. These predictions are incorporated into the force field as

(17)

(18)

There are 19 parameters in Eq. 18, which dictate the relative weights of the different energy terms. We could not combine them naïvely, i.e., let all w_i = 1, because the energy terms are not independent and some interactions are multiply counted. For example, the short-range five-residue correlation energy E₁₅ partly includes the contributions of lower-order correlation energies E_1i (i < 5); the former is also incorporated in the calculations of pairwise interactions. The propensity to regular secondary structure is implemented in different energy terms such as hydrogen bonding, conformational stiffness, and pairwise interactions. Thus in the following, we will first generate a set of nonredundant decoys, and then determine the parameters by maximizing the correlation between the energy and the structural similarity of the decoys to native.

Generation of decoy structures To generate decoys, we selected 30 nonhomologous protein sequences from the PDB (Berman et al., 2000), which cover a variety of lengths (47 ~ 146) and topologies (see proteins marked with in Table 6). We make Monte Carlo runs based on both the SICHO (Kolinski et al., 1998) and CABS force fields using the parallel hyperbolic sampling algorithm (Zhang et al., 2002). To perform the CABS runs, we made a temporary initial estimate of the force field parameters. These simulations start from the native structure. For reasonable force fields, the low temperature replicas stay around the near-native state and the higher temperature replicas move away and generate structures further away from native. If the model force field is not good enough and even low temperature replicas go far away from native, we intermittently stop and restart the simulation from native structures to ensure that a sufficient number of decoys are near native.

TABLE 6 Summary of fold results on 125 benchmark proteins

The decoys are collected from the structure trajectories in all high- and low-temperature replicas. To avoid the overaccumulation of some structure clusters, we introduce a cutoff on the RMSD of structure pairs and ensure that the RMSD of any pair of decoy structures are larger than 3.5 Å. The decoys produced in this way retain their secondary structure and side-chain packing pattern in the low- and middle-temperature replicas. To neglect bad random coil structures present in the high temperature replicas, we remove structures whose radii of gyration are larger than 3N^1/3. The simulation continues until 60,000 decoys are generated for each protein.

Because the force fields are different in the SICHO and CABS models, these two simulations cover different regions of configurational phase space; this is helpful for the divergence of the decoy sets. As shown in Fig. 3 a, the rank of native structure in the decoys produced by the SICHO model is poor if the decoys are evaluated by the SICHO force field; however, if the same decoys are evaluated by the CABS force field, the rank of native structure is much better (Fig. 3 b). Similarly, if the decoys produced by the CABS model simulation are evaluated by CABS force field itself, the rank of native structure is poor (Fig. 3 d); however, if these same decoys are evaluated by the SICHO force field, the rank of native structure is much better (Fig. 3 c). This is a general feature seen in all the decoy sets on the 30 selected proteins; this means that, when the force field used for structure evaluation is different from the force field used for structure generation, it is possible that we can get better identification of native structure than if both force fields are the same. This is understandable because the Monte Carlo simulations always detect the so-called “important phase space” regions that are of low energy. Because of imperfections of the force field, this lowest energy basin usually does not correspond to the native state in most cases (see Fig. 3 e), so the rank of native structure in those decoys produced by the force field itself is poor. Because of the differences in the two force fields, the states in the lowest energy basin of the first force field can be of high energy in the second force field. But the idea is that native structure should be of relatively low energy in a reasonable force field. Therefore, the rank of native structure can be relatively better when ranked by the second force field (Fig. 3 e).

FIGURE 3 Energy versus RMSD of decoys to native structure of protein 1cis_. (a) Decoys generated by Monte Carlo simulations of the SICHO model, energies of decoys are evaluated by the SICHO force field. (b) The same decoys as in a but the energies are evaluated (more ...)

Parameter optimization The aim of our parameter optimization procedure is: i), to maximize the correlation between the energy function of the decoys and the RMSD to the native structure; and ii), to maximize the energy gap between the native state and the ensemble of unfolded states. For the 30 × 60,000 decoy structures, we try to find a set of parameters to minimize the following equation:

(19a)

(19b)

(19c)

(19d)

Here is the RMSD of jth decoy structure of kth training protein. We have a cutoff on the RMSD, i.e., = 4 if RMSD < 4 Å, = 10 if RMSD > 10 Å, because we consider any decoy with a RMSD < 4 Å as good and a RMSD > 10 Å as poor. N_p is the number of undetermined parameters (w_i values) of force fields, is the energy term conjugate to the parameter w_i. denotes the average over the decoys.

The first term G₁, of Eq. 19 b, aims to maximize the correlation coefficient between the RMSD and the total energy. The second term G₂ of Eq. 19 c acts to minimize the χ² between a linear regression (R_k = ηE + b_k) and the energy versus RMSD, where b_k is the individual intercept for the kth training protein, η is the slope of the fit line. Although b_k and η are irrelevant for the determination of the best force field, η decides the scale of the energy function that is related to the temperature range using MC simulations. We will determine η from the simulations. Although both G₁ and G₂ try to enhance the correlation of the energy function to the RMSD from native, the combination of these two terms speeds up the convergence of the optimization procedure and gives better results than when using either one of them alone. Finally, the aim of the third term G₃ is to maximize the relative gap between the native structures and the ensemble of the decoys of all 30 training proteins.

Because the weight parameters w_ri in Eqs. 8 and 17 depend on the results from threading, based on Eq. 19 a we at first optimize the 15 inherent parameters of w_i of the intrinsic force field with the threading-based restraint parameters w_ri = 0. In the second step, we have the 15 w_i values fixed at their optimized values and we optimize the remaining four threading parameters w_ri. To obtain the optimized values of 30 + N_P parameters in Eq. 19, we develop a minimization approach based on the CERN MINUIT package (James, 1998), which can handle and find the global minimum of a generation function of up to 100 variable parameters. To avoid some unphysical subminima and to speed up the optimization processes, we have put a loose physical restriction on each parameter.

In Fig. 4 a, we show an example of the energy versus RMSD correlation for 1fas_. If we simply add all the subenergy items with naïve weight factor w_i = 1, the global minimum of the force field is ~8.5 Å away from native structure and the correlation coefficient of total energy and RMSD is 0.44 (Fig. 4 b). With the optimization of the weight factors, the global minimum state is much closer to native and the energy versus RMSD correlation coefficient equals to 0.69 (Fig. 4 c).

FIGURE 4 The energy versus RMSD for the decoy structures of 1fas_ produced by the CABS model. (a) Correlations of 19 subenergy terms with the RMSD to native. (b) Combined energy with wi = 1. (c) Combined energy with optimized weight parameters.

In Table 2, we show the average correlation coefficients and z-scores of the different energy terms over 30 × 60,000 decoys. It is shown that the combined energy with optimized weight factors has higher correlation coefficients and nativelike recognition capability than the naïve combination of energy and each of the single energy terms alone.

TABLE 2 Summary of the CABS force-field weights

Because the energy barriers can be too high for the simulation to cross, it is well known that the canonical Metropolis protocol usually results in the simulations being trapped in local energy minima in rugged force fields (Newman and Barkema, 1999). In recent work (Zhang et al., 2002), we developed a new parallel hyperbolic sampling (PHS) algorithm to alleviate the problem of “ergodicity breaking.” The point of this algorithm is that the local high-energy barriers are flattened by a nonlinear transformation, i.e.:

(20)

In this work, we will use the PHS algorithm as the energy update protocol for the CABS model. The conformational update is first applied on the C_α chain. Then positions of C_β and SG units are determined accordingly. Five kinds of C_α-chain movements are used in our simulations.

• Movement 1: Basic 2-bond and 3-bond movements (Fig. 5 a), in which a 2- or 3-bond fragment is replaced by a fragment of the same length, but with a new conformation. Because of the limited number of conformations of 2- and 3-bond fragments on the lattice, all basic moves can be prefabricated, i.e., they are calculated only once and then randomly selected during the simulation. With the current lattice, we have 67,272 2-bond fragments and 14,507,376 3-bond fragments.

  FIGURE 5 Schematic diagrams of the movements employed in the Monte Carlo simulations. The Cα-traces before and after movements are denoted by the solid and dashed lines, respectively. (a) A basic prefabricated 3-bond update of the fragment [i, i + (more ...)
• Movement 2: 4-, 5-, and 6-bond movements (Fig. 5 b), which consist of consecutive 2- and 3-bond moves.
• Movement 3: 6- to 12-bond translation (Fig. 5 c), in which a randomly chosen fragment of 6–12 bonds is translated over a small distance.
• Movement 4: Multibond sequence shift (Fig. 5 d), which is performed through a permutation of a randomly chosen 2-bond piece and another randomly chosen 3-bond piece. Because the conformation of the fragment between the permutation points is not modified, the net result of this permutation is a sequence shift along the modeling chain as marked by the arrows in Fig. 5 d. Although the acceptance probability of this movement can be quite low, it can substantially increase the probability of extrusion and resorption of tangled structures and help the simulation get out of some local energy traps.
• Movement 5: Extremity movements (Fig. 5 e), which reconstruct the conformation of the N-or C-terminus through a random walk from a chosen point to the extremity.

In each of the above randomly chosen movements, a geometric restriction on the virtual C_α-C_α bond angles to lie in the range of [65°, 165°] is put on all new conformations. The smaller moves with higher acceptance rates are performed with greater frequency, which lead to a better simulation of the process of the fine repacking of side chains after a larger change of the main chain local geometry.

Because only the energy difference between two conformations is involved in Eq. 20, in each step of updates we only need to calculate the energies of the fragments whose conformation changed to save CPU time. Before any energy computation, the test for excluded volume violation of the C_α and C_βs are always performed, and trial conformations that would lead to steric collisions of chain units are rejected.

Table 3 shows the lowest energies identified by different algorithms using different move sets for the same CPU time and demonstrates how the two aspects of the energy update protocol and movements influence the efficiency of Monte Carlo simulations. The 20 test cases cover protein lengths from 36 to 174 residues. For the same algorithm, the simulations with a more comprehensive move set always do better than these including only basic 2- and 3-bond movements, because the larger moves can cross over local energy barriers more efficiently. When using simple movements, the PHS algorithm does significantly better than the RS algorithm, because the local energy obstacles, which are difficult to surmount by simple movements, are flattened in the PHS simulation. When using combined move sets, there is no obvious difference in the performance of the PHS and RS simulations for small proteins (say, <100 residues). However, for larger proteins, the PHS simulations almost always identify lower energy structures than the RS simulations do. This may mean that the roughness of energy landscape is correlated with protein length. For small proteins, the local energy barriers are not too high and can be surmounted when using a larger set of movements. For large proteins, however, the local energy barriers are still difficult to surmount with the combined movements. So the flattening of energy landscape improves the sampling.

TABLE 3 Lowest energies found in test simulations of different algorithms and move sets

The prediction accuracy of secondary structures has been considerably improved with the utilization of the multiple sequence alignments (Benner and Gerloff, 1991). It was found that the secondary structure information can be extracted from the sequence evolutionary information (Branden and Tooze, 1999). In Table 4, we show the results of secondary structure predictions on 125 test proteins based on the three most-often-used sequence-based predictors: PHD (Rost and Sander, 1994), SAM-T99 (Karplus et al., 1998), and PSIPRED (Jones, 1999). The average prediction accuracy of the single predictor for the 125 proteins fluctuates from 73.4% to 81.0%, depending on the cutoff of the confidence level for α-helix and β-strand assignments. The accuracy of PSIPRED is slightly better than SAM-99 in our test set, and the accuracy of both prediction methods is better than PHD. The highest prediction accuracy comes from the combination of PSIPRED and SAM-T99 results, where two ways of combination of “overlap” and “consensus” are defined as in Table 5. We have done test runs using six sets of highest secondary structure prediction accuracy (see italic bold numbers in Table 4) in our fold simulations, after the optimization of the force field. The tertiary structure prediction results depend on both the accuracy and the coverage of the secondary structure predictions. The overlap set with a cutoff equaling to 5 and 0.49 for PSIPRED and SAM-T99, respectively (see italic bold numbers in Table 4) works the best. This is used in all subsequent simulations.

TABLE 4 Accuracy and coverage of secondary structure prediction by different predictors

TABLE 5 Combinations of two secondary structure predictors

In this section, we will report the results of applying our methodology to a test set of 125 proteins. We first check the folding ability and convergence of our basic force field (without restraints) on 100 small proteins. Then, we use the methodology on the whole set of proteins under the guide of threading-based restraints. Finally, we describe the protocol of selecting structures from the generated trajectories.

In column four of Table 6, we list the folding results of the CABS model without using predicted protein-specific local and tertiary restraints provided by our threading program. If we define a “successful” fold as one in which at least one of the top five clusters has the RMSD to native below 6.5 Å, we can successfully fold 41 cases using the basic force field. There is an obvious bias of fold success to the protein secondary structure class: 21 foldable cases are α-helical proteins, nine are β-sheet proteins, and 11 are mixed α/β proteins. All the successful folds occur on the 100 small proteins of length N < 120 amino acids. The dependence of RMSD on protein size is shown in Fig. 6 a with both testing (denoted as solid circles) and training (denoted as open circles) protein sets. It shows that the folding results have no obvious bias to the training protein set (13 foldable cases belong to the 30 training proteins, compared with 41 foldable cases to 100 proteins in total). This may mean that the training set of 30 proteins is sufficiently large and representative for a general optimization of the force field.

FIGURE 6 (a) RMSD of the best cluster in the top five clusters versus protein length N in the CABS simulations without using protein-specific restraints. The solid circles denote the training proteins that are used in the optimization of force field. The open (more ...)

To fold proteins longer than 120 residues and to improve the yield of small proteins, we exploit the predicted local and tertiary restraints in our force field (see Eqs. 8 and 17). These restraints are collected from consensus substructures found by our threading program PROSPECTOR (Skolnick and Kihara, 2001), where homologous sequences to the query protein are entirely excluded from the data base. Although a portion of the predicted restraints may be incorrect, it is indeed helpful to guide the simulations to near-native states and significantly improve the folding results in the majority of cases. In column five of Table 6, we list the results of the simulations with restraints. There are 83 cases with a RMSD of the best cluster centroid below 6.5 Å to native, all within the top five clusters. Fifty-one successful cases are from the 65-protein set and 32 are from the harder 60-protein set (The trajectories and cluster centroids of all the 125 proteins are available on our website: http://www.bioinformatics.buffalo.edu/abinitio/125).

The improvement using tertiary restraints occurs on both small and large proteins (see Fig. 6 b). For the 100 small proteins of lengths less than 120 residues, the number of foldable cases with restraints increases to 70 (compared to 41 without restraints). Without restraints especially, the program can never fold proteins of lengths longer than 120 residues. Under the guide of restraints, we can fold 13 of the 25 large proteins; none can be folded without predicted side chain contacts. Moreover, within all 83 cases, 33 cases belong to α-helical proteins, 27 cases to β-sheet proteins, and 23 to mixed α/β proteins, which show a considerably reduced folding bias to the secondary structure class, compared to the pure ab initio results.

The effect of restraints on the degree of folding success depends on its accuracy. In Fig. 7, we show the dependence of the fold improvement on the accuracy of the predicted contacts and local distant restraints. There is a strong correlation between the RMSD improvement and the accuracy of long-range contact restraints. This correlation is much less obvious for the local distance restraints, which seems to indicate that the local short-range restraints are less important. This may be due to the fact that the information of short-range correlations has been included due to the relatively high accurate secondary structure prediction and the short-range distance restraints do not provide much additional information. However, our simulations show that appropriate short-range restraints indeed considerably speed up the formation of local structures.

FIGURE 7 RMSD improvement on including the threading-based tertiary and secondary restraints versus the accuracy of the restraints. , where RMSDwo and RMSDw are the RMSD of the best clusters to native structures in the simulations without and with using the threading-based (more ...)

As expected, when the accuracy of contact restraints is too low, a successful fold from a “pure” ab initio simulation can be spoiled by inclusion of poorly predicted restraints. According to Fig. 7, when accuracy of contact restraints is higher than 22% or the ratio of the number of correct restraints to protein length is larger than 0.2, the restraints almost always have a positive effect on folding. To alleviate the negative influence of the bad restraints, we combine the clusters from both simulations as follows: The best cluster is the lowest energy cluster in most of the successful pure ab initio simulations because only a good ab initio force field can fold a protein without restraints. Thus, we take the lowest energy cluster from the pure ab initio simulations and combine it with the four lowest energy clusters from restraint-based simulations. As shown in column six of Table 6, this combination converts all the significant spoiled cases by the inclusion of poorly predicted restraints into successful folds. Moreover, we retain all the successful folding cases in the restraint-based simulation set.

As a comparison, we also made Monte Carlo runs of the SICHO model on the harder subset of proteins with similar CPU times. The results are shown in the histogram in Fig. 8. It should be noted that these 60 proteins represent diverse structure categories, and no protein from this set was used in the training of either the CABS or the SICHO force field. The folding rate is 1/3 for the SICHO model and 1/2 for the CABS model. However, in fairness, the SICHO model has not yet been subjected to the same optimization procedure as done in the CABS model.

FIGURE 8 Comparison of the folding results by the SICHO and CABS models on the 60-nonhomologous-protein set. The shown data are the number of proteins that have their best cluster below a given RMSD threshold versus the RMSD threshold.

In what follows, we first address the issues of how to select the best structure from an individual simulation trajectory. We introduce several quantities that are highly correlated with the likelihood of successful fold selection. We perform five sets of simulations under different restraints, and present an automatic procedure to select the best structures from the multiple simulations by combining appropriate fold selection criteria.

Selecting the best fold from an individual simulation In previous approaches to ab initio structure predictions, the authors usually cluster the generated structures (Shortle et al., 1998; Betancourt and Skolnick, 2001) and choose the cluster with the lowest energy (Kolinski et al., 2001; Kihara et al., 2001). Although clustering has considerable success in selecting the correct folds (Simons et al., 2001; Kihara et al., 2001), this approach can have an inherent contradiction. Although the aim of clustering is to identify the low free energy structures, the selection of the lowest energy structure neglects the configurational entropy, because the structure of lowest energy is not necessary that of lowest free energy. Thus, we consider a combination of the energy and free energy

(21)

Another relevant indicator of the quality of the predicted structures is the normalized structure density of cluster defined as

(22)

In Fig. 9, we show the RMSD to native of all the cluster centroids versus their normalized structure density. There is a strong correlation between the fold quality and the structure density. If we define the best cluster as a cluster of lowest RMSD, most of the best clusters (denoted by solid circles) have higher structure density as compared to the high RMSD clusters. As shown in column 10 of Table 6, by selecting the highest-density cluster, we choose the best fold in 76 of 125 cases and 67 of the chosen clusters have a RMSD below 6.5 Å.

FIGURE 9 RMSD to native of all cluster centroids for 125 proteins versus the normalized structure density. The solid circles denote the best clusters of lowest RMSD to native in each of the 125 proteins.

Indicator of likelihood of success of the folding simulation Now we turn to the issue of how to judge the likelihood of success of a blind simulation. As mentioned above, if the force field is a combination of consistent and reinforcing energy subterms, the resultant landscape tends to have a funnellike shape with deep energy basins, which results in convergent structure clusters in the fold simulation. Recent experimental studies of denatured state showed that this funnellike landscape is a basic and necessary characteristic of real proteins to keep the native structure as a stable and unique state (Shea and Brooks, III, 2001). This funneling characteristic can be quantitatively evaluated by the maximum cluster density D_max, or the maximum multiplicity rate of clusters , where M_max is the multiplicity of the largest cluster. It can also be a represented by the normalized Y-gap between the energy basin of lowest Y and other basins, i.e.:

(23)

In Fig. 10, we show the dependence of the RMSD of the best cluster on D_max, R_max, and L-score, respectively, demonstrating that these parameters can be considered as indicators of the likelihood of success of the simulations.

FIGURE 10 RMSD of the best cluster to native versus different funneling parameters of the energy landscape. (a) The maximum structure density Dmax. (b) The maximum multiplicity Rmax. (c) L-score of energy landscape (defined in Eq. 23).

In Fig. 11, we show the successful folding rate and average RMSD of best clusters versus the threshold values of maximum cluster density. With higher density cutoff, we have higher rate of successful folds and lower average RMSD. For example, for the simulations with D_max > 0.18, 95% of cases (63 of 66 cases) are successfully folded, and the average RMSD is 3.92 Å. This is dramatically better than the overall fold rate 66% (83 of 125) and the overall average RMSD of 5.90 Å. Furthermore, if we select the highest D cluster in these 66 cases of D_max > 0.18, 82% of them (54 cases) have RMSD below 6.5Å.

FIGURE 11 (a) Rate of successful fold (best RMSD < 6.5 Å) versus the cutoff of maximum density (Dmax > Dcut). (b) Average RMSD versus the cutoff of maximum density.

Automatic procedure of selecting top five clusters from multiple simulations To demonstrate the usage of the combination of above-defined parameters, we make five sets of simulations on the 60 hard proteins, each set taking different restraints that were obtained by using different threading procedures and cutoff parameters. On average, there are ~10 clusters for each protein in each individual run. To select the five best clusters for each protein from these roughly 50 clusters, we at first sort the clusters in each simulation according to Y and D, and the different simulation sets according to D_max. Then, we choose the five clusters according to following automatic procedures:

• Select the five clusters of highest D from five sets of simulations.
• If any pair of clusters is of the same fold (<2 Å), displace the cluster selected from lower D_max simulation with the cluster of lowest Y from the simulation of higher D_max.
• Repeat step 2 until five different clusters are chosen.

In column three of Table 7 we show the selection result according to the automatic procedure. Compared with the absolutely best clusters in column four, this procedure allows us to select almost all the best folds in the top five clusters (all 37 successful cases with a RMSD < 6.5Å). Column two shows the selection of five clusters if we choose them just according to cluster energy E in different sets of simulations, which is much worse than that by above combined selection procedure.

TABLE 7 Selection of top five clusters from multiple simulation runs

In this work, we have developed a new ab initio modeling approach to the tertiary protein structure prediction, based on a simplified lattice representation of the C_α, C_β, and center of side group of protein chains. This new lattice description has a high geometric fidelity. The basic energy function consists of general short-range correlations biased to regular and predicted secondary structures, amino acid-dependent short- and long-range interactions derived from the PDB data base, hydrogen bonds, electrostatic interactions, one-body burial interactions, and a general bias to the expected contact order and contact number that depends on protein size and secondary structure. These energy terms from different sources are combined and optimized by a set of 30 × 60,000 nonredundant structure decoys, by maximizing both the correlation of RMSD of decoys to native and their energies, and the relative energy gap between native and decoy ensemble. This combined force field provides a basic working platform for further assembly and optimization of tertiary structures when threading information (i.e., predicted side-chain contact restraints) are available. It has also shown to be able to successfully assemble structures from sparse NMR experimental NOE data (Li et al., 2002). Here, we used the intrinsic platform (without restraints) on the folding experiment of 100 small proteins (<120 amino acids). 41% of them can be successfully folded with the best RMSD of the top five clusters below 6.5 Å. Twenty-one foldable cases are α-helical proteins, nine are β-sheet proteins, and 11 are mixed α/β-proteins. There is no obvious bias to the training set (13/30 foldable cases for training proteins compared to 41/100 foldable cases in total), which demonstrates that the training set of decoys is large enough for a universal derivation of the force field.

The long-range contact prediction and short-range distance prediction are collected from templates found by our threading program PROSPECTOR (Skolnick and Kihara, 2001). These data are incorporated into our CABS force field as loose side-chain pairwise and local distance restraints. It should be mentioned that, even when no template is hit with significant z-scores in the threading program, some useful information could still be extracted from the consensus substructures with weak z-score hits. These threading-based restraints in most cases can significantly improve the folding results, even if the prediction accuracy is low. More specifically, when the accuracy of contact prediction is higher than 22% or the ratio of correctly predicted contact number to protein length is larger than 20%, the effect of restraints on the folding is almost always positive. There is no obvious sensitivity on the accuracy of local short-range distance predictions. This may be because short-range interactions are already dictated by the high accurate secondary structure predictions (the combination of the PSIPRED (Jones, 1999) and SAM-T99 (Karplus et al., 1998) secondary structure predictors) that have been incorporated in our force field, and therefore the short-range restraints do not provide much additional information.

The improvement by including the tertiary restraints occurs for both small- and large-size proteins. For 100 small proteins <120 residues, the program can fold 70 cases with restraints compared to 41 in pure ab initio folding. The intrinsic force field, especially, can never fold proteins >120 residues in length. Under the guide of restraints, however, the program can fold 13 cases of the 25 larger proteins with lengths ranging from 120 to 174 residues. Overall, in the restraint-based simulations, 33 foldable cases belong to α-helical proteins, 27 belong to β-sheet cases, and 23 belong to α/β-proteins, which shows a less obvious bias toward the protein topology category than the pure ab initio simulations.

We found that the structure density of cluster D and combination of energy and free energy Y are more discriminative than the often-used energy E or cluster size M in the selecting of best-folded structures. In the folding of 125 proteins, if we select one cluster according to the lowest E, or biggest M, or lowest Y, or highest D, the numbers of cases that select the lowest RMSD are 61, 67, 73, and 76, respectively. The numbers of cases that have a RMSD below 6.5 Å in the first cluster in these selections are 58, 60, 65, and 67, respectively.

The coherence of the energy terms in the force field and the funnellike characteristics of the energy landscape can be quantitatively evaluated by the maximum cluster density, maximum multiplicity rate, or L-score. There are strong correlations between the best RMSD and these funneling parameters, which demonstrate that the parameters can be used as indicators of the likelihood of success of fold simulations. For the simulations of 125 proteins, if we take a cutoff of maximum density D_max > 0.18, 95% of cases (63 of 66) are successfully folded, which is much larger than the overall folding rate of 66% (83 of 125).

The combination of the discriminative parameters and the indicator parameters of likelihood of success folds can be used for selection of the best structures from multiple simulations that are run using different force fields (e.g., based on different tertiary restraints). In an evaluation of five sets of test simulation runs, by sorting the simulations according to the indicator parameters and sorting the clusters according to the discriminative parameters, we can select almost all the absolute best structures in the top five chosen clusters. This could not be achieved by the selection based on traditional average energy or cluster size. Because our procedures are fully automatic from the trajectory generation to the identification of final structures, these approaches can be applied to large-scale structure predictions. A comprehensive prediction survey of PDB structure data base and the subsequent genome-scale structure predictions based on these approaches are in progress.

This research was supported in part by the Division of General Science of the National Institutes of Health (grant numbers GM-37408 and GM-48835). We are also grateful to Dr. Adrian K. Arakaki for help in the preparation of the figures.