Calculating Acid-Dissociation Constants of Proteins Using the Boundary Element Method

Bridging Eulerian and Lagrangian Poisson-Boltzmann solvers by ESES

The Poisson-Boltzmann (PB) model is a widely used electrostatic model for biomolecular solvation analysis. Formulated as an elliptic interface problem, the PB model can be numerically solved on either Eulerian meshes using finite difference/finite element methods or Lagrangian meshes using boundary element methods. Molecular surface generators, which produce the discretized dielectric interfaces between solutes and solvents, are critical factors in determining the accuracy and efficiency of the PB solvers. In this work, we investigate the utility of the Eulerian Solvent Excluded Surface (ESES) software for rendering conjugated Eulerian and Lagrangian surface representations, which enables us to numerically validate and compare the quality of Eulerian PB solvers, such as the MIBPB solver, and the Lagrangian PB solvers, such as the TABI-PB solver. Furthermore, with the ESES software and its associated PB solvers, we are able to numerically validate an interesting and useful but often neglected source-target symmetric property associated with the linearized PB model.

The electrostatic potential inside and around α-lactalbumin: Fluctuations and mean-field models

The electrostatic potential (EP) generated by the protein α-lactoalbumin in the presence of added salt is computed as a thermal average at a given point in space. With this aim, constant pH Monte Carlo (MC) simulations are performed within the primitive model, namely, the solvent is treated as a continuum dielectric. The study of the thermal and spatial fluctuations of the EP reveals that they are remarkably high inside the protein. The calculations indicate that fluctuations inside the protein are mainly due to the asymmetric distribution of the charge groups, while the charge fluctuations of the titratable groups play a minor role. The computed EP matches very well with the one obtained from the Poisson equation for the average charge density in spherical symmetry. The Tanford–Kirkwood multipole expansion reproduces the simulated angular-averaged potential rather accurately. Surprisingly, two of the simplest mean-field models, the linear Poisson–Boltzmann (PB) equation and Donnan potential, provide good estimations of the average EP in the effective protein surface (surface EP). The linear PB equation predicts a linear relationship between charge and surface EP, which is numerically reproduced only if the small ions within the protein are taken into account. On the other hand, the partition coefficients of the small ions inside and outside the protein predicted by Donnan theory reproduce reasonably well the simulation results.

Computing Protein pKas Using the TABI Poisson–Boltzmann Solver

A common approach to computing protein pKas uses a continuum dielectric model in which the protein is a low dielectric medium with embedded atomic point charges, the solvent is a high dielectric medium with a Boltzmann distribution of ionic charges, and the pKa is related to the electrostatic free energy which is obtained by solving the Poisson–Boltzmann equation. Starting from the model pKa for a titrating residue, the method obtains the intrinsic pKa and then computes the protonation probability for a given pH including site–site interactions. This approach assumes that acid dissociation does not affect protein conformation aside from adding or deleting charges at titratable sites. In this work, we demonstrate our treecode-accelerated boundary integral (TABI) solver for the relevant electrostatic calculations. The pKa computing procedure is enclosed in a convenient Python wrapper which is publicly available at the corresponding author’s website. Predicted results are compared with experimental pKas for several proteins. Among ongoing efforts to improve protein pKa calculations, the advantage of TABI is that it reduces the numerical errors in the electrostatic calculations so that attention can be focused on modeling assumptions.

Computational alanine scanning mutagenesis--an improved methodological approach

Alanine scanning mutagenesis of protein-protein interfacial residues can be applied to a wide variety of protein complexes to understand the structural and energetic characteristics of the hot-spots. Binding free energies have been estimated with reasonable accuracy with empirical methods, such as Molecular Mechanics/Poisson-Boltzmann surface area (MM-PBSA), and with more rigorous computational approaches like Free Energy Perturbation (FEP) and Thermodynamic Integration (TI). The main objective of this work is the development of an improved methodological approach, with less computational cost, that predicts accurately differences in binding free energies between the wild-type and alanine mutated complexes (DeltaDeltaG(binding)). The method was applied to three complexes, and a mean unsigned error of 0.80 kcal/mol was obtained in a set of 46 mutations. The computational method presented here achieved an overall success rate of 80% and an 82% success rate in residues for which alanine mutation causes an increase in the binding free energy > 2.0 kcal/mol (warm- and hot-spots). This fully atomistic computational methodological approach consists in a computational Molecular Dynamics simulation protocol performed in a continuum medium using the Generalized Born model. A set of three different internal dielectric constants, to mimic the different degree of relaxation of the interface when different types of amino acids are mutated for alanine, have to be used for the proteins, depending on the type of amino acid that is mutated. This method permits a systematic scanning mutagenesis of protein-protein interfaces and it is capable of anticipating the experimental results of mutagenesis, thus guiding new experimental investigations.

Nonlocal continuum electrostatic theory predicts surprisingly small energetic penalties for charge burial in proteins

We study the energetics of burying charges, ion pairs, and ionizable groups in a simple protein model using nonlocal continuum electrostatics. Our primary finding is that the nonlocal response leads to markedly reduced solvent screening, comparable to the use of application-specific protein dielectric constants. Employing the same parameters as used in other nonlocal studies, we find that for a sphere of radius 13.4 Å containing a single +1e charge, the nonlocal solvation free energy varies less than 18 kcal/mol as the charge moves from the surface to the center, whereas the difference in the local Poisson model is ∼35 kcal/mol. Because an ion pair (salt bridge) generates a comparatively more rapidly varying Coulomb potential, energetics for salt bridges are even more significantly reduced in the nonlocal model. By varying the central parameter in nonlocal theory, which is an effective length scale associated with correlations between solvent molecules, nonlocal-model energetics can be varied from the standard local results to essentially zero; however, the existence of the reduction in charge-burial penalties is quite robust to variations in the protein dielectric constant and the correlation length. Finally, as a simple exploratory test of the implications of nonlocal response, we calculate glutamate pK(a) shifts and find that using standard protein parameters (ε(protein) = 2-4), nonlocal results match local-model predictions with much higher dielectric constants. Nonlocality may, therefore, be one factor in resolving discrepancies between measured protein dielectric constants and the model parameters often used to match titration experiments. Nonlocal models may hold significant promise to deepen our understanding of macromolecular electrostatics without substantially increasing computational complexity.

Trans-channel interactions in batrachotoxin-modified skeletal muscle sodium channels: voltage-dependent block by cytoplasmic amines, and the influence of mu-conotoxin GIIIA derivatives and permeant ions

External mu-conotoxins and internal amine blockers inhibit each other's block of voltage-gated sodium channels. We explore the basis of this interaction by measuring the shifts in voltage-dependence of channel inhibition by internal amines induced by two mu-conotoxin derivatives with different charge distributions and net charges. Charge changes on the toxin were made at residue 13, which is thought to penetrate most deeply into the channel, making it likely to have the strongest individual interaction with an internal charged ligand. When an R13Q or R13E molecule was bound to the channel, the voltage dependence of diethylammonium (DEA)-block shifted toward more depolarized potentials (23 mV for R13Q, and 16 mV for R13E). An electrostatic model of the repulsion between DEA and the toxin simulated these data, with a distance between residue 13 of the mu-conotoxin and the DEA-binding site of approximately 15 A. Surprisingly, for tetrapropylammonium, the shifts were only 9 mV for R13Q, and 7 mV for R13E. The smaller shifts associated with R13E, the toxin with a smaller net charge, are generally consistent with an electrostatic interaction. However, the smaller shifts observed for tetrapropylammonium than for DEA suggest that other factors must be involved. Two observations indicate that the coupling of permeant ion occupancy of the channel to blocker binding may contribute to the overall amine-toxin interaction: 1), R13Q binding decreases the apparent affinity of sodium for the conducting pore by approximately 4-fold; and 2), increasing external [Na(+)] decreases block by DEA at constant voltage. Thus, even though a number of studies suggest that sodium channels are occupied by no more than one ion most of the time, measurable coupling occurs between permeant ions and toxin or amine blockers. Such interactions likely determine, in part, the strength of trans-channel, amine-conotoxin interactions.

Intermolecular complexation and phase separation in aqueous solutions of oppositely charged biopolymers

Turbidity measurements performed at 450nm were used to follow the process of complex formation, and phase separation in gelatin-agar aqueous solutions. Acid (Type-A) and alkali (Type-B) processed gelatin (polyampholyte) and agar (anionic polyelectrolyte) solutions, both having concentration of 0.1% (w/v) were mixed in various proportions, and the mixture was titrated (with 0.01 M HCl or NaOH) to initiate associative complexation that led to coacervation. The titration profiles clearly established observable transitions in terms of the solution pH corresponding to the first occurrence of turbidity (pH(C), formation of soluble complexes), and a point of turbidity maximum (pH(phi), formation of insoluble complexes). Decreasing the pH beyond pH(phi) drove the system towards precipitation. The values of pH(C) and pH(phi) characterized the initiation of the formation of intermolecular charge neutralized soluble aggregates, and the subsequent formation of microscopic coacervate droplets. These aggregates were characterized by dynamic light scattering. It was found that Type-A and -B gelatin samples formed soluble intermolecular complexes (and coacervates) with agar molecules through electrostatic and patch-binding interactions, respectively.

Assignment of polar states for protein amino acid residues using an interaction cluster decomposition algorithm and its application to high resolution protein structure modeling

We have developed a new method (Independent Cluster Decomposition Algorithm, ICDA) for creating all-atom models of proteins given the heavy-atom coordinates, provided by X-ray crystallography, and the pH. In our method the ionization states of titratable residues, the crystallographic mis-assignment of amide orientations in Asn/Gln, and the orientations of OH/SH groups are addressed under the unified framework of polar states assignment. To address the large number of combinatorial possibilities for the polar hydrogen states of the protein, we have devised a novel algorithm to decompose the system into independent interacting clusters, based on the observation of the crucial interdependence between the short range hydrogen bonding network and polar residue states, thus significantly reducing the computational complexity of the problem and making our algorithm tractable using relatively modest computational resources. We utilize an all atom protein force field (OPLS) and a Generalized Born continuum solvation model, in contrast to the various empirical force fields adopted in most previous studies. We have compared our prediction results with a few well-documented methods in the literature (WHATIF, REDUCE). In addition, as a preliminary attempt to couple our polar state assignment method with real structure predictions, we further validate our method using single side chain prediction, which has been demonstrated to be an effective way of validating structure prediction methods without incurring sampling problems. Comparisons of single side chain prediction results after the application of our polar state prediction method with previous results with default polar state assignments indicate a significant improvement in the single side chain predictions for polar residues.

Rotationally inelastic scattering of OH (Π3∕22, v=0, J=3∕2, f) by HBr (Σ1, v=0, J<4)

Relative state-to-state cross sections of OH molecules in the (2)Pi(32), v=0, J=32, M(J)=32, f state have been determined for transitions up to (2)Pi(32), v=0, J=112, f and (2)Pi(12), v=0, J=72, e states by collisions with HBr molecules ((1)Sigma, v=0, J<4) at 750 cm(-1) collision energy. In order to investigate features of the anisotropy of the OH-HBr potential energy surface, the steric asymmetries, which account for the effect of the OH orientation with respect to the collision partner, have been measured. A comparison with other systems previously studied shows strong similarities with the OH-HCl system.

Theoretical investigations of prostatic acid phosphatase

The phosphotyrosyl protein phosphatase activity of prostatic acid phosphatase (PAP) has been well established. It has also been suggested that PAP partly regulates the activity of growth factor receptors by dephosphorylating the autophosphorylysable tyrosines in them. We studied the binding of the peptides from epidermal growth factor receptor (EGFR) and its homolog (ErbB-2), corresponding to their autophosphorylation sites, to PAP using theoretical modeling and molecular dynamics (MD) simulation methods. Nine different peptides, each with a phosphotyrosine residue, were docked on human PAP. The binding energies of these peptide-PAP complexes were calculated theoretically and compared to experimentally obtained affinities. The peptide Ace--DNLpYYWD--NH2 from ErbB-2(1197-1203) showed the most favorable free energy of binding when estimated theoretically. The results demonstrate that the presence of another tyrosine residue proximate to C-terminal of autophosphorylysable Tyr enhances the binding affinity considerably. The presence of a bulky group instead prevents the binding, as is observed in case of peptide Ace--NLYpYWDQ--NH2 which failed to bind, both in theoretical calculations and experiments. Thus we demonstarted that PAP could potentially bind to EGFR and Erbb-2 and dephosphorylate them. Thus it could be involved in the regulation of the function of such receptors. In addition, complexes of a peptide from AngiotensinII and phosphotyrosine(pY) with human PAP were also modeled. The effects of different protonation states of the titratable active site residues on ligand (pY) binding have also been investigated. For a favorable binding His12 and Asp258 should be neutral, His257 should be positively charged and the phosphate group of the ligand should be in PO(4) (3-) state. Furthermore, the analysis of protein motion as observed during simulations suggests the loop-loop contact in the PAP dimer to be of importance in cooperativity.

26 kDa Endochitinase from Barley Seeds: An Interaction of the Ionizable Side Chains Essential for Catalysis

To explore the structure essential for the catalysis in 26 kDa endochitinase from barley seeds, we calculated theoretical pKa values of the ionizable groups based on the crystal structure, and then the roles of ionizable side chains located near the catalytic residue were examined by site-directed mutagenesis. The pKa value calculated for Arg215, which is located at the bottom of the catalytic cleft, is abnormally high (>20.0), indicating that the guanidyl group may interact strongly with nearby charges. No enzymatic activity was found in the Arg215-mutated chitinase (R215A) produced by the Escherichia coli expression system. The transition temperature of thermal unfolding (T(m)) of R215A was lower than that of the wild type protein by about 6.2 degrees C. In the crystal structure, the Arg215 side chain is in close proximity to the Glu203 side chain, whose theoretical pKa value was found to be abnormally low (-2.4), suggesting that these side chains may interact with each other. Mutation of Glu203 to alanine (E203A) completely eliminated the enzymatic activity and impaired the thermal stability (deltaT(m) = 6.4 degrees C) of the enzyme. Substrate binding ability was also affected by the Glu203 mutation. These data clearly demonstrate that the Arg215 side chain interacts with the Glu203 side chain to stabilize the conformation of the catalytic cleft. A similar interaction network was previously found in chitosanase from Streptomyces sp. N174 [Fukamizo et al. (2000) J. Biol. Chem. 275, 25633-25640]; hence, this type of interaction seems to be at least partly conserved in the catalytic cleft of other glycosyl hydrolases.

Monte Carlo simulation of an inhomogeneous dielectric continuum model for B-DNA

Thermodynamic and structural properties of the counterion atmosphere surrounding B-DNA are calculated by Monte Carlo simulation in a spatially inhomogeneous, but piecewise uniform, dielectric continuum cell model - the "barbarous" model. A boundary element formulation is implemented to study the sensitivity of these properties with respect to perturbations in the location of discontinuous dielectric boundaries relative to fixed and mobile charges. High concentrations are considered corresponding to the liquid crystalline hexagonally ordered phase of DNA. Primitive model results are verified against other simulation reports and a comparison of barbarous model predictions with experimental data is discussed. The internal energy, osmotic coefficient, radial distributions and the population ratio of counterions in the geometrically resolved major and minor grooves are all found to strongly depend on the dielectric boundary position. This suggests that a self-consistent development of the model should consider a free surface problem where the boundary is not specified a priori.

Calculation of affinities of peptides for proteins

Several methodologies were employed to calculate the Gibbs standard free energy of binding for a collection of protein-ligand complexes, where the ligand is a peptide and the protein is representative for various protein families. Almost 40 protein-ligand complexes were employed for a continuum approach, which considers the protein and the peptide at the atomic level, but includes solvent as a polarizable continuum. Five protein-ligand complexes were employed for an all-atom approach that relies on a combination of the double decoupling method with thermodynamic integration and molecular dynamics. These affinities were also computed by means of the linear interaction energy method. Although it generally proved rather difficult to predict the absolute free energies correctly, for some protein families the experimental ranking order was correctly reproduced by the continuum and all-atom approach. Considerable attention has also been given to correctly analyze the affinities of charged peptides, where it is required to judge the effect of one or more ions that are being decoupled in an all-atom approach to preserve electroneutrality. The various methods are further judged upon their merits.

A conserved arginine plays a role in the catalytic cycle of the protein disulphide isomerases

The pK(a) values of the CXXC active-site cysteine residues play a critical role in determining the physiological function of the thioredoxin superfamily. To act as an efficient thiol-disulphide oxidant the thiolate state of the N-terminal cysteine must be stabilised and the thiolate state of the C-terminal cysteine residue destabilised. While increasing the pK(a) value of the C-terminal cysteine residue promotes oxidation of substrates, it has an inhibitory effect on the reoxidation of the enzyme, which is promoted by the formation of a thiolate at this position. Since reoxidation is essential to complete the catalytic cycle, the differential requirement for a high and a low pK(a) value for the C-terminal cysteine residue for different steps in the reaction presents us with a paradox. Here, we report the identification of a conserved arginine residue, located in the loop between beta5 and alpha4 of the catalytic domains of the human protein disulphide isomerase (PDI) family, which is critical for the catalytic function of PDI, ERp57, ERp72 and P5, specifically for reoxidation. An examination of the published NMR structure for the a domain of PDI combined with molecular dynamic studies suggest that the side-chain of this arginine residue moves into and out of the active-site locale and that this has a very marked effect on the pK(a) value of the active-site cysteine residues. This intra-domain motion resolves the apparent dichotomy of the pK(a) requirements for the C-terminal active-site cysteine.

Identification by integrated computer modeling and light scattering studies of an electrostatic serum albumin-hyaluronic acid binding site

Dynamic light scattering and turbidimetry, carried out on solutions of hyaluronic acid (HA) and bovine or human serum albumin (SA) at fixed ionic strength (I), revealed a critical pH corresponding to the onset of HA-SA soluble complex formation. Subsequent reduction of pH below pH(c), corresponding to an increase in protein net positive charge, results in phase separation of the complex. The sensitivity of pH(c) to I indicated the primacy of electrostatic interactions in this process. Since pH(c) was always above the pK(a) of HA, these effects could be attributed to the influence of protein charge. The electrostatic potential around HSA was modeled using DelPhi (MSI) under pH, I conditions corresponding to incipient binding, phase separation, and noninteraction. At all incipient binding conditions (i.e., pH(c), at varying I), an identical region of positive potential 5 A from the protein van der Waals surface appeared. This unique domain intensified with a decrease in pH or I (corresponding to stronger binding), and diminished with an increase in pH or I (i.e., at noninteracting conditions). The size and low curvature of this domain could readily accommodate a 12 nm (decamer) sequence of HA. Simple electrostatic considerations indicate an electrostatic binding energy for the formation of this complex of ca. 1 kT, consistent with the condition of incipient complex formation. We suggest that such weak electrostatic binding may characterize nonspecific interactions for other protein-gylcosaminoglycan pairs.

Theoretical calculation of pKa reveals an important role of Arg205 in the activity and stability of Streptomyces sp. N174 chitosanase

Based on the crystal structure of chitosanase from Streptomyces sp. N174, we have calculated theoretical pK(a) values of the ionizable groups of this protein using a combination of the boundary element method and continuum electrostatics. The pK(a) value obtained for Arg(205), which is located in the catalytic cleft, was abnormally high (>20.0), indicating that the guanidyl group may interact strongly with nearby charges. Chitosanases possessing mutations in this position (R205A, R205H, and R205Y), produced by Streptomyces lividans expression system, were found to have less than 0.3% of the activity of the wild type enzyme and to possess thermal stabilities 4-5 kcal/mol lower than that of the wild type protein. In the crystal structure, the Arg(205) side chain is in close proximity to the Asp(145) side chain (theoretical pK(a), -1.6), which is in turn close to the Arg(190) side chain (theoretical pK(a), 17.7). These theoretical pK(a) values are abnormal, suggesting that both of these residues may participate in the Arg(205) interaction network. Activity and stability experiments using Asp(145)- and Arg(190)-mutated chitosanases (D145A and R190A) provide experimental data supporting the hypothesis derived from the theoretical pK(a) data and prompt the conclusion that Arg(205) forms a strong interaction network with Asp(145) and Arg(190) that stabilizes the catalytic cleft.

Protonation states and pH titration in the photocycle of photoactive yellow protein

Photoactive yellow protein (PYP) undergoes a light-driven cycle of color and protonation states that is part of a mechanism of bacterial phototaxis. This article concerns functionally important protonation states of PYP and the interactions that stabilize them, and changes in the protonation state during the photocycle. In particular, the chromophore pK(a) is known to be shifted down so that the chromophore is negatively charged in the ground state (dark state) even though it is buried in the protein, while nearby Glu46 has an unusually high pK(a). The photocycle involves changes of one or both of these protonation states. Calculations of pK(a) values and protonation states using a semi-macroscopic electrostatic model are presented for the wild-type and three mutants, in both the ground state and the bleached (I(2)) intermediate state. Calculations allowing multiple H-bonding arrangements around the chromophore also have been carried out. In addition, ground-state pK(a) values of the chromophore have been measured by UV-visible spectroscopy for the wild-type and the same three mutants. Because of the unusual protonation states and strong electrostatic interactions, PYP represents a severe test of the ability of theoretical models to yield correct calculations of electrostatic interactions in proteins. Good agreement between experiment and theory can be obtained for the ground state provided the protein interior is assumed to have a relatively low dielectric constant, but only partial agreement between theory and experiment is obtained for the bleached state. We also present a reinterpretation of previously published data on the pH-dependence of the recovery of the ground state from the bleached state. The new analysis implies a pK(a) value of 6.37 for Glu46 in the bleached state, which is consistent with other available experimental data, including data that only became available after this analysis. The new analysis suggests that signal transduction is modulated by the titration properties of the bleached state, which are in turn determined by electrostatic interactions. Overall, the results of this study provide a quantitative picture of the interactions responsible for the unusual protonation states of the chromophore and Glu46, and of protonation changes upon bleaching.

A self-consistent, microenvironment modulated screened coulomb potential approximation to calculate pH-dependent electrostatic effects in proteins

An improved approach is presented for calculating pH-dependent electrostatic effects in proteins using sigmoidally screened Coulomb potentials (SCP). It is hypothesized that a key determinant of seemingly aberrant behavior in pKa shifts is due to the properties of the unique microenvironment around each residue. To help demonstrate this proposal, an approach is developed to characterize the microenvironments using the local hydrophobicity/hydrophilicity around each residue of the protein. The quantitative characterization of the microenvironments shows that the protein is a complex mosaic of differing dielectric regions that provides a physical basis for modifying the dielectric screening functions: in more hydrophobic microenvironments the screening decreases whereas the converse applies to more hydrophilic regions. The approach was applied to seven proteins providing more than 100 measured pKa values and yielded a root mean square deviation of 0.5 between calculated and experimental values. The incorporation of the local hydrophobicity characteristics into the algorithm allowed the resolution of some of the more intractable problems in the calculation of pKa. Thus, the divergent shifts of the pKa of Glu-35 and Asp-66 in hen egg white lysozyme, which are both about 90% buried, was correctly predicted. Mechanistically, the divergence occurs because Glu-35 is in a hydrophobic microenvironment, while Asp-66 is in a hydrophilic microenvironment. Furthermore, because the calculation of the microenvironmental effects takes very little CPU time, the computational speed of the SCP formulation is conserved. Finally, results from different crystal structures of a given protein were compared, and it is shown that the reliability of the calculated pKa values is sufficient to allow identification of conformations that may be more relevant for the solution structure.