|Home | About | Journals | Submit | Contact Us | Français|
We present a novel approach for computing biomolecular interaction binding affinities based on a simple path integral solution of the Fokker-Planck equation. Computing the free energy of protein-ligand interactions can expedite structure-based drug design. Traditionally, the problem is seen through the lens of statistical thermodynamics. The computations can become, however, prohibitively long for the change in the free energy upon binding to be determined accurately. In this work we present a different approach based on a stochastic kinetic formalism. Inspired by Feynman's path integral formulation, we extend the theory to classical interacting systems. The ligand is modeled as a Brownian particle subjected to the effective non-bonding interaction potential of the receptor. This allows the calculation of the relative binding affinities of interacting biomolecules in water to be computed as a function of the ligand's diffusivity and the curvature of the potential surface in the vicinity of the binding minimum. The calculation is thus exceedingly rapid. In test cases, the correlation coefficient between actual and computed free energies is >0.93 for accurate data-sets.
The accumulation of atomic-resolution information of protein-ligand interactions available through public databases of X-ray crystallographic and NMR structures1 has paved the way for the development of several theoretical methods for binding free energy prediction2 based on atomistic representations of the protein-ligand complex. Typically, the theoretical foundation is based on statistical thermodynamics. For an interaction between a protein and a ligand, E and I, interacting in solution to form a complex E*I we can write the biomolecular reaction
At equilibrium we can write for the standard free energy of association
where μoE*I, μoE, μoI are the standard chemical potentials of the complex and the individual species respectively, R is the ideal gas constant, T the temperature and 1/Ki =kon/koff is the binding equilibrium constant.
These macroscopic thermodynamic properties connect to microscopic properties determined by atomistic computer simulations through the classical statistical thermodynamics relationship
where QE*I , QE , QI , QW are the molecular, canonical ensemble partition functions of the complex, the individual species and the solvent respectively. In principle, the partition functions enumerate all the possible microscopic states of the molecules. In practice, the direct calculation of the partition function for as complex a system as a solvated protein is theoretically and computationally unfeasible, because of the configurational integral. For one of the interacting species, say inhibitor I, this integral is
where and are the are the dimensions of the configurational space available to the N molecules of species I and the M molecules of solvent W. is the potential energy of interaction between I and W.
There have been ingenious efforts to approximate the calculation, decomposing the free energy into numerous components that can be tractably computed3-14. Nonetheless, difficulties persist, mainly because the error in the calculations of free energy components is larger than the actual, absolute value of the binding strength. Therefore, although theoretically on firm ground, the statistical thermodynamics approach to use atomistic models of biological molecules to predict the free energy of binding can be fruitful only with the injection of empirically derived corrections or severe simplifying assumptions. We propose the following alternative approach.
Considering the phase orbit of a classical mechanical system in phase space and integrating out the solvent degrees of freedom, the path will resemble the motion of a Brownian particle described by a Langevin equation:
where W(t) is Gaussian noise with average <W(t)>=0 and correlation W(t)W(t') = 2Dδ(t–t'), D is the diffusion of the system in phase space, kB and T are the Boltzmann factor and temperature respectively.
An equivalent Fokker-Planck equation is
The effective potential is
The conditional probability solution for remaining in the same state space position for small time intervals is
Here t is very short, of the order of the average duration of solvent collisions, and n is the number of degrees of freedom.
In a lucid analysis of the Langevin equation, Grooth 19 surmises that the friction coefficient γ = kB T/D = 2 ms f, where ms is the mass of the solvent molecules and f is the number of collisions per second, giving the average time per collision as 2 ms D/kB T. Therefore,
If we consider the reversible reaction between the two areas of phase space representing the bound and unbound states of the ligand given by
where [Ibound]=[E*I], the transitional probability to stay in the bound conformation (denoted by xB) is20
Equations 14 and 15 provide an elegant means for computing the free energy of biomolecular interactions given the bound state structure. Eq 15 has only two variables – the diffusion coefficient of the ligand, D, and the second derivatives of the potential for the bound complex. The latter can be computed quickly given the bound state, as described later.
Diffusion coefficients for organic molecules can also be estimated to high accuracy by semiempirical derivatives of the Stokes-Einstein relation which typically give scaling functions D~V−0.6 where V is the molar volume of the solute.21
Even though, in theory, the right hand side (RHS) of eq 15 can be exactly determined, experimental errors would give rise to empirical coefficients for the two terms in RHS. In practice, therefore, these variables in eq 15 must be trained with a set of protein-ligand complexes with known binding free energies.
Here we present a systematic study of the application of eq 15 towards the prediction of binding affinities of three different enzyme-inhibitor systems – bovine trypsin, β-secretase and aldose reductase. We restrict our study to systems for which experimental inhibitor affinity measures (Ki, IC50) as well as X-ray crystal structures for the bound complexes are available. We find that eq 15 is not only a convenient and accurate means for predicting binding affinities, but also provides useful insights on the effect of the oft-ignored assay conditions on measured binding affinities.
Diffusivities for organic molecules at infinite dilutions can be accurately determined by semiempirical derivatives of the Stokes-Einstein relationship of which the following equation by Wilke-Chang is perhaps the most popular,21
where χ indicates the degree of solvent association, Mw is the molar mass of the solvent (18.015 g mol−1), η is solvent viscosity (mPa s) at temperature T (K) and V is the Le Bas molar volume22 of the solute at its normal boiling point. It has been shown previously23 that the van der Waals volume of the molecule in Å3 (Vvdw) estimated by molecular modeling software can approximate Le Bas volume. Vvdw was therefore calculated in MOE24 using a grid approximation with a spacing of 0.75 Å.
The first and second terms on the RHS of eq 15 (denoted as TermA and TermB, respectively) were computed using DWC for diffusivity D. These computed terms were used to train the linear regression fit, giving an equation of the form:
where A and B are empirical constants and I is the intercept. For the trypsin dataset ln[Ki′] = ln[1 + 4Ki/[ETotal]], where the binding constant, Km, and the total enzyme concentration, [ETotal] were known. The binding affinity data were approximated as ln[1+Ki] or ln[IC50] for the other two sets.
Since χ is a dimensionless empirical parameter, different values have been proposed for χ in the literature for different types of molecules – 2.6 by Wilke and Chang21, 2.9 for organic electrolytes25, 2.26 for nonelectrolytes26 and 1.61 for aromatics27. In this study we used χ =2.6 and introduced a new empirical parameter , such that D = *DWC, to correct for the use of Vvdw instead of Le Bas volume, as well as for the presence of other solutes in the assay buffer. A second normalizing term, λ, was introduced to eliminate the intercept.
The modified eq 15 in three dimensions is therefore,
Note that parameters and λ are not necessary to obtain a useful equation for predicting activities. However, determination of helps to compare the effective diffusivities of the ligands between datasets, while parameter λ nondimensionalizes the probability of eq 10.
For the first derivative of the potential to be zero in Veff (eq 9), the solvated enzyme/ligand complex has to be proximal to its lowest energy conformation. Protein/ligand complexes were therefore minimized prior to computation of the second derivatives. The ligand/receptor complexes were minimized in MOE24 to a RMS gradient of 0.001 using the MMFF94 force field28. The RMS gradient is the product of norm of the potential gradient and the square root of the number of unfixed atoms. Non-bonded interactions were evaluated without any cut-offs. During minimization, the solvent was implicitly modeled through the Generalized-Born model as implemented in MOE. Only residues of the receptor within 5 Å of the ligand were selected for the minimization since we were only interested in minor side-chain modifications in the vicinity of the ligand.
If we treat the inhibitor as a rigid body at the potential minumum, 2U(x)/x2 of eq 15 is the trace of the resultant three dimensional Hessian matrix (k = kxx + kyy + kzz). The components of the Hessian matrix were evaluated from a finite difference of the forces about the potential minimum as follows:
with Δi=Δj=1×10−6 Å, where i and j refer to the three Cartesian dimensions x,y and z, Fi is the force at position i. Forces were evaluated using the Potential function in MOE. The trace of the Hessian is therefore k = kxx + kyy + kzz.
Statistical analysis and linear regression was carried out using JMP29. The leave-one-out cross validation coefficient, R2cv, and the root mean square error of cross validation (scv) were computed using MATLAB30.
TermA and TermB of eq 17 were computed for the entire dataset (Table 1). The linear regression fit for the entire dataset is presented in figure 1. While the R2 of the fit is not very significant (0.51), a close examination of figure 1 presents some remarkable trends. Significant clustering is observed for points belonging to the same protein, with either systematic over- or under prediction within data-points of those sub-sets. We assumed that these differences were due to either the specific details of the assay buffers, or the chemical characteristics of ligands themselves, both of which would affect the ligands’ diffusivities. The Wilke-Chang equation (eq 16) evaluates diffusivities for molecules at infinite dilution, and systematic deviations from the calculated diffusivities are expected that are dependent on the assay buffers. We therefore decided that each subset had to be trained separately for better fit. Also intriguing was the presence of two clusters in the β-secretase data-set, one of which was significantly underpredicted (figure 1) and belonged entirely to data from the work of the same group presented in two different papers31,32. A comparison of the assays between this subset and the rest of the β-secretase complexes revealed that these assays had 10% dimethylsulfoxide (DMSO) in the assay buffer. DMSO is a highly viscous solvent and is bound to have significant effect on the diffusivity of the ligand molecules. Aminabhavi and Gopalakrishna33 report a viscoscity (mPa.s) increase to 1.787 from 0.891 as the mole fraction of water changes from 1 to 0.9. A 10% DMSO buffer, therefore increases the buffer viscoscity by ~30%, sufficient to cause the relative shift observed in figure 1. It is also quite likely that other effects of DMSO, such as altered solubility and hydrophobicity play an important role towards the activity as well. These data points were therefore separated into two, separate subsets for training the linear regression fit.
The statistics of the linear regression fit for the four different subsets belonging to bovine trypsin, β-secretase (with and without DMSO) and aldose reductase are presented in Table 2.
Experimental evaluations of the inhibitor affinity measurements are very sensitive to experimental conditions with typical coefficients of variation of ~20% in the reported analytical measurements. The best scenario therefore involves data obtained by the same laboratory and the use of the same assay. The trypsin data-set, therefore, was of very high quality, since all of the data points used in that set were obtained by the same team. Katz et al 34,35 conducted a high-throughput study involving the determination of over 300 high resolution crystal structures of trypsin bound with small molecules inhibitors (2-(2-phenol-indoles and 2-(2-phenol)-benzamidazoles) over a wide range of pH. Scanning through the available structures led to a data-set of 12 inhibitors with their inhibition constants (Kis). This reduced data-set was due to the fact that all the inhibitors studied here exhibit a parabolic pH dependence with a minimum in Ki and therefore only structures determined at pHs close to Ki(min) were used, to ensure that the selected crystal structures captured the protonation state of the inhibitors without any ambiguity. At pH's away from Ki(min) a mixture of protonation states is expected for the inhibitors thereby complicating the picture. This high resolution data-set involving crystal structures for all the inhibitors ensures that any errors due to in silico mutations and the subsequent molecular mechanics minimization step are avoided. The linear regression fit for this data-set is shown in figure 2.
While there are no outliers in this data-set, the hydroxyl groups of structures 1o36 and 1o3g had to be protonated for a better fit. This was not unreasonable considering that the pKa's of the hydroxyl groups for the free molecules are predicted to be close to the pH of the assays. The activities are predicted to very high accuracy (R2 = 0.93, σ= 0.62, Table 2) in the prediction of ln[Ki′]. A σ= 0.62 implies that the Kis are predicted to the correct order of magnitude, such resolution is rarely seen in computational free energy prediction.
The β-secretase data-set is derived from the work of several independent laboratories each with differing assay conditions, substrates, substrate- and enzyme concentrations and measures of affinity (Ki and IC50s) (Table 1). The LHS of eq 18 has the term 1+Ki/[E] (= IC50/[E*I]) (see supporting information). The concentration of the enzyme used in these assays was not available for most of the data-points and was therefore approximated by ln[1+Ki] or ln[IC50] for consistency. This approximation is expected to increase the noise in the dataset to above that of experimental error. The linear regression fit for the two subsets of β-secretase are presented in figures 3 and and4.4. 1tqf and 2f3e were included in data-set 1 based on the regression fit. It was not possible for us to determine if that was reasonable, as details of their assays36,37 are not readily available. Congreve et al point that the inhibition data for the inhibitor of 2oht should be treated with caution as they exhibited a high hill slope31. Interestingly, 2oht is identified as an outlier in the linear regression, and in fact, its inclusion in the fit significantly reduces the accuracy of the fit (R2 = 0.95, σ = 0.5488, F = 63.43), to below the accuracy obtained by including TermB alone in the fit (R2 = 0.95, σ = 0.5213, F = 140.3626). This underscores the importance of accuracy in the data-set in training the regression fit. For this dataset, TermB dominates in the fit, with TermA offering further resolution to the prediction of activities. In fact the agreement between predicted and actual IC50's is remarkable, a σ = 0.2913 implying that the predicted IC50s are within ~35% of experimental Ki's, with experimental error itself typically of the order of 20%. The domination of TermB in the prediction formula is due to the suppressed values of the diffusivities obtained (Table 2), which appears to be the effect of DMSO in the assay. Since, in TermA the trace of the Hessian matrix is multiplied by D2, its value is very sensitive to the diffusion coefficient. Significantly lower values of D therefore increase the significance of TermB towards the inhibitory activity.
The second data-set has five outliers, 1ym1, 1ym2, 2iso, 2irz and 2iqg. A sixth outlier, the ligand of 2oah, was assayed without separation of its stereoisomers, therefore justifying its exclusion from the fit. Due to the significant reduction in the accuracy of prediction with the inclusion of 2oht in the previous data-set, we chose to exclude 1ym1, 1ym2, 2iqg, 2iso and 2irz from the fit, even though we are aware that the exclusion of these points artificially improves the statistics of the fit. Supporting our choice is the fact that the four of the outliers are from the work of two different papers 38,39, and three of these – 2iqg, 2iso, and 2irz (and 2ph6) use as a substrate a fusion protein containing maltose binding protein at the amino terminal end and 125 amino acids of the amyloid precursor protein amino terminal end. All the other assays use significantly smaller peptide substrate in the assay. The use of a constant χ in the Wilke-Chang equation for the empirical fit could lead to errors in the estimation of D, since its value suggested value varies from 1.61 to 2.9 depending on the chemical nature of the molecule25,27. This is a flaw in our analysis since the specific chemical characteristics of the molecule were not considered in the empirical determination of D. The ligands of 1ym1 and 1ym2, for instance, have very little aromatic character to them while that is not true of the other ligands. Also, conditions specific to those assays may also be responsible for the observed deviation. Even so, the deviation in the predicted ln[Ki′]s for these outliers is within a factor of 2.5, which implies that the predicted activity is correct to an order of magnitude. In general, since the ligands in this study are significantly more complex than those used to drive the empirical fit in the Wilke-Chang equation, extrapolation to more complex molecules may increase the error.
The aldose reductase data-set 40-53 is expected to be most noisy, since the data spans over a decade of work and is obtained from the work of several different laboratories. In fact an examination of the assays shows that the affinity data was obtained from the use of a wide range of substrates (glucose, xylitol, xylose, glyceraldehydes, and benzyl alcohol). The LHS of eq 18 involves [E] whose value depends on [ETotal] and the ratio of [S]/Km, where [S ] is the substrate concentration. Assays using different substrates using varying ratios of [S]/Km are therefore likely to introduce further noise into the data-set beyond that due to experimental error. In fact, for two of the inhibitors, different values of affinity separated by an order of magnitude have been reported by different groups; 44nM 44 and 40054 nM for zenarestat (1iei) and, 354nM and 6054 nM for zopolrestat (1mar). We used the more recent estimates in our analysis for consistency.
Since eq 18 predicts the affinities to very high precision, noisy data is expected to lead to inaccurate empirical prediction formulas for the activity. This is observed in the fit for this data-set (figure 5), where the contribution of the TermA is not captured at all. This data-set has three outliers - 1pwl, 1z8a, 1z89. While we can not explain the inaccuracy in the prediction of 1pwl, an examination of the assay buffers of the ligands of 1z89 and 1z8a, not surprisingly, found 20% DMSO in the assay buffer. Here, (σ = 1.22) activities are predicted to within an order of magnitude, which is close to the expected level of noise in this data-set. The widely varying chemical composition of the ligands in this data-set will also increase errors in the evaluation of the diffusivities. A high sensitivity to the diffusivity is expected in the prediction of the activities due to the nature of eq 18, with the sensitivity being much higher for TermA than TermB, since a D2 features in TermA while TermB involves ln (1/D3).
In general it appears that eq 18 can predict free energies almost to the precision of the available data for a linear regression fit. Of great interest is the effect of the magnitude of the diffusion coefficient on the observed inhibitor affinities. The systematic error observed in the residuals in figure 1 (Table 1) appears to be due to the different correction factor () that was required for different enzyme assays. The correction factor, , was determined by equalizing the coefficients of TermA and TermB in the linear regression fit (Table 2). The value of explains the strong dependence on TermB in data-set 2a. This appears to be due to the lower diffusivities of these molecules in comparison to the trypsin dataset. Larger diffusivities of the ligands in the trypsin data-set reduce the accuracy of prediction with TermB alone in the trypsin dataset (R2 = 0.39, σ=1.55, F=5.84). The decrease in the diffusivities of data-set 2a is most likely due to the increased viscosity of the solution due to addition of DMSO, since viscosity appears in the denominator of Wilke-Chang equation. The sensitivity of the activity of the molecules to their diffusion coefficients highlights the importance of the conditions of the assay towards the observed activity.
An interesting implication of the strong dependence on TermB for subset 2a is that the energetics of the binding are of little consequence for predicting activity for such systems. The activity almost entirely depends on the diffusion coefficient. The Wilke-Chang equation implies that activity will be strongly correlated with molecular volume. Since molecular masses and volumes are strongly correlated (R2 = 0.96 for the entire data-set), the molecular mass would also be highly correlated with activity. This interesting fact has been observed by Kim and Skolnick,55 where a comparison of scoring functions of different docking algorithms led to the observation that molecular mass was as good a predictor of activity as the scoring functions themselves. There, the authors argue that this observation may be due to the fact that most rational drug design involves improvement of efficacy of drug leads through the addition of favorable functional moeties, thus leading to an increase in mass. Equation 18 implies that the reason for the observed correlation is more fundamental since the size of the molecule directly affects its diffusivity. In fact we find that the use of the logarithm of the molecular weight (ln(MW)) and k/MW (where k is the trace of the hessian matrix) resulted in empirical correlations with similar accuracies to those achieved by eq 18.
Apart from the composition of the assay buffer, other factors such as temperature also strongly influence D. For instance, from the Wilke-Chang equation, the ratio of diffusivities of a molecule at infinite dilution at temperatures 310 and 298 would be 1.34. Thus, for a low molecular weight ligand, eq 18 would predict significantly different activities at different temperatures.
One common assumption made to simplify the configurational integral of eq 4 is that it can be approximated by “end-point” calculations that limit the computation to the solvated bound/unbound ligand/receptor alone. Equation 18 implies that such an approximation will capture only a part of the picture, neglecting the contribution of the molecules’ diffusivity towards binding affinity.
Despite the success of eq 18 with the trypsin and β-secretase data-sets and the fairly low computational cost of its evaluation, several problems exist with its usage on a large scale for predicting activities. The first comes from the limitations of the empirical nature of the Wilke-Chang equation for predicting diffusivities. It is not surprising therefore, that the highest accuracy in the prediction of activities was achieved for data-sets 1 and 2a where the ligands were of similar chemical nature. For a data-set with molecules of widely differing properties, such as data-set 3, the use of a single empirical coefficient will most likely lead to significant errors in the prediction of diffusion coefficients. Besides, the success of the Wilke-Chang equation has been demonstrated for a limited number of small organic molecules in solution and it is not clear if its extrapolation to larger and chemically more diverse molecules typical of drug candidates would lead to accurate predictions of diffusion coefficients. It appears therefore that independent methods for determining diffusivities may have to supplement drug design experiments.
A second difficulty of the method outlined here is with the determination of the trace of the Hessian matrix. The trace is very sensitive to the final minimized complex as well as to the minimization protocol used. The correlation coefficient, R2, of the linear regression fit for trypsin data-set decreases significantly when cut-offs (8 Å) were used for evaluating the non-bonded interactions. A Generalized-Born solvation model was necessary during minimization to simulate the effect of the solvent. This is expected since the minimum of Veff is for the solvated complex. Explicit solvation could likely improve the fit even more, but the computational costs would quickly become prohibitive. Careful attention must be paid to the protonation states of the ligands and any titrable protein side-chains. For instance, 1o36 and 1o3g fit only after protonation of a hydroxyl group.
Also interesting is the relationship that results for Veff, where the significant factor for unbinding is not the interaction potential as intuition suggests, but its second derivative. In other words, the unbinding free energy depends not on the depth of the potential well that the ligand finds itself in, but rather on the curvature of this potential in the vicinity of energy minimum. This is a consequence of the over-damping effect of the solvent that appears as an assumption in the Fokker-Plank equation.
Our greatest challenge in this work, was to find data-sets of sufficiently large size (for statistical significance) that met all of our criteria – namely the availability of affinity data together with crystal structures for every complex, as well as the absence of significant conformational rearrangements in the active-site associated with binding. The absence of metal atoms and significant interactions with crystal waters was also considered essential as we were interested in the simplest possible data-sets to test the theory. We believe that it is due to this choice that we obtained reasonably good fits despite significant simplification of the physics. It is quite likely that other rate-limiting steps such as secondary binding events prior to unbinding or conformational rearrangements in the active-site during unbinding would significantly complicate the picture and perhaps demand a less trivial solution.
Two other methods for affinity prediction that were derived involving assumptions, albeit very different from ours, that lead to a simplification of the underlying physics and hence a reduction of the computational expense can be mentioned here. The first of these is the Linear Interaction Energy (LIE) method developed by Aqvist and coworkers56. The LIE method was derived as an approximation to free energy perturbation, where a linear response is assumed for electrostatic interactions, together with an empirical expression for nonpolar effects. The interaction energies are evaluated from the average of molecular dynamics simulation (MD) energies. The second method is the Mining Minima method developed by Gilson and coworkers57 and applies to the direct computation of the configurational integral of a molecule (eq 4) as a sum of the contribution of low energy states using a Monte Carlo (MC) technique. Both methods require conformational sampling through MD or MC techniques. In comparison, our method has the advantage that a single conformation is sufficient for free energy prediction.
In a recent review of docking and scoring functions Leach et al58 point out that the accuracy of scoring functions has reached a plateau and that there is need for a breakthrough to develop. The path-integral formalism presented in this work offers an alternate perspective on the problem of binding free energy prediction.
The practical complexity of the traditional statistical thermodynamics approach implies that the effect of the assay buffer cannot be reasonably accounted for without significantly increasing the already intractably high computational time of estimating the configurational integral of eq 4. The path integral formalism presented here is elegant not only for its high accuracy, but also for its computational expense, which is minimal.
In its current form, eq 18 has great potential for the development of improved scoring functions for docking algorithms. A more extensive verification of eq 18 is certainly necessary, and at this point this is mainly limited by lack of high quality data-sets such as that of bovine trypsin. It remains to be seen how eq 18 performs in real life drug design scenarios, where typically only a single bound complex is available as a template for the prediction of other bound complexes through docking. It would also greatly help if inhibition experiments were supplemented with experimental determination of the ligands’ diffusivities in their respective buffers, which would remove the reliance on empirical relations such as Wilke-Chang equation that was used in this work.
The authors would like to thank Dr. Abdallah Sayyed-Ahmad and Dr. Anthony D Hill for helpful discussions. This work was supported by a grant from NIH (GM 070989). Computational support from the Minnesota Supercomputing Institute (MSI) is gratefully acknowledged. This work was also partially supported by National Computational Science Alliance under MCA04T033.