Variational electrostatic projection (VEP) methods for efficient modeling of the macromolecular electrostatic and solvation environment in activated dynamics simulations

New methods for the calculation of electrostatic interactions between the active dynamical region and surrounding external solvated macromolecular environment in hybrid quantum mechanical/molecular mechanical (QM/MM) simulations are presented. The variational electrostatic projection (VEP) method, and related variational reverse-mapping procedure (VEP-RVM) utilize an expansion in Gaussian surface elements for representation of electrostatic interactions. The use of a discretized surface that surrounds the active dynamical region greatly reduces the number of interactions with the particles of the external environment. The methods are tested on two catalytic RNA systems: the hammerhead and the hairpin ribozymes. It is shown that with surface elements numbering from 302 to 1202 points the direct VEP and VEP-RVM methods are able to obtain relative force errors in the range of 0.5-0.05% and 0.09-0.0001%, respectively, using a 4.0 A projection buffer. These results are encouraging and provide an essential step in the development of new variational macromolecular solvent-boundary methods for QM/MM calculations of enzyme reactions.     

GBr(6): a parameterization-free, accurate, analytical generalized born method

The Poisson-Boltzmann (PB) equation is widely used for modeling electrostatic effects and solvation for macromolecules. The generalized Born (GB) model has been developed to mimic PB results at substantial lower computational cost. Here, we report an analytical GB method that reproduces PB results with high accuracy. The analytical approach builds on previous work of Gallicchio and Levy (J. Comput. Chem. 2004, 25, 479), and incorporates an improvement, proposed by Grycuk (J. Chem. Phys. 2003, 119, 4817), of the Coulomb-field approximation used in most GB methods. Tested against PB results, our GB method has an average unsigned relative error of only 0.6% for a representative set of 55 proteins and of 0.4% and 0.3%, respectively, for folded and unfolded conformations of cytochrome b562 sampled in molecular dynamics simulations. The dependencies of the electrostatic solvation free energy on solute and solvent dielectric constants and on salt concentration are fully accounted for in our method.     

An improved pairwise decomposable finite-difference Poisson-Boltzmann method for computational protein design

Our goal is to develop accurate electrostatic models that can be implemented in current computational protein design protocols. To this end, we improve upon a previously reported pairwise decomposable, finite difference Poisson-Boltzmann (FDPB) model for protein design (Marshall et al., Protein Sci 2005, 14, 1293). The improvement involves placing generic sidechains at positions with unknown amino acid identity and explicitly capturing two-body perturbations to the dielectric environment. We compare the original and improved FDPB methods to standard FDPB calculations in which the dielectric environment is completely determined by protein atoms. The generic sidechain approach yields a two to threefold increase in accuracy per residue or residue pair over the original pairwise FDPB implementation, with no additional computational cost. Distance dependent dielectric and solvent-exclusion models were also compared with standard FDPB energies. The accuracy of the new pairwise FDPB method is shown to be superior to these models, even after reparameterization of the solvent-exclusion model.     

Accurate conformation‐dependent molecular electrostatic potentials for high‐throughput in silico drug discovery

The atom‐centered partial charges‐approximation is commonly used in current molecular modeling tools as a computationally inexpensive alternative to quantum mechanics for modeling electrostatics. Even today, the use of partial charges remains useful despite significant advances in improving the efficiency of ab initio methods. Here, we report on new parameters for the EEM and SFKEEM electronegativity equalization‐based methods for rapidly determining partial charges that will accurately model the electrostatic potential of flexible molecules. The developed parameters cover most pharmaceutically relevant chemistries, and charges obtained using these parameters reproduce the B3LYP/cc‐pVTZ reference electrostatic potential of a set of FDA‐approved drug molecules at best to an average accuracy of 13 ± 4 kJ mol^?1; thus, equipped with these parameters electronegativity equalization‐based methods rival the current best non‐quantum mechanical methods, such as AM1‐BCC, in accuracy, yet incur a lower computational cost. Software implementations of EEM and SFKEEM, including the developed parameters, are included in the conformer‐generation tool BALLOON , available free of charge at http://web.abo.fi/fak/mnf/bkf/research/johnson/software.php . © 2009 Wiley Periodicals, Inc. J Comput Chem, 2010     

Numerical integration techniques for curved-element discretizations of molecule-solvent interfaces

Surface formulations of biophysical modeling problems offer attractive theoretical and computational properties. Numerical simulations based on these formulations usually begin with discretization of the surface under consideration; often, the surface is curved, possessing complicated structure and possibly singularities. Numerical simulations commonly are based on approximate, rather than exact, discretizations of these surfaces. To assess the strength of the dependence of simulation accuracy on the fidelity of surface representation, here methods were developed to model several important surface formulations using exact surface discretizations. Following and refining Zauhar's work [J. Comput.-Aided Mol. Des. 9, 149 (1995)], two classes of curved elements were defined that can exactly discretize the van der Waals, solvent-accessible, and solvent-excluded (molecular) surfaces. Numerical integration techniques are presented that can accurately evaluate nonsingular and singular integrals over these curved surfaces. After validating the exactness of the surface discretizations and demonstrating the correctness of the presented integration methods, a set of calculations are presented that compare the accuracy of approximate, planar-triangle-based discretizations and exact, curved-element-based simulations of surface-generalized-Born (sGB), surface-continuum van der Waals (scvdW), and boundary-element method (BEM) electrostatics problems. Results demonstrate that continuum electrostatic calculations with BEM using curved elements, piecewise-constant basis functions, and centroid collocation are nearly ten times more accurate than planar-triangle BEM for basis sets of comparable size. The sGB and scvdW calculations give exceptional accuracy even for the coarsest obtainable discretized surfaces. The extra accuracy is attributed to the exact representation of the solute-solvent interface; in contrast, commonly used planar-triangle discretizations can only offer improved approximations with increasing discretization and associated increases in computational resources. The results clearly demonstrate that the methods for approximate integration on an exact geometry are far more accurate than exact integration on an approximate geometry. A MATLAB implementation of the presented integration methods and sample data files containing curved-element discretizations of several small molecules are available online as supplemental material.     

Microscopic and semimicroscopic calculations of electrostatic energies in proteins by the POLARIS and ENZYMIX programs

Different microscopic and semimicroscopic approaches for calculations of electrostatic energies in macromolecules are examined. This includes the Protein Dipoles Langevin Dipoles (PDLD) method, the semimicroscopic PDLD (PDLD/S) method, and a free energy perturbation (FEP) method. The incorporation of these approaches in the POLARIS and ENZYMIX modules of the MOLARIS package is described in detail. The PDLD electrostatic calculations are augmented by estimates of the relevant hydrophobic and steric contributions, as well as the effects of the ionic strength and external pH. Determination of the hydrophobic energy involves an approach that considers the modification of the effective surface area of the solute by local field effects. The steric contributions are analyzed in terms of the corresponding reorganization energies. Ionic strength effects are studied by modeling the ionic environment around the given system using a grid of residual charges and evaluating the relevant interaction using Coulomb's law with the dielectric constant of water. The performance of the FEP calculations is significantly enhanced by using special boundary conditions and evaluating the long-range electrostatic contributions using the Local Reaction Field (LRF) model. A diverse set of electrostatic effects are examined, including the solvation energies of charges in proteins and solutions, energetics of ion pairs in proteins and solutions, interaction between surface charges in proteins, and effect of ionic strength on such interactions, as well as electrostatic contributions to binding and catalysis in solvated proteins. Encouraging results are obtained by the microscopic and semimicroscopic approaches and the problems associated with some macroscopic models are illustrated. The PDLD and PDLD/S methods appear to be much faster than the FEP approach and still give reasonable results. In particular, the speed and simplicity of the PDLD/S method make it an effective strategy for calculations of electrostatic free energies in interactive docking studies. Nevertheless, comparing the results of the three approaches can provide a useful estimate of the accuracy of the calculated energies. © 1993 John Wiley & Sons, Inc.     

Efficient calculations of coulombic interactions in biomolecular simulations with periodic boundary conditions

To make improved treatments of electrostatic interactions in biomacromolecular simulations, two possibilities are considered. The first is the famous particle–particle and particle–mesh (PPPM) method developed by Hockney and Eastwood, and the second is a new one developed here in their spirit but by the use of the multipole expansion technique suggested by Ladd. It is then numerically found that the new PPPM method gives more accurate results for a two-particle system at small separation of particles. Preliminary numerical examination of the various computational methods for a single configuration of a model BPTI–water system containing about 24,000 particles indicates that both of the PPPM methods give far more accurate values with reasonable computational cost than do the conventional truncation methods. It is concluded the two PPPM methods are nearly comparable in overall performance for the many-particle systems, although the first method has the drawback that the accuracy in the total electrostatic energy is not high for configurations of charged particles randomly generated. © 1993 John Wiley & Sons, Inc.     

On combining Thole's induced point dipole model with fixed charge distributions in molecular mechanics force fields

The Thole induced point dipole model is combined with three different point charge fitting methods, Merz–Kollman (MK), charges from electrostatic potentials using a grid (CHELPG), and restrained electrostatic potential (RESP), and two multipole algorithms, distributed multipole analysis (DMA) and Gaussian multipole model (GMM), which can be used to describe the electrostatic potential (ESP) around molecules in molecular mechanics force fields. This is done to study how the different methods perform when intramolecular polarizability contributions are self‐consistently removed from the fitting done in the force field parametrization. It is demonstrated that the polarizable versions of the partial charge models provide a good compromise between accuracy and computational efficiency in describing the ESP of small organic molecules undergoing conformational changes. For the point charge models, the inclusion of polarizability reduced the the average root mean square error of ESP over the test set by 4–10%. © 2015 Wiley Periodicals, Inc.     

The ENUF method—Ewald summation based on nonuniform fast Fourier transform: Implementation,parallelization, and application

Computer simulations of model systems are widely used to explore striking phenomena in promising applications spanning from physics, chemistry, biology, to materials science and engineering. The long range electrostatic interactions between charged particles constitute a prominent factor in determining structures and states of model systems. How to efficiently calculate electrostatic interactions in simulation systems subjected to partial or full periodic boundary conditions has been a grand challenging task. In the past decades, a large variety of computational schemes has been proposed, among which the Ewald summation method is the most reliable route to accurately deal with electrostatic interactions between charged particles in simulation systems. In addition, extensive efforts have been done to improve computational efficiencies of the Ewald summation based methods. Representative examples are approaches based on cutoffs, reaction fields, multi-poles, multi-grids, and particle-mesh schemes. We sketched an ENUF method, an abbreviation for the Ewald summation method based on the nonuniform fast Fourier transform technique, and have implemented this method in particle-based simulation packages to calculate electrostatic energies and forces at micro- and mesoscopic levels. Extensive computational studies of conformational properties of polyelectrolytes, dendrimer-membrane complexes, and ionic fluids demonstrated that the ENUF method and its derivatives conserve both energy and momentum to floating point accuracy, and exhibit a computational complexity of with optimal physical parameters. These ENUF based methods are attractive alternatives in molecular simulations where high accuracy and efficiency of simulation methods are needed to accelerate calculations of electrostatic interactions at extended spatiotemporal scales.     

A computationally inexpensive modification of the point dipole electrostatic polarization model for molecular simulations

We present an approximation, which allows reduction of computational resources needed to explicitly incorporate electrostatic polarization into molecular simulations utilizing empirical force fields. The proposed method is employed to compute three-body energies of molecular complexes with dipolar electrostatic probes, gas-phase dimerization energies, and pure liquid properties for five systems that are important in biophysical and organic simulations-water, methanol, methylamine, methanethiol, and acetamide. In all the cases, the three-body energies agreed with high level ab initio data within 0.07 kcal/mol, dimerization energies-within 0.43 kcal/mol (except for the special case of the CH(3)SH), and computed heats of vaporization and densities differed from the experimental results by less than 2%. Moreover, because the presented method allows a significant reduction in computational cost, we were able to carry out the liquid-state calculations with Monte Carlo technique. Comparison with the full-scale point dipole method showed that the computational time was reduced by 3.5 to more than 20 times, depending on the system in hand and on the desired level of the full-scale model accuracy, while the difference in energetic results between the full-scale and the presented approximate model was not great in the most cases. Comparison with the nonpolarizable OPLS-AA force field for all the substances involved and with the polarizable POL3 and q90 models for water and methanol, respectively, demonstrates that the presented technique allows reduction of computational cost with no sacrifice of accuracy. We hope that the proposed method will be of benefit to research employing molecular modeling technique in the biophysical and physical organic chemistry areas.     

Incorporation of local structure into kriging models for the prediction of atomistic properties in the water decamer

Machine learning algorithms have been demonstrated to predict atomistic properties approaching the accuracy of quantum chemical calculations at significantly less computational cost. Difficulties arise, however, when attempting to apply these techniques to large systems, or systems possessing excessive conformational freedom. In this article, the machine learning method kriging is applied to predict both the intra‐atomic and interatomic energies, as well as the electrostatic multipole moments, of the atoms of a water molecule at the center of a 10 water molecule (decamer) cluster. Unlike previous work, where the properties of small water clusters were predicted using a molecular local frame, and where training set inputs (features) were based on atomic index, a variety of feature definitions and coordinate frames are considered here to increase prediction accuracy. It is shown that, for a water molecule at the center of a decamer, no single method of defining features or coordinate schemes is optimal for every property. However, explicitly accounting for the structure of the first solvation shell in the definition of the features of the kriging training set, and centring the coordinate frame on the atom‐of‐interest will, in general, return better predictions than models that apply the standard methods of feature definition, or a molecular coordinate frame. © 2016 The Authors. Journal of Computational Chemistry Published by Wiley Periodicals, Inc.     

MLIMC: Machine Learning-Based Implicit-Solvent Monte Carlo

Monte Carlo (MC) methods are important computational tools for molecular structure optimizations and predictions. When solvent effects are explicitly considered, MC methods become very expensive due to the large degree of freedom associated with the water molecules and mobile ions. Alternatively implicit-solvent MC can largely reduce the computational cost by applying a mean field approximation to solvent effects and meanwhile maintains the atomic detail of the target molecule. The two most popular implicit-solvent models are the Poisson-Boltzmann (PB) model and the Generalized Born (GB) model in a way such that the GB model is an approximation to the PB model but is much faster in simulation time. In this work, we develop a machine learning-based implicit-solvent Monte Carlo (MLIMC) method by combining the advantages of both implicit solvent models in accuracy and efficiency. Specifically, the MLIMC method uses a fast and accurate PB-based machine learning (PBML) scheme to compute the electrostatic solvation free energy at each step. We validate our MLIMC method by using a benzene-water system and a protein-water system. We show that the proposed MLIMC method has great advantages in speed and accuracy for molecular structure optimization and prediction.     

Comparison of the accuracy of approximate methods TrESP and TrCAMM for evaluation of pigment coupling in light-harvesting complexes

Due to the large size of light-harvesting complexes, approximate methods are commonly used for calculating the coupling energies between their constituents. Two approximate methods are studied in this work: TrESP and the recently suggested TrCAMM method, based on the expansion of the one-particle density matrix. The quality of approximation of the electrostatic potential in the framework of these two approaches has been compared for two biological pigments as an example—chlorophyll a and bacteriochlorophyll a. It has been shown that both approaches provide high accuracy of approximation of the matrix elements of the electrostatic potential operator. In addition, it has been demonstrated that symmetrization of the one-particle density matrix in the framework of the TrCAMM method significantly improves the calculation accuracy and makes the computational procedure unambiguous.     

Accuracy assessment of semiempirical molecular electrostatic potential of proteins

 Accurate electrostatic maps of proteins are of great importance in research of protein interaction with ligands, solvent media, drugs, and other biomolecules. The large size of real-life proteins imposes severe limitations on computational methods one can use for obtaining the electrostatic map. Well-known accurate second-order M?ller–Plesset and density functional theory methods are not routinely applicable to systems larger than several hundred atoms. Conventional semiempirical tools, as less resource demanding ones, could be an attractive solution but they do not yield sufficiently accurate calculation results with reference to protein systems, as our analysis demonstrates. The present work performs a thorough analysis of the accuracy issues of the modified neglect of differential overlap type semiempirical Hamiltonians AM1 and PM3 on example of the calculation of the molecular electrostatic potential and the dipole moment of natural amino acids. Real capabilities and limitations of these methods with application to protein modeling are discussed. Received: 26 April 2002 / Accepted: 19 September 2002 / Published online: 14 February 2003     

Electrodiffusion: a continuum modeling framework for biomolecular systems with realistic spatiotemporal resolution

A computational framework is presented for the continuum modeling of cellular biomolecular diffusion influenced by electrostatic driving forces. This framework is developed from a combination of state-of-the-art numerical methods, geometric meshing, and computer visualization tools. In particular, a hybrid of (adaptive) finite element and boundary element methods is adopted to solve the Smoluchowski equation (SE), the Poisson equation (PE), and the Poisson-Nernst-Planck equation (PNPE) in order to describe electrodiffusion processes. The finite element method is used because of its flexibility in modeling irregular geometries and complex boundary conditions. The boundary element method is used due to the convenience of treating the singularities in the source charge distribution and its accurate solution to electrostatic problems on molecular boundaries. Nonsteady-state diffusion can be studied using this framework, with the electric field computed using the densities of charged small molecules and mobile ions in the solvent. A solution for mesh generation for biomolecular systems is supplied, which is an essential component for the finite element and boundary element computations. The uncoupled Smoluchowski equation and Poisson-Boltzmann equation are considered as special cases of the PNPE in the numerical algorithm, and therefore can be solved in this framework as well. Two types of computations are reported in the results: stationary PNPE and time-dependent SE or Nernst-Planck equations solutions. A biological application of the first type is the ionic density distribution around a fragment of DNA determined by the equilibrium PNPE. The stationary PNPE with nonzero flux is also studied for a simple model system, and leads to an observation that the interference on electrostatic field of the substrate charges strongly affects the reaction rate coefficient. The second is a time-dependent diffusion process: the consumption of the neurotransmitter acetylcholine by acetylcholinesterase, determined by the SE and a single uncoupled solution of the Poisson-Boltzmann equation. The electrostatic effects, counterion compensation, spatiotemporal distribution, and diffusion-controlled reaction kinetics are analyzed and different methods are compared.     

Efficient and accurate solvation energy calculation from polarizable continuum models

A new approach is proposed to enhance the efficiency and accuracy for calculation of the long-range electrostatic interaction from implicit solvation models, i.e., the polarizable continuum model (PCM) and its variants, conductorlike PCM/conductorlike screening model and integral equation formalism PCM. In these methods the solvent electrostatics effects are represented by a set of discrete apparent charges distributed on tesserae of the molecular cavity surface embedding the solute. In principle, the accuracy of these methods is improved if the cavity surface is tessellated to finer tesserae; however, the computational time is increased rapidly. We show that such undesired dependency between accuracy and efficiency is a result of the inaccurate treatment of the apparent charge self-contribution to the potential and/or electric field. By taking into account the full effects due to the size and curvature of the segment occupied by each apparent charge, the error in calculated electrostatic solvation free energy is essentially zero for ions (point charge at the center of a sphere) regardless of the degree of tessellation used. For molecules where gradient of apparent charge density is nonzero at the cavity surface, we propose a multiple-sampling technique which significantly lowers the calculated error compared to the original PCM methods, especially when very few numbers of tesserae are used.     

Monte Carlo-based linear Poisson-Boltzmann approach makes accurate salt-dependent solvation free energy predictions possible

The prediction of salt-mediated electrostatic effects with high accuracy is highly desirable since many biological processes where biomolecules such as peptides and proteins are key players can be modulated by adjusting the salt concentration of the cellular milieu. With this goal in mind, we present a novel implicit-solvent based linear Poisson-Boltzmann (PB) solver that provides very accurate nonspecific salt-dependent electrostatic properties of biomolecular systems. To solve the linear PB equation by the Monte Carlo method, we use information from the simulation of random walks in the physical space. Due to inherent properties of the statistical simulation method, we are able to account for subtle geometric features in the biomolecular model, treat continuity and outer boundary conditions and interior point charges exactly, and compute electrostatic properties at different salt concentrations in a single PB calculation. These features of the Monte Carlo-based linear PB formulation make it possible to predict the salt-dependent electrostatic properties of biomolecules with very high accuracy. To illustrate the efficiency of our approach, we compute the salt-dependent electrostatic solvation free energies of arginine-rich RNA-binding peptides and compare these Monte Carlo-based PB predictions with computational results obtained using the more mature deterministic numerical methods.     

Development of spherical gaussian-based point charge models for the evaluation of electrostatic potentials

A set of charge- and dipole moment-conserving point charge models of the electron density of a spherical Gaussian-based wavefunction is described. A particular model may be chosen to satisfy the computational and/or accuracy requirements of a given application. Criteria for selection of models for the evaluation of electrostatic potentials are detailed and tested in an application to the porphyrin molecule.     

Method to assess packing quality of transmembrane alpha-helices in proteins. 1. Parametrization using structural data

Integral membrane proteins (MPs) are pharmaceutical targets of exceptional importance. Modern methods of three-dimensional protein structure determination often fail to supply the fast growing field of structure-based drug design with the requested MPs' structures. That is why computational modeling techniques gain a special importance for these objects. Among the principal difficulties limiting application of these methods is the low quality of the MPs' models built in silico. In this series of two papers we present a computational approach to the assessment of the packing "quality" of transmembrane (TM) alpha-helical domains in proteins. The method is based on the concept of protein environment classes, whereby each amino acid residue is described in terms of its environment polarity and accessibility to the membrane. In the first paper we analyze a nonredundant set of 26 TM alpha-helical domains and compute the residues' propensities to five predefined classes of membrane-protein environments. Here we evaluate the proposed approach only by various test sets, cross-validation protocols and ability of the method to delimit the crystal structure of visual rhodopsin, and a number of its erroneous theoretical models. More advanced validation of the method is given in the second article of this series. We assume that the developed "membrane score" method will be helpful in optimizing computer models of TM domains of MPs, especially G-protein coupled receptors.     

Parametric models to compute tryptophan fluorescence wavelengths from classical protein simulations

Fluorescence spectroscopy is an important method to study protein conformational dynamics and solvation structures. Tryptophan (Trp) residues are the most important and practical intrinsic probes for protein fluorescence due to the variability of their fluorescence wavelengths: Trp residues emit in wavelengths ranging from 308 to 360 nm depending on the local molecular environment. Fluorescence involves electronic transitions, thus its computational modeling is a challenging task. We show that it is possible to predict the wavelength of emission of a Trp residue from classical molecular dynamics simulations by computing the solvent‐accessible surface area or the electrostatic interaction between the indole group and the rest of the system. Linear parametric models are obtained to predict the maximum emission wavelengths with standard errors of the order 5 nm. In a set of 19 proteins with emission wavelengths ranging from 308 to 352 nm, the best model predicts the maximum wavelength of emission with a standard error of 4.89 nm and a quadratic Pearson correlation coefficient of 0.81. These models can be used for the interpretation of fluorescence spectra of proteins with multiple Trp residues, or for which local Trp environmental variability exists and can be probed by classical molecular dynamics simulations. © 2018 Wiley Periodicals, Inc.