Protein Functional Surfaces: Global Shape Matching and Local Spatial Alignments of Ligand Binding Sites

Surface Shape Signatures for Rapid Global Surface Shape Comparison

Unlike protein folds, a common lexicon has not emerged to precisely describe the shapes formed by surfaces yet quantitative measures can be computed to empirically describe shape. Here, we introduce a metric, the SurfaceShapeSignature (SSS), which describes the global shape of a protein surface that can be used for comparison. Adapted from three-dimensional database object retrieval, the method represents the signature of an object as a probability distribution sampled from a shape function measuring global geometric properties[24]. The complexity of shape matching is then reduced to the comparison of two probability distributions. This approach was used to quickly and successfully retrieve and classify complex shapes from three-dimensional databases.

After a protein's binding surfaces has been identified (see Methods), the SSS is constructed by systematically measuring the Euclidean distance between all unique atom pairs for a given surface. This is seen for the nicotinamide-adenine-dinucleotide phosphate (NADP) binding surface from human pathogen S. pyogenes (PDB:2ahr) in Figure 1c. The inter-atomic distances are then sorted to form the shape signatures. The SSS distributions for a selection of heme, nicotinamide adenine dinucleotide (NAD) and adenosine 5'-triphosphate (ATP) binding surfaces are shown in Figure 1d. For reference, SSS distributions for a selection of DNA and metal binding surfaces are also shown. Once the shape distributions for two surfaces are computed, we apply the Kolmogorov-Smirnov (KS) test[25] to compare the probability distributions. The KS test identifies the greatest distance between the observed and expected cumulative frequencies and is bound between zero (identical distributions) and 1 (different distributions).

Advantages and Limitations of Surface Shape Signatures

The primary advantages of this approach are computational robustness and efficiency. In the original work of Osada, the robustness of shape signatures was verified against a variety of transformations including scale, rotation and mirroring and were insensitive to model simplification. These properties are ideal for applications in structural biology as minor conformational changes or small surface perturbations should not mask overall shape similarities. For example, we should be able to detect the similarity between the apo and bound states of a binding site, yet still discern between the two states. Insensitivity to model simplification allows for meaningful comparisons between surfaces comprised of non-trivial atoms count differences. The implementation and execution of this algorithm is computationally straightforward and allows a query surface to be searched against the GPSS library in minutes.

Biochemical functions rely on a combination of shape and chemical compatibility and therefore one should not expect that shape alone could describe this complexity. Figure 1 shows that NAD, ATP, and heme binding surfaces share similar and, in some cases, overlapping distributions. The SSS comparison can convey gross shape similarity but is void of chemical typing information and cannot be used to infer specific functional information. Instead, it provides a fast and robust comparison metric utilized as the entry point into more comprehensive surface shape matching methodology.

Surface Shape Signature Similarity Threshold

The choice of a shape similarity threshold has a significant impact on the efficiency and specificity of surface library searching. Operating under the pretense that ligands with similar shape and molecular weight are more likely to bind to similar pockets, we correlated these properties to SSS KS distances. Our training data is taken from querying cAMP-dependent kinase (PDB:1atp) ATP binding surface against the GPSS library. With a correlation coefficient of 0.458, we see that some degree of surface similarity can be inferred simply by molecular weight (Figure 2a). We highlight a region on the plot (yellow) that corresponds to +/- 100 D from the molecular weight of ATP and identify the similarity distance of 0.22, along the x-axis, in which only 11% of surfaces exist as outliers. Next, ATP was compared to ligands corresponding to the surfaces in our library using the molecular shape matching application ROCS (OpenEye Scientific Software, Inc). ROCS identifies the best superposition of two molecules by optimizing their overlapping volume and reports a normalized Tanimoto score. Tanimoto values greater than 0.7 are generally regarded as having similar shape and are highlighted on our plot (cyan). The Tanimoto scores are correlated to the SSS distance in Figure 2b. We identify a distance of 0.24, in which only 10% are outliers. Since our assumptions about molecular weight and ligand shape similarity are simplistic and surface comparison hope to identify more evolutionary distant relationships, we set our default threshold distance at 0.3. In the benchmarking retrieval experiments presented in Results section, our default threshold excludes less than 1% of true-positive surfaces, all of which can be justified by unique structural incident (e.g. multiple binding pockets, mutation experiments, low resolution structure, crystallographic error).

Local Spatial Surface Residue Alignments for Physicochemical Texture Similarity

The spatial arrangement of localized surface patterns and orientation of side chains are used to assess our evaluations of biochemical function complimentarity. An exhaustive enumeration and search algorithm, SurfaceAlign, is applied to detect groups of spatially conserved amino acid residue sets. In this application, the term "conserved" refers to the identical residues common to two protein surfaces. Our three-dimensional representation of a surface residue is represented by a single point located at the center of mass of all atoms identified as contributing to a particular surface.

The alignment of two surfaces is performed by combinatorially identifying the best superposition of the maximum subset of conserved residues between surfaces. For all conserved residues of each type, we enumerate all combinations and permutations to create unique coordinate sets of common residues between the two surfaces. A visual depiction of an alignment between the heme binding pockets of myoglobin from P. catodon (PDB:1mbn) and structural genomics target hemoglobin alpha-1 (PDB:1xq5) is shown in Figure 3. The geometric dissimilarity of each coordinate combination is evaluated following the methods of Umeyama[26], which is based on singular value decomposition of the correlation matrix of the coordinates to identify the least square rotational matrix, translation vector and the root mean square distance (RMSD). We utilize two variants of the RMSD: the coordinate root mean square distance (cRMSD), for atomic coordinates as represented in our residue model, and the orientation root mean square distance (oRMSD)[18]. The oRMSD is a derivative dissimilarity measure that reduces the effect of outliers on an RMSD value and simulates the conformational flexibility of amino acid side chains.

Statistical Significance of Aligned Surface Distances

RMSD calculations are sensitive to the number of data points compared, making it necessary to assess the statistical significance of raw distance measures. This is accomplished by converting calculated cRMSD and oRMSD values to a probability value (P-value) measuring the likelihood of obtaining a specific RMSD value for a solution with a given number of residues. This allows for the meaningful comparisons between alignments with differing number of common residues. Following the method described by Binkowski et al[18], random surface alignments were performed for alignment solutions of N _res residues, where 3 ≤ N _res ≤ 100. Calculations numbering 10¹⁰ were computed to construct lookup tables associating cRMSD and oRMSD scores to P-values. To minimize the inherent bias in the PDB toward particular protein families, special consideration was employed to utilize a non-redundant surface library consisting of proteins sharing less than 90% whole-protein sequence similarity.

Surface Volumes Overlap of Aligned Surfaces

While the alignment compares key residues important for shared biochemical function, the volume overlap of an alignment provides a comparison of all atoms belonging to a surface. The overlap volume is defined as the volume difference between the superimposed surfaces and is calculated from the formula:V _AB = V _A + V _B - V_A∪B

Where V _A , V _B , and V_A∪Bare the volumes of surface A, B and the superimposed construct AB, respectively. Volumes are calculated using the weighted Delauney triangulation and alpha shape methods [27–29].

A SVOT score is bounded between 0 (representing non-overlapping surfaces) and 1.0 (representing identical surfaces). The SVOT values for the spatial alignment series show in Figure 3d are 0.50, 0.52, 0.73, 0.62, 0.91, and 0.89 when viewed from left (black) to right (white). The best spatial alignment, as judged by RMSD values, does not guarantee the best SVOT score (i.e. SVOT scores are not correlated to RMSD).

The interpretation of the SVOT score is not straightforward for surfaces that have a large volume disparity. Figure 4 shows a large surface pocket on F420-0:gamma-glutamyl ligase homolog from A. fuldgidus (PDB:2g9i, a) where a sub-surface (b, forest green) is highly conserved to the GDP binding surface in GDP-binding protein from B. taurus (PDB:1tad, c). The volume overlap of the superimposed surfaces has a calculated SVOT score of 0.37, suggesting low similarity (Figure 4e, purple). In this case, the SVOT score fails to account for strong sub-surface similarity; hidden by the overall volume disparity.

Where V _a and V _b are the volumes of the surfaces formed only by the conserved residues used in the alignment solution. The lSVOT for the alignment in Figure 4g is 0.71, conveying stronger similarity.

The rSVOT is bounded between 0 and 1 and represents the fraction of a surface utilized in the alignment solution. In Figure 4, the rSVOT value of 0.52 confirms that it is a sub-pocket match. The rSVOT score can be used to automatically detect surfaces that do not have the desired properties for a particular search (e.g. excluding sub-surface matches during a library search).

A Library of Protein Functional Surfaces

An exclusion contact surface has been calculated for every heteroatom molecule associated with a protein's PDB file and organized into the Global Protein Surface Survey (GPSS). The GPSS contains three-dimensional libraries of functionally annotated surfaces from ligand, DNA, metal and peptide binding surfaces and is updated weekly to correspond to PDB deposits. Libraries are publicly accessible through a web browser[32] or via a PyMol[33] plugin. In this study, we utilize only the ligand binding surfaces of the GPSS: 113,921 members representing 5,575 unique ligands (PDB version: November 2006). For this subset of the GPSS, the average number of residues forming a surface is 12 and the average molecular weight of the bound ligands is 305. To reduce redundancy and improve search efficiency, we further limit our search library to a single ligand of each type from each protein deposit. The first ligand of each type, as described in the PDB file, was selected.

Receiver Operator Characteristic Curves

Surface retrieval benchmarking experiments are summarized in a Receiver Operator Characteristic (ROC) curve, where the sensitivity is plotted against its specificity at various significance levels of summed probabilities. In the ROC curve, the x-axis represents the false positive rate, or 1-specificity, which is calculated by as 1-TN/(TN+FP), where TN is the number of true negatives and FP is the number of false positives. The y- axis represents the true positive rate, or sensitivity, and is calculated as TP/(TP+FN), where FN is the number of false negatives. An overall performance measure of a classification test can be calculated by the area under the ROC curve (AUC)[35]. Bound between 0 and 1, an AUC of 1 is indicative of a perfectly accurate classification test, in which all true positives are distinguished from false positive. An AUC of 0.5 corresponds to a random classification test (e.g. a coin flip). The AUC is a combined measure of sensitivity and specificity.

For our SurfaceScreen methodology, the ROC curves measure our ability to accurately identify similar surfaces from a large database. A true-positive data point in our database retrieval experiments is defined when the ligand from the query surface matches the ligand from the corresponding library surface (e.g. retrieving heme binding sites when a heme binding site is used as the query). This definition should be considered conservative, as protein surfaces have the ability to bind multiple ligands and, in some cases, a false-positive prediction may indeed be a biologically relevant hit.

ATP Conformation Classification

The coordinates for all ATP molecules in the PDB were identified and extracted. Multiple occurrences of the ligand in structure deposit were included and treated as unique molecules. A pairwise, least-squares superposition of all remaining molecules was performed and RMSD values recorded in a distance matrix. Complete linkage clustering was applied to the data matrix. For clarity, we chose to discover the minimum number of conformation families that would accurately represent all ATP molecules. A range of cut values was tested and we chose to set the number of clusters to four based on manual visualization and analysis. The four conformations represent the bent, extended, and two intermediary forms of ATP.

Detection of a Convergent Heme Binding Surface

Convergent evolution presents a far more difficult challenge for annotation of proteins of unknown function. Structural genomics target, IsdG from S. aureus[49] (PDB:1xbw), shows no significant sequence similarity and does not contain the conserved N-terminal histidine or the GXXXG motif characteristic present in the heme-monooxygenases family, yet this enzyme displays classical heme-monooxygenase activity[49]. Also, while all known members of heme-oxygenase superfamily are of all α-helical fold, IsdG adopts α+β sandwich with an anti-parallel β-sheet and ferredoxin-like fold and a β-barrel at the dimer interface[49].

The structure of IsdG has a prominent pocket formed between the α-helices and beta sheets (Figure 8a). This is the largest surface pocket identified by the CASTp webserver[50]. Querying this surface against the GPSS library reveals a striking similarity to the heme binding pocket in heme oxygenase (HmuO) from C. diphtheriae (PDB:1iw0, Figure 8b). The SSS distributions have distance of 0.06 (Figure 8d). There are 19 conserved residues between the surfaces that come from diverse regions of the primary sequence (Figure 8c). The surfaces align with cRMSD P-value of 9.84 × 10^-3 (Figure 8ef) and oRMSD P-value of 5.32 × 10^-4 (Figure 8gh). Superposition of the surfaces results in gSVOT of 0.78, lSVOT 0.84, and rSVOT of 0.93 (Figure 8i). The gSVOT overlap is highlighted in Figure 8j. The SurfaceScreen score for the comparison is 1.98, which ranked fourth overall against the search library.

Several other structural homologues to IsdG have subsequently been solved in the structural genomics effort. A clustering of the putative binding sites for four additional enzymes is shown in Figure 9. Surface analysis reveals that the heme binding pocket is well conserved in protein TT1390 from T. thermopilus (PDB:1iuj) and protein BC2969 from B. cereus (PDB:1tz0) suggesting that both T. thermopilus and B. cereus can acquire iron through heme degradation. Despite overall structural similarity, the binding surface is not well conserved in ActVA-Orf66 from S. coelicolor (PDB:1lq9), a protein involved in antibiotic synthesis that is known to bind 6-deoxydihydrokalafungin (6-DHHK)[51]. Expectedly, this surface forms the most distant branch of the clustering. Although all proteins appear to function as monooxygenases they operate on very different substrates suggesting that convergent evolution may be an important driving force to evolve new functions from existing protein scaffolds. In this manner, surface analysis could be used to define a chemical structure space by interpolating between known substrates clustered on each node.

Specific Nucleotide Binding Site Retrieval

ATP is a multifunctional nucleotide associated that has been classified to catalyze 58 different reactions by the Enzyme Commission (EC). In over 300 structural complexes, ATP binding is associated with domains from 45 homologous superfamilies, some sharing less than 8% sequence identity[52]. The nucleotide is quite flexible and adopts a wide range of conformations, some in less than energetically favorable states[52]. To determine the extent that conformational variability exerts on similarity searching, we conducted retrieval experiments with query surfaces binding ATP in diverse conformations: cAMP- dependent kinase (PDB:1atp) protein kinase CK2 from Z. mays (PDB:1a6o)[53], ATP:corrinoid adenosyltransferase from S. typhimurium (PDB:1g5t)[54], PurT-encoded glycinamide ribonucleotide transformylase from E. coli (PDB:1kj8)[55]. The conformations were selected by clustering all ATP molecules by their three-dimensional shape similarity (see Methods) and are shown in Figure 10.

The AUCs calculated using the SurfaceScreen score were 79.1%, 80.1%, 83.0%, and 85.4% for ATP:corrinoid adenosyltransferase, PurT-encoded glycinamide ribonucleotide transformylase, protein kinase CK2 and cAMP-dependent kinase, respectively (Figure 10e). The extended ATP form, which is the most prominent form in the PDB, had the best retrieval rate, while the bent form ATP had the poorest. Overall, the rates underperform compared to the HIV-1 inhibitor and heme binding surface retrievals. Despite the influence of ligand conformation of surface conformations, our method appears rather tolerant to flexible ligands (and their corresponding binding surfaces) albeit at the expense of the specificity seen in more rigid molecules.

It should be noted that the retrieval rates for ATP are especially conservative, as a disproportional number of ATP binding surfaces complexed with other molecules are in the PDB. For example, there are 11 structures of Protein Kinase A from Bos taurus (PDB:1xha,1xh8,1xh7,1xh6,1xh5,1xh4,1veb,1svg,1sve,1svh) which have ATP competitive inhibitors bound. By correcting for protein kinase inhibitor complexes, the AUCs improve approximately 4% across all query surfaces. Unfortunately, there is no automated method to associate these types of complexes with their native cofactors except through literature analysis, which can be uninformative, as every structure deposit does not have a corresponding publication. The authors are developing an automated database that catalogs such natural cofactor/inhibitor relationships between structures in the PDB.

Non-specific Nucleotide Retrieval

In some proteins, ligand binding is not an exclusive event, as some proteins are able to utilize different cofactors to catalyze the same reaction. Casein kinase 2 (CK2) is a highly conserved eukaryotic serine/threonine kinase that plays a key role in various cellular processes and possesses dual-cosubstrate specificity for guanosine-5'-triphosphate (GTP) or ATP[53]. This feature, whose biological significance is not well understood, is exceptional among eukaryotic protein kinases. Querying the CK2 ATP binding surface (Figure 10b), we can retrieve GTP binding surfaces from the GPSS library with AUC of 83.7%, slightly better than ATP. Given that the two molecules differ only in their nucleoside, this result is not surprising, as we have shown ligand shape complimentarity is a strong precursor of overall surface similarity. In a comparison of the surface retrieval rates for all nucleotides binding surfaces in the PDB against the CK2 ATP binding surface, we observe trends which mirror ligand shape similarity: purine derivatives retrieval is better than the pyrimidines and tri-phosphate molecules retrieval is better than di-phosphates which are retrieved better than mono-phosphates. These results suggest that our method may be useful to identify a diverse set of molecular shapes that could potentially bind to a given surface.

Prediction and Validation of a GDP Binding Site

The F₄₂₀ coenzyme plays important roles in archaea and eubacteria in a variety of biochemical reactions (e.g. methanogenesis, the formation of secondary metabolites, the degradation of nitroaromatic compounds, DNA repair)[56]. CofE, a F₄₂₀-0:gamma-glutamyl ligase, is responsible for the last two enzymatic steps in coenzyme F₄₂₀-2 biosynthesis. Belonging to a structurally uncharacterized family of enzymes, CofE from A. fulgidus was solved as a structural genomics target by the Midwest Center for Structural Genomics and found to be of novel fold (PDB:2g9i) (Figure 4a). Solvent accessible cavities were calculated using the CASTp webserver [27–29, 34], and the largest pocket, presumably the F420, GTP and L-glutamate binding pocket, was selected to query against the GPSS ligand surface library. The top-ranking surface was from GDP-binding protein from B. taurus (PDB:1tad, red, Figure 4c). A GDP molecule is posed into the surface based on the superposition from the alignment (Figure 11a, red GDP molecule). Based on this prediction, the protein was co-crystallized with GDP and a model of the complex was determined (Figure 11a, green GDP molecule)[56]. The GDP position had RMSD of 1.0Ǻ from the predicted pose (Figure 11b). The addition of the ligand also improved the resolution of the structure from 2.50Ǻ to 1.35Ǻ and allowed two loop regions to be modeled where no electron density was previously seen (Figure 11a, magenta).

Automated Protein Kinase Classification by ATP Binding Site Comparisons

Within the transferase family, the most well conserved ATP binding surfaces belong to protein kinases. Protein kinases play vital roles in regulating cellular pathways by phosphorylating other proteins. Malfunctioning kinases have been linked to a variety of diseases such as immunodeficiency, endocrine disorders and cancer, making them the target of drug discovery efforts. At their highest level, kinases are divided by the amino acid residues they target (serine/threonine or tyrosine) and further classified by more specific biochemical activity. Functional classification of kinase families has been undertaken by many methodologies utilizing primary sequence, structure, binding sites, pharmacophore profiles and expert manual analysis [57–60]

We apply our comparison methodology to protein kinases to assess our ability to automatically classify them into their functional subfamilies using only ATP binding surfaces. A dataset of 297 protein kinases with bound ligands, including both natural and synthetic molecules, were annotated using a combination of the PDB web query system, EC numbers, KinBase[61] and the Protein Kinase Resource[62]. For consistency, family nomenclature is applied from KinBase.

An all-against-all comparison was performed with results used to populate a distance matrix of SurfaceScreen scores. A dendrogram showing the complete-linkage clustering is shown in Figure 14 where each surface node is color coded by kinase subfamily. Overall, the method shares strong agreement with the annotated classification, as seen by the color banding. CDK2 kinases are the most ordered; with all members perfectly clustered together and distinct sub-grouping separating nucleotide ligands from small compound inhibitors. The CK2 and CAMP families also show divergence between natural ligands and inhibitors. The CAMP groupings are further clustered by the molecular weight of their bound ligands. The mitogen activated protein kinases (MAP) are successfully classified into their sub-families, but on distant nodes in the graph. In all families, we observe differentiation based on the activation state of the kinase.

Surface based classification conveys many similarities to other methods, but has advantages of additional functional insight that could not be automated in other methods. The ability to distinguish between different activation states and to organized surfaces based on ligand types and molecular weight differences could prove useful in developing binding profiles for enhanced specificity in kinase drug discovery.

Binding Surface Similarity of c-Abl Kinase Inhibitors

The surreptitious fusion of the cellular form of Abelson leukemia virus tyrosine kinase (c-Abl) with the breakpoint cluster region (BCR) gene disrupts the internal control mechanism causing increased tyrosine kinase activity[63]. The fusion protein BCR-Abl results in the disease chronic myelogenous leukemia (CML). Five structures of c-Abl proteins can be found in the PDB with two classes of small molecule inhibitors: pyrido [2,3.d]pyrimidine-type (PDB:1m52, 1opk) and 2-phenylaminopyrimidine-type (PDB:1iep, 1fpu, 1opj). Both inhibitor classes bind in the ATP binding, but 2-phenylaminopyrimidine-type bind exclusively in the inactive conformation of the activation loop. The smaller pyrido [2,3.d]pyrimidine-type class are indifferent to activation state, making them more potent but less specific inhibitors[64]. Our clustering accounts for this behavior and groups them in distinct nodes.

The 2-phenylaminopyrimidine-type inhibitor STI-571 (Figure 15a) is an effective inhibitor of c-Abl activity for treatment against CML [64–66]. It has been shown to be specific for tyrosine kinases and also inhibits stem-cell factor receptor kinase c-Kit (PDB:1t46). Results from querying c-Abl (PDB:1opj, Figure 15a) against the GPSS library show cKit is the best scoring non-ABL kinase. This cross reactivity is detected in our cluster (Figure 14b).

A surprising member of this cluster node is serine/threonine kinase p38 mitogen-activated protein (MAP) kianse (PDB:1kv2, Figure 14b). p38 MAP kinases play critical roles in regulation of proinflammatory cytokines such as tumor necrosis factor and interleukin-1 and are a target for many inflammatory diseases[67]. STI-571 is not currently known to inhibit, by design or mechanism, any serine/threonine kinase[68]. The structure of 1kv2, complexed with inhibitor B96, occupies a unique conformation for p38 MAPs, where the highly conserved DFG motif is turned out[67] (Figure 15b). This is the first observation of this state in a serine/threonine kinase, where it is a hallmark for tyrosine kinases (Figure 15a).

A superposition, based on the alignment of the binding surfaces of c-Abl and p38 MAP, shows that the inhibitors bind in similar orientation (Figure 15c). STI-571 can be posed into the p38 MAP surface, based on the alignment of the surfaces, with no steric clashing and preserves the orientation of several polar atoms (Figure 15d). While the conformation of this p38 MAP is unique and presumed to occur infrequently, its existence presents opportunity to explore the use of STI-571 and analogs for additional therapeutic uses. Automated surface classification can also provides important cross-reactivity analysis; where unexpected binding sites similarities could result in undesired side effects.