Finite-difference Poisson-Boltzmann Research Articles

Towards a transferable nonelectrostatic model for continuum solvation: The electrostatic and nonelectrostatic energy correction model.

In this work, we introduce an electrostatic and non-electrostatic (ENE) correction to the solvation energy based on the Solvent-Accessible Surface Area (SASA) of the solute and the solvent static dielectric constant. The proposed correction was developed for neutral solutes in non-aqueous solvents, considering three different implicit solvation models based on a Self-Consistent Reaction Field treatment of solute-solvent mutual polarization using an Apparent Surface Charge formalism, namely the Integral Equation Formalism of the Polarizable Continuum Model using a continuous surface charge scheme (PCM), the Solvation Model based on solute electron density (SMD), and the generalized Finite-Difference Poisson-Boltzmann (FDPB) model. The proposed correction was parametrized on a diverse training set of 4980 solvation data from the Solv@tum database of experimental solvation energies, and validated on the non-aqueous subset of the MNSOL database comprising 2140 solvation energies. The performances of the proposed ENE models with minimal and extended parameters formulations have been analyzed and the latter variant has been further compared to the widely used Cavity, Dispersion, and Solvent structural effects (CDS) non-electrostatic model originally developed for the SMx family of implicit solvation models. Overall, a very good agreement between the computed solvation energies with the ENE correction and the reference experimental data has been found on both the training and test sets for all continuum solvation models considered. Furthermore, results for the ENE correction are on par with the reference CDS non-electrostatic model for both SMD and FDPB electrostatics, but with the advantage of using a lower number of parameters and thus an improved transferability between different electrostatics treatments.

Generalizing Continuum Solvation in Crystal to Nonaqueous Solvents: Implementation, Parametrization, and Application to Molecules and Surfaces.

We present an extension of a generalized finite-difference Poisson-Boltzmann (FDPB) continuum solvation model based on a self-consistent reaction field treatment to nonaqueous solvents. Implementation and reparametrization of the cavitation, dispersion, and structural (CDS) effects nonelectrostatic model are presented in CRYSTAL, with applications to both finite and infinite periodic systems. For neutral finite systems, computed errors with respect to available experimental data on free energies of solvation of 2523 solutes in 91 solvents, as well as 144 transfer energies from water to 14 organic solvents are on par with the reference SM12 solvation model for which the CDS parameters have been developed. Calculations performed on a TiO2 anatase surface and compared to VASPsol data revealed an overall very good agreement of computed solvation energies, surface energies, as well as band structure changes upon solvation in three different solvents, validating the general applicability of the reparametrized FDPB approach to neutral nonperiodic and periodic solutes in aqueous and nonaqueous solvents. For ionic species, while the reparametrized CDS model led to large errors on free energies of solvation of anions, addition of a corrective term based on Abraham's acidity of the solvent significantly improved the accuracy of the proposed continuum solvation model, leading to errors on aqueous pKa of a test set of 83 solutes divided by a factor of 4 compared to the reference solvation model based on density (SMD). Overall, therefore, these encouraging results demonstrate that the generalized FDPB continuum solvation model can be applied to a broad range of solutes in various solvents, ranging from finite neutral or charged solutes to extended periodic surfaces.

Acceleration of Linear Finite-Difference Poisson-Boltzmann Methods on Graphics Processing Units.

Electrostatic interactions play crucial roles in biophysical processes such as protein folding and molecular recognition. Poisson-Boltzmann equation (PBE)-based models have emerged as widely used in modeling these important processes. Though great efforts have been put into developing efficient PBE numerical models, challenges still remain due to the high dimensionality of typical biomolecular systems. In this study, we implemented and analyzed commonly used linear PBE solvers for the ever-improving graphics processing units (GPU) for biomolecular simulations, including both standard and preconditioned conjugate gradient (CG) solvers with several alternative preconditioners. Our implementation utilizes the standard Nvidia CUDA libraries cuSPARSE, cuBLAS, and CUSP. Extensive tests show that good numerical accuracy can be achieved given that the single precision is often used for numerical applications on GPU platforms. The optimal GPU performance was observed with the Jacobi-preconditioned CG solver, with a significant speedup over standard CG solver on CPU in our diversified test cases. Our analysis further shows that different matrix storage formats also considerably affect the efficiency of different linear PBE solvers on GPU, with the diagonal format best suited for our standard finite-difference linear systems. Further efficiency may be possible with matrix-free operations and integrated grid stencil setup specifically tailored for the banded matrices in PBE-specific linear systems.

Computational Tools for Interpreting Ion Channel pH-Dependence.

Activity in many biological systems is mediated by pH, involving proton titratable groups with pKas in the relevant pH range. Experimental analysis of pH-dependence in proteins focusses on particular sidechains, often with mutagenesis of histidine, due to its pKa near to neutral pH. The key question for algorithms that predict pKas is whether they are sufficiently accurate to effectively narrow the search for molecular determinants of pH-dependence. Through analysis of inwardly rectifying potassium (Kir) channels and acid-sensing ion channels (ASICs), mutational effects on pH-dependence are probed, distinguishing between groups described as pH-coupled or pH-sensor. Whereas mutation can lead to a shift in transition pH between open and closed forms for either type of group, only for pH-sensor groups does mutation modulate the amplitude of the transition. It is shown that a hybrid Finite Difference Poisson-Boltzmann (FDPB) – Debye-Hückel continuum electrostatic model can filter mutation candidates, providing enrichment for key pH-coupled and pH-sensor residues in both ASICs and Kir channels, in comparison with application of FDPB alone.

Hemoglobin Bohr Effects: Atomic Origin of the Histidine Residue Contributions

The Bohr effect in hemoglobin, which refers to the dependence of the oxygen affinity on the pH, plays an important role in its cooperativity and physiological function. The dominant contribution to the Bohr effect arises from the difference in the pKa values of His residues of the unliganded (deoxy) and liganded (carbonmonoxy) structures. Using recent high resolution structures, the residue pKa values corresponding to the two structures are calculated. The method is based on determining the electrostatic interactions between residues in the protein, relative to those of the residue in solution, by use of the linearized finite difference Poisson-Boltzmann equation and Monte Carlo sampling of protonation states. Given that good agreement is obtained with the available experimental values for the contribution of His residues in HbA to the Bohr effect, the calculated results are used to determine the atomic origin of the pKa shift between deoxy and carbonmonoxy HbA. The contributions to the pKa shift calculated by means of the linear response approximation show that the salt bridge involving His146 plays an important role in the alkaline Bohr effect, as suggested by Perutz but that other interactions are significant as well. A corresponding analysis is made for the contribution of His143 to the acid Bohr effect for which there is no proposed explanation. The method used is summarized and the program by which it is implemented is described in the Appendix .

Protein Dielectric Constants Determined from NMR Chemical Shift Perturbations

Understanding the connection between protein structure and function requires a quantitative understanding of electrostatic effects. Structure-based electrostatic calculations are essential for this purpose, but their use has been limited by a long-standing discussion on which value to use for the dielectric constants (ε(eff) and ε(p)) required in Coulombic and Poisson-Boltzmann models. The currently used values for ε(eff) and ε(p) are essentially empirical parameters calibrated against thermodynamic properties that are indirect measurements of protein electric fields. We determine optimal values for ε(eff) and ε(p) by measuring protein electric fields in solution using direct detection of NMR chemical shift perturbations (CSPs). We measured CSPs in 14 proteins to get a broad and general characterization of electric fields. Coulomb's law reproduces the measured CSPs optimally with a protein dielectric constant (ε(eff)) from 3 to 13, with an optimal value across all proteins of 6.5. However, when the water-protein interface is treated with finite difference Poisson-Boltzmann calculations, the optimal protein dielectric constant (ε(p)) ranged from 2 to 5 with an optimum of 3. It is striking how similar this value is to the dielectric constant of 2-4 measured for protein powders and how different it is from the ε(p) of 6-20 used in models based on the Poisson-Boltzmann equation when calculating thermodynamic parameters. Because the value of ε(p) = 3 is obtained by analysis of NMR chemical shift perturbations instead of thermodynamic parameters such as pK(a) values, it is likely to describe only the electric field and thus represent a more general, intrinsic, and transferable ε(p) common to most folded proteins.

A Measure of the Broad Substrate Specificity of Enzymes Based on ‘Duplicate’ Catalytic Residues

The ability of an enzyme to select and act upon a specific class of compounds with unerring precision and efficiency is an essential feature of life. Simultaneously, these enzymes often catalyze the reaction of a range of similar substrates of the same class, and also have promiscuous activities on unrelated substrates. Previously, we have established a methodology to quantify promiscuous activities in a wide range of proteins. In the current work, we quantitatively characterize the active site for the ability to catalyze distinct, yet related, substrates (BRASS). A protein with known structure and active site residues provides the framework for computing ‘duplicate’ residues, each of which results in slightly modified replicas of the active site scaffold. Such spatial congruence is supplemented by Finite difference Poisson Boltzmann analysis which filters out electrostatically unfavorable configurations. The congruent configurations are used to compute an index (BrassIndex), which reflects the broad substrate profile of the active site. We identify an acetylhydrolase and a methyltransferase as having the lowest and highest BrassIndex, respectively, from a set of non-homologous proteins extracted from the Catalytic Site Atlas. The acetylhydrolase, a regulatory enzyme, is known to be highly specific for platelet-activating factor. In the methyltransferase (PDB: 1QAM), various combinations of glycine (Gly38/40/42), asparagine (Asn101/11) and glutamic acid (Glu59/36) residues having similar spatial and electrostatic profiles with the specified scaffold (Gly38, Asn101 and Glu59) exemplifies the broad substrate profile such an active site may provide. ‘Duplicate’ residues identified by relaxing the spatial and/or electrostatic constraints can be the target of directed evolution methodologies, like saturation mutagenesis, for modulating the substrate specificity of proteins.

Reducing grid-dependence in finite-difference Poisson-Boltzmann calculations.

Grid dependence in numerical reaction field energies and solvation forces is a well-known limitation in the finite-difference Poisson-Boltzmann methods. In this study we have investigated several numerical strategies to overcome the limitation. Specifically, we have included trimer arc dots during analytical molecular surface generation to improve the convergence of numerical reaction field energies and solvation forces. We have also utilized the level set function to trace the molecular surface implicitly to simplify the numerical mapping of the grid-independent solvent excluded surface. We have further explored to combine the weighted harmonic averaging of boundary dielectrics with a charge-based approach to improve the convergence and stability of numerical reaction field energies and solvation forces. Our test data show that the convergence and stability in both numerical energies and forces can be improved significantly when the combined strategy is applied to either the Poisson equation or the full Poisson-Boltzmann equation.

Similar structures but different mechanisms: Prediction of FABPs–membrane interaction by electrostatic calculation

The role of fatty acid binding proteins as intracellular fatty acid transporters may require their direct interaction with membranes. In this way different mechanisms have been previously characterized through experimental studies suggesting different models for FABPs–membrane association, although the process in which the molecule adsorbs to the membrane remains to be elucidated. To estimate the importance of the electrostatic energy in the FABP–membrane interaction, we computationally modeled the interaction of different FABPs with both anionic and neutral membranes. Free Electrostatic Energy of Binding (dE), was computed using Finite Difference Poisson Boltzmann Equation (FDPB) method as implemented in APBS (Adaptive Poisson Boltzmann Solver). Based on the computational analysis, it is found that recruitment to membranes is facilitated by non-specific electrostatic interactions. Also energetic analysis can quantitatively differentiate among the mechanisms of membrane association proposed and determinate the most energetically favorable configuration for the membrane-associated states of different FABPs. This type of calculations could provide a starting point for further computational or experimental analysis.

On-the-fly Numerical Surface Integration for Finite-Difference Poisson-Boltzmann Methods.

Most implicit solvation models require the definition of a molecular surface as the interface that separates the solute in atomic detail from the solvent approximated as a continuous medium. Commonly used surface definitions include the solvent accessible surface (SAS), the solvent excluded surface (SES), and the van der Waals surface. In this study, we present an efficient numerical algorithm to compute the SES and SAS areas to facilitate the applications of finite-difference Poisson-Boltzmann methods in biomolecular simulations. Different from previous numerical approaches, our algorithm is physics-inspired and intimately coupled to the finite-difference Poisson-Boltzmann methods to fully take advantage of its existing data structures. Our analysis shows that the algorithm can achieve very good agreement with the analytical method in the calculation of the SES and SAS areas. Specifically, in our comprehensive test of 1,555 molecules, the average unsigned relative error is 0.27% in the SES area calculations and 1.05% in the SAS area calculations at the grid spacing of 1/2Å. In addition, a systematic correction analysis can be used to improve the accuracy for the coarse-grid SES area calculations, with the average unsigned relative error in the SES areas reduced to 0.13%. These validation studies indicate that the proposed algorithm can be applied to biomolecules over a broad range of sizes and structures. Finally, the numerical algorithm can also be adapted to evaluate the surface integral of either a vector field or a scalar field defined on the molecular surface for additional solvation energetics and force calculations.

Dielectric boundary force in numerical Poisson–Boltzmann methods: Theory and numerical strategies

Continuum modeling of electrostatic interactions based upon the numerical solutions of the Poisson–Boltzmann equation has been widely adopted in biomolecular applications. To extend their applications to molecular dynamics and energy minimization, robust and efficient methodologies to compute solvation forces must be developed. In this study, we have first reviewed the theory for the computation of dielectric boundary force based on the definition of the Maxwell stress tensor. This is followed by a new formulation of the dielectric boundary force suitable for the finite-difference Poisson–Boltzmann methods. We have validated the new formulation with idealized analytical systems and realistic molecular systems.

PKa predictions with a coupled finite difference Poisson–Boltzmann and Debye–Hückel method

Modeling charge interactions is important for understanding many aspects of biological structure and function, and continuum methods such as Finite Difference Poisson-Boltzmann (FDPB) are commonly employed. Calculations of pH-dependence have identified separate populations; surface groups that can be modeled with a simple Debye-Hückel (DH) model, and buried groups, with stronger resultant interactions that are dependent on detailed conformation. This observation led to the development of a combined FDPB and DH method for pK(a) prediction (termed FD/DH). This study reports application of this method to ionizable groups, including engineered buried charges, in staphylococcal nuclease. The data had been made available to interested research groups before publication of mutant structures and/or pK(a) values. Overall, FD/DH calculations perform as intended with low ΔpK(a) values for surface groups (RMSD between predicted and experimental pK(a) values of 0.74), and much larger ΔpK(a) values for the engineered internal groups, with RMSD = 1.64, where mutant structures were known and RMSD = 1.80, where they were modeled. The weaker resultant interactions of the surface groups are determined mostly by charge-charge interactions. For the buried groups, R(2) for correlation between predicted and measured ΔpK(a) values is 0.74, despite the high RMSDs. Charge-charge interactions are much less important, with the R(2) value for buried group ΔpK(a) values increasing to 0.80 when the term describing charge desolvation alone is used. Engineered charge burial delivers a relatively uniform, nonspecific effect, in terms of pK(a) . How the protein environment adapts in atomic detail to deliver this resultant effect is still an open question.

Exploring a coarse-grained distributive strategy for finite-difference Poisson-Boltzmann calculations.

We have implemented and evaluated a coarse-grained distributive method for finite-difference Poisson–Boltzmann (FDPB) calculations of large biomolecular systems. This method is based on the electrostatic focusing principle of decomposing a large fine-grid FDPB calculation into multiple independent FDPB calculations, each of which focuses on only a small and a specific portion (block) of the large fine grid. We first analyzed the impact of the focusing approximation upon the accuracy of the numerical reaction field energies and found that a reasonable relative accuracy of 10−3 can be achieved when the buffering space is set to be 16 grid points and the block dimension is set to be at least (1/6)3 of the fine-grid dimension, as in the one-block focusing method. The impact upon efficiency of the use of buffering space to maintain enough accuracy was also studied. It was found that an “optimal” multi-block dimension exists for a given computer hardware setup, and this dimension is more or less independent of the solute geometries. A parallel version of thedistributive focusing method was also implemented. Given the proper settings, the distributive method was able to achieve respectable parallel efficiency with tested biomolecular systems on a loosely connected computer cluster.

Roles of Boundary Conditions in DNA Simulations: Analysis of Ion Distributions with the Finite-Difference Poisson-Boltzmann Method

The wide use of lattice-sum strategies in biomolecular simulations has raised many questions on potential artifacts in these strategies. One interesting question is the artifacts in the counterion distributions of highly charged systems. As one would anticipate, Coulombic interactions under the periodic boundary condition may deviate noticeably from those under the free boundary condition in the highly charged systems, significantly influencing their counterion distributions. On the other hand, the electrostatic screening due to water molecules and mobile ions may effectively damp the possible periodic distortions in Coulombic interactions. Therefore, the magnitude of periodicity-induced artifacts in counterion distributions is not straightforward to dissect without detailed analyses. In this study, we have developed a hybrid explicit counterion/implicit salt representation of mobile ions to address this question. We have chosen a well-studied DNA for easy validation of the minimal hybrid ion representation. Our detailed analysis of continuum ion distributions, explicit ion distributions, radial counterion distribution functions, and sequence-dependent counterion distributions, however, indicates that periodicity artifacts are not apparent at the surface of the tested DNA. Nevertheless, influence of boundary conditions does show up starting at the second solvation shell and becomes apparent at the cell boundary.

On removal of charge singularity in Poisson–Boltzmann equation

The Poisson-Boltzmann theory has become widely accepted in modeling electrostatic solvation interactions in biomolecular calculations. However the standard practice of atomic point charges in molecular mechanics force fields introduces singularity into the Poisson-Boltzmann equation. The finite-difference/finite-volume discretization approach to the Poisson-Boltzmann equation alleviates the numerical difficulty associated with the charge singularity but introduces discretization error into the electrostatic potential. Decomposition of the electrostatic potential has been explored to remove the charge singularity explicitly to achieve higher numerical accuracy in the solution of the electrostatic potential. In this study, we propose an efficient method to overcome the charge singularity problem. In our framework, two separate equations for two different potentials in two different regions are solved simultaneously, i.e., the reaction field potential in the solute region and the total potential in the solvent region. The proposed method can be readily implemented with typical finite-difference Poisson-Boltzmann solvers and return the singularity-free reaction field potential with a single run. Test runs on 42 small molecules and 4 large proteins show a very high agreement between the reaction field energies computed by the proposed method and those by the classical finite-difference Poisson-Boltzmann method. It is also interesting to note that the proposed method converges faster than the classical method, though additional time is needed to compute Coulombic potential on the dielectric boundary. The higher precision, accuracy, and efficiency of the proposed method will allow for more robust electrostatic calculations in molecular mechanics simulations of complex biomolecular systems.

Prediction of the most favorable configuration in the ACBP–membrane interaction based on electrostatic calculations

Acyl-CoA binding proteins (ACBPs) are highly conserved 10 kDa cytosolic proteins that bind medium- and long-chain acyl-CoA esters. They act as intracellular carriers of acyl-CoA and play a role in acyl-CoA metabolism, gene regulation, acyl-CoA-mediated cell signaling, transport-mediated lipid synthesis, membrane trafficking and also, ACBPs were indicated as a possible inhibitor of diazepam binding to the GABA-A receptor. To estimate the importance of the non-specific electrostatic energy in the ACBP–membrane interaction, we computationally modeled the interaction of HgACBP with both anionic and neutral membranes. To compute the Free Electrostatic Energy of Binding ( dE), we used the Finite Difference Poisson Boltzmann Equation (FDPB) method as implemented in APBS. In the most energetically favorable orientation, ACBP brings charged residues Lys18 and Lys50 and hydrophobic residues Met46 and Leu47 into membrane surface proximity. This conformation suggests that these four ACBP amino acids are most likely to play a leading role in the ACBP–membrane interaction and ligand intake. Thus, we propose that long range electrostatic forces are the first step in the interaction mechanism between ACBP and membranes.

Theoretical determination of absolute free energy of reduction of plastoquinone-9 in photosystem II and of plastoquinone- n in DMF

We employed static continuum electrostatics and multi-conformation continuum electrostatics (MCCE) methods to determine the reduction potential ( E red 0 ) of PQ-9 in a section of Photosystem II (PSII). Both methods relied on the finite difference Poisson–Boltzmann (FDPB) solution. The static method brings out a E red 0 value (0.01 V) that is close to the experimental one (0.05 V), thereby demonstrating that the surrounding environment critically decides the net free energy change. The E red 0 value obtained from MCCE (0.04 V) is even closer to the observed value, thereby indicating the importance of protein side-chain and proton motions in the electron transfer process. Furthermore, density functional theory-dielectric polarisable continuum model (DFT-DPCM) was employed to calculate the absolute free energy of reduction of plastoquinone- n (PQ- n, where n is the number of isoprenoid units) in N, N dimethyl formamide (DMF) solvent. The DFT-DPCM method produced reduction potential values of −0.59 and −0.65 V for PQ-1 and PQ-9, respectively. These are more or less in agreement with the experimentally reported values of −0.64 and −0.62 V, respectively.

Role of electrostatics on membrane binding, aggregation and destabilization induced by NAD(P)H dehydrogenases. Implication in membrane fusion

Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) is considered a classical glycolytic protein that can promote the fusion of phospholipid vesicles and can also play a vital role on in vivo fusogenic events. However, it is not clear how this redox enzyme, which lack conserved structural or sequence motifs related to membrane fusion, catalyze this process. In order to detect if this ability is present in other NAD(P)H dehydrogenases with available structure, spectroscopic studies were performed to evaluate the capability of alcohol dehydrogenase (ADH), glutamic dehydrogenase (GDH) and sorbitol dehydrogenase (SDH) to bind, aggregate, destabilize and fuse vesicles. Based on finite difference Poisson–Boltzmann calculations (FDPB) the protein–membrane interactions were analyzed. A model for the protein–membrane complex in its minimum free energy of interaction was obtained for each protein and the amino acids involved in the binding processes were suggested. A previously undescribed relationship between membrane destabilization and crevices with high electropositive potential on the protein surface was proposed. The putative implication of the non-specific electrostatics on NAD(P)H dehydrogenases induced membrane fusion is discussed.

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.

Harmonic averaging of smooth permittivity functions in finite-difference Poisson–Boltzmann Electrostatics

Finite-difference Poisson–Boltzmann calculations offer an efficient and accurate means for electrostatic characterization of solvated molecules. However, discretization of charge and permittivity results in sensitive dependence on molecular position and orientation relative to the finite-difference grid. In this article, an improved method for limiting the error associated with discretization of the molecular volume, combining harmonic averaging between grid vertices and Gaussian-based smooth permittivity functions, is presented. While both these methods have the broader result of a smoothly varying permittivity, the Gaussian model represents a fundamental description of the dielectric boundary while harmonic averaging serves to provide information about the permittivity between grid points. Grid positional error is reduced by an order of magnitude in calculations of Born ion solvation energies, small molecule and protein solvation energies, and the solvation energy contribution to a protein-inhibitor complex.