Implicit Solvents for the Polarizable Atomic Multipole AMOEBA Force Field

Abstract

Computational protein design, ab initio protein/RNA folding, and protein-ligand screening can be too computationally demanding for explicit treatment of solvent. For these applications, implicit solvent offers a compelling alternative, which we describe here for the polarizable atomic multipole AMOEBA force field based on three treatments of continuum electrostatics: numerical solutions to the nonlinear and linearized versions of the Poisson-Boltzmann equation (PBE), the domain-decomposition conductor-like screening model (ddCOSMO) approximation to the PBE, and the analytic generalized Kirkwood (GK) approximation. The continuum electrostatics models are combined with a nonpolar estimator based on novel cavitation and dispersion terms. Electrostatic model parameters are numerically optimized using a least-squares style target function based on a library of 103 small-molecule solvation free energy differences. Mean signed errors for the adaptive Poisson-Boltzmann solver (APBS), ddCOSMO, and GK models are 0.05, 0.00, and 0.00 kcal/mol, respectively, while the mean unsigned errors are 0.70, 0.63, and 0.58 kcal/mol, respectively. Validation of the electrostatic response of the resulting implicit solvents, which are available in the Tinker (or Tinker-HP), OpenMM, and Force Field X software packages, is based on comparisons to explicit solvent simulations for a series of proteins and nucleic acids. Overall, the emergence of performative implicit solvent models for polarizable force fields opens the door to their use for folding and design applications.

Figure 1.

The total implicit solvent potential of mean force ΔW_hydration(X) can be formulated using a thermodynamic cycle composed of five steps. Step 1: the solute is decharged in $\mathrm{vacuum}-U_{\mathit{elec}}^{\mathit{vacuum}}(\mathbf{X})$. Step 2: dispersion interactions are removed between the solute and surrounding medium, which has no energetic cost in vacuum. Step 3: a solute-shaped cavity is fonned in water ΔW_cav(X), which is unfavorable and proportional to solvent excluded volume for small solutes. Step 4: favorable solute-solvent dispersion interactions are added ΔW_disp(X). Step 5: the solute is charged in solvent $U_{\mathit{elec}}^{\mathit{water}}(\mathbf{X})$ to yield an overall electrostatic contribution of $\Delta {\text{W}}_{\mathrm{elec}}(\mathbf{X})=U_{\mathit{elec}}^{\mathit{water}}(\mathbf{X})-U_{\mathit{elec}}^{\mathit{vacuum}}(\mathbf{X})$.

Figure 2.

The total dipole moment for the parameterization set of 103 molecules from QM using MP2/6-311G**(2d,2p) are compared to those of the resulting AMOEBA models (red dashed line: Y = 0.9986 · X +0.0011; R² = 0.9999).

Figure 3.

Surface tension is not constant for small solutes, but increases approximately linearly until the effective radius of the solute grows to beyond ~10 Å. For large (flat) solutes, surface tension asymptotes toward the experimental value for a water-vapor interface of 0.103 kcal/mol/Ų. This length scale dependence can be approximately captured by a cavitation free energy difference that switches between using SEV and SASA via either a simple, non-differentiable form (ΔW_sphere(r) given by Equation 6, black dashes) or the smooth form used in this work (ΔW_cav (r) given by Equation 11, green solid line). The asymptotic surface tension of ΔW_cav (r) can be reduced relative to the experimental value (e.g. to 0.08 kcal/mol/Ų, dashed blue line) to capture cavitation for solutes with a large effective radius, but which are more highly curved than a simple sphere (e.g. a DNA double helix, RNA molecule, or protein).

Figure 4.

The pairwise buffered 14-7 potential (U_{14 – 7}) can be decomposed into purely repulsive (U_rep, Equation 13) and attractive (U_WCA, Equation 16) contributions, which are plotted for the AMOEBA water oxygen atom (r_0,ij = 3.405 Å, ε_ij = 0.11 kcal/mol). The cavitation free energy (ΔW_cav) represents the process of growing in the repulsive potential (U_rep) for solute atoms. Subsequently, the dispersion free energy (ΔW_disp) models the process of adding the attractive U_WCA solute-solvent interactions to recover the full U_{14 – 7} potential in the context of an uncharged solute.

Figure 5.

Panel A. Dispersion free energy differences (ΔW_disp, solid black curve) are given by the integral of the attractive WCA potential (U_WCA, dotted blue curve) over solvent for the interaction of two AMOEBA water oxygen atoms. The derivative of ΔW_disp with respect to R₀ (dashed green curve) shows a maximum slightly beyond the minimum energy separation distance (vertical red line) due to the volume element 4πr²dr increasing more quickly than U_WCA approaches zero just beyond r_0,ij. Panel B. Dispersion interactions of an explicit water oxygen with continuum water oxygen (also plotted in Panel A) and hydrogen. The interaction of oxygen with continuum hydrogen is multiplied by 2 (two hydrogen per water). Panel C. Same as Panel B, but for two explicit water hydrogen interacting with continuum water.

Figure 6.

Illustration of the integration limits for the dispersion free energy based on solute van der Waals parameters and separation distance. Panel A: The minimum separation distance r_0,ij for atoms i and j is based on AMOEBA mixing rules and used to determine the beginning of the WCA dispersion integral R₀. When R₀ < r_0,ij, integration for the constant portion of the WCA potential begins at either R₀ or r_ij − ρ_j, whichever is larger, and ends at r_0,ij or r_ij + ρ_i, whichever is smaller. If R₀ > r_0,ij, then only the variable portion of WCA potential factors into the dispersion free energy. Panel B: If atom j is completely engulfed by the sphere defined by R₀, no solvent is blocked and no dispersion energy must be removed. Panel C: When the two atoms overlap or are close together such that r_ij − ρ_j < R₀, integration of the attractive WCA potential begins at R₀ and ends at r_ij + ρ_j (the furthest edge of atom j). Panel D: When atom j is outside the beginning of the integration, integration of the variable portion of the WCA potential begins at r_ij − ρ_j (the closest edge of atom j) and ends at r_ij + ρ_j.

Figure 7.

Panel A: Fit of GK self energies to perfect PB self energies (Y = 1.040X + 0.048, R² = 0.996). Panel B: Fit of GK cross term energies to perfect PB cross tenn energies (Y = 1.059X − 0.089, R² = 0.994). Both self and aggregate cross-term energies are reported for 1424 atoms.

Figure 8.

Shown is a comparison of experimental and APBS solvation free energy differences using either AMOEBA van der Waals R_min radii to describe the solute-solvent boundary (Panel A) or using fit radii (Panel B). When using AMOEBA R_min radii, the linear regression gave Y = 0.9683 · X + 0.3215 with R² = 0.9702. When using fit radii, the linear regression gave Y = 1.0080 · X + 0.1416 with R² = 0.9979.

Figure 9.

Shown is a comparison of experimental and ddCOSMO solvation free energy differences using either AMOEBA van der Waals R_min radii to describe the solute-solvent boundary (Panel A) or using fit radii (Panel B). When using AMOEBA R_min radii, the linear regression gave Y = 1.0017 · X − 0.1465 with R² = 0.9738. When using fit radii, the linear regression gave Y = 1.0001 · X + 0.0015 with R² = 0.9981.

Figure 10.

Shown is a comparison of experimental and GK solvation free energy differences using either AMOEBA van der Waals R_min radii to describe the solute-solvent boundary (Panel A) or using fit radii (Panel B). When using AMOEBA R_min radii, the linear regression gave Y = 0.9245 · X + 0.0360 with R² = 0.9567. When using fit radii, the linear regression gave Y = 0.9999 · X − 0.0046 with R² = 0.9987.

Figure 11.

The validation set includes nine nucleic acids and nine proteins. The nucleic acid set can be further broken down into sets of four RNA helices (2JXQ, 1F5G, 1MIS, and 2L8F), three RNA hair pins (1ZIH, 2KOC, and 1SZY), and two DNA double helices (1D20 and 2HKB). Additional information on individual molecules and simulation conditions can be found in Supplementary Table S10. Coordinate files for each molecule (Tinker “XYZ” format) are provided in the Supporting Information.

Figure 12.

The APBS energy values for the biomolecular validation set are plotted vs. ddCOSMO and GK for both van der Waals radii and fit radii. Panel A. APBS electrostatics energy using van der Waals radii vs ddCOSMO (Y = 1.0523X, R² = 0.9990) and GK (Y = 0.9854X, R² = 0.9999). Panel B. APBS PBE electrostatics energy using fit radii vs ddCOSMO (Y = 0.9884X, R² = 0.9991) and GK (Y = 0.9833X, R² = 0.9991). Full data is available in Supplementary Tables S11 and S13.

Figure 13.

Comparison of dipole moment magnitudes for fixed charge force fields and AMOEBA across environments. Panel A: AMOEBA vacuum dipole moment magnitudes vs those for fixed charge force fields (AMBER ff99SB R²: 0.987, OPLS-AA/L R²: 0.979, CHARMM22/CMAP R²: 0.982). Panel B: AMOEBA explicit solvent dipole moment magnitudes vs those for fixed charge force fields (AMBER ff99SB R²: 0.998, OPLS-AA/L R²: 0.996, CHARMM22/CMAP R²: 0.996). Dotted lines at x = y are to guide the eye.

Figure 14.

Comparison of dipole moment magnitudes for AMOEBA in explicit solvent vs. vacuum, PB, ddCOSMO and GK environments for the validation set (vacuum R²: 0.990, PB/ddCOSMO/GK R²: 0.999 in each case). The dotted line at x = y is to guide the eye.