A Novel Side-Chain Orientation Dependent Potential Derived from Random-Walk Reference State for Protein Fold Selection and Structure Prediction

Background An accurate potential function is essential to attack protein folding and structure prediction problems. The key to developing efficient knowledge-based potential functions is to design reference states that can appropriately counteract generic interactions. The reference states of many knowledge-based distance-dependent atomic potential functions were derived from non-interacting particles such as ideal gas, however, which ignored the inherent sequence connectivity and entropic elasticity of proteins. Methodology We developed a new pair-wise distance-dependent, atomic statistical potential function (RW), using an ideal random-walk chain as reference state, which was optimized on CASP models and then benchmarked on nine structural decoy sets. Second, we incorporated a new side-chain orientation-dependent energy term into RW (RWplus) and found that the side-chain packing orientation specificity can further improve the decoy recognition ability of the statistical potential. Significance RW and RWplus demonstrate a significantly better ability than the best performing pair-wise distance-dependent atomic potential functions in both native and near-native model selections. It has higher energy-RMSD and energy-TM-score correlations compared with other potentials of the same type in real-life structure assembly decoys. When benchmarked with a comprehensive list of publicly available potentials, RW and RWplus shows comparable performance to the state-of-the-art scoring functions, including those combining terms from multiple resources. These data demonstrate the usefulness of random-walk chain as reference states which correctly account for sequence connectivity and entropic elasticity of proteins. It shows potential usefulness in structure recognition and protein folding simulations. The RW and RWplus potentials, as well as the newly generated I-TASSER decoys, are freely available in http://zhanglab.ccmb.med.umich.edu/RW.


Introduction
The basic hypothesis of protein folding theory is that protein structure generally has the lowest Gibbs free energy in the native state [1]. Therefore, an accurate energy function is the key to solve the protein folding and protein structure prediction problems. The commonly used potential function can be divided into two categories [2]. The first is physics based potential [e.g. AMBER [3], CHARMM [4] and GROMOS [5] etc], which can in principle be derived from the laws of physics. Although atomiclevel structure refinement can be achieved with the molecular dynamics (MD) simulations in some isolated instances, no systematic structure improvement has been observed [6,7,8]. The second is knowledge-based potential [e.g. RAPDF [9], KBP [10], DFIRE [11], DOPE [12], OPUS-PSP [13,14], free-rotating chain-based potential [15], or the more composite TASSER/I-TASSER [16,17,18] and ROSETTA [19] potentials], which is derived from the statistical regularities [20] of the solved protein structures in the PDB library [21].
The Knowledge-based potentials include contact potentials [22,23], orientation-dependent potentials [13,14,24], and distance-dependent potentials [9,10,11,12,25,26,27,28]. According to the reference state calculations, the distance-dependent potentials can be further divided into two classes: that using statistical reference states [RAPDF [9] and KBP [10]] and that using analytical reference states [DFIRE [11] and DOPE [12]]. It has been argued that the analytical reference state potential has better performance [11,12]. For example, the DFIRE potential used a reference state derived from a set of uniformly distributed noninteracting points in finite spheres [11]. DOPE [12] later introduced an improved reference state which used noninteracting atoms in a homogeneous sphere with the radius dependent on a sample native structure [12]. Both DOPE and DFIRE were derived from a non-interacting ideal gas reference state and the major difference is that DOPE also takes into account the size effect of proteins.
The Knowledge-based potentials were successfully applied to many areas, including fold recognition [23,29,30,31], ab initio protein structure prediction [16,19,32,33,34], protein structure refinement [13,14,35], structural model assessment [9,10,11,12], protein-protein docking [22] and protein stability prediction [11,22]. Despite the success of the potentials, more accurate accounting of atomic interactions will undoubtedly increase the power of the potentials in each of the application areas. In general, a protein is essentially a continuous sequential chain of the amino acid residues. The reference state, which accounts for the expected number of atom pairs at certain distance when interactions vanish, should correctly reflect and counteract the inherent chain connectivity effect. This feature, however, cannot be captured by the current ideal gas based reference state. Recently, Cheng et al. showed that a more physical reference model, such as free-rotating chain-based potential, could improve the performance of statistical potentials [15]. Aloy and Oliva introduced a method to split the knowledge-based potentials in biologically meaningful terms which allows a better combination of most relevant scoring functions [36]. Rykunov and Fiser performed a systematic comparison of publicly available scoring functions on CASP decoys which shows a critical role of reference state definitions. Based on the observation, the authors developed a residue based potential that employs a shuffled reference state with considering side-chain orientations and demonstrates advantages in structure decoy recognition [37].
In this work, we proposed a new distance-dependent atomic potential using a random-walk ideal chain as the reference state. This reference state was derived from a linear freely-jointed chain model, which can be considered as the segments of an ideal polymer chain performing a random walk (or ''random flight'') in three dimension space. We term the new potential ''RW potential''. The orientation-dependent all-atom potential, such as OPUS-PSP (it used a set of 19 rigid-body blocks extracted from the chemical structures of all 20 amino acid residues), can capture the feature of side-chain packing [13]. In this paper, a new orientation-dependent potential term was also added to RW. 20 vector pairs were defined to describe the side-chain orientation of 20 amino acids. The orientation term was then generated from the orientation specific packing statistics of those vector pairs in a nonredundant high-resolution structural database. The RW potential and the hybrid potential (RWplus) were optimized on CASP models and tested on eight commonly used decoy sets, as well as a new decoy set from real-life I-TASSER structure assembly followed by MD refinements. Detail comparisons with the stateof-the-art potentials demonstrated the advantage of the new reference state of chain connectivity and the side-chain orientation specificity.

Results
We tested our potential in three ways: (1) the ability to select native structures from structural decoys; (2) the ability to select the best models from structural decoys when the native structures are excluded; (3) the correlation between the potential and the similarity (TM-score and RMSD) of the structural decoy to the native.
As a control, we compared the results of RW and RWplus mainly with two frequently used atomic potentials, DFIRE [11] and DOPE [12]. DFIRE was developed by Zhou and Zhou [11] and we calculated the DFIRE score by the DFIRE program, which is provided by the authors (http://sparks.informatics. iupui.edu/download/ddfire_bin.tgz) [38]. DOPE was developed by Shen and Sali [12] and we calculated DOPE scores from the MODELLER-9v7 package (http://salilab.org/modeller). In the end of the section, we presented a comparison of RW and RWplus with all potentials listed in the Rykunov and Fiser benchmark set [37].

Testing on native structure selection
The ability of native structure selection of DFIRE, DOPE, RW and RWplus is tested using eight independent decoy sets (see Methods), where the experimental structures are mixed with other decoys generated by computers. The purpose is to rank the native structure as the lowest energy conformation using automatic scoring. Meanwhile, the significance of the energy of the native structures (E native ) is evaluated by the normalized energy gap between E native and the average energy of all decoys (E aveage ), i.e. Z-scoreẼ native {E average À Á s, where s is the energy deviation of all decoys. The results of RW, RWplus, DFIRE and DOPE on the native structure selections are listed in Table 1. While there are some fluctuations for the selection ability of different potentials among different decoy sets, RWplus potential correctly identified 123 native structures for a total of 168 targets with a success rate of 73%. The RW potential correctly identified 120 native structures for a total of 168 targets with a success rate of 71%. DFIRE and DOPE were successful for 115 and 98 targets, resulting in a total success rate of 68% and 58% respectively. The improvement of RW and RWplus was also reflected by the Z-score of the native structures. DOPE. Among the eight decoy sets, RWplus and RW have the lowest Z-scores for six decoy sets (fisa, fisa_casp3, lmds, lattice_ssfit ROSETTA and I-TASSER). For the remaining two decoy sets (4state_reducred and Moulder), the Z-scores of all potentials are worse than 24.0 and the selections of different potentials are somewhat random. This is mainly due to the quality of the decoy sets, for example, having poorly packed native structures. RW with additional orientation energy term has a consistent better performance than RW. The average Z-scores of RWplus are lower than RW for seven out of eight decoy sets and the successful selection rate of RWplus is 2% higher than RW. This improvement is due to the contribution of the orientation dependent energy term, which cannot be counted by the pairwise distance dependent potential. With orientation energy term, the most-probable side-chain packing patterns in high-resolution experimental structures, such as p-p and cation-p interactions, can be correctly recognized and be assigned lower orientation energy than the less favorite patterns. Thus the RWplus energies of the native structures are lower than RW and average Z-scores values of RWplus are much better.

Selection of best models from I-TASSER and ROSETTA decoys
The ability to identify native structure from structural decoys is only a minimum request to measure the potentials. Although the selection of the native structures has been a common goal of many protein potential developments [9,10,11,12,39], the usefulness of the criterion is limited. First, there are no native structures which are generated from computer simulations, and all computer models of structure predictions have some level of errors. Second, since the experimental structures are usually perfect conformations in many aspect of features (i.e. H-bonding, atomic clashes, secondary structure regularities, rotamer optimizations, electrostatics interactions etc), it is a relatively easy task to pick out the native structure from a set of computer decoys. On some occasions, a simple counting of some special features (e.g. the atomic clashes) may be enough to distinguish the native structures from the roughly generated computer decoys. So, in what follows, we focus on the more challenging and realistic cases of identifying the best decoys from real-time simulations by I-TASSER [34] and ROSETTA [40,41], or examining the correlation of the energy with the quality of decoys (i.e. RMSD and TM-score to the native). In this respect, we do not consider the decoy sets generated by manual variation of the native structures because the quality of the decoys usually has a strong correlation with the radius of gyrations. We also exclude the decoy sets from homologous modeling because the decoys are usually biased to specific templates and the distance to the initial template may be an efficient metric for decoy recognition [18].
We used RMSD and TM-score as the two criteria for assessing the quality of every structural decoy. RMSD is defined as the root mean squared derivation of all Ca pairs of the decoy to the native structure. Because RMSD weights all distances equally, it is insensitive to the global topology for large RMSD of decoys (e.g. a mis-oriented decoy may have a big RMSD although the global topology in the core region is correct). TM-score [42] weights the large distance at a small weight which makes the magnitude of TM-score more sensitive to the topology rather than the outlier of the structures [43,44]. TM-score ranges in (0, 1] where higher values indicate better quality. Table 2 summarizes the result of best model selection by DFIRE, DOPE, RW and RWplus for the TASSER decoys. If we consider the first model as ranked by the lowest energy, the average RMSD of the first models by RW is 5.20 Å which is 0.4 Å and 0.1 Å lower than that by DFIRE and DOPE, respectively. RWplus has the lowest average RMSD 5.19, which is slightly better than RW. The average TM-score of the first model selected by RW is 0.569 which is also higher than that obtained by DFIRE (0.558) and DOPE (0.560) and RWplus has the best average TMscore (0.575). Apparently, none of the methods could select the absolute best structure as the highest rank model in the decoy sets, which has an average RMSD/TM-score = 3.3 Å /0.675. We also consider the quality of the best decoys which are in the top-five and top-ten lowest-energy decoys in Table 2. The selected models by RW are consistently closer to the native structure than those by DFIRE and DOPE.
Despite of the advantage of RW compared with other methods, we found that it could not select models better than those selected by the structure clustering program SPICKER [45], which was designed to identify the most frequently occurred structural state in the simulation. When we cluster the 500 decoys of I-TASSER, where the redundant decoys have been removed, the average RMSD and TM-score for the first model (CLOSC) are 4.99 Å and 0.572, respectively. If we run SPICKER in the original I-TASSER trajectories (i.e. the 12,500-32,000 conformations which include structural redundancy), the RMSD and TM-score for the first model are 4.84 Å and 0.589, respectively. Here, CLOSC in SPICKER is the structure decoy which is the closest to the cluster centroid (COMBO) where the COMBO structure is calculated by averaging all structural decoys in the cluster. Because the cluster identified by SPICKER has the highest multiplicity and partition function n:Z~Ð e {bE(s) v(s)ds where s is the conformation phase space and v is state density, it is actually selecting the state of the lowest Helmholtz free-energy, i.e. F H~{ k B T ln Z. These results show the advantage of selecting models from the lowest of inherent free-energies.
In Table 2, we also present the result of near-native structure selections by DFIRE, DOPE, RW and RWplus for the Table 2. Average RMSD (Å ) and TM-score (in parenthesis) of models selected from I-TASSER and ROSETTA decoy sets. ROSETTA decoys. RW potential consistently selected models closer to the native structure than those by DFIRE and DOPE in the top-five and top-ten lowest-energy decoys, while DFIRE selected the best first models with an average TM-score 0.469, which is slightly better than DOPE (0.466) and RW (0.460).

Correlation between potential score and modeling errors
Except for the ability of selecting good models from structure decoys, another important criterion of assessing the potential development is to examine the correlation of the potential with the similarity of decoys to the native structure [16]. This is to some extent more important to protein folding because it can determine how structure assembly simulations are guided to the near-native states. Certainly, a golf-hole-like potential may be perfect in selecting good models but it is useless in protein folding because it lacks a middle-range funnel in such an energy landscape.
In Table 3, we present the Pearson correlation coefficients between Ca RMSD (and TM-score) and the potential energies as given by DFIRE, DOPE RW and RWplus for the I-TASSER decoys. Overall, RWplus has the best correlation coefficients. RW and DFIRE have comparable correlation coefficients although the average correlation coefficient of RW is slightly higher than that of DFIRE. The correlation coefficients of all three potentials are much higher than DOPE. More specifically, the RWplus potential yields an average energy-RMSD correlation coefficient of 0.53, compared with that of RW (0.  The average energy-RMSD and energy-TM-score correlation coefficients for the ROSETTA decoys are also listed in Table 3. Again, RW and DFIRE have comparable correlation coefficients with the correlation of RW being slightly higher, while both of these are obviously higher than DOPE. RWplus has correlation coefficient between RW and DOPE for the ROSETTA decoys.

Comparison with other potentials in the Rykunov and Fiser benchmark set
A comprehensive benchmarking survey of quality assessment scoring functions relative to a list of other publicly available potentials is shown in Table 4. The data of the potentials were adopted from Rykunov and Fiser [37] who compared the GDT_TS scores of the models recognized by each of the potentials. The model decoys for the 143 protein targets were generated during CASP5-CASP8 experiments. Data in Table 4 are sorted by the average rank of the lowest energy decoy structure according to the GDT_TS score for the decoy set excluding the native structure. To obtain the correct GDT_TS scores and RW and RWplus scores, the models in Rykunov and Fiser decoy set were first cleaned up by removing the remarks and hetero atoms and the residue numbers in the models were then reordered according to the native structures. The data of RWplus and RW were calculated from the cleaned Rykunov and Fiser's decoy set, which can be downloaded from http://zhanglab.ccmb.med.umich.edu/RW/casp_good.tar.gz. The RWplus and RW ranked second and third place respectively and have comparative performance to the best potential QMEAN6 [46] for average rank with and without native structures, which is a composite potential combining six structural descriptors including distance, solvation, torsion, secondary structure predictions [46].
The performance of RWplus and RW varies depending on the presence or absence of the native structures. RWplus outperforms RW for average rank without native structures, but has worse performance for the average rank with native structures. RWplus can correctly select 57 best decoys for models without native structures which is 6 more than RW, whereas RW can correctly select 110 best decoys which is 4 more than RWplus. RW has significant better performance than other pair-wise distance dependent potentials of the same type of statistics, such as DFIRE [11] and DOPE [12], which indicate that the RW reference state, which mimics the entropic elasticity and chain connectivity, are efficient to counteract generic interactions.

Comparison of different pair-wise distance dependent statistical potentials
Most of the atomic statistical potentials in the literature used the same equation with the major difference in the derivation of the reference state. To examine the detailed differences of the overall potentials, we compared in Figure 2 the distance dependence of RW and DFIRE potentials for four representative pairs of atom types in main chain-main chain, main chain-side chain, hydrophobic side chain-hydrophobic side chain, and polar side chain-hydrophobic side chain groups.
For all four pairs, RW potential has a steeper repulsion at short distance than DFIRE, and thus can assign a lower energy to the atom pairs with a favorite distance and give favorable structure lower energy. For example, the Ile C-beta atom-Leu C-beta atom pair has a deeper valley at 6Å and a higher peak at 9Å , which increases the energy gap between good pairs and bad pairs and therefore also increases the sensitivity of the potential to the structural variations. These subtle changes are mainly due to the difference in calculating the reference state where DFIRE considers the reference state as idea gas and RW treats it as a freely-joint chain with chain connectivity. The overall distance dependences of the potentials are qualitatively similar, because

Comparison of reference states of different pair-wise distance dependent potentials
To examine directly the reference states of DFIRE, DOPE and RW potentials, we calculated the ratio of reference state at a distance R to that at a distance cutoff R cut ( = 15Å ) for a protein of 100 amino acids. For DFIRE, the expected number of atom pairs (a, b) in the distance shell R to RzdR [11] is where V is the volume of the ideal gas system and c = 1.61. N a and N b are the number of atoms of type a and b, respectively. For DOPE, the potential is derived from the distance probability density function [12]: where, p a,b,R ð Þ and p ref a,b,R ð Þ are the observed and reference distance probability density function of atom pair (a, b), respectively.
N N obs (a,b)~S P R N obs (a,b,R)T is the average number of observed atom pairs (a, b). Since p ref a,b,R ð Þequals to the normalization function n R,a ð Þ [12] and n R,a ð Þ~6 where a~ffi ffiffiffiffiffiffi ffi 3=5 p R g is the size of the sample protein structure and R g is the radius of gyration, the expected number of atom pairs (a, b) in the same distance shell can be written as where R c is some upper bound of the statistical potential. From Eq. (18), we can obtain the expected number of atom pairs (a, b) in the same distance shell for RW In Figure 3, we present the ratio of reference states at distance R to that at R cut versus R for FIRE, DOPE and RW. It is shown that the RW potential has a lower ratio than DFIRE and DOPE at short distance, whereas at long distance the ratio of RW is similar to that of DFIRE but lower than that of DOPE. This difference makes the RW potential a steeper potential at short distance as showed in Figure 3 and therefore help increase the sensitivity of the potential to the short range interactions.

Estimation of Kuhn length b and distance cutoff R 0
There are two tuning parameters, the Kuhn length b and distance cutoff R 0 , in the RW potential derivation.
The Kuhn length b was introduced to match the scale of the FJC with that of real protein chains. We found that the optimized value of l( = b2), at which RW achieves the best performance, equals to 460; this corresponds to a Kuhn length b = 21.4 Å . The value coincides with the data of the single molecule stretching experiments with atomic force microscope [47] and laser tweezers [48], where the persistent length of the muscle protein titin is between 4 Å and 20 Å [47,48,49], which corresponds to 8-40 Å in the Kuhn length according to the polymer theory [50]. R 0 is the distance cutoff where the atomic pair-wise interaction vanishes. Increasing the cutoff can in principle extract more information from protein structures and improve the accuracy. But the long distance signal may be unstable which therefore, may not be well matched by an analytical equation. By trial and error, we set R 0 = 15.5 Å as the distance cutoff in RW, which is slightly larger than 15 Å used with DFIRE and DOPE.

Conclusion
We have constructed a new transferable distance-dependent, atomic statistical potential RW, using an ideal random-walk chain of a rigid step length as the reference state. Because the ideal chain has no amino acid-specific interactions between the subunits but keeps the sequence continuity, it mimics the generic entropic elasticity and connectivity of polymer protein molecules, which could not be described by other reference states such as ideal gas systems used in DFIRE and DOPE. As a result, the RW potential has a steeper energy at short distances than these analytical potentials, which helps the RW potential to capture strong signals at short-range interactions. This is particularly important since the atomic potential in our modeling is essentially a short-range one. We also combined RW with a side-chain orientation-dependent energy term and built a hybrid potential RWplus. It is found that the orientation energy term does improve the ability of RW in recognizing the native-like structural features. RW and RWplus have been extensively tested on nine sets of structural decoys from manual assembly, threading, homologous modeling, and ab initio simulations. RWplus correctly recognized the native structures in 73% of cases which is 5-15% higher than other state of the art pair-wise statistical methods. The RW potential selected better quality models than other distant-dependent statistical potentials from ROSETTA and I-TASSER simulations. When compared with a comprehensive list of publicly available other potentials, including composite potentials combining terms from multiple resources, RWplus and RW show comparable performance to the currently best quality assessment scoring functions for the decoy selections. The general correlation coefficient between the RW/ RWplus potentials and the RMSD/TM-score is 0.50-0.53 for the I-TASSER decoys which is higher than DFIRE, and significantly higher than DOPE -although the correlation coefficients for the ROSETTA decoys are slightly lower for all potentials. This strong correlation, together with and the decoy recognition power, demonstrates the exciting probability of using the potential in improving the efficiency of protein folding and protein structure refinement algorithms. The corresponding work of employing RW and RWplus to I-TASSER based ab initio protein folding is in progress.

Construction of pair-wise distance dependent potential
A variety of distance-dependent, pair-wise, statistical potentials [9,10,11,12,20] are derived from the inverse of Boltzmann's law: where k is the Boltzmann constant and T is the Kelvin temperature; R is the distance between atoms of atom type a and b; p obs a,b,R ð Þand N obs a,b,R ð Þare the observed probability and number of atom pairs (a, b) within a distance shell R to RzDR respectively; and p exp a,b,R ð Þ and N exp a,b,R ð Þ are the expected probability and number of atom pairs (a, b) in the same distance shell respectively when there is no interactions between atoms. The purpose of N exp a,b,R ð Þis to rule out by normalization the average and generic dependence of atom pairs (a, b) which do not stem from the atom-atom pair interactions. The method of counting N obs a,b,R ð Þis the same among different methods while the method of calculating N exp a,b,R ð Þis what makes one potential differ from another. Because one of the major purposes of the potential u u(a,b,R) is to recognize the correct conformations from the structural decoys generated in the structural modeling simulations where the decoys are from continuous sequences in most cases, the generic feature of the chain connectivity is a major consideration for calculating N exp a,b,R ð Þin our model.
Here, we applied the freely-jointed chain (FJC) model [50,51] to construct a random-walk reference state, which keeps the general chain connectivity but has no long-range interactions between nodes except for the entropy elasticity that is generic in all protein structures. The expected number of atom pairs at a distance shell R for the FJC can be calculated by N exp~ N N obs (a,b)P(R), where P(R) is the probability for the atom pair in a spherical shell with radius between R and R+dR and where N N obs (a,b)S P R N obs (a,b,R)T is the average number of atom pairs of type a and b in a protein molecule.
Consider a linear polymer to be a FJC with n subunits, each of the Kohn length b, which occupy zero volume so that no part of the chain excludes another, i.e. there is no interaction between the subunits (the excluded volume will be discussed later). One can regard the segments of each such chain in an ensemble as performing a random walk in the three-dimensional space. Since the atoms of distance R can be observed in the residue pairs of different order of distances along the chain, we first consider the conformation of FJC in a set of (n+1) position vectors  (Figure 4). Since the bond vectors r I i are independent of each other, the distribution function of the polymer conformation can be written as where Y r where k~D k I D.
In the small k region, it can be approximated as For the FJC with a constant bond length b, we have Sr 2 T Y~b 2 , thus Eq. (13) is a Gaussian function integration which can be explicitly carried out [52]. The probability distribution function can be written as As a function of the end-to-end distance R~DR I D, this probability distribution can be rewritten in the spherical coordinate system: The probability function of distance R for an atom pair with residue number i and i+n is the probability of the end-to-end vector R I being in the spherical shell with radius between R and R+dR if n is less than the protein sequence length N. In contrast, if n is larger than N, the probability function of distance R is zero, i.e. Given all different order of residue distances, the probability of distance R is Because the model developed here has mapped the FJC nodes to protein residues while the potential in Eq. (6) accounts for the interactions of protein atoms, there is no definite correspondence between the Kuhn length b of the FJC model and the residue scale of real proteins. Therefore, we consider l~b 2 as a freely-tuned parameter to match the scale of the FJC with that of a real protein chain. The tuning of this parameter can also partially amend the generic excluded volume interactions of the protein chain which have not been considered in the derivation of the ideal FJC model. Thus, the final statistical potential equation is Suppose u u(a,b,R)~0 at certain distance R0, the potential can be rewritten as where R 0 is the second parameter tuned for identifying the location where the atomic pair-wise interaction vanishes.

Construction of orientation dependent potential
To specify the side-chain packing orientation, we define 20 vector pairs as shown in Figure 5. For each residue type except GLY and ALA, a unique vector pair is defined based on three most representative side-chain atoms. Totally 18 vector pairs are used to represent the orientation of side-chain atoms and 2 vector pairs are used to represent the orientation of main-chain atoms. The relative orientation of two vector pairs (A and B) can be expressed by three variables: two direction vectorR R AB andR R BA and a torsion angle V as shown in Figure 6. A is the vector pair of A 0 A 1 and A 0 A 2 , which represents the orientation of three sidechain atoms A 0 , A 1 and A 2 . B is the vector pair of B 0 B 1 and B 0 B 2 , which represents the orientation of three side-chain atoms B 0 , B 1 and B 2 .R R AB is the direction vector from A 0 to B 0 .R R BA is the direction vector from B 0 to A 0 . V is the torsion angle of A 1 A 0 B 0 B 1 . We coarse-grained the orientation space into 2704 bins for two vector pairs due to the limited amount of available protein structure data and the balance between the number of bins and the available structure data for statistical analysis [13]. As illustrated in Figure 6, the direction vectorR R BA can be coarsegrained into 26 bins based on two parameters h and Q, where h and Q are the spherical angles of vectorR R BA in the reference frame of B 0 B 1 B 2 . The definition of 26 bins is illustrated in Table 5. The direction vectorR R AB can also be coarse-grained into 26 bins in the same way. The torsion angle V is coarse-grained into four bins spanning p/2 radians each. Thus, for two vector pairs, the number of bins is 26|26|4~2704.
To calculate the total orientation-dependent packing energy, we define the packing energy for two vector pairs A and B in relative orientation space using a similar Boltzmann formula as Eq. (6): Here, OAB is the relative orientation between vector types A and If we assume that every two vector pairs have the same random orientation distribution for the reference state, we can calculate the expected number of vector pair (A, B) as: Þ is the expected probability of relative orientation OAB in the reference state. We assume that every two vector pairs have no interactions in the reference state and the three orientation variables (R R AB ,R R BA and V) are independent. They should have random distributions in orientation space and p exp O AB ð Þ can be calculated as: where p randomR R AB and p randomR R BA are the probability of a vector with random orientation in space and p random V ð Þ is the probability of a random torsion angle in four bins spanning p/2 radians each and should be equal to 0.25. p randomR R AB and p randomR R BA can be obtained by calculating the fraction of surface area for each bin in a spherical surface. The probabilities of 26 bins are listed in Table 5 Construction of hybrid distance and orientation dependent potential The hybrid potential E RWplus is composed of a distance dependent energy term E RW and an orientation dependent term E orient . Therefore the total energy can be calculated by the sum of energies of all distance pairs and vector pairs of non-consecutive residues: Here, d A,B ð Þis 1 when vector pairs A and B are in contact (at least there is one heavy atom pair with distance less than 10 Angstrom) and 0 otherwise; w orient is a weight parameter optimized against training decoy sets.
Experimental structure database for potential statistics 1,383 high-resolution experimental structures were used to calculate the statistical potential. The non-redundant protein list was constructed with the PISCES web server [53], with a percentage identity cutoff 20%, a resolution cutoff 1.6 Å , and a R-factor cutoff 0.25 Å . For the RW potential calculation, the distance cutoff is R 0 . The pair distance from 0 to R 0 was divided into bins with a bin width dR = 0.5 Å . A total of 158 residuespecific atom types, same as DOPE [12], were used.

Parameter training
The RW potential is trained on the conformations generated in the CASP7 and CASP8 experiments [54,55]. This training set includes 203 single-domain targets taken from http://predictioncenter. org/download_area/CASP8 and http://predictioncenter.org/ download_area/ CASP7. Only the decoys with full-length structures were considered and those with missed residues were removed for the convenience of potential evaluations. The final decoy set has 300 to 500 models for each target.
The RW and RWplus potential was optimized using the CASP decoys by tuning parameters l and R 0 in Eq. (19) and w orient in Eq. (23). The objective is to maximize the number of correctly selected native structures from decoys and the Pearson's correlation coefficient between the RW potential and the TMscore of the modeling decoys. When l equals to 460, R 0 equals to 15.5 Å and w orient equals to 0.1, we found that the potential has the best performance with an average Pearson's correlation coefficient with TM-score to the native structure of 0.64; due to the difficulty of the CASP decoy set, the native structure was correctly selected in only 77 out of 203 targets.

Testing structural decoy sets
Eight multiple decoy sets, including the 4-state_reduced [56], fisa [57], fisa_casp3 [57], lmds, lattice_ssfit [58], moulder [59], ROSETTA [40] and I-TASSER decoys sets, were used to evaluate the performance of the statistical potential. The first five decoy sets are available through Decoys 'R' Us [60] (http://dd. compbio.washington.edu/). The moulder decoy set by John and Sali is derived by iterative target-template alignment and comparative modeling of 20 target sequences that are remotely related to their template structures [59]; it contains 300 decoy models for each target, based on a wide range of target-template alignment accuracy (http://salilab.org/ decoys).
The ROSETTA decoy set by Baker and coworkers [40,41] contains 20 random models and 100 lowest scoring models from 10,000 decoys, which were generated for 58 small proteins using ROSETTA de novo structure predictions followed by all-atom refinement (http://depts.washington.edu/bakerpg).
The I-TASSER decoy set includes the atomic structure decoys generated for 56 non-homologous small proteins. The backbone structures were first generated by the I-TASSER ab initio modeling by Wu et al. [34], where for each protein target 12,500-32,000 conformations were taken from the trajectories of 3 lowest-temperature replicas of the simulations. Because this raw decoy set may contain redundant structures and some conformations have steric clashes, we select 300-500 nonredundant decoys for each target by iterative structure clustering [45] where one representative conformation is taken from each cluster. The selected reduced decoys are then refined by energy minimization with the OPLS-AA force field [61] using GROMACS 4.0 simulation package [62] for the purpose of removing the steric clashes and regulating secondary structure details. However, the topology of the I-TASSER decoys is not changed by the energy minimization. A full set of I-TASSER decoys is downloadable at http://zhanglab.ccmb.med.umich. edu/decoys.