A polynomial time algorithm for computing the area under a GDT curve

Background Progress in the field of protein three-dimensional structure prediction depends on the development of new and improved algorithms for measuring the quality of protein models. Perhaps the best descriptor of the quality of a protein model is the GDT function that maps each distance cutoff θ to the number of atoms in the protein model that can be fit under the distance θ from the corresponding atoms in the experimentally determined structure. It has long been known that the area under the graph of this function (GDT_A) can serve as a reliable, single numerical measure of the model quality. Unfortunately, while the well-known GDT_TS metric provides a crude approximation of GDT_A, no algorithm currently exists that is capable of computing accurate estimates of GDT_A. Methods We prove that GDT_A is well defined and that it can be approximated by the Riemann sums, using available methods for computing accurate (near-optimal) GDT function values. Results In contrast to the GDT_TS metric, GDT_A is neither insensitive to large nor oversensitive to small changes in model’s coordinates. Moreover, the problem of computing GDT_A is tractable. More specifically, GDT_A can be computed in cubic asymptotic time in the size of the protein model. Conclusions This paper presents the first algorithm capable of computing the near-optimal estimates of the area under the GDT function for a protein model. We believe that the techniques implemented in our algorithm will pave ways for the development of more practical and reliable procedures for estimating 3D model quality.


Background
Advances in the area of protein three-dimensional structure prediction depend on the ability to accurately measure the quality of a protein model. One of the most popular and most reliable measure of the protein model quality is GDT_TS. It is defined as the average value of GDT_P θ computed for four distance cutoffs θ = 2 i , i = 0, 3, where GDT_P θ is the percentage of model residues (represented by their C α atoms) that can be placed under θ ångströms from the corresponding residues in the experimental structure [1,2]. In a "high-accuracy" version of GDT_TS, denoted by GDT_HA, the distance cutoffs are cut in half (θ = 2 i , i = −1, 2) [3]. In both approaches, the underlying assumption is that the experimental (crystallographic or NMR) structure is close to the real (native) structure (which is sometimes not true due to experimental errors).
Several methods exist for computing GDT_TS. The LGA algorithm [4] can estimate GDT_TS quickly, but those estimates deviate from the true GDT_TS values in about 10 % of the cases [5]. Rigorous algorithms for computing GDT_TS have also been developed [6][7][8][9], but they are computationally much more expensive.
The GDT_TS is commonly interpreted as an approximation of the area under the GDT curve, denoted by GDT_A [10][11][12]. Unfortunately, since the measure is approximated using the GDT function values at only several distance cutoffs, the errors in the area approximation are large. As we demonstrate later, GDT_TS is not only overly sensitive to small but also insensitive to large changes in the protein model's coordinates.
In this paper, we present a polynomial time algorithm for computing GDT_A. Our method runs on the order O(n 3 ), where n represents the length of the protein model (and O hides the log factor). The algorithm returns "nearoptimal" GDT_A scores, meaning that the errors in our estimates can be made arbitrary small i.e., smaller than any upfront specified vale. Although our method is theoretical, we believe that its parallel implementations, coupled with carefully designed speed up techniques, can result in a practical and widely used software tool.
The rest of this paper is structured as follows. First, we present three examples that illustrate drawbacks of GDT_TS and advantages of GDT_A. Then, we place our theory on a firm mathematical ground, which enables us to formally define the GDT_A computation problem. Finally, we describe the actual algorithm for GDT_A and provide its running time analysis.

Definition of the GDT function
The GDT function is a mapping that relates each distance cutoff θ to the percentage of model residues that can be placed at distance ≤θ from the corresponding residues in the experimentally determined structure. The graph of a GDT function provides a valuable insight into the quality of a protein model (Fig. 1). More specifically, the closer the graph runs to the horizontal axis (in other words, the smaller the area under the graph), the better the model.
As a single numerical measure of the model quality, GDT_TS is extensively used at CASP to rank different models for the same target [13,14]. Since it represents the average of GDT_P θ at several distance cutoffs, GDT_TS is often viewed as an approximation of the area under the GDT curve (GDT_A) [10][11][12]: However, as we demonstrate below, such a sparse sampling of the values of GDT function compromises the reliability of GDT_TS.
In our first example, we analyze the protein model for the target T0482, submitted by the group TS208 at CASP8 (Fig. 1). The GDT_TS score of this particular model was not even among the best dozen at CASP8, despite the fact that it fits the largest number of residues at distance ≤∼ 4 from the corresponding residues in the experimental structure. In fact, the blue model (Fig. 1a) can be superimposed onto the experimental structure so that all of its residues are at distance ≤ 8 from the residues in the experimental structure ( Fig. 1b), while no such superposition exists for any other model, even for the distance cutoff of 10Å. Interestingly, according to the MAMMOTH algorithm [15], the blue model is the best model for this particular target, while the DALI [16] algorithm ranks it as the second best.
Although it is impossible to tell whether #13 GDT_TS rank is more or less fair than #1 and #2 rank assigned by MAMMOTH and DALI, respectively, it is also not difficult to see that the ranking by the area under the GDT plot (GDT_A) would serve as a good compromise between these extremes. The next two examples illustrate further disadvantages of GDT_TS. As seen in Fig. 2, better GDT_TS scores can be assigned to obviously worse models. Moreover, as demonstrated in Fig. 3, very similar models can have significantly different GDT_TS scores.

Mathematical formalism
Strictly speaking, the GDT function is not well-defined. Zooming into the plot of the model highlighted in Fig. 1a, we see a set of many small vertical segments, meaning that each point on the horizontal axis is mapped to zero or more points on the vertical axis (Fig. 4). On the other hand, the inverse function (mapping each distance cutoffs θ to the percentage of residues in the model structure that can be fit under the distance θ from the corresponding residues in the experimental structure) is obviously well defined. This allows us to define the area under the GDT plot as the complement of the area under the inverse function: where Total_Area represents the area of the rectangular region under consideration (100 × 10). We start our mathematical formalism by first defining a protein structure.
(2) GDT _A = Total_Area − GDT _A Definition 1 A protein structure a is a sequence of points in the three dimensional Euclidean space R 3 The sequence elements a i can represent individual atoms, but it is more typical (in particular in protein structure prediction experiments) to assume that each point a i corresponds to the alpha-carbon atom of the protein's ith amino acid.
In what follows, we formally define the GDT function [17]. For simplicity of presentation, we will modify the codomain of GDT to represent the "fraction of residues" (ranging from 0 to 1) instead of "percentages of residues" (ranging from 0 to 100). We note that this simple rescaling of the ordinate values will have no effects on the results obtained in our study.

Definition 2
Let a = (a 1 , . . . , a n ) be a protein structure consisting of n amino acids, let b = (b 1 , . . . , b n ) be a 3D model of a, and let θ be a positive constant. The Hubbard function (or GDT function) is the function where denotes the Euclidean norm on R 3 and τ is a rigid transformation (a composition of a rotation and a translation).
Proof Since H b is monotony non-decreasing and since the range of H b is a finite subset of (0,1], it follows that H b must be a stepwise function. To complete the proof, we note that the number of steps in H b matches the size of its range, which does not exceed n − 1, where n is the length of b. For simplicity of presentation, from now on (and whenever the model b is implied), we will omit the subscript in H b and denote the Hubbard function only by H.
(3) a = (a 1 , . . . , a n ). Fig. 2 Insensitivity of GDT_TS. This theoretical example shows no sensitivity of GDT_TS to large variations in model quality. Surprisingly, the red model has a better GDT_TS score than the better blue model, even though it is worse by all standards. Notice that, unlike GDT_TS, the GDT_A measure is not skewed by the values at the cutoff points 1, 2, 4 and 8 Å. In fact, the GDT_A score of the blue model is twice as good as that of the red model Fig. 3 Oversensitivity of GDT_TS. a A four helix bundle-like (toy) protein (dashed grey line) along with two of its, almost identical, models (red and blue). A realistic example of such a target protein (PDB ID: 1JM0A) is shown on the right (b). In this example, we assume that the protein and its models are extended to the right to include 100 or more residues. Note that, if d ∈ {1 Å, 2 Å, 4 Å, 8 Å} then the GDT score of the blue model is significantly higher than that of the red model. For instance, if d = 2 Å, then the blue model has the GDT_TS score of about 87.5 since ~50 % all of its residues can be fit at distance ≤1 Å and 100 % under each distance 2, 4 and 8 Å from the corresponding residues in the experimental structure (dashed grey). On the other hand, the GDT_TS score of the red model is only about 75, since only ~50 % of the red model' s residues can be placed under 1 and 2 Å and 100 % under 4 and 8 Å. In fact, no matter how close the red model gets to the blue model, its GDT_TS score will never improve. Note also that the blue and red models have almost identical GDT_A scores, since GDT_A is not sensitive to small coordinate changes

Algorithm for GDT_A
The area under H is the sum of the areas of the rectangular regions (θ i )(θ i − θ i−1 ): where θ 0 = 0 and θ k+1 = θ (Fig. 5). It would be trivial to compute Area had we known all θ i and all function values H(θ i ). Unfortunately, even if we knew the step points θ i , it would be computationally very difficult to compute the function values at them, since the best to date algorithm for computing H(θ i ) runs on the order of O(n 7 ) [7]. Hence, we resort to using the Riemann sums to approximate (instead of to compute exactly) the area under the graph of H.
The following definition and an accompanying theorem can be found in virtually any mathematical analysis textbook.

Theorem 2 Let f be a real, non-decreasing, Riemann integrable function on an interval [a, b]. Then
where is the upper Riemann sum of f the and Observe that, since H b is piecewise continuous, it must be integrable on [0, θ]. Thus, the area under the graph of H is To approximate Area with a Riemann sum, one can define the partition points ǫ, 2ǫ, . . . , mǫ, where m = ⌈θ/ǫ⌉ (Fig. 6) and then compute an estimate Area(ǫ ) of Area as  Fig. 1. What appears to be the graph of a continuous function is, in fact, a set consisting of many separated vertical line segments The error |Area -Area (ǫ)| in the estimate (8) is below 2ǫ. Up to a half of this error is due to the error in the Riemann sum with the remaining error being due to the possible placement of the last partition point mǫ outside the interval [0, θ].
Unfortunately, computing the area estimates according to (8) is still a challenging problem, because (as we mentioned above), there is no computationally effective procedure for finding the function values H(iǫ). To circumvent the problem, we utilize an efficient algorithm capable of computing the lower bound estimates We then compute an estimate Area(ǫ) of Area as Since Area(ǫ) − Area(ǫ) < 2ǫ, it follows that Area(ǫ) is a 4ǫ-approximation of Area. Below we show how to compute all H i 's, and, in turn, Area(ǫ) in time O n 3 logn/ǫ 6 , where n is the length of b. Our algorithm takes advantage of an efficient procedure for computing near optimal GDT_TS values [5].
Let T(b) denotes the image of the model structure b under the transformation T. Denote by MAX(T , θ) the largest fraction of residues from T(b) that are at distance ≤θ from the corresponding residues in the experimental structure a. To find each H i , it is enough to compute a rigid body transformation T i satisfying Denote by T θ a transformation that places a largest subset b θ of residues from b at distance ≤θ from the corresponding residues in the experimental structure. Given T θ , one can easily compute b θ by calculating all n distances between the residues a i and T θ (b i ). Note that

P(T θ , θ) = H(θ).
We approximate the transformation T θ by a so-called "near-optimal" transformation i.e., a transformation that places at least as many residues from the model structure under distance θ + ǫ as the optimal transformation T θ places under the distance θ. From now on, we will use T ǫ θ to denote a "near-optimal" transformation and the corresponding set of residues will be denoted by b ǫ θ . Observe that P T ǫ θ , θ + ǫ ≥ P(T θ , θ ) = H (θ). Building upon any procedure for computing T ǫ θ , one can develop an algorithm for Area(ǫ) by substituting P T ǫ θ i , θ i +ǫ for H i in (10), where θ i = (i − 1)ǫ. Several existing methods can be modified and made suitable for finding T ǫ θ . The most efficient such method relies on the concept of "radial pair" [5]. Theorem 3 Let T 1 and T 2 be two transformations and let (s k , s l ) be a radial pair of S. If �T 1 (s k ) − T 2 (s k )� < ǫ/3 and �T 1 (s l ) − T 2 (s l )� < ǫ/3 then there exists a rotation R around the line through T 1 (s k ) and T 1 (s l ) such that �R T 1 s p − T 2 (s p )� < ǫ, for any s p in S. The rotation R can be found in time O (nlogn), where n is the size of S.
A proof of the above theorem can be found in [5]. The algorithm for finding R is fairly straightforward and it relies on the so-called plane-sweep approach [18].
The Theorem 3 implies that one choice for the nearoptimal transformation T ǫ θ is the transformation R ∘ T, where T is any transformation that maps the points b k and b l from the radial pair (b k , b l ) of b θ to the ǫ/3 neighborhoods of T θ (b k ) and T θ (b l ), respectively, and R is the rotation around the radial axis T (b k )T (b l ) that maps the remaining points from T(b θ ) to the ǫ-neighborhoods of the corresponding points from T θ (b θ ).
In search for a radial pair of b θ , the algorithm in [5] explores all n 2 possible pairs of residues in b. For each candidate radial pair (b k , b l ), the algorithm generates a finite, representative set of transformations that map b k and b l into θ + ǫ/3 neighborhoods of a k and a l , respectively (see the paragraph below for more details). For every such transformation T, a plane-sweep algorithm [18] is used to find a rotation R around the axis T (b k )T (b l ) that maximizes the number of residues from R(T(b)) that can be placed at distance < θ + ǫ from the corresponding residues in a. A finite set of transformations that map the residues b k and b l into the θ + ǫ/3 neighborhoods of a k and a l , respectively, is constructed in such a way to ensure that for at least one of those transformation T, �T (b k ) − T θ (b k )� < ǫ/3 and �T (b l ) − T θ (b l )� < ǫ/3. This can be achieved by partitioning R 3 into small cubes of side length slightly smaller than √ 3ǫ/9 and then collecting the vertices of the cubes that are inside the open ball of radius θ + ǫ/6 around a k (Fig. 7). The elements of this set, denoted by A k , are the candidate points T(b k ). The number of points in A k is O(1/ǫ 3 ) and at least one of them must be at distance < ǫ/6 from T θ (b k ) (Fig. 7). For each point a k ∊ A k , the set A l (a k ) of possible images of b l under T is computed by discretizing the spherical cap S(a, r) and B(a, r) denote the sphere and the open ball in R 3 with center a and radius r, respectively, in such a way that at least one point from A l (a k ) is found at distance < ǫ/3 from T θ (b l ) (Fig. 7). We note that size of A l (a k ) is O(1/ǫ 2 ). Hence, the total number of candidate pairs of points (T(b k ), T(b l )) is O(1/ǫ 5 ).
An obvious to compute T ǫ θ 1 , . . . , T ǫ θ m is to run the just described algorithm m times in succession, for θ = θ 1 , . . . , θ = θ m . However, such an approach results in many unnecessary repeated calculations as the area around a k and the corresponding spherical cap in the neighborhoods of a l are discretized over and over again. Moreover, all transformations T and R, generated and inspected during the procedure for finding T ǫ θ i , are inspected again during the procedure for finding T ǫ θ j , for each j > i. We show that all transformations T ǫ θ 1 , . . . , T ǫ θ m and the corresponding values H 1 , …, H m can be computed, at once, during the procedure of finding the last transformation, namely T ǫ θ m . As demonstrated in the pseudocode above, the transformation T is generated only once for each pair of points a k , a l ∈ A k × A l a k and a sweepplane algorithm for finding R is called only once for each i satisfying �a k − a k � < θ i + ǫ/6 and �a l − a l � < θ i + ǫ/3. The values of H i are updated on the fly.

Running time analysis
To analyze the algorithm's running time, we note that the number of iterations of the first for loop is equal to the number of candidate radial pairs (b k , b l ), which is

Conclusions
Estimating the quality of a protein 3D model is a challenging task. Automatically generated GDT_TS score is helpful as the first raw approximation but this measure is neither sensitive nor selective enough to be exclusively relied upon in ranking different models for the same target. In this paper, we show that using a more accurate approximation of the area under the GDT curve as the criterion of model quality addresses many of the drawbacks of GDT_TS. We also present a rigorous O n 3 algorithm for computing the area under the GDT curve for a given model, where n is the model's length. The area estimate returned by our method is "near-optimal", meaning that the error in the estimate can be made smaller than any upfront specified value.
Despite the cubic asymptotic running time with a relatively large hidden constant, we believe that the techniques presented in this paper can guide a future development of a computationally efficient computer program, in particular since our methodology is amenable to parallel implementations. A heuristic version of the algorithm for estimating the area under the GDT plot can be found at http://bioinfo.cs.uni.edu/GDT_A.html.