Finite Difference Poisson-Boltzmann Method Research Articles

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.

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.

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.

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.

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.

Electrostatic Energy Calculation on the pH-Induced Conformational Change of Influenza Virus Hemagglutinin

The pH-induced conformational change of influenza virus hemagglutinin (HA) has been investigated by calculating the change of electrostatic energy of the fragment of HA2 upon pH change. The average charge and electrostatic free energy are calculated as a function of pH for the fusion peptide (residues 1–20 of HA2) and the polypeptide of residues 54–77 of HA2 by using the finite difference Poisson-Boltzmann method. It is found that as pH decreases from 8 to 5, the electrostatic free energy of the fusogenic state is lowered by ∼2 kcal/\\ and the fusogenic state is less ionized compared to that of the native state for both polypeptides. For the fusion peptide at the fusogenic state, most of ionizable residues are neutral at acidic pH except Glu-11. For the polypeptide of residues 54–77 at the fusogenic state, most of residues except Glu-74 and His-64 are fully charged between pH 5 and pH 8.

Atomic Charge Parameters for the Finite Difference Poisson-Boltzmann Method Using Electronegativity Neutralization.

An optimization of Rappe and Goddard's charge equilibration (QEq) method of assigning atomic partial charges is described. This optimization is designed for fast and accurate calculation of solvation free energies using the finite difference Poisson-Boltzmann (FDPB) method. The optimization is performed against experimental small molecule solvation free energies using the FDPB method and adjusting Rappe and Goddard's atomic electronegativity values. Using a test set of compounds for which experimental solvation energies are available and a rather small number of parameters, very good agreement was obtained with experiment, with a mean unsigned error of about 0.5 kcal/mol. The QEq atomic partial charge assignment method can reflect the effects of the conformational changes and solvent induction on charge distribution in molecules. In the second section of the paper we examined this feature with a study of the alanine dipeptide conformations in water solvent. The different contributions to the energy surface of the dipeptide were examined and compared with the results from fixed CHARMm charge potential, which is widely used for molecular dynamics studies.

Computational Analysis of the Membrane Association of Group IIA Secreted Phospholipases A2: A Differential Role for Electrostatics

Secreted phospholipases A2 (sPLA2's) are enzymes that hydrolyze glycerophospholipids at the sn-2 position, which leads to the production of lipid mediators of many cellular processes. These interfacial enzymes are regulated by their lipid specificity at two levels: membrane binding and substrate recognition. Different sPLA2's utilize different combinations of electrostatic and hydrophobic interactions to adsorb to membrane surfaces, which results in the wide range of membrane binding behaviors observed. Here, the finite difference Poisson Boltzmann (FDPB) method is used to quantitatively analyze the contribution of electrostatic interactions to the membrane association of two highly basic group II sPLA2's: Agkistrodon piscivorus piscivorus (AppD49) sPLA2 and nonpancreatic human group IIA (hGIIA) sPLA2. The calculations predict how membrane binding is affected by ionic strength, membrane composition, substitutions of residues in the enzymes, and the presence of calcium in the active site. In addition, the results provide molecular models for the membrane-associated forms of the enzymes. Furthermore, these models account for (1) changes in orientation and protonation state of both the native and charge reversal forms of the enzymes at the membrane surface and (2) the effect of protein/vesicle aggregation, as observed for hGIIA sPLA2. Importantly, the modeling quantitatively describes the complex membrane binding behaviors of these interfacial enzymes in terms of simple physical forces and provides structural information that is difficult to obtain experimentally. The computational analysis shows that nonspecific electrostatic interactions not only play a major role in recruiting these enzymes to membrane surfaces but also orient the enzymes for productive catalysis at the membrane interface.

How to improve nature: study of the electrostatic properties of the surface of α-lactalbumin

It was recently shown that alpha-lactalbumin interacts with histones and simple models of histone proteins such as positively charged polyamino acids, suggesting that some fundamental aspects of the protein surface electrostatics may come into play. In the present work, the energies of charge-charge interaction in apo- and Ca(2+)-loaded alpha-lactalbumin were calculated using a Tanford-Kirkwood algorithm with either solvent accessibility correction or using a finite difference Poisson-Boltzmann method. The analysis revealed two major regions of alpha-lactalbumin that possessed highly unfavorable electrostatic potentials: (a) the Ca(2+)-binding loop and its neighboring residues and (b) the N-terminal region of the protein. Several individual mutants were prepared to neutralize specific individual surface acidic amino acids at both the N-terminus and Ca(2+)-binding loop of bovine alpha-lactalbumin. These mutants were characterized by intrinsic fluorescence, differential scanning microcalorimetry and circular dichroism. The structural and thermodynamic data agree in every case with the theoretical predictions, confirming that the N-terminal region is very sensitive to changes in charge. For example, desMet D14N mutation destabilizes protein and decreases its calcium affinity. On the other hand, desMet E1M and desMet D37N substitutions increase the thermal stability and calcium affinity. The Met E1Q is characterized by a marked increase in protein stability, whereas desMet E7Q and desMet E11L display a slight increase in calcium affinity and thermal stability. Examination of the unfavorable energy contributed by Glu1 and the energetically favorable consequences of neutralizing this residue suggests that nature may have made an error with bovine alpha-lactalbumin from the viewpoint of stabilizing structure and conformation.

Computational Protein Design Is a Challenge for Implicit Solvation Models

Increasingly complex schemes for representing solvent effects in an implicit fashion are being used in computational analyses of biological macromolecules. These schemes speed up the calculations by orders of magnitude and are assumed to compromise little on essential features of the solvation phenomenon. In this work we examine this assumption. Five implicit solvation models, a surface area-based empirical model, two models that approximate the generalized Born treatment and a finite difference Poisson-Boltzmann method are challenged in situations differing from those where these models were calibrated. These situations are encountered in automatic protein design procedures, whose job is to select sequences, which stabilize a given protein 3D structure, from a large number of alternatives. To this end we evaluate the energetic cost of burying amino acids in thousands of environments with different solvent exposures belonging, respectively, to decoys built with random sequences and to native protein crystal structures. In addition we perform actual sequence design calculations. Except for the crudest surface area-based procedure, all the tested models tend to favor the burial of polar amino acids in the protein interior over nonpolar ones, a behavior that leads to poor performance in protein design calculations. We show, on the other hand, that three of the examined models are nonetheless capable of discriminating between the native fold and many nonnative alternatives, a test commonly used to validate force fields. It is concluded that protein design is a particularly challenging test for implicit solvation models because it requires accurate estimates of the solvation contribution of individual residues. This contrasts with native recognition, which depends less on solvation and more on other nonbonded contributions.

Implementation and testing of stable, fast implicit solvation in molecular dynamics using the smooth-permittivity finite difference Poisson-Boltzmann method.

A fast stable finite difference Poisson-Boltzmann (FDPB) model for implicit solvation in molecular dynamics simulations was developed using the smooth permittivity FDPB method implemented in the OpenEye ZAP libraries. This was interfaced with two widely used molecular dynamics packages, AMBER and CHARMM. Using the CHARMM-ZAP software combination, the implicit solvent model was tested on eight proteins differing in size, structure, and cofactors: calmodulin, horseradish peroxidase (with and without substrate analogue bound), lipid carrier protein, flavodoxin, ubiquitin, cytochrome c, and a de novo designed 3-helix bundle. The stability and accuracy of the implicit solvent simulations was assessed by examining root-mean-squared deviations from crystal structure. This measure was compared with that of a standard explicit water solvent model. In addition we compared experimental and calculated NMR order parameters to obtain a residue level assessment of the accuracy of MD-ZAP for simulating dynamic quantities. Overall, the agreement of the implicit solvent model with experiment was as good as that of explicit water simulations. The implicit solvent method was up to eight times faster than the explicit water simulations, and approximately four times slower than a vacuum simulation (i.e., with no solvent treatment).

A Molecular Dynamics Study of Reovirus Attachment Protein σ1 Reveals Conformational Changes in σ1 Structure

Molecular dynamics simulations were performed using the recently determined crystal structure of the reovirus attachment protein, σ1. These studies were conducted to improve an understanding of two unique features of σ1 structure: the protonation state of Asp 345, which is buried in the σ1 trimer interface, and the flexibility of the protein at a defined region below the receptor-binding head domain. Three copies of aspartic acids Asp 345 and Asp 346 cluster in a solvent-inaccessible and hydrophobic region at the σ1 trimer interface. These residues are hypothesized to mediate conformational changes in σ1 during viral attachment or cell entry. Our results indicate that protonation of Asp 345 is essential to the integrity of the trimeric structure seen by x-ray crystallography, whereas deprotonation induces structural changes that destabilize the trimer interface. This finding was confirmed by electrostatic calculations using the finite difference Poisson-Boltzmann method. Earlier studies show that σ1 can exist in retracted and extended conformations on the viral surface. Since protonated Asp 345 is necessary to form a stable, extended trimer, our results suggest that protonation of Asp 345 may allow for a structural transition from a partially detrimerized molecule to the fully formed trimer seen in the crystal structure. Additional studies were conducted to quantify the previously observed flexibility of σ1 at a defined region below the receptor-binding head domain. Increased mobility was observed for three polar residues (Ser 291, Thr 292, and Ser 293) located within an insertion between the second and third β-spiral repeats of the crystallized portion of the σ1 tail. These amino acids interact with water molecules of the solvent bulk and are responsible for oscillating movement of the head of ∼50° during 5 ns of simulations. This flexibility may facilitate viral attachment and also function in cell entry and disassembly. These findings provide new insights about the conformational dynamics of σ1 that likely underlie the initiation of the reovirus infectious cycle.

A Computational Model for the Electrostatic Sequestration of PI(4,5)P 2 by Membrane-Adsorbed Basic Peptides

The multivalent acidic phospholipid phosphatidylinositol 4,5-bisphosphate (PI(4,5)P 2) plays a key role in many biological processes. Recent studies show that unstructured clusters of basic residues from a number of peripheral proteins can laterally sequester PI(4,5)P 2 in membranes. Specifically, experiments suggest that the basic effector domain of the myristoylated alanine-rich C kinase substrate (MARCKS), or a peptide corresponding to this domain, MARCKS(151–175), sequesters several PI(4,5)P 2 and that this sequestration is due to nonspecific electrostatic interactions. Here, we use the finite difference Poisson-Boltzmann method to test this hypothesis by calculating the electrostatic free energy of lateral sequestration of PI(4,5)P 2 by membrane-adsorbed basic peptides: Lys-7, Lys-13, and FA-MARCKS(151–175), a peptide based on MARCKS(151–175). In agreement with experiments, we find that the electrostatic free energy becomes more favorable when: 1), Lys-13 and FA-MARCKS(151–175) sequester several PI(4,5)P 2; 2), the linear charge density of the basic peptide increases; 3), the mol percent monovalent acidic lipid in the membrane decreases; and 4), the ionic strength of the solution decreases. In addition, the electrostatic sequestration free energy is in excess of the entropic penalty associated with localizing PI(4,5)P 2. Our calculations, thus, provide a structural and quantitative description of the observed interaction of PI(4,5)P 2 with membrane-adsorbed basic sequences.

Protein molecular dynamics with electrostatic force entirely determined by a single Poisson-Boltzmann calculation.

The electrostatic force including the intramolecular Coulombic interactions and the electrostatic contribution of solvation effect were entirely calculated by using the finite difference Poisson-Boltzmann method (FDPB), which was incorporated into the GROMOS96 force field to complete a new finite difference stochastic dynamics procedure (FDSD). Simulations were performed on an insulin dimer. Different relative dielectric constants were successively assigned to the protein interior; a value of 17 was selected as optimal for our system. The simulation data were analyzed and compared with those obtained from 500-ps molecular dynamics (MD) simulation with explicit water and a 500-ps conventional stochastic dynamics (SD) simulation without the mean solvent force. The results indicate that the FDSD method with GROMOS96 force field is suitable to study the dynamics and structure of proteins in solution if used with the optimal protein dielectric constant.

Electrostatics significantly affect the stability of designed homeodomain variants

The role of electrostatic interactions in determining the stability of designed proteins was studied by constructing and analyzing a set of designed variants of the Drosophila engrailed homeodomain. Computational redesign of 29 surface positions results in a 25-fold mutant with moderate stability, similar to the wild-type protein. Incorporating helix dipole and N-capping considerations into the design algorithm by restricting amino acid composition at the helix termini and N-capping positions yields a ninefold mutant of the initial design (a 23-fold mutant of wild-type) that is over 3 kcal mol −1 more stable than the protein resulting from the unbiased design. Four additional proteins were constructed and analyzed to isolate the effects of helix dipole and N-capping interactions in each helix. Based on the results of urea-denaturation experiments and calculations using the finite difference Poisson-Boltzmann method, both classes of interaction are found to increase the stability of the designed proteins significantly. The simple electrostatic model used in the optimization of rotamers by iterative techniques (ORBIT) force-field, which is similar to the electrostatic models used in other protein design force-fields, is unable to predict the experimentally determined stabilities of the designed variants. The helix dipole and N-capping restrictions provide a simple but effective method to incorporate two types of electrostatic interactions that impact protein stability significantly.

Combining Bioinformatics and Biophysics to Understand Protein-Protein and Protein-Ligand Interactions.

INTRODUCTION. The increasing numbers of proteins whose three-dimensional structures have been determined will have major impact on the ability to exploit genomic data. Sequence alignments will become more meaningful, protein structure prediction will become more accurate, and the prediction of protein function will become increasingly refined and precise. Such developments will require that sequence, structure, and physical chemical information be fully integrated and correlated with biological data in as much detail as possible. We have been developing a series of computational tools with the goal of detecting relationships among amino acid sequence, protein structure and protein function. In this context, recent computational advances in using structure to improve sequence alignments, in homology model building and in the calculation of binding affinities will be summarized as will their combined use, with specific application to understanding the principles of protein-protein and protein-ligand interactions. METHODS. Our basic approach involves calculating protein folding and binding free energies as well as contributions of individual amino acids to these free energies, and correlating these energetic contributions with sequence patterns and with physical and chemical properties of the protein. With regard to binding, we are particularly interested in delineating features that may dictate specificity vs. affinity and in predicting binding specificity from sequence alone. Binding free energies are described in terms of electrostatic and hydrophobic interactions. The former are calculated using finite difference Poisson-Boltzmann methods while the latter are generally calculated from free energy-surface area relationships. RESULTS AND DISCUSSIONS. Our approach to the calculation of binding free energies is validated by its ability to select the correct conformation of protein-protein complexes from a large number of alternative complex geometries. The factors that determine binding free energies will be discussed with particular emphasis on understanding how protein interfaces are designed, in a structural sense, to exploit different combinations of electrostatic and hydrophobic interactions to achieve both affinity and specificity. This information is derived in part by calculating binding free energies for different interfaces and correlating these with structural features on the interface and with sequence patterns. In addition, we have created a database of all protein-protein interfaces of known structure. This allows us to carry out statistical analysis of different interfaces which, when combined with binding free energy calculations, provides a

Salt and pH-dependent properties of native and mutant insulin

A fast and effective model for predicting the salt and pH dependent properties of protein complexes is presented. It is based on the formal charge parameter sets of ionizable groups and applied in conjunction with the finite difference Poisson-Boltzmann (FDPB) method to calculate the electrostatic interactions. All simulations were performed on the native 2Zn insulin and its fast-acting mutants such as B9D (B9Ser→Asp), B9E (B9Aer→Glu), B9EB10D (B9Ser→ Glu, B10His→Asp), and B10D (B10His→Asp). The salt and pH dependent properties of these dimers were analyzed from the aspect of electrostatic interaction, and the theoretical basis of the fast-acting behavior of these mutants was explained. It is found that the results agree well with experimental observations.

Future directions in protein function prediction.

New directions in computational methods for the prediction of protein function are discussed. THEMATICS, a method for the location and characterization of the active sites of enzymes, is featured. THEMATICS, for Theoretical Microscopic Titration Curves, is based on well-established finite-difference Poisson-Boltzmann methods for computing the electric field function of a protein. THEMATICS requires only the structure of the subject protein and thus may be applied to proteins that bear no similarity in structure or sequence to any previously characterized protein. The unique features of catalytic sites in proteins are discussed. Discussion of the chemical basis for the predictive powers of THEMATICS is featured in this paper. Some results are given for three illustrative examples: HIV-1 protease, human apurinic/apyrimidinic endonuclease, and human adenosine kinase.