A generalized Kirkwood implicit solvent for the polarizable AMOEBA protein model

Abstract

Computational simulation of biomolecules can provide important insights into protein design, protein-ligand binding interactions, and ab initio biomolecular folding, among other applications. Accurate treatment of the solvent environment is essential in such applications, but the use of explicit solvents can add considerable cost. Implicit treatment of solvent effects using a dielectric continuum model is an attractive alternative to explicit solvation since it is able to describe solvation effects without the inclusion of solvent degrees of freedom. Previously, we described the development and parameterization of implicit solvent models for small molecules. Here, we extend the parameterization of the generalized Kirkwood (GK) implicit solvent model for use with biomolecules described by the AMOEBA force field via the addition of corrections to the calculation of effective radii that account for interstitial spaces that arise within biomolecules. These include element-specific pairwise descreening scale factors, a short-range neck contribution to describe the solvent-excluded space between pairs of nearby atoms, and finally tanh-based rescaling of the overall descreening integral. We then apply the AMOEBA/GK implicit solvent to a set of ten proteins and achieve an average coordinate root mean square deviation for the experimental structures of 2.0 Å across 500 ns simulations. Overall, the continued development of implicit solvent models will help facilitate the simulation of biomolecules on mechanistically relevant timescales.

FIG. 1.

Pictorial representation of effective radii for an arbitrary globular molecule in implicit water. The effective radii, a_i and a_j, for atoms i and j, respectively, are larger than the intrinsic radii for both atoms (ρ_i and ρ_j). Since atom j is more deeply buried within the molecule than atom i, a_j is larger than a_i. The effective radius for ion k is approximately equivalent to its intrinsic radius ρ_k but could be slightly larger due to descreening by the nearby molecule.

FIG. 2.

Shown is the neck region (shaded) between two atoms. A neck is only formed when the atoms are close enough to exclude water, which is represented here by a sphere with radius ρ_w.

FIG. 3.

Shown are the bonding aware neck scaling scheme and associated equations used to calculate an individual atomic S_neck,i scale factor. The scaling factor is reduced as the number of heavy atoms bound to the atom of interest increases. If the atom of interest has no bound heavy atoms, the scale factor is set to 1.0.

FIG. 4.

Demonstration of the bonding aware neck scaling scheme for two ions in the AMOEBA GK implicit solvent. Using a single, fixed scaling factor leads to underestimation of the neck integral value at all separation distances for pairs of ions, while chemically aware scaling concords with reference PB calculations.

FIG. 5.

Comparison of permanent electrostatic energies for biomolecules calculated in APBS and GK using a van der Waals definition of the solute volume. GK values used element specific HCT scaling factors. The dashed gray line is the best fit regression line with slope = 1.0001 and R² = 0.9971. The solid black is x = y to guide the eye.

FIG. 6.

Effective radii for a molecule with a 1.7 Å base radius (carbon atom) without rescaling (dashed line) and with the final tanh rescaling function based on β₀ = 0.9563, β₁ = 0.2578, and β₂ = 0.0810 (solid line) across a range of input volume integrals (Ψ).

FIG. 7.

Comparison of permanent self-energies with perfect effective Kirkwood radii and GK fit radii for proteins. Energies calculated with only the tanh correction have a slope of 1.192 and a R² of 0.997; energies calculated with the full correction (neck and tanh) have an improved slope of 0.987 and a R² of 0.998. The solid black line is x = y to guide the eye.

FIG. 8.

Comparison of 2TRX effective radii calculated in GK without interstitial space corrections (orange plusses) and with interstitial space corrections (blue circles) to perfect effective Kirkwood radii calculated in APBS. Fit GK base radii are used for both series. The first plot (left) shows all effective radii for 2TRX (thioredoxin), while the second plot (right) shows the range of effective radii from 1 to 5 Å. The solid black line is x = y to guide the eye.

FIG. 9.

Electrostatic hydration energy of target proteins was calculated using APBS and GK, including contributions from permanent multipoles and induced dipoles determined via a self-consistent reaction field. APBS used a molecular surface based on fit GK base radii. GK also used fit radii, with effective Born radii computed with (Slope: 0.987, R²: 0.9937) and without (Slope: 1.575, R²: 0.9907) interstitial space corrections.

FIG. 10.

Superposition of the deposited x-ray crystallography or NMR structure (gray) with the lowest-RMSD structure from the largest cluster (green). The time step of the snapshot in the 500 ns trajectory and its RMSD to the experimental structure are displayed. All representative snapshots were taken from the second half of the trajectories.

FIG. 11.

Protein non-terminal backbone RMSDs to the experimental x-ray crystallography or NMR structure across 500 ns production trajectories (Å).

FIG. 12.

Histograms of RMSD values across 500 ns of simulation for AMOEBA/GK implicit solvent are shown for each protein.

FIG. 13.

Comparison of average dipole moment magnitudes across production MD trajectories in a vacuum and GK implicit solvent. Implicit solvent dipole moment magnitudes are reported with (blue circles) and without (orange pluses) interstitial space corrections. The dashed lines are the best fit regression lines for GK dipole moment magnitudes, while the solid black line is x = y to guide the eye.