Continuum Electrostatics Methods Research Articles

Characterizing the water wire in the Gramicidin channel found by Monte Carlo sampling using continuum electrostatics and in molecular dynamics trajectories with conventional or polarizable force fields

Water molecules play a key role in all biochemical processes. They help define the shape of proteins, and they are reactant or product in many reactions and are released as ligands are bound. They facilitate the transfer of protons through transmembrane proton channel, pump and transporter proteins. Continuum electrostatics (CE) force fields used by program Multiconformation CE (MCCE) capture electrostatic interactions in biomolecules with an implicit solvent, which captures the averaged solvent water equilibrium properties. Hybrid CE methods can use explicit water molecules within the protein surrounded by implicit solvent. These hybrid methods permit the study of explicit hydrogen bond networks within the protein and allow analysis of processes such as proton transfer reactions. Yet hybrid CE methods have not been rigorously tested. Here, we present an explicit treatment of water molecules in the Gramicidin A (gA) channel using MCCE and compare the resulting distributions of water molecules and key hydration features against those obtained with explicit solvent Molecular Dynamics (MD) simulations with the nonpolarizable CHARMM36 and polarizable Drude force fields. CHARMM36 leads to an aligned water wire in the channel characterized by a large absolute net water dipole moment; the MCCE and Drude analysis lead to a small net dipole moment as the water molecules change orientation within the channel. The correct orientation is not as yet known, so these calculations identify an open question.

Characterizing the Water Wire in the Gramicidin Channel Found by Monte Carlo Sampling Using Continuum Electrostatics and in Molecular Dynamics Trajectories with Conventional or Polarizable Force Fields

Water molecules play a key role in all biochemical processes. They help define the shape of proteins, and they are reactant or product in many reactions and are released as ligands are bound. They facilitate the transfer of protons through transmembrane proton channel, pump and transporter proteins. Continuum electrostatics (CE) force fields used by program Multiconformation CE (MCCE) capture electrostatic interactions in biomolecules with an implicit solvent, which captures the averaged solvent water equilibrium properties. Hybrid CE methods can use explicit water molecules within the protein surrounded by implicit solvent. These hybrid methods permit the study of explicit hydrogen bond networks within the protein and allow analysis of processes such as proton transfer reactions. Yet hybrid CE methods have not been rigorously tested. Here, we present an explicit treatment of water molecules in the Gramicidin A (gA) channel using MCCE and compare the resulting distributions of water molecules and key hydration features against those obtained with explicit solvent Molecular Dynamics (MD) simulations with the nonpolarizable CHARMM36 and polarizable Drude force fields. CHARMM36 leads to an aligned water wire in the channel characterized by a large absolute net water dipole moment; the MCCE and Drude analysis lead to a small net dipole moment as the water molecules change orientation within the channel. The correct orientation is not as yet known, so these calculations identify an open question.

Simple models for nonpolar solvation: Parameterization and testing.

Implicit solvent models are important for many biomolecular simulations. The polarity of aqueous solvent is essential and qualitatively captured by continuum electrostatics methods like Generalized Born (GB). However, GB does not account for the solvent-induced interactions between exposed hydrophobic sidechains or solute-solvent dispersion interactions. These "nonpolar" effects are often modeled through surface area (SA) energy terms, which lack realism, create mathematical singularities, and have a many-body character. We have explored an alternate, Lazaridis-Karplus (LK) gaussian energy density for nonpolar effects and a dispersion (DI) energy term proposed earlier, associated with GB electrostatics. We parameterized several combinations of GB, SA, LK, and DI energy terms, to reproduce 62 small molecule solvation free energies, 387 protein stability changes due to point mutations, and the structures of 8 protein loops. With optimized parameters, the models all gave similar results, with GBLK and GBDILK giving no performance loss compared to GBSA, and mean errors of 1.7 kcal/mol for the stability changes and 2 Å deviations for the loop conformations. The optimized GBLK model gave poor results in MD of the Trpcage mini-protein, but parameters optimized specifically for MD performed well for Trpcage and three other small proteins. Overall, the LK and DI nonpolar terms are valid alternatives to SA treatments for a range of applications. © 2017 Wiley Periodicals, Inc.

Benchmarking pKa prediction methods for Lys115 in acetoacetate decarboxylase.

Three different pKa prediction methods were used to calculate the pKa of Lys115 in acetoacetate decarboxylase (AADase): the empirical method PROPKA, the multiconformation continuum electrostatics (MCCE) method, and the molecular dynamics/thermodynamic integration (MD/TI) method with implicit solvent. As expected, accurate pKa prediction of Lys115 depends on the protonation patterns of other ionizable groups, especially the nearby Glu76. However, since the prediction methods do not explicitly sample the protonation patterns of nearby residues, this must be done manually. When Glu76 is deprotonated, all three methods give an incorrect pKa value for Lys115. If protonated, Glu76 is used in an MD/TI calculation, the pKa of Lys115 is predicted to be 5.3, which agrees well with the experimental value of 5.9. This result agrees with previous site-directed mutagenesis studies, where the mutation of Glu76 (negative charge when deprotonated) to Gln (neutral) causes no change in Km, suggesting that Glu76 has no effect on the pKa shift of Lys115. Thus, we postulate that the pKa of Glu76 is also shifted so that Glu76 is protonated (neutral) in AADase. Graphical abstract Simulated abundances of protonated species as pH is varied.

Entropy gain due to water release upon ligand binding

Experimental thermodynamic data of the ligand-receptor association showed that the entropy changes upon binding are positive and large enough to be important driving forces of the binding process for a considerable number of ligand-receptor complexes [1]. The expected source behind such an entropy increase is the release of the water molecules from the binding pocket and from around the ligand to get more freedom in the bulk phase. However, this important source of entropy is usually ignored due to the lake of a method to compute it. Therefore, we developed a method to compute the entropic gain due to water release upon ligand binding. Water molecules close to a charged residue are restricted from free rotation due to the interaction between their dipole moments and the electric field. This restriction result in a decrease in their entropy. The loss of entropy can be gained again when the water molecules displaced from area that under the electric field upon the binding process. To compute the loss of entropy due to the electric field, we developed a method based on the knowledge of the electric field which affects the rotation of the water molecules using the continuum electrostatic methods. This provides a possibility for the calculation of the entropic loss from the partition function of the water molecules under an electric field. For monovalent ions, the computed loss in entropy upon solvation correlates well with experimentally measured one. The loss of entropy of water close to hydrophobic molecules (Hydrophobic entropy) cannot be explained by the effect of the electric field from the molecule because the magnitude of the electric field is not strong enough to restrict the mobility of the water molecules. However, the loss in entropy results from the asymmetric interactions of the water molecules with the hydrophobic molecules and the neighbour water molecules. To compute the hydrophobic entropy we developed an empirical method which is derived from a linear relation between the hydration entropy of small hydrophobic molecules and the hydrophobicity of the water molecules around the surface of the small molecule. A training data set of 142 small molecules shows a significant linear relation between the hydrophobic entropy and the computed hydrophobicity.

Optimizing Electrostatic Field Calculations with the Adaptive Poisson–Boltzmann Solver to Predict Electric Fields at Protein–Protein Interfaces. I. Sampling and Focusing

Continuum electrostatics methods are commonly used to calculate electrostatic potentials in proteins and at protein-protein interfaces to aid many types of biophysical studies. Despite their ubiquity throughout the biophysical literature, these calculations are difficult to test against experimental data to determine their accuracy and validity. To address this, we have calculated the Boltzmann-weighted electrostatic field at the midpoint of a nitrile bond placed at a variety of locations on the surface of the protein RalGDS, both in its monomeric form as well as when docked to four different constructs of the protein Rap, and compared the computation results to vibrational absorption energy measurements of the nitrile oscillator. This was done by generating a statistical ensemble of protein structures using enhanced molecular dynamics sampling with the Amber03 force field, followed by solving the linear Poisson-Boltzmann equation for each structure using the Applied Poisson-Boltzmann Solver (APBS) software package. Using a two-stage focusing strategy, we examined numerous second stage box dimensions, grid point densities, box locations, and compared the numerical result to the result obtained from the sum of the numeric reaction field and the analytic Coulomb field. It was found that the reaction field method yielded higher correlation with experiment for the absolute calculation of fields, while the numeric solutions yielded higher correlation with experiment for the relative field calculations. Finer grid spacing typically improved the calculation, although this effect was less pronounced in the reaction field method. These sorts of calculations were also very sensitive to the box location, particularly for the numeric calculations of absolute fields using a 10(3) Å(3) box.

Continuous development of schemes for parallel computing of the electrostatics in biological systems: implementation in DelPhi.

Due to the enormous importance of electrostatics in molecular biology, calculating the electrostatic potential and corresponding energies has become a standard computational approach for the study of biomolecules and nano-objects immersed in water and salt phase or other media. However, the electrostatics of large macromolecules and macromolecular complexes, including nano-objects, may not be obtainable via explicit methods and even the standard continuum electrostatics methods may not be applicable due to high computational time and memory requirements. Here, we report further development of the parallelization scheme reported in our previous work (Li, et al., J. Comput. Chem. 2012, 33, 1960) to include parallelization of the molecular surface and energy calculations components of the algorithm. The parallelization scheme utilizes different approaches such as space domain parallelization, algorithmic parallelization, multithreading, and task scheduling, depending on the quantity being calculated. This allows for efficient use of the computing resources of the corresponding computer cluster. The parallelization scheme is implemented in the popular software DelPhi and results in speedup of several folds. As a demonstration of the efficiency and capability of this methodology, the electrostatic potential, and electric field distributions are calculated for the bovine mitochondrial supercomplex illustrating their complex topology, which cannot be obtained by modeling the supercomplex components alone.

Protein: Ligand Recognition: Simple Models for Electrostatic Effects

Free energy simulations are a powerful tool to study molecular recognition. The most rigorous variants can provide in depth understanding for a particular system, but are not suited for high throughput application to large libraries of compounds. Related, but less expensive methods are increasingly popular, including continuum electrostatic methods like PBSA ("Poisson-Boltzmann Surface Area") and Linear Response or Linear Interaction Energy methods (LRA, LIE). Here, we review the theoretical background of these methods and provide a unified framework. We focus on the electrostatic contributions to the binding free energy, analyzing nonpolar contributions more briefly. The methods reviewed introduce a multi-step pathway for ligand unbinding, with distinct steps that uncharge the bound ligand, then recharge the unbound ligand. They assume that the system responds to the charging/uncharging in a linear way. With this approximation, the free energy can be described by its one or two first derivatives with respect to a progress variable. The methods can then be classified according to which states of the system are actually simulated and the number of free energy derivatives (one or two) that are employed. The analysis should help clarify the relations between several important free energy methods and the approximations they make. It can suggest new ways to test them, and provide routes for their improvement.

Multi-Level and Interleaved Poisson-Boltzmann (MLIPB) Method for Parallel Computing of the Electrostatics and its Application in Delphi

Electrostatic potential and energies originating from the charge distribution within biomolecules have great impact on intra- and inter-molecular interactions. Due to the enormous importance of electrostatics in molecular biology, calculating the electrostatic potential and corresponding energies has become a standard computational approach to study biomolecules and nano-objects immersed in water and salt phase or other media. However, the electrostatics of large macromolecules and macromolecular complexes may not be obtainable via explicit methods and even the standard continuum electrostatics methods may not be applicable due to high computational and memory requirements. Here we will present a new multi-level (parallelization is achieved at different levels of the algorithm) and interleaved (parallelization is implemented by interleaving the computational tasks) method, the MLIPB method, to parallelize standard methods for computing electrostatics potential and energies in the framework of Poisson-Boltzmann equation (PBE). The MLIPB method is implemented in the popular software DelPhi, resulting in speedup of several folds without compromising accuracy or imposing additional assumptions. As a demonstration of efficiency and capability of this methodology, electrostatic potential and energies are calculated to illustrate the pathway of electron transfer between the component complexes and cytochrome c binding in the bovine mitochondrial supercomplex. This work is supported by NIH grant R01GM093937.

Developing hybrid approaches to predict pKavalues of ionizable groups

Accurate predictions of pKa values of titratable groups require taking into account all relevant processes associated with the ionization/deionization. Frequently, however, the ionization does not involve significant structural changes and the dominating effects are purely electrostatic in origin allowing accurate predictions to be made based on the electrostatic energy difference between ionized and neutral forms alone using a static structure and the subtle structural changes be accounted by using dielectric constant larger than two. On another hand, if the change of the charge state is accompanied by a large structural reorganization of the target protein, then the relevant conformational changes have to be explicitly taken into account in the pKa calculations. Here we report a hybrid approach that first predicts the titratable groups whose ionization is expected to cause large conformational changes, termed "problematic" residues, and then applies a special protocol on them, while the rest of the pK(a)s are predicted with rigid backbone approach as implemented in multi-conformation continuum electrostatics (MCCE) method. The backbone representative conformations for "problematic" groups are generated with either molecular dynamics simulations with charged and uncharged amino acid or with ab-initio local segment modeling. The corresponding ensembles are then used to calculate the titration curves of the "problematic" residues and then the results are averaged to obtain the corresponding pKa.

Cooperativity in a Cluster of Carboxylic Groups in the Active Site of a Protein

The active site of staphylococcal nuclease (SNase) has a cluster of four carboxylic residues, two with near normal pKa values, one (Asp-19) with a low pKa of 2.1, and one (Asp-21) with an unusually elevated pKa of 6.5. Crystal structures provide little insight into why Asp-21 absorbs all the repulsive interactions while Asp-19 is actually stabilized in its charged form. To examine this unusual example of cooperativity the pKa values of Asp and Glu residues were measured using NMR spectroscopy in variants in which ionizable groups in the active site were replaced with neutral ones. The low pKa of Asp-19 reflects hydrogen bonds and a favorable Coulomb interaction whereas the elevated pKa of Asp-21 reflects repulsive Coulomb interactions with the other carboxylic groups in the cluster. Even in the variant where all repulsive interactions within the cluster are eliminated, the pKa of Asp-21 is elevated relative to the pKa of Asp in water due to its ability to act as a hydrogen bond donor in the protonated state. The amplification of small perturbations is a hallmark of cooperative systems and is precisely why structure-based calculations of pKa values in clusters of ionizable groups, such as those present in the active sites of many proteins, are extremely challenging. pKa calculations with a variety of continuum electrostatics methods using static and MD-relaxed structures were performed to illustrate the difficulties inherent to pKa calculations in cooperative systems, where calculations can fail owing to amplification of small and unavoidable inaccuracies in structure and in assumptions built into the calculations.

Molecular determinants of the pKa values of Asp and Glu residues in staphylococcal nuclease

Prior computational studies of the acid-unfolding behavior of staphylococcal nuclease (SNase) suggest that the pK(a) values of its carboxylic groups are difficult to reproduce with electrostatics calculations with continuum methods. To examine the molecular determinants of the pK(a) values of carboxylic groups in SNase, the pK(a) values of all 20 Asp and Glu residues were measured with multidimensional and multinuclear NMR spectroscopy in an acid insensitive variant of SNase. The crystal structure of the protein was obtained to describe the microenvironments of the carboxylic groups. Fourteen Asp and Glu residues titrate with relatively normal pK(a) values that are depressed by less than 1.1 units relative to the normal pK(a) of Asp and Glu in water. Only six residues have pK(a) values shifted by more than 1.5 units. Asp-21 has an unusually high pK(a) of 6.5, which is probably the result of interactions with other carboxylic groups at the active site. The most perturbed pK(a) values appear to be governed by hydrogen bonding and not by Coulomb interactions. The pK(a) values calculated with standard continuum electrostatics methods applied to static structures are more depressed than the measured values because Coulomb effects are exaggerated in the calculations. The problems persist even when the protein is treated with the dielectric constant of water. This can be interpreted to imply that structural relaxation is an important determinant of the pK(a) values; however, no major pH-sensitive conformational reorganization of the backbone was detected using NMR spectroscopy.

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.

Investigating the mechanisms of photosynthetic proteins using continuum electrostatics

Computational methods based on continuum electrostatics are widely used in theoretical biochemistry to analyze the function of proteins. Continuum electrostatic methods in combination with quantum chemical and molecular mechanical methods can help to analyze even very complex biochemical systems. In this article, applications of these methods to proteins involved in photosynthesis are reviewed. After giving a short introduction to the basic concepts of the continuum electrostatic model based on the Poisson-Boltzmann equation, we describe the application of this approach to the docking of electron transfer proteins, to the comparison of isofunctional proteins, to the tuning of absorption spectra, to the analysis of the coupling of electron and proton transfer, to the analysis of the effect of membrane potentials on the energetics of membrane proteins, and to the kinetics of charge transfer reactions. Simulations as those reviewed in this article help to analyze molecular mechanisms on the basis of the structure of the protein, guide new experiments, and provide a better and deeper understanding of protein functions.

Electrostatic Effects in a Network of Polar and Ionizable Groups in Staphylococcal Nuclease

His121 and His124 are embedded in a network of polar and ionizable groups on the surface of staphylococcal nuclease. To examine how membership in a network affects the electrostatic properties of ionizable groups, the tautomeric state and the p K a values of these histidines were measured with NMR spectroscopy in the wild-type nuclease and in 13 variants designed to disrupt the network. In the background protein, His121 and His124 titrate with p K a values of 5.2 and 5.6, respectively. In the variants, where the network was disrupted, the p K a values range from 4.03 to 6.46 for His121, and 5.04 to 5.99 for His124. The largest decrease in a p K a was observed when the favorable Coulomb interaction between His121 and Glu75 was eliminated; the largest increase was observed when Tyr91 or Tyr93 was substituted with Ala or Phe. In all variants, the dominant tautomeric state at neutral pH was the N ε2 state. At one level the network behaves as a rigid unit that does not readily reorganize when disrupted: crystal structures of the E75A or E75Q variants show that even when the pivotal Glu75 is removed, the overall configuration of the network was unaffected. On the other hand, a few key hydrogen bonds appear to govern the conformation of the network, and when these bonds are disrupted the network reorganizes. Coulomb interactions within the network report an effective dielectric constant of 20, whereas a dielectric constant of 80 is more consistent with the magnitude of medium to long-range Coulomb interactions in this protein. The data demonstrate that when structures are treated as static, rigid bodies, structure-based p K a calculations with continuum electrostatics method are not useful to treat ionizable groups in cases where p K a values are governed by short-range polar and Coulomb interactions.

Structure‐Based pKa Calculations Using Continuum Electrostatics Methods

Electrostatic free energy is useful for correlating structure with function in proteins in which ionizable groups play essential functional roles. To this end, the pK(a) values of ionizable groups must be known and their molecular determinants must be understood. Structure-based calculations of electrostatic energies and pK(a) values are necessary for this purpose. This unit describes protocols for pK(a) calculations with continuum electrostatics methods based on the numerical solution of the linearized Poisson-Boltzmann equation by the method of finite differences. Critical discussion of key parameters, approximations, and shortcomings of these methods is included. Two protocols are described for calculations with methods modified empirically to maximize agreement between measured and calculated pK(a) values. Applied judiciously, these methods can contribute useful and novel insight into properties of surface ionizable groups in proteins.

Binding of Aminoglycosidic Antibiotics to the Oligonucleotide A-Site Model and 30S Ribosomal Subunit: Poisson−Boltzmann Model, Thermal Denaturation, and Fluorescence Studies

The binding of paromomycin and similar antibiotics to the oligonucleotide A-site model and the small (30S) ribosomal subunit has been studied using continuum electrostatics methods. Crystallographic information from complexes of paromomycin, tobramycin, and Geneticin bound to an A-site oligonucleotide, and paromomycin and streptomycin complexed to the 30S subunit was used as a foundation to develop structures of similar antibiotics in the same ribosomal binding site. Relative binding free energies were calculated by combining the electrostatic contribution, which was obtained by solving the Poisson-Boltzmann equation, with a surface-area-dependent apolar term and contributions from conformational changes. These computed results showed good correlation with the experimental data resulting from fluorescence binding assays and thermal denaturation studies, demonstrating the ability of the Poisson-Boltzmann model to provide insight into the electrostatic mechanisms for aminoglycoside binding and direction for designing more effective antibiotics.

Electrostatic properties and macroscopic electrodiffusion in OmpF porin and mutants

The bacterial porin OmpF found in the outer membrane of E. coli is a wide channel, characterized by its poor selectivity and almost no ion specificity. It has an asymmetric structure, with relatively large entrances and a narrow constriction. By applying continuum electrostatic methods we determine the ionization states of titratable amino acid residues in the protein and calculate self-consistently the electric potential 3-D distribution within the channel. The average electrostatic properties are then represented by an effective fixed charge distribution along the pore which is the input for a macroscopic electrodiffusion model. The theoretical predictions agree with measurements performed under different salt gradients and pH. The sensitivity of reversal potential and conductance to the direction of the salt gradient and the solution pH is captured by the model. The theory is also able to explain the influence of the lipid membrane charge. The same methodology is satisfactorily applied to some OmpF mutants involving slight structural changes but a large number of net charges. The correlation found between atomic structure and ionic selectivity shows that the transport characteristics of wide channels like OmpF and its mutants are mainly regulated by the collective action of a large number of residues, rather than by the specific interactions of residues at particular locations.

Multiconformation continuum electrostatics analysis of the NhaA Na + /H + antiporter of Escherichia coli with functional implications

Sodium proton antiporters are essential enzymes that catalyze the exchange of sodium ions for protons across biological membranes. Protonations and deprotonations of individual amino acid residues and of clusters formed by these residues play an important role in activating these enzymes and in the mechanism of transport. We have used multiconformation continuum electrostatics method to investigate the protonation states of residues in the sodium proton exchanger NhaA from Escherichia coli, the structure of which has been determined recently by x-ray crystallography. Our calculations identify four clusters of electrostatically tightly interacting residues as well as long-range interactions between residues required for activation. The importance of many of these residues has been demonstrated by the characterization of site-directed mutants. A number of residues with extreme pKa values, including several of the "pH sensor," can only undergo protonation/deprotonation reactions subsequent to conformational changes. The results of the calculations provide valuable information on the activation of the antiporter and the role of individual amino acid residues, and provide a solid framework for further experiments.

Molecular mechanisms of pH‐driven conformational transitions of proteins: Insights from continuum electrostatics calculations of acid unfolding

The acid unfolding of staphylococcal nuclease (SNase) is very cooperative (Whitten and García-Moreno, Biochemistry 2000;39:14292-14304). As many as seven hydrogen ions (H+) are bound preferentially by the acid-unfolded state relative to the native (N) state in the pH range 3.2-3.9. To investigate the mechanism of acid unfolding, structure-based pKa calculations were performed with a variety of continuum electrostatic methods. The calculations reproduced successfully the H+ binding properties of the N state between pH 5 and 9, but they systematically overestimated the number of H+ bound upon acid unfolding. The calculated pKa values of all carboxylic residues in the N state were more depressed than they should be. The discrepancy between the observed and the calculated H+ uptake upon acid unfolding was not improved by using high protein dielectric constants, structures relaxed with molecular dynamics, or other empirical modifications implemented previously by others to maximize agreement between measured and calculated pKa values. This suggests an important role for conformational fluctuations of the backbone as important determinants of pKa values of carboxylic groups. Because no global or subglobal conformational changes have been observed previously for SNase under acidic conditions above the acid-unfolding region, these fluctuations must be local. The acid unfolding of SNase does not seem to involve the disruption of the N state by accruement of intramolecular repulsive interactions, nor the protonation of key ion paired carboxylic residues. It is more consistent with modest contributions from many H+ binding groups, with an important role for local conformational fluctuations in the coupling between H+ binding and the global structural transition.