Continuum description of solvent dielectrics in molecular-dynamics simulations of proteins

Linearly scaling and almost Hamiltonian dielectric continuum molecular dynamics simulations through fast multipole expansions

Hamiltonian Dielectric Solvent (HADES) is a recent method [S. Bauer et al., J. Chem. Phys. 140, 104103 (2014)] which enables atomistic Hamiltonian molecular dynamics (MD) simulations of peptides and proteins in dielectric solvent continua. Such simulations become rapidly impractical for large proteins, because the computational effort of HADES scales quadratically with the number N of atoms. If one tries to achieve linear scaling by applying a fast multipole method (FMM) to the computation of the HADES electrostatics, the Hamiltonian character (conservation of total energy, linear, and angular momenta) may get lost. Here, we show that the Hamiltonian character of HADES can be almost completely preserved, if the structure-adapted fast multipole method (SAMM) as recently redesigned by Lorenzen et al. [J. Chem. Theory Comput. 10, 3244-3259 (2014)] is suitably extended and is chosen as the FMM module. By this extension, the HADES/SAMM forces become exact gradients of the HADES/SAMM energy. Their translational and rotational invariance then guarantees (within the limits of numerical accuracy) the exact conservation of the linear and angular momenta. Also, the total energy is essentially conserved-up to residual algorithmic noise, which is caused by the periodically repeated SAMM interaction list updates. These updates entail very small temporal discontinuities of the force description, because the employed SAMM approximations represent deliberately balanced compromises between accuracy and efficiency. The energy-gradient corrected version of SAMM can also be applied, of course, to MD simulations of all-atom solvent-solute systems enclosed by periodic boundary conditions. However, as we demonstrate in passing, this choice does not offer any serious advantages.

Parameterization of the Hamiltonian Dielectric Solvent (HADES) Reaction-Field Method for the Solvation Free Energies of Amino Acid Side-Chain Analogs

Optimization of the Hamiltonian dielectric solvent (HADES) method for biomolecular simulations in a dielectric continuum is presented with the goal of calculating accurate absolute solvation free energies while retaining the model's accuracy in predicting conformational free-energy differences. The solvation free energies of neutral and polar amino acid side-chain analogs calculated by using HADES, which may optionally include nonpolar contributions, were optimized against experimental data to reach a chemical accuracy of about 0.5 kcal mol(-1). The new parameters were evaluated for charged side-chain analogs. The HADES results were compared with explicit-solvent, generalized Born, Poisson-Boltzmann, and QM-based methods. The potentials of mean force (PMFs) between pairs of side-chain analogs obtained by using HADES and explicit-solvent simulations were used to evaluate the effects of the improved parameters optimized for solvation free energies on intermolecular potentials.

Electrostatics of proteins in dielectric solvent continua. I. An accurate and efficient reaction field description

We present a reaction field (RF) method which accurately solves the Poisson equation for proteins embedded in dielectric solvent continua at a computational effort comparable to that of an electrostatics calculation with polarizable molecular mechanics (MM) force fields. The method combines an approach originally suggested by Egwolf and Tavan [J. Chem. Phys. 118, 2039 (2003)] with concepts generalizing the Born solution [Z. Phys. 1, 45 (1920)] for a solvated ion. First, we derive an exact representation according to which the sources of the RF potential and energy are inducible atomic anti-polarization densities and atomic shielding charge distributions. Modeling these atomic densities by Gaussians leads to an approximate representation. Here, the strengths of the Gaussian shielding charge distributions are directly given in terms of the static partial charges as defined, e.g., by standard MM force fields for the various atom types, whereas the strengths of the Gaussian anti-polarization densities are calculated by a self-consistency iteration. The atomic volumes are also described by Gaussians. To account for covalently overlapping atoms, their effective volumes are calculated by another self-consistency procedure, which guarantees that the dielectric function ε(r) is close to one everywhere inside the protein. The Gaussian widths σ(i) of the atoms i are parameters of the RF approximation. The remarkable accuracy of the method is demonstrated by comparison with Kirkwood's analytical solution for a spherical protein [J. Chem. Phys. 2, 351 (1934)] and with computationally expensive grid-based numerical solutions for simple model systems in dielectric continua including a di-peptide (Ac-Ala-NHMe) as modeled by a standard MM force field. The latter example shows how weakly the RF conformational free energy landscape depends on the parameters σ(i). A summarizing discussion highlights the achievements of the new theory and of its approximate solution particularly by comparison with so-called generalized Born methods. A follow-up paper describes how the method enables Hamiltonian, efficient, and accurate MM molecular dynamics simulations of proteins in dielectric solvent continua.

Electrostatics of proteins in dielectric solvent continua. II. Hamiltonian reaction field dynamics

In Paper I of this work [S. Bauer, G. Mathias, and P. Tavan, J. Chem. Phys. 140, 104102 (2014)] we have presented a reaction field (RF) method, which accurately solves the Poisson equation for proteins embedded in dielectric solvent continua at a computational effort comparable to that of polarizable molecular mechanics (MM) force fields. Building upon these results, here we suggest a method for linearly scaling Hamiltonian RF/MM molecular dynamics (MD) simulations, which we call "Hamiltonian dielectric solvent" (HADES). First, we derive analytical expressions for the RF forces acting on the solute atoms. These forces properly account for all those conditions, which have to be self-consistently fulfilled by RF quantities introduced in Paper I. Next we provide details on the implementation, i.e., we show how our RF approach is combined with a fast multipole method and how the self-consistency iterations are accelerated by the use of the so-called direct inversion in the iterative subspace. Finally we demonstrate that the method and its implementation enable Hamiltonian, i.e., energy and momentum conserving HADES-MD, and compare in a sample application on Ac-Ala-NHMe the HADES-MD free energy landscape at 300 K with that obtained in Paper I by scanning of configurations and with one obtained from an explicit solvent simulation.

Interpreting the Coulomb-field approximation for generalized-Born electrostatics using boundary-integral equation theory

The importance of molecular electrostatic interactions in aqueous solution has motivated extensive research into physical models and numerical methods for their estimation. The computational costs associated with simulations that include many explicit water molecules have driven the development of implicit-solvent models, with generalized-Born (GB) models among the most popular of these. In this paper, we analyze a boundary-integral equation interpretation for the Coulomb-field approximation (CFA), which plays a central role in most GB models. This interpretation offers new insights into the nature of the CFA, which traditionally has been assessed using only a single point charge in the solute. The boundary-integral interpretation of the CFA allows the use of multiple point charges, or even continuous charge distributions, leading naturally to methods that eliminate the interpolation inaccuracies associated with the Still equation. This approach, which we call boundary-integral-based electrostatic estimation by the CFA (BIBEE/CFA), is most accurate when the molecular charge distribution generates a smooth normal displacement field at the solute-solvent boundary, and CFA-based GB methods perform similarly. Conversely, both methods are least accurate for charge distributions that give rise to rapidly varying or highly localized normal displacement fields. Supporting this analysis are comparisons of the reaction-potential matrices calculated using GB methods and boundary-element-method (BEM) simulations. An approximation similar to BIBEE/CFA exhibits complementary behavior, with superior accuracy for charge distributions that generate rapidly varying normal fields and poorer accuracy for distributions that produce smooth fields. This approximation, BIBEE by preconditioning (BIBEE/P), essentially generates initial guesses for preconditioned Krylov-subspace iterative BEMs. Thus, iterative refinement of the BIBEE/P results recovers the BEM solution; excellent agreement is obtained in only a few iterations. The boundary-integral-equation framework may also provide a means to derive rigorous results explaining how the empirical correction terms in many modern GB models significantly improve accuracy despite their simple analytical forms.

Electrostatics of proteins in dielectric solvent continua. II. First applications in molecular dynamics simulations

In the preceding paper by Stork and Tavan, [J. Chem. Phys. 126, 165105 (2007)], the authors have reformulated an electrostatic theory which treats proteins surrounded by dielectric solvent continua and approximately solves the associated Poisson equation [B. Egwolf and P. Tavan, J. Chem. Phys. 118, 2039 (2003)]. The resulting solution comprises analytical expressions for the electrostatic reaction field (RF) and potential, which are generated within the protein by the polarization of the surrounding continuum. Here the field and potential are represented in terms of Gaussian RF dipole densities localized at the protein atoms. Quite like in a polarizable force field, also the RF dipole at a given protein atom is induced by the partial charges and RF dipoles at the other atoms. Based on the reformulated theory, the authors have suggested expressions for the RF forces, which obey Newton's third law. Previous continuum approaches, which were also built on solutions of the Poisson equation, used to violate the reactio principle required by this law, and thus were inapplicable to molecular dynamics (MD) simulations. In this paper, the authors suggest a set of techniques by which one can surmount the few remaining hurdles still hampering the application of the theory to MD simulations of soluble proteins and peptides. These techniques comprise the treatment of the RF dipoles within an extended Lagrangian approach and the optimization of the atomic RF polarizabilities. Using the well-studied conformational dynamics of alanine dipeptide as the simplest example, the authors demonstrate the remarkable accuracy and efficiency of the resulting RF-MD approach.

Electrostatics of proteins in dielectric solvent continua. I. Newton’s third law marries qE forces

The authors reformulate and revise an electrostatic theory treating proteins surrounded by dielectric solvent continua [B. Egwolf and P. Tavan, J. Chem. Phys. 118, 2039 (2003)] to make the resulting reaction field (RF) forces compatible with Newton's third law. Such a compatibility is required for their use in molecular dynamics (MD) simulations, in which the proteins are modeled by all-atom molecular mechanics force fields. According to the original theory the RF forces, which are due to the electric field generated by the solvent polarization and act on the partial charges of a protein, i.e., the so-called qE forces, can be quite accurately computed from Gaussian RF dipoles localized at the protein atoms. Using a slightly different approximation scheme also the RF energies of given protein configurations are obtained. However, because the qE forces do not account for the dielectric boundary pressure exerted by the solvent continuum on the protein, they do not obey the principle that actio equals reactio as required by Newton's third law. Therefore, their use in MD simulations is severely hampered. An analysis of the original theory has led the authors now to a reformulation removing the main difficulties. By considering the RF energy, which represents the dominant electrostatic contribution to the free energy of solvation for a given protein configuration, they show that its negative configurational gradient yields mean RF forces obeying the reactio principle. Because the evaluation of these mean forces is computationally much more demanding than that of the qE forces, they derive a suggestion how the qE forces can be modified to obey Newton's third law. Various properties of the thus established theory, particularly issues of accuracy and of computational efficiency, are discussed. A sample application to a MD simulation of a peptide in solution is described in the following paper [M. Stork and P. Tavan, J. Chem. Phys., 126, 165106 (2007).

Continuum description of ionic and dielectric shielding for molecular-dynamics simulations of proteins in solution

We extend our continuum description of solvent dielectrics in molecular-dynamics (MD) simulations, which has provided an efficient and accurate solution of the Poisson equation, to ionic solvents as described by the linearized Poisson-Boltzmann (LPB) equation. We start with the formulation of a general theory for the electrostatics of an arbitrarily shaped molecular system, which consists of partially charged atoms and is embedded in a LPB continuum. This theory represents the reaction field induced by the continuum in terms of charge and dipole densities localized within the molecular system. Because these densities cannot be calculated analytically for systems of arbitrary shape, we introduce an atom-based discretization and a set of carefully designed approximations. This allows us to represent the densities by charges and dipoles located at the atoms. Coupled systems of linear equations determine these multipoles and can be rapidly solved by iteration during a MD simulation. The multipoles yield the reaction field forces and energies. Finally, we scrutinize the quality of our approach by comparisons with an analytical solution restricted to perfectly spherical systems and with results of a finite difference method.