text:
# Anisotropic coarse-grained statistical potentials improve the ability to identify native-like protein structures

## I Introduction

A major goal of proteomics is to determine structures of proteins rapidly. It is impractical to obtain very high resolution structures on a genome wide scale. Furthermore, for processing biological functions it is important to characterize the dynamics and thermodynamics of a network of proteins. From a computational perspective, it is currently impossible to realize these goals using all-atom molecular dynamics simulations Duan and Kollman (1998); Ferrara and Caflisch (2000). Thus, there is an urgent need to develop coarse-grained, yet reliable, models for protein structures. Construction of such models requires the determination of interaction potentials between amino acid residues. The wealth of structural data on a number of proteins in the Protein Data Bank (PDB) Berman et al. (2000) has been a source for obtaining interaction potentials Lee et al. (1999); Miyazawa and Jernigan (1999); Godzik et al. (1995).

The idea of using the frequencies of amino acid pairing to determine potential interaction parameters was first proposed by Tanaka and Scheraga Tanaka and Scheraga (1976). With the exception of a few studies Sippl (1990), most of the “knowledge-based” potentials have been obtained solely in terms of residue-residue contacts. Sippl Sippl (1990) and others have introduced an explicit distance dependence in the database-derived mean force potentials using the Boltzmann formula. This method, known as the “Boltzmann device”, assumes that the known protein structures from the PDB correspond to classical equilibrium states. The side chain - side chain (SC-SC) potentials can be related to distance-dependent probability densities $`f(r)`$ by the relation

$$U^{ij}(r)=kTln\left[\frac{f^{ij}(r)}{f_{ref}(r)}\right]$$
(1)
where $`f^{ij}(r)`$ is the probability density for a side chain of type $`i`$ to be separated by a distance $`r`$ from a side chain of type $`j`$, and the choice of the reference probability density $`f_{ref}(r)`$ is very important Betancourt and Thirumalai (1999); Sippl (1995). Sippl Sippl (1995) suggests that $`r`$ can be the distance between two atoms or some other structural parameter such as dihedral angles. Other statistical potentials developed subsequently using this approach have interpreted $`f(r)`$ as the distance-dependent probability density.

In recent years, a number of studies have evaluated the goodness of statistical potentials based on pairwise additive interactions between residues Betancourt and Thirumalai (1999); Vendruscolo and Domany (1998); Ben-Naim (1997); Tobi and Elber (2000); Tobi et al. (2000); Meller and Elber (2001); Meller et al. (2001); Gatchell et al. (2000). These studies have concluded that pairwise additive potentials, dependent only on the radial distance between residues, are inadequate for structure prediction. One of the major drawbacks of using contact potentials or potentials that only depend on the radial distance is that the relative orientation between the amino acids, which plays a role in the packing of the protein interior, is not taken into account. To account for the dense interior of the folded states of proteins it is crucial to optimize the relative orientation of the participating side chains Richards and Lim (1994). It is the purpose of this paper to examine the relevance of interaction potentials between amino acids that include the orientational dependence.

By analyzing the angular distribution of side chains around amino acids, Bahar and Jernigan Bahar and Jernigan (1996) showed that some residue pairs have specific coordination states with much higher probabilities than those expected from random distributions. Another indication that the relative residue-residue orientations are important is the recent success of the united-residue force field (UNRES) developed by Scheraga and coworkers Liwo et al. (1997a, b, 1998). The SC-SC interactions of the UNRES potentials are parameterized as van der Waals potentials with a pair-specific angular dependence included in the well depth. The UNRES potentials employ a generalized Gay-Berne potential Gay and Berne (1981); Vorobjev (1990a, b) that assumes that the interacting sites are ellipsoids of revolution placed at the center of mass of the side chains. Although the relative side chain orientations are described in a simplified manner, the success of the UNRES potentials demonstrates the importance of including orientational dependence in the side chain interaction potential.

This idea is also supported by recent studies of protein side chain packing Kussell et al. (2001) which found that current potentials can accommodate more than one rotamer for 95% of side chain positions. These studies suggest that interaction potentials that can discriminate between a large number of competing rotamer states are required.

The main goals of this paper are (a) to extract statistical information about the relative residue-residue orientations and distances in proteins using the available high resolution PDB structures and (b) to investigate the utility of the orientational information in enhancing the ability of distance dependent potentials to identify native protein folds. These goals can contribute to a better understanding of the specific residue packing of native protein structures. Motivated by the above mentioned results of Jernigan et al. Bahar and Jernigan (1996) and Scheraga et al. Liwo et al. (1997a), we approached the first goal by defining local reference frames (LRFs) that permit a more precise, quantitative description of the relative orientation of a given pair of side chains. The new LRFs are related more closely to the three-dimensional configuration of the individual side chains, rather than to their relative position with respect to the backbone or to the neighboring residues. We could extract thus novel radial and angular dependent inter-residue statistical potentials that reflect the specific distributions of distances and relative SC-SC orientations as observed in native structures of real proteins. In order to address our second goal we tested the efficacy of the new potentials by using a large database of incorrect models of real proteins Hendlich et al. (1990); Samudrala and Levitt (2000), known as decoys. These decoys are computer-generated models of protein structures, specifically designed for being used in evaluating the capacity of various potential functions to distinguish the native-like conformations from non-native ones.

Before presenting our methods and the results that we obtained, a few clarifications are necessary. First, a complete set of potentials for coarse-grained protein folding simulations should take into account interactions between side chains, interactions between side chains and peptide groups Keskin and Bahar (1998), and energetics dependent on torsional angles that define the backbone structure Liwo et al. (1997a, b, 1998). We only investigate the specific features of SC-SC interactions, as has been done before in studies of statistical potentials that are only contact or distance dependent. Secondly, we constructed our own corresponding set of distance-dependent potentials rather than employing potentials constructed by other authors. In this manner we minimized the possibility that a better parameter fitting, specific computational implementation, or other technical aspects could affect our results. Finally, as in similar studies Thomas and Dill (1996) we do not directly address the issue of the dependence of the results on the training database size or sampling problems. Instead, we used a standard, reproducible approach, by employing the set of non-homologous proteins that was used by Scheraga et al. Liwo et al. (1997a, b, 1998) for similar purposes.

## II Methods

The relative residue-residue orientations are thought to be directly related to the nature of the forces that shape the specific three-dimensional structures of proteins Bahar and Jernigan (1996); Liwo et al. (1997a). However, a quantitative approach to the statistical extraction of these orientational information from high resolution structures is still needed. The necessity of including orientation dependent interactions is established in section II.1 by computing correlation functions that probe orientational order in the PDB structures. It is shown that, for a given native state topology, the orientational packing of side chains may be decomposed into a linear combination of simple cluster geometries. However, the use of simple orientational order parameters is *not sufficient* for discriminating between basic protein architectures such as $`\alpha `$-helix or $`\beta `$-sheet. As such, we are motivated to use more detailed quantitative descriptions of relative residue orientations. In order to achieve this goal, in section II.2 we introduce definitions of amino acid dependent local reference frames (LRFs) that permit a standard description of the relative SC-SC orientations. A method for extracting the radial and angular dependent pair distribution functions is presented in section II.3, emphasizing the limits imposed on the statistical analysis by the accuracy of the available experimental database of protein structures.

### II.1 Measures of Orientational Order. The Dependence of Anisotropy on Native State Topology

The specificity of SC-SC contacts suggests that the relative residue orientations should play a significant role in determining packing in proteins. This idea, supported by the study of Bahar and Jernigan Bahar and Jernigan (1996), builds up on the more general theme of how side chains pack in the native states of proteins Richards and Lim (1994); Kussell et al. (2001). The observation that the interior of protein structures is densely packed Soyer et al. (2000); Levitt et al. (1997); Tsai et al. (1999) raises the question if one can use a simple crystal lattice description (e.g. $`sc`$, $`bcc`$, etc.) for the side chain packing. While there are many qualitative descriptions of residue packing, they are mainly based on the notion that hydrophobic interactions are responsible for globular shapes Thomas and Dill (1996) and that specific polar interactions and hydrogen bonding are important stabilizing factors. Because these studies do not address quantitatively the issue of relative side chain orientations, we used orientational order parameters (OOPs) Steinhardt et al. (1983) to asses their importance.

The Orientational Order Parameters (OOPs) were introduced Steinhardt et al. (1983) to analyze the internal structure of liquids and glasses and have been used as reaction coordinates for studying the structure of a nucleating system van Duijneveldt and Frenkel (1992); ten Wolde et al. (1995, 1996). The OOPs are defined in terms of “bonds” that connect a central particle to its neighbors. For each “bond”, one can compute the corresponding values of the spherical harmonic functions $`Y_{lm}(\theta ,\varphi )`$ that can be used as *local order* parameters

$$Q_{lm}(𝐫)Y_{lm}(\theta (𝐫),\varphi (𝐫))\text{ .}$$
(2)
where $`𝐫`$ is the position vector of the central particle and $`\theta (r)`$ and $`\varphi (r)`$ are the polar and azimuthal angles of vector $`𝐫`$ with respect to any reference frame (see, e.g. van Duijneveldt and Frenkel (1992)). To make the $`Q_{lm}`$ values representative for describing the orientational order of an entire system of particles, the *global* orientational parameters are defined as

$$\overline{Q}_{lm}\frac{1}{N_b}\underset{i=1}{\overset{N_b}{}}Q_{lm}(𝐫_i)$$
(3)
where $`N_b`$ is the total number of bonds in the system.

Both $`Q_{lm}`$ and $`\overline{Q}_{lm}`$ depend on the choice of local reference frame in which the $`Y_{lm}`$ values are calculated. Second- and third-order rotationally invariant combinations can be constructed using

$$Q_l\left(\frac{4\pi }{2l+1}\underset{m=l}{\overset{l}{}}\left|\overline{Q}_{lm}\right|^2\right)^{1/2}$$
(4)
and

$$\begin{array}{c}W_l\hfill \\ _{\stackrel{m_1,m_2,m_3}{m_1+m_2+m_3=0}}\left(\begin{array}{ccc}l& l& l\\ m_1& m_2& m_3\end{array}\right)\overline{Q}_{lm_1}\overline{Q}_{lm_2}\overline{Q}_{lm_3}\hfill \end{array}$$

where the term $`\left(\begin{array}{ccc}l& l& l\\ m_1& m_2& m_3\end{array}\right)`$ is the Wigner-3j symbol Steinhardt et al. (1983). In liquids, the ratio

$$\widehat{W}_lW_l/\left(\underset{m=l}{\overset{l}{}}\left|\overline{Q}_{lm}\right|^2\right)^{3/2}$$
(5)
is not sensitive to the exact definition of neighboring bonds of a particle. The orientational order parameters $`Q_l`$ and $`W_l`$, together with the reduced order parameter $`\widehat{W}_l`$ have specific values for a number of simple cluster geometries such as face-centered-cubic (fcc), hexagonal close-packed (hcp), body-centered-cubic (bcc), simple cubic (sc) and icosahedral (icos). If only the spherical term ($`l=0`$) dominates, then it is clear that orientational potentials are not expected to be important. The emergence of non-zero values of $`Q_l`$ and $`W_l`$ not only signify that side chain orientations are important in dense packing but may also point to the nature of orientational ordering.

All the low-order values ($`l10`$) of the parameters $`Q_l`$, $`W_l`$ and $`\widehat{W}_l`$ are sensitive to the structural details that are specific to the simple cluster geometries mentioned above Steinhardt et al. (1983); van Duijneveldt and Frenkel (1992). However, for practical reasons, monitoring the values of $`Q_4`$, $`Q_6`$, $`\widehat{W}_4`$ and $`\widehat{W}_6`$ suffices to discriminate between the different cluster geometries ten Wolde et al. (1995, 1996). The values of these four OOPs (Fig. 1a) serve as a reference to which the values computed using PDB structures can be compared.

Due to scale differences between $`Q_l`$ and $`\widehat{W}_l`$, the values of $`\widehat{W}_4`$ and $`\widehat{W}_6`$ were multiplied by 10.

We calculated the $`Q_4`$, $`Q_6`$, $`\widehat{W}_4`$ and $`\widehat{W}_6`$ parameters (shown in Fig. 1b) for three sets of 50 proteins from the same family with $`\alpha `$-helical and $`\beta `$-strand architectures listed in Table 1. For each set, we have specifically selected homologue structures to maximize the possibility of finding OOPs that have different, characteristic values for different protein architectures.

If similar OOP values were obtained for the homologue sets of hemoglobins (Hb) and myoglobins (Mb), yet different from the ones measured for immunoglobulins (Ig), it would be a strong indication that the OOPs are a sensitive measure of protein architecture. As shown next, this is not the case and, therefore, the more detailed investigation of the relative orientations of side chains is necessary. The three rows in Fig. 1b correspond to averages calculated: (a) for all the types of side chains (top), (b) only for five highly hydrophobic side chains (middle), and (c) for five hydrophilic residues (bottom). The hydrophobic side chains considered were Ile, Leu, Val, Phe and Met, while the polar residues were Asn, Gln, Ser, Thr and His. In Fig. 1b, the columns correspond to the three sets of proteins that were analyzed. The first set had fifty different immunoglobulin structures with $`\beta `$-sheet architecture (left), the second set consisted of fifty hemoglobins with $`\alpha `$-helical structures (middle) and, finally, the third set comprised fifty different myoglobin structures which also have $`\alpha `$-helical architecture (right).

While the computation of $`Q_l`$ and $`\widehat{W}_l`$ in liquids is straight forward, for their correct calculation in proteins one needs to consider the following points: (1) The size of the protein must be large enough to obtain meaningful average values. Because helices or $`\beta `$-sheets do not have “interior” points, the notion of packing itself is meaningful only for tertiary structures. Hence, the proteins must contain enough tertiary “structure” Richards and Lim (1994); (2) A meaningful cut-off distance must be selected in the identification of near neighbors for a given side chain. The three protein data sets used in this study are sufficiently large that $`Q_l`$ and $`\widehat{W}_l`$ can be easily computed to assess the degree of anisotropy in the relative residue-residue orientations. We used a neighbor “cut-off” distance that is 1.2 times the value of the position of the first peak in the radial distribution function for each structure. This definition ensured that all the side chains in the first coordination shell were counted as neighbors Steinhardt et al. (1983). In the calculations presented here the centers of the residues were taken to be the geometric centers of the heavy atoms in the side chains.

Comparison of the results in Fig. 1a and Fig. 1b shows that a simple “standard” crystallographic type of order is absent in proteins. This shows that the interior of proteins does not have a simple point group symmetry. We calculated the values of the correlation coefficients between the OOPs presented in Fig. 1b and the ones in Fig. 1a. These calculations show that a simple “standard” crystallographic type of order is absent in proteins. The interior of proteins does not have a simple point group symmetry, which would have permitted the development of a simple, quantitative description of side chains packing. In Fig. 2 are shown the correlation coefficients between the theoretical OOPs values (Fig. 1a) and the OOPs (Fig. 1b) calculated for the three sets of proteins with $`\alpha `$-helical and $`\beta `$-strand structures (Table 1) .

The absolute magnitudes of the correlation coefficients are not very meaningful when comparing small data sets, but their *relative* values can give us a quantitative indication if the order parameter values are closer to one type of cluster geometry (e.g. fcc) than to another (e.g. icos). The results are shown in Fig. 2a by averaging:(a) over all the residue types (top, *all*), (b) over the five most hydrophobic residues (middle, *H*) and (c) over the five most strongly polar side chain types (bottom, *P*). Please note that the *all* data represents averages computed over all 20 amino acids, while the *H* and *P* sets correspond to only 5 amino acids each. Therefore, the results for the *all* set are not necessarily averages of the corresponding *H* and *P* values. The three curves in each figure correspond to correlation coefficients calculated for the immunoglobulins with $`\beta `$-sheet architecture (circles), for the set of hemoglobins with $`\alpha `$-helical structures (squares), and for the $`\alpha `$-helical myoglobins (triangles). A different plot of the same correlation coefficients is also shown in Fig. 2b for the values calculated for: (a) the set of immunoglobulins (top, $`\beta `$), (b) the set of hemoglobins (middle, $`\alpha _1`$), and (c) the set of myoglobins (bottom, $`\alpha _2`$). The three curves correspond for each figure to correlation coefficients calculated by averaging over all the residue types (circles), over the five most hydrophobic residues (squares), and over the five most strongly polar side chain types (triangles).

The analysis of the results in Fig. 2 leads to the following conclusions: (1) There is no clear trend towards a single geometry associated to side chain packing. Nevertheless, we observe significant values for $`Q_4`$ and $`Q_6`$ that suggests that high order orientational correlations characteristic of fcc, bcc, or icosahedral (icos) clusters (or a combination of these) are observed in proteins with $`\alpha `$-helical and $`\beta `$-sheet architecture. (2) While there are strong correlations indicating an icosahedral packing of residues (most evident for polar residues in Fig. 2a, bottom), we also find significant contributions from the other types of cluster geometries. For example, for both $`\beta `$-strand (Fig. 2b top) and $`\alpha `$-helical (Fig. 2b middle) architectures, both the hcp and sc geometries make significant contributions. This is in accord with the observations of Bagci et al. Bagci et al. (2002) who showed that there is a general uniform distribution of residues in proteins, two-thirds being approximately fcc packed and one-third occupying random positions. We also observed high correlations between the fcc packing and the orientational order of both hydrophobic and strongly polar residues in some $`\alpha `$-helical molecular structures like in myoglobins (Fig. 2a middle and bottom and Fig. 2b bottom), but there is also a high correlation for the icosahedral character of packing in those cases. (3) There is a strong similarity in the distribution of OOPs between the $`\beta `$-sheet immunoglobulin structures and the $`\alpha `$-helical hemoglobins (Fig. 2b). However, there is considerable difference in the orientational order between the two $`\alpha `$-helical proteins (hemoglobins and myoglobins). This suggests that even within a given fold there can be variations in the precise orientational registry of the side chains. Similarly, it appears that proteins with vastly different architecture can have similar orientational order. A more detailed investigation of orientational order in a large data set of proteins is warranted.

The results presented here show that the bond-orientational order parameters are useful quantitative tools in investigating the orientational order of side chains in proteins. They also demonstrate that, for a given native state topology, the orientational packing of side chains may be decomposed into a linear combination of simple cluster geometries. The importance of orientational degrees of freedom in the packing of side chains reinforces the need for deciphering anisotropic inter-residue potentials from PDB structures. A method to extract such potentials is described next.

### II.2 The Local Reference Frames for Amino Acids

As demonstrated in the previous sections, we need to define local reference frames (LRFs) for all the types of protein side chains in order to build a quantitative description of their relative orientations. The definition of LRFs for SC-SC interactions is illustrated in Fig. 3.

The local $`Oz`$ axis is defined by the positions of the $`C_\alpha `$ and $`C_\beta `$ atoms. The local $`Ox`$ axis is defined by the positions of the $`C_\beta `$ and $`C_\gamma `$ atoms. The right handed nature of the LRFs determines automatically the direction of the $`Oy`$ axis if $`Ox`$ and $`Oz`$ are known. Once the local axes are defined, the polar angles $`\theta `$ and $`\varphi `$ and the inter-residue distances $`r`$ can be used as internal degrees of freedom that describe the relative residue-residue orientations as shown in Fig. 3. These definitions are altered for the following special cases: (1) Gly does not have $`C_\beta `$ and the local $`Oz`$ axis is defined by the bisector of the angle defined by $`N^i`$, $`C_\alpha ^i`$ and $`C^i`$, (2) both Gly and Ala do not have $`C_\gamma `$ and the local $`Ox`$ axis is defined as parallel to the direction defined by the backbone atoms $`N^i`$ and $`C^i`$, (3) for Cys and Ser the corresponding coordinates of the S and O atoms are substituted for the coordinates of the missing $`C_\gamma `$, and (4) the coordinates of the midpoint between the two $`C_\gamma `$ atoms are used to define the direction of the $`Ox`$ axis for Ile and Val. For each amino acid the interaction center ($`S^i`$) is defined as the center of mass of the heavy atoms in the side chain with the exception of Gly for which the position of $`S^i`$ coincides with $`C_\alpha `$. Sample LRFs for Ile and Val are also shown in Fig. 4.

### II.3 Finding the structure-derived potentials

From the definition of the LRFs, it follows that the SC-SC interaction potentials depend on five independent parameters that define the relative positions of two distant side-chains $`i`$ and $`j`$: $`r_{ij}`$, $`\varphi _{ij}`$, $`\theta _{ij}`$, $`\varphi _{ji}`$ and $`\theta _{ji}`$. We assume that these are the only geometrical factors that influence the two-body interactions between two residues $`i`$ and $`j`$ ($`|ji|4`$). One more independent parameter (a torsional angle around $`r_{ij}`$) can be used for a complete description of the relative orientations of the two residues in three dimensions, but it is not employed in this work because its influence is expected to be important only for short range interactions. Accounting for this sixth angular parameter would also increase dramatically the statistical requirements for the protein database that is analyzed. Assuming pairwise interactions, the distance and orientation dependent potential for the residue pair $`ij`$ is

$$U^{ij}=U_{DO}^{ij}(r_{ij},\varphi _{ij},\theta _{ij},\varphi _{ji},\theta _{ji})$$
(6)
We further assume that $`U_{DO}^{ij}`$ can be decomposed as

$`U_{DO}^{ij}(r_{ij},\varphi _{ij},\theta _{ij},\varphi _{ji},\theta _{ji})`$ $`=`$ $`U_{DO}^{ij}(r_{ij},\varphi _{ij},\theta _{ij})`$ (7)
$`+U_{DO}^{ji}(r_{ji},\varphi _{ji},\theta _{ji})`$
We use the $`U_{DO}`$ notation for the statistical potentials that are both *distance* and *orientation* dependent, and the $`U_D`$ notation for potentials that depend solely on inter-residue *distances*. The assumption in Eq. 7 on the separability of potentials is not always valid. For a system of interacting side chain pairs, described by a Boltzmann equilibrium, Eq. 7 is consistent, however, with the probabilistic relation

$`P_{total}^{ij}(r_{ij},\varphi _{ij},\theta _{ij},\varphi _{ji},\theta _{ji})`$ $`=`$ $`P^{ij}(r_{ij},\varphi _{ij},\theta _{ij})`$ (8)
$`\times P^{ji}(r_{ji},\varphi _{ji},\theta _{ji})`$
where $`P_{total}^{ij}(r_{ij},\varphi _{ij},\theta _{ij},\varphi _{ji},\theta _{ji})`$ is the probability to find a pair of interacting side chains $`i`$ and $`j`$ separated by a distance $`r_{ij}=r_{ji}`$ between their interaction centers, and with relative orientations given by the set of ($`\varphi _{ij},\theta _{ij}`$) angles in the $`LRF_i`$ frame, and by the set of ($`\varphi _{ji},\theta _{ji}`$) angles in the $`LRF_j`$ frame (see Fig. 3). Eq. 7 implies, therefore, that the relative orientations of the local reference frames $`LRF_i`$ and $`LRF_j`$ of two interacting side chains $`i`$ and $`j`$ do not depend on each other. As suggested by previous studies Sippl (1995); Thomas and Dill (1996); Bahar and Jernigan (1996), this type of independence could be expected for side chains that are separated by a large enough number of peptide bonds along the backbone. This assumption is reasonable in our analysis because only side chains that are separated by at least five peptide bonds along the protein backbone are considered. This corresponds to residues that are found on a “topological level” $`k5`$ Sippl (1990). Besides this constraint, for simplicity, we are also considering (as in Thomas and Dill (1996)) that the SC-SC interactions are independent on sequence separation (i.e. for $`k5`$, all side chains are on the same topological level).

#### II.3.1 The Boltzmann Device

The Boltzmann device assumes that the known protein structures from protein databases (such as PDB) correspond to classical equilibrium states. The SC-SC potentials can therefore be related to position probability densities $`f(r)`$ (see Eq. 1), where $`r`$ can be the radial distance or the angle between the side-chains Sippl (1995). In many studies, $`f(r)`$ can be replaced by the normalized pair distribution functions $`g(r)`$ such that

$$U_D^{ij}(r)=kTln\left[\frac{g^{ij}(r)}{g_{ref}(r)}\right]$$
(9)
for the distributions depending only on distances. We adopt a more general treatment that defines

$$U_{DO}^{ij}(r,\varphi ,\theta )=kTln\left[\frac{P^{ij}(r,\varphi ,\theta )}{P_{ref}(r,\varphi ,\theta )}\right]$$
(10)
for the distance and orientation-dependent distributions. To be consistent with previous studies, we consider the reference pair distribution functions $`P_{ref}`$ to be the corresponding radial or angular pair distributions that are obtained through an analysis of all twenty residue types. A database of non homologous proteins can be used to estimate the pair distributions and to extract amino acid specific interaction potentials that are consistent with a large set of protein structures.

An important issue that appears when using probability density functions with the Boltzmann device for constructing statistical potentials is “the problem of small data sets”. As noted by Sippl Sippl (1990), dividing the SC-SC pair frequencies by both side chain type and distance intervals results in situations when the available data is too small for conventional statistical procedures. This problem was solved by Sippl by proposing a “sparse data correction” formula that builds the correct probability densities as linear combinations between the measured data and the reference, total probability densities obtained by averaging over all twenty SC types. For the general, orientation-dependent probability densities the sparse data correction can be written as

$`P_{corr}^{ij}(r,\varphi ,\theta )`$ $`=`$ $`{\displaystyle \frac{1}{1+m^{}\sigma }}P_{ref}(r,\varphi ,\theta )`$ (11)
$`+{\displaystyle \frac{m^{}\sigma }{1+m^{}\sigma }}P^{ij}(r,\varphi ,\theta )`$
where $`P^{ij}`$ are the actual probability densities obtained from the database for the $`ij`$ pair of side chains, $`P_{corr}^{ij}`$ are the corrected probabilities, and $`P_{ref}`$ is the reference probability density. A modification introduced by the orientational dependence in our case is that the number of measurements $`m`$ becomes $`m^{}=m/sin(\theta _k)`$, as $`k`$ equiangular intervals are used for the $`\theta `$ angle. This is necessary for accounting for the azimuthal dependence of volume elements in spherical coordinates. The $`\sigma `$ parameter is a constant that controls how many actual measurements $`m^{}`$ must be observed so that both the actual probabilities and the reference would have equal weights. As in other studies, we used the value $`\sigma =1/50`$ Sippl (1990); Hendlich et al. (1990); Thomas and Dill (1996).

#### II.3.2 Residue-Residue Interaction Probability Maps

Using our proposed definitions of residue-specific local reference frames (introduced in section II.2) and the histogram extraction procedure (described in detail in Appendix A), we constructed radial and orientational SC-SC interaction probability maps. All the combinations of side chain types were investigated. From the normalized angular histograms, the resulting probability density maps for residue-residue interactions were calculated. For purpose of illustration, the data shown in all the probability density maps presented in this paper were obtained by averaging over the entire distance range for which the radial and angular histograms were constructed (2Å to 26Å), as explained above.

The orientational probability maps for interactions between Ile (one of the most hydrophobic residues) and Arg (highly hydrophilic) are shown in Fig. 4 together with a schematic representation of their associated LRFs. The LRFs in Fig. 4 are shown in the proximity of the $`C_\alpha `$, to illustrate the computational process of constructing them for actual structures. However, for calculations the LRFs are translated such that their origin corresponds to the geometrical center of the heavy atoms in the side chains (i.e. the interaction centers $`S_i`$). The orientational anisotropy is stronger around Arg than around Ile (Fig. 4), which is possibly due to the relatively longer Arg side chain. Lys, another strong hydrophilic residue, presents a relative orientational interaction probability map similar to Arg but with an even stronger anisotropic character. The similarity between the Ile-Ile and Ile-Arg maps and between the Arg-Ile and Arg-Arg is largely a reflection of the fact that the statistical data collected here does not depend on the relative orientation of the second side chain ($`SC_j`$) in the reference system of the first SC ($`LRF_i`$). This finding may be a consequence of the assumptions behind Eq. 8.

The most noticeable difference between the maps calculated in the LRFs of Ile and Arg is that the relative amplitudes of the interaction probabilities at the “poles” are reversed. In other words, the most preferred interaction loci for Ile is around its “north pole” (i.e. towards the positive $`Oz`$ axis, away from the local backbone) while for Arg is in the proximity of its “south pole”. To investigate if this feature is due to the hydrophobic/hydrophilic character, we computed average orientational interaction probability maps for three groups of residues: (a) five highly hydrophobic (H) residues Ile, Val, Leu, Phe, Met, (b) eight less hydrophobic (LH) residues Ala, Gly, Cys, Trp, Tyr, Pro, Thr, Ser, and (c) seven highly hydrophilic ($`P`$) residues His, Glu, Asn, Gln, Asp, Lys and Arg. This classification was suggested by various hydropathy scales Moran et al. (1994); Eisenberg (1984). The computed average maps are presented in Fig. 5 together with the maps for the average virtual residue “all”.

The “all-all” map shows that, for any residue there is a slightly higher probability to find more residues along the attached negative $`Oz`$ axis in their respective LRFs than in the “northern hemisphere”. The definition of the LRFs was made such that the positive $`Oz`$ direction points away from the nearest backbone atoms. Due to the finite size and the relatively compact three-dimensional shapes of most proteins analyzed here, we do expect to find less residues in that direction for higher SC-SC distances and, therefore, at larger distances from the backbone. However, the “H-all” map (Fig. 5) shows that the highly hydrophobic amino acids present the reverse situation: a higher average probability to find other residues at the “north pole” of their side chains. This arises because such residues are mostly found in the “interior” of proteins, “protected” from water. The same arguments apply to the observation that, on average, the highly hydrophilic residues have preferential coordination loci shown on the “P-all” map in the proximity of their “south pole”, as it is expected that their positive $`Oz`$ axis points more often toward the water molecules. It is reassuring that the average over the less hydrophobic residues, shown in the “LH-all” map of Fig. 5, does not reflect the existence of significant orientational preferences.

## III Results: How important is the information on relative residue-residue orientations?

We performed tests using the distance and orientation dependent statistical potentials ($`U_{DO}`$) as scoring functions and we assessed their performance on a standard database of decoys developed by Samudrala and Levitt Samudrala and Levitt (2000). The results are compared with respect to the performance of potentials dependent solely on distance ($`U_D`$). Similar decoy tests have been shown to be useful in analyzing the ability of various potential schemes to correctly recognize the native state Gatchell et al. (2000); Lee et al. (2001). There are two main aspects that need to be taken into consideration when using decoy protein structures for tests: (1) The energy of the native state should be as low as possible, and (2) The decoy structure that has the lowest root mean square $`C_\alpha `$ distance (RMSD) from the native state should also have a low energy Gatchell et al. (2000); Lee et al. (2001); Lee and Kollman (2001). For studying these aspects from a quantitative point of view, we use the Z-score for both the energy ($`Z_E`$) and the RMSD deviations of the decoy structures ($`Z_{RMSD}`$). The definition of these two types of Z-scores is illustrated in Fig. 6 where the distance-dependent energies ($`U_D`$) that we computed for the set of 654 decoys of Calbindin (3icb), is plotted as a function of the $`C_\alpha `$ RMSD, calculated with respect to the native state. The Z-score of a statistical quantity $`x`$ is

$$Z=\frac{x\overline{x}}{\sigma }$$
(12)
where $`\sigma `$ is the standard deviation and $`\overline{x}`$ is the mean of the distribution of $`x`$ values. In our case, we use $`x`$ both as the energy and the RMSD of the set of decoys. In all cases, a good scoring function will lead to a negative Z-score. In Fig. 6, we show the radial energy term computed using our method for the Calbidin decoy set (3icb in Samudrala and Levitt (2000)). The high similarity between this plot and the energies depicted in Fig. 1 of Samudrala and Levitt Samudrala and Levitt (2000) confirms that the present radial potentials are essentially the same as those used by other groups. In this figure, the circle shows the position of the native state and the diamond shows the position of the decoy structure that has the lowest energy. For a potential (or scoring function) that is efficient in discriminating the native state of a protein from a set of decoys it is expected that (1) the native state (circle) corresponds to the lowest interaction energy, and (2) the decoy with the lowest energy (diamond) should have as small an RMSD as possible. Both criteria are important and we find that both the $`Z_E`$ and $`Z_{RMSD}`$ scores defined here are useful, quantitative measures of the performance of the potential energy function. These Z-scores are proportional to the distance of the corresponding parameter of interest (depicted in Fig. 6 as interrupted lines) from the mean values of their distributions (solid lines), in units of standard deviations $`\sigma `$.

Due to the large number and diversity of decoy sets that are employed in this paper, we do not repeat the specific description of the methods used for generating each decoy set and of their names (e.g. “single”,“lmds”, “fisa”, etc.). That information is provided by the “Decoys ‘R’ Us” database (http://dd.stanford.edu) and by the corresponding publications Samudrala and Levitt (2000). Details are provided for the decoys that are relevant to the results obtained for our tests.

In Fig. 7 are compared total statistical potential values obtained for distance-only dependent potentials ($`U_D`$, squares), and for distance and orientation dependent potentials ($`U_{DO}`$, circles) for decoy structures in the “single” sets (“misfold” and “pdb\_error”) of the “Decoys ‘R’ Us” database Samudrala and Levitt (2000).

The “misfold” set contains, for each native protein, an alternative structure generated by placing the same sequences on different folds with the same number of residues (e.g. the “1bp2on2paz” notation means that the “1bp2” sequence has been placed on the “2paz” fold Samudrala and Levitt (2000)). The incorrect structures in the “pdb\_error” set are experimental structures that have been substantially re-refined or found to contain errors. The points united by dotted lines (filled symbols) correspond to incorrect structures, while the continuous lines (open symbols) join scores obtained for correct states (i.e. “native” for the “misfold” set, and “latest” for the “pdb\_error” set). The arrow in Fig. 7 emphasizes the case where the $`U_D`$ score fails to identify the native state, while $`U_{DO}`$ succeeds. Overall, the “single” decoy sets proved to be relative easy tests for both the $`U_D`$ and $`U_{DO}`$ scores. The orientation dependent $`U_{DO}`$ scores consistently result in better (i.e. smaller) values for the correct protein structures in all cases studied.

In all the tests performed, some of the decoy sets that were not appropriate for analysis were eliminated. For example, non-standard side chains were present in the decoys or in the native structures, for which no $`U_{DO}`$ potentials were constructed. Another reason for eliminating a few structures was the absence of enough heavy atom coordinates in some of their large side chains. This prevented the construction of LRFs for those side chain and, therefore, the correct estimation of their $`U_{DO}`$ potentials.

Much more challenging decoy tests than the ones presented above for the “single” decoys can be performed using the “multiple” decoy sets from the “Decoys ‘R’ Us” database Samudrala and Levitt (2000). The main goal of these tests is to identify the native structure from a set of many (i.e. tens, hundreds, or even thousands) of decoys. We discuss next the results of these tests using the $`Z_E`$ and $`C_\alpha `$ RMSD Z-scores defined as shown in Fig. 6. From the “multiple” decoy sets available in the “Decoys ‘R’ Us” database, the “lattice” set was eliminated from the analysis presented here because the new $`U_{DO}`$ statistical potentials used in this paper were constructed using a training set of real, non-homologous protein structures. The specific constraints imposed by the tetrahedral lattice topology on side chain conformations for the “lattice” decoys affect more negatively the performance of the $`U_{DO}`$ potentials than the performance of the distance-only dependent scores, $`U_D`$.

In Fig. 8 are shown energy ($`Z_E`$) and $`C_\alpha `$ RMSD Z-scores ($`Z_{RMSD}`$) calculated for the multiple decoy sets “lmds”, “fisa\_casp3”, “fisa” and “4state” Samudrala and Levitt (2000). Z-scores calculated for the distance-only dependent statistical potentials $`U_D`$ (dark bars) and for the new $`U_{DO}`$ potential (white bars) values are compared. In Fig. 8, and in all the subsequent figures presenting Z-scores, more negative values are better, and the cases in which $`U_{DO}`$ gives better results than using $`U_D`$ are emphasized using arrows. It is observed that in a large number of cases the inclusion of orientational information improves the performance of both $`Z_E`$ and $`Z_{RMSD}`$ scores. An especially interesting case is the one of the “lmds” set in which all-atom distance dependent scores were shown to perform poorly Samudrala and Levitt (2000); Samudrala and Moult (1998). In this case, when both the $`Z_E`$ and $`Z_{RMSD}`$ Z-scores are considered, the new distance- and orientation-dependent potentials $`U_{DO}`$ performed much better than the distance-only dependent $`U_D`$ in a majority of cases.

In Figs. 9-11 are shown the corresponding results for the “hg\_structal”, “ig\_structal\_hires” and “ig\_structal” decoy sets of the “Decoys ‘R’ Us” database. For the “hg\_structal” decoy sets of globins, improvements are observed in about a third of the test cases. However, for both the “ig\_structal\_hires” and “ig\_structal” decoy sets of immunoglobulins, more than 70% of tests present an improvement in the performance the energy Z-scores when $`U_{DO}`$ potentials are employed ($`Z_E(U_{DO})`$).

The results for all the sets of “multiple” decoy tests are summarized in Table 2. As mentioned above, particularly strong improvements are observed when using $`U_{DO}`$ for the “lmds” set, for which distance-dependent scoring functions were reported before to have a weak performance Samudrala and Levitt (2000); Samudrala and Moult (1998). The very good energy Z-scores ($`Z_E(U_{DO})<Z_E(U_D)`$ for more than $`70\%`$ of the decoy sets) that were obtained for the large decoy sets of immunoglobulins (“ig\_structal” and “ig\_structal\_hires”) suggest that considering the orientational dependence ($`U_{DO}`$) in the cases of proteins with a significant $`\beta `$-sheet architecture it is more important than for proteins that are mainly $`\alpha `$-helical (e.g. the “hg\_structal” sets). At the same time, the fact that $`U_{DO}`$ improves both Z-scores in about a third of the “hg\_structal” decoy tests, could be explained by the globular shape of these predominantly $`\alpha `$-helical structures. These results agree to the rather intuitive observation that details about the relative side chain orientations could play a more significant role in compact, globular proteins than in proteins with relatively extended chains. Overall, it is shown that, for a majority of test cases, the $`Z_E`$ and the $`Z_{RMSD}`$ Z-scores are improved by including the relative SC-SC orientations in the structure analysis. This demonstrates that the performance of the statistical potentials in discriminating the native state is significantly enhanced by the inclusion of orientation dependent information.

The results of the Z-score calculations for the decoy sets analyzed in this paper show that there are cases where the more detailed, distance- and orientation-dependent potentials ($`U_{DO}`$) do not provide a better performance in discriminating the native state as compared to using simple distance-dependent potentials ($`U_{DO}`$). Since both the $`U_{DO}`$ and $`U_D`$ potentials employed here were extracted from the same protein database and for the same number of radial bins, it is expected that the $`U_{DO}`$ values, while more detailed, are at the same time more prone to statistical errors than the $`U_D`$ potentials. The sparse data correction procedure ensures that reasonable values are employed for situations when little or no specific statistical data on relative residue-residue distances and orientations is available. However, this procedure alone can only ensure that realistic, accurate statistical potentials are obtained when a very large training database of non-homologue protein structures is employed. The relatively reduced size of the database of structures employed here (and in Liwo et al. (1997a)) for extracting the statistical data could be an important factor responsible for the reduced $`U_{DO}`$ performance in some cases.

The type of protein architecture and the number of residues are also factors that could influence the relative performance of the $`U_{DO}`$ and $`U_D`$ potentials. A useful quantitative parameter that is directly related to the type of protein architecture is the contact order (CO) Fiebig and Dill (1993); Plaxco et al. (1998); Dima and Thirumalai (2002) of a protein structure. The CO definition employed here is

$$CO=\frac{1}{N_c}\underset{<i,j>}{\overset{N_c}{}}\mathrm{\Delta }S_{i,j}$$
(13)
where $`N_c`$ is the total number of contacts and $`\mathrm{\Delta }S_{i,j}`$ is the sequence separation, in residues, between contacting residues $`i`$ and $`j`$. Two residues are considered to be in contact when $`|ij|>1`$ and any of the heavy atoms of residue $`i`$ is within 3.75 Å of any of the heavy atoms of residue $`j`$ Shen et al. (2002). When the CO is small the protein structure presents mostly local contacts (i.e. $`|ij|<6`$), and when CO is larger, the contacts are nonlocal Fiebig and Dill (1993); Dima and Thirumalai (2002). As shown in Fig. 12a we observed an inverse proportionality between the CO values calculated for all the native states of the decoys analyzed in this paper and the fraction of local contacts for various protein architectures. Based on this observation, we investigated the relative performance of the $`U_{DO}`$ and $`U_D`$ potentials as a function of CO values (Fig. 12b) and sequence length (Fig. 12c). In Figs. 12b and 12c is shown the dependence of the difference $`\mathrm{\Delta }Z_E=Z_E(U_{DO})Z_E(U_D)`$ for the energy Z-scores on contact order values (CO) and sequence length of the native state (N). Negative $`\mathrm{\Delta }Z_E`$ values correspond to better performing scores for the distance- and orientation-dependent potentials ($`U_{DO}`$). The lines represent linear fits. While the trends are not very strong, it is observed that, for the energy Z-scores the novel $`U_{DO}`$ statistical potentials perform better then the $`U_D`$ potentials for longer proteins, presenting large contact orders, such as the $`\beta `$-sheet structures from the “ig\_structal” and “ig\_structal\_hires” sets. It is observed that the inclusion of the information on relative SC-SC orientations makes the current version of the $`U_{DO}`$ statistical potentials more useful than $`U_D`$ for relatively large, globular proteins ($`N>100`$ residues), with a significant content of $`\beta `$-sheet structure and nonlocal contacts. At the same time, it is also observed that even for relatively short ($`N<100`$ residues) and $`\alpha `$-helical proteins (low CO values), both the $`Z_E`$ and $`Z_{RMSD}`$-scores can be improved when using $`U_{DO}`$ potentials but for a smaller number of cases. These results suggest that details about side chain - backbone interactions should be included in statistical potentials for short or $`\alpha `$-helical proteins with a high content of local contacts.

## IV Conclusions

Our results demonstrate that it is possible to improve the performance of energy based scoring functions by using statistical information extracted from the relative residue-residue orientations. Calculations of orientational order parameters for investigating the residue packing in proteins showed that architecture dependent packing is best described by a linear combination of simple cluster geometries (fcc, icos, hcp and bcc). This result reinforces the need for extracting orientation dependent potentials using PDB structures. We have defined a novel local reference system for each residue for quantitatively describing the relative three-dimensional orientations in a manner that is independent of the neighboring amino acids. Arguments based on experimental resolution of protein structures were used to derive the optimal bin size employed for collecting statistical data with both angular and distance dependence. The orientational information has both common and complementary significance as compared to the information that can be extracted from the relative residue-residue distances alone. The tests that we performed on a standard data base of artificially generated decoy structures suggest that this complementarity can be very important for a large class of protein structures. The novel distance and orientation dependent statistical potentials were shown to present an enhanced ability to recognize native-like protein folds, especially for larger structures with high contact orders (e.g. immunoglobulins). They should find use in constructing the next generation of coarse grained off-lattice protein simulations. These new potentials could also be instrumental in developing more efficient computational methods for structure prediction on much larger scales than it is currently possible, addressing one of the major goals of proteomics. To achieve this it will be important to also include the relative orientations between side chains and the protein backbone Keskin and Bahar (1998).

## V Acknowledgments

We thank Ruxandra I. Dima, Alan van Giessen and Tom Keyes for useful discussions. This work was supported by the National Institutes of Health R01 NS41356-01 (JES and DT), the National Science Foundation CHE-9975494 (JES), and CHE-0209340 (DT). The authors are grateful to the Center for Scientific Computing and Visualization at Boston University for computational resources. The data visualization was carried out using VMD Humphrey et al. (1996) and Matlab (The Mathworks, Inc., Natick, MA).

*

## Appendix A Histogram Extraction

Previous studies suggest that only protein structures that are determined with a resolution of 2Å or better should be used in the computation of $`g^{ij}(r)`$ Liwo et al. (1997a) from protein databases Berman et al. (2000). A resolution of 2Å for protein structures that are determined by X-ray crystallography corresponds to an accuracy of +/- 0.2Å in atomic positions Moody and Wilkinson (1990). The 2Å resolution is often good enough to accurately assign hydrogen bonds and to allow for a limited interpretation of solvent structure.

The high resolution structural data are binned for the construction of either radial or orientation dependent pair distributions $`g^{ij}(r)`$ and $`g^{ij}(\varphi ,\theta )`$. It is important, therefore, to identify the minimal radial and angular bin sizes that ensure a certain high level of confidence that the information has been analyzed correctly and that no important artifacts due to the limitations of the experimental methods have been overlooked. In building a statistical distribution, the data that is collected from protein structures Berman et al. (2000) consists of inter-atomic distances. We expect that such data has an accuracy of $`dR=+/0.4`$Å.

For the radial dependence of the data, let the bin size be $`L`$. In the process of assigning a certain pair distance to a certain bin, if that value is exactly at the bin boundary, we have a probability of correct assignment $`P_{ca}=0.5`$. On the other hand, if the value is exactly in the middle of the bin, this probability is maximum (if $`L2dR`$ than $`P_{ca}=1`$). We estimate that the probability of correct radial assignment depends on the bin size $`L`$ and the accuracy $`dR`$ as

$$P(L,dR)=2/L_0^{L/2}P_{ca}(x,dR)𝑑x$$
(14)
This estimate is based on the assumption that the probability of correct assignment to a certain bin of a certain pair distance value is directly proportional to the difference $`x`$ between that value and the nearest bin boundary, or that $`P_{ca}(x,dR)=a+bx`$. Using the above assumed conditions for the values of $`P_{ca}`$ at the bin middle and on the boundary, we find that

$`P_{ca}(x,dR)`$ $`=`$ $`(1+x/dR)/2\text{ if }x<dR`$
$`P_{ca}(x,dR)`$ $`=`$ $`1\text{ if }xdR`$ (15)
Therefore, for the average probability of correct assignment we find that

$`P(L,dR)`$ $`=`$ $`1{\displaystyle \frac{dR}{2L}}\text{ if }L2dR\text{ and}`$
$`P(L,dR)`$ $`=`$ $`0.5+{\displaystyle \frac{L}{8dR}}\text{ if }L<2dR.`$ (16)
From this result we can infer that for $`dR=0.4`$Å, if we choose a radial bin size $`L=0.5`$Å, the confidence of correct bin assignment is about 65.6%. In this work we employ a radial bin width of $`L=1.2`$Å. This bin width is at least twice larger than the values employed by other authors. However, it ensures that we have a confidence level of at least 83.3% that the radial bin assignment is correct.<sup>2</sup><sup>2</sup>2If a Gaussian probability of correct assignment, $`P_{ca}(x,dR)`$, is used instead of the linear dependence assumed above, with the same constraints ($`P_{ca}(0,dR)=0.5`$ and $`P_{ca}(dR,dR)=1`$) for the same value $`dR=0.4`$Å, the confidence level for a bin size $`L=1.2`$Å is estimated to be 87.3%, and only 71.5% for $`L=0.5`$Å.

Similar considerations are used for studying the angular orientation pair distributions. Consider the polar coordinate system and the division of the $`\theta `$ and $`\varphi `$ values. We assume that we have $`N`$ bins in $`\theta `$ and $`2N`$ bins in $`\varphi `$. The solid angle corresponding to one angular bin is $`d\mathrm{\Omega }=2\pi /N^2`$. This solid angle is proportional to the ratio between the area element that insures a certain high value for the confidence level of correct assignment, and the square of the corresponding minimum pair distance (radius) at which we expect that confidence level. Let us require a confidence level of 80% at a distance of 5Å. At larger distances the confidence level will be higher, but it will decrease at closer ranges. Since $`d\mathrm{\Omega }L^2/(R+dR)^2`$, and $`L=1.2`$Å, $`R=5.0`$Å and $`dR=0.4`$Å, we have $`d\mathrm{\Omega }0.22^2=2\pi /N^2`$. This reasoning leads us to employ the nearest even value of $`N=12`$ for the number of bins that divide the angular interval for $`\theta `$. Previous studies Bahar and Jernigan (1996) used a smaller value ($`N=3`$).

The use of a larger radial bin size (1.2Å vs. 0.5Å) than previously suggested, is needed to ensure a high confidence level (at least 83.3%) of correct radial bin assignment. Here, we employ twenty radial distance bins with $`L=1.2`$Å for distances starting at $`2`$Å. When we collect statistical data on the relative 3-dimensional orientation of pairs of residues, we obtain a similarly high confidence level by using $`N=12`$ for $`\theta `$ and $`N=24`$ for $`\varphi `$.

All the calculations presented in this section are based on the assumption that all protein structures analyzed have a resolution of 2Å or better. If proteins of different structural resolution are analyzed, the optimal bin size can be estimated accordingly. These arguments assure that optimal radial and angular bin sizes (i.e. values that are small enough to provide a good resolution, yet large enough to correspond to a high statistical confidence level) are employed.