INTRODUCTION

TOP ABSTRACT INTRODUCTION THEORETICAL BASIS AND... RESULTS AND DISCUSSION CONCLUSIONS APPENDIX REFERENCES

The electrostatic interaction between two charges in a dielectric medium decreases rapidly as the separation between the two entities is increased. This is a result of the dielectric medium's intrinsic ability to shield the electrostatic field arising from the charge. There is compelling experimental and theoretical evidence that the general functional form of this screening is sigmoidal. Early experiments on organic acids (Debye and Pauling, 1925; Schwarzenbach, 1936; Webb, 1926) were collected together and presented as mean values in tabular form (Conway et al., 1951) that produced a smooth, sigmoidal screening as a function of separation with an asymptotic value of about 80. Experiments on water as a dielectric medium indicate that the electrostatic field, at loci only two solvation layers (6 Å) away from a charge, is dramatically damped, almost approaching the bulk dielectric screening of  = 80 (Harvey and Hoekstra, 1972; Pennock and Schwan, 1969; Takashima and Schwan, 1965). When two charges are in close proximity, however, the dielectric screening rapidly diminishes approaching the vacuum value because there is no intervening dielectric. These conditions are simultaneously satisfied with a radially dependent sigmoidal screening function. Additionally, it has been demonstrated that the generalized Born-surface area continuum solvation model also gives rise to sigmoidal screening (Mehler, 1996b). Thus, the general sigmoidal screening form applies to both high and low dielectric media. Theoretically, the Lorentz-Debye-Sack (LDS) theory of polar solvation also leads to sigmoidally screened electrostatic fields (Debye, 1929; Lorentz, 1952; Sack, 1926; Sack, 1927). LDS theory has been used to develop radial dielectric permittivity profiles to describe ion and dipole field sources (Ehrenson, 1989) and to derive an expression for the hydration energies of spherical ions (Bucher and Porter, 1986) from the integral form of the Born Equation (Born, 1920). Subsequently, Ehrenson (1989) showed that by using Bötcher's analytic expression (Böttcher, 1938) of Onsager's theory (Onsager, 1936), the reaction field corrected LDS theory for both dipolar and asymmetric ion field sources also yielded sigmoidal screening. (For a review of LDS theory and its relation to protein electrostatics, see Mehler, 1996a.) In various applications (Collura et al., 1994; Mehler, 1996a; Mehler and Solmajer, 1991), it has been shown that sigmoidal screening leads to the reliable prediction of properties dependent on electrostatic effects in proteins, and therefore is an attractive formulation for the development of implicit solvent models (Guarnieri et al., 1998; Hassan et al., 1999).

Recently, a self-consistent approach based on using sigmoidally screened Coulomb potentials (SCP) for calculating pH-dependent properties in proteins has been derived and applied to several systems (Mehler, 1996b). The results indicated that the reliability of the method is similar to that achieved by the finite difference Poisson-Boltzmann (FDPB), but requires substantially less computing effort. Interestingly, although these various methods predict a majority of the pK_a shifts correctly, they were consistent in the occasional dramatic errors on singularly important residues. The divergent pK_a shifts of the two titratable acids of hen egg white lysozyme, Glu-35 and Asp-66, are particularly problematic. The divergent pK_a shifts cannot be simply explained by the buried fraction (Demchuk and Wade, 1996), since the two titratable groups are nearly completely buried. Thus, it appears more likely that the unique local environment around each residue is a key determinant for the pK_a of Glu-35 to be two pK units higher than its water value, whereas in Asp-66 the pK_a is more than two units lower. Ponnuswamy et al. (1980), showed that at a detailed level, the local environmental hydrophobicity that surrounds each residue is distinct and can be quite different from all other residues, implying that a different screening would be required for each site (for PB calculations, each site would require its own internal permittivity value) as has been pointed out by Warshel (King et al., 1991).

In this work we hypothesize that a key source of single residue aberrations in pK_a is a result of the special protein microenvironment that exists around the singular residue. To test this hypothesis, a strategy has been developed to quantitatively characterize the hydrophilicity or hydrophobicity of the microenvironment around each titratable group, and then use this property to modulate the components of the electrostatic free energy operative in proteins. The advantage of this approach is that by using physicochemical properties independent of the parameters required to calculate the electrostatic effects, one might gain additional insight into the relation between protein structure and the forces that control its function. Accounting for every local microenvironment in the entire protein structure could become prohibitively costly. Therefore, to preserve the computational speed of the SCP, the quantitative calculation of the microenvironments has been implemented in a way that requires virtually no additional CPU time.

In the next section of this paper the method for quantifying the microenvironments is presented, and the dielectric screening function used previously (Mehler, 1996a; Mehler and Eichele, 1984) is reformulated to account for the local variation in the dielectric properties of proteins. The algorithm is applied to calculating the pK_a in bovine pancreatic trypsin inhibitor (BPTI), triclinic hen egg white lysozyme (HEWL), ribonuclease A (RnaseA), ribonuclease T1 (RnaseT), ribonuclease HI (RnaseHI), calbindinD9k (Cab) and the B1 immunoglobulin G-binding domain of protein G (ProtG). These seven proteins provide a data set of 103 measured pK_a values that have been used to develop the parametrization of the algorithm. Subsequently, the method is tested on turkey ovomucoid inhibitor domain 3, HIV protease, and other crystal forms of some of the proteins in the test set. Finally, the results obtained here are compared with other calculations and potential improvements for the quantification of the microenvironments are considered.

	THEORETICAL BASIS AND CHARACTERIZATION OF MICROENVIRONMENTS

TOP ABSTRACT INTRODUCTION THEORETICAL BASIS AND... RESULTS AND DISCUSSION CONCLUSIONS APPENDIX REFERENCES

Calculation of the electrostatic free energy

The calculation of the pK_a of titratable groups in proteins is greatly simplified by considering the thermodynamic cycle (Bashford and Karplus, 1990; Bashford and Karplus, 1991; Warshel, 1981)

(1)

where (p) and (s) refer to protein or solvent, respectively. Thus, the pK_a of group A in the protein can be calculated from its pK_a in model solvent (water) and the additional changes in free energy that arise when the titratable group is transferred from water into the protein. Here, only the electrostatic contributions are considered and pK_a(p) is obtained from the relation

(2)

where w_A^int is the interaction free electrostatic energy of the charged group in the field of all the other groups in the protein, and w_A^tr is the change in self-energy on transferring the group from water into the protein.

In Mehler (1996b), it was shown that with an SCP of the form

(3)

where q_j is the charge at point r_j, r_j = |r  r_j| and D(r) is a screening function, the total electrostatic free energy of a system consisting of M groups could be expressed as the sum of the interaction term and a transfer contribution as

(4)

where A may refer to subgroups of each amino acid residue and other groups contributing to the total charge of the system. As discussed in the introduction, both experimental results and theoretical considerations indicated that D(r) should have a sigmoidal form (for details, see the Appendix) as originally proposed by Mehler and Eichele (1984).

For the calculation of the pK_a, it is not necessary to evaluate the total electrostatic free energy, but only the contribution from each titratable group in the field due to all the other groups in the protein. In most pK_a calculations the full titration charge is placed at one or more fixed positions in the protonatable moiety of the N titratable residues, and the equilibrium charge state is calculated from the distribution of the charge over the 2^N ionization microstates (Bashford and Karplus, 1991; Beroza et al., 1991; Gilson, 1993). In the procedure developed in Mehler (1996b), the assignment of the ionization charge over the atoms of the protonatable moiety was determined variationally to allow it to respond to the protein environment. To reduce computing time, the direct determination of the 2^N charge microstates was bypassed by coupling the total variational ionization charge of each group, at the given pH, to the Henderson-Hasselbalch equation. In this paper the method presented in Mehler (1996b) is still used, but has been modified; the details of the modified approach are discussed in the Appendix.

The calculations reported in Mehler (1996b) were based on a single screening function that had been parametrized earlier on the basis of pK_a shifts in bifunctional organic acids and bases (D_wds defined in Mehler and Eichele, 1984). Thus, the only parameter introduced in Mehler (1996b) was a scaling of the interaction energy of nearby charges. To improve and generalize the approach followed in Mehler (1996b), it is hypothesized that the screening for groups deeply buried in the protein should be less than for solvent exposed groups. This reflects the commonly accepted view that the protein interior is more hydrophobic than water. To implement this hypothesis, two new screening curves were constructed by redefining D₀ and  (see Appendix, Eq. A1) so that, for small separations, one curve would provide more screening than D_wds, used in Mehler (1996b), whereas the other would provide less screening. These two screening functions, D₁ and D₂, as well as the screening function in water, D_w, are shown in Fig. 1 (D₃ will be discussed below) and the values of D₀ and  are given in Table I. Although these two screening functions appear to be very similar, for small values of r (see inset, Fig. 1) the energy is quite sensitive to small changes in the screening. Thus, at r = 2 Å, the difference between the two curves for unit charges amounts to nearly 4.5 pK units, while, at r = 3 Å, the difference is 1.4 pK units. In contrast, at r = 8 Å, the difference is only 0.09 pK units and continues to decrease with increasing r. Thus, for small separations between groups, the contributions to the interaction energy are quite sensitive to which screening function is used. This is also the case for the transfer energy, since the Born radii, R_a in Eq. A2, lie between about 1.2 Å and 2.6 Å. Which screening function is used for a given titratable group depends on the fraction of the group's solvent accessible surface area that is buried in the protein, i.e., the buried fraction. When these two screening functions are incorporated into the calculation, there was some improvement in the pK_a values, but the pK_a of Glu-35 in lysozyme was still too small by more than one pK unit.

View larger version (22K): [in this window] [in a new window]   FIGURE 1   Screening functions defined in Table 1. Solid line, Dw; short dashed line, D1; long dashed line, D2; and dots, D3. The inset is a detail from 0.5 to 3 Å. (See the Appendix for details).

Quantification of the microenvironment

Analysis of the local environments around Glu-35 and Asp-66 in lysozyme shows that for the former the nearest neighbors are largely hydrophobic residues, whereas for the latter residue the nearest neighbors are much more hydrophilic. The increased hydrophobicity around Glu-35 has consequences for both components of the electrostatic free energy: the energy transfer term will become more positive (more unfavorable), which can be effected by a decrease in R^p (see Eq. A6). The screening of the Coulomb potential should also decrease, especially for interactions with nearby atoms so that they are larger in magnitude relative to a less hydrophobic region. However, this will be partially offset by the fairly small partial charges on atoms in hydrophobic groups. Since w^int is usually negative, whereas w^tr will be positive for hydrophobic regions, the resulting shift in pK_a will depend on the subtle balance between these two quantities. The transfer of Asp-66 from water to a relatively hydrophilic environment, where D_p(R^p) is larger, will result in a smaller transfer energy (less unfavorable; it could be negative if the local environment is more hydrophilic than water), and the Coulombic screening should be greater. These considerations clearly suggest that the relevant quantity for characterizing the local environment, here called the microenvironment, is the hydrophobicity or hydrophilicity (here we use Hpy to refer to these properties and their values in particular cases) of the residues or their fragments near the given titratable group.

There are a large number of Hpy scales available for evaluating the hydrophobicity of amino acid residues (Cornette et al., 1987). Generally these scales have been used to evaluate the hydrophobicity of a particular residue, or at least its side chain. Here, the hydrophobicity of the particular residue is not required, but instead, because the response of the residue to changes in the local environment is the quantity of interest, it is a quantitative measure of the hydrophobicity/hydrophilicity of the microenvironment that is needed. A hydrophobicity scale based on atoms or small fragments is most appropriate for the algorithm proposed here because it allows the microenvironment to be defined in a simple, distance dependent way. In addition, this type of scale can account for the fact that side chains generally do not have a single Hpy value along their entire length, nor do all the atoms of a side chain necessarily belong to a given microenvironment. A further important advantage of a fragment-based scale is that the presence of coenzymes, prosthetic groups, ligands, or nucleic acids can all be accounted for with a consistent scale. From these considerations, the Rekker Fragmental Hydrophobic Constants (Rekker, 1977; Rekker, 1979) were selected to evaluate the hydrophobicities of the microenvironments in the present calculations.

The Rekker scale was the first fragmental system developed and is based on assigning hydrophobicity parameters and correction factors to small functional fragments to calculate log P values (for details, see Rekker, 1977). This scale has been widely used in the pharmaceutical industry as a means to calculate partition coefficients of potential drug candidates. Here, we use it not to obtain the hydrophobicity/hydrophilicity of a particular residue, but to obtain this quantity for the local protein environment. This is done by (1) identifying fragments of all the residues that are in the neighborhood of the titratable residue, (2) assigning the Rekker coefficients to these fragments, (3) summing them to determine the local hydrophobicity/hydrophilicity of the region around the titratable residue, and (4) assigning an appropriate screening for this environment. Three other fragmental systems have been proposed, of which two are based on functional fragments (Leo et al., 1975; Suzuki and Kudo, 1990) and one is atom based (Ghose and Crippen, 1986). All the fragmental systems that have been proposed to date are calibrated on the basis of partition coefficients (log P). The parameters are therefore completely independent of protein structural information. Thus, their successful use to describe local protein environments will give some indication of the universality of the hydrophobic concept and the transferability of scales developed from different physical concepts. A recent report compared the four methods and showed that the Rekker system was one of the most accurate available (Mannhold et al., 1995).

In this initial application the original values of the hydrophobic fragmental constants and correction factors were used (Rekker, 1977). For each fragment the values were partitioned to the atoms belonging to that fragment, and to simplify the calculations, net charge on a specific group was not accounted for. Then the microenvironment can be defined on the basis of a distance criterion, and only atoms satisfying it are included in the calculation of the total Hpy around the given group. The initial set of values that have been used in the calculations are tabulated in Table 2. Note that the water value was estimated from the fragmental constants of aliphatic OH and H_neg (Rekker, 1977).

The microenvironment of any group is defined as the region that lies within a sphere of radius 4.25 Å from any (nonhydrogenic) atom belonging to the group. This radius was chosen because it essentially includes the first shell around each atom, but does not extend to atoms beyond this. In addition, this value is in reasonable agreement with the criterion used by Ponnuswamy et al. (1980). In analogy with the formula used to calculate log P from the fragmental hydrophobic constants (Rekker, 1977), the Hpy value for each microenvironment of a group is calculated from the formula

(5)

In Eq. 5, RFHC_b is the contribution of atom b to the fragment's fragmental hydrophobic constant as given in Table 2, N_A is the number of atoms in group A and
to ensure that each atom in the microenvironment is counted only once. Schematic representations of the microenvironments around Glu-25 and Asp-66 in lysozyme are presented in Fig. 2. These diagrams show that in most cases only some atoms of a given amino acid residue are located within 4.25 Å of one or more atoms of the group in question. Moreover, these arrangements of the atoms in space clearly show how the protein architecture has evolved to create the hydrophobic microenvironment around Glu-35 and the hydrophilic microenvironment around Asp-66.

View larger version (21K): [in this window] [in a new window]   FIGURE 2   Residues defining the microenvironments of Glu-35 and Asp-66 in HEWL. The number under each residue is its contribution to the total hydrophobicity or hydrophilicity. Atoms within 4.25 Å of any atom of the titratable moiety (shown in bold) are labeled. See text for details.

In addition to the Hpy of the protein microenvironments, the Hpy in the solvent are needed to calculate the gain or loss of hydrophilicity in the solvent to protein transfer process. These were calculated by immersing each titratable residue in an equilibrated water droplet and carrying out 100 ps of dynamics after heating to 300 K. The residues were capped with neutral fragments at the termini. From the interval 70-100 ps, 31 structures were extracted, an Hpy calculated for each structure using the water fragment values given in Table 2, and the final value was obtained from the average of the Hpy values calculated from the 31 structures using Eq. 5. The results for the protonatable residues are given in Table 3 and the atoms belonging to the titratable groups, as defined in the PAR19 parameter set of CHARMM (Brooks et al., 1983), are also shown. It should be noted that in analogy to Rekker's (1977) formula for calculating log P from the hydrophobic fragmental constants, Hpy as defined in Eq. 5 is an extensive quantity. Therefore, in an environment that is uniformly hydrophobic or hydrophilic its magnitude will tend to increase with increasing size of the group. In addition, basic groups will tend to orient the waters with the dipole pointing away from the titration site (oxygen pointing toward the group), which will tend to increase the hydrophilicity of the microenvironment, whereas for the acidic groups the dipole will tend to orient toward the titration site (hydrogens nearer the acidic group) tending to increase the hydrophobicity.

Three measures for characterizing the hydrophobicity of the microenvironments are defined: first is Hpy_A defined in Eq. 5, second is the total hydrophobicity of the microenvironment due to the fraction of the group buried in the protein and the fraction exposed to solvent

(6)

where Hpy_A⁰ is the value of Hpy in the solvent, and finally, a measure of the hydrophobicity of a given group's microenvironment relative to what it would be if the group was totally immersed in water, i.e.,

(7)

These quantities have been calculated for all proteins in the data set, and average values, calculated from the seven proteins in the data set for each type of titratable residue, are given in Table 4. It is of interest that the buried fractions and microenvironments are quite different for the different residue types. The reciprocal of rHpy is an indication of the increase in the environmental hydrophobicity on transferring the titratable group from solvent to its final position in the protein. Thus, for the lysines, this increase is only about 1.5, which is not surprising considering the small value of . Histidine is the second most buried residue, but on the average, its microenvironment appears to be quite hydrophilic, and the increase in hydrophobicity is less than two. The titratable moiety of Tyr is the most buried with the most hydrophobic microenvironments. In an initial more extended study of the microenvironments of all groups in a sample of 13 globular proteins comprising a total of 7132 groups, it was found that and had values of 0.81 and 0.43, respectively (Mehler, unpublished results). Therefore, the average increase in hydrophobicity is about 2.3, which compares well with the value of 2.6 obtained in an earlier study by Ponnuswamy et al. (1980) using very different methods and assumptions.

The buried fractions (BF) and microenvironment properties of each titratable moiety in lysozyme are given in Table 5. It is seen that only a few values of Hpy are greater than zero. However, as will be discussed below, it was found that titratable groups with a BF less than 0.7 should be considered as solvent exposed for the pK_a calculations. Considering groups with BF > 0.7, only Glu-35 has a positive value of Hpy. The values of THpy or rHpy provide additional insight into the change in the hydrophobicity of the microenvironment when the group is transferred from water to protein. In HEWL, THpy is always negative and the highest value is 0.59 for Glu-35; the next largest value is 1.99 for Tyr-53. Comparison of the microenvironment of Glu-35 with the average values and standard deviations for Glu (Table 4) also clearly exhibits the unusual nature of the local environment around Glu-35.

Considering rHpy, the smaller the value the more hydrophobic the protein microenvironment of the group relative to water. Thus, the value of Glu-35 is the smallest (it implies that its protein microenvironment is around 20 times more hydrophobic than water!) in HEWL, and the next value is five times larger. Comparison of the microenvironments of Glu-35 and Asp-66 shows how the structural characteristics (see Fig. 2) translate into the quantitative description of the differences in hydrophobicity between the two microenvironments provided by the three hydrophobicity parameters. In HEWL, only Asp-52 has an Hpy value slightly more negative than Asp-66. Thus, the protein contribution to the microenvironments of these two residues are the most hydrophilic in HEWL, while that of Glu-35 is the most hydrophobic. It should finally be noted that the values of rHpy can be negative, indicating extreme hydrophobicity, or greater than one, indicating a microenvironment more hydrophilic than water.

Parametrization and the microenvironment

The hydrophobicity measures of the microenvironments in lysozyme (Fig. 2, Table 5) reveal a complex mosaic of different dielectric regions that present determinants of electrostatic properties that are quite independent of the buried fraction. In several important cases a totally buried residue, which usually is considered to exist in a low dielectric region, is shown to actually be found in a quite hydrophilic microenvironment. Thus, transfer of such a residue from water to protein is, in fact, only a small perturbation because the screening in such a microenvironment of the protein will be close to that in water. Nevertheless the microenvironments around the titratable groups of proteins show that they are generally more hydrophobic than in water (compare Tables 3, 4, and 5).

The values of D₀ and  used to define D₁ and D₂ were determined as described above to account for the expected differences in dielectric behavior depending on the extent of burial, but the explicit values were not calibrated against the pK_a of the protein data set. To incorporate these two screening functions into the algorithm, a switching value had to be determined, which was accomplished by assigning D₁ to all protonatable groups with a BF below a given threshold and D₂ for the others. A value of 0.7 was found to give the lowest root mean square deviation (rmsd) between calculated and experimental pK_a values in the data set and is used for switching between D₁ (BF < 0.7) and D₂ (BF  0.7) for all subsequent calculations. The calculations carried out in this way do not consider the explicit effects of the microenvironments and therefore are labeled as "menv-free" calculations in the following discussions.

From Fig. 2 and Tables 3, 4, and 5, it is apparent that for individual residues there can be considerable variation around the mean values of the properties characterizing the microenvironments. It seems reasonable, therefore, to make use of this to identify microenvironments that exhibit especially large deviations from the mean, and use the physical implications of these deviant regions to modulate the electrostatic properties of the titratable groups in question. To incorporate this idea into the algorithm, it is necessary to define quantitatively when the properties of a microenvironment deviate enough from the mean values to require modification of the menv-free assignments of the screening functions given above. These parameters were determined in the following way. The microenvironments of groups for which the calculated pK_a showed the largest errors were analyzed to rationalize on physical-chemical grounds the sources of the errors and the changes in screening that might lead to improvement. However, only reduction in overall rmsd was used to decide if a particular combination of screening functions was to replace the default assignments depending only on the BF value. This approach was taken because of the observation from the comparison of results using different crystal forms that a particular large error may be due to a conformation that is not relevant for the solution structure (see below) or that the sources of a particular error may be due to one or more inadequacies of the method. For these reasons, and because the results indicate that the quantification of the microenvironments will need further extension and refinement, exhaustive optimization was not carried out. Modification of the default screening was applied to titratable groups with BF  0.7, where the hydrophobicity or hydrophilicity of the microenvironments was far from the mean value; the border region with 0.7  BF  0.8 was treated differently from groups with BF > 0.8. For groups with very small buried fractions (BF < 0.3) and microenvironments very close to water, a third screening function was introduced (D₃ in Table 1 and Fig. 1) that is much closer to D_w than D₁ or D₂. Table 6 gives the thresholds for which the default conditions are modified for calculating the electrostatic energies, and they are discussed in the following paragraphs.

The protonatable moiety of Glu-35 in lysozyme is a buried titratable group in a very hydrophobic environment (see Fig. 2 and Table 5). According to Rekker (1977), Table IX, 1, p. 301), amino acid residues for which the fragmental hydrophobic constants sum to a relative hydrophobicity greater than one, are considered lipophilic, and this value has been used here to define an extremely hydrophobic microenvironment. As mentioned above, only a few buried tyrosines are in more hydrophobic microenvironments than Glu-35, but their pK_a values have not been measured. The effect of very hydrophobic environments is to reduce screening and to make the transfer of the titratable group from water to the protein more unfavorable. Since the pK_a has been measured only for Glu-35, any parameter adjustment is arbitrary, so that in these cases (see Table 6) w^tr has been scaled to increase the energy penalty for transferring a charge from solvent to the protein interior.

The average value of Hpy for the microenvironments surrounding the titratable residues in the seven proteins of the data set is about 1.5, but, as noted (see Table 4), the variations around this mean are considerable. For Asp and Glu when BF > 0.7, the microenvironments of a fairly large number of them are considerably more hydrophilic than 1.5. For such cases, the screening D₂ may imply a microenvironment that is too hydrophobic so that the large, unfavorable transfer energy results in pK_a values that are too high for the acids and too low for the bases. For 14 of the 53 aspartic and glutamic acids in the data set, this was the case, and the overestimation of these pK_a values was corrected by using D₁ instead of D₂ for calculating w^tr, because this implies transfer to a more hydrophilic microenvironment. For tyrosine and the bases the number of cases meeting the conditions of burial with BF greater than 0.7 and surrounded by very hydrophilic microenvironments was too small to assess any trends. It was found empirically, however, that using D₁ instead of D₂ for calculating w^int gave a small improvement in the overall rmsd. This result may be fortuitous because it appears to be counterintuitive. At the same time the response of any titratable group to its environment is determined by the balance between w^int and w^tr, so that, as the differences in properties given in Table 4 suggest, this response may depend on residue type. Finally, the third screening function, D₃, is only used when the properties of the microenvironment indicate conditions that are very close to pure solvent.

Computational details

The calculations reported below were carried out using protein coordinates obtained from the protein data bank (Bernstein et al., 1977) as follows: BPTI:4PTI* (Marquart et al., 1983), HEWL:2LZT* (Ramanadham et al., 1981) and 1HEL (Wilson et al., 1992), RnaseA:3RN3* (Howlin et al., 1989), RnaseT:3RNT* (Kostrewa et al., 1989), RnaseHI:2RN2* (Katayanagi et al., 1992) and 1RNH (Yang et al., 1990), Cab:3ICB* (Szebenyi and Moffat, 1986), protG:1PGA* and 1PGB (Gallagher et al., 1994), turkey ovomucoid inhibitor domain 3(OMT3):1PPF (Bode et al., 1986), and hiv protease (1HIV) (Thanki et al., 1992), where the starred coordinate sets were used in the optimization procedure. The preparation of the coordinates was carried out in the same way as in I, and the values of the model pK_a used here are N-term, 7.5; C-term, 3.8; His, 6.3; Glu, 4.4; Asp, 4.0; Tyr, 10.0; Lys, 10.4; and Arg, 12.0. The partial charges and group structures of the amino acid residues are those defined in the PAR19 topology file of CHARMM (Brooks et al., 1983), whereas the partial charges of the neutral forms of the titratable groups were taken from the PAR22 (MacKerell et al., 1992; MacKerell et al., 1995) topology files. The (BF)_a of Eq. A6 were evaluated from the solvent-accessible surface areas calculated with CHARMM using a probe radius of 1.0 Å, and the group average (defined by CHARMM) was used to scale all the w_a^tr for the given group (see Table 5 for the HEWL values). The Born radii were determined in the same way as described in the appendix of Mehler (1996b) except that atomic values were calculated, and the value obtained for D_w was used in both self energy terms in Eq. A6. For charges closer than 1.9 Å, w^int was scaled by a factor of 0.35. This scaling of nearby charges is similar to that used in Mehler (1996b) and compensates for the breakdown of the point charge model for nearby charges. However, due to the modified variational procedure (see Appendix) the cutoff distance could be reduced from 2.8 Å in Mehler (1996b) to the present value. Therefore, the scaling is essentially only used where uncompensated steric clashes leave two charges too close together. In addition, ionic strength is taken into account using a Debye screening as in Mehler (1996b). Finally, the polar protons of histidine are treated symmetrically in assigning the charge of the neutral species.

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION THEORETICAL BASIS AND... RESULTS AND DISCUSSION CONCLUSIONS APPENDIX REFERENCES

Calculation of pK_a

The calculated pK_a for the test suite of seven proteins are presented in Tables 7-13. To evaluate the effect on the pK_a of introducing the dependence of the electrostatics on the microenvironment, the pK_1/2 calculated with the screening functions D₁ and D₂ only, i.e., menv-free are also tabulated. Calculations were carried out in steps of 0.25 pH units and the pK_1/2 of each titratable group was estimated by linear interpolation between the two pH values where the fraction of charged species crossed the 0.5 value. The ionic strength was adjusted to be as close as possible to the experimental conditions and is 0.1 except where stated otherwise. For most of the proteins in the data set the experimental pK_a values to be used for calculating rmsd are clear from the quoted reference, but in some cases more than one choice is possible, and these will be noted as required for clarity. Finally, the Null model introduced by Antosiewicz et al., (1994) assumes that the pK_a values in the protein are the same as in the model (solvent), and it was noted that the rmsd for this model is usually quite small, since most of the pK_a values do not shift much from their model values. Because of this the rmsd resulting from the calculations provide only a partial measure of reliability. Of equal or greater importance are the number and size of the large errors, and these will also be considered in assessing the reliability of the approach.

Bovine pancreatic trypsin inhibitor

Calbindin D_9k

Hen egg lysozyme

View larger version (16K): [in this window] [in a new window]   FIGURE 3   Glu-35 of HEWL titration curve. Closed circles, menv-free calculation; open circles, includes microenvironment.

View larger version (19K): [in this window] [in a new window]   FIGURE 4   Energy components of Glu-35 as a function of pH. Closed symbols, transfer energy; open symbols, interaction energy; diamonds, menv-free calculation; circles, complete calculation.

Immunoglobulin G-binding domain B1 of Protein G

Ribonuclease A and Ribonuclease T1

Ribonuclease HI

Turkey ovomucoid third domain and HIV protease

Assessment of the results

Table 15 summarizes the reliability of the SCP approach using the description of microenvironments developed in this paper. Incorporation of the modulating effects of the titration site microenvironments on the electrostatic free energy leads to an improvement of about 15% in the rmsd relative to the menv-free calculation for the entire data set of 103 measured pK_a values. It should be realized, however, that with two different screening functions some accounting of environment, based on the buried fraction only, has already been included in the menv-free calculations. Therefore, the rmsd is smaller than that obtained by Antosiewics et al. (1994) from their FDPB calculations, but closer to the overall rmsd obtained by Demchuk and Wade (1996). More significant is the decrease in the larger errors due to inclusion of the environmental effects: at the error threshold of 0.5 pK unit, the decrease is 17% and increases to 33% at a threshold of 1.0 pK unit. A scatter plot of the pK_a values from both menv-free and microenvironment calculations is given in Fig. 5. The regression line has a slope of 0.98 and an intercept of 0.21. The correlation coefficient was 0.99. The scatter plot shows that accounting for the microenvironment induces shifts in the points so that they lie closer to the ideal line pK_calc = pK_exp.

View larger version (18K): [in this window] [in a new window]   FIGURE 5   Scatter plot of calculated versus experimental pKa values. Open triangles, menv-free results; closed circles, complete calculation.

An error analysis according to residue type is given in Table 16. With the present data set the largest errors were found in the pK_a values of Asp and Glu, with the bases having smaller errors; for Tyr and the termini not enough data is available to be statistically meaningful. The error trends from the present calculations are somewhat different from the trends noted in Antosiewicz et al. (1996), where the rmsd from all the residues were about the same.

The mean error of Asp shows that the calculated pK_a are too large. This trend is clearly seen in Fig. 6, where the experimental (Tanford and Roxby, 1972) and calculated titration curves of lysozyme have been plotted. The overestimation of the acidic pK_a values results in the slope of the calculated curve being too steep in the pH region 2-4. In the pH interval 4-6, the titration curve is primarily controlled by His-15 and Glu-35 that both titrate in this interval. The pK_1/2 points of the calculated titration curves of these two residues (see inset, Fig. 6) are shifted relative to the measured values (Bartik et al., 1994) such that their equilibrium charges approximately cancel in the pH interval 4.5-5.5 and lead to a flatter slope in this portion of the calculated titration curve. Inspection of the slope of the His-15 experimental titration curve suggests that it is steeper than the calculated slope (Bartik, K., C. Redfield, and C. M. Dobson, private communication). Finally, from a pH of 6 to 11, the experimental and calculated curves are in good agreement.

View larger version (21K): [in this window] [in a new window]   FIGURE 6   Calculated (open circles) and experimental (closed circles) titration curves of lysozyme. Inset shows calculated titration curves of His-15 and Glu-35.

Limitations of pK_a calculations using a single structure

In their analysis of pK_a values calculated with the FDPB algorithm, Antosiewicz et al. (1996) included average values obtained from NMR structures, and they noted a trend toward improvements in the results. The use of simulation to incorporate protein flexibility (relaxation effects) into the calculation of protein electrostatics has been addressed by Warshel and collaborators using linear response theory (Sham et al., 1997; Sham et al., 1998). Other approaches have been reported (Alexov and Gunner, 1997; Beroza and Case, 1996), and these authors also found some improvements in the calculated values. However, at present not enough systems have been tested to come to any general conclusions. Variability in the conformation of side chains was already pointed out by Bashford and Karplus (1990) and others (Yang et al., 1993; You and Bashford, 1995) in pK calculations using the trigonal and tetragonal crystal forms of HEWL. These differences are already reflected in the buried fractions of the titratable residues and in their microenvironments, which can show variations that, in some cases, are large enough to change the screening assignments for the calculation of the electrostatic free energy components.

Oda et al. (1993)