A fast and simple method to calculate protonation states in proteins.

A simple model for electrostatic interactions in proteins, based on a distance and position dependent screening of the electrostatic potential, is presented. It is applied in conjunction with a Monte Carlo algorithm to calculate pK(alpha) values of ionizable groups in proteins. The purpose is to furnish a simple, fast, and sufficiently accurate model to be incorporated into molecular dynamic simulations. This will allow for dynamic protonation calculations and for coupling between changes in structure and protonation state during the simulation. The best method of calculating protonation states available today is based on solving the linearized Poisson-Boltzmann equation on a finite difference grid. However, this model consumes far too much computer time to be a practical alternative. Tests are reported for fixed structures on bacteriorhodopsin, lysozyme, myoglobin, and calbindin. The studies include comparisons with Poisson-Boltzmann calculations with dielectric constants 4 and 20 inside the protein, a model with uniform dielectric constant 80 and distance-dependent dielectric models. The accuracy is comparable to that of Poisson-Boltzmann calculations with dielectric constant 20, and it is considerably better than that with epsilon = 4. The time to calculate the protonation at one pH value is at least 100 times less than that of a Poisson-Boltzmann calculation. Proteins 1999;36:474-483.     

Electrostatic effects on the alpha-helix and beta-strand formation of BPTI(16-36) studied by Monte Carlo simulated annealing.

The electrostatic effects on the secondary structure forming tendencies of a peptide fragment with residues 16-36 of bovine pancreatic trypsin inhibitor, BPTI(16-36), are studied using Monte Carlo simulated annealing simulations. We consider three dielectric functions epsilon(r) of distance r: constant dielectric function (epsilon = 2; strong electrostatic interactions) and sigmoidal functions varying from epsilon(0) = 2 to epsilon(infinity) = 47 (intermediate) and to epsilon(infinity) = 78 (weak). Simulations with epsilon = 2 suggest that this peptide exhibits a significant propensity for beta-strand formations in accordance with a beta-sheet structure of the relevant segment in native BPTI. The tendency for alpha-helix formations becomes almost comparable with that of beta-strands in the simulation with epsilon(infinity) = 47, and there appears no appreciable conformational propensity for this case. Finally, the results with epsilon(infinity) = 78 generate low-energy conformations with conspicuous alpha-helices. These findings suggest the possibility that the change in electrostatic interactions can be the key factor for the conformational transitions of peptides between alpha-helix and beta-sheet that have recently been observed in experiments. These changes in electrostatic interactions can arise from those in various environmental factors such as conformations of the rest of the protein molecule and solvent conditions.     

Comparison of calculation and experiment implicates significant electrostatic contributions to the binding stability of barnase and barstar

The contributions of electrostatic interactions to the binding stability of barnase and barstar were studied by the Poisson-Boltzmann model with three different protocols: a), the dielectric boundary specified as the van der Waals (vdW) surface of the protein along with a protein dielectric constant (epsilon (p)) of 4; b), the dielectric boundary specified as the molecular (i.e., solvent-exclusion (SE)) surface along with epsilon (p) = 4; and c), "SE + epsilon (p) = 20." The "vdW + epsilon (p) = 4" and "SE + epsilon (p) = 20" protocols predicted an overall electrostatic stabilization whereas the "SE + epsilon (p) = 4" protocol predicted an overall electrostatic destabilization. The "vdW + epsilon (p) = 4" protocol was most consistent with experiment. It quantitatively reproduced the observed effects of 17 mutations neutralizing charged residues lining the binding interface and the measured coupling energies of six charge pairs across the interface and reasonably rationalized the experimental ionic strength and pH dependences of the binding constant. In contrast, the "SE + epsilon (p) = 4" protocol predicted significantly larger coupling energies of charge pairs whereas the "SE + epsilon (p) = 20" protocol did not predict any pH dependence. This study calls for further scrutiny of the different Poisson-Boltzmann protocols and demonstrates potential danger in drawing conclusions on electrostatic contributions based on a particular calculation protocol.     

Role of electrostatic interaction in the stability of the hexamer of constitutive phycocyanin from Fremyella diplosiphon

The role of electrostatic interactions between the charges carried by the titratable groups in the different aggregates of constitutive phycocyanin from Fremyella diplosiphon has been studied by using a simple theoretical approach based on the modified Tanford-Kirkwood model. The electrostatic potential has also been calculated by means of a numerical solution of the linearized Poisson-Boltzmann equations using the finite-differences technique. The pH dependence of the electrostatic contribution to free energy suggests an electrostatic stabilization of the alpha- and beta-subunits as well as of the (alpha beta)-monomer over a broad pH interval. The charge distributions in the individual alpha- and beta-subunits produce electrostatic complementarity and promote the assembly of the subunits to the (alpha beta)-monomer, as well as of the monomers to the larger trimeric and hexameric aggregates. Trimer-trimer electrostatic interactions exhibit strong pH dependence, predicting an association/dissociation equilibrium with a midpoint at pH 6. The electrostatic trimer-trimer interactions correspond to the steric fit, suggesting that electrostatic interactions may initially help to orient the trimers during aggregation. The distribution of the electrostatic potential of the monomers and of the higher aggregates suggests that it plays an important role also in phycocyanin-linker protein binding.     

A parameterized, continuum electrostatic model for predicting protein pKa values

Recognizing the limits of trying to achieve chemical accuracy for pK(a) calculations with a purely electrostatic model, we include empirical corrections into the Poisson-Boltzmann solver macroscopic electrostatics with atomic detail (Bashford, Biochemistry 1990;29:10219-10225), to improve the reliability and accuracy of the model. The total number of parameters is kept to a minimum to maximize the robustness of the model for compounds outside of the fitting dataset. The parameters are based on: (a) the electrostatic interaction between functional groups close to the titratable site, (b) the electrostatic work required to desolvate the residue, and (c) the site-to-site interactions. These interactions are straightforward to calculate once the electrostatic field has been solved for each residue using the linearized Poisson-Boltzmann equation and are assumed to be linearly related to the intrinsic pK(a). Two hundred and eighty-six residues from 30 proteins are used to determine the empirical parameters, which result in a root mean square error (RMSE) of 0.70 for the entire set. Eight proteins with 46 experimentally known values were excluded from the parameterization to test the model. This test set had a RMSE of 1.08. We show that the parameterized model improves the results over other models, although like other models the error is strongly correlated with the degree to which a residue is buried. The parameters themselves indicate that local effects are most important for determining the pK(a), whereas site-to-site interactions are found to be less significant.     

Direct proton magnetic resonance determination of the pKa of the active center histidine in thiolsubtilisin

The serine proteases constitute a group of endopeptidases whose members owe their catalytic activity to the presence of a catalytic triad of amino acids consisting of a serine, a histidine and an aspartate. The pK(a) values for this histidine have been determined for several cases in which there is a negative charge installed at the serine to mimic the oxyanionic intermediate and related transition state for the catalytic pathway. Instances from this laboratory include (1) replacement of the serine by a cysteine in subtilisin to create a thiolate; (2) formation of monoisopropylphosphoryl-Ser 195 monoanionic phosphodiesters (in trypsin and chymotrypsin, Ser 221 in subtilisins); and (3) tetrahedral boronates formed with peptide boronic acids. The nuclear magnetic resonance (NMR) signals pertinent to this histidine, or signals indirectly reflecting the state of ionization of this histidine, have been used effectively to monitor changes in the active center ionization state. In every case studied, there is elevation of the pK(a) at the histidine when the negative charge is installed at the serine position. Herein is reported the first NMR measurement of the active center His 63 pK(a) in thiolsubtilisin Carlsberg; it is elevated by 3 units compared with the parent enzyme. Using a numerical solution (finite difference) of the Poisson-Boltzmann equation, a protein dielectric constant of 4 provides a good estimate of the experimentally observed pK(a) elevations. Very significantly, a very low protein dielectric constant (epsilon(p) = 3-5) is required in all of the comparisons, and for all three enzymes used (chymotrypsin, trypsin, and subtilisin). Finally, we discuss why the electrostatic perturbation sensed at His of the active center is more amplified by a negative charge on the Ser side than the same charge on the Asp side. A plausible explanation is that the positive charge on the imidazolium ring of the His is localized, with the N(delta 1) carrying a smaller fraction, the N(epsilon 2) carrying the bulk of the positive charge.     

Brownian dynamics simulations of ion atmospheres around polyalanine and B-DNA: effects of biomolecular dielectric

We have extended an earlier Brownian dynamics simulation algorithm for simulating the structural dynamics of ions around biomolecules to accommodate dielectric inhomogeneity. The electrostatic environment of a biomolecule immersed in water was obtained by numerically solving the Poisson equation with the biomolecule treated as a low dielectric region and the solvent treated as a high dielectric region. Instead of using the mean-field type approximations of ion interactions as in the Poisson-Boltzmann model, the ions were treated explicitly by allowing them to evolve dynamically under the electrostatic field of the biomolecule. This model thus accounts for ion-ion correlations and the finite-size effects of the ions. For a 13-residue alpha-helical polyalanine and a 12-base-pair bp B-form DNA, we found that the choice of the dielectric constant of the biomolecule has much larger effects on the mean ionic structure around the biomolecule than on the fluctuational and dynamical properties of the ions surrounding the biomolecule.     

Electrostatic interactions in the association of proteins: an analysis of the thrombin-hirudin complex.

The role of electrostatic interactions in stabilization of the thrombin-hirudin complex has been investigated by means of two macroscopic approaches: the modified Tanford-Kirkwood model and the finite-difference method for numerical solution of the Poisson-Boltzmann equations. The electrostatic potentials around the thrombin and hirudin molecules were asymmetric and complementary, and it is suggested that these fields influence the initial orientation in the process of the complex formation. The change of the electrostatic binding energy due to mutation of acidic residues in hirudin has been calculated and compared with experimentally determined changes in binding energy. In general, the change in electrostatic binding energy for a particular mutation calculated by the modified Tanford-Kirkwood approach agreed well with the experimentally observed change. The finite-difference approach tended to overestimate changes in binding energy when the mutated residues were involved in short-range electrostatic interactions. Decreases in binding energy caused by mutations of amino acids that do not make any direct ionic interactions (e.g., Glu 61 and Glu 62 of hirudin) can be explained in terms of the interaction of these charges with the positive electrostatic potential of thrombin. Differences between the calculated and observed changes in binding energy are discussed in terms of the crystal structure of the thrombin-hirudin complex.     

Calculations of free-energy contributions to protein-RNA complex stabilization

The problem of calculating binding affinities of protein-RNA complexes is addressed by analyzing a computational strategy of modeling electrostatic free energies based on a nonlinear Poisson-Boltzmann (NLPB) model and linear response approximation (LRA). The underlying idea is to treat binding as a two-step process. Solutions to the NLPB equation calculate free energies arising from electronic polarizability and the LRA is constructed from molecular dynamics simulations to model reorganization free energies due to conformational transitions. By implementing a consistency condition of requiring the NLPB model to reproduce the solute-solvent free-energy transitions determined by the LRA, a "macromolecule dielectric constant" (epsilon(m)) for treating reorganization is obtained. The applicability of this hybrid approach was evaluated by calculating the absolute free energy of binding and free-energy changes for amino acid substitutions in the complex between the U1A spliceosomal protein and its cognate RNA hairpin. Depending on the residue substitution, epsilon(m) varied from 3 to 18, and reflected dipolar reorientation not included in the polarization modeled by epsilon(m) = 2. Although the changes in binding affinities from substitutions modeled strictly at the implicit level by the NLPB equation with epsilon(m) = 4 reproduced the experimental values with good overall agreement, substitutions problematic to this simple treatment showed significant improvement when solved by the NLPB-LRA approach.     

Focusing of electric fields in the active site of Cu-Zn superoxide dismutase: effects of ionic strength and amino-acid modification

In this paper we report the implementation of a finite-difference algorithm which solves the linearized Poisson-Boltzmann equation for molecules of arbitrary shape and charge distribution and which includes the screening effects of electrolytes. The microcoding of the algorithm on an ST-100 array processor allows us to obtain electrostatic potential maps in and around a protein, including the effects of ionic strength, in about 30 minutes. We have applied the algorithm to a dimer of the protein Cu-Zn superoxide dismutase (SOD) and compared our results to those obtained from uniform dielectric models based on coulombic potentials. We find that both the shape of the protein-solvent boundary and the ionic strength of the solvent have a profound effect on the potentials in the solvent. For the case of SOD, the cluster of positive charge at the bottom of the active site channel produces a strongly enhanced positive potential due to the focusing of field lines in the channel-a result that cannot be obtained with any uniform dielectric model. The remainder of the protein is surrounded by a weak negative potential. The electrostatic potential of the enzyme seems designed to provide a large cross-sectional area for productive collisions. Based on the ionic strength dependence of the size of the positive potential region emanating from the active site and the repulsive negative potential barrier surrounding the protein, we are able to suggest an explanation for the ionic strength dependence of the activity of the native and chemically modified forms of the enzyme.     

pH-dependent processes in proteins

Recent improvements in the understanding of electrostatic interactions in proteins serve as a focus for the general topic of pH-dependent processes in proteins. The general importance of pH-dependent processes is first set out in terms of hydrogen ion equilibria, stability, ligand interactions, assembly, dynamics, and events in related molecular systems. The development of various theoretical treatments includes various formalisms in addition to the solvent interface model developed by Shire et al. as an extension of the Tanford-Kirkwood treatment. A number of detailed applications of the model are presented and future potentialities are sketched.     

Protein Electrostatics: Rapid Multigrid-Based Newton Algorithm for Solution of the Full Nonlinear Poisson-Boltzmann Equation

Abstract

A new method for solving the full nonlinear Poisson-Boltzmann equation is outlined. This method is robust and efficient, and uses a combination of the multigrid and inexact Newton algorithms. The novelty of this approach lies in the appropriate combination of the two methods, neither of which by themselves are capable of solving the nonlinear problem accurately. Features of the Poisson-Boltzmann equation are fully exploited by each component of the hybrid algorithm to provide robustness and speed. The advantages inherent in this method increase with the size of the problem. The efficacy of the method is illustrated by calculations of the electrostatic potential around the enzyme Superoxide Dismutase. The CPU time required to solve the full nonlinear equation is less than half that needed for a conjugate gradient solution of the corresponding linearized Poisson-Boltzmann equation. The solutions reveal that the field around the active sites is significantly reduced as compared to that obtained by solving the corresponding linearized Poisson-Boltzmann equation. This new method for the nonlinear Poisson-Boltzmann equation will enable fast and accurate solutions of large protein electrostatics problems.     

Dissecting structural and electrostatic interactions of charged groups in alpha-sarcin. An NMR study of some mutants involving the catalytic residues

The cytotoxic ribonuclease alpha-sarcin is the best characterized member of the ribotoxin family. Ribotoxins share a common structural core, catalytic residues, and active site topology with members of the broader family of nontoxic microbial extracellular RNases. They are, however, much more specific in their biological action. To shed light on the highly specific alpha-sarcin activity, we have evaluated the structural and electrostatic interactions of its charged groups, by combining the structural and pK(a) characterization by NMR of several variants with theoretical calculations based on the Tanford-Kirkwood and Poisson-Boltzmann models. The NMR data reveal that the global conformation of wild-type alpha-sarcin is preserved in the H50Q, E96Q, H137Q, and H50/137Q variants, and that His137 is involved in an H-bond that is crucial in maintaining the active site structure and in reinforcing the stability of the enzyme. The loss of this H-bond in the H137Q and H50/137Q variants modifies the local structure of the active site. The pK(a) values of active site groups H50, E96, and H137 in the four variants have been determined by two-dimensional NMR. The catalytic dyad of E96 and H137 is not sensitive to charge replacements, since their pK(a) values vary less than +/-0.3 pH unit with respect to those of the wild type. On the contrary, the pK(a) of His50 undergoes drastic changes when compared to its value in the intact protein. These amount to an increase of 0.5 pH unit or a decrease of 1.1 pH units depending on whether a positive or negative charge is substituted at the active site. The main determinants of the pK(a) values of most of the charged groups in alpha-sarcin have been established by considering the NMR results in conjunction with those derived from theoretical pK(a) calculations. With regard to the active site residues, the H50 pK(a) is chiefly influenced by electrostatic interactions with E96 and H137, whereas the effect of the low dielectric constant and the interaction with R121 appear to be the main determinants of the altered pK(a) value of E96 and H137. Charge-charge interactions and an increased level of burial perturb the pK(a) values of the active site residues of alpha-sarcin, which can account for its reduced ribonucleolytic activity and its high specificity.     

The effect of protein relaxation on charge-charge interactions and dielectric constants of proteins.

The effect of the reorganization of the protein polar groups on charge-charge interaction and the corresponding effective dielectric constant (epsilon(eff)) is examined by the semimicroscopic version of the Protein Dipole Langevin Dipoles (PDLD/S) method within the framework of the Linear Response Approximation (LRA). This is done by evaluating the interactions between ionized residues in the reaction center of Rhodobacter sphaeroides, while taking into account the protein reorganization energy. It is found that an explicit consideration of the protein relaxation leads to a significant increase in epsilon(eff) and that semimicroscopic models that do not take this relaxation into account force one to use a large value for the so-called "protein dielectric constant," epsilon(p), of the Poisson-Boltzmann model or for the corresponding epsilon(in) in the PDLD/S model. An additional increase in epsilon(eff) is expected from the reorganization of ionized residues and from changes in the degree of water penetration. This finding provides further support for the idea that epsilon(in) (or epsilon(p)) represents contributions that are not considered explicitly. The present study also provides a systematic illustration of the nature of epsilon(eff), supporting our previously reported view that charge-charge interactions correspond to a large value of this "dielectric constant," even in protein interiors. It is also pointed out that epsilon(eff) for the interaction between ionizable groups in proteins is very different from the effective dielectric constant, epsilon'(eff), that determines the free energy of ion pairs in proteins (epsilon'(eff) reflects the effect of preoriented protein dipoles). Finally, the problems associated with the search for a general epsilon(in) are discussed. It is clarified that the epsilon(in) that reproduces the effect of protein relaxation on charge-charge interaction is not equal to the epsilon(in) that reproduces the corresponding effect upon formation of individual charges. This reflects fundamental inconsistencies in attempts to cast microscopic concepts in a macroscopic model. Thus one should either use a large epsilon(in) for charge-charge interactions and a small epsilon(in) for charge-dipole interactions or consider the protein relaxation microscopically.     

Can a simple function for the dielectric response model electrostatic effects in globular proteins?

The relationship between the effective dielectric constant that models the electrostatic effect from a charged side chain in a protein was evaluated both experimentally and theoretically. Experimental values were obtained from the shifts in pKa that resulted from point mutations of side chains in subtilisin. Theoretical values were obtained from an iterative solution to Poisson's equation that considers the dielectric response of the protein and the solvent together with charge positions. There is no simple relationship between the effective dielectric constant and the distance from the charge responsible for the interactions. For some charge positions a linear but not a direct proportional relationship of the effective dielectric with distance of separation was observed. Thus, simple models such as a linear distance-dependence for the dielectric response are not suitable to evaluate electrostatic effects in proteins.     

Modeling salt-mediated electrostatics of macromolecules: the discrete surface charge optimization algorithm and its application to the nucleosome

Much progress has been achieved on quantitative assessment of electrostatic interactions on the all-atom level by molecular mechanics and dynamics, as well as on the macroscopic level by models of continuum solvation. Bridging of the two representations-an area of active research-is necessary for studying integrated functions of large systems of biological importance. Following perspectives of both discrete (N-body) interaction and continuum solvation, we present a new algorithm, DiSCO (Discrete Surface Charge Optimization), for economically describing the electrostatic field predicted by Poisson-Boltzmann theory using a discrete set of Debye-Hückel charges distributed on a virtual surface enclosing the macromolecule. The procedure in DiSCO relies on the linear behavior of the Poisson-Boltzmann equation in the far zone; thus contributions from a number of molecules may be superimposed, and the electrostatic potential, or equivalently the electrostatic field, may be quickly and efficiently approximated by the summation of contributions from the set of charges. The desired accuracy of this approximation is achieved by minimizing the difference between the Poisson-Boltzmann electrostatic field and that produced by the linearized Debye-Hückel approximation using our truncated Newton optimization package. DiSCO is applied here to describe the salt-dependent electrostatic environment of the nucleosome core particle in terms of several hundred surface charges. This representation forms the basis for modeling-by dynamic simulations (or Monte Carlo)-the folding of chromatin. DiSCO can be applied more generally to many macromolecular systems whose size and complexity warrant a model resolution between the all-atom and macroscopic levels.     

Comparison of dielectric response models for simulating electrostatic effects in proteins

Two recent approaches for calculating pK shifts in proteins are compared. The first of these uses Coulomb's law with a distance-dependent dielectric permittivity, epsilon (r), to model the screening effects of the environment, and the second uses a finite difference approach to solve Poisson's equation. It is shown that an explicit form of epsilon (r) which has been fitted to experimentally determined values of the dielectric permittivity in a range from 1 to 21 A can be approximated by a linear form in the functionally significant range of charge separations of approximately 3-10 A, but for distances greater than 10 A the effective permittivity is strongly nonlinear. A statistical analysis of the errors in calculated pK shifts due to electrostatic interactions between charges with separations greater than 10 A shows that there are only marginal differences in reliability between using Coulomb's law with an appropriate form of epsilon (r) or the finite difference approach for solving Poisson's equation. Thus it is concluded that pK shifts can be calculated just as well, and with considerably less effort, using Coulomb's law.     

Changes in stability upon charge reversal and neutralization substitution in staphylococcal nuclease are dominated by favorable electrostatic effects

Single site mutations that reverse or neutralize a surface charge were made at 22 ionizable residues in staphylococcal nuclease. Unfolding free energies were obtained by guanidine hydrochloride denaturation. These data, in conjunction with previously obtained stabilities of the corresponding alanine mutants, unequivocally show that the dominant contribution to stability for virtually all of the wild-type side chains examined is the electrostatic effect associated with each residue's charged group. With only a few exceptions, these charges stabilize the native state, with an average loss of 0.5 kcal/mol of stability upon neutralization of a charge. When the charge is reversed, the average destabilization is doubled. Structure-based calculations of electrostatic free energy with the continuum method based on the finite difference solution to the linearized Poisson-Boltzmann equation reproduce the observed energetics when the polarizability in the protein interior is represented with a dielectric constant of 20. However, in some cases, large differences are found, giving insight into possible areas for improvement of the calculations. In particular, it appears that the assumptions made in the calculations about the absence of electrostatic interactions in the denatured state and the energetic consequences of dynamic fluctuations in the native state will have to be further explored.