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@1374375 on May 9, 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.91 kcal/mol, 95% CI [0.71, 1.13] for SMIRNOFF99Frosst and 1.68 kcal/mol, 95% CI [1.51, 1.84] 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 Initiative. 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. Predictions of protein-ligand absolute binding free energies have achieved root mean square errors around 1-2 kcal/mol for “well-behaved” systems [1,2,3], with deviations an order of magnitude larger for some protein families [4]. Relative free energy calculations on a series of congeneric ligands, using proprietary methods, have also achieved root mean square errors ~1 kcal/mol [5,6]. A variety of techniques for computing absolute binding free energies in host-guest systems have shown accuracy of ~1 kcal/mol, as highlighted in the recent SAMPL5 and SAMPL6 blind challenges [1,7]. These techniques include both quantum and classical dynamics, employing a range of energy and solvation models, with some techniques having knowledge-based steps, docking, or clustering [10,11,12,13,14,15,8,9].

Give the above paragraph another look.

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 [16], GAFF v2.1, and SMIRNOFF99Frosst [17,18].

SMIRNOFF99Frosst, released in late 2018, is the first force field produced by the Open Force Field Initiative. SMIRNOFF99Frosst is derived from AMBER parm99 [19] and Merck’s parm@Frosst [20]. 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 [21]. 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 322 lines of 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 [17]. 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 (on both absolute errors and correlation) and GAFF v2.1 (on absolute errors), predicting experimental binding free energies, enthalpies, and entropies. We then characterize the differences in host conformations sampled 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) [22] and computationally in [23]. As in Henriksen, et al. [23], 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. [22], and the ΔG difference between each stereoisomer is expected to be < 0.1 kcal/mol [24].

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 [25,26] 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 [27] and Joung-Cheatham monovalent ion parameters [28] in each simulation set.

GAFF v1.7 bond, angle, torsion, and Lennard-Jones parameters were applied using the tleap program distributed with AmberTools16.

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, the toplogy and coordinates for the host-guest complex was combined with the solvent and ions, which retained their TIP3P water parameters and Joung-Cheatham ion parameters, respectively. This was accomplished by the ParmEd program [29], which enables saving the OpenMM system created by the Open Force Field Toolkit in AMBER-format .prmtop and .inpcrd files. 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 provided in the literature [12,30,31,32]. 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 [12] 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 with the GAFF v1.7 force field) and AMBER 18 (calculations with the GAFF v2.1 and SMIRNOFF99Frosst force fields) molecular dynamics software [33]. Each window for each system was independently solvated and simulated. The host-guest simulations using GAFF v1.7 were previously published in Henriksen, et al. [23] and are described in additional detail therein. 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 create a rectangular periodic box, reducing the amount of solvent required for the calculation. 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 [22]. In the GAFF simulations, hydrogen mass repartitioning [34] was used to adjust the mass of hydrogen by a factor of 3, keeping the total molecular weight of each molecule constant, 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. An isotropic analytic correction to the Lennard-Jones interactions is applied beyond the cutoff distance [35]. Production NPT simulations were run for a minimum of 2.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 uncertainty in the work done in each λ window. 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.

Excluding the first and last window, the average window length was 11.8 ns and 5.39 ns for GAFF v1.7 and GAFF v2.1 simulations, respectively. In the interest of simplicity, 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 the work done by each restraint in each simulation window, \(\left\langle \frac{\partial U}{\partial \lambda_i} \right\rangle_\text{SEM}\), was estimated using blocking analysis [36]. Based on these values, resampling with replacement was used to create 100,000 hypothetical \(\frac{\partial U}{\partial \lambda}\) curves. The standard deviation in the integrated work, \(\int \frac{\partial U}{\partial \lambda_i}\,\mathrm{d}\lambda_i\), of the 100,000 curves was used to determine the standard error of the mean of ΔG. The standard error of the mean of ΔH was computed from the standard error of the mean of the total potential energy in each end point window, estimated using blocking analysis, added in quadrature. The standard error of the mean −TΔS was calculated using the uncertainties in ΔG and ΔH added 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 uncertainties 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 [16]) 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

Binding free energies

Despite having far fewer numerical parameters, SMIRNOFF99Frosst does about as well as GAFF v1.7 and better than GAFF v2.1 on predicting ΔG compared to values measured using ITC or NMR. SMIRNOFF99Frosst has an overall deviation from experiment of 0.91 kcal/mol, 95% CI [0.71, 1.13] on binding free energies across the 43 host-guest systems, compared to 0.88 kcal/mol, 95% CI [0.71, 1.08] for GAFF v1.7 and 1.68 kcal/mol, 95% CI [1.51, 1.85] for GAFF v2.1 (Figures 3, 16, Table 2, Table 7, Table 8, and Table 9) On the whole, GAFF v1.7 agrees well with SMIRNOFF99Frosst (Figure 20); the overall root mean squared error (RMSE) and mean signed error (MSE) between the methods is 0.80 kcal/mol, 95% CI [0.58, 1.04] and -0.47 kcal/mol, 95% CI [-0.67, -0.28], respectively. Both SMIRNOFF99Frosst and GAFF v1.7 systematically underestimate the ΔG for cyclic alcohols. GAFF v2.1 significantly overestimates the binding of all compounds by more than 1.5 kcal/mol, with a mean signed error of -1.56 kcal/mol, 95% CI [-1.74, -1.37]. Despite this, GAFF v2.1 has a strikingly strong correlation with experiment across all functional group classes. This may be traced to differences in host conformations sampled by GAFF v2.1, which indicate a more consistently open cyclodextrin “pocket” for guests to bind (Figure 14).

Binding enthalpies and entropies

On both binding enthalpies and entropies, SMIRNOFF99Frosst and GAFF v1.7 agree reasonably well with each other (ΔH RMSE = 1.65 kcal/mol, 95% CI [1.32, 2.02], −TΔS RMSE = 1.59 kcal/mol, 95% CI [1.24, 1.96]) (Figure 20). The deviations between SMIRNOFF99Frosst and GAFF v2.1 are higher for ΔH (RMSE = 2.96 kcal/mol, 95% CI [2.52, 3.42]) and lower for −TΔS (RMSE = 1.55 kcal/mol, 95% CI [1.25, 1.85]) On binding enthalpies, SMIRNOFF99Frosst agrees the best with experiment (RMSE = 1.85 kcal/mol, 95% CI [1.40, 2.29]), followed by GAFF v2.1 (RMSE = 2.21 kcal/mol, 95% CI [1.77, 2.66]), and then GAFF v1.7 (RMSE = 2.54 kcal/mol, 95% CI [2.08, 2.99]), In some cases, GAFF v1.7 underestimates ΔH by over 3 kcal/mol and up to 5 kcal/mol (b-chp).

On binding entropies, GAFF v2.1 has the lowest RMSE compared to experiment (RMSE = 1.47 kcal/mol, 95% CI [1.09, 1.99]), followed by SMIRNOFF99Frosst (RMSE = 1.90 kcal/mol, 95% CI [1.49, 2.32]), and GAFF v1.7 (RMSE = 2.21 kcal/mol, 95% CI [1.74, 2.68]) (Figure 16).
All force fields perform poorly predicting −TΔS for carboxylate guests. It is worth noting that all force fields predict a much smaller entropic component of binding for a-coc by 3-5 kcal/mol, which does not easily fit inside the primary cavity of cyclodextrin (Figure 4).

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 (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 (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, 95% CI [-0.54, 0.30] and RMSE = 0.76 kcal/mol, 95% CI [0.43, 1.12]) 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 binds 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, 95% CI [0.22, 1.21] and RMSE = 1.07 kcal/mol, 95% CI [0.66, 1.58]) (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, 95% CI [-0.73, -0.01] and RMSE = 0.87 kcal/mol, 95% CI [0.58, 1.16]) (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 (Figures 24, 25, and 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. Note, that hydroxyl hydrogens are assigned σ = 0 Å and ε = 0 kcal/mol in both GAFF v1.7 and SMIRNOFF99Frosst v1.0.5, but later versions of SMIRNOFF99Frosst adopt small σ and ε values based on a similiar atom type in parm@Frosst [37,38,39]. 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 (Figure 2, GAFF atom types h2-c3-c3-oh, SMIRKS pattern [#1:1]-[#6X4:2]-[#6X4:3]-[#8X2:4]), 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). Partly this is due to the addition of dihedral terms with a barrier height of exactly 0.00 kcal/mol in GAFF, which are used to override wildcard parameters that might match the same atom types. For example, GAFF v2.1 applies a three term energy function to the atom types c3-os-c3-c3, whereas SMIRNOFF99Frosst employs a two term energy function for the SMIRKS pattern ([#6X4:1]-[#6X4:2]-[#8X2H0:3]-[#6X4:4]), but only the terms with periodicity 2 and 3 have nonzero barrier heights in GAFF v2.1. SMIRNOFF99Frosst uses two nonzero terms to model the potential barrier for the SMIRKS pattern [#6X4:1]-[#6X4:2]-[#8X2H1:3]-[#1:4] yet GAFF v2.1 applies a single term with a barrier height of exactly 0.00 kcal/mol for the atom types c3-c3-oh-ho.

I am considering removing this table. I think it is difficult to interpret: there are cases where things are missing, things are zero, and duplicate parameters applied to different sets of atoms and listed multiple times.

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 (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

Code and data availability

• GitHub repository used to convert AMBER input files from GAFF force field to SMIRNOFF99Frosst.
• GitHub repository for setting up the attach-pull-release calculations using paprika version 0.0.3.
• GitHub repository for analyzing the simulations and generating the plots in this manuscript.
• GitHub repository for the Open Force Field group containing the toolkit and force field XML file.

Author contributions

Conceptualization, DRS, NMH, JDC, MKG; Methodology, DRS, NMH; Software, DRS, NMH; Formal Analysis, DRS, NMH, JDC, MKG; Investigation, DRS, NMH; Resources, MKG, JDC; Data Curation, DRS, NMH; Writing-Original Draft, DRS, NMH; Writing - Review and Editing, DRS, NMH, JDC, MKG; Visualization, DRS; Supervision, JDC, MKG; Project Administration, MKG; Funding Acquisition, MKG.

Acknowledgments

Disclosures

The authors declare the following competing financial interest(s): MKG has an equity interest in and is a cofounder and scientific advisor of VeraChem LLC.

List of abbreviations

APR, attach-pull-release; CD, cyclodextrin; GAFF, Generalized AMBER Force Field

References