- Email: [email protected]
Effect of ligand torsion number on the AutoDock mediated prediction of protein-ligand binding affinity
G Model JIEC 4893 No. of Pages 7
Journal of Industrial and Engineering Chemistry xxx (2019) xxx–xxx
Contents lists available at ScienceDirect
Journal of Industrial and Engineering Chemistry journal homepage: www.elsevier.com/locate/jiec
Effect of ligand torsion number on the AutoDock mediated prediction of protein-ligand binding affinity Dinesh Kumar Sriramulua , Sangwook Wub,c,* , Sun-Gu Leea,* a b c
Department of Polymer Science and Chemical Engineering, Pusan National University, Busan, 609-735, Republic of Korea Department of Physics, Pukyong National University, Busan, 608-737, Republic of Korea Pharmcadd, IntelliumCenter, Centum 6-ro, Haeundae-gu, Busan, 48059, Korea
A R T I C L E I N F O
A B S T R A C T
Article history: Received 27 September 2019 Received in revised form 2 December 2019 Accepted 7 December 2019 Available online xxx
Molecular docking simulation is a useful tool in the prediction of protein-ligand binding affinity on a large scale and has great potential in various application fields such as virtual screening of potential drug molecules. However, the reliability of molecular docking is still weak in the estimation of ligandbinding free energy, which limits the applicability of molecular docking simulation. Ligand torsion number is related to the flexibility of ligand and generally incorporated as a crucial variable in the thermodynamic function of binding free energy. In this study, we investigated how the ligand torsion number has influence on the binding affinity prediction of AutoDock, a popular molecular docking simulation tool. The pKd values of various protein-ligands were estimated by using the binding free energy function of AutoDock and compared with their experimental pKd values. The torsion number dependent comparison showed that the predicted binding affinities were mostly underestimated in the complexes of higher torsion numbers, whereas the underestimated and overestimated cases were relatively balanced at relatively lower torsion numbers. A new weight factor for torsion-free energy term of binding energy function was determined and introduced to make correction to the underestimation of binding affinity of ligands with high torsion numbers. It is expected that the torsion number dependent deviation pattern of AutoDock and its correction strategy are useful in the largescale validation of protein-ligand binding affinity. © 2019 The Korean Society of Industrial and Engineering Chemistry. Published by Elsevier B.V. All rights reserved.
Keywords: Protein-ligand interactions Binding energy AutoDock Torsional free energy Weight factor Torsion number
Introduction Proteins realize their biological functions through their direct interaction with ligands, nucleic acids and other proteins [1–3]. Among these interactions, protein-ligand interactions are involved in various crucial biological processes such as enzyme reaction, signal transduction, and immune response [4]. Understanding protein-ligand interactions at a molecular level is deeply related to a broad area of applications. For instance, it plays an essential role in the design of the protein-ligand binding site as well as in structure-based ligand design. Predicting the affinity of proteinligand binding based on molecular structures is very useful and integral in such applications. Theoretically, the prediction can be achieved by estimating a free energy change that takes in the contribution of enthalpic and entropic changes in a protein-ligand
* Corresponding authors. E-mail addresses: [email protected] (S. Wu), [email protected] (S.-G. Lee).
binding [5]. Although it is still challenging to predict the thermodynamic terms correctly, various methodologies have been developed based on the thermodynamic concept and used in a variety of applications [6]. A molecular docking simulation is a typical approach employed in the study of protein-ligand interactions. In general, several ligand poses are generated inside the binding pocket of the receptor molecule, and the poses are evaluated with the help of a scoring function in molecular docking simulation [7]. The approach has been widely used in the study of the ligand-binding site, ligand orientation, and protein-ligand molecular interactions [8]. Also, scoring functions for molecular docking simulation are mostly devised on the basis of thermodynamic energy terms in protein-ligand binding, and therefore the approach is sometimes utilized in the estimation of protein-ligand binding affinity [9–11]. Primarily, it is popularly used in the estimation on a large scale due to its advantage in the aspect of computational time and capacity. However, there is a significant deviation between predicted and experimental binding affinities, which limits the reliability of docking-based estimation of binding affinity.
https://doi.org/10.1016/j.jiec.2019.12.009 1226-086X/© 2019 The Korean Society of Industrial and Engineering Chemistry. Published by Elsevier B.V. All rights reserved.
Please cite this article in press as: D.K. Sriramulu, et al., Effect of ligand torsion number on the AutoDock mediated prediction of protein-ligand binding affinity, J. Ind. Eng. Chem. (2019), https://doi.org/10.1016/j.jiec.2019.12.009
G Model JIEC 4893 No. of Pages 7
2
D.K. Sriramulu et al. / Journal of Industrial and Engineering Chemistry xxx (2019) xxx–xxx
To overcome the limitation of docking-based affinity prediction, it is essential to understand the factors which induce the deviation of predicted affinity from experimental affinity. The number of rotational or torsion bonds in a ligand is recognized to be an essential variable in the estimation of ligand-binding free energy [12]. It is related to the flexibility of ligand and has a significant influence on the entropic term of free energy. The effect of torsion number on ligand-binding free energy is generally incorporated into the binding energy functions for molecular docking simulations. Despite the importance of torsion number in the estimation of binding energy, it has not been studied systematically how the number of torsion bonds affects the docking-based affinity prediction. There have been only a few studies on the effect of torsion number which focused on the prediction efficiency of ligand-binding pose [13,14]. Here, we estimated the binding affinities of various protein-ligand complexes by using a molecular docking simulation method and analyzed the effect of torsion number on the prediction of protein-ligand binding affinity. AutoDock (version 4.2.6) [15] is a popular molecular docking simulation tool that has been used in a variety of protein-ligand interaction studies. An essential feature of AutoDock is that it's scoring function is a binding free energy function composed of thermodynamically meaningful energy terms, i.e., intermolecular energy term and torsion energy term. Intermolecular binding energy term includes the contributions from van der Waal’s interaction, hydrogen bonds, electrostatic potential, and desolvation energy in the protein-ligand interaction. The torsion energy term incorporates the free energy released due to the rotational bonds present in the ligand [16]. In this study, various proteinligand complexes taken from a database were docked, and their binding affinities were estimated by using the default binding energy function of AutoDock. The predicted binding affinities were compared with experimental values, and their deviations were analyzed depending on the number of torsion bonds in the ligand. The torsion number effect was further investigated for intermolecular energy term and torsion energy term. From these analyses, we were able to find out a specific torsion number dependent deviation pattern in the binding affinity prediction. Finally, an attempt was made to reduce the deviation by modifying a weight factor in the torsion energy term of the binding free energy function. Methods Dataset used in the study The protein-ligand complexes were taken from different databases. A total of 172 protein-ligand complexes in Supplementary Table 1 were extracted from the PDBbind [17] database and used in the study on the torsion number effect. For the study on validation of modified energy function, ten complexes in Supplementary Table 2 were selected from BindingDB [18] and BindingMoad [19] in a way that those structures are exclusively available only in these databases. In the selection of protein-ligand complexes from the databases, the protein structure of a single chain (no multiple or multi-chain structures) and availability of experimental Kd values were used as a criterion to avoid complexity in analysis. All these structures were downloaded by the RCSB [20].
region to that of its crystal structure and hence avoid binding to the other areas. Protein and ligand coordinates were separated for each complex. With the help of AutoDock tools [15] (ADT, version 1.5.6), these protein and ligand structures were processed to a format recognized by the ADT (*.pdbqt files) by adding all hydrogen atoms, Gasteiger charges and merging the non-polar hydrogen atoms. The hydrogens are added to the receptor and ligand structures using Babel in Autodock [15]. The number of torsions for the ligand are determined by the Autotors, an inbuilt tool in the ADT. The torsional root for the ligand molecule is usually a rigid part of the ligand molecule from which branches with rotatable bonds emerge [23]. The Autodock often tries to identify the root for the molecule unless the user specifies it. In this work, root was automatically selected by the Autodock. The grid box was positioned based on the information from the PDBsum and CastP servers. The atom specific affinity maps for all ligand atom types, electrostatic and desolvation potentials were computed with the help of Autogrid (version 4.2.6) [15]. The docking simulation for every proteinligand pair consisted of 100 iterations following Genetic Algorithm methods. The population size was set to 150 with 25 million energy evaluations using a Lamarckian Genetic Algorithm local search method in AutoDock (version 4.2.6) [15]. All the other parameters like the number of generations, the maximum number of top individuals to survive to the next generation, mutation rate and the crossover rate were set to their default value. The total minimum energy resulting from the 100 docking runs were used as the binding energy for analysis. Conversion of energy terms to affinity terms The minimum binding energy (DGbinding) obtained from docking was converted to a dissociation equilibrium constant (Kd) based on their relation, DGbinding ¼ RT lnK d [3] and the Kd was converted to pKd, used as a binding affinity in this study. Eventually, the DG
equation used in the conversion was pKd = log10(e RT ). In this calculation, 1.987 cal K1 mol1 and 298 K were used as R (gas constant) and T (temperature) values. The intermolecular energy term (DGintermolecular) and torsion energy term (DGtor) obtained from molecular docking were converted to the pKinter and pKtor corresponding to the pKd term for binding free energy in the same manner. Results Dataset collection, molecular docking, and prediction of binding affinity A dataset of 172 protein-ligand complexes was taken from the PDBbind database. Molecular docking was performed for these 172 pairs with the help of AutoDock. In the docking, the ligand-binding site in the PDB structure of the complex was re-targeted. The lowest binding energy values from the several dockings were taken, converted to pKd values, and then compared with their corresponding experimental pKd values. Supplementary Table 1 shows the selected protein-ligand complexes, their experimental pKd values and predicted pKd values. Effect of torsion number on experimental and predicted binding affinities
Molecular docking Molecular docking of the protein-ligand complexes was performed with the help of AutoDock. The ligand-binding site of each complex was identified by PDBsum [21] and CASTp [22] servers. This process was done to make the ligand bind in a similar
To investigate the effect of torsion number present in a ligand on the prediction of ligand-binding affinity, the experimental and predicted pKd values were plotted against the torsion number of a ligand in Fig. 1(a) and (b), respectively. In the figures, the average values were also plotted for respective torsion numbers to
Please cite this article in press as: D.K. Sriramulu, et al., Effect of ligand torsion number on the AutoDock mediated prediction of protein-ligand binding affinity, J. Ind. Eng. Chem. (2019), https://doi.org/10.1016/j.jiec.2019.12.009
G Model JIEC 4893 No. of Pages 7
D.K. Sriramulu et al. / Journal of Industrial and Engineering Chemistry xxx (2019) xxx–xxx
3
Fig. 1. The effect of the torsion number on experimental and predicted binding affinities. (a) Distribution of experimental binding affinity (pKd(experimental)) values based on the torsional number of the ligand. The experimental binding affinity (Kd) values were obtained from the PDBbind database and were converted to their corresponding pKd values. (b) Distribution of predicted binding affinity (pKd(predicted)) values based on the torsional number of the ligand. The binding energy (DGbinding) values were obtained from the docking run for and were converted to their corresponding pKd values. (c) Distribution of difference between experimental and predicted binding affinity values based on the torsional number of the ligand. The values are calculated by subtracting the predicted values from the experimental values. If the difference between these values are greater than zero, then it indicates that the prediction is underestimated and if the difference is less than zero, then the prediction is overestimated. In the figures, the red dot indicates the average of data for each torsion number (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article).
understand their torsion number dependent patterns. The averages of experimental pKd values were not significantly varied and remained in a similar range as the torsion number increased (Fig. 1(a)). On the other hand, the averages of predicted pKd values were overall lower and not uniformly distributed compared to those of the experimental values. The average predicted affinity values showed a pattern of decrease as the ligand torsion number increased (Fig. 1(b)). To interpret the torsion number effect on the deviation between predicted affinity and experimental affinity more clearly, the difference in the affinity values was calculated by subtracting predicted pKd values from experimental pKd values and plotted against the torsion number (Fig. 1(c)). In the plot, the difference of zero means a perfect correlation between the experimental and predicted affinity values. The positive and negative values of difference indicate that the predicted affinity is underestimated and overestimated, respectively, compared to the experimental affinity. This result suggests that the predicted affinities are overall underestimated. For example, assuming that the difference in the range of 1 to +1 is the case that predicted values and experimental values are close, 124 protein-ligand complexes were out of the range, among which 27 complexes were in the overestimated region while the rest 97 complexes were in the underestimated region. Meanwhile, the plot exhibited a different pattern depending on the torsion number. The underestimated prediction became predominant in the complexes of higher torsion number, whereas the underestimated and overestimated cases were relatively balanced in the complexes of relatively lower torsion number. This pattern could be confirmed from the estimation of the average
of the differences for the complexes of torsion number <15 and 15, estimated to be 0.85 (2.33) and 2.22(2.4), respectively. The average of the difference values in each torsion number also confirmed this pattern. These results indicate that the predicted affinities for the ligands with higher torsion numbers are biased toward being underestimated. In other words, as the torsion number of the ligand increases, the predicted binding affinity values are patterned to be lower than their corresponding experimental values. Effect of torsion number on intermolecular and torsion energy terms The free energy of binding (DGbinding) in AutoDock is calculated as the sum of free energy of intermolecular energy term (DGintermolecular) and the free energy of the torsion term (DGtorsion) (Eq. (1)).
DGbinding = DGintermolecular + DGtorsion
(1)
In the equation, the intermolecular energy term includes the contributions of van der Waal’s interaction, hydrogen bonds, electrostatic potential, and desolvation energy term. The torsion energy term is estimated simply by multiplying a weight factor (Wtor) and torsion number of ligand (Ntor) (Eq. (2)).
DGtorsion =Wtor * Ntor
(2)
Therefore, the torsion energy term would be directly affected by ligand torsion number, while the effect of torsion number on the intermolecular term is unclear. Here, to understand how the intermolecular energy term and torsion energy term influence the identified deviation pattern of predicted binding affinity, i.e.,
Please cite this article in press as: D.K. Sriramulu, et al., Effect of ligand torsion number on the AutoDock mediated prediction of protein-ligand binding affinity, J. Ind. Eng. Chem. (2019), https://doi.org/10.1016/j.jiec.2019.12.009
G Model JIEC 4893 No. of Pages 7
4
D.K. Sriramulu et al. / Journal of Industrial and Engineering Chemistry xxx (2019) xxx–xxx
the bias toward the underestimation of predicted pKd at higher torsion numbers, the effects of torsion number on intermolecular energy term and torsion energy term were separately investigated. As the binding free energy was converted to the binding affinity term, i.e., pKd, the intermolecular energy term, and torsion energy term were converted to corresponding terms as described in the Methods section, designated as pKinter and pKtor, respectively. Now pKd can be defined as a sum of pKinter and pKtor. Fig. 2(a) shows the pKinter values of the dataset and the average of the values in each torsion number plotted against the torsion number. The plot roughly indicates a pattern of a slight increase in pKinter values as the torsion number increases, but the torsion number dependency was fragile. On the other hand, the plot of pKtor against torsion number showed an apparent torsion number dependency, a linear decrease proportional to the torsion number (Fig. 2(b)). This result is reasonable, considering that the torsion energy term is calculated by multiplying the constant weight factor (0.29) and torsion number of the ligand. In the plot, the decrease was determined to be 0.21 per torsion number, and the average pKtor for the entire dataset was estimated to be 2.18. The above results suggest that the thermodynamic strategy in the prediction of binding affinity which adds the intermolecular term and torsion term is roughly effective. On average, the estimated pKinter values in Fig. 2(a) were higher than the experimental pKd values in Fig. 1(a), and therefore the addition of negative pKtor values to the pKinter values could be useful to find a predicted binding affinity close to an experimental affinity. Meanwhile, the torsion number dependencies of pKinter and pKtor values in Fig. 2 indicate that there are two possibilities in the underestimation-biased prediction of affinity for the ligands with higher torsion number observed in Fig. 1; the torsion number dependent increase of pKinter values is too weak, or the torsion number dependent decrease of pKtor values is too substantial. This finding also implies that the bias of underestimation would be corrected if the current torsion number dependency of the intermolecular term and torsion term can be modified. Modification of torsion-free energy term for correcting the underestimation-biased deviation Our next study focused on reducing the underestimation-biased deviation of predicted binding affinities for the higher torsion number ligands. Theoretically, the deviation of predicted affinity would be caused by the incorrect estimation of the intermolecular energy term and torsion energy term. Here, we assumed that the
underestimation-biased deviation was created by the higher torsion number dependency of torsion energy term just because the torsion energy term includes the torsion number (Ntor) as a direct variable. Therefore, it was hypothesized that the modification of weight factor (Wtor) value (0.29), could weaken the current torsion number dependency and eventually relieve the bias toward underestimation. To check whether such approach can be useful, we computed an effective Wtor to improve the correlation between the predicted binding energy and the experimental binding energy for each complex based on the Eq. (3). Effective Wtor = (DGexperimental DGintermolecular)/Ntor
(3)
Fig. 3 shows the effective Wtor values for the entire dataset plotted against the torsion number. It was observed that the effective Wtor values were significantly varying. At lower torsion number ligand molecules, there was no bias in distribution. The effective Wtor values were found to be equally distributed above and below the default value used in AutoDock. This result suggests that it is not easy to achieve a modified Wtor value that can be universally applied to the correction of the deviations of ligands with lower torsion numbers. On the other hand, as the torsion number increases, the variability of effective Wtor values decreased significantly, and the effective Wtor values were found to be mostly lower than the default value used in AutoDock (0.29). This fact indicates that the underestimation-biased deviation of predicted affinities for the higher torsion number ligands can be solved by employing a modified Wtor value lower than current 0.29. To test the above hypothesis, the binding affinities for the ligands with torsion number of >12 (48 complexes) were estimated by employing an α Wtor (0.29) * Ntor instead of Wtor (0.29) * Ntor in the estimation of DGtorsion. The value of modification factor (α) was determined by dividing 0.129, an average of the effective Wtor values for ligands with torsion number of >12 by 0.29, the default Wtor. To check the effect of modification, the difference between the experimental and modified predicted affinity values were plotted against the torsion number and compared with the differences before the change (Fig. 4(a)). The plot clearly shows that the use of the modification factor could relieve the underestimation-biased deviation problem of high torsion number ligands. The average difference between the two affinities after the modification was estimated to be 0.01 (2.2), whereas the average difference before modification was 2.19 (2.27). The frequency distribution of difference also confirmed that the bias was clearly reduced by using the modification factor (Fig. 4(b)).
Fig. 2. The effect of the torsion number on the intermolecular and torsional energy terms. (a) Distribution of the pKintermolecular [log10(DGintermolecular)] of 172 ligands based on the corresponding torsional number. The intermolecular energy values obtained from 172 individual docking results are converted to their respective pKd values. The red dot indicates the average of data for each torsion number. (b) Distribution of the pKtorsion[log10(DGtorsion)] of 172 ligands based on the corresponding torsional number. The torsional free energy values obtained from 172 individual docking results are converted to their respective pKd values (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article).
Please cite this article in press as: D.K. Sriramulu, et al., Effect of ligand torsion number on the AutoDock mediated prediction of protein-ligand binding affinity, J. Ind. Eng. Chem. (2019), https://doi.org/10.1016/j.jiec.2019.12.009
G Model JIEC 4893 No. of Pages 7
D.K. Sriramulu et al. / Journal of Industrial and Engineering Chemistry xxx (2019) xxx–xxx
5
Fig. 3. Distribution of effective Wtor values plotted against ligand torsion numbers. The effective Wtor values were calculated based on the experimental binding affinity, predicted intermolecular energy value and the torsional number of the ligand. The dotted line in red indicates the default value used by the AutoDock. In order to show a better representation of the plot, two entries were removed as these two values were too deviated from the rest of the values (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article).
Fig. 4. Comparison of difference between experimental and predicted binding affinity values for high torsion number ligand complexes before and after use of modification factor. (a) Distribution of difference values dependent on ligand torsion number. The values were calculated by subtracting the predicted affinities from the experimental affinities for the 48 complexes with >12 of ligand torsion number in the test dataset. (b) Frequency distribution of difference values. The frequency distribution ratio was calculated by dividing the number of complexes in the given range of difference by total number of complexes. In this plot, the frequencies in the difference range of 0–2 were same before and after the use of modification factor.
Application of modification factor to the complexes in the validation dataset Finally, the strategy of using the modification factor was validated by applying it for the binding affinity prediction of the complexes out of the test dataset used above. 17 protein structures with the ligand torsion number >12 were taken from BindingDB and BindingMoad databases, designated to be ‘validation dataset,’ and their binding affinities were calculated based on the default Wtor (0.29) and with the use of the modification factor(α) value determined above. Supplementary Table 2 shows the 17 proteinligand complexes, their experimental binding affinities, and predicted binding affinities. Based on the data, the differences between the experimental pKd and the predicted binding pKd were calculated and plotted against the torsion number (Fig. 5(a)). The difference values for the prediction based on default Wtor were all
positive, except for one, and their average was identified to be 2.45 with a standard deviation of 1.54. These results confirmed that the use of default Wtor in the prediction of binding affinity resulted in the biased pattern of underestimation for higher torsion numbers. On the other hand, the prediction based on the ‘α’ factor didn’t show the biased pattern of difference. The average difference was determined to be 0.31 with a standard deviation of 1.45. This data implicates that the prediction accuracy could also be somewhat improved by relieving the bias of underestimation. The frequency distribution of difference also confirmed that the bias was reduced by using the modification factor in the prediction of binding affinities for the validation dataset (Fig. 5(b)). Discussion Prediction of protein-ligand binding affinity is of vital importance in the field of structure-based drug discovery. Currently,
Please cite this article in press as: D.K. Sriramulu, et al., Effect of ligand torsion number on the AutoDock mediated prediction of protein-ligand binding affinity, J. Ind. Eng. Chem. (2019), https://doi.org/10.1016/j.jiec.2019.12.009
G Model JIEC 4893 No. of Pages 7
6
D.K. Sriramulu et al. / Journal of Industrial and Engineering Chemistry xxx (2019) xxx–xxx
Fig. 5. Comparison of difference between experimental and predicted binding affinity values for the validation dataset before and after use of modification factor. (a) Distribution of difference values dependent on ligand torsion number. The values were calculated by subtracting the predicted affinities from the experimental affinities for the 17 complexes with >12 of ligand torsion number in the validation dataset obtained from the BindingDB and BindingMoad databases. (b) Frequency distribution of difference values. The frequency distribution ratio was calculated by dividing the number of complexes in the given range of difference by total number of complexes.
Molecular Dynamics (MD) simulation or MD-related approaches are recognized to be promising in the prediction, but their computational time, cost, and complexity could be a limitation in the large-scale estimation of binding affinity. The molecular docking approach has been popularly used in the prediction of protein-ligand interactions due to its availability for the screening of ligands on a large scale. However, the predictability of binding affinity is very weak because of too many simplifications in its evaluation. It may be meaningful to attempt to solve the problems in the docking-based affinity prediction. In this study, we have investigated how the torsion number of ligands influences the affinity prediction of AutoDock. It was identified that the pattern of predicted affinity is biased toward underestimation as the torsion number increases. Based on the discovered pattern, we attempted to correct the torsion-free energy term of AutoDock and demonstrated that the underestimation-biased prediction for high torsion ligands could be relieved by using a modification factor. It is expected that the discovered deviation pattern and its correction strategy are useful in the design of the large-scale screening process of ligands using molecular docking approach. The underestimation-biased prediction pattern for high torsion ligands is an only ligand-dependent pattern, and the correction strategy would be useful regardless of target proteins. To relieve the underestimation-biased deviation at higher torsion numbers, we have used the strategy using a modified factor (α), which can weaken the torsion number dependency of torsion energy term. However, this doesn’t mean that the predicted intermolecular energy term is correct, and only the predicted torsion-free energy term is incorrect. Theoretically, the contribution of intermolecular energy term to the total binding energy is more significant as the torsion number decreases. As shown in Fig. 1(c), the deviation of the predicted affinity from the experimental values was significantly high for the ligands with low torsion numbers, and the predicted binding affinities were not biased toward overestimation or underestimation. This pattern indicates that there are significant deviations in the prediction of intermolecular energy term. On the other hand, the torsion term effect on the total binding free energy becomes more significant as the torsion number increases. Therefore, the deviation of predicted affinity for the lower torsion numbers might be caused by the incorrect estimation of torsion energy term in addition to the inaccurate evaluation of intermolecular energy term. Meanwhile, the underestimation-biased trend of higher torsion number is presumed to be contributed mainly by the torsion energy term effect. If the intermolecular term effect is more significant than the
torsion term, the deviation would have followed the non-biased pattern observed in low torsion numbers. If the underestimation-biased trend of predicted affinity for high torsion numbers is mainly from the torsion energy term as hypothesized above, the modification factor may have a physical meaning related to the torsion energy term. Based on the effective Wtor values in Fig. 3, the values for the ligands with high torsion numbers are mostly positive values and lower than the currently used Wtor, 0.29. The positive values indicate that the torsion energy term is positive, and the torsion energy before ligandbinding is higher than after ligand-binding, as understood currently. On the other hand, the lower values than 0.29 suggest that the torsion energy difference between before binding and after binding is not so more significant than expected. There are two possibilities responsible for making such an assumption. First is that the ligands with high torsion numbers are more compact than expected when it is surrounded by water molecules in the solution. The other possibility is that the ligands with high torsional number may have far more flexible conformations than expected inside the binding pocket. Of course, there is a possibility that this kind of interpretation of the physical meaning of the torsional energy term is not meaningful because many other important factors, such as protein flexibility were not considered in our study. The Wtor value of Autodock has been generally used regardless of ligand torsion number. In this study, we found that there was an underestimation-biased deviation pattern for the ligands with a relatively higher torsion number when the Wtor was used and demonstrated the idea that the use of modification factor could relieve the bias in the affinity prediction of high torsion number ligands. To apply the approach more practically, further studies may be needed. For instance, more detailed statistical analyses and optimization studies need to be performed to achieve a more optimized modification factor. In addition, the use of modification factor is not effective to low torsion ligands. For the ligands with torsion number of 412, the average of the effective Wtor values was determined to be 0.21, closer to the default Wtor value, and there was no clear deviation pattern in the prediction of binding affinities. As mentioned above, the prediction of more accurate intermolecular energy term is necessarily needed to reduce the deviation of affinity prediction for low torsion ligands. Declaration of interests None.
Please cite this article in press as: D.K. Sriramulu, et al., Effect of ligand torsion number on the AutoDock mediated prediction of protein-ligand binding affinity, J. Ind. Eng. Chem. (2019), https://doi.org/10.1016/j.jiec.2019.12.009
G Model JIEC 4893 No. of Pages 7
D.K. Sriramulu et al. / Journal of Industrial and Engineering Chemistry xxx (2019) xxx–xxx
Acknowledgements This research was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea Government (MSIT) (No. 2018R1A2B6001670). Appendix A. Supplementary data Supplementary material related to this article can be found, in the online version, at doi:https://doi.org/10.1016/j.jiec.2019.12.009. References [1] C. Bissantz, B. Kuhn, M. Stahl, J. Med. Chem. 53 (2010) 5061. [2] J.L. Medina-Franco, O. Mendez-Lucio, K. Martinez-Mayorga, Adv. Protein Chem. Struct. Biol. 96 (2014) 1. [3] Xing Du, Yi Li, Yuan-Ling Xia, Shi-Meng Ai, Jing Liang, Peng Sang, Xing-Lai Ji, S.Q. Liu, Int. J. Mol. Sci. 17 (2016) 17. [4] W.J, H.N. Motlagh, J. Li, V.J. Hilser, Nature 508 (2014) 331. [5] G. Klebe, H.J. Bohm, J. Recept. Signal Transduct. Res. 17 (1997) 459. [6] I.A. Guedes, C.S. de Magalhaes, L.E. Dardenne, Biophys. Rev. 6 (2014) 75. [7] P. Ferrara, H. Gohlke, D.J. Price, G. Klebe, C.L. Brooks III, J. Med. Chem. 47 (2004) 3032. [8] D. Plewczynski, M. Lazniewski, R. Augustyniak, K. Ginalski, J. Comput. Chem. 32 (2011) 742.
7
[9] R.D. Smith, J.B. Dunbar Jr., P.M. Ung, E.X. Esposito, C.Y. Yang, S. Wang, H.A. Carlson, J. Chem. Inf. Model 51 (2011) 2115. [10] Z. Gaieb, S. Liu, S. Gathiaka, M. Chiu, H. Yang, C. Shao, V.A. Feher, W.P. Walters, B. Kuhn, M.G. Rudolph, S.K. Burley, M.K. Gilson, R.E. Amaro, J Comput. Aided Mol. Des. 32 (2018) 1. [11] N.S. Pagadala, K. Syed, T.J, Biophys. Rev. 9 (2017) 91. [12] Z. Wang, Y. Kang, D. Li, H. Sun, X. Dong, X. Yao, L. Xu, S. Chang, Y. Li, T. Hou, J. Phys. Chem. B 122 (2018) 2544. [13] Z. Wang, H. Sun, X. Yao, D. Li, L. Xu, Y. Li, S. Tian, T. Hou, Phys. Chem. Chem. Phys. 18 (2016) 12964. [14] S. Mukherjee, T.E. Balius, R.C. Rizzo, J. Chem. Inf. Model 50 (2010) 1986. [15] G.M. Morris, R. Huey, W. Lindstrom, M.F. Sanner, R.K. Belew, D.S. Goodsell, A.J. Olson, J. Comput. Chem. 30 (2009) 2785. [16] R. Huey, G.M. Morris, A.J. Olson, D.S. Goodsell, J. Comput. Chem. 28 (2007) 1145. [17] Renxiao Wang, Xueliang Fang, A. Yipin Lu, S. Wang, J. Med. Chem. 47 (2004) 2977. [18] L.T. Gilson, M.K. Baitaluk, M. Nicola, G. Hwang, L. Chong, J. Nucleic Acids Res. 4 (2016) D1045. [19] tL. Hu, M.L. Benson, R.D. Smith, M.G. Lerner, H. Carlson, Proteins 60 (2005) 333. [20] H.M. Berman, J. Westbrook, Z. Feng, G. Gilliland, T.N. Bhat, H. Weissig, I.N. Shindyalov, P.E. Bourne, Nucleic Acids Res. 28 (2000) 235. [21] R.A. Laskowski, E.G. Hutchinson, A.D. Michie, A.C. Wallace, M.L. Jones, J.M. Thornton, Trends Biochem. Sci. 22 (1997) 488. [22] W. Tian, C. Chen, X. Lei, J. Zhao, J. Liang, Nucleic Acids Res. 46 (2018) W363. [23] D.S.G. Garrett, M. Morris, Ruth Huey, Arthur J. Olson, (1996).
Please cite this article in press as: D.K. Sriramulu, et al., Effect of ligand torsion number on the AutoDock mediated prediction of protein-ligand binding affinity, J. Ind. Eng. Chem. (2019), https://doi.org/10.1016/j.jiec.2019.12.009
Journal of Industrial and Engineering Chemistry xxx (2019) xxx–xxx
Contents lists available at ScienceDirect
Journal of Industrial and Engineering Chemistry journal homepage: www.elsevier.com/locate/jiec
Effect of ligand torsion number on the AutoDock mediated prediction of protein-ligand binding affinity Dinesh Kumar Sriramulua , Sangwook Wub,c,* , Sun-Gu Leea,* a b c
Department of Polymer Science and Chemical Engineering, Pusan National University, Busan, 609-735, Republic of Korea Department of Physics, Pukyong National University, Busan, 608-737, Republic of Korea Pharmcadd, IntelliumCenter, Centum 6-ro, Haeundae-gu, Busan, 48059, Korea
A R T I C L E I N F O
A B S T R A C T
Article history: Received 27 September 2019 Received in revised form 2 December 2019 Accepted 7 December 2019 Available online xxx
Molecular docking simulation is a useful tool in the prediction of protein-ligand binding affinity on a large scale and has great potential in various application fields such as virtual screening of potential drug molecules. However, the reliability of molecular docking is still weak in the estimation of ligandbinding free energy, which limits the applicability of molecular docking simulation. Ligand torsion number is related to the flexibility of ligand and generally incorporated as a crucial variable in the thermodynamic function of binding free energy. In this study, we investigated how the ligand torsion number has influence on the binding affinity prediction of AutoDock, a popular molecular docking simulation tool. The pKd values of various protein-ligands were estimated by using the binding free energy function of AutoDock and compared with their experimental pKd values. The torsion number dependent comparison showed that the predicted binding affinities were mostly underestimated in the complexes of higher torsion numbers, whereas the underestimated and overestimated cases were relatively balanced at relatively lower torsion numbers. A new weight factor for torsion-free energy term of binding energy function was determined and introduced to make correction to the underestimation of binding affinity of ligands with high torsion numbers. It is expected that the torsion number dependent deviation pattern of AutoDock and its correction strategy are useful in the largescale validation of protein-ligand binding affinity. © 2019 The Korean Society of Industrial and Engineering Chemistry. Published by Elsevier B.V. All rights reserved.
Keywords: Protein-ligand interactions Binding energy AutoDock Torsional free energy Weight factor Torsion number
Introduction Proteins realize their biological functions through their direct interaction with ligands, nucleic acids and other proteins [1–3]. Among these interactions, protein-ligand interactions are involved in various crucial biological processes such as enzyme reaction, signal transduction, and immune response [4]. Understanding protein-ligand interactions at a molecular level is deeply related to a broad area of applications. For instance, it plays an essential role in the design of the protein-ligand binding site as well as in structure-based ligand design. Predicting the affinity of proteinligand binding based on molecular structures is very useful and integral in such applications. Theoretically, the prediction can be achieved by estimating a free energy change that takes in the contribution of enthalpic and entropic changes in a protein-ligand
* Corresponding authors. E-mail addresses: [email protected] (S. Wu), [email protected] (S.-G. Lee).
binding [5]. Although it is still challenging to predict the thermodynamic terms correctly, various methodologies have been developed based on the thermodynamic concept and used in a variety of applications [6]. A molecular docking simulation is a typical approach employed in the study of protein-ligand interactions. In general, several ligand poses are generated inside the binding pocket of the receptor molecule, and the poses are evaluated with the help of a scoring function in molecular docking simulation [7]. The approach has been widely used in the study of the ligand-binding site, ligand orientation, and protein-ligand molecular interactions [8]. Also, scoring functions for molecular docking simulation are mostly devised on the basis of thermodynamic energy terms in protein-ligand binding, and therefore the approach is sometimes utilized in the estimation of protein-ligand binding affinity [9–11]. Primarily, it is popularly used in the estimation on a large scale due to its advantage in the aspect of computational time and capacity. However, there is a significant deviation between predicted and experimental binding affinities, which limits the reliability of docking-based estimation of binding affinity.
https://doi.org/10.1016/j.jiec.2019.12.009 1226-086X/© 2019 The Korean Society of Industrial and Engineering Chemistry. Published by Elsevier B.V. All rights reserved.
Please cite this article in press as: D.K. Sriramulu, et al., Effect of ligand torsion number on the AutoDock mediated prediction of protein-ligand binding affinity, J. Ind. Eng. Chem. (2019), https://doi.org/10.1016/j.jiec.2019.12.009
G Model JIEC 4893 No. of Pages 7
2
D.K. Sriramulu et al. / Journal of Industrial and Engineering Chemistry xxx (2019) xxx–xxx
To overcome the limitation of docking-based affinity prediction, it is essential to understand the factors which induce the deviation of predicted affinity from experimental affinity. The number of rotational or torsion bonds in a ligand is recognized to be an essential variable in the estimation of ligand-binding free energy [12]. It is related to the flexibility of ligand and has a significant influence on the entropic term of free energy. The effect of torsion number on ligand-binding free energy is generally incorporated into the binding energy functions for molecular docking simulations. Despite the importance of torsion number in the estimation of binding energy, it has not been studied systematically how the number of torsion bonds affects the docking-based affinity prediction. There have been only a few studies on the effect of torsion number which focused on the prediction efficiency of ligand-binding pose [13,14]. Here, we estimated the binding affinities of various protein-ligand complexes by using a molecular docking simulation method and analyzed the effect of torsion number on the prediction of protein-ligand binding affinity. AutoDock (version 4.2.6) [15] is a popular molecular docking simulation tool that has been used in a variety of protein-ligand interaction studies. An essential feature of AutoDock is that it's scoring function is a binding free energy function composed of thermodynamically meaningful energy terms, i.e., intermolecular energy term and torsion energy term. Intermolecular binding energy term includes the contributions from van der Waal’s interaction, hydrogen bonds, electrostatic potential, and desolvation energy in the protein-ligand interaction. The torsion energy term incorporates the free energy released due to the rotational bonds present in the ligand [16]. In this study, various proteinligand complexes taken from a database were docked, and their binding affinities were estimated by using the default binding energy function of AutoDock. The predicted binding affinities were compared with experimental values, and their deviations were analyzed depending on the number of torsion bonds in the ligand. The torsion number effect was further investigated for intermolecular energy term and torsion energy term. From these analyses, we were able to find out a specific torsion number dependent deviation pattern in the binding affinity prediction. Finally, an attempt was made to reduce the deviation by modifying a weight factor in the torsion energy term of the binding free energy function. Methods Dataset used in the study The protein-ligand complexes were taken from different databases. A total of 172 protein-ligand complexes in Supplementary Table 1 were extracted from the PDBbind [17] database and used in the study on the torsion number effect. For the study on validation of modified energy function, ten complexes in Supplementary Table 2 were selected from BindingDB [18] and BindingMoad [19] in a way that those structures are exclusively available only in these databases. In the selection of protein-ligand complexes from the databases, the protein structure of a single chain (no multiple or multi-chain structures) and availability of experimental Kd values were used as a criterion to avoid complexity in analysis. All these structures were downloaded by the RCSB [20].
region to that of its crystal structure and hence avoid binding to the other areas. Protein and ligand coordinates were separated for each complex. With the help of AutoDock tools [15] (ADT, version 1.5.6), these protein and ligand structures were processed to a format recognized by the ADT (*.pdbqt files) by adding all hydrogen atoms, Gasteiger charges and merging the non-polar hydrogen atoms. The hydrogens are added to the receptor and ligand structures using Babel in Autodock [15]. The number of torsions for the ligand are determined by the Autotors, an inbuilt tool in the ADT. The torsional root for the ligand molecule is usually a rigid part of the ligand molecule from which branches with rotatable bonds emerge [23]. The Autodock often tries to identify the root for the molecule unless the user specifies it. In this work, root was automatically selected by the Autodock. The grid box was positioned based on the information from the PDBsum and CastP servers. The atom specific affinity maps for all ligand atom types, electrostatic and desolvation potentials were computed with the help of Autogrid (version 4.2.6) [15]. The docking simulation for every proteinligand pair consisted of 100 iterations following Genetic Algorithm methods. The population size was set to 150 with 25 million energy evaluations using a Lamarckian Genetic Algorithm local search method in AutoDock (version 4.2.6) [15]. All the other parameters like the number of generations, the maximum number of top individuals to survive to the next generation, mutation rate and the crossover rate were set to their default value. The total minimum energy resulting from the 100 docking runs were used as the binding energy for analysis. Conversion of energy terms to affinity terms The minimum binding energy (DGbinding) obtained from docking was converted to a dissociation equilibrium constant (Kd) based on their relation, DGbinding ¼ RT lnK d [3] and the Kd was converted to pKd, used as a binding affinity in this study. Eventually, the DG
equation used in the conversion was pKd = log10(e RT ). In this calculation, 1.987 cal K1 mol1 and 298 K were used as R (gas constant) and T (temperature) values. The intermolecular energy term (DGintermolecular) and torsion energy term (DGtor) obtained from molecular docking were converted to the pKinter and pKtor corresponding to the pKd term for binding free energy in the same manner. Results Dataset collection, molecular docking, and prediction of binding affinity A dataset of 172 protein-ligand complexes was taken from the PDBbind database. Molecular docking was performed for these 172 pairs with the help of AutoDock. In the docking, the ligand-binding site in the PDB structure of the complex was re-targeted. The lowest binding energy values from the several dockings were taken, converted to pKd values, and then compared with their corresponding experimental pKd values. Supplementary Table 1 shows the selected protein-ligand complexes, their experimental pKd values and predicted pKd values. Effect of torsion number on experimental and predicted binding affinities
Molecular docking Molecular docking of the protein-ligand complexes was performed with the help of AutoDock. The ligand-binding site of each complex was identified by PDBsum [21] and CASTp [22] servers. This process was done to make the ligand bind in a similar
To investigate the effect of torsion number present in a ligand on the prediction of ligand-binding affinity, the experimental and predicted pKd values were plotted against the torsion number of a ligand in Fig. 1(a) and (b), respectively. In the figures, the average values were also plotted for respective torsion numbers to
Please cite this article in press as: D.K. Sriramulu, et al., Effect of ligand torsion number on the AutoDock mediated prediction of protein-ligand binding affinity, J. Ind. Eng. Chem. (2019), https://doi.org/10.1016/j.jiec.2019.12.009
G Model JIEC 4893 No. of Pages 7
D.K. Sriramulu et al. / Journal of Industrial and Engineering Chemistry xxx (2019) xxx–xxx
3
Fig. 1. The effect of the torsion number on experimental and predicted binding affinities. (a) Distribution of experimental binding affinity (pKd(experimental)) values based on the torsional number of the ligand. The experimental binding affinity (Kd) values were obtained from the PDBbind database and were converted to their corresponding pKd values. (b) Distribution of predicted binding affinity (pKd(predicted)) values based on the torsional number of the ligand. The binding energy (DGbinding) values were obtained from the docking run for and were converted to their corresponding pKd values. (c) Distribution of difference between experimental and predicted binding affinity values based on the torsional number of the ligand. The values are calculated by subtracting the predicted values from the experimental values. If the difference between these values are greater than zero, then it indicates that the prediction is underestimated and if the difference is less than zero, then the prediction is overestimated. In the figures, the red dot indicates the average of data for each torsion number (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article).
understand their torsion number dependent patterns. The averages of experimental pKd values were not significantly varied and remained in a similar range as the torsion number increased (Fig. 1(a)). On the other hand, the averages of predicted pKd values were overall lower and not uniformly distributed compared to those of the experimental values. The average predicted affinity values showed a pattern of decrease as the ligand torsion number increased (Fig. 1(b)). To interpret the torsion number effect on the deviation between predicted affinity and experimental affinity more clearly, the difference in the affinity values was calculated by subtracting predicted pKd values from experimental pKd values and plotted against the torsion number (Fig. 1(c)). In the plot, the difference of zero means a perfect correlation between the experimental and predicted affinity values. The positive and negative values of difference indicate that the predicted affinity is underestimated and overestimated, respectively, compared to the experimental affinity. This result suggests that the predicted affinities are overall underestimated. For example, assuming that the difference in the range of 1 to +1 is the case that predicted values and experimental values are close, 124 protein-ligand complexes were out of the range, among which 27 complexes were in the overestimated region while the rest 97 complexes were in the underestimated region. Meanwhile, the plot exhibited a different pattern depending on the torsion number. The underestimated prediction became predominant in the complexes of higher torsion number, whereas the underestimated and overestimated cases were relatively balanced in the complexes of relatively lower torsion number. This pattern could be confirmed from the estimation of the average
of the differences for the complexes of torsion number <15 and 15, estimated to be 0.85 (2.33) and 2.22(2.4), respectively. The average of the difference values in each torsion number also confirmed this pattern. These results indicate that the predicted affinities for the ligands with higher torsion numbers are biased toward being underestimated. In other words, as the torsion number of the ligand increases, the predicted binding affinity values are patterned to be lower than their corresponding experimental values. Effect of torsion number on intermolecular and torsion energy terms The free energy of binding (DGbinding) in AutoDock is calculated as the sum of free energy of intermolecular energy term (DGintermolecular) and the free energy of the torsion term (DGtorsion) (Eq. (1)).
DGbinding = DGintermolecular + DGtorsion
(1)
In the equation, the intermolecular energy term includes the contributions of van der Waal’s interaction, hydrogen bonds, electrostatic potential, and desolvation energy term. The torsion energy term is estimated simply by multiplying a weight factor (Wtor) and torsion number of ligand (Ntor) (Eq. (2)).
DGtorsion =Wtor * Ntor
(2)
Therefore, the torsion energy term would be directly affected by ligand torsion number, while the effect of torsion number on the intermolecular term is unclear. Here, to understand how the intermolecular energy term and torsion energy term influence the identified deviation pattern of predicted binding affinity, i.e.,
Please cite this article in press as: D.K. Sriramulu, et al., Effect of ligand torsion number on the AutoDock mediated prediction of protein-ligand binding affinity, J. Ind. Eng. Chem. (2019), https://doi.org/10.1016/j.jiec.2019.12.009
G Model JIEC 4893 No. of Pages 7
4
D.K. Sriramulu et al. / Journal of Industrial and Engineering Chemistry xxx (2019) xxx–xxx
the bias toward the underestimation of predicted pKd at higher torsion numbers, the effects of torsion number on intermolecular energy term and torsion energy term were separately investigated. As the binding free energy was converted to the binding affinity term, i.e., pKd, the intermolecular energy term, and torsion energy term were converted to corresponding terms as described in the Methods section, designated as pKinter and pKtor, respectively. Now pKd can be defined as a sum of pKinter and pKtor. Fig. 2(a) shows the pKinter values of the dataset and the average of the values in each torsion number plotted against the torsion number. The plot roughly indicates a pattern of a slight increase in pKinter values as the torsion number increases, but the torsion number dependency was fragile. On the other hand, the plot of pKtor against torsion number showed an apparent torsion number dependency, a linear decrease proportional to the torsion number (Fig. 2(b)). This result is reasonable, considering that the torsion energy term is calculated by multiplying the constant weight factor (0.29) and torsion number of the ligand. In the plot, the decrease was determined to be 0.21 per torsion number, and the average pKtor for the entire dataset was estimated to be 2.18. The above results suggest that the thermodynamic strategy in the prediction of binding affinity which adds the intermolecular term and torsion term is roughly effective. On average, the estimated pKinter values in Fig. 2(a) were higher than the experimental pKd values in Fig. 1(a), and therefore the addition of negative pKtor values to the pKinter values could be useful to find a predicted binding affinity close to an experimental affinity. Meanwhile, the torsion number dependencies of pKinter and pKtor values in Fig. 2 indicate that there are two possibilities in the underestimation-biased prediction of affinity for the ligands with higher torsion number observed in Fig. 1; the torsion number dependent increase of pKinter values is too weak, or the torsion number dependent decrease of pKtor values is too substantial. This finding also implies that the bias of underestimation would be corrected if the current torsion number dependency of the intermolecular term and torsion term can be modified. Modification of torsion-free energy term for correcting the underestimation-biased deviation Our next study focused on reducing the underestimation-biased deviation of predicted binding affinities for the higher torsion number ligands. Theoretically, the deviation of predicted affinity would be caused by the incorrect estimation of the intermolecular energy term and torsion energy term. Here, we assumed that the
underestimation-biased deviation was created by the higher torsion number dependency of torsion energy term just because the torsion energy term includes the torsion number (Ntor) as a direct variable. Therefore, it was hypothesized that the modification of weight factor (Wtor) value (0.29), could weaken the current torsion number dependency and eventually relieve the bias toward underestimation. To check whether such approach can be useful, we computed an effective Wtor to improve the correlation between the predicted binding energy and the experimental binding energy for each complex based on the Eq. (3). Effective Wtor = (DGexperimental DGintermolecular)/Ntor
(3)
Fig. 3 shows the effective Wtor values for the entire dataset plotted against the torsion number. It was observed that the effective Wtor values were significantly varying. At lower torsion number ligand molecules, there was no bias in distribution. The effective Wtor values were found to be equally distributed above and below the default value used in AutoDock. This result suggests that it is not easy to achieve a modified Wtor value that can be universally applied to the correction of the deviations of ligands with lower torsion numbers. On the other hand, as the torsion number increases, the variability of effective Wtor values decreased significantly, and the effective Wtor values were found to be mostly lower than the default value used in AutoDock (0.29). This fact indicates that the underestimation-biased deviation of predicted affinities for the higher torsion number ligands can be solved by employing a modified Wtor value lower than current 0.29. To test the above hypothesis, the binding affinities for the ligands with torsion number of >12 (48 complexes) were estimated by employing an α Wtor (0.29) * Ntor instead of Wtor (0.29) * Ntor in the estimation of DGtorsion. The value of modification factor (α) was determined by dividing 0.129, an average of the effective Wtor values for ligands with torsion number of >12 by 0.29, the default Wtor. To check the effect of modification, the difference between the experimental and modified predicted affinity values were plotted against the torsion number and compared with the differences before the change (Fig. 4(a)). The plot clearly shows that the use of the modification factor could relieve the underestimation-biased deviation problem of high torsion number ligands. The average difference between the two affinities after the modification was estimated to be 0.01 (2.2), whereas the average difference before modification was 2.19 (2.27). The frequency distribution of difference also confirmed that the bias was clearly reduced by using the modification factor (Fig. 4(b)).
Fig. 2. The effect of the torsion number on the intermolecular and torsional energy terms. (a) Distribution of the pKintermolecular [log10(DGintermolecular)] of 172 ligands based on the corresponding torsional number. The intermolecular energy values obtained from 172 individual docking results are converted to their respective pKd values. The red dot indicates the average of data for each torsion number. (b) Distribution of the pKtorsion[log10(DGtorsion)] of 172 ligands based on the corresponding torsional number. The torsional free energy values obtained from 172 individual docking results are converted to their respective pKd values (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article).
Please cite this article in press as: D.K. Sriramulu, et al., Effect of ligand torsion number on the AutoDock mediated prediction of protein-ligand binding affinity, J. Ind. Eng. Chem. (2019), https://doi.org/10.1016/j.jiec.2019.12.009
G Model JIEC 4893 No. of Pages 7
D.K. Sriramulu et al. / Journal of Industrial and Engineering Chemistry xxx (2019) xxx–xxx
5
Fig. 3. Distribution of effective Wtor values plotted against ligand torsion numbers. The effective Wtor values were calculated based on the experimental binding affinity, predicted intermolecular energy value and the torsional number of the ligand. The dotted line in red indicates the default value used by the AutoDock. In order to show a better representation of the plot, two entries were removed as these two values were too deviated from the rest of the values (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article).
Fig. 4. Comparison of difference between experimental and predicted binding affinity values for high torsion number ligand complexes before and after use of modification factor. (a) Distribution of difference values dependent on ligand torsion number. The values were calculated by subtracting the predicted affinities from the experimental affinities for the 48 complexes with >12 of ligand torsion number in the test dataset. (b) Frequency distribution of difference values. The frequency distribution ratio was calculated by dividing the number of complexes in the given range of difference by total number of complexes. In this plot, the frequencies in the difference range of 0–2 were same before and after the use of modification factor.
Application of modification factor to the complexes in the validation dataset Finally, the strategy of using the modification factor was validated by applying it for the binding affinity prediction of the complexes out of the test dataset used above. 17 protein structures with the ligand torsion number >12 were taken from BindingDB and BindingMoad databases, designated to be ‘validation dataset,’ and their binding affinities were calculated based on the default Wtor (0.29) and with the use of the modification factor(α) value determined above. Supplementary Table 2 shows the 17 proteinligand complexes, their experimental binding affinities, and predicted binding affinities. Based on the data, the differences between the experimental pKd and the predicted binding pKd were calculated and plotted against the torsion number (Fig. 5(a)). The difference values for the prediction based on default Wtor were all
positive, except for one, and their average was identified to be 2.45 with a standard deviation of 1.54. These results confirmed that the use of default Wtor in the prediction of binding affinity resulted in the biased pattern of underestimation for higher torsion numbers. On the other hand, the prediction based on the ‘α’ factor didn’t show the biased pattern of difference. The average difference was determined to be 0.31 with a standard deviation of 1.45. This data implicates that the prediction accuracy could also be somewhat improved by relieving the bias of underestimation. The frequency distribution of difference also confirmed that the bias was reduced by using the modification factor in the prediction of binding affinities for the validation dataset (Fig. 5(b)). Discussion Prediction of protein-ligand binding affinity is of vital importance in the field of structure-based drug discovery. Currently,
Please cite this article in press as: D.K. Sriramulu, et al., Effect of ligand torsion number on the AutoDock mediated prediction of protein-ligand binding affinity, J. Ind. Eng. Chem. (2019), https://doi.org/10.1016/j.jiec.2019.12.009
G Model JIEC 4893 No. of Pages 7
6
D.K. Sriramulu et al. / Journal of Industrial and Engineering Chemistry xxx (2019) xxx–xxx
Fig. 5. Comparison of difference between experimental and predicted binding affinity values for the validation dataset before and after use of modification factor. (a) Distribution of difference values dependent on ligand torsion number. The values were calculated by subtracting the predicted affinities from the experimental affinities for the 17 complexes with >12 of ligand torsion number in the validation dataset obtained from the BindingDB and BindingMoad databases. (b) Frequency distribution of difference values. The frequency distribution ratio was calculated by dividing the number of complexes in the given range of difference by total number of complexes.
Molecular Dynamics (MD) simulation or MD-related approaches are recognized to be promising in the prediction, but their computational time, cost, and complexity could be a limitation in the large-scale estimation of binding affinity. The molecular docking approach has been popularly used in the prediction of protein-ligand interactions due to its availability for the screening of ligands on a large scale. However, the predictability of binding affinity is very weak because of too many simplifications in its evaluation. It may be meaningful to attempt to solve the problems in the docking-based affinity prediction. In this study, we have investigated how the torsion number of ligands influences the affinity prediction of AutoDock. It was identified that the pattern of predicted affinity is biased toward underestimation as the torsion number increases. Based on the discovered pattern, we attempted to correct the torsion-free energy term of AutoDock and demonstrated that the underestimation-biased prediction for high torsion ligands could be relieved by using a modification factor. It is expected that the discovered deviation pattern and its correction strategy are useful in the design of the large-scale screening process of ligands using molecular docking approach. The underestimation-biased prediction pattern for high torsion ligands is an only ligand-dependent pattern, and the correction strategy would be useful regardless of target proteins. To relieve the underestimation-biased deviation at higher torsion numbers, we have used the strategy using a modified factor (α), which can weaken the torsion number dependency of torsion energy term. However, this doesn’t mean that the predicted intermolecular energy term is correct, and only the predicted torsion-free energy term is incorrect. Theoretically, the contribution of intermolecular energy term to the total binding energy is more significant as the torsion number decreases. As shown in Fig. 1(c), the deviation of the predicted affinity from the experimental values was significantly high for the ligands with low torsion numbers, and the predicted binding affinities were not biased toward overestimation or underestimation. This pattern indicates that there are significant deviations in the prediction of intermolecular energy term. On the other hand, the torsion term effect on the total binding free energy becomes more significant as the torsion number increases. Therefore, the deviation of predicted affinity for the lower torsion numbers might be caused by the incorrect estimation of torsion energy term in addition to the inaccurate evaluation of intermolecular energy term. Meanwhile, the underestimation-biased trend of higher torsion number is presumed to be contributed mainly by the torsion energy term effect. If the intermolecular term effect is more significant than the
torsion term, the deviation would have followed the non-biased pattern observed in low torsion numbers. If the underestimation-biased trend of predicted affinity for high torsion numbers is mainly from the torsion energy term as hypothesized above, the modification factor may have a physical meaning related to the torsion energy term. Based on the effective Wtor values in Fig. 3, the values for the ligands with high torsion numbers are mostly positive values and lower than the currently used Wtor, 0.29. The positive values indicate that the torsion energy term is positive, and the torsion energy before ligandbinding is higher than after ligand-binding, as understood currently. On the other hand, the lower values than 0.29 suggest that the torsion energy difference between before binding and after binding is not so more significant than expected. There are two possibilities responsible for making such an assumption. First is that the ligands with high torsion numbers are more compact than expected when it is surrounded by water molecules in the solution. The other possibility is that the ligands with high torsional number may have far more flexible conformations than expected inside the binding pocket. Of course, there is a possibility that this kind of interpretation of the physical meaning of the torsional energy term is not meaningful because many other important factors, such as protein flexibility were not considered in our study. The Wtor value of Autodock has been generally used regardless of ligand torsion number. In this study, we found that there was an underestimation-biased deviation pattern for the ligands with a relatively higher torsion number when the Wtor was used and demonstrated the idea that the use of modification factor could relieve the bias in the affinity prediction of high torsion number ligands. To apply the approach more practically, further studies may be needed. For instance, more detailed statistical analyses and optimization studies need to be performed to achieve a more optimized modification factor. In addition, the use of modification factor is not effective to low torsion ligands. For the ligands with torsion number of 412, the average of the effective Wtor values was determined to be 0.21, closer to the default Wtor value, and there was no clear deviation pattern in the prediction of binding affinities. As mentioned above, the prediction of more accurate intermolecular energy term is necessarily needed to reduce the deviation of affinity prediction for low torsion ligands. Declaration of interests None.
Please cite this article in press as: D.K. Sriramulu, et al., Effect of ligand torsion number on the AutoDock mediated prediction of protein-ligand binding affinity, J. Ind. Eng. Chem. (2019), https://doi.org/10.1016/j.jiec.2019.12.009
G Model JIEC 4893 No. of Pages 7
D.K. Sriramulu et al. / Journal of Industrial and Engineering Chemistry xxx (2019) xxx–xxx
Acknowledgements This research was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea Government (MSIT) (No. 2018R1A2B6001670). Appendix A. Supplementary data Supplementary material related to this article can be found, in the online version, at doi:https://doi.org/10.1016/j.jiec.2019.12.009. References [1] C. Bissantz, B. Kuhn, M. Stahl, J. Med. Chem. 53 (2010) 5061. [2] J.L. Medina-Franco, O. Mendez-Lucio, K. Martinez-Mayorga, Adv. Protein Chem. Struct. Biol. 96 (2014) 1. [3] Xing Du, Yi Li, Yuan-Ling Xia, Shi-Meng Ai, Jing Liang, Peng Sang, Xing-Lai Ji, S.Q. Liu, Int. J. Mol. Sci. 17 (2016) 17. [4] W.J, H.N. Motlagh, J. Li, V.J. Hilser, Nature 508 (2014) 331. [5] G. Klebe, H.J. Bohm, J. Recept. Signal Transduct. Res. 17 (1997) 459. [6] I.A. Guedes, C.S. de Magalhaes, L.E. Dardenne, Biophys. Rev. 6 (2014) 75. [7] P. Ferrara, H. Gohlke, D.J. Price, G. Klebe, C.L. Brooks III, J. Med. Chem. 47 (2004) 3032. [8] D. Plewczynski, M. Lazniewski, R. Augustyniak, K. Ginalski, J. Comput. Chem. 32 (2011) 742.
7
[9] R.D. Smith, J.B. Dunbar Jr., P.M. Ung, E.X. Esposito, C.Y. Yang, S. Wang, H.A. Carlson, J. Chem. Inf. Model 51 (2011) 2115. [10] Z. Gaieb, S. Liu, S. Gathiaka, M. Chiu, H. Yang, C. Shao, V.A. Feher, W.P. Walters, B. Kuhn, M.G. Rudolph, S.K. Burley, M.K. Gilson, R.E. Amaro, J Comput. Aided Mol. Des. 32 (2018) 1. [11] N.S. Pagadala, K. Syed, T.J, Biophys. Rev. 9 (2017) 91. [12] Z. Wang, Y. Kang, D. Li, H. Sun, X. Dong, X. Yao, L. Xu, S. Chang, Y. Li, T. Hou, J. Phys. Chem. B 122 (2018) 2544. [13] Z. Wang, H. Sun, X. Yao, D. Li, L. Xu, Y. Li, S. Tian, T. Hou, Phys. Chem. Chem. Phys. 18 (2016) 12964. [14] S. Mukherjee, T.E. Balius, R.C. Rizzo, J. Chem. Inf. Model 50 (2010) 1986. [15] G.M. Morris, R. Huey, W. Lindstrom, M.F. Sanner, R.K. Belew, D.S. Goodsell, A.J. Olson, J. Comput. Chem. 30 (2009) 2785. [16] R. Huey, G.M. Morris, A.J. Olson, D.S. Goodsell, J. Comput. Chem. 28 (2007) 1145. [17] Renxiao Wang, Xueliang Fang, A. Yipin Lu, S. Wang, J. Med. Chem. 47 (2004) 2977. [18] L.T. Gilson, M.K. Baitaluk, M. Nicola, G. Hwang, L. Chong, J. Nucleic Acids Res. 4 (2016) D1045. [19] tL. Hu, M.L. Benson, R.D. Smith, M.G. Lerner, H. Carlson, Proteins 60 (2005) 333. [20] H.M. Berman, J. Westbrook, Z. Feng, G. Gilliland, T.N. Bhat, H. Weissig, I.N. Shindyalov, P.E. Bourne, Nucleic Acids Res. 28 (2000) 235. [21] R.A. Laskowski, E.G. Hutchinson, A.D. Michie, A.C. Wallace, M.L. Jones, J.M. Thornton, Trends Biochem. Sci. 22 (1997) 488. [22] W. Tian, C. Chen, X. Lei, J. Zhao, J. Liang, Nucleic Acids Res. 46 (2018) W363. [23] D.S.G. Garrett, M. Morris, Ruth Huey, Arthur J. Olson, (1996).
Please cite this article in press as: D.K. Sriramulu, et al., Effect of ligand torsion number on the AutoDock mediated prediction of protein-ligand binding affinity, J. Ind. Eng. Chem. (2019), https://doi.org/10.1016/j.jiec.2019.12.009