Importance of vdW and long-range exchange interactions to DFT-predicted docking energies between plumbagin and cyclodextrins

TL;DR

This study highlights the importance of van der Waals and long-range exchange interactions in density functional theory calculations for accurately predicting docking energies between plumbagin and cyclodextrins.

Contribution

It demonstrates that incorporating GD3 dispersion correction and sufficient long-range exchange enhances DFT prediction reliability for molecular docking energies.

Findings

GD3 dispersion correction improves DFT accuracy

Long-range exchange is crucial for reliable predictions

Semi-empirical methods qualitatively reproduce docking stabilization

Abstract

We calculated the docking energies between plumbagin and cyclodextrins, using density functional theory (DFT) with several functionals and some semi-empirical methods. Our DFT results revealed that GD3 dispersion force correction significantly improves the reliability of prediction. Also sufficient amount of long-range exchange is important to make it reliable further, agreeing with the previous work on argon dimer. In the semi-empirical methods, PM6 and PM7 qualitatively reproduce the stabilization by docking , yet under- and over-estimating the docking energies by ~10 kcal/mol, respectively.

Tables

Table 1. Table 1: List of exchange-correlation functionals we tested. We examined the reliability of common functionals, B3LYP, M06L, and M06-2X, and their relatives with GD3 and/or CAM corrections for predicting the docking energies between plumbagin and cyclodextrins.

Plain	GD3	CAM	CAM+GD3
B3LYP	CAM-B3LYP	B3LYP-GD3	CAM-B3LYP-GD3
M06L	M06L-GD3	–	–
M06-2X	M06-2X-GD3	–	–

Full text

\pubvolume

xx \issuenum1 \articlenumber5

\historyReceived: date; Accepted: date; Published: date \Title Importance of vdW and long-range exchange interactions to DFT-predicted docking energies between plumbagin and cyclodextrins

\Author Tom Ichibha 1,**, Ornin Srihakulung1,2,3*, Guo Chao1, Adie Tri Hanindriyo4, Luckhana Lawtrakul2, Kenta Hongo5,6,7,8, and Ryo Maezono1,8, \AuthorNames Tom Ichibha, Ornin Srihakulung, Guo Chao, Adie Tri Hanindriyo, Luckhana Lawtrakul, Kenta Hongo, and Ryo Maezono

\corresCorrespondence: ichibha@icloud.com

1 Introduction

Biopharmaceuticals are manufactured, extracted, and semi-synthesized from biological sources. Compared to chemically synthesized pharmaceuticals, they have high potency at low dose and possess distinctive medical properties.Mitragotri et al. (2014); George and Abraham (2006) On the other hand, they often lack physical and chemical stability, which poses problems for long-lasting storage, oral ingestion, and so on. One of the most promising ways to solve the problem is host-guest docking technology.George and Abraham (2006); Mitragotri et al. (2014) The biopharmaceutical molecules (guest) are combined with the carrier molecules (host) and are stabilized both physically and chemically. It also allows one to control where and how fast the biopharmaceutical is released from the carrier and absorbed in human/animal body.Vyas et al. (2008)

The binding energy of the docking is of great relevance to the stability and the release rate/timing of the biopharmaceuticals. Density functional theory is expected to be able to provide theoretical prediction with high reliability beyond those of semi-empirical methods.Christensen et al. (2016) Yet, even for DFT, it is difficult to make a reliable energy estimation since multiple non-covalent forces such as hydrogen bond or dispersion force are complexly intertwined to realize the docking interaction.Ye et al. (2016) In other words, special treatments for the non-covalent forces are needed. One example is long-range correction,Iikura et al. (2001); Yanai et al. (2004) which enhances the proportion of the exact exchange term for long-range interactions, which can improve the description of van der Waals forces.Tsuneda (2012) Another example is using a family of Minnesota functionals, whose parameter training-set includes weak-force interaction systems and has an improved accuracy to describe such systems.Zhao and Truhlar (2007) However, they cannot describe the asymptotic decline of van der Waals forces in proportion to $R^{-6}$ ($R$: inter-atomic distance), since they do not explicitly contain dispersion interactions by its construction and also their internal parameters are optimized only around equilibrium geometries.Hongo and Maezono (2017) This decline can be described by Grimmes’s dispersion correction (GD3),Grimme et al. (2010) which adds an empirical function akin to the Lennard-Jones potentialsGrimme (2011) and systematically improves the description of van der Waals systems. Grimme et al. (2010); Burns et al. (2011)

For several systems, especially for host-guest docking of biopharmaceutical compounds, system size is often an insurmountable obstacle for conventional Kohn-Sham(KS)-DFT applications. On the other hand, the recent progress of DFT out of the conventional KS-DFT has shown promise to realize protein-scale DFT calculation. Romero-Muñiz et al. (2018); Mohr et al. (2015) Therefore, it is a significantly meaningful task to study what exchange-correlation functional works for the host-guest docking in order to prepare for similar large scale calculations.

In this work, we tested eight functionals listed in Table 1, targeting the binding energy between plumbagin and cyclodextrins. This system is a representative example of the host-guest docking of biopharmaceutical compounds, providing a reasonably-sized system for the conventional KS-DFT calculation. We concluded the co-existence of vdW correction and sufficient amount of long-range exchange is essential for having a quantitatively reliable prediction. In addition to DFT, we applied semi-empirical methods PM3, PM6, and PM7,Christensen et al. (2016) and found PM6 and PM7 reproduce the stabilization of docking qualitatively. Nevertheless, they under- and over-estimate the docking energy by $\sim$10 kcal/mol, compared to the best DFT predictions.

2 System

Plumbagin is an organic molecule including two benzene rings,Oommen et al. (1999) which is reported to be effective against prostate cancer.Aziz et al. (2008); Abedinpour et al. (2017) However, this molecule cannot exist for a long time under normal atmospheric conditions as oxidation and degradation can cause losses of up to 63.8% in a single month. Suthanurak et al. (2010) Docking plumbagin molecules within cyclodextrin has been considered as an effective method to extend its short shelf-life.

Cyclodextrin (CD) is a circular molecule of glucose units as shown in Figure 1. It is broadly used as a carrier of pharmaceuticals due to the following merits (other than stabilization): Vyas et al. (2008)

• •

The ring size is adjustable to the size of guest molecule by changing the number of glucose units $n$ ($\geq$6). For $n$ equaling 6, 7, and 8, it is called $\alpha$-, $\beta$-, and $\gamma$-CD, respectively.

• •

Docking with CD improves drug solubility or dissolution, which is essential for drugs with poor water solubility.

• •

The release rate/timing is controllable by assigning different functional groups at $R_{1}$ and $R_{2}$ in Figure 1.

We selected $\beta$-CD (BCD) as the host molecule, since it has been experimentally used for docking with plumbagin,Oommen et al. (1999) and we obtained as well the docking energies for plain BCD and its two variants, Methyl-BCD (MBCD) and Hydroxy Propyl-BCD (HPBCD), shown in Figure 1.

3 Methods

We obtained the docking structures between plumbagin and BCDs from docking analysis with Lamarckian algorithm using AutoDock 4.2.6.Morris et al. (1998) This algorithm is often used to predict the ligand arrangements of protein systems. A set of genes representing the ligand arrangements are optimized to get energetically stable structures. Each of the genes consists of translations, orientations, and conformations of the ligands. We optimized the arrangement of plumbagin as the ‘ligand’ of BCDs.

The molecular structures of plumbagin and BCDs are taken from the entries: PVVAQS01, BCDEXD03, BOYFOK04, and KOYYUS in the Cambridge Structural Database.Groom et al. (2016) We optimized their structures using semi-empirical PM7 before doing docking analysis. In the docking analysis, the translation of plumbagin is discretized on a 50$\times$38$\times$24 grid with spacings of 0.375 $\AA$. We run 100 iteration for 150 of initial population of genes. At the end of each iteration, we selected only one gene with the lowest energy to "survive" to the next iteration. The energies were calculated with empirical force field, where electrostatic interaction was given based on Gasteriger charges. Gasteiger and Marsili (1978) The other input parameters were set to be the default values of Autodock 4.2.6.

We performed DFT calculations using Gaussian09/16. Frisch et al. ; Frisch et al. (2016) We used the 6-31G++($d$,$p$) basis set, since a family of 6-31G basis sets are often used for docking systems. Bachrach (2013); Bachrach and Andrews (2014); Ye et al. (2016); Bouhadiba et al. (2017); Deka and Bhattacharyya (2017) We corrected the basis set superposition error by the counterpoise method.Simon et al. (1996) In addition to the common functionals B3LYP, M06-2X, and M06L, we compared the reliability of CAM-B3LYP with long-range exchange correctionYanai et al. (2004) and B3LYP-GD3, M06-2X-GD3, M06L-GD3, and CAM-B3LYP-GD3 with Grimme’s dispersion correction,Grimme et al. (2010) as shown in Table 1. GD3 correction introduces a pair-wise function akin to the Lennard-Jones potential to the original functional. Thus, the corrected functional is no longer within the framework of density functionals.

Semi-empirical methods are generally used to predict the structures of large-size systems represented by proteins.Christensen et al. (2016) Its cost and reliability is lower than DFT but higher than the methods based on practical potentials. This method performs one-electron integrals only, so the scaling of the calculation cost is kept from rising above $\mathcal{O}(N)$ ($N$: number of electrons).Christensen et al. (2016) We employed three kinds of semi-empirical methodsChristensen et al. (2016) and calculated the docking energies using Gaussian09/16.Frisch et al. ; Frisch et al. (2016) These methods are in the orders of PM3, PM6, and PM7: PM6 improves on PM3 by refining the core-core interaction term and introducing the $d$-type basis function.Christensen et al. (2016) PM7 method further improves on this by introducing corrections for dispersion and hydrogen bonding.Stewart (2013)

4 Results and Discussions

We found that the structures obtained by the docking analysis are classified into two types of conformations as shown in Figure 2. In type-I (II), the hydroxyl phenolic (methyl quinone) group of plumbagin is placed around narrow-side of the cavity in BCDs. Thus, we took both types of conformations for our benchmark.

Figure 3 shows the predictions of the docking energies between plumbagin and BCDs for type-I and II by DFT with the functionals listed in Table 1 and the semi-empirical methods. First, in comparing the DFT results from functionals without GD3 correction, we observe that M06L and M06-2X reproduce the stabilization by docking qualitatively, while B3LYP and CAM-B3LYP do not. This is recognized to be originating from the artifact of (CAM-)B3LYP, which are optimized just for covalent systems, resulting in a poor description for non-covalent forces.Marom et al. (2011) Similar problems have been observed for the stacking of b-DNAKen et al. (2018) and the bonding of NiPc dimer.Marom et al. (2011)

DFT results for functionals with GD3 correction show that all functionals, with the notable exception of M06L-GD3, give similar predictions. It is also most surprising that the predictions from CAM-B3LYP-GD3 and M06-2X-GD3 show almost the same tendencies (see Figure 3b). It is unlikely to be a product of coincidence that two functionals based on different design concepts accidentally converge to the same results. Rather, it seems natural to conclude that the two functionals predicted very close values to the true docking energy, and resulting in the coincidence. Here, the next important question is why these functionals work better than M06L-GD3 and B3LYP-GD3. This would be attributed to their having larger amount of exact exchange used to describe long-range interactions: Kamiya et al. established in their benchmark calculations on the bonding of argon dimer that the balanced evaluation of van der Waals correlation and long-range exchange is significant to reliably describe the van der Waals interaction.Kamiya et al. (2002) They observed as well that argon dimer is over-bound when the long-range exchange is lacking.Kamiya et al. (2002) It is in agreement with our results, which shows M06L-GD3 without exact exchange predicts much higher docking energy than the other GD3 functionals.

Finally, we discuss here the results given by PM3, PM6, and PM7 semi-empirical methods. First, looking at PM3 and PM6, only PM6 reproduces the stabilization by docking with BCDs. This would be attributed to the refinement of the core-core interaction in PM6 compared to PM3. PM7 also reproduces the stabilization for all six patterns. However, we observe that the binding energies predicted by PM6 (PM7) are smaller (larger) by $\sim$10 kcal/mol than those produced from the more rigorous ab initio methods (DFT with CAM-B3LYP-GD3 and M06-2X-GD3). The failures of PM6 and PM7 are due to not explicitly containing the dispersion force correction and lacking exact exchange, respectively.

5 Conclusion

We have investigated various types of functionals which give a reliable estimation of the binding energy of docking between plumbagin and BCD, MBCD, and HPBCD, which are representative host-guest docking systems in the biopharmaceutical field. Comparing the predictions of non-GD3 functionals, we find functionals M06L and M06-2X reproduce the stabilization by host-guest docking while B3LYP and CAM-B3LYP do not. This is due to B3LYP and CAM-B3LYP having been optimized just for covalent systems. Among functionals with GD3 corrections, we concluded that CAM-B3LYP-GD3 and M06-2X-GD3 predict the binding energies very reliably, since they give surprisingly similar predictions, which cannot be considered to be just accidental. This would be attributed to both functionals possess sufficient amount of long-range exchange to properly describe non-covalent forces. Lastly, from the semi-empirical methods, PM6 and PM7 reproduce the stabilization by host-guest docking. However, each method under- and over-estimates the binding energy by $\sim$10kcal/mol, respectively.

6 Acknowledgments

The computation in this work has been performed using the facilities of the Research Center for Advanced Computing Infrastructure (RCACI) at JAIST. T.I. is grateful for financial suport from Grant-in-Aid for JSPS Research Fellow (18J12653). O.S. is grateful for financially supported by SIIT-JAIST Dual degree scholarship from Thailand’s National Electronics and Computer Technology Center (NECTEC), Sirindhorn International Institute of Technology (SIIT) and Japan Advanced Institute of Science and Technology (JAIST). K.H. is grateful for financial support from a KAKENHI grant (JP17K17762), a Grant-in-Aid for Scientific Research on Innovative Areas “Mixed Anion” project (JP16H06439) from MEXT, PRESTO (JPMJPR16NA) and the Materials research by Information Integration Initiative (MI2I) project of the Support Program for Starting Up Innovation Hub from Japan Science and Technology Agency (JST). R.M. is grateful for financial supports from MEXT-KAKENHI (17H05478 and 16KK0097), from Toyota Motor Corporation, from I-O DATA Foundation, and from the Air Force Office of Scientific Research (AFOSR-AOARD/FA2386-17-1-4049). R.M. and K.H. are also grateful to financial supports from MEXT-FLAGSHIP2020 (hp170269, hp170220).

7 Author Contributions

Data curation, Tom Ichibha, Ornin Srihakulung, Guo Chao, and Ryo Maezono; Investigation, Tom Ichibha, Ornin Srihakulung and Guo Chao; Supervision, Luckhana Lawtrakul, Kenta Hongo and Ryo Maezono; Validation, Tom Ichibha, Ornin Srihakulung, Guo Chao and Adie Tri Hanindriyo; Visualization, Ornin Srihakulung; Writing - original draft, Tom Ichibha; Writing - review & editing, Tom Ichibha, Adie Tri Hanindriyo, Luckhana Lawtrakul, Kenta Hongo and Ryo Maezono