INTRODUCTION

TOP ABSTRACT INTRODUCTION THEORETICAL METHODS RESULTS AND DISCUSSION REFERENCES

Electrostatic interactions play important roles in the stability of proteins, as illustrated by protein unfolding at extreme pH values, yet quantitative understanding of these roles has proven elusive due to a number of factors such as the strength and long-range nature of these interactions, strong mediation by solvent, and interference of nonelectrostatic effects. Over 20 years ago Perutz (Perutz and Raidt, 1975; Perutz, 1978) noted that salt bridges could potentially increase folding stability. Experimental studies on charge mutations have led to inconsistent conclusions on the contributions of charge-charge interactions to protein stability, and efforts to introduce stabilizing salt bridges have met with mixed success (Anderson et al., 1990; Dao-pin et al., 1991; Sali et al., 1991; Marqusee and Sauer, 1994; Waldburger et al., 1995; Meeker et al., 1996; Tissot et al., 1996; Spek et al., 1998; Vetriani et al., 1998; Huang et al., 1998; Ogasahara et al., 1998; Grimsley et al., 1999; Merz et al., 1999; Ramos et al., 1999; Giletto and Pace, 1999; Loladze et al., 1999; Perl et al., 2000; Pace, 2000; Spector et al., 2000; Strop and Mayo, 2000; Takano et al., 2000; Burkhard et al., 2000; Shaw et al., 2001; Perl and Schmid, 2001; Sanchez-Ruiz and Makhatadze, 2001; Olson et al., 2001; Kammerer et al., 2001; Loladze and Makhatadze, 2002). In contrast, in a number of theoretical studies based on continuum electrostatics, the view appears to have emerged that overall electrostatic interactions destabilize or marginally stabilize proteins and protein complexes (Novotny and Sharp, 1992; Hendsch and Tidor, 1994; Elcock, 1998; Elcock et al., 1999; Sheinerman et al., 2000; Lee and Tidor, 2001). The main argument is that the desolvation cost for bringing two charges together upon protein folding or complex formation is so large that it may more than offset the energetic contribution of the charge-charge interaction. We have recently noted that the desolvation cost calculated by continuum electrostatics is very sensitive to the definition of the boundary between the high solvent dielectric and the low protein dielectric (Vijayakumar and Zhou, 2001). The large desolvation cost calculated by others using the solvent-exclusion (SE) surface as the dielectric boundary is reduced substantially when we used the van der Waals (vdW) surface. A priori, neither surface is preferred and the choice must be resolved by testing against experiment. We found that the vdW surface gave much better agreement between calculated and experimental effects of 12 charge mutations on the folding stability of barnase.

In this paper we present Poisson-Boltzmann (PB) calculations for the effects of charge reversal at five solvent-exposed sites (K16E, R119E, K135E, K147E, and R154E) and charge neutralization and proton titration of the H31-D70 semi-buried salt bridge on the stability of T4 lysozyme. We adopt the protocol of using vdW surface as the dielectric boundary, as established in our earlier study on charge mutations in barnase. We also explicitly account for residual charge-charge interactions in the unfolded state. These residual interactions have been shown to be important for accounting for the pH dependence of the unfolding free energy (Elcock, 1999; Pace et al., 2000; Zhou, 2002).

T4 lysozyme has a net charge of +9e at neutral pH (a total of 47 ionizable groups with 10 Asp, 8 Glu, 13 Lys, 13 Arg, 1 His, and the N- and C-terminals). This large net positive charge might signal significant repulsion between like charges and reversing some of the positive charges may stabilize the protein. This was the motivation of Dao-pin et al. (1991) for studying the effects of charge reversal at the five solvent-exposed sites. The electrostatic contributions of the five charge reversal mutations to G_unfold from our PB calculations range from 1.6 to 1.3 kcal/mol. This indicates that the variable effects of charge reversal observed by Dao-pin et al. are not unexpected. The variability of the effects arises because a charge reversal may stabilize the unfolded state not as much as (as in the cases of K16E and R119E), as much as (in the case of K135E), or even more than (in the case of K147E) the folded state, or it may destabilize the folded state while stabilizing the unfolded state (in the case of R154E).

The H31N, D70N, and H31N/D70N mutations were found to destabilize the protein by 2.9, 1.3, and 1.6 kcal/mol, and the pK_a values of H31 and D70 were shifted to 9.4 and 0.6, respectively. These results are in good accord with experimental data of Anderson et al. (1990). In contrast, if the SE surface were used, the H31N/D70N mutant would be found to more stable than the wild-type protein by 1.3 kcal/mol. These and additional calculation results for 27 charge mutational on five other proteins lead us to conclude that, while solvent-exposed charges may not reliably contribute to protein stability, semi-buried salt bridges can provide significant stabilization.

	THEORETICAL METHODS

TOP ABSTRACT INTRODUCTION THEORETICAL METHODS RESULTS AND DISCUSSION REFERENCES

We calculated the electrostatic contributions of charge mutations to the folding stability of T4 lysozyme. The calculated results for G_el, the electrostatic component of the change G_unfold in unfolding free energy, were compared with experimental results for G_unfold. G_el may be conveniently viewed to be composed of two terms:

(1)

For the first term G, the unfolded state is assumed to be devoid of any charge-charge interactions. Thus the unfolded state is simply modeled as the residue under mutation alone exposed to the solvent. The second term accounts for residual charge-charge interactions in the unfolded state. G was obtained from the dielectric continuum model for the unfolded protein and for the isolated mutation residue (a primitive representation of the unfolded state) (Vijayakumar and Zhou, 2001). G was obtained from the Gaussian-chain model (Zhou, 2002).

Generation of mutant structures

The locations of the five exposed charged residues and the H31-D70 salt bridge are shown in Fig. 1. For the purpose of calculating G, it is important that mutant structures have minimal differences from the wild-type structure. That is, structural changes are isolated to the mutated side chains alone. Otherwise, differences in other parts of the protein will overwhelm the mutated residue in contributing to G. Our standard procedure for preparing protein structures for continuum electrostatic calculations are as follows (Vijayakumar and Zhou, 2001). First, hydrogens were added to the x-ray structure of wild-type T4 lysozyme (PDB entry 3lzm; Weaver and Matthews, 1987) by using the program InsightII (Molecular Simulations, Inc.). Mutations were generated in InsightII and optimized by energy minimizing side-chain atoms beyond C. The AMBER force field (Weiner et al., 1984) was used for the minimization.

View larger version (52K): [in this window] [in a new window]   FIGURE 1   Locations of five solvent-exposed charges and the H31-D70 salt bridge in T4 lysozyme.

For the five solvent-exposed charged residues, mutant structures have been determined by x-ray crystallography (Dao-pin et al., 1991). These mutant structures provide opportunities for modifying the standard procedure. For the R119E, K135E, and K147E mutants, structural changes are limited to just the mutated side chains. The conformations of these three mutated side chains were transplanted to the wild-type protein as the starting conformations for optimization. For the K16E and R154E mutants, there are concurrent changes in side chains not mutated. For these two mutants, we simply reverted to the standard procedure. The H31N and D70N mutants were generated by the standard procedure, whereas the H31N/D70N double mutant was generated by making a single D70N mutation on the H31N mutant.

Calculation of G

The calculation G followed the same protocol as established in our earlier work on barnase charge mutations (Vijayakumar and Zhou, 2001). Briefly, the PB equation was solved by the UHBD program (Madura et al., 1995), with the vdW surface used by deselecting the "nmap 1.5, nsph 500" option. The electrostatic potential  was calculated first from a 100 × 100 × 100 grid with 1.5 Å spacing centered at the geometric center of the wild-type protein. This was followed by a 140 × 140 × 140 grid with 0.5 Å spacing at the same center. A final round of focusing at the CB atom of a mutated side chain was introduced on a 60 × 60 × 60 grid with 0.25 Å spacing. The electrostatic energies of the folded protein and the isolated mutation residue were calculated by

(2)

where q_i are the partial charges. G was obtained by taking the difference in G_el between the unfolded state (as represented by the isolated mutation residue) and the folded state, and then the difference of these differences between the mutant and the wild type.

As in our earlier work, Amber charges and radii were used. Asp, Glu, and the C-terminal were unprotonated, whereas the lone His [H31, with a pK_a measured at 9.1 (Anderson et al., 1990)], the N-terminal, Lys, and Arg were protonated. For the five exposed charged residues, the ionic strength was 50 mM and the temperature was 65°C. For mutations on the H31-D70 salt bridge, the ionic strength was 100 mM and the temperature was 10°C. The protein dielectric constant was 4 and the solvent dielectric constant was set be that of water at the particular temperature.

Calculation of interaction energy in the folded state

G can be decomposed into three terms (Vijayakumar and Zhou, 2001):

(3)

where G_solv is the difference between the mutant and the wild-type protein in the changes in desolvation cost upon unfolding, G'_el is the difference in the electrostatic energies when the mutated residue is completely discharged, and G is the difference in the interaction energies of the mutated residue with the rest of the protein. The first term can be calculated by discharging the whole protein except for the mutated residue, while the second term can be obtained by discharging the mutated residue. The third term can be obtained by calculating the electrostatic potential of the rest of the protein and then multiplying the partial charges of the mutated residue (Eq. 2 without the factor of 1/2). Alternatively, it can also be obtained by multiplying the potential of the mutated residue with the partial charges of the rest of the protein. In the latter procedure, one may also obtain information about the interaction energy between the mutated residue and each partial charge of the rest of the protein.

Calculation of pK_a values

The pK_a values of ionizable groups in a protein are determined by the electrostatic energies of protonation. In general, protonations of different ionizable groups are coupled. However, the proton titration of H31 (with a measured pK_a of 9.1) occurs in a pH range in which all acidic groups are deprotonated and the other basic groups are still protonated. This suggests that the pK_a of H31 can be calculated from the electrostatic energy of protonation while fixing the protonation states of all other ionizable groups. This approach was used previously (Vijayakumar and Zhou, 2001) to accurately calculate the pK_a of D93 of barnase. The ionizable group was assigned appropriate partial charges in both the deprotonated and protonated states.

If the protonation of an ionizable group is viewed as a mutation, its pK_a can be calculated as

(4)

where pK₀ is the pK_a of a model compound, and Z₀ and Z₁ represent the unprotonated and protonated forms, respectively, of the ionizable group. An implicit assumption in Eq. 4 is that the isolated ionizable group takes the model compound pK_a (i.e., pK₀). The pK_a values of H31 in the wild-type protein and the D70N mutant, and D70 in the wild-type protein and the H31N mutant, were calculated using Eq. 4. In calculating the pK_a values of D70, the protonation states of all ionizable groups were again fixed. Because the assignment of protonation states of the other groups requires an initial guess for the pK_a of the group under investigation, the calculation of this pK_a may be viewed as a confirmation of the initial guess. The model-compound pK_a values were 4.0 for Asp and 6.3 for His.

Calculation of G

In the unfolded state, the Gaussian-chain model predicts the interaction energy of the mth ionizable group by Zhou (2002)

(5)

where x_i are the protonation states and x_i0 are their values when the groups are neutral, and W_mi are the strengths of interaction as specified by the Gaussian-chain model. G is the difference in G between the mutant and the wild-type protein.

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION THEORETICAL METHODS RESULTS AND DISCUSSION REFERENCES

Electrostatic contributions of reversing solvent-exposed charges

As shown in Table 1, the electrostatic contributions of the five charge reversal mutations to G_unfold were found to range from 1.6 to 1.3 kcal/mol. The variability of the effects arises because a charge reversal may stabilize the unfolded state not as much as (as in the cases of K16E and R119E), as much as (in the case of K135E), or even more than (in the case of K147E) the folded state, or it may destabilize the folded state while stabilizing the unfolded state (in the case of R154E). Such variability is consistent with the experimental observations of Dao-pin et al. (1991).

Quantitative agreement with the experimental results for G_unfold are reasonable for K16E, K147E, and R154E. For R119E, G_el indicates that the charge reversal should increase the folding stability by 1.3 kcal/mol, but experimentally R119E was found to have a marginal effect on the stability. For K135E, G_el indicates that the charge reversal should have a marginal effect; but experimentally, K135E was found to decrease the stability by 1 kcal/mol. A perfect match between calculated G_el and experimental G_unfold is not expected. The measured G_unfold will likely have nonelectrostatic contributions (arising, e.g., from nonpolar interactions and effects of side-chain entropy). The continuum electrostatic model obviously has its own limitations. For all five mutations, inclusion of residual charge-charge interactions in the unfolded state brings calculated G_el into much closer agreement with experimental G_unfold.

The destabilization of the folded state by R154E is due to the interaction with D127. In the wild-type protein, both the NE and NH1 atoms of R154 are 5.1 Å away from the OD2 of D127. This favorable interaction is changed to a repulsive one upon the R154E mutation.

Electrostatic contribution of the H31-D70 salt bridge

As shown in Table 1, the H31N, D70N, and H31N/D70N mutations were found to destabilize the protein by 2.9, 1.3, and 1.6 kcal/mol. These results are in good accord with experimental data of Anderson et al. (1990), who observed a decrease of 7°C in the melting temperature at pH 7 for all the three mutants. This decrease in melting temperature corresponds to a decrease in G_unfold of ~2.5 kcal/mol. The coupling energy for this salt bridge,

(6)

is 2.63 kcal/mol, agreeing with the experimental result of 2.5 kcal/mol. In general, the coupling energy is less corrupted by nonelectrostatic effects, thus comparison between electrostatic calculation and experiment in this case is the fairest.

The 2.9 kcal/mol value of G by the mutation H31N consists of 0.9 kcal/mol in G_solv, 0.5 kcal/mol in G'_el, and 4.3 in G. The interaction energy almost exclusively comes from the interaction of residue 31 with D70. The 1.3 kcal/mol value of G by the mutation D70N consists of 1.0 kcal/mol in G_solv, 0.1 kcal/mol in G'_el, and 2.4 in G. This interaction energy almost exclusively comes from the interaction of residue 70 with H31. The nearly twofold difference in G between the two mutants comes about because Asn is a much closer substitute for Asp than for His.

In contrast to the results on residual electrostatic effects found for the five exposed charges, such interactions in the unfolded state are negligible for H31 and D70. This difference has to do with the distribution of charges along the sequence. In the Gaussian-chain model for the unfolded state, charge-charge interactions are dominated by charges close to the sequence. For example, the strong stabilization of the unfolded state by K147E reflects the fact that the five nearest ionizable groups (K135, R137, R145, R148, and R154) along the sequence are all positively charged. However, for both H31 and D70, charges on the two sides along the sequence have opposite signs and interactions with them are nearly canceled (D20 and E22 versus K35 and K43 for the former and D61, E62, E64, K65 versus D72, R76, R80, and K83 for the latter). Of course, because of the large sequence separation between H31 and D70, the coupling between them in the unfolded state is expected to be weak and indeed found to be negligible.

The pK_a values of H31 and D70 in the wild-type protein were calculated to be 9.4 and 0.6, respectively. These are in excellent agreement with the measured values of 9.1 and 0.5 (Anderson et al., 1990, 1993). Upon the D70N mutation, the pK_a of H31 was found to decrease to 6.4, close to the experimental value of 6.9 (Anderson et al., 1990). The pK_a shift of H31 to a normal value upon the D70N mutation validates the earlier result that electrostatic interactions with H31 are dominated by the D70 residue. Upon the H31N mutation, the pK_a of D70 was found to increase to 3.1, again consistent with experiment (Anderson et al., 1990).

Critical importance of the choice of dielectric boundary

The results reported above were obtained by using the vdW surface as the dielectric boundary, a protocol that we have previously found to give the best predictions for the effects of charge mutations in barnase (Vijayakumar and Zhou, 2001). Here again we found that the common practice of using the SE surface yields unsatisfactory results. Using the SE surface, G for the H31N, D70N, and H31N/D70N mutations were found to be 4.1, 0.1, and 1.3 kcal/mol, respectively. The prediction that the H31N/D70N mutant is more stable than the wild-type protein by 1.3 kcal/mol is in stark contrast to the experimental finding and the vdW-surface calculation result.

The main reason for the wrongly predicted higher stability of the H31N/D70N mutant by using the SE surface is the excessively high cost for desolvating H31 and D70 (7.5 versus 1.7 kcal/mol predicted by using the vdW surface). Meanwhile, the predicted coupling energy, 5.3 kcal/mol, is twice as large as the experimental and vdW-surface calculation results. These findings confirm the shortcomings of using the SE surface as noted in our previous study of charge mutations in barnase (Vijayakumar and Zhou, 2001).

The difference between the vdW and SE surfaces consists of crevices not accessible to a spherical probe (with a radius of 1.4 Å). For a semi-buried charged residue, the crevices around neighboring residues add up and significantly change the accessibility of the charge. For six charged side chains in barnase, we have shown that, on average, 90% of the vdW surfaces are exposed, but only 30% of the SE surfaces are exposed (Vijayakumar and Zhou, 2001). This significantly increased exposure of the vdW surface accounts for the decreased desolvation cost. The vdW surface can leave small holes in the protein interior, but we did not find the presence of these holes to have any consequence on calculation results. For wild-type T4 lysozyme, three small holes were found. We filled these holes with dummy atoms (with radii of 0.3-0.5 Å). The solvation energy of the protein was unchanged.

As might be expected, however, the choice of the dielectric boundary is much less important for the solvent-exposed charged residues. Using the SE surface, G was found to be 1.31, 3.07, 0.47, 0.03, and 1.24 kcal/mol, respectively, for K16E, R119E, K135E, K147E, and R154E. These (except for R119E) are similar to the results, listed in Table 1, obtained using the vdW surface.

Application of calculation protocol to other proteins

We appear to have established an electrostatic calculation protocol that reasonably predicts effects of charge-charge interactions in proteins. As illustrations of the robustness of the protocol, we have studied 27 additional charge mutations on five other proteins. These include neutralizations of six semi-buried aspartates and glutamates involved in salt-bridge and hydrogen-bonding interactions in human lysozyme, charge reversals of five solvent-exposed aspartates and glutamates in ribonuclease Sa (net charge on the wild-type protein is 7e), mutations that eliminate three of the 12 differences between the sequences of Bacillus subtilis cold shock protein B and the thermophilic Bacillus caldolyticus cold shock protein, and alanine substitutions of three residues (D14, R17, and S77) forming a semi-buried salt-bridge/hydrogen-bonding network in  repressor.

Comparison of calculated and experimental results on G for these proteins is presented in Table 1. Overall, the agreement is reasonable. However, there are a number of overestimates of the effects of charge mutations (D76N on human lysozyme, D25K and E74K on ribonuclease Sa, R17A and R17A/S77A on  repressor). These can be partly attributed to the fact that, in the present protocol, residues around the mutation site are not allowed to relax. In particular, optimizations of residues around D76N in human lysozyme and R17A in  repressor might lower the effects of the mutations; however, accurate molecular modeling of these optimizations is difficult.

The calculation results for the charge reversals of the five solvent-exposed residues in ribonuclease Sa again indicate that the effects of such charged residues are variable. The net charge of 7e on the wild-type protein should provide a generally favorable environment for the positive charge resulting from a D or E to K mutation, yet E41K is destabilizing because of the higher desolvation cost of E41K and the loss of the favorable interaction between E41 and R40 (OE2 to NE distance at 5.8 Å). We do not have an explanation for the experimentally observed destabilizing effect of the D17K mutation (Shaw et al., 2001).

The calculation results for the five mutations on Bacillus subtilis cold shock protein B are worth noting. Most of the other mutations studied are destabilizing, whereas four of the five mutations on this protein are stabilizing, in total agreement with experiment (Perl and Schmid, 2001). Experimentally, the E3R/E66L double mutation was found to contribute 3.4 kcal/mol toward the 3.8 kcal/mol difference in stability between the mesophilic and thermophilic proteins. The calculated G for the E3R/E66L mutant was also 3.4 kcal/mol. Main contributions to this G are 1.8 kcal/mol from the elimination of the desolvation cost for E66 (E3 and R3 have the desolvation cost), 1.2 kcal/mol from better interactions of R3 than E3 with the rest of the protein (the net charge on the wild-type protein is 6e), and 0.3 kcal/mol from the elimination of the electrostatic repulsion between E3 and E66.

When the SE surface was used instead, agreement with experiment deteriorated significantly. Specifically, relative to the vdW-surface results, the magnitudes of calculated effects of charge mutations were further increased by 0.7, 1.0, 0.7, 4.8, 1.9, 2.1, and 2.2 kcal/mol for the human lysozyme D67N, ribonuclease Sa D17K, D25K, E41K, and E74K, and Bacillus subtilis cold shock protein B E66L and E3R/E66L mutations, respectively. Moreover, the  repressor D14A, D14A/R17A, D14A/A77, and D14A/R17A/S7A mutations were incorrectly predicted to stabilize the protein. Again, the fault primarily lies in overestimated desolvation cost. For example, the changes in desolvation cost upon the ribonuclease Sa E41K, Bacillus subtilis cold shock protein B E66L, and  repressor D14A mutations were 1.0, 1.9, and 2.0 kcal/mol, respectively, according to the vdW surface. These became 4.1, 4.6, and 7.3 kcal/mol, respectively, according to the SE surface.

We also studied mutations of Asp-76 in ribonuclease T₁. This residue is buried and forms four hydrogen bonds with three polar side chains and a buried, conserved water molecule. Giletto and Pace (1999) investigated the contributions of Asp-76 to the stability of ribonuclease T₁ by mutating it to Asn, Ser, and Ala. Both urea and thermal unfolding showed that the mutants were less stable by ~3.5 kcal/mol. Our calculation found G_el to be 4.0, 3.6, and 3.9 kcal/mol for D76N, D76S, and D76A, respectively, in close agreement with experiment. The pK_a of Asp-76 was predicted to 1.1, also in good agreement with the experimental value of 0.5 (Giletto and Pace, 1999). One of the three polar side chains hydrogen-bonded to Asp-76 is from Tyr-11. Model building suggests that a mutation of Tyr-11 to Arg may introduce a salt bridge with Asp-76. The Y11R mutation is predicted to stabilize ribonuclease T₁ by 2.4 kcal/mol. Whether this mutation indeed stabilizes the protein awaits experimental test.

As another application, we calculated the pK_a of the single histidine, His-68, in ubiquitin. The pK_a was found to be downshifted to 5.0, mainly because of an unfavorable interaction between His-68 while protonated and Lys-6. The calculated result is in reasonable agreement with the experimental value of 5.9 measured by NMR (Ibarra-Molero et al., 1999).

Contrasting roles of semi-buried salt bridges and exposed charged residues

Both of our earlier studies of semi-buried salt bridges in barnase and the present study of the H31-D70 salt bridge in T4 lysozyme and the D14-R17 salt bridge in  repressor show that they make significant contributions to the folding stability. However, we now have shown that the effects of exposed charges exhibit variability. Recognizing the different roles of semi-buried salt bridges and exposed charges is very important. This recognition may help reconcile some of the conflicting reports regarding the contributions of electrostatic interactions to protein stability.

In conclusion, that the use of the SE surface as the dielectric boundary may lead to overestimated desolvation cost suggests a reexamination of the popular view that electrostatic interactions are generally destabilizing. While solvent-exposed charges may not reliably contribute to protein stability, semi-buried salt bridges can provide significant stabilization.