Binding thermodynamics of host-guest systems with SMIRNOFF99Frosst from the Open Force Field Group

This manuscript (permalink) was automatically generated from slochower/smirnoff-host-guest-manuscript@ff89e76 on April 24, 2019.

Authors

• David R. Slochower
  0000-0003-3928-5050 · slochower · drslochower
  Skaggs School of Pharmacy and Pharmaceutical Sciences, University of California, San Diego, La Jolla, CA 92093, USA

• Niel M. Henriksen
  0000-0002-7916-0757 · nhenriksen
  Atomwise, Inc., San Francisco, CA 94105, USA

• John D. Chodera
  0000-0003-0542-119X · jchodera · jchodera
  Computational and Systems Biology Program, Sloan Kettering Institute, Memorial Sloan Kettering Cancer Center, New York, NY 10065

• David L. Mobley
  0000-0002-1083-5533 · davidlmobley · davidlmobley
  Department of Pharmaceutical Sciences and Department of Chemistry, University of California, Irvine, CA 92697, USA

• Michael K. Gilson
  0000-0002-3375-1738 · michaelkgilson
  Skaggs School of Pharmacy and Pharmaceutical Sciences, University of California, San Diego, La Jolla, CA 92093, USA

Note, I only added authors who reviewed the outline on Google Docs. We can revisit this.

Abstract

Designing ligands that bind their target with high affinity and specificity is a key step in small-molecule drug discovery. Yet accurate predictions of protein-ligand binding free energies are difficult and errors in the calculations can be traced to challenges adequately sampling conformational space, ambiguous protonation states, and other causes. Noncovalent complexes between a cavity-containing host molecule and drug-like guest molecules have emerged as powerful tools for modeling protein-ligand binding. Due to their small size and extensive experimental characterization, calculations of host-guest binding free energies, enthalpies, and entropies offer an opportunity to directly probe, and ultimately optimize, force fields.

The Open Force Field Initiative aims to create a modern, open software infrastructure for automatically generating and validating force fields using high-quality data sets. The first force field to arise out of this effort, named SMIRNOFF99Frosst, has one tenth the number of parameters of a typical general small molecule force field, such as GAFF, yet predicts binding thermodynamics that are on average, at least as accurate. Here, we report the results of free energy calculations on 43 α and β-cyclodextrin host-guest pairs for which experimental data are available. Our calculations were performed using the attach-pull-release method as implemented in the open source package, pAPRika. On binding free energies, the root mean square error of the predictions relative to experiment is 0.82 kcal/mol for SMIRNOFF99Frosst and 1.58 kcal/mol for GAFF version 2.1. These results suggest significant room for improvement in force fields, and will help create a transparent and robust method of evaluating future candidate parameter sets from the Open Force Field Group. Improving the performance of force fields for predicting binding affinities will help reduce the timescale and cost required to generate drug candidates.

Introduction

Accurate predictions of protein-ligand binding free energies are a key goal of computational chemistry. Despite this, calculations of protein-ligand binding thermodynamics involve a number of challenging choices, including specifying the protonation state of ionizable residues, adding hydrogens or otherwise adjusting the initial protein structure, and placing the ligand in the binding pocket, for which there is no consensus in the computational chemistry community. Recently, calculations on the reversible work of ligand binding to bromodomain proteins was completed [1], with deviations from experiment in the 1-2 kcal/mol range and uncertainties on the order of ~0.5 kcal/mol. Here, we report the calculation of binding free energies, enthalpies, and entropies of drug-like guest molecules to α- and β-cyclodextrin host molecules, converged to within ~0.1 kcal/mol, using the attach-pull-release method. These calculations, which are easier to sample and have been experimentally characterized using a variety of methods, offer an opportunity to benchmark—and ultimately optimize—new and existing force fields. We compare the predictions of three force fields: GAFF v1.7 [2], GAFF v2.1, and SMIRNOFF99Frosst [3,4].

SMIRNOFF99Frosst, released in late 2018, is the first force field produced by the Open Force Field Initiative. SMIRNOFF99Frosst is derived from AMBER parm99 [5] and Merck’s parm@Frosst [6]. Instead of relying on atom types to assign force field parameters to compounds, which is the procedure followed by the tleap program that parameterizes molecules in AmberTools, SMIRNOFF99Frosst and the Open Force Field Toolkit use the local chemical environment of each atom to apply force field parameters using SMIRKS strings [7]. This process simplifies and effectively uncouples the parameters for each term in the force field. That is, the addition of a new Lennard-Jones parameter does not require the addition of new bonded, angle, and dihedral parameters involving the same atom. These factors lead to a much more lean force field specification; there are over 3000 lines of parameters in GAFF v1.7, over 6000 lines of parameters in GAFF v2.1, and just ~300 parameters in SMIRNOFF99Frosst version 1.0.5.

Thus far, SMIRNOFF99Frosst has been tested on hydration free energies of 642 small molecules, and the densities and dielectric constants of 45 pure organic liquids [3]. Here we benchmark SMIRNOFF99Frosst using noncovalent binding thermodynamics using two flexible host molecules and thirty three guests containing three different functional group moieties. We first show that SMIRNOFF99Frosst does about as well as the conventional force fields, GAFF v1.7 and GAFF v2.1, predicting experimental binding free energies, enthalpies, and entropies. We then characterize the conformational differences produced by SMIRNOFF99Frosst compared to the other force fields.

Methods

Choice of host-guest systems

In this study, we report the binding thermodynamics of 43 host-guest complexes (Figure 1 and Table 1) computed using three different force fields. The complexes consist of either α- or β-cyclodextrin as host molecules and a series of small molecule guests containing ammonium, carboxylate, or cyclic alcohol functional groups. Cyclodextrins are cyclic polymers consisting of six (αCD) or seven (βCD) glucose monomers in the shape of a truncated cone. The equilibrium constants and standard molar enthalpies of binding for these 43 complexes have been measured using isothermal titration calorimetry (ITC) and nuclear magnetic resonance spectroscopy (NMR) [8] and computationally in [9]. As in Henriksen, et al. [9], only a single stereoisomer was considered for the 1-methylammonium guests because it was not clear whether a mixture or a pure solution was used in Rekharsky, et al. [8], and the ΔG difference between each stereoisomer is expected to be < 0.1 kcal/mol [10].

Application of force field parameters

We sought to compare force fields directly, and as such, attempted to minimize additional differences between the simulations with each force field. In all simulations, we applied AM1-BCC [11,12] partial atomic charges to both the host and guest molecules using the antechamber program. The host charges were calculated using a single glucose molecule with methoxy caps on the O1 and O4 alcohols (Figure 2); each glucose monomer in the cyclodextrin polymer has identical charges. We used TIP3P water [13] and Joung-Cheatham monovalent ion parameters [14] in each simulation set.

GAFF v1.7 bond, angle, torsion, and Lennard-Jones parameters were applied using the tleap program distributed with AmberTools16. These simulations were performed as part of Henriksen, et al. [9] and are described in additional detail therein.

GAFF v2.1 parameters were applied in an identical manner to the GAFF v1.7 parameters, using the tleap program distributed with AmberTools18 and substituting leaprc.gaff for leaprc.gaff2 in the tleap input file. In GAFF v2.1, the bond and angle parameters have been updated to reproduce small molecule geometries obtained from high-level quantum mechanical calculations. The force constants for the bond and angle parameters were tuned to reproduce the vibrational spectra of over 600 molecules. The torsion parameters were optimized to reproduce the rotational potential energy surface of 400 model compounds. Finally, the Lennard-Jones coefficients were redeveloped to reproduce interaction energies and pure liquid properties. See #7.

To apply SMIRNOFF99Frosst parameters, we followed a multistep process, beginning with the prepared GAFF v1.7 files. The host and guest molecules were parameterized with the Open Force Field Toolkit version 0.0.3, SMIRNOFF version 1.0, and SMIRNOFF99Frosst version 1.0.5. Once parameterized with SMIRNOFF99Frosst, ParmEd [15] was used to combine the host and guest with the solvent and ions, which retained their TIP3P water parameters and Joung-Cheatham ion parameters, respectively. See #8.

Thermodynamic calculation

We used the attach-pull-release (APR) method as implemented in the open source package pAPRika version 0.0.3, to calculate absolute binding free energies. A complete description of the APR method has been characterized in the literature [16,17,18,19]. The attachment and release phases consisted of 15 independent windows. During the attachment phase, the force constants on the host and guest are scaled by a \(\lambda\) parameter that goes from \(\lambda = 0\), at which point all restraints are turned off, to \(\lambda = 1\), at which point all restraints are at their maximum force constant. The \(\lambda\) windows are more densely spaced where the force constant is smaller to improve sampling along highly curved regions of the potential of mean force. Conformational restraints were applied between neighboring glucose units of the cyclodextrin to limit the incursion of monomers into the host cavity. These restraints were applied along the pseudodihedrals \(\ce{{O5}_n-{C1}_n-{O1}_n-{C4}_{n+1}}\) and \(\ce{{C1}_n-{O1}_n-{C4}_{n+1}-{C5}_{n+1}}\) to improve convergence and sampling of the bound state (Figure 2 for atom names). To further improve sampling, we applied a hard wall restraint that confined the guest molecule to within a sphere of 12.3 and 13.5 Å of αCD and βCD, respectively, during the bound state. The release phase is the conceptual reverse of the attach phase, in which the conformational restraints on the host are gradually turned off (\(\lambda =1 \rightarrow 0\)) in the absence of the guest. This explicit release phase is performed once for αCD and once for βCD. Finally, an analytic correction is performed to compute the work of moving the guest from the restricted volume enforced by the APR restraints to standard state at 1 M concentration.

The pulling phase consisted of 45 independent, equally spaced windows. During the pulling phase, the \(\lambda\) parameter represents the target value of a distance restraint with constant force constant. This target distance is increased uniformly in 45 increments of 0.4 Å, yielding windows that separate the host and guest by 18 Å over the course of the calculation.

Due to the asymmetry of the primary and secondary alcohols of cyclodextrin (Figure 17), and the lack of symmetry of the small molecule guests, there are generally two distinct binding poses that do not interconvert during the simulation timescale. To account for this effect, we separately compute the binding free energy and enthalpy for each orientation [16] and combine the results to produce a single value for each host-guest combination using the following equation:

\[\begin{equation} \Delta G = -RT \ln(\exp(-\beta \Delta G_\text{primary}) + \exp(-\beta \Delta G_\text{secondary})). \end{equation}\]

In this manuscript, we refer to calculations whereby the guest is pulled out of the primary face of cyclodextrin with a -p suffix and calculations whereby the guest is pulled out of the secondary orientation with a -s suffix.

Thermodynamic integration was used to compute the binding free energy (ΔG). The binding enthalpy (ΔH) was computed as the difference in mean potential energy of the bound state (in the absence of any restraints) and the unbound state (where the guest is held far away from the host, but the conformational restraints on the host are disabled). The binding entropy (ΔS) was computed by subtraction using ΔG and ΔH.

Simulations

Simulations were performed with the pmemd.cuda module of AMBER 16 (calculations of the GAFF v1.7 force field) and AMBER 18 (calculations of the GAFF v2.1 and SMIRNOFF99Frosst force fields) molecular dynamics software [20]. Each window for each system was independently solvated and simulated.

Solvation consisted of 2000 TIP3P waters for the αCD systems and 2210 waters for the βCD systems in an orthorhombic box. The host and guest were oriented via non-interacting dummy atoms along the simulation box \(z\) axis to minimize the amount of solvent required. Each simulation contained enough Na⁺ or Cl^- ions to neutralize the host-guest complex and an additional 50 mM NaCl to match the experimental conditions in [8]. In the GAFF simulations, hydrogen mass repartitioning [21] was used to adjust the mass of hydrogen atoms to 3.024 Da by transferring mass from bonded heavy atoms, enabling a simulation timestep of 4 fs. Equilibration consisted of 50,000 steps of energy minimization, 100 ps of heating from 0 to 300 K, followed by 2000 ps of additional NPT simulations. A Langevin thermostat, the Monte Carlo barostat, a nonbonded cutoff of 9 Å and default PME parameters, were used for the NPT simulations. Production NPT simulations were run for a minimum of 5 ns and maximum of 50 ns per window, except for the windows used to calculate the enthalpy, which were simulated for 1 μs. In the GAFF v1.7 and GAFF v2.1 simulations, the exact length of each window’s simulation was determined by the restraint coordinate uncertainty. For restraint energy \(U\) in \(\lambda\) window \(i\), each window (except for the windows used to calculate ΔH) was simulated until the value \(w(\lambda_i)\):

\[\begin{equation} w(\lambda_i) = \begin{cases} \left\langle \frac{\partial U}{\partial \lambda_i} \right\rangle_\text{SEM} \cdot \frac{\lambda_{i+1}}{2} & i = 0 \\ \left\langle \frac{\partial U}{\partial \lambda_i} \right\rangle_\text{SEM} \cdot \frac{\lambda_{i+1} - \lambda_{i-1}}{2} & i \in [1, N-1] \\ \left\langle \frac{\partial U}{\partial \lambda_i} \right\rangle_\text{SEM} \cdot \frac{1 - \lambda_{i-1}}{2} & i = N \\ \end{cases} \end{equation}\]

fell below a threshold of 0.02 kcal/mol during the attach phase and 0.1 kcal/mol during the pull phase. In contrast, SMIRNOFF99Frosst production NPT simulations were run for 10 ns per window, except for the first and last window which ran for 1 μs to calculate ΔH.

Statistics

The uncertainty in each simulation window was estimated using blocking analysis [22] and propagated into an uncertainty for ΔG after integrating 100,000 curves and an uncertainty for ΔH by adding in quadrature.

For each force field, we computed the root mean squared error (RMSE), mean signed error (MSE), the coefficient of determination (R²), Kendall’s rank correlation coefficient (τ), and the slope and intercept of the computed properties relative to the experimental values. We also test how the force fields compare to each other using these statistics. Comparisons with experiment have 43 data points, for the 43 unique host-guest complexes listed in Table 1; comparisons between force fields have 86 data points, representing the two separate calculations performed for each host-guest pair.

The overall RMSE and R² values for each comparison are reported as the mean with the 95% confidence interval in brackets. These values are estimated by using the uncertanties assigned to each data point to create 100,000 hypothetical graphs through resampling with replacement. The R² values for each functional group subset is also reported in the bottom right corner in each graph.

Results

This results section is organized as follows. We first present a comparison of SMIRNOFF99Frosst and two iterations of the General AMBER Force Field (GAFF [2]) on predicting binding free energies (ΔG) and binding enthalpies (ΔH) of small molecule guests to α-cyclodextrin (αCD) and β-cyclodextrin (βCD). We then detail how the conformational preferences of guest molecules changes between force fields and finally we summarize the parameter differences between SMIRNOFF99Frosst and GAFF along with the effects of the parameter differences.

Comparison with experimental binding free energies, enthalpies, and entropies

Despite having far fewer numerical parameters, SMIRNOFF99Frosst does about as well as, or even better than, GAFF v1.7 on predicting ΔG and ΔH, compared to experimental values measured using ITC or NMR. SMIRNOFF99Frosst has an overall deviation from experiment under 1 kcal/mol on binding free energies and under 2 kcal/mol on binding enthalpies and entropies across the 43 host-guest systems (Figure 3, Figure 16, Table 2, Table 7, Table 8, and Table 9).

GAFF v1.7 agrees well with SMIRNOFF99Frosst (Figure 20); the overall root mean squared error (RMSE) between the methods is 0.88 kcal/mol. Compared to experiment, GAFF v1.7 has RMSE values of 0.88 kcal/mol, 2.54 kcal/mol, and 2.21 kcal/mol on ΔG, ΔH, and −TΔS. Both SMIRNOFF99Frosst and GAFF v1.7 systematically underestimate the ΔG and ΔH for cyclic alcohols. In some cases, GAFF v1.7 underestimates ΔH by over 3 kcal/mol and up to 5 kcal/mol (b-chp).

The predictions made by GAFF v2.1 exhibit significant differences from those of SMIRNOFF99Frosst (RMSE=1.90 kcal/mol Figure 20) or GAFF v1.7 (RMSE=2.1 kcal/mol, Figure 20). GAFF v2.1 has a large negative systematic deviation from the experimental values in both ΔG and ΔH, yet strikingly strong correlations with the experimental values across all three functional group classes. Furthermore, GAFF v2.1 has the lowest −TΔS RMSE of any force field (Figure 16). Compared to SMIRNOFF99Frosst and GAFF v1.7, GAFF v2.1 uniformly overestimates both the binding free energy and binding enthalpy, with slopes > 1 in both cases (Table 2).

I’m thinking it might be cleaner and more palatable to separate the comparisions out by ΔG first, then ΔH, and then −TΔS to break things up a bit. I’m made a mush of the results, trying to touch upon the key differences without covering everything, for example, the second paragraph in this section begins by implicitly using ΔG for “agrees well” but then discusses ΔH in the final sentence.

Guest preferences for binding in the primary or secondary cyclodextrin cavity

The asymmetry of the hosts and the guests leads to two distinct bound states for each host-guest pair: one where the primary functional group of the guest interacts with the primary alcohols of the host and a second conformation where the primary functional group of the guest interacts with the secondary alcohols (Figure 17).

The difference in binding free energy between either orientation (ΔΔG_orientation) can be large, around ~2 kcal/mol for SMIRNOFF99Frosst and GAFF v1.7 and ~5 kcal/mol for GAFF v2.1. SMIRNOFF99Frosst predicts the largest ΔΔG_orientation for the ammonium-containing butylamine and pentylamine with αCD (Figure 4), with the primary orientation being more favorable. GAFF v1.7 predicts a large ΔΔG_orientation for the cyclic alcohols cyclooctanol and cycloheptanol, with the secondary orientation having a more favorable ΔG. This effect is even more apparent with GAFF v2.1 where the ΔΔG_orientation for a-chp and a-coc is greater than 4 kcal/mol. This effect is due, at least in part, to sampling challenges in the bound state for very large guests(Figure 4, bottom right), especially in the narrow primary cavity of the smaller α-cyclodextrin.

Guest preferences for αCD and βCD

Ten guests in the data set bind both αCD and βCD. These ten guests show different patterns of binding between the two host molecules. For example, SMIRNOFF99Frosts underestimates the binding free energy of cyclooctanol both orientations by the same amount (Figure 5). However, despite underestimating the binding free energy of cyclobutanol for αCD by ~0.7 kcal/mol, the binding affinity prediction for βCD is very slightly overestimated by ~0.3 kcal/mol (Table 7). Very similar patterns are observed for GAFF v1.7 and both force fields appear to perform better overall on binding affinities to αCD compared to βCD (Figure 23).

As is the case with the difference in binding free energy between guest orientations, the difference in binding free energy between host molecules using GAFF v2.1 is large. There does not appear to be a clear difference in the accuracy of the predictions for αCD versus βCD (Figure 23).

Trends by guest functional group

SMIRNOFF99Frosst does a good job (MSE = -0.10 kcal/mol and RMSE = 0.76 kcal/mol) estimating the binding free energy of ammonium-containing guests to both αCD and βCD (Figure 6). Shorter chain molecules bind less strongly and the same guest will bind more strongly to αCD than βCD.

I need to include the functional group tables here. I will put the SMIRNOFF99Frosst vs. Experimental tables here and the others in the SI.

SMIRNOFF99Frosst performs reasonably on cyclic alcohols (MSE = 0.70 kcal/mol and RMSE = 1.07 kcal/mol) (Figure 7). The predictions for αCD are uniformly underestimated while those for βCD are mostly under-predicted. The predicted ΔG for cyclooctanol with αCD is particularly poor due to a poor fit in the bound state (Figure 4).

Should say something here about cylooctanol predictions and experiment being poor because it does not fit well, especially in αCD

The binding affinity of carboxylate guests to both αCD and βCD is well characterized by SMIRNOFF99Frosst (MSE = -0.36 kcal/mol and RMSE = 0.87 kcal/mol) (Figure 8).

In all cases, GAFF v1.7 tends to predict slightly weaker binding than SMIRNOFF99Frosst whereas GAFF v2.1 predicts much stronger binding for these compounds (Figure 24, Figure 25, and Figure 26).

Differences in force field parameters between SMIRNOFF99Frosst and GAFF

Next, we summarize the parameter differences among SMIRNOFF99Frosst, a descendant of parm99; GAFF v1.7 (released circa March 2015 according to gaff.dat distributed with AMBER16); and GAFF v2.1 (which is under active development) on the parameters applied to αCD.

The σ and ε parameters are identical between SMIRNOFF99Frosst and GAFF v1.7. Compared to GAFF v2.1, SMIRNOFF99Frosst has deeper well depths for oxygens and decreased σ values for the hydroxyl hydrogens (Figure 9).

Lennard-Jones

Bonded parameters

Compared to GAFF v1.7, SMIRNOFF99Frosst tends to have slightly larger bond force constants, except for the O–H hydroxyl bond force constant, which is much stronger. In GAFF v2.1, the O–H hydroxyl bond force constant is consistent with SMIRNOFF99Frosst, but the carbon-oxygen bond constants are weaker. Equilibrium bond lengths are very similar (Figure 27).

Angle parameters

Relative to GAFF v1.7 and GAFF v2.1, SMIRNOFF99Frosst has fewer unique angle parameters applied to αCD; several distinct parameters appear to be compressed into a single force constant, around 50 kcal/mol/rad² (Figure 11). These parameters correspond to \(\ce{C-C-C}\), \(\ce{C-O-C}\), and \(\ce{O-C-O}\) angles. The \(\ce{C-C-C}\) angles are primarily around the ring of the glucose monomer. The \(\ce{C-O-C}\) angles are both around the ring and between monomers (e.g., \(\ce{C1-O1-C4}\) and \(\ce{C1-O5-C5}\)). Weaker force constants for these parameters may lead to increased flexibility.

Dihedral parameters

The dihedral parameters between SMIRNOFF99Frosst and GAFF v1.7 are extremely similar (where differences occur, they are in the second or third decimal place), with the exception of the \(\ce{H1-C1-C2-O2}\) parameter, for which SMIRNOFF99Frosst applies a dihedral with periodicity = 1 and GAFF v1.7 applies a dihedral with a periodicity of 3 (Table 3 and Figure 12).

The dihedral parameters in GAFF v2.1 differ from those in SMIRNOFF99Frosst, in a number of ways. There are several dihedrals that have a different number of terms in either force field (Table 4).

While this table is accurate, it might not be the best representation because some missing terms happen to correspond to a 0.00 kcal/mol force constant in GAFF v2.1. It might be more useful to exclude those.

In other cases, SMIRNOFF99Frosst and GAFF v2.1 have disagreements on the barrier height after matching the periodicity and phase for a given dihedrals. It is notable that GAFF v2.1 does not have drastically higher force constants for any of the dihedrals, yet GAFF v2.1 produces much more rigid structures (Table 5, Figure 14). The dihedral differences between neighboring glucose monomers demonstrate that SMIRNOFF99Frosst, not GAFF v2.1 has higher force constants (Table 6).

Structural consequences of the force field parameter differences

In both SMIRNOFF99Frosst and GAFF v1.7, the average RMSD of βCD is between 2 and 2.5 Å over 43 μs of simulation. GAFF v2.1 is significanly more rigid, with an average RMSD less than 1.0 Å from the initial structure (Figure 14).

The “flip” pseudodihedral O2_n–C1_n–C4_n+1–O3_n+1 characterizes the orientation of glucose monomers relative to their neighbors. This dihedral is tightly distributed in GAFF v2.1, with all seven dihedrals having a Gaussian-like distribution, centered around -10 degreees (Figure 15,a). In contrast, simulations with both SMIRNOFF99Frosst and GAFF v1.7 report a multipeaked distribution for the dihedral, with a small amount of spread among the individual angles. At any given time point, SMIRNOFF99Frosst adopts a variety of individual pseudodihedral conformations, leading to many conformations with at least one glucose monomer inside the cyclodextrin cavity and distortion of the overall shape of the host binding pocket (Figure 15,b-c). Each pseudodihedral in GAFF v2.1 has a tight distribution; neighboring pseudodihedrals are negatively correlated with each other and positively correlated with the dihedrals on the opposite side of the ring (Figure 15,d).

Discussion

It is striking that SMIRNOFF99Frosst does similarly to both GAFF force fields, despite having fewer parameters. The necessity of adding each “type” of parameter for each atom type means that many parameters in GAFF are duplicates. But it is not clear which parameters can be pruned. SMIRNOFF99Frosst is a terse repesentation of a GAFF-like force field that is a good starting point for future development and optimization. How are we going to take this work forward?

Supporting Information

References