Absolute binding free energy calculations of CBClip host–guest systems in the SAMPL5 blind challenge

Outline of the binding free energy calculations of CBClip and guest molecules

The overall workflow of our SAMPL5 blind binding free energy calculations is shown in Fig. 2. First, we generated the putative binding poses of CBClip and guest molecules in two ways: (1) by slowly creating a guest molecule inside CBClip in the gas phase with a harmonic restraint to keep the guest near the center of mass of CBClip (referred to as “-MD” sets) and (2) by performing docking simulations using the GalaxyDock program (referred to as “-dock” sets), which finds the putative binding poses of protein–ligand, protein–protein and host–guest complexes via highly efficient global optimization [44–47] of the AutoDock4 scoring function [45, 48, 49]. We obtained the initial conformations for the MD simulations from scratch and did not use the conformations provided by the organizers. To obtain the initial conformations for the MD simulations, CBClip and guests were centered at the origin and rotated so that the first principle axis of inertia matches with X-axis using the “COOR ORIENT” command in CHARMM. All initial complex structures were solvated and neutralized with explicit water molecules and counter ions followed by brief molecular dynamics simulations to equilibrate. During the MD simulations, harmonic restraints were imposed to the complexes to keep the initial conformation. The last snapshots of the equilibrium runs were used as the input structures for absolute binding free energy calculations based on the double decoupling scheme [50, 51], which requires alchemical free energy calculations of the guest molecules in the bound and free states.

Free energy calculations were carried out in two ways: (1) a combination of Hamiltonian replica change (HREM) [52–54] and Bennett acceptance ratio (BAR) [11, 18, 55, 56] and (2) thermodynamic integration (TI) [57, 58].After the binding free energy calculations, if a guest molecule has one or more titratable groups whose pKa values are in the range of experimental condition, we performed constant-pH calculations of the bound and free states to consider the effect of possible protonation states [59]. Based on the combinations of the two docking and two free energy calculation methods, we submitted four combinations of prediction schemes: ‘TI-dock’, ‘BAR-dock’, ‘TI-MD’ and ‘BAR-MD’. In addition, we submitted the set that consists of the minimum free energy predictions from the four sets: ‘All’.

Obtaining the initial bound conformations via docking simulation

With the refined CBClip force field parameters, we performed 10 ns MD simulations of the host solvated with explicit water. Using the ART-2 clustering algorithm [60] implemented in CHARMM, we clustered the sampled CBClip conformations from the MD simulation and selected the four largest clusters for following docking simulations.

For each of the four host conformations obtained from MD simulations, up to three bound conformations were selected after host–guest docking with each guest molecule. The docking method, called GalaxyDock-HG, performs global optimization of the AutoDock4 energy [61] with the conformational space annealing (CSA) algorithm [44–49, 62], as illustrated in Fig. 3. The AutoDock4 energy function consists of van der Waals energy, directional hydrogen bond energy, Coulomb electrostatic energy, and desolvation free energy as follows:

{A_ij, B_ij} and {C_ij, D_ij} are parameters for the van der Waals energy and hydrogen bond energy, respectively. h(t) is the weight factor for hydrogen bond to describe its directionality. q_i are partial charges, and ε(r) is a distance-dependent dielectric constant. {S_i, V_i, σ} are desolvation energy parameters. AutoDock4 energy parameters were used in the all-hydrogen topology except that the atomic partial charges of CGenFF were used.

GalaxyDock-HG was developed based on the GalaxyDock protein–ligand docking program [48] with the following modifications: (1) the potential energy was evaluated in the continuous space instead of interpolating energy values at the grid points because both host and guest molecules are treated flexibly and the host–guest systems were small enough to run global optimization efficiently in the continuous space, and (2) the initial set of conformations for CSA, called initial bank, was generated by randomly perturbing initial structures instead of using the geometry-based pre-docking method of GalaxyDock because more intensive pre-docking did not seem necessary for the small systems. More specifically, the initial structures were generated by randomly perturbing the flexible torsion angles of the host molecule first and then by randomly placing a guest molecule onto each host conformation by perturbing the values of translational, rotational, and torsional degrees of freedom for the guest molecule. A total of 100 conformations were selected as the initial bank after local energy minimization, and the bank was evolved by the CSA algorithm as implemented in GalaxyDock. The same docking run parameters were used as GalaxyDock except that the number of seed conformations used to generate trial conformations was increased from 25 to 50 because the bank size was increased from 50 to 100.

Performance of GalaxyDock-HG was tested on a set of 68 host–guest complex structures selected from those generated by Wang and Pang [63]. Successful docking results with guest-RMSD<2 Å were obtained in 64.7 %/76.5 % of the test set complexes when top 1/top 3 complex structures were selected by clustering the final bank conformations and selecting the lowest-energy conformations among the lowest-energy representatives of the clusters. These success rates were higher than those when the conformations were selected from the final bank conformations by energy only (64.7 %/70.6 %) or by the cluster size (52.9 %/76.5 %). Therefore, we selected up to three host–guest docking conformations for each of the initial host conformation based on energy after clustering, which resulted in 12 conformations for each guest. Among the 12 conformations, the top two scored conformations were used as initial docked conformations for the absolute binding free energy calculations.

Absolute binding free energy calculation

Binding affinities of all guest-CBClip complexes were calculated using the double-decoupling scheme [50, 51] (Fig. 4). This scheme allows one to calculate the absolute binding free energy of a given conformation of a complex. The key to obtaining well-converged results with this scheme is to restrain the relative geometry of a ligand to prevent it from being separated from the host when guest-host interactions become very weak. The thermodynamic cycle shows that the absolute binding free energy of a complex is the difference between the decoupling free energy of a ligand in the bound state and the solvation free of a ligand in the monomeric state:

where $\mathrm{\Delta }{G}_{\text{rest}-\text{on}}^C$ and $\mathrm{\Delta }{G}_{\text{rest}-\text{off}}^C$ are the free energies of turning on and turning off restraints, $\mathrm{\Delta }{G}_{\text{elec}}^C$ and $\mathrm{\Delta }{G}_{\text{vdW}}^C$ are the free energies of decoupling the electrostatic and the vdW interactions of a guest in the bound state, and $\mathrm{\Delta }{G}_{\text{elec}}^M$ and $\mathrm{\Delta }{G}_{\text{vdW}}^M$ are those of a guest in the monomeric state.

To calculate the free energies of a complex, six geometric restraints, one distance, two angle, and three dihedral restraints, are gradually turned on between a guest and a host. With the presence of restraints, the electrostatics and vdW free energies between the guest and host are calculated by decoupling the interactions. The free energy of removing restraints $\mathrm{\Delta }{G}_{\text{rest}-\text{off}}^C$ was calculated by the following formula:

where V is the volume of the unit box, r₀ is a reference distance of a distance restraint, ϑ_A and ϑ_B are the reference angles of angle restraints, K values are the force constants of the six restraints.

All simulations were performed with CHARMM 39b2 [64]. The force field parameters of all ligands were generated with the CGENFF server [65, 66]. The dihedral parameters and atomic partial charges with high penalty scores were further optimized using the standard CGENFF procedure [67, 68]. All complex structures were solvated in a cubic box whose dimensions were 40 × 40 × 40 ų with TIP3P water molecules.

To calculate the free energy contributions of the electrostatic and van der Waals (vdW) interactions, we performed a Hamiltonian replica exchange (HREM) [52–54] simulation and post processed the results using Bennett’s acceptance ratio (BAR) [11, 18, 55, 56] method. For the HREM simulation, we used 32 replicas. The electrostatic interactions were scaled down linearly using 12 replicas. The vdW interactions were gradually turned off based on the serial insertion scheme [69] (i.e., without using soft-core potentials). Twenty nonlinear scaling factors for vdW interactions were used to make exchange ratios between replicas uniform: 0.995, 0.985, 0.960, 0.930, 0.900, 0.850, 0.800, 0.750, 0.700, 0.650, 0.600, 0.550, 0.500, 0.450, 0.390, 0.330, 0.270, 0.220, 0.180, and 0.0. All HREM simulations were performed for 1 ns, which makes 32 ns in total. The free energy calculations of unbound ligands were performed using the same procedures as the complexed calculations. All simulations were repeated three times and averages and standard deviations were calculated from the three trials.

To compare the performances of different free energy calculation methods, we also performed thermodynamic integration calculations. To calculate the free energy cost of turning on restraints between the host and the ligand, we performed thermodynamic integration calculations using 19 windows with lambda values of 0.0005, 0.002, 0.004, 0.00625, 0.01, 0.01875, 0.03, 0.05, 0.075, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, and 1.0. The electrostatic and van der Waals free energy values were calculated using 15 windows with lambda values of 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.675, 0.725, 0.775, 0.825, 0.875, 0.925, 0.975, 1.0. For each window, equilibrium and production runs were performed for 20 and 100 ps, respectively, which amounts to 6 ns of sampling time in total for each target. The λ values were determined from a series of short TI calculations with G1 so that the variances of ${\langle \frac{\mathrm{\partial }U(\mathit{x},\lambda )}{\mathrm{\partial }\lambda }\rangle }_{\lambda }$ remain lower than 1 kcal/mol, which are specified by the DIFFLC variable in the log of the PERT command in CHARMM.

After the SAMPL5 submission, we performed additional longer TI calculations using the docked conformations (“-dock” set) to check the convergence of the TI calculations. For the TI-long simulations, we extended the simulation time of each step by a factor of five, which amounts to 29.5 ns of simulation time in total. Production simulations were run in the NVT ensemble using the Nose–Hoover thermostat [70–72]. All simulations were performed at temperature 300 K. Direct space nonbonded interactions were truncated with a 14.0 Å cutoff, whereas long-range electrostatics were handled with the PME method [73]. SHAKE constraints were applied to bonds involving hydrogen, and the simulation time step was set to 1 fs [74].

Free energy methods

The Hamiltonian replica exchange method (HREM) uses multiple replicas to sample a phase space simultaneously, while replicas are swapped periodically [52]. Each replica corresponds to a different environmental condition. It is well known that proper exchanges between Hamiltonians enhance sampling efficiency while preserving the Boltzmann distribution [53, 54, 75, 76]. In HREM, the ratio between the probability of exchanging x in the mth replica and y in nth replica W(x, E_m; y, E_n) and that of the reverse process W(y, E_m; x, E_n) must follow the relation: $\frac{W(\mathit{x},E_m;\mathit{y},E_n)}{W(\mathit{y},E_m;\mathit{x},E_m)}=exp(-\beta \mathrm{\Delta })$ , where Δ = E_m(y) − E_m(x) + E_n(x) − E_n(y) and this condition can be satisfied if the Metropolis-type criterion is used for exchanges; if Δ>0, W(x, E_m; y, E_n) = exp(−βΔ), otherwise W(x, E_m; y, E_n) = 1.

TI calculates the free energy difference between two states 0 and 1 by using the derivatives of the potential energy with respect to a state variable λ [6, 58], $\frac{\mathrm{\partial }U(\mathit{x},\lambda )}{\mathrm{\partial }\lambda }$. Using integration, the free energy difference between state 0 and 1 is $\mathrm{\Delta }{G}^{TI}=\underset{0}{\overset{1}{\int }}{\langle \frac{\mathrm{\partial }U(\mathit{x},\lambda )}{\mathrm{\partial }\lambda }\rangle }_{\lambda }$ . In practice, multiple simulations are performed at different λ values to evaluate ${\langle \frac{\mathrm{\partial }U(\mathit{x},\lambda )}{d\lambda }\rangle }_{\lambda }$. The total free energy difference is obtained from the obtained ensemble averages of the derivatives using numerical integration schemes, such as the trapezoidal rule or Simpson’s rule [77, 78]. With the trapezoidal rule, the free energy difference is obtained from the following equation:

where $g(\lambda )={\langle \frac{\mathrm{\partial }U(\mathit{x},\lambda )}{d\lambda }\rangle }_{\lambda }$ and n is the total number of intermediate states between state 0 and 1. By using Simpson’s rule with the generalization of Brun that allows to use non-equidistant λ values [79], the integration can be performed using the following equation:

In contrast to TI, BAR uses information from a pair of simulations. From the ensembles sampled obtained with two simulations, the free energy difference is calculated according to

where f is the Fermi function and C is a constant that can be obtained from an iterative procedure [80, 81]. Here, for the van der Waals free energy calculations using BAR, we used the serial insertion scheme [69]. In the serial insertion scheme, van der Waals interactions are turned on or off on an “atom by atom” basis. The advantage of this scheme is that it does not require using soft-core potentials preventing from using fast energy calculation routines.

Parameterization of CBClip

The initial force field parameters of CBClip were generated using the ParamChem automation server for the CHARMM general force field (CGenFF) [65, 66]. To sample possible different conformations of apo-CBClip, we performed multiple 10 ns MD simulations with explicit water starting from the initial conformation provided by the organizers (Fig. 5B). In all MD simulations, the two aromatic walls of the CBClip opened and formed a rather planar conformation mostly within 2 ns (Fig. 5C). In addition, with the opened conformation, the CBClip–guest complexes did not remain stable and become separated in the equilibrium MD simulations, preventing us from performing absolute binding free energy calculations. Thus, we assumed that the opened CBClip conformation is an artifact of the default CGenFF parameters.

We hypothesized that this opened conformation is formed due to the weak force constants of the dihedral angle potentials of 1-3-diazepine rings in CBClip, which were assigned high penalty scores by ParamChem, indicating the possible need to reparameterize. To validate the quality of these dihedral angle parameters, we compared the potential energy surfaces of the 1-3-diazepine rings with the CGenFF parameters and quantum mechanical (QM) calculation at the RIMP2/cc-pvdz level (Fig. 5D). The comparison shows that the initial CGenFF parameters favor the open conformation by about 2 kcal/mol (Fig. 5E). However, the QM calculation favors the closed conformation and has a higher energy barrier than the CGenFF model, suggesting more accurate modeling of the host molecule is needed for the binding affinity calculations.

We optimized the force field parameters of CBClip according to the principle of CGenFF parameterization [67, 68]. The atomic partial charges were optimized by targeting the water interaction computed at the HF/6–31 g* level with model compounds including glycoluril, hexamethyl-glycoluril, and acylic cucurbit[2]uril. The dihedral parameters for the 1-3-diazepine ring in CBClip were optimized using a substituted glycoluril (SMILES: [C@H]12[C@@H]3NC(=O)N1Cc1c(CN2C(=O)N3)cccc1) as model compound, targeting 2D QM potential energy surfaces computed. With the optimized parameters, the closed conformation is energetically more favorable than the open conformation by 2 kcal/mol, with a barrier height of about 11 kcal/mol (Fig. 5F), suggesting this large-scale conformational change might be of importance in binding free energy calculations.

Considering multiple protonation states using constant-pH simulation

We assumed that the pH of the solvent is the same as the experimental condition pH = 7.4. The protonation states of guests were determined by comparing predicted pKa values and pH. The pKa values were predicted by using the same approach that we used for the hydration free energy calculations in SAMPL4 [18]. If the predicted pKa value of a guest differs from experimental pH by more than 2.0, we assumed that the protonation state of a guest is fixed. Otherwise, we performed constant-pH simulations to consider multiple protonation states.

For guest molecules G1, G2 and G3, we performed the constant-pH simulation of the complex to take into account of the effect pH on the binding affinity using the 2D-HREM-EDS method [59, 75, 76, 82]. The pKa values of G1 and G2 are 8.9 and 9.1, respectively. The pKa values of G3 are 7.95 and 8.82, respectively. The deprotonation free energies of G1, G2, and G3 were calculated using HREM combined with a BAR calculation.

For all constant-pH simulations, replica exchanges between EDS Hamiltonians and different pH conditions were attempted every 100 steps, alternating between exchanging EDS s values and pH values. The same nonbonded parameters were used as in case of HREM calculations, and simulations were 1 ns long. The simulations were performed at pH 5.4, 7.4 and 9.4 and the EDS temperatures of 0, 7500, 10,000 and 12,000 were used. After the constant-pH simulations of the three complexes were finished, the pKa shifts of the three guests are calculated by fitting the deprotonated fractions to the Hill equation. By using the obtained pKa values, a binding free energy estimated with fixed charges at a given pH is adjusted by using the following equation [59]:

where $\mathrm{\Delta }{G}_{\mathrm{b}}^0$ and $\mathrm{\Delta }{G}_{\mathrm{b}}^{+}$ are the binding free energies obtained with deprotonated and protonated ligands.

Overall performances

Among the five submissions, the TI-dock submission, which was obtained from the docking and TI calculations, shows the best agreement with the experimental data in terms of the average unsigned error (AUE), 2.94 kcal/mol, and the root mean squared deviation (RMSD) measures, 3.41 kcal/mol (Table 1). The BAR calculations (‘BAR-dock’) starting from the same docked conformations with the ‘TI-dock’ submissions resulted in higher AUE and RMSD values (3.89 kcal/mol and 4.62 kcal/mol). The submissions with initial complex conformations from MD simulations, ‘TI-MD’ and ‘BAR-MD’, resulted in worse RMSD than TI-dock. Opposite to the cases of the ‘dock’ submissions, ‘BAR-MD’ showed better results than ‘TI-MD’, which indicates that neither ‘TI’ or ‘BAR’ did not outperform the other consistently for the CBClip challenge. We also submitted the “All” set, which consists of the minimum value of all predictions. This set is based on the assumption that if all free energy calculations are accurate and one docked conformation has lower binding free energy than the other; the contribution of the lowest binding free energy conformation will contribute most to the final free energy. However, the RMSD of ‘All’ is higher than that of ‘TI-dock’, which indicates that the ‘All’ set may suffer from the overestimation of binding free energy due to the neglect of the relative free energies of the involved conformations.

To compare different protocols statistically, we estimated the error bounds of the AUE and RMSD values by performing a bootstrapping calculation [83]. We assumed that individual free energy terms, $\mathrm{\Delta }{G}_{\text{rest}-\text{on}}^C,\mathrm{\Delta }{G}_{\text{elec}}^C$, and $\mathrm{\Delta }{G}_{\text{vdW}}^C$, obey the Gaussian distribution whose mean and standard deviation are obtained from the three repetitions. Based on the Gaussian distributions, random values were sampled 1000 times. The standard deviations of the AUE and RMSD values were calculated by comparing the random and the experimental values (Table 1). The boostrapping analysis indicates that the results of all protocols are different from each other in a statistically significant way.

Comparison with other submissions

Our best submission, ‘TI-dock’, ranked as the top method among all submissions with the RMSD and AUE metrics. With the same metrics, our other submissions, ‘All’, ‘BAR-MD’, ‘BAR-dock’, and ‘TI-MD’ sets ranked the 3rd, 4th, 6th, and 7th best submissions, which shows that overall performance of our methods are better than other methods. The second best result was submitted by the movable type method devised by the Mertz group [84]. While our submissions were consistently ranked near the top by RMSD and AUE metrics, our results did not correlate very strongly with experiment. Other submissions, especially those by the Michel group, were able to rank order their predictions very well.

Dependency on initial docked conformations

The large difference between the ‘dock’ and the ‘MD’ results indicates that absolute free energy calculations using the double decoupling scheme heavily depend on initial docked conformations because the scheme restraints the relative geometry of ligand. If the native conformation is significantly different from the initial conformation, it can be only accessible during the early stages of calculating $\mathrm{\Delta }{G}_{\text{rest}-\text{on}}^C$, when the restraints are weak. To avoid this problem, one may use very weak restraints. However, if the restraints are too weak, conformational sampling becomes more extensive. Moreover, a weak restraint may slow down the convergence of calculations and the calculation may become unstable due to severe steric clashes when vdW interactions are almost turned off. Thus, with weak restraints, it would be necessary to perform very long free energy simulations. On the other hand, if restraints are strong, a small conformational basin near an initial conformation will be sampled, which makes the free energy calculation converge faster. However, the contributions of conformations, which differ from the initial structure will be ignored, which biases the free energy calculation [51]. Thus, starting from an accurate docked conformation is one of the most essential steps to obtain an accurate binding free energy using the double decoupling scheme.

Comparison of TI and BAR

To obtain more details on our results, we plotted the computed and the experimental binding free energies of each guest molecule (Fig. 6). Based the docked initial conformations (Fig. 6A), the predictions for G1, G2, G4 and G7 show consistently good agreement with the experimental values regardless of the free energy calculation method. For G8, G9 and G10, the TI calculations result in good or reasonable agreement with experiment, but the BAR calculations significantly underestimate the values. For G3, both predictions made with TI and BAR overestimate the binding free energies by 6 ~8 kcal/mol. Among 10 guest molecules, BAR yielded a better prediction than TI only for G5.

When the initial docked poses are obtained with MD, the patterns of estimated binding free energies change substantially (Fig. 6B). Overall, the BAR calculations are in better agreement with experiment than the TI calculations. Both TI and BAR estimated the binding free energies of G1, G4, G5, and G9 reasonably well within an error of 3 kcal/mol. However, TI significantly underestimated the binding free energies of G6 and G10 compared to BAR, which is opposite to the ‘TI-dock’ and ‘BAR-dock’ results. The results of four ligands, G2, G3, G7, and G8, deviate from the experimental values for both TI and BAR. Especially the results of G2 and G7 deteriorate compared to the ‘TI-dock’ and ‘BAR-dock’ results, which indicates that the predicted docked poses of these molecules may have played an important role for the accuracy of the binding free energy calculation. Similar to the ‘TI-dock’ and ‘BAR-dock’ results, the binding of G3 to CBClip was highly overestimated with TI and BAR.

For a more detailed comparison of TI and BAR, the contributions of the free energy terms used in the double decoupling scheme of the ‘TI-dock’ results are listed in Table 2. The major discrepancy between the two methods is observed for the free energy contributions of the vdW interactions in the complex and the monomeric states, $\mathrm{\Delta }{G}_{\text{vdW}}^C$ and $\mathrm{\Delta }{G}_{\text{vdW}}^M$. For all complexes and monomers, BAR exhibits lower vdW free energy values than TI. The average $\mathrm{\Delta }{G}_{\text{vdW}}^C$ and $\mathrm{\Delta }{G}_{\text{vdW}}^M$ values estimated with BAR was smaller than those from TI by 2.7 kcal/mol and 3.5 kcal/-mol, respectively. The difference between the vdW free energy values, ~0.83 kcal/mol, accounts for about 75 % of the average total free energy differences $\langle \mathrm{\Delta }{G}_{\text{total},\mathrm{i}}^{\text{TI}}-\mathrm{\Delta }{G}_{\text{total},\mathrm{i}}^{\text{BAR}}\rangle$, which is −1.1 kcal/mol. In contrast, the average difference between the electrostatic free energy values is −0.27 kcal/mol.

One possible cause of the difference between TI and BAR is that the numerical quadrature in the TI calculations based on the trapezoidal rule is not converged as we have observed in SAMPL3 [11, 77]. To identify whether this is the case, we re-calculated the vdW free energies from the TI calculations using Simpson’s rule (Table 3). When the TI results are integrated with Simpson’s rule, the vdW free energies changed to some degree. The average absolute differences between $\mathrm{\Delta }{G}_{\text{vdW}}^C$ values integrated with the trapezoidal rule and Simpson’s rule results for long and short TI calculations are 1.4 and 1.2 kcal/mol, respectively. These results indicate that the error from numerical quadrature is observed in the longer calculations. However, the numerical quadrature error appears not to be the major source of the difference between the TI with soft-core potential calculations and the HREM–BAR with serial-insertion calculations.

Another possibility is that the TI and BAR calculations sampled different conformational spaces and were not fully converged. The major differences in the protocols for the TI and BAR calculations are the use of soft-core potentials (TI) versus serial insertion (BAR) for the van der Waals free energy calculation, and the use of enhanced sampling with HREM for the BAR calculations, which leads to sampling of different conformational spaces. We observed that the discrepancy between the TI and BAR calculations remains even for the longer TI calculations. This suggests that the TI and BAR calculations are not fully converged and the use of soft-core potential may lead to significantly different conformational space sampling, which hinders the convergence of calculations. More systematic investigations are necessary to identify the source of the difference between the convergence of conventional TI calculations with soft-core potentials and free energy calculations using the serial insertion scheme without soft-core potentials [69].

Convergence of TI calculations

The free energy components of the longer (30 ns) TI calculations are listed in Table 4 to compare the convergence of the longer and shorter TI calculations. The largest differences are observed in G3, G7 and G8 and the longer TI calculations resulted in closer agreements with experiments. The absolute errors of (G3, G7, G8) are reduced from (5.6, 1.8, 3.5) kcal/mol to (0.6, 0.3, 1.7) kcal/mol leading to the average AUE of 2.2 kcal/mol for the longer TI calculations. This result suggests that the free energy calculations of these guests were not properly converged with the short TI calculations. The slower convergences of these guests may be attributed to the fact that they have more flexible bonds than the other guests. The differences between the other guests are less or equal than 1.1 kcal/mol.

Effect of the accuracy of force field parameters

One possible source of error in calculating host–guest binding free energies is the accuracy of force field parameters. The automatic CGENFF procedure provides an estimate of the quality of the parameters with empirical penalty scores based on the chemical similarity between a given molecule and reference molecules with pre-calculated parameters. If the penalty score of a guest is correlated with the accuracy of binding free energy, one may obtain the approximate reliabilities of free energy calculations when experimental data are unavailable. To identify how the accuracy of force field parameters affects the accuracy of free energy calculations, we investigated the correlation between the penalty scores of the CGENFF parameters of guests and the average AUE values of four submissions (Fig. 7). When all results are considered, there is little correlation between two values with an R² value of 0.09. However, if the results of G3 and G9 are excluded, the correlation becomes much stronger; the R² value increases to 0.86. These results suggest that a free energy calculation of a molecule with a lower penalty score is likely to be more accurate than a calculation with a higher penalty score. However, since we investigated the correlation only with 10 guests of SAMPL5, more extensive investigations are necessary. For example, the correlation coefficient between the CGenFF penalty score and AUE for hydration free energies in SAMPL4 was 0.4.

Effect of multiple protonation states

From constant-pH simulations, we obtained the population distributions of different protonation states of G1, G2, and G3. For G1 and G2, the guests in the bound state remained fully protonated for over 99 % of the simulation time at the experimental pH of 7.4. Thus, we did not apply any corrections due to multiple protonation states. From the constant-pH simulation of the CBClip-G3 complex, it was identified that G3 was in the fully protonated state for 98.6 % and in the singly protonated state for 1.4 % of trajectory, which corresponds to a pKa shift of 9.6. Based on the constant-pH result, the total binding free energy was shifted by −0.16 kcal/mol using Eq. 2. Thus, the influence of the protonated states was rather small for the systems under study.

What went right and what went wrong?

From the four different sets of simulations, the relatively small and rigid molecules, G1, G4, and G9 were predicted accurately. The conformations of G1 obtained with docking and MD simulations are similar to each other (Fig. 8A). The benzene ring of G1 is stocked with the two naphthalene sidewalls of CBClip in a parallel fashion. The two ligand conformations of G4 are the mirror images of each other (Fig. 8B). Thus, they adopt similar interaction patterns with CBClip and resulted in similar binding free energies. However, the two docked conformations of G9 are significantly different from each other although their calculated binding free energies are similar (Fig. 8C). The docking simulation placed G9 between the two sidewalls of CBClip in a parallel fashion. However, G9 is dangling around a sulfonate group via electrostatic interaction. The similar binding free energies of the two conformations suggest that the G9-CBClip complex may adopt diverse docking poses. The sulfonate arms of CBClip may play a role not only in enhancing solubility but also in attracting ligands [40].

The binding free energies of G3, G8, and G10, which are bulky and flexible molecules with many rotatable bonds, were poorly calculated. From both docking and MD simulations, G3 binds to one of the sulfonate arms of CBClip via electrostatic interactions (Fig. 8D). The overestimation of binding free energy of G3 indicates that the electrostatic interactions are significantly exaggerated. On the other hand, G8 is weakly bound to CBClip via interactions with the sulfonate groups and does not form a stable complex (Fig. 8E). The significantly underestimated binding free energy of G8 suggests that the docked conformations are incorrect and more extensive conformational sampling is necessary. The aromatic groups of G10 are stacked in between two sidewalls in the two docked conformations, which is similar to the G9-CBClip complex conformation obtained with the docking simulation (Fig. 8F). When the individual binding free energy components of G9 and G10 are compared, the electrostatic free energies in the complex state ( $\mathrm{\Delta }{G}_{\text{elec}}^C$) differ the most, which suggests that electrostatic interactions between G10 and CBClip may be underestimated.