A conformational potential, which will be used for fold recognition, is represented as the sum of coarse-grained long-range E^l and short-range E^s potentials. The long-range potential has two terms, a contact energy E^c reflecting contact frequencies in crystal structures and a repulsive energy E^r to penalize overly dense packing

The short-range potential is a secondary structure potential based on peptide dihedral angles. All of these potentials are estimated as potentials of mean force from the observed distributions of residue-residue contacts and of peptide dihedral angles at the residue level in crystal structures of proteins. In the following, energy is represented in k_BT units, where k_B is the Boltzmann constant and T is temperature.

The total contact energy is defined here as the sum of all pairwise energies between residues,

where e^c(r_i,r_j) is the contact energy between the ith and jth residues, and r_i represents all the atomic positions of the ith residue. The pairwise energy potential is represented as the sum of two terms, one of which is the usual contact potential^2,3,4 and the other is a potential of mean force for relative orientations between contacting residues that is evaluated here from the statistical distribution of relative orientations,

where ^c(r_i,r_j) represents the degree of contact between the ith and jth residues, e is the contact energy for residues of types a_i and a_j in contact, and e^o_aiaj(r_i,r_j) is the orientational energy for the relative direction and rotation between amino acids of type a_i and a_j contact; a_i means the amino acid type of the ith residue. Here, it should be noted that the radial distance between residues is described by specifying whether or not these residues are in contact with each other, and that orientational interactions are assumed only for residues that are in contact with each other.

^c(r_i,r_j) takes a value one for residues that are completely in contact, the value zero for residues that are too far from each other, and values between one and zero for residues whose distance is intermediate between those two extremes, about 6.5 Å between geometric centers of their side chain heavy atoms. Previously, this function was defined as a step function for simplicity. Here, it is defined as a switching function as follows; in the equation below to define residue contacts, r_i means the position vector of a geometric center of side chain heavy atoms or the C atom for GLY,

where S_w is a switching function, and r is the van der Waals radius of a residue of type a which is estimated from the average volume V_a occupied by a residue of type a in protein structures with the packing density of hard sphere ; V_a are those calculated in Refs. 46 and 47 and listed in Ref. 2. A critical distance to define a residue-residue contact is about 6.5 Å, but it is taken to be larger for bulky residues.

Pairwise contact energies are defined as the sum of collapse energy e and a residue-type dependent term e; r means an average residue here.

The energies e for all pairs of the 20 types of residues were recalculated⁴⁴ from 2129 protein species representatives of the SCOP⁴⁸ Release 1.53 with the sampling method³ and with the parameters evaluated in Miyazawa and Jernigan⁴ to correct these values estimated by the Bethe approximation; actually, the estimates of contact energies corrected for the Bethe approximation are divided by 0.263 defined in Eq. (34) of that paper⁴ and used as the values of e. In other words, the intrinsic pairwise interaction energies e_ij are corrected relative to the hydrophobic energies e_ir, and the hydrophobic energies are not corrected at all; see that paper⁴ for the exact definitions of e_ij and e_ir. This scheme is employed, so that all the energy potentials in Eq. (1) have magnitudes estimated as the potential of mean force from observed distributions by assuming a Boltzmann distribution.

The collapse energy e is essential for a protein to fold by canceling out the large conformational entropy of extended conformations but it is difficult to estimate.^2,3 The value –2.55 originally estimated^2,3 for e is used here; as a result, the contact energy e takes a negative value for all amino acid pairs except for LYS-LYS pair.

In the representation of the relative location of one residue with respect to another three degrees of translational freedom and three degrees of rotational freedom must be specified. Here, distances between residues are described by specifying whether or not those residues are in contact with each other. Thus, for contacting residue pairs, two degrees of directional freedom and three degrees of rotational freedom are needed to represent those relative locations. Let us use polar angles (,) and Euler angles (,,) to describe the direction and rotation of one residue relative to another, respectively. A local coordinate system fixed on each residue will be defined later. The potential of mean force for residue orientations is defined as

where f_aa(,,,,) is a probability density function for a residue of type a at the orientation (,,,,) relative to the residue of type a; it satisfies

An obvious relationship between the Euler angles exists for the distribution of residue orientations between f_aa and f_aa:

The relationship in respect to the polar angles (,) is not simple, but (,) can be uniquely calculated from (,,,,). Thus, in principle, f_aa and f_aa must be equal to each other:

However, in the present statistical estimation of the probability density, the relationship above would be approximately satisfied. Therefore, the potential is evaluated in the form of Eq. (11).

The second and the fourth terms in Eq. (11), each of which is the orientational entropy in k_B units, are calculated as

Here it is important to note that this term represents a reference state such that the expected value of the orientational energy for each type of contacting residue pair in the native structures is equal to zero. Thus, this orientational potential represents simply the suitability of a relative orientation between contacting residues, but does not represent at all whether a contact between residues is favorable or not. The latter is supposed to be represented in the present scheme by the usual contact energy e. The reference distribution of residue-residue orientations for these orientational potentials is the uniform distribution, and not the overall distribution for all types of amino acid pairs employed by others.^33,34,35 Therefore, for residue pairs whose distributions coincide with the overall distribution, the latter potentials give always no preference but the present potentials give a preference. This is a desirable behavior for orientational potentials, because such an overall distribution of residue-residue orientations would not be an intrinsic characteristic of non-native conformations but rather of native structures of proteins.

Instead of directly evaluating the frequency distributions of relative residue-residue orientations in angular space, we estimate it with a series expansion in spherical harmonics functions. The use of spherical harmonics functions to represent orientational distributions of residue-residue pairs was attempted by Onizuka et al.³³ and Buchete et al.^34,35 The probability density is expanded as follows in the series of spherical harmonics functions which makes a complete orthonormal system with the (,,,,) variables.

g is represented as

where Y is the normalized spherical harmonics function, P is the associated Legendre function; P with m_p = 0 is the Legendre polynomial. Then, the coefficients in the expansion of Eq. (18) can be calculated from the observed density distribution by

Thus, the coefficient of the first constant term in Eq. (18) that corresponds to the uniform distribution is obvious;

Buchete et al.^34,35 employed spherical harmonics functions only for smoothing the frequency distributions of residue-residue relative orientations observed in angular coordinate space. However, to estimate the expansion coefficients, the formal representation of an observed probability function with the function can be used,³³ that is,

and then, the expansion coefficients are calculated as

where (_µ,_µ,_µ,_µ,_µ) is a set of angles observed for the contact µ between residue types a and a, and w_µ is a weight for this contact. The summations in the equations above are over all contacts of amino acid types a versus a. A contact between amino acid types a and a is counted as one half of a contact for a versus a and another half for a versus a; N_aa + N_aa is equal to the actual number of contacts between amino acid types a and a. Thus, a weight w_µ is equal to 0.5w^c, where w^c is a sampling weight for each protein that is described in the section "Datasets of protein structures used." In Eq. (24), residues are regarded to be in contact if the geometric centers of side chains or C atoms for GLY are within 6.5 Å.

The sample size limits the frequencies of those modes whose expansion coefficients can be reliably estimated. High order terms are less reliably estimated than the low order terms. Bayesian statistical analysis suggests using "pseudo counts" for expected occurrences of residue pairs.^8,49 As a result, the expansion coefficients of the observed distribution are estimated as follows:

where is taken to be

in order to reduce statistical errors resulting from the small size of samples; in Eq. (31) is a parameter to be optimized. Equation (31) means that more samples are required to determine higher frequency modes. In Eq. (27), the first term becomes more effective than the second term in the limit of small numbers of N_aa, and inversely the second term becomes more effective than the first term in the limit of large numbers of N_aa.

Then, higher order terms in Eq. (18), which tend to reflect artificial contributions from the small size of samples, are ignored by evaluating only the lower order terms with

and

where O_cutoff is a cutoff value for expansion terms.

In order to reduce the number of expansion terms, we choose only terms in the expansion whose coefficients have absolute values larger than a certain cutoff value. Thus, the probability density function is evaluated as

where H is the Heaviside step function which takes a value of one for zero and positive values of the argument and is otherwise zero. Finally the estimate of the probability density f_aa(,,,,) is cut off at sufficiently low and high values in such a way that its logarithm takes a value within an appropriate range; for example, –7–ln f_aa(,,,,) + ln(cg₀₀₀₀₀)1.

The orientational entropy defined by Eq. (17) is evaluated with the observed probability distribution of Eq. (24).

A repulsive potential used here is the one described in details in Ref. 3 to prevent packing at overly high densities; it consists of a hard core repulsion e^hc, an excess contact energy e, and a repulsive packing potential e,

where S_w is defined by Eq. (7). The repulsive packing potentials e for the 20 types of residues are estimated from the observed distributions of the numbers of contacting residues in dense regions of protein structures by assuming a Boltzmann distribution.³   N(a_i,n) is the observed number of residues of type a_i that are surrounded by n residues in the database of protein structures.   q is a coordination number, which is defined as the maximum feasible number of contacting residues around a residue, for the amino acid of type a_i.    in Eq. (40) is a small value (= 10^–6) that is added to avoid the divergence of the logarithm function. The observed distribution N(a_i,n) used here is one⁴⁴ compiled from 2129 protein species representatives of the SCOP⁴⁸ Release 1.53 with our sampling method.³

The short-range potential is evaluated here by the sum of dihedral angle dependent energies e(_i,_i) over all residues:

For this secondary structure potential, a 10° mesh over (,) space is used to count frequencies of amino acids observed in protein native structures, and this intraresidue potential e for each amino acid type a is evaluated as

where N_a(,) is the number of amino acids of type a at (,) observed in protein native structures, and N_a is their sum over the entire (,) space, that is, the number of amino acids of type a. The second term is a constant term that corresponds to a reference energy, so that the (,) energy expected for each type of residue in the native structures is equal to zero.

The observed distribution N_a(,) used here is one⁴⁴ compiled from 2129 protein species representatives of the SCOP⁴⁸ Release 1.53 with the sampling method³ used to reduce the weights of contributions of structures having high sequence identity.

To estimate the orientational potential, proteins each of which represent a different protein fold were collected. Release 1.61 of the SCOP database⁴⁸ was used for the classification of protein folds. Representatives of species are the first entries in the protein lists for each species in SCOP; if these first proteins in the lists are not appropriate (see below) to use, for the present purpose, then the second ones are chosen. These species are all those belonging to the protein classes 1–5; that is, classes of all , all , /, + , and multidomain proteins. Classes of membrane and cell surface proteins, small proteins, peptides, and designed proteins are not used. Proteins whose structures⁵⁰ were determined by NMR or having stated resolutions worse than 2.5 Å are removed to assure that the quality of proteins used is high. Also, proteins whose coordinate sets consist either of only C atoms, or include many unknown residues, or lack many atoms or residues, are removed. In addition, proteins shorter than 50 residues are also removed. As a result, the set of species representatives includes 4435 protein domains; this dataset is named here as dataset A.

The recognition power of the orientational potentials for the protein native structures is evaluated by using decoy sets, "Decoys'R'Us." ³⁹ To avoid a bias, orientational potentials to be tested are compiled from a dataset of protein structures, in which native proteins included in the decoy sets are removed; the total number of proteins is reduced to 4369; this dataset is named dataset B.

Also, to remove sampling biases that result from sequence similarities among these representative proteins, a sampling weight for each protein is determined by the sampling method based on a sequence identity matrix between sequences, which is described in Ref. 3. In other words, each of the structures having similar sequences is sampled with a weight less than 1. As a result, the 4435 protein sequences of the dataset A correspond to the effective number, 3522, of sequences and include the effective number, 1 467 302, of residue-residue contacts. The 4369 sequences in the other protein dataset B corresponds to the effective number, 3506, of sequences and include the effective number, 1 463 806, of contacts. The orientational distributions of contacting residues are evaluated in the multimeric state of the complete protein structure for each protein domain.