GBr6NL: a generalized Born method for accurately reproducing solvation energy of the nonlinear Poisson-Boltzmann equation

The nonlinear Poisson-Boltzmann (NLPB) equation can provide accurate modeling of electrostatic effects for nucleic acids and highly charged proteins. Generalized Born methods have been developed to mimic the linearized Poisson-Boltzmann (LPB) equation at substantially reduced cost. The computer time for solving the NLPB equation is approximately fivefold longer than for the LPB equation, thus presenting an even greater obstacle. Here we present the first generalized Born method, GBr(6)NL, for mimicking the NLPB equation. GBr(6)NL is adapted from GBr(6), a generalized Born method recently developed to reproduce the solvation energy of the LPB equation [Tjong and Zhou, J. Phys. Chem. B 111, 3055 (2007)]. Salt effects predicted by GBr(6)NL on 55 proteins overall deviate from NLPB counterparts by 0.5 kcal/mol from ionic strengths from 10 to 1000 mM, which is approximately 10% of the average magnitudes of the salt effects. GBr(6)NL predictions for the salts effects on the electrostatic interaction energies of two protein:RNA complexes are very promising.     

Secondary structure bias in generalized Born solvent models: comparison of conformational ensembles and free energy of solvent polarization from explicit and implicit solvation

The effects of the use of three generalized Born (GB) implicit solvent models on the thermodynamics of a simple polyalanine peptide are studied via comparing several hundred nanoseconds of well-converged replica exchange molecular dynamics (REMD) simulations using explicit TIP3P solvent to REMD simulations with the GB solvent models. It is found that when compared to REMD simulations using TIP3P the GB REMD simulations contain significant differences in secondary structure populations, most notably an overabundance of alpha-helical secondary structure. This discrepancy is explored via comparison of the differences in the electrostatic component of the free energy of solvation (DeltaDeltaG(pol)) between TIP3P (via thermodynamic Integration calculations), the GB models, and an implicit solvent model based on the Poisson equation (PE). The electrostatic components of the solvation free energies are calculated using each solvent model for four representative conformations of Ala10. Since the PE model is found to have the best performance with respect to reproducing TIP3P DeltaDeltaG(pol) values, effective Born radii from the GB models are compared to effective Born radii calculated with PE (so-called perfect radii), and significant and numerous deviations in GB radii from perfect radii are found in all GB models. The effect of these deviations on the solvation free energy is discussed, and it is shown that even when perfect radii are used the agreement of GB with TIP3P DeltaDeltaG(pol) values does not improve. This suggests a limit to the optimization of the effective Born radius calculation and that future efforts to improve the accuracy of GB models must extend beyond such optimizations.     

Generalized born model with a simple smoothing function

Based on recent developments in generalized Born (GB) theory that employ rapid volume integration schemes (M. S. Lee, F. R. Salabury, Jr., and C. L. Brooks III, J Chem Phys 2002, 116, 10606) we have recast the calculation of the self-electrostatic solvation energy to utilize a simple smoothing function at the dielectric boundary. The present GB model is formulated in this manner to provide consistency with the Poisson-Boltzmann (PB) theory previously developed to yield numerically stable electrostatic solvation forces based on finite-difference methods (W. Im, D. Beglov, and B. Roux, Comp Phys Commun 1998, 111, 59). Our comparisons show that the present GB model is indeed an efficient and accurate approach to reproduce corresponding PB solvation energies and forces. With only two adjustable parameters--a(0) to modulate the Coulomb field term, and a(1) to include a correction term beyond Coulomb field--the PB solvation energies are reproduced within 1% error on average for a variety of proteins. Detailed analysis shows that the PB energy can be reproduced within 2% absolute error with a confidence of about 95%. In addition, the solvent-exposed surface area of a biomolecule, as commonly used in calculations of the nonpolar solvation energy, can be calculated accurately and efficiently using the simple smoothing function and the volume integration method. Our implicit solvent GB calculations are about 4.5 times slower than the corresponding vacuum calculations. Using the simple smoothing function makes the present GB model roughly three times faster than GB models, which attempt to mimic the Lee-Richards molecular volume.     

Application of the frozen atom approximation to the GB/SA continuum model for solvation free energy

The generalized Born/surface area (GB/SA) continuum model for solvation free energy is a fast and accurate alternative to using discrete water molecules in molecular simulations of solvated systems. However, computational studies of large solvated molecular systems such as enzyme-ligand complexes can still be computationally expensive even with continuum solvation methods simply because of the large number of atoms in the solute molecules. Because in such systems often only a relatively small portion of the system such as the ligand binding site is under study, it becomes less attractive to calculate energies and derivatives for all atoms in the system. To curtail computation while still maintaining high energetic accuracy, atoms distant from the site of interest are often frozen; that is, their coordinates are made invariant. Such frozen atoms do not require energetic and derivative updates during the course of a simulation. Herein we describe methodology and results for applying the frozen atom approach to both the generalized Born (GB) and the solvent accessible surface area (SASA) parts of the GB/SA continuum model for solvation free energy. For strictly pairwise energetic terms, such as the Coulombic and van-der-Waals energies, contributions from pairs of frozen atoms can be ignored. This leaves energetic differences unaffected for conformations that vary only in the positions of nonfrozen atoms. Due to the nonlocal nature of the GB analytical form, however, excluding such pairs from a GB calculation leads to unacceptable inaccuracies. To apply a frozen-atom scheme to GB calculations, a buffer region within the frozen-atom zone is generated based on a user-definable cutoff distance from the nonfrozen atoms. Certain pairwise interactions between frozen atoms in the buffer region are retained in the GB computation. This allows high accuracy in conformational GB comparisons to be maintained while achieving significant savings in computational time compared to the full (nonfrozen) calculation. A similar approach for using a buffer region of frozen atoms is taken for the SASA calculation. The SASA calculation is local in nature, and thus exact SASA energies are maintained. With a buffer region of 8 A for the frozen-atom cases, excellent agreement in differences in energies for three different conformations of cytochrome P450 with a bound camphor ligand are obtained with respect to the nonfrozen cases. For various minimization protocols, simulations run 2 to 10.5 times faster and memory usage is reduced by a factor of 1.5 to 5. Application of the frozen atom method for GB/SA calculations thus can render computationally tractable biologically and medically important simulations such as those used to study ligand-receptor binding conformations and energies in a solvated environment.     

Second derivatives in generalized Born theory

Generalized Born solvation models offer a popular method of including electrostatic aspects of solvation free energies within an analytical model that depends only upon atomic coordinates, charges, and dielectric radii. Here, we describe how second derivatives with respect to Cartesian coordinates can be computed in an efficient manner that can be distributed over multiple processors. This approach makes possible a variety of new methods of analysis for these implicit solvation models. We illustrate three of these methods here: the use of Newton-Raphson optimization to obtain precise minima in solution; normal mode analysis to compute solvation effects on the mechanical properties of DNA; and the calculation of configurational entropies in the MM/GBSA model. An implementation of these ideas, using the Amber generalized Born model, is available in the nucleic acid builder (NAB) code, and we present examples for proteins with up to 45,000 atoms. The code has been implemented for parallel computers using both the OpenMP and MPI environments, and good parallel scaling is seen with as many as 144 OpenMP processing threads or MPI processing tasks.     

The treatment of solvation by a generalized Born model and a self-consistent charge-density functional theory-based tight-binding method

We present a model to calculate the free energies of solvation of small organic compounds as well as large biomolecules. This model is based on a generalized Born (GB) model and a self-consistent charge-density functional theory-based tight-binding (SCC-DFTB) method with the nonelectrostatic contributions to the free energy of solvation modeled in terms of solvent-accessible surface areas (SA). The parametrization of the SCC-DFTB/GBSA model has been based on 60 neutral and six ionic molecules composed of H, C, N, O, and S, and spanning a wide range of chemical groups. Effective atomic radii as parameters have been obtained through Monte Carlo Simulated Annealing optimization in the parameter space to minimize the differences between the calculated and experimental free energies of solvation. The standard error in the free energies of solvation calculated by the final model is 1.11 kcal mol(-1). We also calculated the free energies of solvation for these molecules using a conductor-like screening model (COSMO) in combination with different levels of theory (AM1, SCC-DFTB, and B3LYP/6-31G*) and compared the results with SCC-DFTB/GBSA. To assess the efficiency of our model for large biomolecules, we calculated the free energy of solvation for a HIV protease-inhibitor complex containing 3,204 atoms using the SCC-DFTB/GBSA and the SCC-DFTB/COSMO models, separately. The computed relative free energies of solvation are comparable, while the SCC-DFTB/GBSA model is three to four times more efficient, in terms of computational cost.     

Effective Born radii in the generalized Born approximation: the importance of being perfect

Generalized Born (GB) models provide, for many applications, an accurate and computationally facile estimate of the electrostatic contribution to aqueous solvation. The GB models involve two main types of approximations relative to the Poisson equation (PE) theory on which they are based. First, the self-energy contributions of individual atoms are estimated and expressed as "effective Born radii." Next, the atom-pair contributions are estimated by an analytical function f(GB) that depends upon the effective Born radii and interatomic distance of the atom pairs. Here, the relative impacts of these approximations are investigated by calculating "perfect" effective Born radii from PE theory, and enquiring as to how well the atom-pairwise energy terms from a GB model using these perfect radii in the standard f(GB) function duplicate the equivalent terms from PE theory. In tests on several biological macromolecules, the use of these perfect radii greatly increases the accuracy of the atom-pair terms; that is, the standard form of f(GB) performs quite well. The remaining small error has a systematic and a random component. The latter cannot be removed without significantly increasing the complexity of the GB model, but an alternative choice of f(GB) can reduce the systematic part. A molecular dynamics simulation using a perfect-radii GB model compares favorably with simulations using conventional GB, even though the radii remain fixed in the former. These results quantify, for the GB field, the importance of getting the effective Born radii right; indeed, with perfect radii, the GB model gives a very good approximation to the underlying PE theory for a variety of biomacromolecular types and conformations.     

FACTS: Fast analytical continuum treatment of solvation

An efficient method for calculating the free energy of solvation of a (macro)molecule embedded in a continuum solvent is presented. It is based on the fully analytical evaluation of the volume and spatial symmetry of the solvent that is displaced from around a solute atom by its neighboring atoms. The two measures of solvent displacement are combined in empirical equations to approximate the atomic (or self) electrostatic solvation energy and the solvent accessible surface area. The former directly yields the effective Born radius, which is used in the generalized Born (GB) formula to calculate the solvent-screened electrostatic interaction energy. A comparison with finite-difference Poisson data shows that atomic solvation energies, pair interaction energies, and their sums are evaluated with a precision comparable to the most accurate GB implementations. Furthermore, solvation energies of a large set of protein conformations have an error of only 1.5%. The solvent accessible surface area is used to approximate the nonpolar contribution to solvation. The empirical approach, called FACTS (Fast Analytical Continuum Treatment of Solvation), is only four times slower than using the vacuum energy in molecular dynamics simulations of proteins. Notably, the folded state of structured peptides and proteins is stable at room temperature in 100-ns molecular dynamics simulations using FACTS and the CHARMM force field.     

New analytic approximation to the standard molecular volume definition and its application to generalized Born calculations

In a recent article (Lee, M. S.; Salsbury, F. R. Jr.; Brooks, C. L., III. J Chem Phys 2002, 116, 10606), we demonstrated that generalized Born (GB) theory provides a good approximation to Poisson electrostatic solvation energy calculations if one uses the same definitions of molecular volume for each. In this work, we present a new and improved analytic method for reproducing the Lee-Richards molecular volume, which is the most common volume definition for Poisson calculations. Overall, 1% errors are achieved for absolute solvation energies of a large set of proteins and relative solvation energies of protein conformations. We also introduce an accurate SASA approximation that uses the same machinery employed by our GB method and requires a small addition of computational cost. The combined methodology is shown to yield an efficient and accurate implicit solvent representation for simulations of biopolymers.     

VBSM: a solvation model based on valence bond theory

A new solvation model, called VBSM, is presented. The model combines valence bond (VB) theory with parameters determined for the SM6 solvation model (Kelly, C. P.; Cramer, C. J.; Truhlar, D. G. J. Chem. Theo. Comp. 2005, 1, 1133-1152). VBSM, like SM6, is based on the generalized Born (GB) approximation for bulk electrostatics and atomic surface tensions to account for cavitation, dispersion, and solvent structure (CDS). The solvation free energy of VBSM includes (i) a self-consistent polarization term obtained by using VB atomic charges in a GB reaction field with a VB self-consistent field procedure that minimizes the total energy of the system with respect to the valence bond orbitals and (ii) a geometry-dependent CDS term to account for deviations from bulk-electrostatic solvation. Test calculations for a few systems show that the liquid-phase partial atomic charges obtained by VBSM are in good agreement with liquid-phase charges obtained by charge model CM4 (Kelly, C. P.; Cramer, C. J.; Truhlar, D. G. J. Chem. Theo. Comp. 2005, 1, 1133-1152). Free energies of solvation are calculated for two prototype test cases, namely, for the degenerate S(N)2 reaction of Cl(-) with CH(3)Cl in water and for a Menshutkin reaction in water. These calculations show that the VBSM method provides a practical alternative to single-configuration self-consistent field theory for solvent effects in molecules and chemical reactions.     

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.     

Generalized Born model with a simple, robust molecular volume correction

Generalized Born (GB) models provide a computationally efficient means of representing the electrostatic effects of solvent and are widely used, especially in molecular dynamics (MD). A class of particularly fast GB models is based on integration over an interior volume approximated as a pairwise union of atom spheres-effectively, the interior is defined by a van der Waals rather than Lee-Richards molecular surface. The approximation is computationally efficient, but if uncorrected, allows for high dielectric (water) regions smaller than a water molecule between atoms, leading to decreased accuracy. Here, an earlier pairwise GB model is extended by a simple analytic correction term that largely alleviates the problem by correctly describing the solvent-excluded volume of each pair of atoms. The correction term introduces a free energy barrier to the separation of non-bonded atoms. This free energy barrier is seen in explicit solvent and Lee-Richards molecular surface implicit solvent calculations, but has been absent from earlier pairwise GB models. When used in MD, the correction term yields protein hydrogen bond length distributions and polypeptide conformational ensembles that are in better agreement with explicit solvent results than earlier pairwise models. The robustness and simplicity of the correction preserves the efficiency of the pairwise GB models while making them a better approximation to reality.     

Performance assessment of semiempirical molecular orbital methods in describing halogen bonding: quantum mechanical and quantum mechanical/molecular mechanical-molecular dynamics study

The performance of semiempirical molecular-orbital methods--MNDO, MNDO-d, AM1, RM1, PM3 and PM6--in describing halogen bonding was evaluated, and the results were compared with molecular mechanical (MM) and quantum mechanical (QM) data. Three types of performance were assessed: (1) geometrical optimizations and binding energy calculations for 27 halogen-containing molecules complexed with various Lewis bases (Two of the tested methods, AM1 and RM1, gave results that agree with the QM data.); (2) charge distribution calculations for halobenzene molecules, determined by calculating the solvation free energies of the molecules relative to benzene in explicit and implicit generalized Born (GB) solvents (None of the methods gave results that agree with the experimental data.); and (3) appropriateness of the semiempirical methods in the hybrid quantum-mechanical/molecular-mechanical (QM/MM) scheme, investigated by studying the molecular inhibition of CK2 protein by eight halobenzimidazole and -benzotriazole derivatives using hybrid QM/MM molecular-dynamics (MD) simulations with the inhibitor described at the QM level by the AM1 method and the rest of the system described at the MM level. The pure MM approach with inclusion of an extra point of positive charge on the halogen atom approach gave better results than the hybrid QM/MM approach involving the AM1 method. Also, in comparison with the pure MM-GBSA (generalized Born surface area) binding energies and experimental data, the calculated QM/MM-GBSA binding energies of the inhibitors were improved by replacing the G(GB,QM/MM) solvation term with the corresponding G(GB,MM) term.     

Implicit solvation based on generalized Born theory in different dielectric environments

In this paper we are investigating the effect of the dielectric environment on atomic Born radii used in generalized Born (GB) methods. Motivated by the Kirkwood expression for the reaction field of a single off-center charge in a spherical cavity, we are proposing extended formalisms for the calculation of Born radii as a function of external and internal dielectric constants. We demonstrate that reaction field energies calculated from environmentally dependent Born radii lead to much improved agreement with Poisson-Boltzmann solutions for low dielectric external environments, such as biological membranes or organic solvent, compared to previous methods where the calculation of Born radii does not depend on the environment. We also examine how this new approach can be applied for the calculation of transfer free energies from vacuum to a given external dielectric for a system with an internal dielectric larger than one. This has not been possible with standard GB theory but is relevant when scoring minimized or average structures with implicit solvent.     

Improved methods for semiempirical solvation models

We present improved algorithms for the SMx (x = 1, 1a, 2, 3) solvation models presented previously [see the overview in C. J. Cramer and D. G. Truhlar, J. Comp.-Aided Mol. Design, 6 , 629 (1992)]. These models estimate the free energy of solvation by augmenting a semiempirical Hartree-Fock calculation on the solute with the generalized Born (GB) model for electric polarization of the solvent and a surface tension term based on solvent-accessible surface area. This article presents three improvements in the algorithms used to carry out such calculations, namely (1) an analytical accessible surface area algorithm, (2) a more efficient radial integration scheme for the dielectric screening computation in the GB model, and (3) a damping algorithm for updating the GB contribution to the Fock update during the iterations to achieve a self-consistent field. Improvements (1) and (2) decrease the computer time, and improvement (3) leads to more stable convergence. Improvement (2) removes a small systematic numerical error that was explicitly absorbed into the parameterization in the SMx models. Therefore, we have adjusted the parameters for one of the previous models to yield essentially identical performance as was obtained originally while simultaneously taking advantage of improvement (2). The resulting model is called SM2.1. The fact that we obtain similar results after removing the systematic quadrature bias attests to the robustness of the original parameterization. © 1995 by John Wiley & Sons, Inc.     

Accurate and efficient generalized born model based on solvent accessibility: derivation and application for LogP octanol/water prediction and flexible peptide docking

A novel method for fast and accurate evaluation of the generalized Born radii in macromolecular solvation electrostatics calculations is proposed, based on the solvent accessibility of the first two solvation layers around an atom. The reverse generalized Born radii calculated by the method have correlation coefficient of 98.7% and RMSD of 0.031 A(-1) with the values obtained using a precise but significantly slower numerical boundary element solution. The method is applied to derive an estimate of the free solvation energy difference between octanol and water and to predict LogP octanol-water. A nine-parameter model is optimized on an 81 compound training set and applied to predict LogP(ow) for an external evaluation set of 19 drug molecules with RMSD of 0.9. The new GB approximation is also tested in Monte Carlo docking simulations of the fully flexible p53 peptide fragment to MDM2. The best energy solution found in the simulations has RMSD of 2.8 A to the X-ray structure.     

New Born radii deriving method for Generalized Born model

Here we report a method to calculate Born radii, an important parameter used in a Generalized Born model. Traditional methods to derive Born radii are mostly based on a complicated formula, while our method is easier and more direct. Atoms are classified according to their atom type, and the Born radii of each type are obtained by fitting to experimental solvation free energy. The SMARTS language is used for the exact definition of atoms types, and Ullmann's subgraph isomorphism algorithm is used to deduce the environment. A generic algorithm is used for the parameter fitting because of its efficiency in searching a huge phase space, and its results are then optimized by using the conjugate gradient method. The final parameter set is fitting from a training set containing 357 molecules and is tested using a test set of 44 small organic molecules, and the average error is 0.58 kcal/mol for 36 neutral molecules and is 1.67 kcal/mol for 8 ions. The model is further tested under organic molecules, biopolymers, and a protein-inhibitor complex and yields reliable results in all these cases. This method can be used to accelerate molecular docking calculations.     

Electrostatics of ligand binding: parametrization of the generalized Born model and comparison with the Poisson-Boltzmann approach

An accurate and fast evaluation of the electrostatics in ligand-protein interactions is crucial for computer-aided drug design. The pairwise generalized Born (GB) model, a fast analytical method originally developed for studying the solvation of organic molecules, has been widely applied to macromolecular systems, including ligand-protein complexes. However, this model involves several empirical scaling parameters, which have been optimized for the solvation of organic molecules, peptides, and nucleic acids but not for energetics of ligand binding. Studies have shown that a good solvation energy does not guarantee a correct model of solvent-mediated interactions. Thus, in this study, we have used the Poisson-Boltzmann (PB) approach as a reference to optimize the GB model for studies of ligand-protein interactions. Specifically, we have employed the pairwise descreening approximation proposed by Hawkins et al.(1) for GB calculations and DelPhi for PB calculations. The AMBER all-atom force field parameters have been used in this work. Seventeen protein-ligand complexes have been used as a training database, and a set of atomic descreening parameters has been selected with which the pairwise GB model and the PB model yield comparable results on atomic Born radii, the electrostatic component of free energies of ligand binding, and desolvation energies of the ligands and proteins. The energetics of the 15 test complexes calculated with the GB model using this set of parameters also agrees well with the energetics calculated with the PB method. This is the first time that the GB model has been parametrized and thoroughly compared with the PB model for the electrostatics of ligand binding.     

Calculation of absolute ligand binding free energy to a ribosome-targeting protein as a function of solvent model

A comparative analysis is provided of the effect of different solvent models on the calculation of a potential of mean force (PMF) for determining the absolute binding affinity of the small molecule inhibitor pteroic acid bound to ricin toxin A-chain (RTA). Solvent models include the distance-dependent dielectric constant, several different generalized Born (GB) approximations, and a hybrid explicit/GB-based implicit solvent model. We found that the simpler approximation of dielectric screening and a GB model, with Born radii fitted to a switching-window dielectric-boundary surface Poisson solvent model, severely overpredicted the binding affinity as compared to the experimental value, estimated to range from -4.4 to -6.0 kcal/mol. In contrast, GB models that are parametrized to fit the Lee-Richards molecular surface performed much better, predicting binding free energy within 1-3 kcal/mol of experimental estimates. However, the predicted free-energy profiles of these GB models displayed alternative binding modes not observed in the crystal structure. Finally, the most rigorous and computationally costly approach in this work, which used a hybrid explicit/implicit solvent model, correctly determined a binding funnel in the PMF near the crystallographic bound state and predicted an absolute binding affinity that was 2 kcal/mol more favorable than the estimated experimental binding affinity.     

Continuum solvation models in the linear interaction energy method

The linear interaction energy (LIE) method in combination with two different continuum solvent models has been applied to calculate protein-ligand binding free energies for a set of inhibitors against the malarial aspartic protease plasmepsin II. Ligand-water interaction energies are calculated from both Poisson-Boltzmann (PB) and Generalized Born (GB) continuum models using snapshots from explicit solvent simulations of the ligand and protein-ligand complex. These are compared to explicit solvent calculations, and we find close agreement between the explicit water and PB solvation models. The GB model overestimates the change in solvation energy, and this is caused by consistent underestimation of the effective Born radii in the protein-ligand complex. The explicit solvent LIE calculations and LIE-PB, with our standard parametrization, reproduce absolute experimental binding free energies with an average unsigned error of 0.5 and 0.7 kcal/mol, respectively. The LIE-GB method, however, requires a constant offset to approach the same level of accuracy.