Estimating affinities of calcium ions to proteins

Abstract

Ca²⁺-ions have a range of affinities to different proteins, depending on the various functions of these proteins. This makes the determination of Ca²⁺-protein affinities an interesting subject for functional studies. We have investigated the performance of two methods – Fold-X and AutoDock vina – in the prediction of Ca²⁺-protein affinities. Both methods, although based on different energy functions, showed virtually the same correlation with experimental affinities. Guided by insight from experiment, we further derived a simple linear model based on the solvent accessible surface of Ca²⁺ that had practically the same performance in terms of absolute errors as the more complex docking methods.

Keywords:

• metal ions
• binding
• free energy
• crystal structure
• solvent accessible surface

Introduction

Calcium ions (in the following termed Ca²⁺-ions or Ca²⁺) are important signalling agents that mediate a large number of intra- and extra-cellular processes,¹ for instance blood clotting, neurotransmission, or muscle contraction. Many of these processes involve proteins that bind Ca²⁺ more or less transiently, and accordingly, with a wide range of Ca²⁺-protein affinities.² In order to study the functional mechanisms of these proteins, it is desirable to determine these affinities. For some proteins these affinities have been determined experimentally,³ but to our experience these data are not available for most of the Ca²⁺-binding proteins. Given that the Protein Data Bank (PDB⁴) currently contains about five thousand structures of such proteins, it would be attractive for mechanistic studies to have a method at hand to quickly estimate Ca²⁺-protein affinities based on the structures of the corresponding complexes.

Theoretically, it should be possible to compute the Ca²⁺-protein affinity by free energy techniques based on physical models, eg, by pulling the ion out of its pocket with a series of umbrella potentials in molecular dynamics simulations and integrating over the potential of mean force,⁵ or by the Molecular Mechanics – Poisson-Boltzmann/Surface Area (MM-PB/SA) method.⁶ Despite their indisputable potential, these techniques have in the case of Ca²⁺-protein affinity so far not been demonstrated quantitative agreement with experiment, and, in addition, they are relatively costly in terms of computational resources. An alternative would be an estimation using empirical approaches. An early attempt towards a fast estimation of Ca²⁺-protein affinity was the work by Boguta et al⁷ who related secondary structure information with Ca²⁺-protein affinity. They found that for some proteins these relations could be used to classify their affinities, while for other proteins their scheme was less successful.⁸ More recently, Schymkowitz et al⁹ have published Fold-X, a method and empirical force field developed for, amongst other things, the fast prediction of the binding sites and affinities of metal ions, including Ca²⁺ and its rival Mg²⁺. Since Fold-X is a kind of docking method for special ligands, it would be interesting to compare the predictive performance of this method with that of a non-specialized state-of-the-art docking method in order to assess the advancement achieved by the special parametrization of protein-metal-ion interactions in Fold-X.

The importance of Ca²⁺-protein binding has over the years led to large body of experimental work from which a qualitative picture has emerged of the factors that govern Ca²⁺-protein affinity.^2,3 It is particularly notable that Ca²⁺- binding usually seems to be dominated by a gain in entropy, probably due to the release of water molecules from the solvation shell of the ion; in other words: the less ligands of the protein-bound Ca²⁺ are water molecules (and the more are functional groups of the protein), the tighter the binding. Such qualitative models may also be helpful for guiding the development of methods for affinity estimation.

In what follows we will address some of the points raised above, namely, we will compare the correlation of experimental affinities with affinities predicted by Fold-X and the state-of-the-art docking method AutoDock vina.¹⁰ The surprising finding of this comparison has prompted further study to possibly identify more simple computational models with the same power and speed in estimation of Ca²⁺-protein affinities. We will show that such models can be found and that they reflect knowledge attained by experimental work.

Materials and methods

We compared two docking methods for their ability to estimate Ca²⁺-protein affinities, Fold-X⁹ and AutoDock vina.¹⁰ Fold-X has been published as a method for the prediction of positions of metal-ions on proteins, and for the prediction of affinities between proteins and metal-ions. The energetic model underlying Fold-X has an ad hoc form with a number of parameters that have been fitted to experimental data.⁹ We obtained Fold-X versions 2.5.2 and 3.0b3 as executables for Linux from the respective server, including documentation. Despite considerable efforts we were not able to generate estimates of affinities of metal-ions to proteins using the commands described in the documentation of the software; we suspect that the option in those versions of the software is dysfunctional. Fortunately, the authors of Fold-X have offered as Table 4 in their supplementary material to Ref⁹ a list of 48 Ca²⁺-binding pockets in 19 X-ray structures, mostly with experimentally determined affinities, and affinities predicted with an earlier version of Fold-X. Hence, we took these data (“Fold-X dataset”) as basis for the comparison, specifically the columns “experimental energy” and “predicted energy” of Table 4 in Ref.⁹ For five of the binding pockets, two experimental energies were given; since the differences between the first and second energies were relatively small, only the first value was considered in each case. For ease of comparison with Ref ⁹ we give all affinities in units of kcal/mol (1 kcal/mol = 4.1868 kJ/mol).

The 19 X-ray structures of the Fold-X dataset were retrieved from the Protein Data-bank (PDB⁴) for comparative analysis with AutoDock vina.¹⁰ AutoDock vina version 1.0.3 for Linux was downloaded as executable from the website of its authors.

Usually, some information is missing from X-ray structures that is needed for energy calculations. Most importantly, this is the case for hydrogen positions, including also hydrogen bond networks. Related to this is the possibility to optimize X-ray structures by flipping carbonyl-oxygens and –NH groups (both groups have similar electron densities). Finally, sometimes X-ray structures contain atomic overlaps that can be removed relatively easily. To see whether by considering these effects, affinity predictions can be improved, we used three different protocols with AutoDock vina.

In the first protocol, PDB files were prepared with the AutoDockTool suite of AutoDock 4¹¹ by removing water molecules, non-standard residues, and alternate positions of residues. Then polar hydrogens were added to the protein in standard orientation without rotational optimization, and Gasteiger charges¹² were computed for protein atoms by AutoDock. Finally, all Ca²⁺-ions in the protein were redocked with AutoDock vina using default parameters, except for the “search space”, ie, the volume in which the optimal docking position is searched for. Since we were interested in affinities at the crystallographically determined positions of the Ca²⁺-ions, and not in finding optimal positions, we restricted the search space to the minimum allowed by AutoDock vina around the crystallographically determined positions, namely a cube of 1 Å (0.1 nm) length in x-, y-, and z-directions around the crystallographic positions of the Ca²⁺-ions. This protocol was applied to all Ca²⁺-ions in the Fold-X dataset.

In the second protocol we introduced in the preparation of the protein structures a further step in which the “reduce” method¹³ was used to optimize positions of hydrogen atoms around crystallographically determined heavy-atom positions, including also potential flips of amide-groups in side chains of asparagine and glutamine. Otherwise this second protocol had the same elements as the first, including also the restrained docking of Ca²⁺-ions with AutoDock vina.

The third protocol was similar to the first one, but additionally optimized hydrogen positions with pdb2pqr,¹⁴ including debumping to avoid steric clashes and optimization of the hydrogen bonding network.

For computations of solvent accessible surfaces (SAS) we used MSMS version 2.5.7¹⁵ with standard atomic radii for protein atoms and 1 Å for Ca²⁺. All statistical analyses were carried out with R.¹⁶

Results and discussion

Comparison of Fold-X and AutoDock vina

We first analyzed the correlation of Ca²⁺-affinity predictions by Fold-X and AutoDock vina with experimental data in the Fold-X dataset. To this end, linear models were fitted using the least-squares algorithm. The Fold-X predictions had a value of Pearson correlation coefficient r of 0.67 with a least-square fitted straight line ΔG_exp = 3.4221 kcal/mol + 0.5854 ΔG_pred (Figure 1). The predictions by AutoDock vina (using the first protocol described in “Materials and methods”) had r = 0.71 with a best-fit line ΔG_exp = −2.393 kcal/mol + 7.309 ΔG_pred (Fig. 2). Both above correlation coefficients are different from zero with high significance according to t-tests (p ≈ 10⁻⁶). Although the range of predictions by Fold- X has a much better overlap with the range of experimental affinities, the linear correlation as given by r is slightly worse than that obtained with the AutoDock vina predictions. However, closer inspection of the data shows that the ranking of r-values is mainly based on a single value, the outlier in the lower left corner in Fig. 2. This outlier is given in the Fold-X dataset with a negative experimental affinity of −1.2 kcal/mol, though without reference to an experimental source of that value. The corresponding structure in the Fold-X dataset is that of a calmodulin of Paramecium tetraurelia¹⁷ (PDB entry 1exr), and the conspicuous affinity value probably refers to an unusual fifth Ca²⁺-ion bound to a pocket that, according to the crystallographers, probably had been created by crystal contacts and thus may be without functional relevance. In the Fold-X dataset there is no prediction given for this pocket, while AutoDock vina produces the mentioned outlier. We can interpret this complex as representative of an extremely weakly bound Ca²⁺, and if we replace the experimental value of −1.2 kcal/mol by 0 kcal/mol (ie, zero affinity) the correlation coefficient r of AutoDock vina and experiment is practically unchanged at 0.71. If we omit this experimental value altogether, the value of r drops to 0.63, which is somewhat lower than r = 0.67 of Fold-X with experiment. If we consider the 95% confidence intervals based on a t-test for r we obtain [0.46, 0.81] for Fold-X vs experiment and [0.41, 0.78] for AutoDock vina vs experiment. Thus, the correlation coefficients between Fold-X predictions and experiment, and AutoDock vina predictions and experiment are virtually equal. This is astonishing since Fold-X had been at least partially calibrated with the same dataset, as mentioned in Ref,⁸ while AutoDock vina has probably not been specifically developed to solve the problem for which we have employed it here.

Figure 1 Correlation of Fold-X predictions and experiment. The straight line is a least-square fit between experimentally determined Ca²⁺-protein affinities (ΔG_exp) and affinities predicted with Fold-X (ΔG_pred). Pearson correlation coefficient is r = 0.67. Data.⁸

Figure 2 Correlation of AutoDock vina predictions and experiment. The dashed straight line is the least-square fit between experimentally determined Ca²⁺-protein affinities (ΔG_exp) and affinities predicted with AutoDock vina (ΔG_pred). Pearson correlation coefficient is r = 0.71. If the outlier in the lower-left is dropped, r decreases to 0.63 (solid line). Experimental affinities.⁸

The two other protocols used in conjunction with AutoDock vina (see “Materials and methods”) had little effect and did not improve the correlation with experiment (r = 0.71 and r = 0.69 for the second and third protocol, respectively). This may be due to the fact that these protocols mainly affect hydrogen positions, while Ca²⁺-binding pockets usually are dominated by anionic groups with few protons.

A simple model for estimating Ca²⁺–protein affinities

The fact that Fold-X and AutoDock vina did show the same correlation with experiment could be a consequence of the two underlying energy models capturing the same dominating cause of Ca²⁺-protein affinity. If this hypothesis is true, we should find a strongly decreased correlation of the two models with experiment after removing the contribution of that dominating cause out of the data.

A candidate for such a dominating effect mentioned in the introduction is the entropy gain due to water molecules that are released from the first hydration shell of Ca²⁺ on binding of the ion to the protein. In other words, the more water molecules are still attached to the protein-bound Ca²⁺, the lower the entropy gain and thus, the lower the affinity. This argument suggests an avenue to a computational test of the above hypothesis of a dominating factor: if we assume that the number of water molecules attached to the protein-bound Ca²⁺ is proportional to the solvent accessible surface (SAS) of the ion, we should expect a negative linear correlation of experimental free energy of binding ΔG_exp and SAS. The part of ΔG_exp not explained by the correlation with SAS is then contained in the residuals e_i = αSAS_i + β −ΔG_exp,i, with α and β slope and intercept, respectively, of the least-squares fitted linear model, and SAS_i and ΔG_exp,i the SAS and experimental affinity, respectively, of the ith Ca²⁺ in the Fold-X dataset. If the above hypothesis of a dominating factor is true, there should be a much smaller correlation between the predicted ΔG_FoldX,i (and ΔG_vina,i) and e_i as compared to the correlation between ΔG_FoldX,i (and ΔG_vina,i) and ΔG_exp,i. In the following we carry out this partial correlation analysis.

As solvent probe radius we first assumed 1.4 Å, a value that is frequently used to model a molecular “water-sphere”. With this sphere the distribution of SAS values was strongly skewed with a peak at the lowest SAS values. In fact, in thirteen of the binding pockets in the Fold-X dataset, Ca²⁺ was not accessible at all (SAS = 0 Ų). An Anderson-Darling test rejected with high significance that the SAS values are normally distributed. Therefore, the correlation with ΔG_exp was not tested with Pearson correlation coefficient r but with Spearman rank correlation ρ. We found that ρ = −0.52 was significantly different from zero (significance level 0.05, P = 2 · 10⁻⁴) (see Table 1). This correlation dropped only slightly when the outlier discussed above was omitted.

Table 1 Spearman rank correlation ρ of four models with experimental affinities ΔG_exp and residuals e

Modela	ρ(ΔGexp, Model)b	ρ(eΔGexp, SAS1.4, Model)c	ρ(eΔGexp, SAS0.5, Model)d
SAS1.4	−0.52 (2 · 10^−4)
SAS0.5	−0.52 (2 · 10^−4)
ΔGFoldX	0.66 (2 · 10^−6)	0.52 (4 · 10^−4)	0.39 (0.01)
ΔGvina	0.64 (3 · 10^−6)	0.50 (4 · 10^−4)	0.34 (0.02)

a Notes: four linear least squares fit models with input variables SAS with probe radius 1.4 Å, SAS with probe radius of 0.5 Å (see Eq. (1)$\mathrm{\Delta }{G}_{exp}\approx -0.63·SA{S}_{0.5}+11.95$(1) ), affinities computed with Fold-X, affinities computed with AutoDock vina;

b Spearman rank correlation ρ of experimental affinities ΔG_exp with predictions of fitted models;

c ρ of residuals of model based on fit with SAS^1.4 with Fold-X and AutoDock vina, respectively;

d ρ of residuals of model based on fit with SAS^0.5 with Fold-X and AutoDock vina. Numbers in parentheses are p-values for the null-hypothesis that ρ = 0.

The correlation of ΔG_FoldX with the residuals e_i of the least-squares fitted linear model ΔG_exp(SAS) was somewhat lower (ρ = 0.52) than the correlation of ΔG_FoldX with ΔG_exp, but remained highly significant (P = 4 · 10⁻⁴). The same was true for ΔG_vina (ρ = 0.50, P = 4 · 10⁻⁴). In view of our hypothesis this means that there is an effect on the affinity that can be formulated in terms of SAS, but it may not be dominating affinity.

Our argument has so far neglected the fact that the crystal structures are results of an averaging process. A protein at ambient temperatures explores many conformations, so that the solvent accessibility of Ca²⁺-ions computed for the crystal structure may not reflect its true accessibility. Since Ca²⁺-ions in proteins are often surrounded by a tightly packed first co- ordination shell and thus have minimum solvent accessibility, we expect that conformational flexibility could perturb that packing and thus lead on average to a higher solvent accessibility. An approach that takes mobility into account could be to simulate the molecular dynamics and compute the accessibility as thermodynamic average. As this is computationally expensive, and we were more interested in a fast approximation, we tried to find a faster alternative that works in a similar direction. In a sense, the higher accessibility due to the protein flexibility can be mimicked by using a probe with a smaller radius. We therefore carried out the partial correlation analysis described above with a series of smaller probe radii between 1.4 Å and 0.3 Å. The highest correlations of SAS with ΔG_exp were obtained with 0.4 Å and 0.5 Å. As the numerically safe minimum probe radius in MSMS (see “Materials and methods”) is 0.5 Å, we completed our analysis with this value (see Table 1 and Figure 3).

Figure 3 Correlation of SAS (probe radius 0.5 Å) and experiment. The straight line is the least-squares fit between experimentally determined Ca²⁺-protein affinities (ΔG_exp) and the solvent accessible surfaces (SAS) of the Ca²⁺-ions. Spearman rank correlation coefficient ρ is −0.53 (P = 2 · 10⁻⁴ for null-hypothesis ρ = 0). Experimental affinities.⁸

While the correlation ρ of SAS with ΔG_exp changed only marginally from −0.52 to −0.53 and the p-value remained constant, the correlation of the residuals with ΔG_FoldX and ΔG_vina, respectively, dropped more strongly to 0.39 and 0.34 with p-values of 0.01 and 0.02 indicating no longer highly significant correlation. Omission of the discussed outlier does not change the picture; conversely, one could argue that a fully solvent exposed Ca²⁺ with ΔG = 0 (which is approximately the case for the outlier) should be included in the data to represent the limiting case of no binding. Overall we can conclude that SAS with probe radius of 0.5 Å (SAS_0.5) indeed models a dominating effect on affinity.

For a simple model to predict ΔG from SAS we did a least-squares fit of SAS_0.5 and ΔG_exp and found

$\mathrm{\Delta }{G}_{exp}\approx -0.63·SA{S}_{0.5}+11.95$(1)

with SAS_0.5 in Ų and ΔG_exp in kcal/mol. Following the above argument we replaced in the fit the outlier by the point given by theoretical values for an unbound Ca²⁺, ie, SAS_0.5 = 4π(r_ca−ion + r_probe)² = 4π1.5² = 18.85 and ΔG = 0, and enforced inclusion of this point (18.85, 0) in the linear model. The mean of the absolute errors of this model on the Fold-X dataset was 1.9 kcal/mol, which is not much larger than the mean absolute error of 1.8 kcal/mol of Fold-X itself against ΔG_exp.

We assessed the robustness of the model of Eq. (1)$\mathrm{\Delta }{G}_{exp}\approx -0.63·SA{S}_{0.5}+11.95$(1) in a leave-one-out test: Each of the ΔG_exp values was left out of the fitting procedure, and then a model was derived from the other values. The ΔG_exp of the Ca²⁺-ion left out was then predicted by applying the new model to the SAS_0.5 value of the left-out Ca²⁺-ion. This was iterated over all ΔG values. The resulting mean of absolute errors was 2.0 kcal/mol.

Thus, we have achieved our goal of a simple and fast computational procedure that allows an estimation of Ca²⁺- protein affinities based on the structure of the complexes. Judged from the numerical experiments described above, the accuracy of the method should be high enough to classify Ca²⁺-binding pockets into weakly or strongly binding. The accuracy is limited by several factors, of which we mention two: First, as pointed out in Ref ⁸ the experimental data on which the model has been based may in part not satisfy modern standards. Second, the simple model of Eq. (1)$\mathrm{\Delta }{G}_{exp}\approx -0.63·SA{S}_{0.5}+11.95$(1) does completely neglect that the binding of calcium is often accompanied by global re-arrangements of protein conformation that also affect the free energy of binding.²

Finally, we can speculate how Fold-X and AutoDock vina with their different energy functions could nevertheless capture the effect expressed in terms of SAS that we mainly attribute to the entropy gain due to release of water molecules bound to the solvated Ca²⁺. Neither the energy function of Fold-X nor that of AutoDock vina contains a term that explicitly takes into account this physical effect. However, both Fold-X and AutoDock vina evaluate the Ca²⁺-protein affinity essentially by estimating interactions of Ca²⁺ with the atoms of the protein lining the binding pocket. According to our experience these pockets are dominated by anionic groups and groups with negative partial charges. The direct interaction of Ca²⁺ with such negative groups is in fact taken into account by Fold-X and AutoDock vina, and this may be the cause of the apparent correlation with SAS: the more negative groups are around a Ca²⁺, the lower the predicted affinity due to direct interaction, but also the lower the SAS of that Ca²⁺, because each of the neighboring groups will supplant water molecules and lead to their release.

Acknowledgments/disclosures

The authors report no conflicts of interest in this work. Funding by BMBF grant number 01EZ0933 is gratefully acknowledged.