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.
Fig. 1. CASP 8 example. aRed lines represent GDT plots of different 3D models of the target protein T0482 submitted during the CASP8
experiment. The model submitted by the group TS208 is represented in blue. b Structural alignment of the model submitted by the group TS208 (blue) and the experimental structure (red) over 67 residues evaluated by the CASP assessors
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]–12]:
(1)
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 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 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.
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
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:
(2)
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.
Fig. 4. A closer look at the GDT_TS function. Zooming into the GDT plot of the model highlighted in Fig. 1. What appears to be the graph of a continuous function is, in fact, a set consisting
of many separated vertical line segments
Definition 1
A protein structure a is a sequence of points in the three dimensional Euclidean space
(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 be a protein structure consisting of n amino acids, let be a 3D model of a, and let be a positive constant. The Hubbard function (or GDT function) is the function , defined by , where denotes the Euclidean norm on and ? is a rigid transformation (a composition of a rotation and a translation).
Theorem 1
H b is a stepwise function with finitely many steps, 1 ? k ? n ? 1.
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.
Algorithm for GDT_A
The area under H is the sum of the areas of the rectangular regions (? i
)(? i
 ? ? i?1
):
(4)
where ?0
 = 0 and (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(n7
) 7]. Hence, we resort to using the Riemann sums to approximate (instead of to compute
exactly) the area under the graph of H.
Fig. 5. The general shape of the Hubbard function. Notice that the values ? i
along with the function values H b
(? i
), , uniquely determine the area under the graph of H b
. At the biannual CASP experiment,
The following definition and an accompanying theorem can be found in virtually any
mathematical analysis textbook.
Definition 3
If is a function then where a = x0
  x1
  ···  x n
 = b is the partition of the interval [a, b] and v i
denotes the supremum of f over [x i?1
, x i
], is called the upper Riemann sum of f on [a, b].
Theorem 2
Let f be a real, non-decreasing, Riemann integrable function on an interval [a, b]. Then
(5)
where
(6)
is the upper Riemann sum of f the and
Observe that, since H b
is piecewise continuous, it must be integrable on . Thus, the area under the graph of H is
(7)
To approximate Area with a Riemann sum, one can define the partition points , where (Fig. 6) and then compute an estimate Area() of Area as
(8)
Fig. 6. Approximation of the Hubbard function by the Riemann sum. Area() is the sum of the areas of all rectangular regions
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 .
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 H i
of H(i), satisfying , . We then compute an estimate of Area as
(9)
Since , it follows that is a -approximation of Area. Below we show how to compute all H i
’s, and, in turn, in time , 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 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 to denote a “near-optimal†transformation and the corresponding set of residues will
be denoted by . Observe that .
Building upon any procedure for computing , one can develop an algorithm for by substituting for H i
in (10), where . Several existing methods can be modified and made suitable for finding . The most efficient such method relies on the concept of “radial pair†5].
Definition 4
Let be a set of points in the three-dimensional Euclidean space. An ordered pair of points
is called a radial pair of S if s j
is the furthest point from s i
among all points in S.
Theorem 3
Let T1 and T2be two transformations and let be a radial pair of S. Ifand then there exists a rotation R around the line through T1
(s k
) andsuch that, 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 near-optimal transformation 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 neighborhoods of T ?
(b k
) and T ?
(b l
), respectively, and R is the rotation around the radial axis 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 n2
possible pairs of residues in b. For each candidate radial pair , the algorithm generates a finite, representative set of transformations that map
b k
and b l
into 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 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 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, and . This can be achieved by partitioning into small cubes of side length slightly smaller than and then collecting the vertices of the cubes that are inside the open ball of radius
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 and at least one of them must be at distance 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 , where S(a, r) and B(a, r) denote the sphere and the open ball in 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 from T ?
(b l
) (Fig. 7). We note that size of A l
(a k
) is . Hence, the total number of candidate pairs of points (T(b k
), T(b l
)) is .
Fig. 7. Discretizing the space of rigid body transformations. 2D illustration of A k
(the set of the vertices of the squares shown on the left) and the set A l
(a k
) generated for one a k
 ? A k
An obvious to compute is to run the just described algorithm m times in succession, for . 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 , are inspected again during the procedure for finding , for each j  i.
We show that all transformations and the corresponding values H1
, …, H m
can be computed, at once, during the procedure of finding the last transformation,
namely . As demonstrated in the pseudocode above, the transformation T is generated only once for each pair of points and a sweep-plane algorithm for finding R is called only once for each i satisfying and . 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 The number of iterations of the second for loop matches the number of pairs of grid points around a k
and a l
, which is . Each one of iterations of the third for loop calls a O(nlogn) plane-sweep procedure to compute an optimal rotation and (if needed) to update the
value H i
. Hence, the asymptotic time complexity of the three nested for loops is