Search tips
Search criteria 


Logo of scirepAboutEditorial BoardFor AuthorsScientific Reports
Sci Rep. 2017; 7: 11652.
Published online 2017 September 14. doi:  10.1038/s41598-017-10412-z
PMCID: PMC5599562

Global organization of a binding site network gives insight into evolution and structure-function relationships of proteins


The global organization of protein binding sites is analyzed by constructing a weighted network of binding sites based on their structural similarities and detecting communities of structurally similar binding sites based on the minimum description length principle. The analysis reveals that there are two central binding site communities that play the roles of the network hubs of smaller peripheral communities. The sizes of communities follow a power-law distribution, which indicates that the binding sites included in larger communities may be older and have been evolutionary structural scaffolds of more recent ones. Structurally similar binding sites in the same community bind to diverse ligands promiscuously and they are also embedded in diverse domain structures. Understanding the general principles of binding site interplay will pave the way for improved drug design and protein design.


Ligand binding sites are responsible for various biological processes of proteins such as signal transduction and enzymatic activity. Understanding the characteristics of binding sites is essential in drug discovery and protein engineering. Still, there are many open questions that should be addressed to better understand binding sites. Are there common motifs of binding site structures? How much are binding sites similar or different quantitatively? How do binding sites evolve? Do they follow the same or different evolutionary pathways from protein domain structures? How do binding site similarities relate to known global structural, that is, fold similarities? Exploiting the vast structural data available in the Protein Data Bank (PDB)1, many binding site structure comparison algorithms and databases were developed210. However, previous studies were mainly focused on finding individual similar structures, and did not explore the global organization of binding site similarities4, 69. There have been studies to find clusters of known binding sites6, 7, 1116. However, most of them focused on certain protein families and did not address the global and evolutionary relationships between the representative binding site structures.

In contrast to binding site structures, protein domains have been investigated extensively since domains have long been considered as basic protein structural units that are stable, function and evolve17. There are many approaches to classify the protein domain universe systematically, such as SCOP18 and CATH19, that provide global hierarchical organizations of domain structures. To visualize the complex global relationships between domains the whole protein domain space has been projected onto two- and three-dimensional maps2023. Alternatively, the domain space has been represented as a network in which nodes are protein domains, and connections are drawn between domains that have similar folds24, 25. The three-dimensional projection of the domain space shows that all domain structures can be clustered into four clusters, which approximately correspond to the four SCOP classes, i.e., all alpha, all beta, alpha + beta, and alpha/beta classes2022. Another study proposes that the alpha/beta class forms a densely populated and functionally diverse core region of the protein domain universe23. A recent network-based approach shows that alpha/beta domains form a large connected domain network, whereas all-alpha, all-beta, and alpha + beta domains form smaller and disconnected networks24.

It is still an open question whether the domain is the optimal level to describe and classify protein structures26. Several studies suggest the existence of highly conserved and frequently occurring subdomain level motifs2730. Assuming that subdomain level local structures, such as binding sites, are the basic units of protein structure evolution, additional questions can be asked. What is the process by which the binding sites evolve and how does this process differ from the evolution of domain structures? Although convergent31 and divergent32 evolution of binding sites has been described, a mechanism that can encompass both evolutionary pathways has not been found yet. To help find answers to these questions, we analyzed global relationships between domains and binding sites.

In this study, we generated a global weighted network of all known binding sites based on their structural similarities calculated using ProBiS24, 10. Since the generated networks have large numbers of nodes and edges, they are incomprehensible at first glance. To reveal the hidden global organization of the networks, we reduced its complexity using a community detection approach that finds a subset of nodes that are more densely connected than the rest of the network. This approach, used in a protein-protein network improves the accuracy of protein function prediction3335.

We found that the size distribution of binding site communities follow a power-law distribution, implying that the largest communities may be the motifs of the most ancient binding sites and they many have served as structural scaffolds for other binding sites. We also found that functional diversity of proteins is independent of the binding site community size, indicating that binding sites are tightly coupled to protein function. We believe that the global binding site network communities can contribute to development of new approaches to function prediction, drug discoveries, and binding site design.

Results and Discussion

We generated a weighted binding site network from an all-to-all similarity comparison of binding site structures extracted from non-redundant sets of protein structures. In the binding site network, nodes are binding site structures and similar sites are connected by edges. Figure 1 presents the workflow in a sequence of 5 steps: 1) we generated the non-redundant protein structure databases using sequence identity cutoffs of 40, 70, and 90%; 2) non-redundant binding site structures were extracted from the databases; 3) structural similarities between all binding sites were calculated using the ProBiS2, 3; 4) the similarity scores were normalized into z-scores and the pairs of binding site pairs whose z-score is higher than a pre-defined threshold were connected by edges and 5) structurally similar binding sites were classified into communities to reduce the complexity of the network. Using the identified binding site communities, we investigated their size distribution, functional enrichment, and their relationships with ligands and domain structures.

Figure 1
Flow diagram for binding site community analysis.

Binding site network is continuous

To identify the continuity of the binding site network, we constructed networks with z-score thresholds from 1.5 to 4.0 with an interval of 0.5, calculated the fraction of binding sites included in the largest connected network and counted the number of separate networks (Fig. S1). We identified that most binding sites are connected as a single continuous network until z-score reaches 2.5. At this z-score, only the 3% highest similarities are considered and 53% of the nodes belong to the largest connected network. This suggests that binding site structure space is continuous until this level of statistical significance. When more stringent z-scores are used, the network becomes scattered into small and disconnected networks. When z-score is 3.0, the largest connected network includes only 9.5% of the nodes and all binding sites are divided into 1,176 separate networks. For further analysis, we used the largest connected network generated with a sequence identity cutoff of 70% and a z-score threshold of 1.5, which contains 15,789 binding sites (99.8% of all non-redundant binding sites) and 547,305 edges.

Community structure of the binding-site network reveals the global structure of binding sites

To reduce the complexity of the network, we detected the binding sites communities, i.e., groups of binding sites that are more closely related to each other than to the rest, using the Infomap algorithm36, 37. This procedure leads to the weighted network, in which nodes correspond to binding site communities and the widths of edges represent their structural similarity (Fig. 2). This community network reveals that larger communities are located at the core region of the network and smaller communities are at peripheral region. The two largest communities (C1.HEM.CLA and C2.CIT.AKG, which are most enriched with citric acid and heme molecules) are connected with most communities indicating that they play roles of network hubs. The largest communities show significant structural similarities between themselves. The 30 highest edge weights are observed between the 31 largest communities (Fig. S2). Only a few significant inter-community similarities are not connected with C1.HEM.CLA or C2.CIT.AKG.

Figure 2
Binding site community network. The 39 highest similarities between binding site communities and associated 20 binding site communities using ProBiS are displayed. A node corresponds to a binding site community and its size is proportional to the number ...

The network of binding site communities reveals the relationships between binding sites. It is noticeable that C2.CIT.AKG has a strong connection with C4.IPE.UDP and many connections with other communities, indicating that the structure of C2.CIT.AKG may be close to the structural scaffold of many binding sites. It is also identified that a group of four communities, which interact with ligands containing adenosine, C6.NAD.NAP, C9.FAD.NAD, C8.SAH.SAM, and C11.NAD.NAP, are closely inter-connected. This suggests that these communities may be structurally diversified due to functional reasons although they interact with similar ligands. Thus, comparing the difference between these communities may provide insight into how protein structures have evolved to make similar ligands interact with various proteins in a distinctive way.

To test the robustness of our findings, we generated another binding site network using a different binding site comparison program G-LoSA38 and performed the community detection on the network (Fig. S3). The largest binding site community is most enriched with heme molecules. It is also identical that the five largest communities are also the most strongly connected. To verify that the similarity between community structures obtained with ProBiS and G-LoSA are statistically significant, we calculated the normalized mutual information (NMI) between two community structures and compared it with the values obtained with randomly permutated communities. The NMI value of the two original community structures was 0.337 while that of the random permutations was 0.1226 ± 0.0012, resulting in a P-value < 0.0001, which demonstrates that the community structures are independent of the choice of a binding site comparison algorithm. We also iterated the identical analysis with networks constructed with a different z-value (z = 2.0) (Fig. S4) and sequence identity thresholds, 40% and 90% (Fig. S5), and obtained similar community structures to those shown in Fig. 2.

Our analysis is the first study that identified the similarities between the communities of binding sites compared to previous PDB-wide binding site analyses13, 14. The most recent comprehensive global comparison of binding sites was reported by Gao and Skolnick13. They performed a PDB-wide clustering of about 20,000 non-redundant binding site structures using the APoc method8, and suggested that all binding site structures may be categorized into around 1,000 shapes. The APoc method normalizes a similarity score based on the size of a query binding-site structure, which makes the score asymmetric and dependent on a query-template definition. Thus, only the pairs that show significant similarities in both directions are used, i.e., only binding sites with similar sizes are considered, leading to tightly inter-connected communities. Similarly, Kinjo and Nakamura also performed a global similarity analysis of binding sites14. By performing all-to-all similarity comparison of 180,000 binding sites known by June 2008, they constructed 11,532 separate networks consisting of highly similar binding sites both functionally and structurally. The large number of separated networks is due to the use of an exceedingly stringent similarity criterion to connect binding sites (P-value < 10−15). In contrast to previous studies, we used the z-score of the ProBiS scoring function4, which normalizes a similarity score based on the number of aligned residues, resulting in scores independent of the binding site size. Thus, the binding site communities detected in this study consist of more remotely and weakly related binding sites. In addition, our network reveals relative similarities between binding site communities, which may shed light on the evolution of binding site structures.

Sizes of binding site communities follow a power-law distribution

The size distribution of binding site communities follows a power-law distribution. The linearity of the complementary cumulative distribution function39 clearly shows that the community sizes are distributed by a power-law: f(k) = 15953k −2.4, where f(k) is the number of communities with a size of k (Fig. 3A and B). The power law distribution of the communities indicates that a few extremely large communities coexist with many small ones. The 30 communities shown in Fig. 2 include 47.7% of all binding sites (Fig. 3C), thus, a few largest communities and their direct neighbors represent the majority of them. The generative models explaining the power-law distribution are generally based on the constant birth rate assumption, i.e., the probabilities to be duplicated or emerge a new community are identical for all binding sites regardless of their structure or function4043. The only dominant factor that determines the future community size is its current size. Thus, the power-law distribution suggests that binding sites may have evolved based on a simple and universal mechanism, not constrained by a structural or functional necessity44, 45.

Figure 3
Size distributions of binding site communities. (A) The frequency of binding site communities of size k, (B) the complementary cumulative distribution function (cdf) of community sizes P(k), and (C) the cumulative fraction of binding sites included in ...

The power-law distribution of the binding site communities also implies that the largest communities may be the most ancient binding site structures and they might have been structural scaffolds for the evolution of other binding sites. A power-law distribution is widely found in genomics44, 46 and it has been successfully explained by simple generative models41, 44, which indicate that a larger community is generally older than a smaller one. The analyses of protein domain structures show that alpha/beta domains are located at denser regions, i.e., have more structural neighbors, and, generally, are older than other folds2224. To further test this hypothesis, we calculated the enrichment of taxonomies of binding sites47. If the hypothesis is valid, binding sites from older organisms should be mostly found in large communities. Our analysis indicates this is the case. Binding site structures from cyanobacteria are enriched in only the first two largest communities (Table S5). C1.HEM.CLA is enriched with a DNA binding and cytochrome-c oxidase activity and C2.CIT.AKG is highly enriched with transition metal ion (Zn2+, Fe2+, Fe3+, Mn2+ and Co2+) binding functions. This is consistent with a domain-based proteome study48 and phylogenetic analyses49, 50. Whereas binding sites from human proteins are enriched not only in large communities (C3.GDP.ADP and C7.ANP.ATP) but also in many small communities (C13, C51, C52, C90, C100, C183, C189, and C241) (Table S6). Interestingly, C3.GDP.ADP is highly enriched with ATPase activity and an ABC transporter function. This fact is consistent with the recent domain structure analysis, which found that the ABC transporters might be the oldest aerobic metabolic enzymes51. Thus, the taxonomy analysis also supports the hypothesis that larger communities may be more ancient and have been evolutionary structural scaffolds of smaller ones.

The binding sites that belong to small communities with less than 15 binding sites, which do not follow the power-law distribution mainly bind to glycans or buffering agents. The most frequently found ligand in these small communities is N-acetyl-D-glucosamine (NAG), which is one of the most common building blocks of glycans. Because N-glycosylation process involves the formation of a covalent bond between a glycan and the residue at a glycosylation site, NAG binding site residues may play little role in its binding, which makes these residues non-specific binding partners. Following NAG, 2-methyl-2,4-pentanediol (MPD) and 2-(N-morpholino) ethanesulfonic acid (MES) are the next most frequently found ligands in the small communities; these are used as buffering or precipitation agents. This suggests that the binding sites in the small communities may not be true binding sites and the power-law fitting may allow us to discriminate between specific and non-specific interactions.

Binding site communities are functionally specific

To study the relationships between binding site structures and functions and to functionally characterize each community, we performed the functional enrichment analysis by mapping the gene ontology (GO) annotations52 onto communities. We found that the largest binding site communities are associated with different functions and have few overlapped functions despite their significant structural similarities. The enriched functions of C1.HEM.CLA detected with a P-value threshold of 0.001 are compared with those of the other top 30 largest communities. Among them, only 9 communities have at least one common functional annotation. Further, the all-to-all comparison between the top 30 communities showed that the average number of shared functional annotation is only 0.47. These results indicate that a binding-site structure is strongly coupled with the function of a protein and our binding-site centric classification of proteins is thus very function specific.

For example, the two tightly connected communities C6.NAD.NAP and C9.FAD.NAD have NAD as their major ligands and are enriched with oxidoreductase activity related terms. However, C6.NAD.NAP is enriched with NADH-specific enoyl-ACP reductase activity (GO:0004318) while C9 is enriched with dihydrolipoyl dehydrogenase activity (GO:0004791). Overall, the two communities share only 4 enriched GO terms although C6 and C9 have 24 and 18 enriched terms, respectively. This indicates that our community analysis also differentiates the binding sites that bind to the same ligand but have different functions.

Next, we investigated the functional diversity of proteins associated with binding site communities. It is known that the domain structures located at the core region of the protein structures map tend to have more diverse functions than the proteins at peripheral regions23. We measured the functional diversity by obtaining the average number of distinct GO molecular function (N MF) and biological process (N BP) terms associated with each protein included in each community23. Our results show a different pattern from the results of domain structure analysis; there is little correlation between the centrality of a binding site community and the functional diversity of its members (Fig. 4A). The functional diversity of proteins is almost uniform regardless of its community size, and is 4.9, which is the average functional diversity of all proteins in the network.

Figure 4
Shannon information (entropy) values of the ligand/domain compositions and the functional diversity of binding site communities The x-axes represent the community size using a log-scale. The y-axis of (A) represents the functional diversity of the communities. ...

Binding site communities are promiscuous

Do similar binding sites bind to similar ligands? To quantify the ligand specificity of binding site communities, we calculated the ligand composition Shannon information S (entropy) of a binding site community. The entropy is calculated using the following formula:

An external file that holds a picture, illustration, etc.
Object name is 41598_2017_10412_Article_IEq5.gif
, where i is a ligand or a domain index and p i is the fraction of ligand in a community. A large value of S indicates that a binding site in the community binds to diverse ligands, and S close to zero, indicates that it binds to a few specific ligands. We found an almost linear relationship between the size and the ligand composition Shannon information of the binding site communities (Fig. 4B). The linear relationship is close to the theoretical maximum Shannon information when all ligands are different, i.e., S = log N comm, where N comm is the size of a community. We also calculated the variance of distances between ligands in a community to identify whether the same relationship is valid when the chemical similarities of ligands are considered (Fig. 4C). The variances of ligand similarities of binding site communities are compared with those of the same number of randomly selected ligands. A lower variance indicates that the included ligands are more similar to each other. The analysis shows that the majority of detected communities have large variances approaching random distribution. On the other hand, only a few communities have significantly lower variances than randomly selected sets, e.g., eleven communities have variances lower than 0.1. These results indicate that similar binding sites included in the same community bind to diverse ligands and the diversity of ligands increases with the size of communities, which is consistent with previous studies13, 54.

This high ligand promiscuity of binding site communities may be explained by the constructive neutral evolution scenario44, 45, 55, in which a system evolves via the accumulation of irreversible dependencies between related parts of the system, not adaptation. The scenario is known to be consistent with the power-law distribution of community44, 45. If binding sites in larger communities are more ancient than ones in smaller communities, more diverse ligands could have been tested against them, which may have resulted in more binding partners than recently evolved sites.

Similar binding sites are found in diverse domain structures

Are similar binding sites embedded in similar domain structures? To answer this question, we calculated domain composition Shannon information S values by investigating to which CATH domain a binding site belongs. Overall, the domain composition S domain of a binding site community increases almost linearly with the logarithm of its size, similar to the ligand composition S ligand (Fig. 4D). In other words, a binding site structure from a larger community is associated with more diverse backbone structures. However, S domain values deviate more from the theoretical maximum compared to S ligand values, which indicates that a binding site structure depends weakly on its domain structure.

The fact that binding site structures depend only weakly on their corresponding domain structures indicates that the protein structure segments may have evolved independently. In other words, the unit of protein structure evolution could be smaller than a domain2729, 56. If a binding site and its backbone structure were strongly coupled, S domain should be significantly lower than the theoretical maximum and almost completely independent of the community size. Halabi et al.29 suggest the concept of protein sectors, in which a whole domain can be divided into the subgroups of residues that are structurally adjacent and have distinct functional roles. They show that the S1A serine protease domain consists of three functionally independent groups of connected residues that have the roles of ligand specificity, protein stability and catalytic core, respectively29. The protein sector analyses of the PDZ, PAS, SH2, and SH3 domain families show that the ligand binding sites of these domains are detected as independent sectors29.

The protein sector hypothesis could explain the differences between the functional diversity distributions of the domain network and our binding site network. It was identified that the core of the domain network, mostly alpha/beta classes of SCOP, has high functional diversity, and the peripheral region has low functional density23, 24. However, the functional diversity of binding sites is almost uniform across the binding site network, which suggests that binding sites are more important in determining the function of a protein than its domain structure. If a domain structure consists of independent sectors, the combinations of functionally distinct sectors can achieve high functional diversity of a domain. A binding site that recognizes a specific ligand can form a domain with various sectors performing distinct enzymatic activities, and the sectors in the domain communicate via an allosteric mechanism. In contrast, the domains with lower functional diversity may consist of sectors having similar functions.


It is reasonable to propose that the global binding site network and detecting its community structure can contribute to development of new methods for structure-based function inference, drug discovery, and binding site design. High functional specificity of binding site communities proposes that more accurate protein function predictions will be possible using the knowledge on similar binding sites and the enriched functions of the communities than by using protein domain structures34, 35. The enhanced prediction accuracy may be more pronounced when a target protein adopts an alpha/beta structure, since this fold has high functional promiscuity23, which can lead to many false positive predictions. Also, knowing which binding residues are conserved in a community may significantly improve the success rate of drug discovery and protein design.

Our analysis identified that interactions between ligands and binding pockets are generally promiscuous. This explains the “off-target” effects of drugs because one specific ligand may bind to different binding pockets13, 57. This also suggests that a conventional strategy to find new drug candidates that considers the nearest neighbors of a specific ligand may be inefficient. Although binding site-ligand interactions are generally promiscuous, our analysis also indicates that different binding site communities are functionally distinct. Based on this observation, one alternative approach to overcome the limitation of promiscuous interaction is to extract useful information from the ensemble of related binding sites and their ligands. For example, certain interactions between a given ligand and related binding sites in a community may be conserved, while no such interactions exist between the same ligand and binding sites in another community. If this is the case, such difference may be used to screen or design new more selective drug candidates. Similar approaches, using an ensemble of related ligand-binding site interactions or considering remote homologs of a target binding site, have been suggested previously5863. We believe that the network of binding site communities will provide a basis for exploring more efficient computational approaches for drug discovery and design.

Materials and Methods

ProBiS algorithm

To define surface residues, solvent accessible area was defined by rolling a spherical probe of a radius of 1.4 Å on the van der Waals surface of protein atoms. Residues located up to 4 Å below this surface were considered as surface residues. Binding site residues were determined as those surface residues separated by <5 Å from any ligand atom. Ligands were defined as those HET codes in the PDB file having >7 heavy atoms, and not being a modified residue denoted by a MODRES code. Thus, surface binding sites as well as binding sites deep inside the proteins that are connected with the exterior by a channel with radius of at least 1.4 Å were considered. Ligand binding sites that are completely buried within a single protein chain are very rare in our experience and were not considered in this study; nevertheless, binding sites buried by two or more chains were still considered, as each protein chain was considered separately. Metal ions are not considered as ligands in this study. Nevertheless, metal ion binding sites are considered if they are in vicinity (<5 Å) of a small molecule ligand.

A detected binding site surface patch was represented as a graph consisting of a set of vertices and edges connecting them. Vertices in a graph represent the functional groups of surface amino acid residues. A functional group is characterized with 5 groups based on its physicochemical property: hydrogen bond donor, acceptor, mixed donor/acceptor, aromatic, and aliphatic64. If two vertices are separated by <15 Å, they are connected by an edge.

For a given pair of binding site graphs, their product graph was constructed. For two graphs, G 1 and G 2, the vertex set of the Cartesian product graph is defined as: HV(G1) × V(G2) = {(uv)|u ∈ V(G1) and v ∈ (G2)}. If the physicochemical properties of two vertices of the initial graphs, u and v, are different, the corresponding pair is not considered to generate the product graph. Two vertices of a product graph, (u 1,v 1) and (u 2,v 2), are connected by an edge if and only if the distances between u 1 and v 1 and between u 2 and v 2 differ by <2 Å. The generated product graph is an approximate representation of all possible superposition of two structures. From the generated product graph, we detected its maximum clique2, 65, i.e., the largest complete sub-graph of a graph where all vertices are connected to each other.

For each pair of aligned binding sites, we calculated its similarity score using structural similarity and evolutionary similarity4. For each pair of binding site alignment, four criteria were used to calculate the local alignment score: (1) surface angle, (2) surface patch RMSD, (3) surface patch size, and (4) an E-value calculated with the Karlin-Altschul equation66. For each surface patch, a surface vector originating from the geometric center of the patch and pointing to the perpendicular direction of the surface was generated. If the angle between a pair of surface vectors was larger than 90° or the number of aligned vertices were smaller than 10, the pair was discarded. For the remaining pairs, the alignment scores, al score, were calculated as follows:

An external file that holds a picture, illustration, etc.
Object name is 41598_2017_10412_Article_IEq7.gif
, where RMSD is the surface patch RMSD between pairs of superimposed vertices, n vert is the number of aligned vertices, and e value is the alignment expectation value4. The raw alignment scores were normalized to z-scores as follows:
An external file that holds a picture, illustration, etc.
Object name is 41598_2017_10412_Article_IEq8.gif
, where μ and σ are the average and the standard deviation of the alignment scores of all pairs, and the values of μ and σ are 2.0 and 2.2, respectively.

Constructing a binding site network and detecting communities

We constructed an initial non-redundant protein single chain database from a Nov 2013 PDB release1 using a sequence identify cutoff of 95% resulting in 42,282 unique protein chains. Based on the non-redundant protein chain database, an all-to-all comparison between the surfaces of the non-redundant chains was performed using ProBiS24. All details on the ProBiS algorithm are reported in SI Methods. From the statistics of similarity scores from the all-to-all comparison, we calculated the z-scores of all similarity scores. A pair of binding sites is connected by an edge if the z-score of the pair is higher than a threshold value and the z-score of the edge is represented by its weight. If the z-score of a pair is lower than a threshold, the pair is considered to be disconnected. The G-LoSA method was used on the same binding site set. A community detection analysis was performed using the Infomap method based on the minimum description length principle36, 37. This method detects the community structure of a network by minimizing the amount of information to describe the itinerary of a random walker on the network. One of the advantages of the Infomap method is that it yields how strongly communities are related, which is interpreted as the similarity between a pair of binding site communities in this study. The weights of inter-community edges are proportional to the probability of a random walker selecting the path. A higher inter-community probability indicates higher structural similarity between binding site communities.

Power-law fitting

We fitted the size distribution of binding site communities to a power-law distribution using the powerlaw Python package67. The optimal minimum threshold of the cluster size for fitting was determined by the algorithm to yield the minimum Kolmogorov-Smirnov distance between the data and the fit68. The sizes of communities and their frequencies are fitted to a power-law distribution, p(k) = Pr(Kk) =  Ckα, where Pr(K) is the probability and p(k) probability density to form a community with k binding sites; α is a scaling parameter and C is a normalization constant. The fitting showed that the communities with more than 14 binding sites, k min = 15, follow a power-law distribution with α = 2.4067, 68. The total number of communities with more than 14 sites is 236. The most common way to identify a power-law distribution is to check the linearity of the cumulative distribution (cdf), P(k) = Pr(K > k), on a log-log plot68.

Ligand and domain composition Shannon information calculation

The ligand and domain composition Shannon information values were calculated using the following formula: S = −∑ipi ln pi, where i is a ligand or a domain index and p i is the fraction of ligand or domain i relative to the total number of binding sites in a community. For ligand composition Shannon information calculation, we used the ligand IDs defined in the PDB database. Ligands are considered as different if they have different ligand IDs in the PDB. To identify the promiscuity of a binding site community, we also calculated the variance of ligand similarities. The variance of ligands in a community C is calculated with the following equation:

An external file that holds a picture, illustration, etc.
Object name is 41598_2017_10412_Article_IEq11.gif
, where T ij is the Tanimoto coefficient53 between ligands i and j.

For domain composition Shannon information calculation, the chain IDs of residues participating in forming a binding site were identified and their domain IDs defined by the CATH database19, 69 were used. If a binding site consists of residues from multiple distinct chains, the binding site was considered to be located at the interface of those chains. A binding site located at the interface of a homo-oligomer was considered to be different from that of the corresponding monomer.

Functional enrichment analysis

To identify which functional annotations are significantly enriched in a community, the P-values of GO terms52 associated with proteins included in a community were calculated. In this study, we used the MF and BP terms of GO annotations. P-values were calculated using a hypergeometric distribution:

An external file that holds a picture, illustration, etc.
Object name is 41598_2017_10412_Article_Equ1.gif

where n is the total number of proteins in the network with known GO terms, m is the size of a community, f is the total number of proteins in the network associated with a function of interest, and k is the number of proteins in the community with the function of interest.

In other words, the proteins in the network with known associated GO terms were used as the background protein list. The GO annotations of the proteins were adopted from the SIFTS database47. Since GO terms are hierarchical, only the most specific terms of each protein were used23. For each community, P-values were adjusted for multiple comparisons using the Benjamini-Hochberg procedure with a false discovery rate level of 0.0170, indicating that less than 1% of identified enriched functional annotations are expected to be false positives.

Data availability

All relevant data are available from the authors upon request.

Electronic supplementary material

Dataset 1(3.3M, zip)

Dataset 2(519K, xlsx)

Dataset 3(39K, xlsx)

Dataset 4(60K, xlsx)

Dataset 5(49K, xlsx)

Dataset 6(14K, xlsx)


The authors thank Richard Pastor and Mark Knepper for helpful discussions. J.L. and B.R.B. are supported by the Intramural Research Program of the National Institutes of Health, National Heart, Lung and Blood Institute (NHLBI) under Project No. Z01 HL001050. J.K. and D.J. are supported by Grant Nos. J1-6743, L7-8269 and P1-0002 of the Slovenian Research Agency. J.K. acknowledges the Fulbright Scholarship at the NHLBI, National Institutes of Health. This work utilized the high-performance computational capabilities at the National Institutes of Health, Bethesda, MD (NHLBI LoBoS cluster).

Author Contributions

Author Contributions

J.L. conceived and designed the experiments. J.L. and J.K. performed the experiments and analyzed the data. J.L., J.K., D.J., and B.R.B. wrote the paper.


Competing Interests

The authors declare that they have no competing interests.


Juyong Lee and Janez Konc contributed equally to this work.

Electronic supplementary material

Supplementary information accompanies this paper at doi:10.1038/s41598-017-10412-z

Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


1. Berman HM, et al. The Protein Data Bank. Nucleic Acids Res. 2000;28:235–242. doi: 10.1093/nar/28.1.235. [PMC free article] [PubMed] [Cross Ref]
2. Konc J, Janežič D. ProBiS algorithm for detection of structurally similar protein binding sites by local structural alignment. Bioinformatics. 2010;26:1160–1168. doi: 10.1093/bioinformatics/btq100. [PMC free article] [PubMed] [Cross Ref]
3. Konc J, Depolli M, Trobec R, Rozman K, Janežič D. Parallel-ProBiS: Fast parallel algorithm for local structural comparison of protein structures and binding sites. J. Comput. Chem. 2012;33:2199–2203. doi: 10.1002/jcc.23048. [PubMed] [Cross Ref]
4. Konc J, Česnik T, Konc JT, Penca M, Janežič D. ProBiS-database: Precalculated binding site similarities and local pairwise alignments of PDB structures. J. Chem. Inf. Model. 2012;52:604–612. doi: 10.1021/ci2005687. [PMC free article] [PubMed] [Cross Ref]
5. Konc J, Janežič D. ProBiS-ligands: A web server for prediction of ligands by examination of protein binding sites. Nucleic Acids Res. 2014;42:215–220. doi: 10.1093/nar/gku460. [PMC free article] [PubMed] [Cross Ref]
6. Kufareva I, Ilatovskiy AV, Abagyan R. Pocketome: An encyclopedia of small-molecule binding sites in 4D. Nucleic Acids Res. 2012;40:535–540. doi: 10.1093/nar/gkr825. [PMC free article] [PubMed] [Cross Ref]
7. Ito JI, Tabei Y, Shimizu K, Tomii K, Tsuda K. PDB-scale analysis of known and putative ligand-binding sites with structural sketches. Proteins Struct. Funct. Bioinforma. 2012;80:747–763. doi: 10.1002/prot.23232. [PubMed] [Cross Ref]
8. Gao M, Skolnick J. APoc: Large-scale identification of similar protein pockets. Bioinformatics. 2013;29:597–604. doi: 10.1093/bioinformatics/btt024. [PMC free article] [PubMed] [Cross Ref]
9. Nisius B, Sha F, Gohlke H. Structure-based computational analysis of protein binding sites for function and druggability prediction. J. Biotechnol. 2012;159:123–134. doi: 10.1016/j.jbiotec.2011.12.005. [PubMed] [Cross Ref]
10. Konc J, et al. ProBiS-CHARMMing: Web Interface for Prediction and Optimization of Ligands in Protein Binding Sites. J. Chem. Inf. Model. 2015;55:2308–2314. doi: 10.1021/acs.jcim.5b00534. [PubMed] [Cross Ref]
11. Xie L, Bourne PE. Detecting evolutionary relationships across existing fold space, using sequence order-independent profile-profile alignments. Proc. Natl. Acad. Sci. USA. 2008;105:5441–5446. doi: 10.1073/pnas.0704422105. [PubMed] [Cross Ref]
12. Xie L, Xie L, Bourne PE. A unified statistical model to support local sequence order independent similarity searching for ligand-binding sites and its application to genome-based drug discovery. Bioinformatics. 2009;25:305–312. doi: 10.1093/bioinformatics/btp220. [PMC free article] [PubMed] [Cross Ref]
13. Gao M, Skolnick J. A Comprehensive Survey of Small-Molecule Binding Pockets in Proteins. PLoS Comput. Biol. 2013;9:e1003302. doi: 10.1371/journal.pcbi.1003302. [PMC free article] [PubMed] [Cross Ref]
14. Kinjo AR, Nakamura H. Comprehensive Structural Classification of Ligand-Binding Motifs in Proteins. Structure. 2009;17:234–246. doi: 10.1016/j.str.2008.11.009. [PubMed] [Cross Ref]
15. Zhang Z, Grigorov MG. Similarity networks of protein binding sites. Proteins Struct. Funct. Genet. 2006;62:470–478. doi: 10.1002/prot.20752. [PubMed] [Cross Ref]
16. Park K, Kim D. Binding similarity network of ligand. Proteins. 2008;71:960–71. doi: 10.1002/prot.21780. [PubMed] [Cross Ref]
17. Chothia C. Proteins. One thousand families for the molecular biologist. Nature. 1992;357:543–544. [PubMed]
18. Andreeva A, et al. SCOP database in 2004: refinements integrate structure and sequence family data. Nucleic Acids Res. 2004;32:D226–229. doi: 10.1093/nar/gkh039. [PMC free article] [PubMed] [Cross Ref]
19. Sillitoe I, et al. CATH: comprehensive structural and functional annotations for genome sequences. Nucleic Acids Res. 2015;43:D376–D381. doi: 10.1093/nar/gku947. [PMC free article] [PubMed] [Cross Ref]
20. Hou J, Sims GE, Zhang C, Kim S-H. A global representation of the protein fold space. Proc. Natl. Acad. Sci. USA. 2003;100:2386–2390. doi: 10.1073/pnas.2628030100. [PubMed] [Cross Ref]
21. Hou J, Jun S-R, Zhang C, Kim S-H. Global mapping of the protein structure space and application in structure-based inference of protein function. Proc. Natl. Acad. Sci. USA. 2005;102:3651–3656. doi: 10.1073/pnas.0409772102. [PubMed] [Cross Ref]
22. Choi I-G, Kim S-H. Evolution of protein structural classes and protein sequence families. Proc. Natl. Acad. Sci. USA. 2006;103:14056–14061. doi: 10.1073/pnas.0606239103. [PubMed] [Cross Ref]
23. Osadchy M, Kolodny R. Maps of protein structure space reveal a fundamental relationship between protein structure and function. Proc. Natl. Acad. Sci. USA. 2011;108:12301–12306. doi: 10.1073/pnas.1102727108. [PubMed] [Cross Ref]
24. Nepomnyachiy S, Ben-Tal N, Kolodny R. Global view of the protein universe. Proc. Natl. Acad. Sci. USA. 2014;111:11691–11696. doi: 10.1073/pnas.1403395111. [PubMed] [Cross Ref]
25. Pascual-García A, Abia D, Ortiz ÁR, Bastolla U. Cross-over between discrete and continuous protein structure space: Insights into automatic classification and networks of protein structures. PLoS Comput. Biol. 2009;5:e1000331. doi: 10.1371/journal.pcbi.1000331. [PMC free article] [PubMed] [Cross Ref]
26. Valas RE, Yang S, Bourne PE. Nothing about protein structure classification makes sense except in the light of evolution. Curr. Opin. Struct. Biol. Biol. 2009;19:329–334. doi: 10.1016/ [PMC free article] [PubMed] [Cross Ref]
27. Szustakowski JD, Kasif S, Weng Z. Less is more: Towards an optimal universal description of protein folds. Bioinformatics. 2005;21:66–71. doi: 10.1093/bioinformatics/bti1111. [PubMed] [Cross Ref]
28. Friedberg I, Godzik A. Connecting the protein structure universe by using sparse recurring fragments. Structure. 2005;13:1213–1224. doi: 10.1016/j.str.2005.05.009. [PubMed] [Cross Ref]
29. Halabi N, Rivoire O, Leibler S, Ranganathan R. Protein sectors: evolutionary units of three-dimensional structure. Cell. 2009;138:774–786. doi: 10.1016/j.cell.2009.07.038. [PMC free article] [PubMed] [Cross Ref]
30. Lupas AN, Ponting CP, Russell RB. On the evolution of protein folds: are similar motifs in different protein folds the result of convergence, insertion, or relics of an ancient peptide world? J. Struct. Biol. 2001;134:191–203. doi: 10.1006/jsbi.2001.4393. [PubMed] [Cross Ref]
31. Gherardini PF, Wass MN, Helmer-Citterich M, Sternberg MJE. Convergent Evolution of Enzyme Active Sites Is not a Rare Phenomenon. J. Mol. Biol. 2007;372:817–845. doi: 10.1016/j.jmb.2007.06.017. [PubMed] [Cross Ref]
32. Horvath, M. M., Wang, X., Resnick, M. a. & Bell, D. a. Divergent evolution of human p53 binding sites: Cell cycle versus apoptosis. PLoS Genet. 3, 1284–1295 (2007). [PMC free article] [PubMed]
33. Lee J, Gross SP, Lee J. Modularity optimization by conformational space annealing. Phys. Rev. E - Stat. Nonlinear, Soft Matter Phys. 2012;85:56702. doi: 10.1103/PhysRevE.85.056702. [PubMed] [Cross Ref]
34. Lee J, Lee J. Hidden information revealed by optimal community structure from a protein-complex bipartite network improves protein function prediction. PLoS One. 2013;8:e60372. doi: 10.1371/journal.pone.0060372. [PMC free article] [PubMed] [Cross Ref]
35. Lee J, Gross SP, Lee J. Improved network community structure improves function prediction. Sci. Rep. 2013;3:2197. doi: 10.1038/srep02197. [PMC free article] [PubMed] [Cross Ref]
36. Rosvall M, Bergstrom CT. An information-theoretic framework for resolving community structure in complex networks. Proc. Natl. Acad. Sci. USA. 2007;104:7327–31. doi: 10.1073/pnas.0611034104. [PubMed] [Cross Ref]
37. Rosvall M, Bergstrom CT. Maps of random walks on complex networks reveal community structure. Proc. Natl. Acad. Sci. USA. 2008;105:1118–1123. doi: 10.1073/pnas.0706851105. [PubMed] [Cross Ref]
38. Lee HS, Im W. Identification of Ligand Templates using Local Structure Alignment for Structure-Based Drug Design. J. Chem. Inf. Model. 2012;52:2784–2795. doi: 10.1021/ci300178e. [PMC free article] [PubMed] [Cross Ref]
39. Newman, M. E. J. Networks. An introduction. Oxford University Press (2010).
40. Newman M. Power laws, Pareto distributions and Zipf’s law. Contemp. Phys. 2005;45:323–351. doi: 10.1080/00107510500052444. [Cross Ref]
41. Karev GP, Wolf YI, Rzhetsky AY, Berezovskaya FS, Koonin EV. Birth and death of protein domains: a simple model of evolution explains power law behavior. BMC Evol. Biol. 2002;2:18. doi: 10.1186/1471-2148-2-18. [PMC free article] [PubMed] [Cross Ref]
42. Koonin EV, Wolf YI, Karev GP. The structure of the protein universe and genome evolution. Nature. 2002;420:218–223. doi: 10.1038/nature01256. [PubMed] [Cross Ref]
43. Barabasi A-L, Albert R. Emergence of scaling in random networks. Science (80-.). 1999;286:509–512. doi: 10.1126/science.286.5439.509. [PubMed] [Cross Ref]
44. Koonin EV. Are there laws of genome evolution? PLoS Comput. Biol. 2011;7:e1002173. doi: 10.1371/journal.pcbi.1002173. [PMC free article] [PubMed] [Cross Ref]
45. Gray MW, Lukes J, Archibald JM, Keeling PJ, Doolittle WF. Cell biology. Irremediable complexity? Science. 2010;330:920–921. [PubMed]
46. Luscombe NM, Qian J, Zhang Z, Johnson T, Gerstein M. The dominance of the population by a selected few: power-law behaviour applies to a wide variety of genomic properties. Genome Biol. 2002;3:research0040.1–0040.7. doi: 10.1186/gb-2002-3-8-research0040. [PMC free article] [PubMed] [Cross Ref]
47. Velankar S, et al. SIFTS: Structure Integration with Function, Taxonomy and Sequences resource. Nucleic Acids Res. 2013;41:483–489. doi: 10.1093/nar/gks1258. [PMC free article] [PubMed] [Cross Ref]
48. Dupont CL, Yang S, Palenik B, Bourne PE. Modern proteomes contain putative imprints of ancient shifts in trace metal geochemistry. Proc. Natl. Acad. Sci. USA. 2006;103:17822–17827. doi: 10.1073/pnas.0605798103. [PubMed] [Cross Ref]
49. Dupont CL, Butcher A, Valas RE, Bourne PE, Caetano-Anollés G. History of biological metal utilization inferred through phylogenomic analysis of protein structures. Proc. Natl. Acad. Sci. USA. 2010;107:10567–10572. doi: 10.1073/pnas.0912491107. [PubMed] [Cross Ref]
50. David La, Alm EJ. Rapid evolutionary innovation during an Archaean genetic expansion. Nature. 2011;469:93–96. doi: 10.1038/nature09649. [PubMed] [Cross Ref]
51. Kim KM, et al. Protein domain structure uncovers the origin of aerobic metabolism and the rise of planetary oxygen. Structure. 2012;20:67–76. doi: 10.1016/j.str.2011.11.003. [PubMed] [Cross Ref]
52. Gene T, Consortium O, Gene T, Go O. Gene Ontology Consortium: going forward. Nucleic Acids Res. 2014;43:D1049–D1056. [PMC free article] [PubMed]
53. Willett P, Barnard JM, Downs GM. Chemical Similarity Searching. J. Chem. Inf. Comput. Sci. 1998;38:983–996. doi: 10.1021/ci9800211. [Cross Ref]
54. Ma B, Shatsky M, Wolfson HJ, Nussinov R. Multiple diverse ligands binding at a single protein site: a matter of pre-existing populations. Protein Sci. 2002;11:184–197. doi: 10.1110/ps.21302. [PubMed] [Cross Ref]
55. Stoltzfus A. On the possibility of constructive neutral evolution. J. Mol. Evol. 1999;49:169–181. doi: 10.1007/PL00006540. [PubMed] [Cross Ref]
56. Friedberg I, Godzik A. Fragnostic: Walking through protein structure space. Nucleic Acids Res. 2005;33:249–251. doi: 10.1093/nar/gki363. [PMC free article] [PubMed] [Cross Ref]
57. Zhou H, Gao M, Skolnick J. Comprehensive prediction of drug-protein interactions and side effects for the human proteome. Sci. Rep. 2015;5:11090. doi: 10.1038/srep11090. [PMC free article] [PubMed] [Cross Ref]
58. Skolnick J, Gao M, Roy A, Srinivasan B, Zhou H. Implications of the small number of distinct ligand binding pockets in proteins for drug discovery, evolution and biochemical function. Bioorg. Med. Chem. Lett. 2015;25:1163–1170. doi: 10.1016/j.bmcl.2015.01.059. [PMC free article] [PubMed] [Cross Ref]
59. Jian J-W, et al. Predicting Ligand Binding Sites on Protein Surfaces by 3-Dimensional Probability Density Distributions of Interacting Atoms. PLoS One. 2016;11:e0160315. doi: 10.1371/journal.pone.0160315. [PMC free article] [PubMed] [Cross Ref]
60. Salentin S, Haupt VJ, Daminelli S, Schroeder M. Polypharmacology rescored: Protein-ligand interaction profiles for remote binding site similarity assessment. Prog. Biophys. Mol. Biol. 2014;116:174–186. doi: 10.1016/j.pbiomolbio.2014.05.006. [PubMed] [Cross Ref]
61. Marsh, L. Strong Ligand-Protein Interactions Derived from Diffuse Ligand Interactions with Loose Binding Sites. Biomed Res. Int. 2015, (2015). [PMC free article] [PubMed]
62. Tan, Z., Chaudhai, R. & Zhang, S. Polypharmacology in Drug Development: A Minireview of Current Technologies. ChemMedChem 1211–1218 doi:10.1002/cmdc.201600067 (2016). [PubMed]
63. Duran-frigola M, et al. Detecting similar binding pockets to enable systems polypharmacology. PLoS Comput. Biol. 2017;13:e1005522. doi: 10.1371/journal.pcbi.1005522. [PMC free article] [PubMed] [Cross Ref]
64. Schmitt S, Kuhn D, Klebe G. A new method to detect related function among proteins independent of sequence and fold homology. J. Mol. Biol. 2002;323:387–406. doi: 10.1016/S0022-2836(02)00811-2. [PubMed] [Cross Ref]
65. Konc J, Janežič D. An improved branch and bound algorithm for the maximum clique problem. MATCH Commun. Math. Comput. Chem. 2007;58:569–590.
66. Karlin S, Altschul SF. Methods for assessing the statistical significance of molecular sequence features by using general scoring schemes. Proc. Natl. Acad. Sci. USA. 1990;87:2264–2268. doi: 10.1073/pnas.87.6.2264. [PubMed] [Cross Ref]
67. Alstott, J., Bullmore, E. & Plenz, D. Powerlaw: A python package for analysis of heavy-tailed distributions. PLoS One9 (2014). [PMC free article] [PubMed]
68. Clauset A, Shalizi CR, Newman MEJ. Power-law distributions in empirical data. SIAM Rev. 2007;51:661–703. doi: 10.1137/070710111. [Cross Ref]
69. Orengo CA, et al. CATH–a hierarchic classification of protein domain structures. Structure. 1997;5:1093–1108. doi: 10.1016/S0969-2126(97)00260-8. [PubMed] [Cross Ref]
70. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. R. Statisitical Soc. Ser. B. 1995;57:289–300.

Articles from Scientific Reports are provided here courtesy of Nature Publishing Group