Introduction

Tremendous advances in technological capabilities have enabled computational approaches to be applied to discern a broad range of physical, chemical, and biological phenomena across scales in molecular science [1,2,3,4,5,6]. With emphasis on molecular design, computational approaches have found great utility towards innovation in drug discovery. Considering the time and cost of the drug pipeline, from the discovery process to market, in silico biophysical methods serve an important role in expediting and reducing the cost of the discovery process, facilitating the identification, optimization, and refinement of potential drug candidates and providing comprehensive insight into the mechanism of action and structure–property relationships at the atomic level that are ultimately critical to a drug’s efficacy [7,8,9,10,11,12].

In computational strategies towards structure-based design, an important step is the prediction of probable conformations of a ligand bound to the host. To identify better possible candidate binding modes, they can be ranked via scoring functions and further evaluated via molecular simulation and free energy calculations. From free energy calculations, selectivity profiles may be constructed not only to determine binding affinities but also to provide understanding into how the ligand recognizes its host.

Because of the complexity that occurs in ligand-bound protein systems, relatively smaller representative models such as polymer-based host–guest systems are used to assess free energy methods [13,14,15,16,17,18]. Although host structures selected to represent proteins are typically much smaller than proteins, they are large enough to possess a cavity or binding pocket that allows non-covalent binding of multiple guest molecules. The advantage of using host–guest systems for assessing free energy methods is that they tend to be more rigid and symmetric than proteins, which results in fewer conformations that need to be sampled [19,20,21,22,23]. Even in the representation of proteins by more simplistic models, modeling binding free energies for these smaller models is challenging since no clear “best” computational chemistry approach has been identified; efforts are needed to better resolve strategies towards predictions of binding free energies. Statistical Assessment of the Modeling of Proteins and Ligands (SAMPL) blind challenges provide a unique platform to validate available methods and stimulate the development of new methods for quantitative predictions [13, 16, 18, 24,25,26]. In these challenges, binding affinities and other physicochemical properties are predicted, using computational models without the benefit of insight from experiment; they are then later compared to unpublished experimental measurements that allow the comparison of different computational prediction methods.

For the prediction of free energies, there are many methods available such as free energy perturbation [27], replica exchange free energy perturbation [28], and thermodynamic integration [29]. However, these methods are computationally demanding (memory, disk space, CPU time) as they converge poorly for large systems and require dividing the model into multiple intermediate steps or multiple configurations for accurate predictions [30,31,32]. In contrast, end state free energy methods are path-independent and do not require sampling of multiple configurations. These methods offer a simpler, computationally less costly approach to predict the binding free energy [33,34,35,36,37,38]. Moreover, implicit solvation models can be used to further reduce the computational cost while predicting binding free energies similar to these more sophisticated and computationally demanding methods.

While classical molecular dynamics (MD) methods are commonly used to investigate host–guest interactions, molecular mechanics (MM) force fields result in a limited treatment of effects resulting from polarization, charge transfer, and many body effects which can impact the description of properties such as binding free energies [9, 39,40,41,42,43]. To better account for these effects, quantum mechanical (QM) approaches, which are more costly, are commonly used in drug discovery research [9, 44], and have been used in previous SAMPL competitions [17, 45,46,47]. For example, in the SAMPL5 competition for host–guest binding, Caldarau et al. [45] used density functional theory (DFT) with an added dispersion correction (DFT-D3) and the wavefunction-based domain-based local pair-natural orbital coupled cluster (DLPNO-CCSD(T)) method to predict the binding energies for octa-acid (OA) host–guest systems. In this approach, they used TPSS-D3/def2-SVP optimized structures and host structures are constrained during MD simulations to reduce the flexibility of the host and limit the structural distortions resulting from the repulsion between the negative charge of the ligands and the large negative charge of the OA hosts. This approach yielded binding energies approximately 12.0 kcal/mol greater than the experimental binding affinities, with a low correlation coefficient (r2 ≈ 0), and a statistically insignificant Kendall’s rank correlation coefficient (τ ≤ 0.20) for all attempts for the host–guest systems in the SAMPL5 blind challenge due to incorrect representative structures, not sampling enough conformational binding positions for ligands, and thermochemical corrections that yielded up to a 7.2 kcal/mol difference depending on the method of choice. This performance demonstrates the limited sampling capabilities of current QM methods compared to MD methods, obtained representative structures, as well as thermodynamic and solvation corrections.

Contrary to this, in the SAMPL4 competition for host–guest binding, Mikulskis et al. [47] were successful with both MM- and QM-based approaches for OA hosts with mean absolute deviations (MADs) less than 2.0 kcal/mol. Their MM approach, which utilized free energy perturbation (FEP) calculations, yielded MADs of approximately 1.0 kcal/mol while their QM approaches with DFT-D3 optimized structures yielded MADs of approximately 1.0–2.0 kcal/mol depending on the implementation of a solvent in the calculations, i.e. no solvent, implicit solvent, or a combined implicit-explicit solvent. However, the combination of FEP and DFT-D3 did not yield favorable results due to the large difference between the MM and DFT potential energy functions. Sure et al. [46] provided another successful attempt at using DFT-D3 for the SAMPL4 competition for host–guest binding of a macrocyclic cucurbit[7]uril host by optimizing the geometry at the TPSS-D3/def2-TZVP level of theory after pre-optimizing possible binding scenarios with the HF-3c semiempirical method. These optimizations were followed by single point calculations using PW6B95-D3/def2-QZVP with the g- and f-functions for non-hydrogen and hydrogen atoms removed, respectively, with the COSMO-RS implicit solvent model, which yielded an MAD of 2.0 ± 0.5 kcal/mol. These two studies highlight that for the SAMPL4 competition, host–guest structure optimization and higher-level MM-based approaches like FEP can be vital in characterizing correct binding interactions at the QM level.

In this work, efforts in MD and QM methods are combined to predict binding affinities for fourteen ligands to a macrocyclic cucurbit[8]uril host [19, 21, 22, 48] and eight ligands to two variants of the OA deep-cavity cavitands [20, 23]. Using MD simulations to obtain representative structures, MM- and QM-based methods are utilized to predict binding free energies. Within the QM methods, the use of a resolution-of-the-identity (RI) approximation designed for larger molecules [49], Grimme’s D3 atom-pairwise dispersion corrections with Becke-Johnson damping [50], and truncated correlation consistent basis sets for the hydrogen atoms [51] are evaluated to probe how different electronic structure approaches that reduce the computational cost contribute to predicting binding affinities. Insights into what strategies are more favorable for host guest-binding will help to build a framework for predicting host–guest binding affinities using QM approaches.

Methods

System preparation

The initial structures for the guest molecules, shown in Figs. 1 and 2, and the three host molecules, shown in Fig. 3, cucurbit[8]uril (CB8), octa-acid (OA), and tetramethyl octa-acid (TEMOA), that were issued with the SAMPL6 challenge dataset were used to generate the host–guest systems. The CB8 molecule has no formal charge whereas the octa-acids (OA/TEMOA) have eight deprotonated carboxylic acid groups and thus a formal charge of − 8. Even though OA and TEMOA are water-soluble structurally similar deep-cavity cavitands, the TEMOA host has four methyl groups in place of four hydrogen atoms present in the OA host located on the upper rim of the cavitand that enclose the hydrophobic binding pocket.

Fig. 1
figure 1

The guest molecules for the cucurbit[8]uril (CB8)

Full size image
Fig. 2
figure 2

The guest molecules for the octa-acid (OA) and tetramethyl octa-acid (TEMOA) hosts

Full size image
Fig. 3
figure 3

The host molecules: cucurbit[8]uril (CB8), octa-acid (OA), and tetramethyl octa-acid (TEMOA)

Full size image

Initial binding poses of guest molecules binding to the host were generated using the docking feature implemented in Molecular Operating Environment (MOE) v2016.08 [52]. The London ΔG scoring function [53] was used to estimate ligand placement in the pocket. The top 100 poses given by the London ΔG scoring function for each host–guest complex were refined to a list of ten poses by rescoring the flexible receptor and ligand conformation using the GBVI/WSA ΔG scoring algorithm [52, 53]. Among these ten poses, those with minor structural differences were discarded from the list of ten poses and the chemically relevant poses with the highest GBVI/WSA ΔG scores were selected for further investigation. The chosen host–guest poses were minimized under the AMBER10: Extended Hückel Theory (EHT) potential implemented in MOE, which employs Amber ff10 and EHT bonded parameters [54,55,56].

To generate force field parameters, the AM1-BCC scheme [57] was used for generating partial charges for the guest and host molecules using the Antechamber suite. As the guest molecule CB8-G13 contains a platinum atom, the Mulliken charges were calculated from a geometry optimization using B3LYP [58, 59] in conjunction with the 6-31G(d) basis set [60] and the effective core potential basis set Lan2L2DZ [61] for the platinum atom. All electronic structure calculations were performed in Gaussian 16 [62]. The host–guest systems were further prepared for simulation using the Leap module of AmberTools  [63] under the General Amber Force Field (GAFF) [56]. Each system was neutralized with counterions using parameters from Joung and Cheatham [64], and solvated in a 14.0 Å cube of TIP4P-Ew water [65] beyond the solute. To mimic the ionic strength of the experimental buffers, additional counter ions were added to create a buffer of 150 mM sodium chloride for the CB8 complexes and 60 mM sodium chloride for the OA/TEMOA complexes.

Simulation protocol

The host–guest systems were relaxed using NVT ensembles over six minimization procedures with decreasing restraints on the host of 500.0, 200.0, 20.0, 10.0, 5.0 kcal/mol (Å2), and then were heated to 300 K over 30 ps. The temperature was maintained at 300 K using Langevin dynamics and the pressure was coupled to 1 atm using isotropic position scaling. Atomistic molecular dynamics simulations were performed for 10 ns in triplicate to account for randomized parameters that affect the MD trajectories. Nonbonded interactions were truncated with a 10.0 Å cutoff, whereas long-range electrostatics were handled with the particle-mesh Ewald (PME) method. Bonds involving hydrogen were constrained using SHAKE, and the simulation time step was set to 2 fs. All simulations were performed with AMBER16.7 [63].

The binding free energies were calculated with MMPBSA approach using the built-in PBSA-solver [66]. The internal and external dielectric constants were set to 1.0 and 80.0, respectively. The solvent accessible surface area (SASA) was determined with the default LCPO method using the modified Bondi atomic radii. Calculations for solute entropic contributions were not considered. For each system, the binding free energy was determined using the final 100 frames from the simulation.

Clusters were formed using the density-based spatial clustering of applications with noise (DBSCAN) algorithm based on two parameters, which are epsilon (Eps) and the minimum number of points in an Eps-neighborhood (MinPts). MinPts was set to 4 and the Eps value for DBSCAN was determined from the threshold point of a sorted 4-dist graph [67]. The cluster conformation representing the greatest number of frames from the MD simulations was used for further analyses. Additional QSAR (quantitative structure–activity relationship) calculations were performed on each guest molecule to determine the van der Waals volume each molecule occupies by using the connection table approximation descriptor in MOE (see Tables S1 and S2 in supplemental information).

Quantum mechanical methods

The individual structures generated from the clustering of MD trajectories, shown in Figs. 4, 5 and 6, for each host–guest complex were used for all quantum chemical calculations. The host and guest molecules were analyzed with the same geometry as from the complex. The thermal corrections for all molecules were calculated at the HF/6-31G(d) level of theory in Gaussian 16 and the vibrational contributions were scaled by 0.8953 [68]. Single point energies were obtained using ORCA 4.0 [69] with the B3PW91 density functional [58, 70, 71] since B3PW91 has been shown to properly treat long-range covalent interactions. In the treatment of the exact exchange in the functional, the RIJCOSX approximation [49] was used with the def2 auxiliary basis set [72] to reduce the computational cost associated with the number of atoms in the host–guest complex since the RIJCOSX approximation has been shown to be five times as efficient for molecules of similar size to the host–guest systems. To mimic the aqueous solution, the SMD implicit solvation model [73] was used with water (ε = 78.4) as the implicit solvent. Grimme’s D3 dispersion correction with Becke-Johnson damping was used to investigate long-range covalent interactions as the inclusion of D3 dispersion improves intermolecular interaction energies predicted with DFT [46, 50, 74, 75].

Fig. 4
figure 4

The structures of the CB8 guest molecules inside the binding pocket. These structures are generated from the clustering analysis

Full size image
Fig. 5
figure 5

The structures of the OA guest molecules inside the binding pocket. These structures are generated from the clustering analysis

Full size image
Fig. 6
figure 6

The structures of the TEMOA guest molecules inside the binding pocket. These structures are generated from the clustering analysis

Full size image
Table 1 The binding free energies in kcal/mol for the CB8 host–guest systems
Full size table
Table 2 The binding free energies in kcal/mol for the OA host–guest systems
Full size table
Table 3 The binding free energies in kcal/mol for the TEMOA host–guest systems
Full size table

The correlation consistent basis set family (cc-pVnZ) [76] was used for all single point calculations since these basis sets were developed to exhibit convergence behavior to the complete basis set (CBS) limit for wavefunction-based methods through extrapolation [77,78,79,80]. Knowing the CBS limit, which removes basis set incompleteness error, the error for the property of interest, i.e. binding free energy, only corresponds to the intrinsic error of the chosen QM method. Therefore, to extrapolate to the Kohn–Sham limit for DFT methods, analogous to the CBS limit for wavefunction-based methods, the cc-pVnZ basis sets were used (n = D, T) with the following extrapolation scheme proposed by Jensen

$$E\left( {{l_{\hbox{max} }}} \right)={E_{CBS}}+A\left( {{l_{\hbox{max} }}+1} \right){e^{ - B\sqrt {{n_s}} }}$$
(1)

where lmax is the maximum angular momentum function in the basis set and ns is the number of s functions in the basis set [81]. The B-parameter was set to 5.5 in agreement with Jensen for use as a two-point extrapolation scheme. Due to the abundance of weak molecular interactions in biomolecules, the calculated binding energies were counterpoise-corrected before the extrapolations were performed on each host, guest, and host–guest complex [82, 83].

Additional electronic structure modeling techniques were applied to the CB8 host–guest systems to examine the impact of various approximations on the binding free energy. Targeting reduction in computational time, the correlation consistent basis sets were truncated via the removal of higher angular momentum basis functions for hydrogen atoms. This has been shown to reduce the computational time by approximately 42.9% and 57.8% when removing 1 d function from the cc-pVTZ basis set, denoted as cc-pVTZ(− 1d), and 2 d functions and 1 f function from the cc-pVQZ basis set, denoted as cc-pVQZ(− 1f2d), respectively, and yielded the results closest to the atomization energies generated with the full basis sets at the complete basis set limit [51].

Binding free energies calculated with and without the use of the resolution-of-the-identity (RI) approximation were examined to gauge how the RI approximation, which leads to a reduction in CPU time, affects the accuracy. To characterize the ionic strength of the solution used in experiment, the dielectric constant for the implicit water solvent was also altered from 78.4 for pure water to 76.4 given the concentration of the sodium chloride solution used in the MD simulations and the experimentally determined relation between the concentration of an ionic solution and the dielectric constant [84].

Results

The binding free energies submitted as part of the SAMPL6 competition are shown in Tables 1, 2 and 3 for CB8, OA, and TEMOA host–guest systems, respectively. For each host–guest complex, statistical measurements were used to gauge the effectiveness of each of the three methods, which are MMPBSA, RI-B3PW91-D3, and RI-B3PW91, in predicting experimental binding free energies. These include the mean absolute error (MAE), the root mean square error (RMSE), Kendall’s Tau (τ) rank correlation coefficient, which measures how well a method ranked calculated binding free energies relative to experimental binding free energies where τ values closer to one correspond to increased qualitative accuracy of the prediction, and the correlation coefficient (r2). To demonstrate there is no correlation in ranking between the calculated binding free energies and the experimental binding free energies, τ values are compared against τcrit, a cutoff value obtained through a table of critical values generated by Monte Carlo simulations of a τ distribution, which is similar to the normal Z distribution, used to reject the null hypothesis [85, 86].

Cucurbit[8]uril (CB8)

The binding free energy predictions for the CB8 host with the three methods submitted were compared to experiment (Table 1). The predicted values were significantly more negative than experimental binding free energies with an MAE of 16.69, 33.58, and 15.54 kcal/mol for MMPBSA, RI-B3PW91-D3, and RI-B3PW91, respectively.

When the binding affinities of the guests to CB8 are ranked from the lowest to the highest binding affinity, MMPBSA did not correctly rank any of the systems but predicted CB8-G12 to have a stronger binding affinity relative to the other complexes, which correlates to experiment well. RI-B3PW91-D3 correctly ranked CB8-G2 as the tenth strongest bound host–guest complex and predicted that CB8-G12 was more tightly bound relative to the other CB8 host–guest systems. RI-B3PW91 correctly ranked CB8-G6, CB8-G2, CB8-G1, and CB8-G3 as fifth, tenth, eleventh, and fourteenth, respectively, while the remaining systems were ranked incorrectly. Unlike both MMPBSA and RI-B3PW91-D3, RI-B3PW91 predicted CB8-G12 to have a lower binding affinity relative to the other CB8 host–guest systems.

Octa-acid (OA)

The three sets of submitted binding free energy predictions for OA are reported in Table 2. All values predicted using MMPBSA were significantly more negative than experimental measurements with an MAE of 6.8 ± 0.2 kcal/mol. When ranking the binding affinities of the guest to the host from lowest to highest binding affinity, MMPBSA correctly placed OA-G2, OA-G4, OA-G6, OA-G5 as first, second, sixth, and eighth, respectively. The other systems were not ranked correctly; OA-G0, OA-G1, OA-G7 and OA-G3 ranked third, fourth, fifth, and seventh, respectively, whereas experimentally ranked fourth, seventh, third, and fifth, respectively.

For RI-B3PW91-D3 and RI-B3PW91, the binding free energy predicted for OA-G2 was determined as a statistical outlier with 99% confidence, visualized in Fig. 8, using Dixon’s Q-Test [87]. When the statistical outlier (OA-G2) was excluded from the RI-B3PW91-D3 set, the MAE, RMSE, Kendall’s Tau (τ), and the correlation coefficient (r2) increased from 35.46 to 38.39 kcal/mol, 36.41 to 38.52 kcal/mol, 0.29 to 0.71, and 0.44 to 0.52, respectively. When the binding free energy for OA-G2 was excluded from the set of binding free energies obtained with RI-B3PW91, the MAE, RMSE, and r2 decreased from 17.87 to 12.85 kcal/mol, 22.51 to 13.39 kcal/mol, and 0.60 to 0.03, respectively, as shown in Table 2. In Fig. 7b, the statistical outlier was removed, which improved and worsened the linear regression model comparing experiment to RI-B3PW91-D3 and RI-B3PW91, respectively. With the exclusion of OA-G2, ranking the binding affinities from lowest to highest, RI-B3PW91-D3 correctly ranked OA-G4, OA-G1, and OA-G5, as first, sixth, and seventh, respectively, while RI-B3PW91 did not correctly ranked any of the systems.

Fig. 7
figure 7

Plots for calculated results in Tables 1, 2 and 3 versus experimental results in kcal/mol for (a) CB8, (b) OA, and (c) TEMOA for MMPBSA (blue), RI-B3PW91-D3 (black), and RI-B3PW91 (green). The dashed lines in each corresponding color refers to the best fit line where the statistical outlier (OA-G2) for RI-B3PW91 and RI-B3PW91-D3 is removed for b and c. The dashed gray line is the y = x line

Full size image
Fig. 8
figure 8

Error plots from experimental results in kcal/mol for (a) CB8 (b) OA, and (c) TEMOA for MMPBSA (blue), RI-B3PW91-D3 (black), and RI-B3PW91 (green) for the submitted results from Tables 1, 2 and 3

Full size image

Tetramethyl octa-acid (TEMOA)

TEMOA is structurally different from OA because of the substitution of four hydrogens around the portal to the binding pocket of OA with four methyl groups. While the same guests bound to TEMOA and OA with similar binding energies, G7 weakly binds to TEMOA relative to the other guests whereas it binds stronger to OA experimentally. Binding free energy predictions using the submitted methods for the TEMOA host are reported in Table 3. Similar to OA, all three methods overestimated the binding free energies relative to experiment. RI-B3PW91-D3 overestimated the binding free energies with an MAE of 38.83 kcal/mol. Of the three methods considered, the MMPBSA method yielded better binding free energies, both quantitatively (MAE of 5.9 ± 0.2 kcal/mol) and qualitatively (τ = 0.79), than the QM-based calculations. MMPBSA ranked TEMOA-G0 and TEMOA-G1 as the third and fourth strongest bound complexes, respectively. Additionally, MMPBSA predicted that TEMOA-G4 and TEMOA-G2 were the most tightly bound complexes while TEMOA-G7 and TEMOA-G5 were the most loosely bound complexes. RI-B3PW91-D3 correctly predicted that TEMOA-G4, TEMOA-G2, and TEMOA-G3 were the first, second, and fifth most tightly bound complexes, respectively. Like MMPBSA, RI-B3PW91-D3 predicted that TEMOA-G5 was a weakly bound host–guest complex relative to the other TEMOA host–guest systems. RI-B3PW91 correctly predicted TEMOA-G0 as the third strongest bound host–guest complex and yielded the lowest deviation from experiment (0.41 kcal/mol) for TEMOA-G2.

Quantum mechanical calculations

The CB8 host–guest systems were used to probe approaches for improving the binding free energy prediction. Specifically, the effects of (1) utilizing truncated correlation consistent basis sets as opposed to standard correlation consistent basis sets; (2) utilizing traditional DFT calculations (neglecting the RI approximation); and (3) modifying the dielectric constant used in the continuum solvation model to reflect the ionic strength of the solution used in experiment were examined.

As shown in Tables 1, 2 and 3, for CB8, OA without the statistical outlier (OA-G2), and TEMOA, the MAE, and RMSE increased by approximately 19.4, 25.5, and 32.6 kcal/mol when using Grimme’s D3 dispersion with RI-B3PW91, respectively, away from experiment. However, when using Grimme’s D3 dispersion, the τ value decreases from 0.05 to − 0.14 for CB8 but increases from − 0.05 to 0.71 when the statistical outlier is removed for OA and increases from − 0.14 to 0.57 for TEMOA. This shows the importance of using a dispersion correction for qualitative ranking of binding affinities.

The binding free energies as a result of utilizing truncated basis sets individually and extrapolated to the Kohn–Sham limit with a two-point extrapolation using cc-pVDZ and cc-pVTZ (cc-pV∞Z[D,T]) and a three-point extrapolation using cc-pVDZ and truncated triple and quadruple correlation consistent basis sets, cc-pVTZ(− 1d) and cc-pVQZ(− 1f2d), denoted as cc(0,− 1,− 2), are reported in Tables 4 and 5, respectively.

Table 4 The binding free energies for CB8 complexes in kcal/mol with various schemes involving not using the RI approximation, changing the dielectric constant of the implicit solvent with the truncated correlation consistent basis sets for hydrogen
Full size table
Table 5 The binding free energies for the CB8 complexes in kcal/mol with various schemes involving not using the RI approximation, changing the dielectric constant of the implicit solvent, and two options for basis set choice when extrapolating to the Kohn–Sham limit
Full size table

For the CB8 complexes in Table 4, using standard DFT (B3PW91-D3) yielded a MAE of 35.33 kcal/mol and 34.19 kcal/mol with cc-pVTZ(− 1d) and cc-pVQZ(− 1f2d), respectively, while RI-DFT (RI-B3PW91-D3) yielded a MAE of 34.81 and 34.18 kcal/mol for cc-pVTZ(− 1d) and cc-pVQZ(− 1f2d), respectively. When changing ε from 78.4 for pure water to 76.4 to account for the ionic strength of the solution (RI-B3PW91-D3 (ε = 76.4)), all metrics (MAE, RMSE, τ, and r2) used to gauge the method’s predictive qualities for the binding free energies did not significantly change with respect to the binding free energies predicted in pure water (RI-B3PW91-D3 (ε = 78.4)).

Table 5 shows the predicted binding free energies for B3PW91-D3 (ε = 78.4), RI-B3PW91-D3 (ε = 78.4), and RI-B3PW91-D3 (ε = 76.4) at the Kohn–Sham limit using cc-pV∞Z[D,T], a two-point extrapolation using cc-pVDZ and cc-pVTZ, and cc(0,− 1,− 2), a three-point extrapolation using cc-pVDZ, cc-pVTZ(− 1d) and cc-pVQZ(− 1f2d) for the CB8 complexes. Using the cc(0,− 1,− 2) basis set choice for extrapolation, the binding free energies predicted by RI-B3PW91-D3 (ε = 78.4) and RI-B3PW91-D3 (ε = 76.4) lowered the MAE by approximately 1.6 kcal/mol and 4.2 kcal/mol, respectively, in regards to using the cc-pV∞Z[D,T] scheme.

Discussion

Calculating end-state binding free energies with MMPBSA is relatively fast and simple but results of a loss in accuracy and reliance compared to other free energy methods. It has been known that various factors affect the performance of the MMPBSA method such as the force field, solute dielectric constant, as well as sampling [33]. In our model, we employed the AM1-BCC partial charge scheme for the guest and host molecules for use with the GAFF force field to increase computational efficiency. GAFF was designed to use partial charges calculated from the restrained electrostatic potential fit (RESP) method [56, 88]. Although, the AM1-BCC scheme was parameterized to reproduce RESP charges, this may only be appropriate for the guest molecules rather than the larger host molecules. The interactions between the host and guest molecules may have been overestimated or underestimated as a result using the AM1-BCC charge scheme, hence the binding affinity predictions may be improved by using the RESP charge model.

Submission analysis

For the methods submitted to the SAMPL6 competition, using RI-B3PW91-D3 yielded higher τ values for OA and TEMOA than using RI-B3PW91 for predicting binding free energies. Since there are eight guests that are bound to OA and TEMOA, τcrit for α = 0.05 is 0.57 for 8 data points. Only MMPBSA correlates with experiment (|τ| > τcrit), as the τ values are 0.64, 0.29, and − 0.21 for MMPBSA, RI-B3PW91-D3, and RI-B3PW91, respectively. However, after removing the statistical outlier, OA-G2, from the dataset, τ increases from 0.29 to 0.71, which implies that RI-B3PW91-D3 also correlates with experiment. As shown in Table 2, RI-B3PW91-D3 ranked the binding free energies more correctly than MMPBSA when the outlier is excluded. For TEMOA, both MMPBSA and RI-B3PW91-D3 correlate with experiment with τ values of 0.79 and 0.57, respectively, which are greater than τcrit.

As shown in Fig. 7a, there is no correlation between experimental and predicted binding free energies for the CB8 host–guest systems. This is supported by r2 ≈ 0 and τ values of − 0.19, − 0.14, 0.12 for MMPBSA, RI-B3PW91-D3, and RI-B3PW91, respectively, which are smaller in magnitude than τcrit for α = 0.05 for 14 data points, which is 0.36. This also shows an inconsistency when using Grimme’s dispersion correction, which may be due to the abundance of N and O atoms present in the CB8 host and empirical descriptors for those atoms. For all sets of the host–guest systems, RI-B3PW91 had a lower MAE and RMSE than RI-B3PW91-D3 by approximately 19.4–32.6 kcal/mol, but as a tradeoff, resulted in qualitatively better predictions of the binding affinities (Fig. 8). This implies that using a dispersion correction overbinds the guest to the host but is needed for proper ranking.

To estimate the relative performance of the methods, the mean signed error (MSE) was used to offset the calculated binding free energies. After the removal of MSE from the MMPBSA and RI-B3PW91-D3 predicted binding free energies for OA and TEMOA, the MAE and the RMSE values are recalculated to estimate the performance of methods in relative terms as shown in Table 6. This correction improved the MAE and RMSE for MMPBSA by 6.8 and 5.9 kcal/mol for OA and TEMOA, respectively. The correction improved the RI-B3PW91-D3 MAE and RMSE by 38.39 and 38.83 kcal/mol for OA without the OA-G2 outlier and TEMOA, respectively.

Table 6 The predicted binding energies for OA and TEMOA using MMPBSA and RI-B3PW91 after the removal of mean signed error (MSE)
Full size table

Impact of truncated basis sets

For the QM calculations, the subset of the CB8 host–guest systems was chosen because the size of these systems is smaller compared to the octa-acid host–guest systems investigated. While using the RI approximation, lowering ε from 78.4 for pure water to 76.4 to account for the ionic strength of the solution increased the MAE by 0.56 kcal/mol. However, altering the dielectric constant from 78.4 to 76.4 to account for the ionic strength of the solution lowered the MAE from 34.85 to 30.65 kcal/mol for the three-point extrapolation with truncated triple-ζ and quadruple-ζ correlation consistent basis sets, yet for RI-B3PW91-D3 (ε = 78.4), the MAE only decreased from 34.29 to 32.69 kcal/mol (Table 5). Therefore, factors that can change the dielectric constant should be considered when using implicit solvent models for binding free energy predictions.

The use of the cc(0,− 1,− 2) basis set scheme lowered the MAE for CB8 complexes by 1.60 kcal/mol relative to using cc-pV∞Z[D,T] (Table 5) for RI-B3PW91-D3 (ε = 78.4). In contrast, when using truncated basis sets and standard basis sets for binding free energies (Table 4), the MAE decreased by 0.51 kcal/mol for the CB8 complexes when using cc-pVTZ as opposed to cc-pVTZ(− 1d) for RI-B3PW91-D3 (ε = 78.4). The MAE decreased by 0.31 kcal/mol when increasing the basis set quality of truncated basis sets for RI-B3PW91-D3 (ε = 78.4). Therefore, within the RI approximation, the decrease in MAE when using cc-pVQZ(− 1f2d) highlights the importance of using higher quality basis sets when extrapolating to the Kohn–Sham limit.

For predictions without the RI approximation, the binding free energies determined using B3PW91-D3/cc-pVTZ yielded a decrease in the MAE by 7.65 kcal/mol relative to B3PW91-D3/cc-pVTZ(− 1d) as shown in Table 4. This is believed to be a result from including the four-center two-electron electron repulsion integrals removed via the RI approximation and the need for additional polarization when describing interactions with hydrogens between the host and the guest. This effect also contributes to the increase of 3.25 kcal/mol in the MAE between B3PW91-D3/cc-pV∞Z[D,T] and B3PW91-D3/cc(0,− 1,− 2). However, as shown in Table 5, when employing truncated basis sets (cc(0,− 1,− 2)), binding free energy predictions when using RI-B3PW91-D3 (ε = 76.4) are more positive and yield a MAE of 0.28 kcal/mol lower than B3PW91-D3 (ε = 78.4). This illustrates that within the RI approximation, changing the dielectric constant is as beneficial to predicting binding free energies as utilizing standard DFT, which is more computationally costly than RI-DFT.

For the CB8-G6 host–guest complex, which was one of the smaller systems in the set of host–guest systems, the number of basis functions decreased from 4016 to 3696 with the truncation of 1 d basis function from the cc-pVTZ basis set for hydrogen and decreased from 7640 to 6872 with the truncation of 1 f and 2 d basis functions from the cc-pVQZ basis set for hydrogen. Since DFT scales approximately N3 to N5 depending on the complexity of the functional where N is the number of basis functions, truncated basis sets become a practical option for further decreasing the computational cost while improving the quantitative prediction of binding free energies for these host–guest systems as truncating 1 d basis function from cc-pVTZ only affected the binding energy predicted with cc-pVTZ by ≤ 0.06 kcal/mol as shown in Table 4 for RI-B3PW91-D3.

Impact of the extrapolation scheme B-parameter

Another factor that can account for the large deviations between host–guest binding energies is the parameter used to fit Eq. 1 for two-point extrapolations. The value of 5.5 proposed by Jensen for the B-parameter, which was used for atoms and diatomics, caused the extrapolation curve to converge at a very rapid rate and is reflected in the predictions for the CB8 complexes, as the binding affinities in Table 1 are identical to those predicted with the cc-pVTZ basis set with the respective method in Table 4. Also, when using the three-point extrapolations with truncated basis sets for the CB8 complexes, the B-parameter yielded an average value of 0.37 (Table S3). Therefore, the value of 0.37 for the B-parameter was applied to two-point extrapolations with cc-pVDZ and cc-pVTZ to gauge how changing the B-parameter affects the extrapolated binding free energies (Table 7). The results from using 0.37 as the B-parameter in a two-point extrapolation show that the MAE decreased by 0.84 and 0.42 kcal/mol for the CB8 and TEMOA complexes, respectively. The MAE did not change for the OA complexes. Setting the B-parameter to 0.37 did not change the τ values for CB8 and OA complexes, however, did increase the τ value from 0.57 to 0.71 for TEMOA.

Table 7 The predicted binding energies when using different values for B in Eq. 1 for two-point extrapolations using cc-pVDZ and cc-pVTZ with RI-B3PW91-D3
Full size table

In addition to applying 0.37 for the B-parameter to predict binding free energies for all host–guest systems using two-point extrapolations with cc-pVDZ and cc-pVTZ, the value of the B-parameter was optimized to the value of 0.12 via minimizing the MAE and was applied (Table 7). For the CB8 host–guest systems, shifting the B-parameter from 5.5 to 0.12 had a noticeable impact on the MAE, which decreased from 34.29 to 29.84 kcal/mol for RI-BWPW91-D3. A similar effect was observed for TEMOA with a decrease in the MAE of 5.07 kcal/mol. There is no notable change in MAE, RMSE, or τ for the OA complexes with the change in the B-parameter. Furthermore, τ increases from 0.57 to 0.93 when the B-parameter is changed from 5.5 to 0.12 for TEMOA with RI-B3PW91-D3, which provides more evidence that dispersion-corrected functionals should be used for qualitative predictions of binding free energies since |τ| > τcrit. The observed trends imply that the value of the B-parameter should be reoptimized when using Eq. 1 for macromolecules.

Compared to other submissions employing QM methods in the SAMPL6 host–guest binding challenge, our approach yielded quantitatively poorer predictions that may have resulted from the approximations considered in this work. In our approach, only a single conformational state of the guest binding to the host system was considered. Additionally, the representative structures of the individual host–guest systems obtained from clustering the MD trajectories were not optimized with QM methods and is reflected in our model chemistries.

Impact of representative geometries

The cause of OA-G2 being a statistical outlier is suspected to be from the orientation of the substituted cyclohexene ring relative to the OA host (Fig. 5). Comparing OA-G2 and TEMOA-G2 in Figs. 5 and 6, where the only difference is the four methyl groups on the host, the structure of the OA-G2 complex has a smaller binding pocket than the TEMOA-G2 complex. While the experimental data suggests that G2 has a stronger binding affinity towards OA than TEMOA, MMPBSA suggests the opposite. More sampling of representative structures would aid in depicting whether the anomalous binding behavior of OA-G2 correlates with the positive binding free energies predicted with DFT.

Although the only difference between CB8-G6 and CB8-G7 was the expansion of the ring for the guest by one CH2 group, the predicted binding affinities for the CB8-G6 and CB8-G7 complexes differed by approximately 17.0 kcal/mol. This may be due to the binding poses of CB8-G6 and CB8-G7 complexes, as G6 bound in a perpendicular fashion inside the binding pocket relative to the host whereas G7 bound in a parallel fashion inside the binding pocket. This would affect nearby electrostatic interactions and why for B3PW91-D3 (ε = 78.4), RI-B3PW91-D3 (ε = 78.4), and RI-B3PW91-D3 (ε = 76.4), there was a 3.00 kcal/mol difference in the change of binding energies between CB8-G6 and CB8-G7 when improving basis set quality via the basis set scheme used for extrapolation (Table 5). Ergo, more sampling of chemically relevant structures or enhanced sampling methods can provide a more robust depiction of the host–guest binding environment.

The volumes of guest molecules for OA and CB8 molecules were compared to each other. The volumes of the guests bound to CB8 are larger than those bound to OA and TEMOA as shown in Tables S1 and S2. The guests CB8-G0, CB8-G1, CB8-G2, CB8-G3, CB8-G4, and CB8-G12 are among the largest ligands for this year’s competition with volumes of 462, 518, 432, 468, 817, and 553 Å3, respectively. These values are more than twice the average volume of OA guests and the absolute error between the experimental and the predicted binding free energies for the larger CB8 guests are among the highest for all our methods (MMPBSA, RI-B3W91-D3 and RI-B3PW91) as shown in Fig. 8. The MMPBSA and RI-B3PW91-D3 methods have a definite correlation with the experiment based on ranking the binding affinities of the octa-acid guest molecules, which were smaller in volume on average compared to the CB8 guests. This correlation is evident from the τ values of 0.64, 0.79, 0.71, 0.57 for MMPBSA (OA), MMPBSA (TEMOA), RI-B3PW91-D3 (OA without OA-G2 outlier), and RI-B3PW91-D3 (TEMOA), respectively.

However, these two methods do not correlate to the CB8 binding free energies since the τ values are − 0.19 and − 0.14 for MMPBSA and RI-B3PW91-D3, respectively. This may result from insufficient sampling as the CB8 guests are larger molecules with higher conformational flexibility. For example, the size of CB8-G4 does not allow the guest to fit entirely into the binding cavity. As a result, most of the CB8-G4 molecule is weakly bound to the host from outside of the binding pocket and only one of the three triethyl amines within the guest can fit into the pocket as shown in Fig. 4. Each triethyl amine group could bind to the host from inside the binding cavity, which would result in alternative binding conformations and affect the overall binding free energy. To better understand binding free energies of these large structures, more sampling of the different binding modes is needed to generate weighted averages based on the thermodynamic stability of predicted poses.

The results for OA and TEMOA systems illustrate that MMPBSA and RI-B3PW91-D3 methods can be used to qualitatively rank binding energies of small molecules. Among those two methods, MMPBSA is computationally less expensive, but RI-B3PW91-D3 predicted the relative binding affinities better for OA and TEMOA host–guest systems. However, the MAE and the corresponding error plots (Fig. 8) indicate that both methods overestimated the binding free energies. The MAE reported for the OA and TEMOA complexes state that MMPBSA and RI-B3PW91-D3 predict overbinding by 6.8 and 35.5 kcal/mol, respectively, for OA complexes and 5.9 and 38.8 kcal/mol, respectively, for TEMOA complexes. For all systems, the MMPBSA method was the best approach overall in terms of quantitative predictions.

Conclusions

When implementing DFT for predicting host–guest binding affinities, the use of Grimme’s D3 dispersion correction was essential for qualitatively predicting the binding free energies for the OA and TEMOA systems even though the MAE exceeded 35.0 kcal/mol for both the OA and TEMOA systems. When using implicit solvent models, factors that can change the dielectric constant, such as the ionic strength of the solution, are relevant for predicting binding free energies, as lowering the dielectric constant lowered the MAE. While RI-B3PW91-D3 reduced the computational cost relative to B3PW91-D3, B3PW91-D3 yielded a lower MAE. To attain more quantitatively favorable results, using cc-pVQZ(− 1f2d) for hydrogen atoms reduces the computational cost relative to using cc-pVQZ while simultaneously providing a better standard for extrapolating to the Kohn–Sham limit than only utilizing cc-pVDZ and cc-pVTZ for extrapolations. Also, truncating 1 d basis function for hydrogen atoms had a very small effect on predicted binding free energies obtained with cc-pVTZ, indicating that truncated basis sets are a viable option to reduce the computational cost while yielding near-identical binding free energies. With the extrapolation scheme utilized, the B-parameter should be revised for macromolecules since reducing the value of the B-parameter from the proposed 5.5 to 0.12 reduced the MAE while providing extrapolated binding energies that were in alignment with those predicted using quadruple-ζ level basis sets. Sampling of different binding poses becomes pertinent for future investigations as binding orientation in the pocket affected the predicted binding free energies by approximately 17.0 kcal/mol when using RI-B3PW91-D3 for guests that only differed by one CH2 group.

All methods presented predict overbinding character for these host–guest systems except for RI-B3PW91 for CB8 host–guest systems. MMPBSA and RI-B3PW91-D3 worked well at ranking binding affinities for smaller guests regardless of the size of the host. The CB8 guest molecules with a larger van der Waals volume yielded poor prediction of binding free energy due to their higher conformational flexibility, which can complicate predicting binding poses. To better understand binding free energies of these large structures, enhanced sampling methods can be used, and multiple host–guest binding poses can be sampled.