From Isotropic to Anisotropic Side Chain Representations: Comparison of Three Models for Residue Contact Estimation

The criterion to determine residue contact is a fundamental problem in deriving knowledge-based mean-force potential energy calculations for protein structures. A frequently used criterion is to require the side chain center-to-center distance or the -to- atom distance to be within a pre-determined cutoff distance. However, the spatially anisotropic nature of the side chain determines that it is challenging to identify the contact pairs. This study compares three side chain contact models: the Atom Distance criteria (ADC) model, the Isotropic Sphere Side chain (ISS) model and the Anisotropic Ellipsoid Side chain (AES) model using 424 high resolution protein structures in the Protein Data Bank. The results indicate that the ADC model is the most accurate and ISS is the worst. The AES model eliminates about 95% of the incorrectly counted contact-pairs in the ISS model. Algorithm analysis shows that AES model is the most computational intensive while ADC model has moderate computational cost. We derived a dataset of the mis-estimated contact pairs by AES model. The most misjudged pairs are Arg-Glu, Arg-Asp and Arg-Tyr. Such a dataset can be useful for developing the improved AES model by incorporating the pair-specific information for the cutoff distance.


Introduction
The accurate identification of inter-residue contact is a crucial step in the understanding of protein structure. The residue contacts observed in crystal structures of globular proteins are generally considered the intrinsic inter-residue interactions. Based on this commonly accepted assumption, structures from the Protein Data Bank (PDB) [1] have been used to elucidate twobody residue contact and packing potentials since 1970s [2]. Miyazawa and Jernigan developed the theory of effective interresidue energy from protein crystal structures [3,4] based on the Bethe Approximation [5,6,7,8] and quasi-chemical approximation [9,10,11,12,13]. Applying Boltzmann's law, Sippl proposed an approach to yield mean force potential of residue interactions as a function of distance [14]. In addition to the residue-distancedependence studies [15,16], the effect of relative orientations on contact energy has been investigated [17,18]. In order to include the influence of multi-residue interactions and local environmental dependence, development of tri-residue [19,20,21], four-residue [22,23,24,25,26,27,28,29,30] and secondary structure-related energy [31] have been the focus of recent research.
The contact models can be classified into two broad categories: the all atom model and the reduced representation model. In the all atom model, a pairs of residues are considered in contact if any two non-hydrogen side chain atoms (NHSA) from residues i,j are within a specified cutoff distance [56,57,58,59]. This model is expected to have accurate determination of the contact pairs [15,47,60,61]. The drawback is that it requires the knowledge of location of all the atoms on the side chains, and that is computationally expensive in structure prediction. Popular reduced representation of the side chains have been proposed through the use of C a atom [62,63], C b atoms, the centroid of amino acid and the centroid of side chain. Models with all atom main chain backbone and a single united atom for side chain have been proposed [53]. Another advanced model has been proposed with hydrogen bonds and flexible ellipsoidal side chains [64,65,66]. However, a more accurate description is required to capture atom-atom interactions in detail.
Two residue side chains are considered to be in contact if the side chain center or the C a atom distance is less than a specified, pre-determined threshold distance [2,3,4,14,25]. The influence of a residue over surrounding medium can be effectively characterized at a limit distance [67,68,69,70]. A cutoff distance of 8.0 Å has been used in multi-body potentials [27], folding rate of proteins [71] and protein stability [72,73]. Other various contacting distances have also been used for protein-folding studies. A commonly used side chain center distance threshold is 6.5 Å [3,4,18]. Spatial contact is considered to exist if C a atom pair or C b atom pair distance is less than 7.0 Å [74,75,76,77].
In order to avoid the drawbacks of arbitrary cutoff distance, Yang, et. al. proposed the parameter-free elastic network model (pfENM) to improve the estimation of B-factor [78]. Although the artificial cutoff distance are not as perfect as we expected, these convenient criteria still find their applications in many fields, especially in protein structure networks [79,80,81]. Cutoff distance is crucial to the contact degree distribution function describing the network behavior.
It is challenging to find a cutoff distance due to the variation in sizes, the preferred orientations and the anisotropy nature of the side chains. However, it can be more and more accurate as the mechanism of the contact is more and more understood. This study compares the following three models: the surface-to-surface model based on the side chain Atom Distance Criteria (ADC), the Isotropic Sphere Side chain (ISS) model and the Anisotropic Ellipsoid Side chain (AES) model using a dataset of 424 high resolution proteins from the PDB. We derive a dataset and illustrate the pairs that were wrongly estimated using the AES model for future study to improve the AES model.

Results and Discussion
It is known that side chains tend to have preferred orientations and exist as certain energetically favorable rotamers [82,83,84,85,86,87,88,89]. This anisotropic nature of the side chains proposes challenges in the determination of the contact. In general, there are side chain overlaps (as defined by van der Waals radii) in experimentally determined NMR and crystallographic protein structures [90,91]. But the number of steric clashes is low. Side chain overlaps defined by covalent radii are even less. We investigated the overlaps in the high resolution PDB structures using three side chain models. The dataset was used in our previous work [49] and it includes 424 protein structures with single-chain, higher than 1.5 Å in resolution, less than 30% sequence identity structures from the PDB that are determined using X-ray crystallography technique [1]. Some PDBs with missed NHSAs are excluded.

Residue-contact distribution for the ADC model
In the ADC model, two amino acids overlap if any pair of atoms, one from each amino acid, is within the overlap cutoff distance. Two non-overlapping amino acids are in contact if any pair of atoms, one from each amino acid, is within the contact distance.
We calculated the overall residue contact degree l n r ð Þ and the overlap degree l ovl n r ð Þ for 424 high-resolution PDBs based on ADC model. l n r ð Þ denotes the total number of contacts among all n r residues in a protein. The overlap degree l ovl n r ð Þ means the total number of side chain overlaps for a protein with n r residues. Since it is not expected to find large number of residue overlaps in the test dataset, the lower the l ovl n r ð Þ, the more accurate contact model. The total number of residues falling within the contact distance of residue i is recorded as the contact degree n cnt i ð Þ.
A ij~1 if residue i and j are in contact 0 otherwise Here n r is the total number of residues in the protein. Residue i and its surrounding neighbors form a residue-contact cluster. This cluster is related to residue i and contains n cls i ð Þ~n cnt i ð Þz1 residues, in which the residue immediately before and after i on the protein sequence are excluded.
An overall residue contact degree l n r ð Þ can be described by the size of the contact network.
l n r ð Þ provides an intuitive understanding of the compactness and overall connectivity. In the same way, overlap degree l ovl n r ð Þ can be defined to describe the side chain overlaps.
The ADC model reveals a linear relation between the contact degree l n r ð Þ and the protein length n r (Figure 1 (A)). Linear fitting formula l n r ð Þ~kn r zb reveals a spontaneous collapse character of protein structures in different sizes. If all data points are matched simultaneously, the linear relationship can be described by l n r ð Þ~0:6n r {8:79 (the lower matching line in Figure 1 (A)) with a confidence R 2 = 0.73. R 2 is the fitting coefficient of determination. The data fitting is facilitated by the Matlab fit function [92]. The relatively low fitting confidence is the result of the deviation of some 'escaping' data points.
An interesting observation is that the ADC model has the ability to classify protein structures. The 'escaping' data points, depicted in Figure 1 (A), constitute another group and have a distinct linear slope. Thus, l n r ð Þ is divided into two separate groups with obviously different slopes. Here we use the linear slope k as a criterion. To separate the data into two groups, we used the line that fit all the data points as a reference, where k all~0 :6 and b all~{ 8:79. The slope of the data point i is calculated as k i~l i n ri ð Þ{b all n ri . If k i wk all z0:2, the data point i is placed in another group. When all the 'escaping' data points are fitted as a separate group, the linear regression is l n r ð Þ~1:3n r {8:81 with a confidence R 2 = 0.93 (the upper matching line in Figure 1 (A)).
Detailed evaluation suggests that proteins with steeper increasing slopes are highly compact and can be considered dense-core proteins [93]. These well-packed structures can roughly be classified into three categories: (1) nearly-perfect globular proteins with short and flexible secondary structures; (2) proteins composed of a bundle of tightly packed alpha helixes; (3) proteins composed of curly b sheets. Figure 1 (B) shows the distribution of the overlaps in the test data set. The fact that very few proteins have overlaps in the dataset suggests that ADC is an accurate side chain model that can be used to reflect the anisotropic effect of the residue side chains. In fact, the largest number of overlaps is 4 in one protein among the entire dataset. In addition, the sparse and random data distributions suggest that systematic misinterpretation of side chain overlaps is avoided in the ADC model.
The residue contact number depends on the r gap , which is included in the definition of residue contact cutoff distance r ij cnt (Equation (9)). Here r gap is a gap distance between atom surfaces in ADC criteria (see Methods section). In order to investigate the influence of r gap on overall residue-contact distributions, l n r ð Þ and l ovl n r ð Þ are calculated for r gap~0 .5 Å , 1.0 Å , 1.5 Å and 2.0 Å , respectively. The contact cutoff distance r ij cnt increases as the r gap increases. Residue pairs with larger distance, which were considered as separated-residues pair, are included as contact pairs. Not all these extra contact pairs reflect intrinsic residue interactions. An appropriate r gap value is required to eliminate unexpected contacts. In the ADC fixed length model section, the derivation of the optimal r gap is presented. Since the overlap cutoff distance is independent of the r gap , the number of residue overlaps remain unchanged as the r gap increases.

Residue-contact distribution for the ISS model
The ISS model uses a sphere to represent the side chain. This simple model can result in spurious side chain overlap as shown in  In this study, we used the geometry center x 0 of heavy-atoms as the center of side chain.
Here x i is the coordinate of atom i. n is the number of heavy atoms.
In the ISS method, the error caused by neglecting side chain anisotropy can also be observed in the distributions of l ovl n r ð Þ, which increases linearly with respect to protein size ( Figure 1(D)).
It is not surprising that the bulky side chains, such as Trp and Arg, lead to more significantly spurious overlaps. For example, phenyl rings ( Figure 2(A)) prefer to form parallel or vertical orientations, and the spherical representation can overestimate the size of it. We also observed that the two lines fitted in the contact degree ( Figure 1(A)) appear as one line (Figure 1(C)) in a linear regression of l n r ð Þ~1:9n r {51 with a confidence R 2 = 0.99. The sensitivity to anisotropy and ability to discriminate among different structure packings are lost in the ISS model. The l n r ð Þ behavior of the ISS model with a surface-gap distance r gap = 0.0 Å appears to be equivalent to that in the ADC model with an atom surface-gap distance of r gap = 1.0 Å or r gap = 1.5 Å .
We also calculated the distributions of l n r ð Þ and l ovl n r ð Þ with r gap~0 .5, 1.0, 1.5 and 2.0 Å respectively for the ISS model (data not shown). The increase in r gap leads to a simultaneous increase in the residue-contact cutoff distance r ij cnt . As a consequence, more residue side chain pairs are considered to be within the contact range. However, l ovl n r ð Þ distributions do not change with different r gap for the reason that the overlap cutoff distance r ij ovl is independent of r gap (see Equation (15)). For all types of r gap , ISS model has significantly more overlaps than ADC model.

Residue-contact distribution for the AES model
In AES model, the residue side chain is represented as an ellipsoid with anisotropic radii in three principal dimensions. An ellipsoid collision-detection algorithm [94,95] was used to determine the side chain contact and overlap. With the anisotropic radii in three principal dimensions, the AES model is much more accurate than the ISS model. Although the AES still has false positive determination of overlaps, the number of misjudgements is less than 5% of that in the ISS model. Figure 1(E)-(F) show the distribution of l n r ð Þ and l ovl n r ð Þ calculated by the AES model. The l ovl 2n r ratio for the AES (Figure 1(F)) is much less than that in the ISS model ( Figure 1(D)). With the ellipsoid overlap criterion, more than 95% of the false residue overlaps in the ISS model have been avoided, and most of the l ovl are less than 10. This number appear to be less than that in previous works [96].
Although the side chain anisotropy is taken into consideration to some extent, the 20 types of side chain conformations are still encompassed in a quadratic surface. The bulky side chain volume will lead to an underestimation of residue side chain distances. As a result, many closely packed residue pairs are mistakenly judged as overlaps. The anisotropic ellipsoidal radii help to improve accuracy in discriminating contact-residue pairs from separateresidue pairs in the AES model. However, the AES model still encounters difficulties in assessing the difference between overlap and contact. Part of residue contacts are taken as overlaps.
We analyzed contact determination algorithms with the three models (see method section). The ISS model is the least computationally intensive, followed by ADC and then AES. The accuracy rank is ADC, AES, ISS from high to low respectively. Algorithm accuracy/computing-cost ratio suggests that the ADC model is a cost effective model with the best accuracy. AES model is the most computational method among the three because the detection of ellipsoid collision algorithm is the most intensive step and needs to be improved in the future.
An analysis of the number of overlaps determined using the three models show that less than 50% of the total pairs of ADC contact are correctly predicted by ISS model. Whereas most of the 424 proteins have more than 95% ADC contacts shared by AES model. Figure 3 shows the overlap number distribution for the 20620 pairs of residues using the AES model. It appears that AES model is successful in determination of contact for most of the pairs. However, AES fails in most of the pairs involving Arg-Glu, Arg-Asp and Arg-Tyr. Arg-Glu pair is one of the most frequently seen false positive overlaps due to their large side chains. Figure 4 (A) illustrated the ellipsoids calculated for an Glu-Arg pair. It appears that the overlapping volume is not much in this case. In another example of Asp-Arg (Figure 4(C)), the false positive overlap involves quite a lot of overlapping volume. Figure 4 show the all-atom side chain positions of residue pair Glu260-Arg286 in 1IO0 (PDB ID) and residue pair Asp49-Arg51 in 1C7K (PDB ID). It is possible to develop an improved AES model that involves pair-specific and relative orientation dependent distance criteria for more accurate representation of the side chains.

Pair-specific contact cutoff distance
A popular contact cutoff is 5 Å between two NHSA atoms. We investigate if this threshold is a good estimation for all pairs of residues in this section. In theory, the cutoff distance in the ADC model should depend on the specific radii of atoms that are in contact and the atom surface gap distance r gap . This is because the two residues can interact through different pairs of atoms. For example, the minimal distance of Val-Phe may occur either between atom pairs CG1-CE1 or CG2-CZ. The contact/overlap distances are distributions, rather than a single value (such as 5 Å ). The relation between r ij cnt distributions and the 5 Å model is an interesting topic. In addition, the method of how to estimate an appropriate r gap value for the cutoff distance is discussed in this section.
The interactions between two residues have preferred distances and orientations, rather than a random packing. For any two residues, the most frequently occurring residue distance is considered its major contact distance. As the two residues depart from, or come close to, each other, the occurrence probability decreases gradually and the interaction energy becomes relatively unstable till the contact distance increases to the upper limit, the cutoff distance. Figure 5 depicts the contact-distance distribution of Val-Phe pair. The histogram of packing distances between Val and Phe is a Gaussian-shaped function with a peak close to 3.9 Å . The peak position, i.e. the preferred contact distance, is almost independent of the cutoff distance. Only when the gap distance r gap is beyond 3 Å will a second peak arise gradually. From the linearly increasing manner and peak-position shift, we ascertained that this peak ('non-local contact') is the result of the increase of cutoff distance, rather than an intrinsic interaction between Val and Phe. We further investigated all the pairs involving Val. The peak distributions of all residue pairs containing Val confirm the steady behavior of residue contact ( Figure 6). As the gap distance r gap increases, the peak positions corresponding to the preferred Val-XXX contact distance remain constant. While the positions of 'non-local contact' peak increase linearly with respect to r gap .
The peak distributions allow us to set the cutoff distance for residue contacts. Between the two contact peaks, there is a low occurrence valley close to 5 Å ( Figure 5). The valley position provides a rough estimation of the cutoff distance. We determined the cutoff distance for all the 210 pairs of residues using the position of the valleys ( Table 1). The popular cutoff distance of r ij cnt = 5 Å appears to be effective in most of the pairs. However, some residue pairs such as Gly-XXX have complicated distributions with multiple valleys. In such cases, a larger cutoff distance will be chosen as the optimal value such that all the preferred contact distances (the stable peaks) can be included. A proper estimation of the gap distance r gap is crucial in the ADC model. The optimal r gap is expected to cover the most intrinsic contact distance between two residues, especially the major preferred contact distance. Meanwhile, the optimal r gap must be small enough to exclude the fake ''preferred contact distance'' (the linearly increased peak position in Figure 6). The appropriate r gap values are estimated statistically from the optimal cutoff distance r ij cnt . First, r ij cnt are selected from the valley positions under different gap-distance conditions ( Figure 5). The valley position does not shift drastically when the gap distance r gap increases. This stable behavior aids us to identify a statistically optimal r ij cnt . Then, from among all gap distances, there should be one critical value at which all contact-pair distances are less than the optimal r ij cnt . The critical value is the appropriate gap distance r gap .
Take Val-Phe for example, Figure 5 shows the valley position near 5 Å , i.e. the optimal cutoff distance r ij cnt = 5 Å . When the r gap varies from 0 to 8 Å , the positions of occurrence bin extend from 4 to 12 Å along the horizontal axis. At the critical gap-distance value r gap &1 Å , all bin positions are lower than 5 Å . In such cases, the optimal gap distance is r gap &1 Å . Although small fluctuations do occur, we note that the optimal gap distance for all 210 residue pair types is around 1 Å .
Instead of the fixed value of 5 Å , the optimal cutoff distances in Table 1 provide pair-specific cutoff distances. The ADC model uses the cutoff distance r ij cnt , which depends on the specific atom pairs between two side chains (see Equation (9). For the same type of residue pair, such as Val-Phe, the minimal side chain distance may occur between different pairs of atoms and hence the cutoff contact distances are usually different. The maximal and minimal cutoff-distance variations can be seen in Figure 5. No matter what cutoff distances are used, either the popular 5-Å criterion or the optimal ones in Table 1, single-value cutoff distances can hardly deal with various atom-to-atom contact cases. The 5 Å may be a good choice for two atoms surrounded by hydrogen atoms. However, it may be too large for two atoms that have no hydrogen atoms attached to them. A fixed, large cutoff-distance value is more convenient for most residue side chain contact, but overestimation can happen when the cutoff distance is adopted for heavy atom pairs in more compactly packed side chains.

Conclusion
The influence of residue side chain anisotropy has been studied for three side chain contact models. The atom distance criteria (ADC) contact model shows high accuracy in the determination of residue contact and overlaps. Protein structures can be classified as closely packed and loosely packed groups with the help of two different linear fit of l n r ð Þ by ADC method. The isotropic sphere side chain (ISS) model has systematically misjudgements in determining of both residue contact and overlaps. The residue surface distances are underestimated and more side chain overlaps emerged. With the different radii in three principle directions, anisotropic ellipsoid side chain (AES) model is more accurate than ISS in determining residue contact. The number of misjudgement is less than 5% of that in ISS method. However, AES need much more computations than ISS model. Based on the algorithm accuracy and complexity analysis, ADC model is recommended as the best all-atom side chain contact determination method. And AES is the most promising coarse-grained method.

Atom distance criteria (ADC) for residue side chain contact and overlap
Two atoms can be considered in 'contact' when they are in close interaction. Van der Waals interaction, a commonly employed close interaction, decreases rapidly with the distance between atoms. Residue contact can be defined when any two non-hydrogen side chain atoms (NHSA) from two residues are in the range of van der Waals interaction. This interaction-based contact definition are usually converted to distance-based contact definition by a cutoff distance of van der Waals interaction [56,57,58]. A popular cutoff distance between atoms from different residues is 5 Å [97]. We discuss the atom distance criteria in details in this section.
X-ray crystallography can barely resolve hydrogen atoms in most protein crystals. As a consequence, hydrogen atoms are absent in most PDB files. Thus the influence of hydrogen atoms that are attached at the NHSA has to be approximately included in determining residue-contact relations. we define the contact radius R cnt as in the following (Figure 2 (C)).

R cnt~Rvdw zd H d H : ð7Þ
Where R vdw is the van der Waal's radius of the side chain atom; Atom interactions are confined to a limited range, such as the contact radius of an atom. If the distance r ij between atom i and j satisfies the criterion r ij vr ij cnt , the atoms are considered to be in contact. Other than a predetermined fixed value, the atom-contact cutoff distance r ij cnt is calculated based on side chain atom-surface distance, which reflects the anisotropy in side chain orientations.  for different pairs of side chains, or for the same pair of side chain with different torsion angles. Generally speaking, the atom distance between different residues cannot be less than the sum of covalent radii (the disulfide bond is about 2.05 Å in length, which is almost equal to the sum of covalent radii of the S atom). Overlaps happen when one atom is within the range of the covalent volume of other atoms, i.e., r ij vr ij ovl .
Here R i ovl is the 'overlap radius' for residue side chain atom i.
R cov is the covalent radius of the atom; d H is a constant value as defined in the contact radius; and R H cov is the covalent radius of the hydrogen atom.

Isotropic sphere side chain (ISS) contact model
In many coarse-grained protein structure models, the isotropic sphere is used as a simplification of residue side chains [2,3,4,14,25]. The sphere model depends on two parameters: center position and radius. The geometry or mass center of heavyatom collections is usually chosen as the sphere-shaped side chain center (Figure 2 (D)). Although the radius of gyration, R g , is commonly used in describing the size of the residue side chain, some atoms will be located outside the range of R g . In the present study, the side chain radius is scaled to envelop all atoms. In order to determine the contact and overlap relationships between residues, effective 'contact radius' and 'overlap radius' have been proposed for sphere-shaped side chains.
The effective 'contact radius' R cnt for isotropic sphere side chain is defined as: Here r Imax is the maximal radius of all r i . The r i denotes the distance between the side chain atom i and the center of the sphere. The atom index corresponding to the maximal radius is the I max . The R Imax vdw is the van der Waals radius of atom I max ; d H and d H have the same definitions as in the ADC model. The effective 'overlap radius' is defined as: Here R Imax cov is the covalent radius of atom I max . The R H cov is the covalent radius of the hydrogen atom. Based on the contact and overlap radii of isotropic sphere side chain model, two cutoff distances are proposed: Where r ij is the distance between the center of the sphere of residue i and j; r gap is a gap distance between spherical surfaces, which can be adjusted to provide some flexibility to residueattraction interactions.

Anisotropic ellipsoid side chain (AES) contact model
Although the ISS model works well for equi-axial, spheroidal atom systems, the radius-of-gyration techniques do not retain three-dimensional anisotropic properties with regard to side chain orientations. A more general ellipsoid side chain model is proposed in this work. The residue side chain is simulated as ellipsoids with three principal axes for arbitrarily shaped atom clusters. The orientations of resulting ellipsoids are then used to study relative positions of the residue side chain. All residue NHSAs are used to calculate three ellipsoidal radii.
The principal radii of the best-fit ellipsoid are along the transformed Cartesian coordinates axes. Principal axes can be obtained from the diagonalization of matrix M. Where Þ are the coordinates of atom i.The major and minor radii (known as the principal radius [99]) are determined directly from the Eigen values l 1 ,l 2 ,l 3 ð Þof M.
Here D M ð Þ denotes the diagonalization of M based on the cyclic Jacobi method [100]; l 1 §l 2 §l 3 are three eigen values; and r 3 §r 2 §r 1 are the major and two minor semi-axes of the best-fit ellipsoid, respectively. If the atomic mass is assumed to be uniform and the side chain to have a unit mass, the eigen values are the principal moments of inertia for the ellipsoid side chain models I 1 ,I 2 ,I 3 ð Þ .
The principal radii of ellipsoid side chain are [99]: The ellipsoid orientation vector v, with respect to the reference coordinate can also be obtained from the eigen vectors [98]. Using the ellipsoid side chain model, many spurious side chain overlaps can be avoided (see Results). Some studies have been reported with regard to the detection of ellipsoid overlap [94,95]. In this study, we apply the algorithm to residue side chain contact determinations. For two ellipsoids A: X T AX~0 and B: X T BX~0, the solution of characteristic equation det lAzB ð Þ0 is used as a simple algebraic condition for the separation of the ellipsoids. Here X~x 1 ,x 2 ,x 3 ,w ð Þ T , where w is the 4 th dimension that represents the constant term in the ellipsoid formula; A and B are 4|4 real, symmetric matrices. The interiors of two ellipsoids are represented by X T AXv0 and X T BXv0. Then, Matrix A and B can be constructed with the ellipsoid principal direction vectors v i (v' i for B) and principal axis radii r i (r' i for B) [94].
Here r 3 §r 2 §r 1 are the major and two minor semi-axes of the best-fit ellipsoid; r gap is a gap distance between ellipsoid surfaces, which represents the decay zone of attraction interaction. If the residue side chain distances clear the first two screening conditions, there still are three possibilities: overlap, contact or separated. According to the ellipsoid collision conditions, overlap can easily be sorted out. The problem is with regard to how contact from separated cases can be discriminated. In a similar manner as in the ADC and ISS models, 'contact radius' are proposed for the ellipsoid side chain model. With this contact radius, the ellipsoids can be scaled to include all atoms represented by van der Waals radius. Then, the ellipsoid collision conditions are checked for the enlarged ellipsoids. If the scaled ellipsoids are still separated, the two residue side chains are separated. Otherwise, the residue pair is said to be in contact.

Analysis of algorithm complexity
For a protein chain with N amino acids, the number of NHSA of residue k is m k ð Þ. The ADC side chain contact algorithm needs to calculate the distances between two heavy atoms i andj. Let t d be the complexity of distance operation x i {x j . The total complexity of ADC for a whole protein chain is: In the case of the ISS model, there are three main steps. The bulky spherical centers have to be estimated first. Then the sphere radii are determined. Finally, the distance between two side chains is calculated and checked. If the geometrical center is considered the side chain center x 0 , the coordinate-averaging operation x i will be involved in calculations for residue k.
Here x i is the location of the non-hydrogen atom i. The largest atom-to-center distance r max i in the side chain k is utilized as isotropic sphere radius. In order to determine the r max i , distance operation x i {x 0 k kis carried out for all m k ð Þ atoms. Finally, the distances between any two side chain centers are calculated and checked.
If t a is the approximate complexity of each add operation for x i m k ð Þ , then m k ð Þ : t a is the complexity of the coordinate-averaging x i . Let t d be the complexity of the distance operation x i {x j , and the total complexity of the ISS model will be: Although residues have different rotamers, the side chain radius will not change too much for such conformational isomers. To simplify the process, the same type of residues is assumed to have the same radius. As a consequence, the calculation of radii is only necessary for 20 types of amino acid, rather than for all the N residues. The complexity can be re-written as: Here m m 20~1 20 The complexity of the AES algorithm is comprised of several aspects: the creation of moment of inertia matrix M, the diagonalization of M and calculation of principal semi-radii, and the determination of residue contact according to ellipsoid collision conditions.
The elements of M are calculated from the relative positions of side chain atom to side chain center. In the same way as in the ISS algorithm, side chain-center calculation complexity is derived by P N k~1 m k ð Þ : t a . Here t a is the approximate complexity related to coordinate-averaging operations. Considering the symmetry of the 3|3 matrix M, relative position estimations require m k ð Þ partial distance operations (only two coordinate axes are used) for each matrix element. Let t d be the complexity of distance operation. Matrix creation has a complexity as: The direct diagonalization of matrix M results in an algorithm complexity t diag , which covers the computing cost of principal radii vectors. The total complexity for the whole protein chain is T M2~N : t diag .
The ellipsoid-collision conditions are based on the solution of the characteristic equation det lAzB ð Þ0. The constructions of A and B need products of three 4|4 matrices and the complexity is 2 : 2 : t A . Here t A is the matrices multiplication complexity. The number of solutions can be obtained by solving the characteristic equation, which has a complexity independent of protein size and total atom number. This complexity is represented as t det . The ellipsoid collision complexity for N residues is T M34 From all the above analysis, the total complexity of the AES algorithm for an entire protein is: m m : t a z6N : m m : t d zN : t diag z4N : When the protein size N increases to a very large value (N& m m), T AES has an asymptotic approximation as N 2 2 t det .
For large proteins (N& m m), the asymptotic approximation complexity of the ADC, ISS and AES algorithms are N 2 m m 2 2 t d , N 2 2 z20 m m 20 : t d and N 2 2 t det , respectively; m m is number of NHSA with respect to all the N residues in a protein; and m m 20 is the average number of NHSA with respect to 20 types of amino acids. The difference between m m and m m 20 is trivial. Thus, an obvious fact is that the ADC model needs a significantly larger number of computations than the ISS model. The AES appears less complex than the ADC model. However, the t det is much larger than t d . If t det can be written as t det~c : t d , the AES complexity will be N 2 2 ct d . The average number of NHSA usually satisfies m m&5. As a consequence, ADC complexity is around N 2 2 : 25t d . When t det §25t d , AES complexity exceeds that of the ADC model. In current algorithms, the complexity t det for solving a fourth-order equation det lAzB ð Þ0 and determining the number of different solutions is significantly greater than 25t d . Stated briefly, the AES is currently the most computationally intensive algorithm model.