Combining the polarizable Drude force field with a continuum electrostatic Poisson-Boltzmann implicit solvation model

Abstract

In this work, we have combined the polarizable force field based on the classical Drude oscillator with a continuum Poisson-Boltzmann/solvent-accessible surface area (PB/SASA) model. In practice, the positions of the Drude particles experiencing the solvent reaction field arising from the fixed charges and induced polarization of the solute must be optimized in a self-consistent manner. Here, we parameterized the model to reproduce experimental solvation free energies of a set of small molecules. The model reproduces well-experimental solvation free energies of 70 molecules, yielding a root mean square difference of 0.8 kcal/mol versus 2.5 kcal/mol for the CHARMM36 additive force field. The polarization work associated with the solute transfer from the gas-phase to the polar solvent, a term neglected in the framework of additive force fields, was found to make a large contribution to the total solvation free energy, comparable to the polar solute-solvent solvation contribution. The Drude PB/SASA also reproduces well the electronic polarization from the explicit solvent simulations of a small protein, BPTI. Model validation was based on comparisons with the experimental relative binding free energies of 371 single alanine mutations. With the Drude PB/SASA model the root mean square deviation between the predicted and experimental relative binding free energies is 3.35 kcal/mol, lower than 5.11 kcal/mol computed with the CHARMM36 additive force field. Overall, the results indicate that the main limitation of the Drude PB/SASA model is the inability of the SASA term to accurately capture non-polar solvation effects. © 2018 Wiley Periodicals, Inc.

Keywords: CHARMM; Drude force field; Poisson-Boltzmann continuum solvation model; binding free energy; electronic polarization; implicit solvent model; molecular dynamics; protein-protein interactions.

Figure 1

Experimental solvation free energies vs. computed solvation free energies. Calculations were done using the C36 force field, shown in the top panels and Drude force field shown in the bottom panels. Left panels: the computed ΔG_solv with the initial vdW radii. Right panel: the computed ΔG_solv with the radii from the work of Nina et al for C36 and radii optimized in this work for the Drude calculations. The SASA tension coefficient was set to 5 cal/mol/Ų for the calculations with the vdW radii and 7.64 cal/mol/Ų for the calculations with the optimized set of radii and the Drude force field. Points for charged molecules are in red; points for polar molecules are in blue; and points for non-polar ones are in black. The gray line shows the perfect match between experimental and computed solvation free energies.

Figure 2

Induced polarization in implicit and explicit models. Left: Displacement of the Drude particles in the BPTI protein in implicit PB and explicit solvents relative to their positions in the protein in vacuum. Right: Angle between displacements of Drude particles induced by solvent interactions in the implicit PB and explicit solvents.

Figure 3

Computed and experimental binding free energy differences for protein mutants. Computations were done using (left) CHARMM36 and (right) Drude force fields. In both cases the protein dielectric constant of 1 was used.

Scheme 1

Thermodynamic cycle used to compute binding free energies in solvent. Superscripts show the medium in which the electronic polarization (position of Drude particles) of the protein atoms is computed. S indicates aqueous solution.