Email updates

Keep up to date with the latest news and content from BMC Structural Biology and BioMed Central.

Open Access Methodology article

Fast dynamics perturbation analysis for prediction of protein functional sites

Dengming Ming12, Judith D Cohn13 and Michael E Wall134*

Author Affiliations

1 Computer, Computational, and Statistical Scienes Division, Los Alamos National Laboratory, Los Alamos, New Mexico, USA

2 School of Life Sciences, Nanjing University, Nanjing, Jiangsu Province, China

3 Bioscience Division, Los Alamos National Laboratory, Los Alamos, New Mexico, USA

4 Center for Nonlinear Studies, Los Alamos National Laboratory, Los Alamos, New Mexico, USA

For all author emails, please log on.

BMC Structural Biology 2008, 8:5  doi:10.1186/1472-6807-8-5


The electronic version of this article is the complete one and can be found online at: http://www.biomedcentral.com/1472-6807/8/5


Received:18 October 2007
Accepted:30 January 2008
Published:30 January 2008

© 2008 Ming et al; licensee BioMed Central Ltd.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Background

We present a fast version of the dynamics perturbation analysis (DPA) algorithm to predict functional sites in protein structures. The original DPA algorithm finds regions in proteins where interactions cause a large change in the protein conformational distribution, as measured using the relative entropy Dx. Such regions are associated with functional sites.

Results

The Fast DPA algorithm, which accelerates DPA calculations, is motivated by an empirical observation that Dx in a normal-modes model is highly correlated with an entropic term that only depends on the eigenvalues of the normal modes. The eigenvalues are accurately estimated using first-order perturbation theory, resulting in a N-fold reduction in the overall computational requirements of the algorithm, where N is the number of residues in the protein. The performance of the original and Fast DPA algorithms was compared using protein structures from a standard small-molecule docking test set. For nominal implementations of each algorithm, top-ranked Fast DPA predictions overlapped the true binding site 94% of the time, compared to 87% of the time for original DPA. In addition, per-protein recall statistics (fraction of binding-site residues that are among predicted residues) were slightly better for Fast DPA. On the other hand, per-protein precision statistics (fraction of predicted residues that are among binding-site residues) were slightly better using original DPA. Overall, the performance of Fast DPA in predicting ligand-binding-site residues was comparable to that of the original DPA algorithm.

Conclusion

Compared to the original DPA algorithm, the decreased run time with comparable performance makes Fast DPA well-suited for implementation on a web server and for high-throughput analysis.

Background

Prediction of protein functional sites is a key aspect of protein function prediction [1], and can be an important step in identifying small-molecule interactions for drug discovery [2]. It can also potentially be used as a pre-processing step to reduce the search space in computational docking algorithms. There are many methods to predict functional sites–here we emphasize those that make use of analysis of protein structure and dynamics. Existing protein structure analysis methods are based on diverse principles, including: association of functional sites with surface clefts that have extreme values of volume [3-6] or other shape descriptors [7-11]; identifying spatial clusters of methyl probes that exhibit energetically favorable interactions with the protein [12]; association of functional sites with charged surface residues either in unfavorable electrostatic environments [13] or with anomalous predicted pH titration curves [14]; identifying spatial clusters of residues whose diversity appears to be correlated with changes in protein function [15,16]; defining structural features (e.g. motifs) associated with functional sites [17-22]; identifying residues that are on average close to other residues in the protein (closeness centrality) [23-25]; and machine-learning prediction of functional sites/residues using sequence, structure, and chemical features from training sets [26-28]. Principles of methods that consider protein dynamics include association of functional sites with: hinge regions [29,30]; regions where the harmonic vibrations are largely determined by high-frequency modes [31]; intrinsically disordered regions that are highly mobile in the absence of a molecular interaction partner [32]; and residues where mutations cause a large change in the couplings of local perturbations to remote, local changes in the distribution of folded vs. unfolded states of the protein [33]. Information from complementary methods may be integrated for functional site prediction [34,35].

We recently developed an additional approach to prediction of protein functional sites that is based on analysis of protein dynamics [36-39]. To help motivate the approach, we note that cellular functions are regulated by molecular interactions that alter protein activity. To enable such control, protein activity, and therefore protein conformational distributions, must be susceptible to alteration by molecular interactions at functional sites. In other words, protein activity should be controllable by allosteric effects (allostery).

Weber [40] recognized the importance of considering changes in the full conformational distribution to understand allostery, as opposed to considering mechanistic changes among discrete, well-defined structural states in earlier models due to Monod, Wyman, and Changeux [41]; and Koshland, Nemethy, and Filmer [42]. Weber's perspective is well-aligned with more recent emphases on the need to consider allostery from a global thermodynamic/statistical perspective [43,36-39,33,45]. It is also well-aligned with modern rate theories based on the control of protein activity by dynamical transitions among conformational substates [46], as originally suggested by spectroscopic assays of ligand-binding at low-temperature [47,48].

Given the above considerations, we hypothesized that protein functional sites might tend to evolve at control points where interactions cause a large change in the protein conformational distribution [36]. To test this hypothesis, we developed a method called dynamics perturbation analysis (DPA) to quantify changes in protein conformational distributions due to molecular interactions [36,37], examined 305 protein structures from the GOLD [49] docking test set [38], and found that interactions at small-molecule binding sites cause a relatively large change in protein vibrations.

Motivated by these results, we developed a DPA-based algorithm that successfully predicts small-molecule binding sites at locations where interactions cause a large change in protein vibrations [38]. This method was evaluated in Ref. [38] using 305 proteins in the GOLD [49] docking test set of protein-ligand structures. For the test, only the top-ranked functional site was selected and was used to predict the location of the ligand-binding site. This is a relatively strict requirement; in other published methods for predicting functional sites [11], performance often is evaluated by allowing for any of several predicted functional sites to overlap a known ligand-binding site. The method produced at least one predicted functional site for 287 of the 305 proteins in the test set. In 87% of cases (250 proteins), at least one predicted residue was in the ligand-binding site. The recall of binding-site residues (percentage of binding-site residues found among the predicted residues) was at least 30% for 80% of cases, and was at least 50% for 76% of the cases. The precision of the predicted residues (percentage of predicted residues found among the binding-site residues) was at least 30% for 68% of the cases, and was at least 50% for 44% of the cases. The statistical significance of the overlaps was assessed using a null model in which surface residues were randomly selected. Using the null model, a P-value was calculated to evaluate predictions for the 250 proteins in which at least one predicted residue was in the ligand-binding site. The P-value estimated the probability of obtaining a precision at least as high as the observed precision by randomly selecting surface residues [38]. For 87% of the cases, the P-value was 10-3 or smaller, indicating a statistically significant overlap. The performance of the DPA method compared favorably to that of a cleft analysis method for predicting ligand-binding residues.

The original DPA algorithm is a highly innovative approach that performs well. However, the computational requirements limit the utility of the original method. For example, it takes about an hour to analyze a 150-residue protein domain using DPA, and the method doesn't scale well to larger systems. Here, we report an improved algorithm based on use of first-order perturbation theory that will facilitate the use of DPA in high-throughput scenarios and increase its utility, e.g., for web server applications. The algorithm, called Fast DPA, enables a dramatic decrease in the time required to predict protein functional sites, with performance that is comparable to the original DPA algorithm.

Methods

Dynamics perturbation analysis

Our overall approach for predicting functional sites is based on a method called dynamics perturbation analysis (DPA) [36,38,37]. In DPA, a protein is decorated with M surface points that interact with neighboring protein atoms, as illustrated for Protein Data Bank entry 1JEF [50] in Fig. 1. The protein conformational distribution P(x) is calculated in the absence of any surface points, and M protein conformational distributions P(m)(x) are calculated for the protein interacting with each point m. The conformational distributions are calculated using a coarse-grained model of molecular vibrations, and the distributions P(m)(x) are calculated from models of the protein in complex with each surface point. The relative entropy, or Kullback-Leibler divergence [51], <a onClick="popup('http://www.biomedcentral.com/1472-6807/8/5/mathml/M1','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1472-6807/8/5/mathml/M1">View MathML</a> between P(x) and P(m)(x) is calculated for each point m, and is used as a measure of the change in the protein conformational distribution upon interacting with point m:

thumbnailFigure 1. Application of Dynamics Perturbation Analysis (DPA) to predict protein functional sites. Left. In this example, the surface of lysozyme (PDB entry 1JEF [50], yellow cartoon) is decorated with test points (533 spheres at a density of 1 point per Å2), and the degree to which the test points individually perturb the protein conformational distribution is calculated (temperature-coded coloring of the spheres). A tri-NAG molecule (purple wireframe) binds in the active site. Warm-colored spheres indicate where the perturbation is large. Center. Points where the perturbation is largest are selected and clustered (green spheres). Right. Cα atoms within 6 Å of the DPA cluster are selected, and the associated residues define the predicted functional site (16 residues). For comparison, Cα atoms within 6 Å of the tri-NAG are selected; we use the associated residues to define the actual functional site (7 residues). The overlapping residues (6 residues) are shown in orange; there are 10 predicted residues that do not exactly match the functional site (green), and there is 1 functional site residue that is not among the predicted residues (purple, in the helix on the right hand side).

<a onClick="popup('http://www.biomedcentral.com/1472-6807/8/5/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1472-6807/8/5/mathml/M2">View MathML</a>

(1)

In the present case (unlike in other useful biological applications [52-56]), the relative entropy is not just an ad hoc measure; rather, it has real biophysical significance [39,57]: <a onClick="popup('http://www.biomedcentral.com/1472-6807/8/5/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1472-6807/8/5/mathml/M3">View MathML</a>, where T is the temperature and kB is Boltzmann's constant, is the free energy required to change the protein conformational distribution from an equilibrium distribution P(x) to a non-equilibrium distribution P(m)(x).

Thus far, DPA calculations have most often been performed using a simple model of protein vibrations–the elastic network model (ENM) [58-61]. In the ENM, Cα atoms are extracted from an atomic model of a protein, and an interaction network is generated by connecting springs between all atom pairs (i, j) separated by a distance less than or equal to a cutoff distance rc. Each spring has the same force constant γ, is aligned with the separation between the connected atoms, and has an equilibrium length equal to the distance dij between the atoms in the initial model. Thus, the potential energy is given by <a onClick="popup('http://www.biomedcentral.com/1472-6807/8/5/mathml/M4','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1472-6807/8/5/mathml/M4">View MathML</a>, where εij = 1 if atoms i and j are connected, and εij = 0 otherwise. The interaction between the protein and a surface point m is modeled by connecting springs of force constant γs between the surface point and all protein atoms within a cutoff distance rs of the surface point. The protein coordinates are not modified in modeling the interaction. The dynamics are defined using normal mode analysis of the model. In this model, the reference distribution P(x) is given by

<a onClick="popup('http://www.biomedcentral.com/1472-6807/8/5/mathml/M5','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1472-6807/8/5/mathml/M5">View MathML</a>

(2)

In Eq. (2), N is the number of atoms in the protein; x0 is the equilibrium structure; and λi and vi are the ith eigenvalue and eigenvector of the Hessian <a onClick="popup('http://www.biomedcentral.com/1472-6807/8/5/mathml/M6','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1472-6807/8/5/mathml/M6">View MathML</a>. The perturbed distribution P(m)(x) is similar to Eq. (2), but substituting the eigenvalues and eigenvectors <a onClick="popup('http://www.biomedcentral.com/1472-6807/8/5/mathml/M7','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1472-6807/8/5/mathml/M7">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1472-6807/8/5/mathml/M8','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1472-6807/8/5/mathml/M8">View MathML</a> of the pseudo-Hessian <a onClick="popup('http://www.biomedcentral.com/1472-6807/8/5/mathml/M9','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1472-6807/8/5/mathml/M9">View MathML</a> for λi and vi. <a onClick="popup('http://www.biomedcentral.com/1472-6807/8/5/mathml/M9','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1472-6807/8/5/mathml/M9">View MathML</a> is derived from the full Hessian H(m) for the protein model in the presence of the surface point m:

<a onClick="popup('http://www.biomedcentral.com/1472-6807/8/5/mathml/M10','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1472-6807/8/5/mathml/M10">View MathML</a>

(3)

The sub-matrix <a onClick="popup('http://www.biomedcentral.com/1472-6807/8/5/mathml/M11','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1472-6807/8/5/mathml/M11">View MathML</a> couples the protein coordinates, the sub-matrix <a onClick="popup('http://www.biomedcentral.com/1472-6807/8/5/mathml/M12','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1472-6807/8/5/mathml/M12">View MathML</a> couples the test-point coordinates, and the submatrix G(m) couples the protein to the test point. In terms of these matrices, <a onClick="popup('http://www.biomedcentral.com/1472-6807/8/5/mathml/M9','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1472-6807/8/5/mathml/M9">View MathML</a> is given by [37]

<a onClick="popup('http://www.biomedcentral.com/1472-6807/8/5/mathml/M13','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1472-6807/8/5/mathml/M13">View MathML</a>

(4)

Using expressions for P(x) and P(m)(x), Eq. (1) becomes [36,37]

<a onClick="popup('http://www.biomedcentral.com/1472-6807/8/5/mathml/M14','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1472-6807/8/5/mathml/M14">View MathML</a>

(5)

The first six modes involve zero eigenvalues and are ignored in the sums. Equation (5) is the central equation that enables DPA.

To use DPA to predict functional sites, we make use of the fact that, empirically, the distribution of y = <a onClick="popup('http://www.biomedcentral.com/1472-6807/8/5/mathml/M1','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1472-6807/8/5/mathml/M1">View MathML</a> values on the surface of a protein calculated using Eq. (5) is observed to obey an extreme value distribution (Fig. 2),

thumbnailFigure 2. Distribution of <a onClick="popup('http://www.biomedcentral.com/1472-6807/8/5/mathml/M1','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1472-6807/8/5/mathml/M1">View MathML</a> values for 4859 points on the surface of lysozyme 1JEF (the number of points was increased in this case to evaluate the fit). The distribution is well-fit by an extreme value distribution (Eq. (6)) with parameters μ = 23.07 and β = 8.45 (solid line). By examining the cumulative distribution (dashed line), the fit is used to find surface points that lie within the upper 96% of the distribution; these points are used to predict functional sites.

<a onClick="popup('http://www.biomedcentral.com/1472-6807/8/5/mathml/M15','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1472-6807/8/5/mathml/M15">View MathML</a>

(6)

First, DPA is performed on a protein and the distribution of <a onClick="popup('http://www.biomedcentral.com/1472-6807/8/5/mathml/M1','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1472-6807/8/5/mathml/M1">View MathML</a> values is modeled using Eq. (6). Points with <a onClick="popup('http://www.biomedcentral.com/1472-6807/8/5/mathml/M1','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1472-6807/8/5/mathml/M1">View MathML</a> values in the upper 96% of the modeled distribution are selected and are spatially clustered. The clusters are ranked according to the mean value of <a onClick="popup('http://www.biomedcentral.com/1472-6807/8/5/mathml/M1','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1472-6807/8/5/mathml/M1">View MathML</a> within the cluster, and all clusters are considered to be potentially associated with a functional site. Finally, residues in the neighborhood of the clusters are selected and form the basis for functional site predictions.

Fast dynamics perturbation analysis

Fast DPA is based on a simple empirical observation: for dynamics defined by normal modes, the total value of Dx in Eq. (5) is highly correlated with just the first (entropic) term,

<a onClick="popup('http://www.biomedcentral.com/1472-6807/8/5/mathml/M16','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1472-6807/8/5/mathml/M16">View MathML</a>

(7)

Hereafter we refer to <a onClick="popup('http://www.biomedcentral.com/1472-6807/8/5/mathml/M17','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1472-6807/8/5/mathml/M17">View MathML</a> simply as <a onClick="popup('http://www.biomedcentral.com/1472-6807/8/5/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1472-6807/8/5/mathml/M18">View MathML</a>. Observation of this correlation motivates the use of <a onClick="popup('http://www.biomedcentral.com/1472-6807/8/5/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1472-6807/8/5/mathml/M18">View MathML</a> as a surrogate for Dx in DPA, and, because <a onClick="popup('http://www.biomedcentral.com/1472-6807/8/5/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1472-6807/8/5/mathml/M18">View MathML</a> only involves eigenvalues, creates an avenue for accelerating DPA. The acceleration arises because the eigenvalues of the normal models of the protein in the presence of test points are well-approximated using first order perturbation theory. In this approximation, the pseudo-Hessian <a onClick="popup('http://www.biomedcentral.com/1472-6807/8/5/mathml/M9','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1472-6807/8/5/mathml/M9">View MathML</a> of the protein in the presence of point m is written as the Hessian H of the protein in the absence of the ligand plus a perturbation term δ<a onClick="popup('http://www.biomedcentral.com/1472-6807/8/5/mathml/M9','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1472-6807/8/5/mathml/M9">View MathML</a>:

<a onClick="popup('http://www.biomedcentral.com/1472-6807/8/5/mathml/M19','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1472-6807/8/5/mathml/M19">View MathML</a>

(8)

where the expression for <a onClick="popup('http://www.biomedcentral.com/1472-6807/8/5/mathml/M9','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1472-6807/8/5/mathml/M9">View MathML</a> is as in previous studies [37,38]. To estimate the eigenvalues of <a onClick="popup('http://www.biomedcentral.com/1472-6807/8/5/mathml/M9','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1472-6807/8/5/mathml/M9">View MathML</a>, we use the canonical first-order perturbation theory expression,

<a onClick="popup('http://www.biomedcentral.com/1472-6807/8/5/mathml/M20','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1472-6807/8/5/mathml/M20">View MathML</a>

(9)

where λi is the ith eigenvalue of H.

The Fast DPA algorithm is the same as the original DPA algorithm, except instead of using values of Dx, the analysis is based on values of <a onClick="popup('http://www.biomedcentral.com/1472-6807/8/5/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1472-6807/8/5/mathml/M18">View MathML</a> estimated using perturbation theory. (It is possible to evaluate all terms in Eq. (5) using first-order perturbation theory, but doing so would not accelerate the method because the computational cost is comparable to that of solving the full eigenvalue problem in original DPA.)

Implementation of Fast DPA

Our implementation of DPA and Fast DPA here follows our previous implementation of DPA for functional site prediction [38]. Given an input PDB structure, MSMS [62] was run with a 1.5 Å probe radius and a triangulation density of 1 vertex per Å2 to generate test points on the surface of the protein. As when using original DPA to predict functional sites, perturbations were calculated using every other point in the MSMS output (we also tried using every point, but this led to decreased performance in the precision measures). The cutoff rc for interactions between protein Cα atoms was 8.5 Å. For some proteins, this cutoff yielded more than six zero-frequency modes, indicating that the network of springs was too sparse (for example, if only one spring connects two domains, then free rotations about the spring yield two additional zero-frequency modes). In these cases, the connectivity of the elastic network model was increased by incrementing rc in 1 Å steps until the additional zero-frequency modes were eliminated. The cutoff rs for interactions between a test point and the protein was 14 Å, and the interaction strength between a test point and protein atoms was γs = 12γ, or 12 times the strength of the interaction between two protein atoms. Results are independent of the value of γ.

Implementation of functional site prediction using DPA

To predict functional sites, the distribution of y = <a onClick="popup('http://www.biomedcentral.com/1472-6807/8/5/mathml/M1','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1472-6807/8/5/mathml/M1">View MathML</a> values was fit using Eq. (6). Points with <a onClick="popup('http://www.biomedcentral.com/1472-6807/8/5/mathml/M1','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1472-6807/8/5/mathml/M1">View MathML</a> values in the upper 96% of the distribution were selected and spatially clustered using the OPTICS algorithm [63] with a distance threshold of 6 Å and a minimum of 3 points per cluster. Cα atoms within 6 Å of any point in a cluster were selected and were used to define predicted functional sites. The sites were ranked according to the mean value of <a onClick="popup('http://www.biomedcentral.com/1472-6807/8/5/mathml/M1','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1472-6807/8/5/mathml/M1">View MathML</a> within the corresponding cluster of points. Only the top-ranked predicted site was used for the evaluation of performance described below.

Results and Discussion

Results that motivate Fast DPA

To motivate the use of <a onClick="popup('http://www.biomedcentral.com/1472-6807/8/5/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1472-6807/8/5/mathml/M18">View MathML</a> instead of Dx for DPA, we analyzed proteins from the GOLD test set. We found that Dx is highly correlated with <a onClick="popup('http://www.biomedcentral.com/1472-6807/8/5/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1472-6807/8/5/mathml/M18">View MathML</a> for these cases; Fig. 3 illustrates the agreement for four proteins. This is not a trivial result mathematically (see Eqs. (5) and (7))–it means that <a onClick="popup('http://www.biomedcentral.com/1472-6807/8/5/mathml/M21','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1472-6807/8/5/mathml/M21">View MathML</a> is highly correlated with <a onClick="popup('http://www.biomedcentral.com/1472-6807/8/5/mathml/M22','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1472-6807/8/5/mathml/M22">View MathML</a>.

thumbnailFigure 3. Values of Dx (y-axis) and <a onClick="popup('http://www.biomedcentral.com/1472-6807/8/5/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1472-6807/8/5/mathml/M18">View MathML</a> (x-axis) calculated using original DPA are plotted for four PDB entries (values of the Pearson correlation, C, between the two, are listed here parenthetically): a) 1AEC [65], from an actinidin-E-64 complex (C = 0.988); b) 1FKI [66], from a FKBP complex (0.989); c) 1JEF [50], from a lysozyme complex (0.992); and d) 1STP [67], from a biotin complex (0.989).

To motivate the use of perturbation theory to estimate <a onClick="popup('http://www.biomedcentral.com/1472-6807/8/5/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1472-6807/8/5/mathml/M18">View MathML</a>, we compared the true eigenvalues to those estimated using perturbation theory for proteins in the GOLD test set. Because in our model the strength of the spring that connects the test points to the protein is 12 times the strength of the spring that connects protein atoms to each other (Methods), it was not obvious that first-order perturbation theory would yield reasonable estimates of eigenvalues. However, we had hoped for success based on the fact that we were only adding a single test point to the model, compared to, typically, O(100) protein Cα atoms. As illustrated for lysozyme in Fig. 4, we did find that Eq. (9) approximates well the true eigenvalues obtained by diagonalization of H(m). Finally, we found that Dx calculated using original DPA was highly correlated with <a onClick="popup('http://www.biomedcentral.com/1472-6807/8/5/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1472-6807/8/5/mathml/M18">View MathML</a> calculated using Fast DPA, as illustrated for four proteins in Fig. 5.

thumbnailFigure 4. Eigenvalues (used for calculation of <a onClick="popup('http://www.biomedcentral.com/1472-6807/8/5/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1472-6807/8/5/mathml/M18">View MathML</a>) that are estimated using perturbation theory (filled triangles) are a good approximation to the true eigenvalues of a lysozyme elastic network model (open circles).

thumbnailFigure 5. Values of Dx calculated using original DPA (y-axis) and <a onClick="popup('http://www.biomedcentral.com/1472-6807/8/5/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1472-6807/8/5/mathml/M18">View MathML</a> calculated using Fast DPA (x-axis) are plotted for four PDB entries (values of the Pearson correlation between the two are listed here parenthetically): a) 1AEC (0.981); b) 1FKI (0.982); c) 1JEF (0.981); d) 1STP (0.980).

Evaluation of Fast DPA for prediction of functional sites

The above results motivated us to develop the Fast DPA algorithm for prediction of protein functional sites (Methods). Through use of first-order perturbation theory, Fast DPA replaces matrix diagonalization by matrix-vector multiplication for each test point (Eq. (9)). Because matrix diagonalization requires O(N3) operations, and matrix-vector multiplication requires O(N2) operations, we expected Fast DPA to run N-fold faster than the original DPA. We found this to be the case (Fig. 6): the original DPA scales roughly as N3.45, while fast DPA scales roughly as N2.29, yielding a factor of N1.16 decrease in the time required to perform Fast DPA vs. DPA (here, N is the number of residues in the protein).

thumbnailFigure 6. Comparison of run times for DPA (upwards-pointing triangles) vs. Fast DPA (downwards-pointing triangles) for various protein sizes. The inset shows the ratio of run times for various protein sizes.

Because Dx calculated using original DPA and <a onClick="popup('http://www.biomedcentral.com/1472-6807/8/5/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1472-6807/8/5/mathml/M18">View MathML</a> calculated using Fast DPA are highly correlated (Fig. 5), we expected the performance of Fast DPA in predicting functional site residues to be comparable to that of the original DPA. We analyzed the performance of the algorithm on the 305-protein GOLD [49] test set, which was used to evaluate the original DPA algorithm [38]. Each prediction has an associated recall (fraction of residues in the binding site that are among those in the rank-1 prediction) and precision (fraction of rank-1 predicted residues that are among those in the binding site). To evaluate performance statistically, we use (1) the fraction of binding sites for which the recall is greater than or equal to a minimum value, and (2) the fraction of fraction of rank-1 predictions for which the precision is greater than or equal to a minimum value.

Figure 7 compares the performance of Fast DPA using different thresholds of the extreme value distribution, and is equivalent to Fig. 8 in [38]. The nominal threshold of 0.96 indicated in this figure is equivalent to that chosen for original DPA. Fig. 8 compares the performance of Fast DPA with original DPA for different thresholds. When the threshold is 0.96 or smaller, the recall statistics of Fast DPA tend to be better, and the precision statistics of original DPA tend to be better. When the threshold is 0.97 or higher, original DPA outperforms Fast DPA in both precision and recall statistics.

thumbnailFigure 7. Comparison of Fast DPA performance using different thresholds of the extreme value distribution (Eq. (6)). The y-axis is either the fraction of proteins for which a prediction is made (squares), the fraction of binding sites with a recall of at least 0.5 (circles), or the fraction of predictions with a precision of at least 0.5 (triangles). The threshold is indicated on the x-axis; the 0.96 threshold used for Figs. 9 and 10 is indicated using a vertical dashed line.

thumbnailFigure 8. Comparison of Fast DPA vs. original DPA precision and recall statistics at different thresholds of the extreme value distribution (Eq. (6)). The curves are similar to precision-recall curves: the y-axis is the fraction of binding sites with a recall of at least 0.5, and the x-axis is the fraction of binding sites with a precision of at least 0.5. Fast DPA values are indicated using open squares, and original DPA is indicated using filled squares. Points corresponding to a threshold of 0.96 are indicated using arrows.

At the nominal threshold value of 0.96, the performance of Fast DPA is comparable to that of original DPA. At this threshold, original DPA yielded 287 rank-1 predictions for the test set (rate of 94%), whereas Fast DPA yielded 267 rank-1 predictions (rate of 86%) (Table 1). However, Fast DPA makes 251 predictions that have at least one residue that overlaps the binding site, while original DPA makes 250 such predictions, yielding a higher rate of locating binding sites for rank-1 Fast DPA predictions (94%) than for original DPA (87%) (Table 1). The recall statistics tend to be a bit better for Fast DPA (Table 1, Fig. 9), and the precision statistics tend to be better for original DPA (Table 1, Fig. 10).

Table 1. Performance statistics for Fast DPA and original DPA using a threshold of 0.96

thumbnailFigure 9. Comparison of recall of binding-site residues using DPA vs. Fast DPA for 287 (number of predictions using DPA) or 267 (number of predictions using Fast DPA) proteins in the 305-protein GOLD test set. The y-axis indicates the fraction of proteins with a recall at least as high as the value on the x-axis (y-values should be read from the top of each step).

thumbnailFigure 10. Comparison of precision of predicted residues using DPA vs. Fast DPA (see also Fig. 9). The y-axis indicates the fraction of proteins with a precision at least as high as the value on the x-axis (y-values should be read from the top of each step).

Conclusion

Use of Fast DPA enables functional site predictions to be performed N-fold faster than original DPA, with comparable performance in predicting residues in functional sites. The acceleration will facilitate optimization of Fast DPA for functional site predictions. Calculations that once took hours using DPA now may be performed in a matter of minutes, making practical the use of DPA via a web server. Indeed, high-throughput analysis using Fast DPA has already produced over 60,000 predicted functional sites for about 50,000 protein domains in the SCOP [64] database (J.D. Cohn, D. Ming, and M.E. Wall, in preparation). These predictions will provide a rich source of information for developing hypotheses concerning mechanisms of protein function.

Authors' contributions

DM implemented the Fast DPA algorithm, tested its performance, and helped to draft the manuscript. JC provided assistance with databases and automation. MW conceived of the study, coordinated the work, and drafted the manuscript. All authors read and approved the final manuscript.

Acknowledgements

Supported by the US Department of Energy through contract DE-AC52-06NA25396. We thank James Faeder for reading the manuscript.

References

  1. Ofran Y, Punta M, Schneider R, Rost B: Beyond annotation transfer by homology: novel protein-function prediction methods to assist drug discovery.

    Drug Discov Today 2005, 10(21):1475-1482. PubMed Abstract | Publisher Full Text OpenURL

  2. Campbell SJ, Gold ND, Jackson RM, Westhead DR: Ligand binding: functional site location, similarity and docking.

    Curr Opin Struct Biol 2003, 13(3):389-395. PubMed Abstract | Publisher Full Text OpenURL

  3. Dundas J, Ouyang Z, Tseng J, Binkowski A, Turpaz Y, Liang J: CASTp: computed atlas of surface topography of proteins with structural and topographical mapping of functionally annotated residues.

    Nucleic Acids Res 2006, 34(Web Server issue):W116-8. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  4. Hendlich M, Rippmann F, Barnickel G: LIGSITE: automatic and efficient detection of potential small molecule-binding sites in proteins.

    J Mol Graph Model 1997, 15(6):359-63, 389. PubMed Abstract | Publisher Full Text OpenURL

  5. Glaser F, Morris RJ, Najmanovich RJ, Laskowski RA, Thornton JM: A method for localizing ligand binding pockets in protein structures.

    Proteins 2006, 62(2):479-488. PubMed Abstract | Publisher Full Text OpenURL

  6. Laskowski RA: SURFNET: a program for visualizing molecular surfaces, cavities, and intermolecular interactions.

    J Mol Graph 1995, 13(5):323-30, 307-8. PubMed Abstract OpenURL

  7. Coleman RG, Burr MA, Souvaine DL, Cheng AC: An intuitive approach to measuring protein surface curvature.

    Proteins 2005, 61(4):1068-1074. PubMed Abstract | Publisher Full Text OpenURL

  8. Coleman RG, Sharp KA: Travel depth, a new shape descriptor for macromolecules: application to ligand binding.

    J Mol Biol 2006, 362(3):441-458. PubMed Abstract | Publisher Full Text OpenURL

  9. Hendrix DK, Kuntz ID: Surface solid angle-based site points for molecular docking.

    Pac Symp Biocomput 1998, 317-326. PubMed Abstract OpenURL

  10. Nayal M, Honig B: On the nature of cavities on protein surfaces: application to the identification of drug-binding sites.

    Proteins 2006, 63(4):892-906. PubMed Abstract | Publisher Full Text OpenURL

  11. Xie L, Bourne PE: A robust and efficient algorithm for the shape description of protein structures and its application in predicting ligand binding sites.

    BMC Bioinformatics 2007, 8 Suppl 4:S9. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  12. Laurie AT, Jackson RM: Q-SiteFinder: an energy-based method for the prediction of protein-ligand binding sites.

    Bioinformatics 2005, 21(9):1908-1916. PubMed Abstract | Publisher Full Text OpenURL

  13. Elcock AH: Prediction of functionally important residues based solely on the computed energetics of protein structure.

    J Mol Biol 2001, 312(4):885-896. PubMed Abstract | Publisher Full Text OpenURL

  14. Ondrechen MJ, Clifton JG, Ringe D: THEMATICS: a simple computational predictor of enzyme function from structure.

    Proc Natl Acad Sci U S A 2001, 98(22):12473-12478. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  15. Lichtarge O, Bourne HR, Cohen FE: An evolutionary trace method defines binding surfaces common to protein families.

    J Mol Biol 1996, 257(2):342-358. PubMed Abstract | Publisher Full Text OpenURL

  16. Yao H, Mihalek I, Lichtarge O: Rank information: a structure-independent measure of evolutionary trace quality that improves identification of protein functional sites.

    Proteins 2006, 65(1):111-123. PubMed Abstract | Publisher Full Text OpenURL

  17. Wallace AC, Borkakoti N, Thornton JM: TESS: a geometric hashing algorithm for deriving 3D coordinate templates for searching structural databases. Application to enzyme active sites.

    Protein Sci 1997, 6(11):2308-2323. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  18. Shulman-Peleg A, Nussinov R, Wolfson HJ: Recognition of functional sites in protein structures.

    J Mol Biol 2004, 339(3):607-633. PubMed Abstract | Publisher Full Text OpenURL

  19. Stark A, Russell RB: Annotation in three dimensions. PINTS: Patterns in Non-homologous Tertiary Structures.

    Nucleic Acids Res 2003, 31(13):3341-3344. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  20. Stark A, Shkumatov A, Russell RB: Finding functional sites in structural genomics proteins.

    Structure 2004, 12(8):1405-1412. PubMed Abstract | Publisher Full Text OpenURL

  21. Liang MP, Brutlag DL, Altman RB: Automated construction of structural motifs for predicting functional sites on protein structures.

    Pac Symp Biocomput 2003, 204-215. PubMed Abstract OpenURL

  22. Barker JA, Thornton JM: An algorithm for constraint-based structural template matching: application to 3D templates with statistical analysis.

    Bioinformatics 2003, 19(13):1644-1649. PubMed Abstract | Publisher Full Text OpenURL

  23. Amitai G, Shemesh A, Sitbon E, Shklar M, Netanely D, Venger I, Pietrokovski S: Network analysis of protein structures identifies functional residues.

    J Mol Biol 2004, 344(4):1135-1146. PubMed Abstract | Publisher Full Text OpenURL

  24. Thibert B, Bredesen DE, del Rio G: Improved prediction of critical residues for protein function based on network and phylogenetic analyses.

    BMC Bioinformatics 2005, 6:213. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  25. Chea E, Livesay DR: How accurate and statistically robust are catalytic site predictions based on closeness centrality?

    BMC Bioinformatics 2007, 8:153. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  26. Ofran Y, Rost B: ISIS: interaction sites identified from sequence.

    Bioinformatics 2007, 23(2):e13-6. PubMed Abstract | Publisher Full Text OpenURL

  27. Gutteridge A, Bartlett GJ, Thornton JM: Using a neural network and spatial clustering to predict the location of active sites in enzymes.

    J Mol Biol 2003, 330(4):719-734. PubMed Abstract | Publisher Full Text OpenURL

  28. Wei L, Altman RB: Recognizing complex, asymmetric functional sites in protein structures using a Bayesian scoring function.

    J Bioinform Comput Biol 2003, 1(1):119-138. PubMed Abstract | Publisher Full Text OpenURL

  29. Ma B, Wolfson HJ, Nussinov R: Protein functional epitopes: hot spots, dynamics and combinatorial libraries.

    Curr Opin Struct Biol 2001, 11(3):364-369. PubMed Abstract | Publisher Full Text OpenURL

  30. Yang LW, Bahar I: Coupling between catalytic site and collective dynamics: a requirement for mechanochemical activity of enzymes.

    Structure 2005, 13(6):893-904. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  31. Haliloglu T, Keskin O, Ma B, Nussinov R: How similar are protein folding and protein binding nuclei? Examination of vibrational motions of energy hot spots and conserved residues.

    Biophys J 2005, 88(3):1552-1559. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  32. Radivojac P, Iakoucheva LM, Oldfield CJ, Obradovic Z, Uversky VN, Dunker AK: Intrinsic disorder and functional proteomics.

    Biophys J 2007, 92(5):1439-1456. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  33. Liu T, Whitten ST, Hilser VJ: Functional residues serve a dominant role in mediating the cooperativity of the protein ensemble.

    Proc Natl Acad Sci U S A 2007, 104(11):4347-4352. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  34. Rossi A, Marti-Renom MA, Sali A: Localization of binding sites in protein structures by optimization of a composite scoring function.

    Protein Sci 2006, 15(10):2366-2380. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  35. Petrova NV, Wu CH: Prediction of catalytic residues using Support Vector Machine with selected protein sequence and structural properties.

    BMC Bioinformatics 2006, 7:312. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  36. Ming D, Wall ME: Quantifying allosteric effects in proteins.

    Proteins 2005, 59(4):697-707. PubMed Abstract | Publisher Full Text OpenURL

  37. Ming D, Wall ME: Allostery in a coarse-grained model of protein dynamics.

    Phys Rev Lett 2005, 95:198103. PubMed Abstract | Publisher Full Text OpenURL

  38. Ming D, Wall ME: Interactions in native binding sites cause a large change in protein dynamics.

    J Mol Biol 2006, 358:213-223. PubMed Abstract | Publisher Full Text OpenURL

  39. Wall ME: Ligand binding, protein fluctuations, and allosteric free energy.

    AIP Conf Proc 2006, 851:16-33. OpenURL

  40. Weber G: Ligand binding and internal equilibria in proteins.

    Biochemistry 1972, 11(5):864-878. PubMed Abstract OpenURL

  41. Monod J, Wyman J, Changeux JP: On the nature of allosteric transitions: a plausible model.

    J Mol Biol 1965, 12:88-118. PubMed Abstract OpenURL

  42. Koshland DE Jr., Nemethy G, Filmer D: Comparison of experimental binding data and theoretical models in proteins containing subunits.

    Biochemistry 1966, 5(1):365-385. PubMed Abstract OpenURL

  43. Pan H, Lee JC, Hilser VJ: Binding sites in Escherichia coli dihydrofolate reductase communicate by modulating the conformational ensemble.

    Proc Natl Acad Sci U S A 2000, 97(22):12020-12025. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  44. Gunasekaran K, Ma B, Nussinov R: Is allostery an intrinsic property of all dynamic proteins?

    Proteins 2004, 57:433-443. PubMed Abstract | Publisher Full Text OpenURL

  45. Hilser VJ, Thompson EB: Intrinsic disorder as a mechanism to optimize allosteric coupling in proteins.

    Proc Natl Acad Sci U S A 2007, 104(20):8311-8315. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  46. Frauenfelder H, Wolynes PG: Rate theories and puzzles of hemeprotein kinetics.

    Science 1985, 229(4711):337-345. PubMed Abstract | Publisher Full Text OpenURL

  47. Austin RH, Beeson K, Eisenstein L, Frauenfelder H, Gunsalus IC, Marshall VP: Dynamics of carbon monoxide binding by heme proteins.

    Science 1973, 181(99):541-543. PubMed Abstract | Publisher Full Text OpenURL

  48. Austin RH, Beeson KW, Eisenstein L, Frauenfelder H, Gunsalus IC: Dynamics of ligand binding to myoglobin.

    Biochemistry 1975, 14(24):5355-5373. PubMed Abstract OpenURL

  49. Jones G, Willett P, Glen RC, Leach AR, Taylor R: Development and validation of a genetic algorithm for flexible docking.

    J Mol Biol 1997, 267(3):727-748. PubMed Abstract | Publisher Full Text OpenURL

  50. Harata K, Muraki M: X-ray structure of turkey-egg lysozyme complex with tri-N-acetylchitotriose. Lack of binding ability at subsite A.

    Acta Crystallogr D Biol Crystallogr 1997, 53(Pt 6):650-657. PubMed Abstract | Publisher Full Text OpenURL

  51. Kullback S, Leibler RA: On information and sufficiency.

    Annals of Math Stats 1951, 22:79-86. OpenURL

  52. del Sol Mesa A, Pazos F, Valencia A: Automatic methods for predicting functionally important residues.

    J Mol Biol 2003, 326(4):1289-1302. PubMed Abstract | Publisher Full Text OpenURL

  53. Liu X, Zhang LM, Guan S, Zheng WM: Distances and classification of amino acids for different protein secondary structures.

    Phys Rev E Stat Nonlin Soft Matter Phys 2003, 67(5 Pt 1):51927. OpenURL

  54. Igarashi Y, Aoki KF, Mamitsuka H, Kuma K, Kanehisa M: The evolutionary repertoires of the eukaryotic-type ABC transporters in terms of the phylogeny of ATP-binding domains in eukaryotes and prokaryotes.

    Mol Biol Evol 2004, 21(11):2149-2160. PubMed Abstract | Publisher Full Text OpenURL

  55. Bhasi K, Zhang L, Brazeau D, Zhang A, Ramanathan M: Information-theoretic identification of predictive SNPs and supervised visualization of genome-wide association studies.

    Nucleic Acids Res 2006, 34(14):e101. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  56. Sterner B, Singh R, Berger B: Predicting and Annotating Catalytic Residues: An Information Theoretic Approach.

    J Comput Biol 2007, 14:1058-1073. PubMed Abstract | Publisher Full Text OpenURL

  57. Qian H: Relative entropy: free energy associated with equilibrium fluctuations and nonequilibrium deviations.

    Phys Rev E 2001, 63(4 Pt 1):42103. OpenURL

  58. Tirion MM: Large amplitude elastic motions in proteins from a single-parameter, atomic analysis.

    Physical Review Letters 1996, 77(9):1905-1908. PubMed Abstract | Publisher Full Text OpenURL

  59. Bahar I, Atilgan AR, Erman B: Direct evaluation of thermal fluctuations in proteins using a single-parameter harmonic potential.

    Fold Des 1997, 2(3):173-181. PubMed Abstract OpenURL

  60. Hinsen K: Analysis of domain motions by approximate normal mode calculations.

    Proteins 1998, 33(3):417-429. PubMed Abstract | Publisher Full Text OpenURL

  61. Atilgan AR, Durell SR, Jernigan RL, Demirel MC, Keskin O, Bahar I: Anisotropy of fluctuation dynamics of proteins with an elastic network model.

    Biophys J 2001, 80(1):505-515. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  62. Sanner MF, Olson AJ, Spehner JC: Reduced surface: an efficient way to compute molecular surfaces.

    Biopolymers 1996, 38(3):305-320. PubMed Abstract OpenURL

  63. Ankerst M, Breunig MM, Kriegel HP, Sander J: OPTICS: ordering points to identify the clustering structure. In Proceedings of the ACM SIGMON International Conference on Management of Data. Volume 28. Philadelphia, PA ; 1999::49-60. OpenURL

  64. Murzin AG, Brenner SE, Hubbard T, Chothia C: SCOP: a structural classification of proteins database for the investigation of sequences and structures.

    J Mol Biol 1995, 247(4):536-540. PubMed Abstract | Publisher Full Text OpenURL

  65. Varughese KI, Su Y, Cromwell D, Hasnain S, Xuong NH: Crystal structure of an actinidin-E-64 complex.

    Biochemistry 1992, 31(22):5172-5176. PubMed Abstract OpenURL

  66. Holt DA, Luengo JI, Yamashita DS, Oh HJ, Konialian AL, Yen HK, Rozamus LW, Brandt M, Bossard MJ, Levy MA, Eggleston DS, Liang J, Schultz LW, Stout TJ, Clardy J: Design, synthesis, and kinetic evaluation of high-affinity FKBP ligands and the X-ray crystal-structures of their complexes with FKBP12.

    J Am Chem Soc 1993, 115:9925-9938. OpenURL

  67. Weber PC, Ohlendorf DH, Wendoloski JJ, Salemme FR: Structural origins of high-affinity biotin binding to streptavidin.

    Science 1989, 243(4887):85-88. PubMed Abstract | Publisher Full Text OpenURL