pH replica-exchange method based on discrete protonation states

We propose a new algorithm for obtaining proton titration curves of ionizable residues. The algorithm is a pH replica-exchange method (PHREM), which is based on the constant pH algorithm of Mongan et al. (J Comput Chem 2004;25:2038-2048). In the original replica-exchange method, simulations of different replicas are performed at different temperatures, and the temperatures are exchanged between the replicas. In our PHREM, simulations of different replicas are performed at different pH values, and the pHs are exchanged between the replicas. The PHREM was applied to a blocked amino acid and to two protein systems (snake cardiotoxin and turkey ovomucoid third domain), in conjunction with a generalized Born implicit solvent. The performance and accuracy of this algorithm and the original constant pH method (PHMD) were compared. For a single set of simulations at different pHs, the use of PHREM yields more accurate Hill coefficients of titratable residues. By performing multiple sets of constant pH simulations started with different initial states, the accuracy of predicted pK(a) values and Hill coefficients obtained with PHREM and PHMD methods becomes comparable. However, the PHREM algorithm exhibits better samplings of the protonation states of titratable residues and less scatter of the titration points and thus better precision of measured pK(a) values and Hill coefficients. In addition, PHREM exhibits faster convergence of individual simulations than the original constant pH algorithm.     

Charge effects on folded and unfolded proteins

We develop a theory for the effects of charge on the stabilization of globular proteins. The folding process is modeled as occurring through a fictitious intermediate state along a two-part thermodynamic pathway in which the molecule (i) increases its density and then (ii) rearranges its ionic groups to the protein surface. The equilibrium for the binding of protons in salt solutions is assumed to be driven by the electrical potential due to the charge distribution, in addition to the intrinsic binding affinity and bulk proton concentration. The potential is calculated for inside and outside a porous sphere model of the protein using the Poisson-Boltzmann relation, wherein the interior dielectric constant is taken to be a linear function of the chain density. The model predicts the slope of the titration curves for native myoglobin in agreement with experiments by Breslow and Gurd (1962). From the similar experiments on the unfolded state, and from the experiments of Privalov et al. (1986) on the intrinsic viscosity of the unfolded molecules, the theory shows that the unfolded state has a much higher density than a chain in a theta solvent and that the density increases with ionic strength. In addition, from the free energy of proton binding to the protein, we also calculate the electrostatic contributions to protein stability, a major contribution deriving from changes in ionization. We consider the example of the stability of myoglobin as a function of pH, ionic strength, and ionic groups buried in the native protein structure. We show that although maximum stability of most proteins should occur at their isoelectric point, the burial of nontitratable groups should lead to maximum stabilities at pH values other than the isoelectric point.     

A fast and accurate computational approach to protein ionization

We report a very fast and accurate physics-based method to calculate pH-dependent electrostatic effects in protein molecules and to predict the pK values of individual sites of titration. In addition, a CHARMm-based algorithm is included to construct and refine the spatial coordinates of all hydrogen atoms at a given pH. The present method combines electrostatic energy calculations based on the Generalized Born approximation with an iterative mobile clustering approach to calculate the equilibria of proton binding to multiple titration sites in protein molecules. The use of the GBIM (Generalized Born with Implicit Membrane) CHARMm module makes it possible to model not only water-soluble proteins but membrane proteins as well. The method includes a novel algorithm for preliminary refinement of hydrogen coordinates. Another difference from existing approaches is that, instead of monopeptides, a set of relaxed pentapeptide structures are used as model compounds. Tests on a set of 24 proteins demonstrate the high accuracy of the method. On average, the RMSD between predicted and experimental pK values is close to 0.5 pK units on this data set, and the accuracy is achieved at very low computational cost. The pH-dependent assignment of hydrogen atoms also shows very good agreement with protonation states and hydrogen-bond network observed in neutron-diffraction structures. The method is implemented as a computational protocol in Accelrys Discovery Studio and provides a fast and easy way to study the effect of pH on many important mechanisms such as enzyme catalysis, ligand binding, protein-protein interactions, and protein stability.     

Prediction of titration properties of structures of a protein derived from molecular dynamics trajectories.

This paper explores the dependence of the molecular dynamics (MD) trajectory of a protein molecule on the titration state assigned to the molecule. Four 100-ps MD trajectories of bovine pancreatic trypsin inhibitor (BPTI) were generated, starting from two different structures, each of which was held in two different charge states. The two starting structures were the X-ray crystal structure and one of the solution structures determined by NMR, and the charge states differed only in the ionization state of N terminus. Although it is evident that the MD simulations were too short to sample fully the equilibrium distribution of structures in each case, standard Poisson-Boltzmann titration state analysis of the resulting configurations shows general agreement between the overall titration behavior of the protein and the charge state assumed during MD simulation: at pH 7, the total net charge of the protein resulting from the titration analysis is consistently lower for the protein with the N terminus assumed to be neutral than for the protein with the N terminus assumed to be charged. For most of the ionizable residues, the differences in the calculated pKaS among the four trajectories are statistically negligible and remain in good agreement with the data obtained by crystal structure titration and by experiment. The exceptions include the N terminus, which responds directly to the change of its imposed charge; the C terminus, which in the NMR structure interacts strongly with the former; and a few other residues (Arg 1, Glu 7, Tyr 35, and Arg 42) whose pKaS reflect the initial structure and the limited trajectory lengths. This study illustrates the importance of the careful assignment of protonation states at the start of MD simulations and points to the need for simulation methods that allow for the variation of the protonation state in the calculation of equilibrium properties.     

Measurement of the mitochondrial membrane potential and pH gradient from the redox poise of the hemes of the bc1 complex

The redox potentials of the hemes of the mitochondrial bc(1) complex are dependent on the proton-motive force due to the energy transduction. This allows the membrane potential and pH gradient components to be calculated from the oxidation state of the hemes measured with multi-wavelength cell spectroscopy. Oxidation states were measured in living RAW 264.7 cells under varying electron flux and membrane potential obtained by a combination of oligomycin and titration with a proton ionophore. A stochastic model of bc(1) turnover was used to confirm that the membrane potential and redox potential of the ubiquinone pool could be measured from the redox poise of the b-hemes under physiological conditions assuming the redox couples are in equilibrium. The pH gradient was then calculated from the difference in redox potentials of cytochrome c and ubiquinone pool using the stochastic model to evaluate the ΔG of the bc(1) complex. The technique allows absolute quantification of the membrane potential, pH gradient, and proton-motive force without the need for genetic manipulation or exogenous compounds.     

Hydrophobic potential of mean force as a solvation function for protein structure prediction

We have developed a solvation function that combines a Generalized Born model for polarization of protein charge by the high dielectric solvent, with a hydrophobic potential of mean force (HPMF) as a model for hydrophobic interaction, to aid in the discrimination of native structures from other misfolded states in protein structure prediction. We find that our energy function outperforms other reported scoring functions in terms of correct native ranking for 91% of proteins and low Z scores for a variety of decoy sets, including the challenging Rosetta decoys. This work shows that the stabilizing effect of hydrophobic exposure to aqueous solvent that defines the HPMF hydration physics is an apparent improvement over solvent-accessible surface area models that penalize hydrophobic exposure. Decoys generated by thermal sampling around the native-state basin reveal a potentially important role for side-chain entropy in the future development of even more accurate free energy surfaces.     

Evaluation of linked protonation effects in protein binding reactions using isothermal titration calorimetry.

A theoretical development in the evaluation of proton linkage in protein binding reactions by isothermal titration calorimetry (ITC) is presented. For a system in which binding is linked to protonation of an ionizable group on a protein, we show that by performing experiments as a function of pH in buffers with varying ionization enthalpy, one can determine the pK(a)'s of the group responsible for the proton linkage in the free and the liganded states, the protonation enthalpy for this group in these states, as well as the intrinsic energetics for ligand binding (delta H(o), delta S(o), and delta C(p)). Determination of intrinsic energetics in this fashion allows for comparison with energetics calculated empirically from structural information. It is shown that in addition to variation of the ligand binding constant with pH, the observed binding enthalpy and heat capacity change can undergo extreme deviations from their intrinsic values, depending upon pH and buffer conditions.     

Incorporating protein conformational flexibility into the calculation of pH-dependent protein properties.

A method for combining calculations of residue pKa's with changes in the position of polar hydrogens has been developed. The Boltzmann distributions of proton positions in hydroxyls and neutral titratable residues are found in the same Monte Carlo sampling procedure that determines the amino acid ionization states at each pH. Electrostatic, Lennard-Jones potentials, and torsion angle energies are considered at each proton position. Many acidic and basic residues are found to have significant electrostatic interactions with either a water- or hydroxyl-containing side chain. Protonation state changes are coupled to reorientation of the neighboring hydroxyl dipoles, resulting in smaller free energy differences between neutral and ionized residues than when the protein is held rigid. Multiconformation pH titration gives better agreement with the experimental pKa's for triclinic hen egg lysozyme than conventional rigid protein calculations. The hydroxyl motion significantly increases the protein dielectric response, making it sensitive to the composition of the local protein structure. More than one conformer per residue is often found at a given pH, providing information about the distribution of low-energy lysozyme structures.     

Probing protein electrostatic interactions through temperature/reduction potential profiles

 The change in the equilibrium reduction potentials of the iron-sulfur proteins, Pyrococcus furiosus rubredoxin and P. furiosus ferredoxin, and heme protein, horse cytochrome c, has been calculated as a function of temperature using a numerical solution to the Poisson-Boltzman equation. Working curves for different internal dielectric constants were generated to best reproduce experimental observation. Based on a comparison of the experimental and simulated change in reduction potential with temperature, it is concluded that the dielectric constant of proteins is temperature-dependent and varies from protein to protein. For example, the temperature-dependent reduction potential of cytochrome c can only be simulated using a different temperature-dependent dielectric constant for each oxidation state, but this was not the case for rubredoxin or ferredoxin. The role of changes in ionization states of cytochrome c at alkaline pHs, where the reduction potential is known to be pH-dependent at room temperature, is also discussed in terms of electrostatic interaction energies as a function of temperature. It appears that temperature/reduction potential profiles may provide a direct method for measuring relative changes in internal protein dielectric constants. Received: 29 April 1996 / Accepted: 1 August 1996     

The calculated free energy effects of 5-methyl cytosine on the B to Z transition in DNA

We have examined the free energy effects of 5-methylation of cytosine on the B in equilibrium Z conformational equilibrium in DNA. Free energy differences were calculated using the free energy perturbation approach, which uses an easily derived equation from classical statistical mechanics to relate the free energy difference between two states to the ensemble average of the potential energy difference between the states. Calculations were carried both in explicit solvent and (for comparison) in vacuo. The free energy values obtained for the explicit solvent systems are total free energies, with contributions from all parts of the system (solvent + solute), and so are relevant to the B in equilibrium Z transitions observed under real (physiological) conditions. We calculate that in solution, methylation makes the B in equilibrium Z transition more favorable by about -0.4 kcal/mole base pair (bp) in free energy. This value compares well with approximate experimentally derived values of about -0.3 kcal/mole-bp. We also discuss a method for determining the free energy difference between conformational states poorly maintained by a potential energy model. Finally, the effects of methylation on the melting temperature of DNA are examined.     

Development of an equation for potentiometric titration of polyelectrolytes using a periodic lattice model. A numerical analysis of interactions among ionizable groups

A numerical method is presented for analysing the potentiometric titration behavior of linear polyelectrolytes. A polyelectrolyte molecule is treated as a one-dimensional lattice containing a large number of lattice points, each of which has an identical ionizable group. In this method, the polyelectrolyte model lattice is divided into identical repeating unit cell systems with a finite number of ionizable groups to calculate the thermodynamic partition function of the system of polyelectrolyte solution. The electrostatic interaction between ionized groups adopted in the present study is given by the Debye-Hückel type screened Coulombic potential. The titration behavior as well as several thermodynamic quantities is derived by making a canonical ensemble summation of all states in a unit cell system under an appropriate periodic boundary condition. This method serves as a model of the two-step ionization often appearing in polyions with strong neighboring interactions such as homo- and copolymers of maleic acid. Several characteristics of the titration behavior of these polyelectrolytes are well reproduced by using a lattice model with pairwise intervals, including the effects of hydrogen bond formation and change in dielectric constant of the medium around an ionizable group. In addition, this method is valid for a more detailed analysis of the titration behavior of polyelectrolytes with various kinds of arrangements of charged groups.     

The influence of a transmembrane pH gradient on protonation probabilities of bacteriorhodopsin: the structural basis of the back-pressure effect

Bacteriorhodopsin pumps protons across a membrane using the energy of light. The proton pumping is inhibited when the transmembrane proton gradient that the protein generates becomes larger than four pH units. This phenomenon is known as the back-pressure effect. Here, we investigate the structural basis of this effect by predicting the influence of a transmembrane pH gradient on the titration behavior of bacteriorhodopsin. For this purpose we introduce a method that accounts for a pH gradient in protonation probability calculations. The method considers that in a transmembrane protein, which is exposed to two different aqueous phases, each titratable residue is accessible for protons from one side of the membrane depending on its hydrogen-bond pattern. This method is applied to several ground-state structures of bacteriorhodopsin, which residues already present complicated titration behaviors in the absence of a proton gradient. Our calculations show that a pH gradient across the membrane influences in a non-trivial manner the protonation probabilities of six titratable residues which are known to participate in the proton transfer: D85, D96, D115, E194, E204, and the Schiff base. The residues connected to one side of the membrane are influenced by the pH on the other side because of their long-range electrostatic interactions within the protein. In particular, D115 senses the pH at the cytoplasmic side of the membrane and transmits this information to D85 and the Schiff base. We propose that the strong electrostatic interactions found between D85, D115, and the Schiff base as well as the interplay of their respective protonation states under the influence of a transmembrane pH gradient are responsible for the back-pressure effect on bacteriorhodopsin.     

A potentiometric method of quantitative analysis of rifampicin

A potentiometric procedure for assay of rifampicin was developed. The procedure implies titration of rifampicin as a monofunctional acid by sodium hydroxide solution (0.1 mol/l) in 75 per cent aqueous methanol. The constant ionic strength of the solution is provided by addition of KCl until its concentration is 0.1 mol/l, the titrant concentration being 10 times higher than the antibiotic concentration in the solution. This provides a precise determination of the concentration ionization constant of the antibiotic as a monofunctional acid (pKa 7.33 +/- 0.01) and an insignificant dilution of the antibiotic solution during the titration promoting precise and reproducible results. The procedure error is 0.20 per cent. The variation coefficient is 0.27 per cent.     

Can current force fields reproduce ring puckering in 2-O-sulfo-α-l-iduronic acid? A molecular dynamics simulation study

The monosaccharide 2-O-sulfo-α-l-iduronic acid (IdoA2S) is one of the major components of glycosaminoglycans. The ability of molecular mechanics force fields to reproduce ring-puckering conformational equilibrium is important for the successful prediction of the free energies of interaction of these carbohydrates with proteins. Here we report unconstrained molecular dynamics simulations of IdoA2S monosaccharide that were carried out to investigate the ability of commonly used force fields to reproduce its ring conformational flexibility in aqueous solution. In particular, the distribution of ring conformer populations of IdoA2S was determined. The GROMOS96 force field with the SPC/E water potential can predict successfully the dominant skew-boat to chair conformational transition of the IdoA2S monosaccharide in aqueous solution. On the other hand, the GLYCAM06 force field with the TIP3P water potential sampled transitional conformations between the boat and chair forms. Simulations using the GROMOS96 force field showed no pseudorotational equilibrium fluctuations and hence no inter-conversion between the boat and twist boat ring conformers. Calculations of theoretical proton NMR coupling constants showed that the GROMOS96 force field can predict the skew-boat to chair conformational ratio in good agreement with the experiment, whereas GLYCAM06 shows worse agreement. The omega rotamer distribution about the C5-C6 bond was predicted by both force fields to have torsions around 10°, 190°, and 360°.     

Design of an Efficient Inhibitor for the Influenza A Virus M2 Ion Channel

Influenza A virus is capable of rapidly infecting large human populations, warranting the development of novel drugs to efficiently inhibit virus replication. A transmembrane ion channel formed by the M2 protein plays an important role in influenza virus replication. A reasonable approach to designing an effective antivirus drug is constructing a molecule that binds in the M2 transmembrane proton channel, blocks H+ proton diffusion through the channel, and thus the influenza A virus cycle. The known anti-influenza drugs amantadine and rimantadine have a weak effect on influenza A virus replication. A new class of positively charged molecules, diazabicyclooctane derivatives with a constant charge of +2, was proposed to block proton diffusion through the M2 ion channel. Molecular dynamics simulations were performed to study the temperature fluctuations in the M2 structure, and ionization states of histidine residues were established at physiological pH values. Two types of diazabicyclooctane derivatives were analyzed for binding with the M2 ion channel. An optimal structure was determined for a blocker to most efficiently bind with the M2 ion channel and block proton diffusion. The new molecule is advantageous over amantadine and rimantadine in having a positive charge of +2, which creates a positive electrostatic potential barrier to proton transport through the M2 ion channel in addition to a steric barrier.     

Potentiometric titration of poly-S-carboxyethyl-L-cysteine

Potentiometric titration curves have been determined for aqueous solutions of poly-S-carboxyethyl-L -cysteine, which is subject to the β-coil transition by a change in pH. Reversibility and time dependence of the titration curves are examined by different methods in order to establish the conditions for obtaining equilibrium curves. The β-coil transition is manifest, at some region on the equilibrium titration curve, if pH – log (α/1 – α) is plotted against α. Assuming a value, 4.00, for pK_int, the free-energy change for the β-coil transition of uncharged polymer has been evaluated from the extrapolation of the observed titration curves and is found to depend on the ionic strength and polymer concentration. The Henderson-Hasselbach plot of the titration curve yields clearer distinction between the β-form and random coil, and it permits estimation of the content of β-form at, a given pH. Comparison of the conformational titration curve with the circular dichroic measurements leads to a value of ?10,000° for [θ]₂₂₃ for the pure β-structure. Precipitation which occurs at low degrees of ionization and, especially, at high ionic strength does not reveal any discontinuous change of the titration curve, which suggests that, the degree of ionization of the precipitated β-form is not very different from that in solution.