Journal of Theoretical and Computational Chemistry

Electrostatic solvation modeling based upon the Poisson–Boltzmann equation is widely used in studies of biomolecular structures and functions. This manuscript provides a thorough review of published efforts to adapt the numerical Poisson–Boltzmann methods to molecular simulations so that these methods can be extended to biomolecular studies involving conformational fluctuation and/or dynamics. We first review the fundamental works on how to define the electrostatic free energy and the Maxwell stress tensor. These topics are followed by three different strategies in developing algorithms to compute electrostatic forces and how to improve their numerical performance. Finally procedures are also presented in detail on how to discretize these algorithms for numerical calculations. Given the pioneer works reviewed here, further developmental efforts will be on how to balance efficiency and accuracy in these theoretical sound approaches — two important issues in applying any numerical algorithms for routine biomolecular applications. Even if not reviewed here, more advanced numerical solvers are certainly necessary to achieve higher accuracy than the widely used classical methods to improve the overall performance of the numerical Poisson–Boltzmann methods.

Methods for computing electrostatic interactions often account implicitly for the solvent, due to the much smaller number of degrees of freedom involved. In the Poisson–Boltzmann (PB) approach the electrostatic potential is obtained by solving the Poisson–Boltzmann equation (PBE), where the solvent region is modeled as a homogeneous medium with a high dielectric constant. PB however is not exempt of problems. It does not take into account for example the sizes of the ions in the atmosphere surrounding the solute, nor does it take into account the inhomogeneous dielectric response of water due to the presence of a highly charged surface. In this paper we review two major modifications of PB that circumvent these problems, namely the size-modified PB (SMPB) equation and the Dipolar Poisson–Boltzmann Langevin (DPBL) model. In SMPB, steric effects between ions are accounted for with a lattice gas model. In DPBL, 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. This model results in a dielectric profile that transits smoothly from the solute to the solvent region as well as in a variable solvent density that depends on the charges of the solute. We show successful applications of the DPBL formalism to computing the solvation free energies of isolated ions in water. Further developments of more accurately modified PB models are discussed.

Traditional implicit methods for modeling electrostatics in biomolecules use a two-dielectric approach: a biomolecule is assigned low dielectric constant while the water phase is considered as a high dielectric constant medium. However, such an approach treats the biomolecule–water interface as a sharp dielectric border between two homogeneous dielectric media and does not account for inhomogeneous dielectric properties of the macromolecule as well. Recently we reported a new development, a smooth Gaussian-based dielectric function which treats the entire system, the solute and the water phase, as inhomogeneous dielectric medium (Li L, Li C, Zhang Z, Alexor E, On the dielectric constant of proteins: Smooth dielectric function for macromolecular modeling and its implementation in Delphi, J Chem Theory Comput9(4):2126–2136, 2013). Here we examine various aspects of the modeling of polar solvation energy in such inhomogeneous systems in terms of the solute–water boundary and the inhomogeneity of the solute in the absence of water surrounding. The smooth Gaussian-based dielectric function is implemented in the DelPhi finite-difference program, and therefore the sensitivity of the results with respect to the grid parameters is investigated, and it is shown that the calculated polar solvation energy is almost grid independent. Furthermore, the results are compared with the standard two-media model and it is demonstrated that on average, the standard method overestimates the magnitude of the polar solvation energy by a factor 2.5. Lastly, the possibility of the solute to have local dielectric constant larger than of a bulk water is investigated in a benchmarking test against experimentally determined set of pKa's and it is speculated that side chain rearrangements could result in local dielectric constant larger than 80.

Protein phosphorylation is a very common post-translational modification, catalyzed by kinases, for signaling and regulation. Phosphotyrosines frequently target SH2 domains. The spleen tyrosine kinase (Syk) is critical for tyrosine phosphorylation of multiple proteins and for regulation of important pathways. Phosphorylation of both Y342 and Y346 in Syk linker B is required for optimal signaling. The SH2 domains of Vav1 and PLC-γ both bind this doubly phosphorylated motif. Here we used a recently developed method to calculate the effects of Y342 and Y346 phosphorylation on the rate constants of a peptide from Syk linker B binding to the SH2 domains of Vav1 and PLC-γ. The predicted effects agree well with experimental observations. Moreover, we found that the same doubly phosphorylated peptide binds the two SH2 domains via distinct mechanism, with apparent rigid docking for Vav1 SH2 and dock-and-coalesce for PLC-γ SH2.

Equilibrium peptide conformations in solution, especially in the presence of salts, has been of interest for several decades. The fundamental interactions that determine the dominant peptide conformations in solution have been experimentally and computationally probed; however, a unified understanding has not yet emerged. In a previous study, we performed metadynamics simulations on the heptapeptide AEAAAEA in Sodium Chloride (NaCl) and Potassium Chloride (KCl) solutions at concentrations ranging from 0.5–2.0 M. Using a three-dimensional collective variable coordinate system, we computed the free energy landscapes in each saline environment as well as in pure water. We found that the presence of Na⁺ and K⁺ ions induces some changes in the stability of the conformers that define the state space, but does not alter the overall energetics between conformers and does not favor helical conformations. We investigate here, how the presence of salts (NaCl and KCl) affects the glutamic–glutamic interaction and its consequences on the stability of each equilibrium conformation. We perform this study through fixed backbone simulations for the most populated conformations identified in our previous work: the α-helix, 3₁₀-helix, π-helix, the extended polyproline II (PPII) and 2.5₁-helix conformations. It was found that for each conformation, there exists stable substates determined by the glutamic acid side-chains distance and orientation, and that Na⁺ and K⁺ cations (de)stabilize preferentially each conformation. It was also found that intramolecular single water mediated hydrogen bonds play a crucial role in the observed (de) stabilization of each equilibrium conformation.

In this paper, molecular dynamics (MD) simulations were performed for a number of benchmark proteins using both the standard assisted model building with energy refinement (AMBER) charge and the dynamically adjusted polarized protein-specific charge (DPPC) from quantum fragment calculations to provide accurate electrostatic interactions. Our result shows that proteins' dynamic structures drifted away from the native structures in simulations under standard (nonpolarizable) AMBER force field. For comparison, proteins' native structures were dynamically stable after a long time simulation under DPPC. The free energy landscape reveals that the native structure is the lowest energy conformation under DPPC, while it is not under standard AMBER charge. To further investigate the polarization effect on the stability of native structures of proteins, we restarted from some decoy structures generated from simulations using standard AMBER charges and then carried out further MD simulation using DPPC to refine those structures. Our study shows that the native structures from these decoy structures can be mostly recovered using DPPC and that the dynamic structures with the highest population in cluster analysis are in close agreement with the corresponding native structures. The current study demonstrates the importance of electronic polarization of protein in stabilizing the native structure.

We investigate the influence of three common definitions of the solute/solvent dielectric boundary (DB) on the accuracy of the electrostatic solvation energy ΔG_el computed within the Poisson–Boltzmann (PB) and the generalized Born (GB) models of implicit solvation. The test structures include small molecules, peptides and small proteins; explicit solvent ΔG_el are used as accuracy reference. For common atomic radii sets BONDI, PARSE (and ZAP9 for small molecules) the use of van der Waals (vdW) DB results, on average, in considerably larger errors in ΔG_el than the molecular surface (MS) DB. The optimal probe radius ρ_w for which the MS DB yields the most accurate ΔG_el varies considerably between structure types. The solvent accessible surface (SAS) DB becomes optimal at ρ_w ~ 0.2Å(exact value is sensitive to the structure and atomic radii), at which point the average accuracy of ΔG_el is comparable to that of the MS-based boundary. The geometric equivalence of SAS to vdW surface based on the same atomic radii uniformly increased by ρ_w gives the corresponding optimal vdW DB. For small molecules, the optimal vdW DB based on BONDI +0.2 Å radii can yield ΔG_el estimates at least as accurate as those based on the optimal MS DB. Also, in small molecules, pairwise charge–charge interactions computed with the optimal vdW DB are virtually equal to those computed with the MS DB, suggesting that in this case the two boundaries are practically equivalent by the electrostatic energy criteria. In structures other than small molecules, the optimal vdW and MS dielectric boundaries are not equivalent: the respective pairwise electrostatic interactions in the presence of solvent can differ by upto 5 kcal/mol for individual atomic pairs in small proteins, even when the total ΔG_el are equal. For small proteins, the average decrease in pairwise electrostatic interactions resulting from the switch from optimal MS to optimal vdW DB definition can be mimicked within the MS DB definition by doubling of the solute dielectric constant. However, the use of the higher interior dielectric does not eliminate the large individual deviations between pairwise interactions computed within the two DB definitions. It is argued that while the MS-based definition of the DB is more physically correct in some types of practical calculations, the choice is not so clear in some other common scenarios.

A large fraction of proteins function as homodimers, but it is not always clear why the dimerization is important for functionality since frequently each monomer possesses a distinctive active site. Recent work (PLoS Computational Biology9(2):e1002924) indicates that homodimerization may be important for forming an electrostatic funnel in the spermine synthase homodimer which guides changed substrates toward the active centers. This prompted us to investigate the electrostatic properties of a large set of homodimeric proteins and resulted in an observation that in a vast majority of the cases the dimerization indeed results in specific electrostatic features, although not necessarily in an electrostatic funnel. It is demonstrated that the electrostatic dipole moment of the dimer is predominantly perpendicular to the axis connecting the centers of the mass of the monomers. In addition, the surface points with highest potential are located in the proximity of the interfacial plane of the homodimeric complexes. These findings indicate that frequent homodimerization provides specific electrostatic features needed for the function of proteins.

Continuum electrostatics modeling of solvation based on the Poisson–Boltzmann (PB) equation has gained wide acceptance in biomolecular applications such as energetic analysis and structural visualization. Successful application of the PB solvent models requires careful calibration of the solvation parameters. Extensive testing and validation is also important to ensure accuracy in their applications. Limitation in the continuum modeling of solvation is also a known issue in certain biomolecular applications. Growing interest in membrane systems has further spurred developmental efforts to allow inclusion of membrane in the PB solvent models. Despite their past successes due to careful parameterization, algorithm development and parallel implementation, there is still much to be done to improve their transferability from the small molecular systems upon which they were developed and validated to complex macromolecular systems as advances in technology continue to push forward, providing ever greater computational resources to researchers to study more interesting biological systems of higher complexity.

The homotetramer pyranose 2-oxidase (P2O) from Tetrametes multicolor contains flavin adenine dinucleotide (FAD) as a cofactor, and displays two conformers with different transient fluorescence spectra and lifetimes (ca. 0.1 ps and 360 ps). The ultrashort lifetimes of isoalloxazine (Iso) are ascribed to the photoinduced electron transfer (ET) from Trp168 to the excited Iso. Here, the structural heterogeneity among the four subunits in solution was studied by means of molecular dynamics simulation (MDS). The ET donor–acceptor distances in crystal and solution were compared. The distribution of the H-bond distances between Iso and the surrounding amino acids revealed appreciable differences among the four subunits. The structural fluctuations in two distant places were examined for the Iso-P and Iso-Q distances (where P and Q are Trp or Tyr) with the correlation coefficients between Iso-P and Iso-Q distances, revealing cooperative motions even though P and Q were more than 1 nm apart and located in different subunits. Moreover, distributions of the distances between Iso and its closest ionic amino acids markedly differed among the four subunits. Electrostatic (ES) energies between the Iso anion and the ionic amino acids in the entire protein were obtained using a static dielectric constant of 1. The ES energy in each subunit was strongly influenced by the other subunits, whilst the distributions of the ES energies greatly differed among the four subunits. This heterogeneous distribution of the ES energy between subunits may contribute to the large differences in the experimentally detected ET rates.