INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Protein stability and function depend on interactions among charged and dipolar groups. All proteins have ionized residues, many bind charged substrates, and changes in charge distribution often occur in the active site with reaction. Ionization state changes determine the pH dependence of protein function (Mellor et al., 1993; Dyson et al., 1997; Alexov et al., 2000) and stability (Hendsch et al., 1996; Spassov et al., 1997; Elcock and McCammon, 1998). Proton transfers between buried residues generate transmembrane proton gradients (Ferguson-Miller and Babcock, 1996; Alexov and Gunner, 1999; Luecke et al., 1999; Sham et al., 1999). Surface-charged residue distribution is critical for protein-protein association (Sheinerman et al., 2000). Computational tools that improve estimates of the energetics of the charge distribution within proteins will thus help us understand how protein function is derived from structure.

The analysis of protein charges and dipoles starts with water as a reference solvent (Warshel and Russell, 1984; Gilson and Honig, 1988; Gunner and Alexov, 2000). Water dipoles, being both polar and polarizable, reorient when a charge is added. Ions and dipoles have favorable interactions with water with its high dielectric constant, which are diminished in protein. In contrast, proteins are highly polar but not very polarizable molecules. Each backbone amide has a dipole greater than water, and more than 50% of all side chains are polar or ionizable. However, proteins have limited, nonuniform flexibility producing a heterogeneous response to charges.

The equilibrium distribution of side chain ionization states and position represents the balance of many interactions and so is hard to calculate. Methods to calculate the free energies of protein charge distributions are often tested by their ability to reproduce experimental pK_as (Bashford and Karplus, 1990; Beroza et al., 1991; Takahashi et al., 1992; Yang and Honig, 1993; Antosiewicz et al., 1994, 1996; Demchuk and Wade, 1996; Alexov and Gunner, 1997; Sham et al., 1997; Havranek and Harbury, 1999; Mehler and Guarnieri, 1999; Nielsen and Vriend, 2001). These benchmark calculations are not ideal because they include many values for surface residues, which have little interaction with protein. Also, many pK_as are not at pHs where the protein is functional or the structure determined, so calculations must also model pH-dependent structural perturbations. However, with more than 100 measured pK_as, calculations can try to match many values with one parameter set and methodology. As reliable computational methods become available, proteins where electrostatic forces play important functional roles can be analyzed (Alexov and Gunner, 1999; Demchuk et al., 2000; Dillet et al., 2000).

The Poisson-Boltzmann equation of continuum electrostatics provides a powerful framework for the analysis of electrostatic interactions in proteins (Warwicker and Watson, 1982; Gilson et al., 1985). These calculations take the distribution of dielectric constants, fixed charges, and mobile ions and return the electrostatic potential. Partial charges are placed at atoms in the protein structure. Water, with a dielectric constant of 80, is the major source of the reaction field (solvation) energy of protein dipoles and charges. Water also screens interatomic interactions in a manner that depends on the distances between protein atoms and the solvent. The Poisson-Boltzmann formalism incorporates the large impact of water easily and reasonably accurately. However, a central problem is how to assign the protein dielectric constant. The dielectric constant of a dried protein is measured to be near 4 (Takashima and Schwan, 1965). Molecular dynamics simulations show the protein inner core is not very polarizable with an effective dielectric constant near 4, whereas the solvent exposed flexible exterior has a value near 30 (Simonson and Brooks, 1996). Calculations using one, low value for protein (4) yield a poor match to benchmark pK_as, but there is significant improvement as _prot is increased to 20 (Antosiewicz et al., 1994, 1996; Demchuk and Wade, 1996).

The dielectric constant measures the response of a medium to changes in charge. In protein calculations the dielectric constant is a parameter that averages many kinds of reorganization so that they need not be included explicitly (Sham et al., 1997, 1998; Simonson, 1998, 2001; Gunner and Alexov, 2000; Schutz and Warshel, 2001). For example, changes in electronic polarization, backbone or side chain position, residue protonation state, and ion binding can occur as the pH is changed or a reaction proceeds. Electronic polarization may be fairly uniform throughout the protein, and the backbone motions may be small when there is no significant structural change. Thus, their contribution to the electrostatic energy may be estimated from an averaged dielectric constant rather than by explicit atomic calculations. However, much of the protein response is dependent on the local environment.

Several methods that remain based on continuum electrostatics add a nonuniform protein response to charges (for reviews, see Schutz and Warshel, 2001; Simonson, 2001). Different dielectric constants have been assigned to different amino acids (Voges and Karshikoff, 1998) or to different regions of the protein (Demchuk and Wade, 1996; Rocchia et al., 2001). Other methods consider explicit protein motions. Early methods for sampling multiple protein conformations either averaged interactions between different side chain atomic positions (You and Bashford, 1995) or averaged the ionization states calculated with different structures (Ripoll et al., 1996; Zhou and Vijayakumar, 1997). Hybrid methods simultaneously calculate ionization states of acidic and basic residues and the conformation states of surface side chains (Beroza and Case, 1996), hydroxyls and waters (Alexov and Gunner, 1997), and polar and ionizable side chains (Alexov and Gunner, 1999).

Other approaches include: adding explicit protein relaxation into the protein dipole Langevin dipole methodology with a linear response approximation (Sham et al., 1997, 1998); the use of a nonuniform response parameterized on fragment hydrophobic constants, which implicitly draws on the relationship between a group's polarity and its ability to partition into water (Mehler and Guarnieri, 1999); as well as an iterative mobile cluster approach for calculation of multiple-site ligand binding to flexible macromolecules (Spassov and Bashford, 1999).

The multiconformation continuum electrostatic (MCCE) method allows multiple positions of hydroxyl and water protons, side chain rotamers, and ligands in the calculation of the pH dependence of the ionization equilibria of titratable groups (Alexov and Gunner, 1997, 1999). An analysis of the pH titrations of 166 residues in 12 proteins is reported here. The benchmark root mean square errors of these titrations are comparable with those found with other approaches (Antosiewicz et al., 1994, 1996; Mehler and Guarnieri, 1999). Because improvement in pK_a calculations depends on the protein changing side chain position with pH, MCCE highlights motions that may be important for maintaining protein structure and function.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

The MCCE method

MCCE calculates the equilibrium conformation and ionization states of protein side chains, buried waters, ions, and ligands at a defined pH. Preselected choices for atomic positions and ionization states for side chains and ligands are used. Each specified side chain or ligand position and charge is called a conformer. Look-up tables are calculated of electrostatic and nonelectrostatic conformer self- and pair-wise interactions. Protein microstates have one conformer for each side chain, and ligand. Monte Carlo sampling establishes each conformer's probability in the Boltzmann distribution at a given pH or E_h (Alexov and Gunner, 1997, 1999). The method can compare mutated with wild-type proteins (Alexov et al., 2000) or protein with reactant or product bound (Alexov and Gunner, 1999). The MCCE procedure is divided into three stages: 1) selection of conformers (steps 1-3 below); 2) calculation of energy look up tables (steps 4-5); 3) Monte Carlo sampling and calculation of pK_as (steps 6-7). Early descriptions of the method can be found in Alexov and Gunner (1997, 1999).

Step 1: generate protein data file, striped of surface waters, complete with all protons except those on hydroxyls or waters

Step 2: identify residues with strong electrostatic pair-wise interactions

Step 3: generate composite protein structure

Step 4: electrostatic self- (reaction field) and pair-wise interactions

(1)

jk'

View larger version (15K): [in this window] [in a new window]   FIGURE 1   Dielectric boundary used for calculation of (A) pair-wise interactions and (B) reaction field energy. A2 is the only conformer with partial charges in the DelPhi calculation. Speckled region at prot = 4. Conformers within dashed lines and the surrounding area at  = 80.

Step 5: self- (torsion) and pair-wise (Lennard-Jones) nonelectrostatic energies for all conformers

Step 6: Monte Carlo sampling to determine the Boltzmann averaged conformer occupancy

(2)

Step 7: calculate residue pK_a and n values from conformer occupancies as a function of pH

(3)

Timing of the MCCE calculations

The analysis of Bovine pancreatic trypsin inhibitor (BPTI) (58 residues and 173 conformations) takes 3 h, whereas RNase H (155 residues and 608 conformers) takes 8 h on a single processor of an R10,000, 180 MHz SGI Origin 200. Most the time is spent in the DelPhi calculations. Monte Carlo sampling to determine all conformer occupancies at 15 pHs takes 20 min for BPTI and 1.5 h for RNase H.

Single conformer continuum electrostatics (SCCE) calculations

SCCE and MCCE calculations use the same protein structures with the same parameters for atomic charge and radius. In SCCE, waters are deleted and hydroxyl protons are placed in the most occupied MCCE position at pH 7. The most populated form of the neutral acidic residue at low pH is the only available neutral conformer. _prot is 4.

Atomic coordinates and experimental pK_a values

MCCE was used to calculate pK_as in 12 proteins (Table 1). The 166 measured values are from ¹H -¹³C heteronuclear resonance experiments unless otherwise specified. The experimental and calculated pK_a for each residue are reported in the supplementary material (SM I-X).

Barnase (SM-I) uses crystal structures 1A2P (Martin et al., 1999), 1B20 (Buckle et al., 1993), and the 20 NMR model structures (1BNR) (Bycroft et al., 1991). Experimental values were measured at low ionic strength (Oliveberg et al., 1995). The pK_a of His-18 was obtained from fluorescence measurements (Loewenthal et al., 1993).

BPTI (SM-II) uses crystal structure (4PTI) (Marquart et al., 1983) and 20 NMR structures (1PIT) (Berndt et al., 1992). pK_a from NMR experiments performed at 41°C were corrected to 25°C using the method described in March et al. (1982).

Intestinal bovine calcium-binding protein (calbindin, CabD) (SM-III) consists of crystal structure in presence of Ca²⁺ (3ICB) (Szebenyi and Moffat, 1986), 24 NMR structures for the holo-protein (1CDN) (Akke et al., 1995), and 33 NMR structures for the apoprotein (1CLB) (Skelton et al., 1995). pK_as of Lys with and without calcium is from Kesvatera et al. (1996). The acidic pK_as in the apoenzyme are from Kesvatera et al. (1999).

Rat T-lymphocyte adhesion glycoprotein (CD2) (SM-IV) uses crystal structure 1HNG (Jones et al., 1992). pK_as of 14 acidic residues determined at low (0.1 mM) and high (1.2 mM) salt concentration were averaged (Chen et al., 2000).

Hen egg-white lysozyme (HEWL) (SM-V) uses triclinic crystal structures (2LZT) (Ramanadham et al., 1981), 1B0D (Vaney et al., 1996), tetragonal crystal structure 1HEL (Wilson et al., 1992), 50 NMR structures (1E8L) (Schwalbe et al., 2001), and the averaged NMR structure (1HWA) (Smith et al., 1993). pK_as of basic residues are from Kuramitsu and Hamaguchi (1980), and values for acidic residues at 100 mM ionic strength are from Bartik et al. (1994).

Turkey ovomucoid inhibitor domain 3 (OMTKY3) (SM-VI) consists of OMTKY3 complexed with human leukocyte elastase (1PPF) (Bode et al., 1986), Streptomyces griseus proteinase B (3SGB) (Read et al., 1983), or -chymotrypsin (1CHO) (Fujinaga et al., 1987). The proteases are removed. There are 50 NMR structures (1OMU) (Hoogstraten et al., 1995). Acidic residue pK_as are measured at low (10 mM) and high ionic strength (1 M) (Schaller and Roberston, 1995). Basic residue pK_as are measured at 15 mM and 500 mM ionic strengths (Forsyth et al., 1998).

B1-immunoglobulin G-binding domain of protein G (ProtG) (SM-VII) uses crystal structure (1PGA) (Gallagher et al., 1994), 60 NMR structures (1GB1) (Gronenborn et al., 1991), and minimized NMR structure (2GB1) (Gronenborn et al., 1991). pK_as at 100 mM ionic strength are from Khare et al. (1997).

Ribonuclease A (RNase A) (SM-VIII) uses the crystal structure with sulfate (3RN3) (Howlin et al., 1989), ion free protein (7RSA) (Wlodawer et al., 1988), and 32 NMR structures (2AAS) (Santoro et al., 1993). pK_as at 200 mM NaCl (Rico et al., 1991; #4211). pK_as of His-12, -48, -105, and -119 in the absence of phosphate are an average of values from Ruterjans and Witzel (1969), Matthews and Westmoreland (1973), Walters and Allerhand (1980), and Quirk and Raines (1999); His-112 and His-119 pK_as determined at 20 mM and 100 mM phosphate (Meadows et al., 1969; Cohen et al., 1973). Sulfate is treated as phosphate.

Ribonuclease Hi (RNase H) (SM-IX) uses crystal structures 2RN2 (Katayanagi et al., 1992) and 1RNH (Yang et al., 1990) and Mg²⁺ containing form (1RDD) (Katayanagi et al., 1993) and 8 NMR derived structures (1RCH) (Yamazaki et al., 1997). pK_as are from Oda et al. (1993, 1994).

Ribonuclease T1 (RNase T) (SM-X) uses the crystal structure with vanadate (3RNT) (Kostrewa et al., 1989). Vanadate was treated as phosphate. Measured pK_as are from Inagaki et al. (1981) and McNutt et al. (1990).

Human Immunodeficiency Virus Type-1 Protease (HIV-1) uses the crystal structure (1HIV) (Thanki et al., 1992) with dihydroethylene-containing inhibitor, which is removed. pK_as are of the aspartyl dyad from inhibitor free enzyme kinetic studies (Hyland et al., 1991; Ido et al., 1991).

For residues in Barnase, OMTKY3, BPTI, Calbindin, and RNaseHi, there is more than one pK_a reported. The value used is either an average or the pK_a of the site closest to the atom losing the proton. Overall, the uncertainties are small ( 0.3 pK units). The n values (Eq. 3) were estimated from the published titration curves.

Even in these well-characterized proteins, only 166 of the 374 ionizable residues (Asp, Glu, His, Tyr, Lys, and Arg) have measured pK_as (Table 1). There are no Arg pK_as reported. Residues titrating in the same pH region influence each other. Thus, errors in the titration of Arg may influence Lys or Tyr. Likewise, errors in titration of missing Asp or Glu can impact the titration of the reported residues.

The following have incomplete titration curves that miss the plateau at low or high pH: Lys-34, Asp-7, and Asp-27 in OMTKY3, Lys-25 in CabD, Tyr-3, Tyr-33, Tyr-45, and Lys-13 in ProtG, Asp-93 and Glu-73 in Barnase, Tyr-23 in BPTI, Asp-48 and Asp-66 in HEWL, and His-114 in RNase H. Generally titration curves justify the reported values. There are eight residues where only a limiting value is available. This is taken as the true pK_a, which may overstate the error of the calculations (see supplementary material).

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Comparison of MCCE calculated and measured pK_as

The first pK_a calculations using detailed protein structures and energies derived from continuum electrostatics held the protein rigid, allowed only changes in residue ionization states, and used a low protein dielectric constant (_prot = 4) (Bashford and Karplus, 1990; Lim et al., 1991; Yang et al., 1993; Honig and Nicholls, 1995). SCCE does not match experimentally determined pK_as very well (Antosiewicz et al., 1994; Demchuk and Wade, 1996). However, SCCE provides good benchmarks for comparing different methods. SCCE root mean square (RMS) errors at _prot of 4 (Table 2) are similar to those found previously. MCCE and SCCE calculations with _prot = 4 are compared in Fig. 2 and Table 2. Experimental pK_as for each residue and values calculated with different Protein Data Bank files and _prot values are in the supplementary material (SM-I-X). For the 166 residues, the average absolute error decreases from 0.46 (SCCE) to 0.11 pH units (MCCE). The RMS error is reduced from 2.29 to 0.81 pH units. Less than 10% of the calculated pK_as differ from experiment by more than 1.5 pH units. With SCCE 40% of the pK_as are in error by this amount and 7% have errors greater than 4 pH units (Table 2).

View larger version (17K): [in this window] [in a new window]   FIGURE 2   Calculated versus measured pKas for (a) single conformer (SCCE) and (b) multiconformer (MCCE) continuum electrostatic calculations. prot = 4. Asp (), Glu (), Ctr (), His (), Ntr (), Tyr (), and Lys ().

Eighty-one residues are on the surface with no special interactions. However, 35 experimental pK_as are shifted by >1.0 pK unit from isolated amino acids (Table 3). Thirty-six residues are buried and have lost >3 pK units (4.08 kcal/mol) reaction field (solvation) energy. Fifty-three have ion pair interactions of >3 pK units. Thus, 85 residues are in one or more perturbed categories with many in more than one class. Surface, noninteracting residues are easiest to model. Their SCCE RMS error is only 1.57, and the MCCE error is 0.67 pH units. Of the unperturbed residues His have the largest RMS error because of His-124 in RNAse H (see below). Perturbed residues have larger errors with SCCE (RMS error 2.58 pH units). MCCE provides significant improvement (RMS error 1.09), but these residues pose the greatest challenge for calculation (Table 3). The largest errors come from residues in strong interactions (RMS error 1.17) as do some of the greatest improvements. One difficulty is that many ion pairs involve Arg for which there are no measured pK_as.

pK_as calculated with different structures of the same protein

Calculations should provide values that are the same for all nominally identical structures of the same protein. The average standard deviation with SCCE of the 32 residues in lysozyme structures (2LZT and 1B0D) is 1.4 pH units and only 0.5 pH units with MCCE (SM-V). Thus, MCCE conformational flexibility greatly diminishes the dependence on the initial structure. In addition, the RMS error for the 19 residues with known pK_as is cut in half.

The pK_as derived from crystal and solution structures were compared for 126 residues in eight proteins (Table 4). MCCE uses a rigid backbone and nonpolar side chains. The sampling of conformational space is increased in an ensemble of NMR structures, where each has somewhat different backbone and side chain positions. Thus, the average RMS variation of residue pK_as is 1.5 pH units when different NMR coordinate files are compared. However, the RMS error of the averaged pK_as is 0.83 pH units, smaller than for the crystal structures (0.94). Similar improvement was found in single conformer studies with _prot of 20 (Antosiewicz et al., 1996). Calculations with single, deposited average NMR structures have slightly larger errors (RMS error 1.02 pH units). However, the improvement with analysis of ensembles of NMR structures comes at the expense of more calculations. Here, 271 NMR structures were analyzed for the comparison of eight proteins.

Comparison of experimental and calculated n values

The Hill coefficient, n, in Eq. 3 measures how much residue titrations influence each other. Neighboring residues titrating in the same pH range have smaller n values. Experimental n values were obtained for 94 residues in CabD, Barnase, OMTKY3, CD2d1, and RNaseH, either from the published values or from the shape of published titration curves. Even in these small proteins, half of the residues have n values smaller than 0.85, whereas only two have values significantly larger than 1. The MCCE derived values (Eq. 3) are well correlated with those determined experimentally (Table 5).

The correspondence between calculated titration curves and MCCE site occupancies become worse as n decreases. This is because Eq. 3 provides only an approximation for the behavior of individual residues, missing the complex pH dependence of clusters of interacting residues (Onufriev et al., 2001).

Calculation of water and ion binding

Solvent exposed waters are removed, whereas buried waters are treated in the same detail as side chains with conformations allowing rotation about the oxygen. Each water has an extra conformer in bulk solvent with no interactions with the protein. The reaction field energy of free, solvated water controls the occupancy of interior sites. The water/hexane partition coefficient suggests a 3.56 pK units (4.84 kcal/mol) reaction field energy for an isolated water in bulk water (Wolfenden and Radzicka, 1994). For water to be bound, the remaining reaction field energy plus added interactions with the protein must be more favorable than this. There are 38 buried, crystallographic waters (Table 1). With a solution reaction field energy of 3.56 pK units, the average water site occupancy is only 33% (pH 7). The fraction increases to 77% as the penalty for removing the water from bulk was decreased to 0.87 pK units. Similar results were found in earlier MCCE calculations (Alexov and Gunner, 1999). In the group of proteins studied here only Asp-66 in HEWL has a pK_a that depends on bound water (see below). The indifference to water in this data set can be seen in the similarity of the pK_as derived from crystal structures and those using NMR structures where no waters are included. However, in other systems, bound water can have a large impact on the stability of buried charges (Gibas and Subramaniam, 1996; Garcia-Moreno et al., 1997; Alexov and Gunner, 1999).

Analysis of salt-bridges and corrections for strong electrostatic interactions

It is often difficult to obtain the correct pK_as for residues in ion pairs (Sindelar et al., 1998). Continuum electrostatics calculations often over-stabilize the ionization of charges, lowering acid pK_as and raising those of bases. MCCE with no correction yields pK_as for the 98 acids 0.59 pH units too low and for the 43 bases 0.54 pH units too high. An empirical function (SOFT) (Eqs. 4a and 4b, Fig. 3) was introduced to weaken strong pair-wise electrostatic interactions with other conformers (G) and with backbone and nonconformer containing side chains (G_{pol, i}).

(4a)

(4b)

With full strength interactions the global RMS error is 1.68 pH units (Tables 2 and 6), whereas it is 0.86 with the SOFT function. The average errors for the acids and bases are now 0.16 and 0.07 pH units, respectively. This function has little effect on the pK_a of most groups, but for 20 residues the error is decreased by >1.5 pK units (Fig. 4).

View larger version (17K): [in this window] [in a new window]   FIGURE 3   SOFT function applied to reduce strong pair-wise () conformer-conformer (Eq. 4a); () conformer-polar group (Eq. 4b) electrostatic interactions. Dashed line, initial G; dotted line, energy halved.

View larger version (20K): [in this window] [in a new window]   FIGURE 4   Comparison of the error in the pKas calculated with and without the SOFT function (Eq. 4). () acids; () bases. Along heavy line pKas are unchanged. Residues in the shaded region have larger errors with the SOFT function. Those in the white region have their errors reduced.

The SOFT function can be viewed as increasing the effective dielectric constant but only for relatively close neighbors (Fig. 3). The largest change halves the interaction energy, equivalent to an _prot of 8. There are several reasons that this may be necessary. Interaction energies are very sensitive to small changes in position when groups are close together. If structures report the closest allowable separation for salt-bridges as the time-averaged position then the interactions will be too large. Small excursions, weakening electrostatic interactions are not taken into account. In addition, MCCE uses Lennard-Jones parameters (SM-XI) that yield good values for hydrogen bonds. These may underestimate the cost of bringing ion pairs together, overestimating the electrostatic interactions. Lastly, the available conformers may not be optimal for stabilizing the charge-dipole state of the ion pair at low and high pH. Thus, the SOFT function represents a pragmatic solution to the problem of obtaining good values for ion pairs. However, future implementations of MCCE will aim to use enhanced local conformer flexibility to render this function unnecessary.

MCCE calculations with different protein dielectric constants

Residue pK_as were calculated at _prot of 4, 8, and 20 (Tables 2 and 3). With SCCE, there is marked improvement as _prot is increased (Antosiewicz et al., 1996; Demchuk and Wade, 1996). However, MCCE calculations are only moderately sensitive to _prot. The RMS error at _prot = 8 (0.69) is slightly smaller than with _prot = 4 (0.84) or _prot = 20 (0.70). Here, site titration is stabilized by either explicit conformation change at small _prot or by the averaged dielectric response as _prot increases. Fewer conformational changes are seen as _prot is increased.

With larger _prot, buried charges have a smaller G_rxn and smaller pair-wise interactions. Because most pair-wise interactions are favorable, changes in these two terms tend to offset each other. However, explicit conformational flexibility can stabilize a charge more than _prot of 20. For example, the MCCE pK_a of Asp-7 in OMTKY3 increases from 2.8, the correct value, to 5.0 as _prot increases from 4 to 20. In other cases the ionized form is too stable with _prot of 20 indicating the high dielectric constant either overestimates the local protein flexibility or weakens unfavorable pair-wise interactions.

Conformational flexibility in the MCCE calculations

In the proteins used here, 55% of the residues have multiple conformers either because they are polar (23%) or ionizable (32%). Overall 62% of the available conformers have some occupancy and 48% are occupied by more than 10% during the titration. Focusing on changes in polar residues (Asn, Gln, Ser, Thr, and waters), 17% do not change conformation between pH 0 and 14, 65% show small changes of 1% to 20%, whereas only 4% change occupancy by >50% (Fig. 5). Most motions are hydroxyl reorientation or rotation of Asn or Gln, which would not be seen by x-ray crystallography. However, these rotations rearrange hydrogen bond networks and so can have significant impact on calculated pK_as (Hooft et al., 1996; Nielsen and Vriend, 2001). In addition, several residues take heavy atom conformations different from that found in the protein data file at some point in the pH titration, and these have significant impact on calculated pK_as (see below).

View larger version (23K): [in this window] [in a new window]   FIGURE 5   Number of polar (Ser, Thr, Asn, Gln) conformers that have changed occupancy by the given percentage as the pH is changed from 0 to 14.

MCCE analysis of individual proteins

The benchmark pK_a calculations show the MCCE method yields as good fit to data as other current methods (Antosiewicz et al., 1996; Mehler and Guarnieri, 1999) (see below). The unique advantage of this method is that it follows explicit movement of side chains and waters during a titration and so highlights important conformational changes (Alexov and Gunner, 1999).

Lysozyme

Polar side chain motions

Small motions yield significant changes

CD2d1

Rotamer changes on the protein surface

View larger version (18K): [in this window] [in a new window]   FIGURE 6   (Top) Position of occupied conformers of interacting residues on the binding surface of rat CD2d1 (1HNG). The order of ionization is Glu-28, Glu-29, and then Glu-41 as the pH is raised. Glu-29 moves from neutral conformer A to ionized position B. The B conformer has more favorable interactions with Lys-43 and Arg-31. At low pH, the neutral Glu-41 is found with the proton at position A (78%) and B (18%). (Bottom) Conformer positions in the active site of RNaseA (3RN3). His-119, His-12, and the doubly ionized PO4 shown. The imidazole nitrogens (ND1 and NE2) in darker gray. His-12* and His-119*, found in the deposited structure, differ from His-12 and His-119 by rotation around the CBCG bond. His-119 positions A and B are included in the deposited structure file. At low pH ionized His-12* predominates (70% occupied) and His-119B* and His-119B are both found. Between pH 3 and 5, a small amount of ionized His-119A is occupied. Both His-119 proton tautomers (ND1 or NE2 protonated) are occupied. PO4 dissociates from the protein as His-12 looses its proton.

Calculation of pK_as in mutated proteins

Prot G

Accurate pK_as of ion pairs can rely on side chain motions

OMTKY3

Problems calculating pK_as of residues in clusters

RNaseA

Selection of His tautomers

Ion binding

Unstable calculated pK_as for residues in a buried ion pair