Incorporating dipolar solvents with variable density in Poisson-Boltzmann electrostatics

Abstract

We describe a new way to calculate the electrostatic properties of macromolecules that goes beyond the classical Poisson-Boltzmann treatment with only a small extra CPU cost. The solvent region is no longer modeled as a homogeneous dielectric media but rather as an assembly of self-orienting interacting dipoles of variable density. The method effectively unifies both the Poisson-centric view and the Langevin Dipole model. The model results in a variable dielectric constant epsilon(r) in the solvent region and also in a variable solvent density rho(r) that depends on the nature of the closest exposed solute atoms. The model was calibrated using small molecules and ions solvation data with only two adjustable parameters, namely the size and dipolar moment of the solvent. Hydrophobicity scales derived from the solvent density profiles agree very well with independently derived hydrophobicity scales, both at the atomic or residue level. Dimerization interfaces in homodimeric proteins or lipid-binding regions in membrane proteins clearly appear as poorly solvated patches on the solute accessible surface. Comparison of the thermally averaged solvent density of this model with the one derived from molecular dynamics simulations shows qualitative agreement on a coarse-grained level. Because this calculation is much more rapid than that from molecular dynamics, applications of a density-profile-based solvation energy to the identification of the true structure among a set of decoys become computationally feasible. Various possible improvements of the model are discussed, as well as extensions of the formalism to treat mixtures of dipolar solvents of different sizes.

FIGURE 1

The PBL model and its lattice representation. The atmosphere around the solute includes free ions of charge +ze or −ze and dipoles representing the water. Each lattice point is associated with a spinlike variable indicating the presence of a free ion (S_j) or a spherical dipole (d_j) of the same size. The grid size (a) used to place these entities and ensure exclusion effects is not necessarily the same as the grid size (h) used for discretizing the free energy equation. The interface between the solute interior and the solvent lattice is defined as either the accessible surface (using R_probe = 1.4 Å) or the molecular surface (using R_shrink = 1.0 Å).

FIGURE 2

Calibration of the magnitude p₀ and size a of the solvent dipole using small molecules and ions. PB versus PBL energy for ions of different charges (top left) and radii (top right) for different values of p₀ with a/2 = 1.4 Å. Same thing for the 16 uncharged amino acids (bottom left), and the remaining four charged amino acids together with a set of charged small molecules (bottom right). This includes acetate, dimethylphosphate, ethoxyde, ethylthiolate, guanidinium, imidazolium, methylammonium, methylthiolate, methoxyde, and tetramethylammonium.

FIGURE 3

Comparison of electrostatic potentials and fields for the PB and PBL methods using 2ACS. Correlation between the PB and GBPL electrostatic potentials inside the protein solute (a) and outside, i.e., in the solvent region (b). (c) Ratio of the magnitudes of the electric field E found with either PBE or PBLE, as a function of the distance to the solute. This is especially important as it is the electric field that drives the solvent density (see Eq. 9). Here we have a/2 = 1.2 Å and p₀ = 1.85D; the same kind of curves was observed with increasing values of p₀ up to 4.8D and/or choosing a/2 = 1.4 Å.

FIGURE 4

Ion condensation. Number of negative (left panel) and positive (right panel) ions placed as a function of threshold density ρ_thres for molecules of different net charges. The result of the same calculation with p₀ = 0 is indicated with a dotted line. (a) Aldose reductase (2ACS, net charge = −2); (b) halophilic alcohol dehydrogenase (2B5V, net charge = −33); and (c) B-DNA (9BNA, 12 bp, net charge = −24).

FIGURE 5

Dielectric properties of the solvent (2ACS). (a) Correlation through space of as a function of the distance away from the solute's surface, for lattice points separated by increasing distances (|r − r′| = 1,2,3,4 lattice points; in this case the grid spacing is 0.7 Å). (b) Profile of the dielectric function in the solvent region as a function of the distance to the center of gravity of the solute for p₀ = 1.85D, p₀ = 2.35D, or p₀ = 4.8D.

FIGURE 6

Radial solvent density profiles for different atom types at p₀ = 1.85D. This was obtained with a solute surface defined as the accessible surface area. Here the value a/2 = 1.2 Å was used; increasing this value to a/2 = 1.4 Å will lower the peaks of the solvent density, while increasing p₀ will increase them.

FIGURE 7

Hydrophobicity scales. (Top) Density profile-derived hydrophobicity scales for the 20 amino acids using either the maximum value of the curve (left) or the integral of its squared derivative (right). (Bottom) Mean values of the atomic surface tension coefficients with standard deviations at p₀ = 1.85D (left) and comparison with the Eisenberg and MacLachlan (61) surface tension coefficients (right). In the latter case, both the p₀ = 1.85D or p₀ = 4.8D curves are shown.

FIGURE 8

Condensed water molecules in the first hydration shell. (a) Number of water molecules placed, as a function of the density threshold (upper curve). The lower curve is the derivative of the upper curve together with a tentative fit of the lower curve with a Gaussian function. (b) Number of water molecules placed using different methods: a constant hydration value (h = 0.4 g/g), a density threshold criterion (1.2 or 1.5 or the density set to the maximum derivative of the N_w = f(ρ_thres) curve in (a)), integration of the excess density and a layer of constant size (1 Å) around the solute (null hypothesis). The dependence on the size (number of amino acids) of the solute proteins is shown.

FIGURE 9

Lysozyme decoys derived from circular permutations of the sequence. (a) Gradient-based partial solvent energy F₁. (b) Solvent energy based on the proportion of poorly solvated SAS. (c) Z-scores derived from local atomic surface tension coefficients. In all cases the correct structure is the last one and its energy is indicated by a continuous line. The mean value of the energy is materialized by a continuous line bracketed by two dotted lines (±1 SD).

FIGURE 10

1CTF decoys. (a) Gradient-based partial solvent energy F₁ as a function of the decoy structure number. The correct structure is the last one and its energy is indicated by a continuous line. The mean value of the energy is materialized by a continuous line bracketed by two dotted lines (±1 SD). (b) Partial solvent energy F₁ as a function of the RMSD with the correct model. The native structure energy (RMSD = 0, circle) is shown with a dotted line.