Modeling Charge Interactions and Redox Properties in DsbA*

Accurate prediction of charge interactions in macro-molecules presents a significant challenge for computational biology. A model for the low Cys 30 p K a and oxidizing power of DsbA

DsbA from Escherichia coli is a periplasmic protein of 189 amino acids that catalyzes the introduction of disulfide bonds (1,2). The oxidizing power of DsbA derives from the low Cys 30 pK a of about 3.5 (3,4), exhibiting a ⌬pK a with respect to normal of around Ϫ5. DsbA is a valuable link to theory because models must account for the low Cys 30 pK a and thereby provide a molecular basis for oxidizing power and physiological function. Continuum electrostatics has become the most common method for calculating macromolecular pK a s (5), but it is not yet consistently accurate, with discussion revolving around the choice of macromolecular relative dielectric, ⑀ p (6 -8). Variation in ⑀ p reflects the difficulty of reproducing microscopic solvation effects in a continuum model, particularly where an ionizable group is buried in the macromolecule, so that part of the high relative dielectric (⑀ s ) solvation shell is swapped for the lower ⑀ p environment (see Fig. 1).
A continuum electrostatic model has been presented for the redox potential difference between E. coli DsbA and E. coli thioredoxin (9), using structural homology and differencing to circumvent changes in charge burial (see Fig. 1). The low thiolate pK a and oxidizing power of DsbA were suggested to arise from several sources including His 32 and Gln 97 side chains and summed peptide dipole contributions. Residues Glu 37 and Glu 38 were each predicted, if deprotonated, to move the Cys 30 pK a by about 0.4 in the more reducing direction. The P34H mutation in thioredoxin supports the model implication for DsbA His 32 (10). Redox potential and Cys 30 pK a measurements have been reported for DsbA mutations Glu 37 and Glu 38 and for the deletion mutant ⌬38 -40, which relates to assessment of the peptide dipoles contribution (11). In contrast to the reported disparity with the predictions (11), the current article will demonstrate that data for ⌬38 -40 are in line with the model. The deleted region is just one part of an implicated section of backbone, and its individual effect, although visible, is not large. The effects of Glu 37 and Glu 38 depend on their protonation states, and it will be shown that mutant measurements (11) and the atomic structure (12) are consistent with an elevated pK a for one of these residues.
Whereas lower ⑀ p gives a reasonable match to experiment in the Cys 30 thiolate difference calculations, it generates large discrepancies in unmodified full pK a calculations. A recently introduced modification (8), empirically accounting for changes in water entropy upon charge burial, is shown to perform qualitatively well when applied across the ionizable groups of DsbA.

EXPERIMENTAL PROCEDURES
Reduced DsbA has been modeled previously (9) from the crystallographic structure of the oxidized molecule at 2-Å resolution (12) by breaking the Cys 30 -Cys 33 disulfide bond and torsioning Cys 30 to maintain van der Waals contact with Cys 33 . The modeled reduced configuration is similar to that of the homologous Cys 32 in the nmr structures of reduced E. coli and human thioredoxins (13,14). In the absence of experimentally determined atomic structures for the DsbA mutants, E37Q and E38Q were assumed to be isosteric with WT, 1 and the deletion mutant ⌬38 -40 could be made with a C ␣ -C ␣ link between WT residues 37 and 41, accompanied by only minor conformational rearrangement upon regularization of this region in the program QUANTA, with the CHARMm force field (Molecular Simulations Inc., Waltham, MA). Neither the WT nor ⌬38 -40 proteins were subjected to extensive energy minimization, so that the structures would remain close to experiment and also match each other away from the deletion site. Calculated differences between WT and the mutants are based on the assumption of minimal structural alteration, which is consistent with stereochemical and hydrogen bonding considerations. The efficacy of such conformational modeling can be assessed when the mutant and reduced WT structures are determined experimentally.
Charge interactions were calculated with FD solutions to the Poisson-Boltzmann equation (15,16), implemented in the program FD-CALC, using ⑀ s ϭ 80 and ⑀ p ranging from 4 to 80. Calculations of pK a differences between DsbA WT and mutants (see Fig. 1) used reported ionizable charge assignment (9) for groups other than Cys 30 . In addition, calculations of the electrostatic free energy, differenced between reduced and oxidized forms and between WT and ⌬38 -40 mutant, ⌬⌬⌬G ϭ (⌬G WT,RED Ϫ ⌬G WT,OX ) Ϫ (⌬G MUT,RED Ϫ ⌬G MUT,OX ), are compared with ratios of redox equilibrium measurements (11), ϪRT ln(K eq WT / K eq MUT ), from 0 to 1 M added NaCl, with T ϭ 298.15 K. It is assumed that the structural changes associated with both disulfide bond (Cys 30 -Cys 33 ) breakage and mutation are confined to separate localized regions and that these localized effects will cancel between WT and mutant (for disulfide bond breakage) and between oxidized and reduced (for the mutation). The calculated ⌬⌬⌬G value is therefore the remaining long range (electrostatic) interaction between the mutation and the Cys 30 thiolate. Ratios of K eq at each added NaCl concentration will remove any ionic strength dependence of the glutathione redox potential. Ionic * This work was supported by Biotechnology and Biological Sciences Research Council of the United Kingdom. The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked "advertisement" in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.
strength dependence was incorporated into the FD computations without a Stern layer. Full pK a calculations ( Fig. 1) used a statistical treatment of interacting ionizable groups (17), extended with a Monte Carlo method for computations with large numbers of such groups (18). This method used 10,000 Monte Carlo steps and a modification that allows for multiple site transitions for pairs that are are coupled by an interaction equivalent to more than 2 pK a units (18). Partial charges (24) and ionizable group (free amino acid) pK a s (6) were allocated. Ionizable residues included were Asp, Glu, His, Arg, Lys, Cys, and the amino-terminal group, whereas the carboxyl-terminal residue (189) is missing from the coordinates, and tyrosines have been omitted from the pK a calculations. Modification to account empirically for solvent entropic change upon amino acid transfer to protein is ⌬E t ϭ ⌬V s E s , where ⌬V s is the fractional change in first hydration shell volume (calculated from the FD grids) and E s is a free energy contribution associated with water ordering for a complete first hydration shell (8). Fitting to a range of experimental pK a s, with ⑀ p ϭ 4, gives values of E s that correspond to about 6 ordered water molecules in the hydration shell of a single charge center group and about 2 for a double charge center group (8). Although there may be further detailed variation between ionizable group types within the charge center groupings, these values are starting points for overall estimates of pH titration curves. Some calculations were made using a Debye-Hü ckel model with a uniform dielectric, ⑀ eff ϭ 50, in place of the FD procedure. The higher dielectric and the neglect of counterion exclusion from the protein interior in this method (19) reduces the size of electrostatic interactions.

RESULTS AND DISCUSSION
Difference Calculations between WT and Mutants E37Q, E38Q, and E37Q/E38Q-The pK a s of these residues in WT DsbA are currently unknown, and full pK a calculations are not sufficiently reliable to provide detailed estimates, largely due to the charge burial term. The more reliable charge-charge estimates yield interactions of about 2 kJ/mol between a deprotonated Glu 37 or Glu 38 and the Cys 30 thiolate. In the WT structure at pH 6.5 (12), Glu 37 and Glu 38 carboxylates approach within 3 Å (Fig. 2), strongly suggesting that they share a proton and that one of the pK a s will be elevated above neutral pH, so that just one carboxylate to thiolate interaction should be counted for comparison to redox equilibrium measurements at pH 7. The calculated WT to double mutant difference would be this single interaction, approximating no conformational change and limited glutamine side chain charge effects. With regard to the single mutants, Glu 38 lies toward the protein exterior and Glu 37 facing the protein interior (Fig. 2), so that Glu 37 is likely to be buried within the E38Q mutant whereas Glu 38 will be solvent-exposed in the E37Q mutant. These considerations would be consistent with deprotonation of Glu 38 in the E37Q mutant but neutralization of Glu 37 in the E38Q mutant. The calculated ⌬⌬Gs, WT to (E37Q, E38Q, E37Q/ E38Q) mutants, for this hypothesis with the ⑀ p ϭ 4 model are (0, Ϫ2, Ϫ2) kJ/mol compared with measured values of (0.6, Ϫ2.0, Ϫ1.5) in 10 mM sodium phosphate (11). This reasonable agreement, constructed upon a plausible hypothesis for Glu 37 and Glu 38 pK a s, demonstrates both the requirement for more accurate full pK a calculations and the success of the published model (9) in suggesting a route toward engineering more oxidizing DsbA molecules at neutral pH. Comparisons at the acidic pH of the Cys 30 pK a are omitted for this set of mutants because it is likely that both Glu 37 and Glu 38 will be protonated in this pH region.
Difference Calculations and Salt Dependence for Redox Equilibria of WT versus ⌬38 -40 Mutant-The deleted residues 38 -40 are within a larger polypeptide region (Fig. 2), which in total is predicted to contribute to thiolate stabilization in the low ⑀ p model. It is possible to make the Glu 37 -His 41 link without substantial disruption to the rest of the protein. Following the discussion in the previous section, Glu 37 becomes solventexposed in the modeled deletion mutant so that, by analogy with WT Glu 38 , it is likely to be deprotonated at the neutral pH of K eq measurements. In difference calculations between WT and ⌬38 -40 it is assumed that the thiolate interactions to the FIG. 1. Full pK a calculations and charge burial within a mac-romolecule. The pK a of a schematically drawn cysteine in the WT protein (pK a WT , top left) is derived from the pK a of the free amino acid (pK a AA , top right) and the electrostatic energy difference between ionization in WT and free amino acid (⌬⌬G WT:AA ϭ ⌬G WT Ϫ ⌬G AA ), pK a WT ϭ pK a AA ϩ ⌬pK a WT:AA , with ⌬pK a WT:AA ϭ (1/2.303RT)⌬⌬G WT:AA , where R is the universal gas constant and T the absolute temperature. The terms ⌬G WT and ⌬G AA contain components from the Born (self) energy (22) and from charge-charge interactions (17). It is assumed that electrostatic contributions associated with neutral cysteine are negligible, focusing on the thiolate interactions in protein and free amino acid. Cysteine pK a in a mutant (pK a MUT ) is also derived from differencing with the free amino acid (lower half of the figure). Significant errors in full pK a calculations may be associated with modeling charge burial (8), denoted by first hydration (hyd) shell occlusion between free amino acid and protein in this figure, because the Born energy is highly dependent on the difference between ⑀ p and ⑀ s . Models with low ⑀ p tend to overestimate the cost of charge burial (6). The figure indicates how differencing between WT and mutant (or related) proteins with structurally similar thiolate environments, ⌬pK a WT:MUT ϭ pK a WT Ϫ pK a MUT , circumvents the hydration shell changes that arise from comparison with the free amino acid (9,23). Suggested mechanisms for full pK a calculations that reduce the cost of charge burial include the use of relatively high ⑀ p (6) and the use of lower ⑀ p together with an empirical estimate of the favorable entropic contribution associated with water liberation from the first hydration shell (8). Other ionizable groups are set at normal neutral pH values except for His 32 and Glu 24 , which are in the vicinity of the active site and may have pK a s around neutral pH (9). These ionizations were set to ϩ0.5 and Ϫ0.5, respectively. For the WT versus ⌬38 -40 mutant calculations, both His 32 and Glu 24 contributions almost cancel on differencing.
Calculated Cys 30 thiolate contributions to ⌬⌬⌬G are compared with experimental values derived from K eq ratios (11) for WT versus ⌬38 -40 over a range of added NaCl concentrations ( Table I). The approximations in the calculations (modeled ⌬38 -40 conformation, ionizable charge assignment, and K eq contribution from Cys 30 thiolate alone included) combined with the small ⌬⌬⌬G values (thermal energy or less) suggest that qualitative rather than quantitative comparisons should be made. It can be seen that the ⑀ p ϭ 4 model is by far the closest to experiment, indicating that charge-charge interactions through a low protein dielectric are important in DsbA and that higher ⑀ p values underestimate these interactions. Within the ⑀ p ϭ 4 model, the listing of total and non-ionizable interactions shows the importance of partial charge stabilization of the Cys 30 thiolate. With the modeled ⌬38 -40 mutant structurally homologous to WT DsbA, these partial charge interactions can be attributed to WT residues 38 -40. The measured effects (11) are therefore consistent with earlier predictions, which estimated interaction energy between the Cys 30 thiolate and cumulative peptide dipoles over residues 25-43 at about Ϫ20 kJ/mol (9). The measured upward Cys 30 pK a shift of 0.5 for the WT to ⌬38 -40 mutation (11) compares with a value of 0.3 by ⑀ p ϭ 4 calculation, demonstrating that these relatively small shifts are roughly in line with the prediction that extensive charge interactions in DsbA play a significant role in generating the low Cys 30 pK a and oxidizing power. Residues 38 -40 of E. coli DsbA are missing in Vibrio cholerae DsbA, but a proline causes the same overall kink in the protein backbone (20). The higher thiolate pK a for V. cholerae compared with E. coli DsbA (21) is consistent with the results for the E. coli DsbA ⌬38 -40 mutant, but more detailed assessment must await full difference calculations between the two WT proteins.
Measurements of folding stabilities for oxidized and reduced E. coli DsbA WT and mutants in guanidinium chloride (11) have not been used for comparison because the key determinant of stability in these experiments is the folding transition cooperativity (rather than the transition midpoint) from 1.5-2.5 M guanidinium chloride. Matching computations would therefore be required to account for relatively small differences in ionic strength variation at these high denaturing salt concentrations, which is beyond the scope of current methods. In regard to discussions of the link between redox potential, Cys 30 pK a , and reduced/oxidized protein stability, the current calculations are consistent with such a link, with qualitative agreement between the low ⑀ p model for Cys 30 interactions and K eq measurements. The remaining discrepancy, such as underestimation of experimental values with ⑀ p ϭ 4, could signal the breakdown of the various assumptions and/or modeling insufficiency. For example, choice of ⑀ p within the lower range (typically 2-4) remains an open question, and ⑀ p Ͻ 4 would yield higher calculated values. The large difference in scale between calculated values of ⌬⌬⌬G (Table I) for ⑀ p ϭ 4 and higher ⑀ p models is due to the largely through-protein nature of the (38 -40)/thiolate interactions, suggesting that measurements with ⌬38 -40 provide a sensitive test of protein dielectric modeling. Single-site mutations that are predicted to yield Ͼ0.5 pK a shift relative to WT Cys 30 while preserving active site stereochemistry are removal of the His 32 ionizable group or removal of the Gln 97 side chain amide (9).
Full pK a Calculations, pH Titration Curves, and Protein Dielectric-The various ⑀ p models are now employed in full pK a  Oxidizing Power in DsbA 2503 calculations for DsbA WT, thereby introducing the charge burial term (Fig. 1). Extensive studies of computed versus measured pK a s in a range of proteins have revealed two important factors. There exists a subset of amino acids with large pK a shifts (often linked to function) and a much larger set with small pK a shifts tending toward protein stabilization (6). When ⑀ p is varied to give the best match to experiment, the result tends toward a higher value (e.g. ⑀ p ϭ 20), which yields the larger set of small pK a shifts (6). The Cys 30 thiolate of DsbA is an excellent example of the subset of large pK a shifts. Although pH titrations of the remaining ionizable groups of DsbA have not been measured, it is appropriate to make a qualitative study of the effect of ⑀ p variation on the overall form of the pH dependence and to ask whether any of the available models are capable of generating a large and stabilizing ⌬pK a for Cys 30 together with an overall set of small ⌬pK a s that tend toward protein stabilization. Also included is the ⑀ p ϭ 4 model with the suggested empirical modification to account for hydration entropy change upon charge burial in full pK a calculations (8). Distributions of ⌬pK a are shown for the various computational models (Fig. 3). The E s parameter in the modification for single charge center groups has been adjusted to match experimental pK a values for the DsbA Cys 30 and thioredoxin Cys 32 thiolates (8), so that reproduction of this match for the DsbA Cys 30 pK a in the modified ⑀ p ϭ 4 model is expected. However, application of the E s modification is much more than a device for fit to experiment, because it represents a key part of solvation energetics (solvation entropy), and the derived values fall within the range of measured ionic hydration numbers. It is important to ask whether other models can generate the same agreement for Cys 30 and also to analyze the overall distribution of ⌬pK a s. Modification of the higher ⑀ p models, with a term accounting for the favorable hydration entropy contribution on charge burial, would have a much smaller effect than with a lower ⑀ p model, because the magnitude of the modification must not exceed that of the unfavorable Born term to avoid a model that favors general charge group burial. In other words, the higher ⑀ p models cannot escape from a general underestimation of electrostatic interactions that leaves the Cys 30 ⌬pK a close to zero and in contrast with experiment.
With regard to the overall distribution of ⌬pK a s, the unmodified ⑀ p ϭ 4 model shows a large spread, with the extension toward significant destabilization that is characteristic of charge burial in such a model. Application of the modification gives a range of ⌬pK a s tending toward the moderate overall stabilization that is the basis of success in the higher ⑀ p and ⑀ eff ϭ 50 models (Fig. 3). One of the largest ⌬pK a shifts upon modification of the ⑀ p ϭ 4 model is for Cys 30 , arising from the significant thiolate burial within DsbA. Whereas the modified ⑀ p ϭ 4 model can target both the overall ⌬pK a distribution and specific large values, such as that of Cys 30 , the cost of reducing all charge interactions in the higher ⑀ p and ⑀ eff ϭ 50 models is likely to be the omission of those larger ⌬pK a s that may be of functional interest. The presence of large calculated ⌬pK a s other than that of Cys 30 in this qualitative application of the modified ⑀ p ϭ 4 model does not necessarily indicate model breakdown, because they include residues which by various indications may have significantly altered pK a s, such as Glu 24 , Cys 33 (3), and Glu 37 /Glu 38 .
This article shows that a low ⑀ p continuum electrostatic model is consistent with pK a shifts and redox equilibrium measurements for DsbA WT versus mutants (11), reinforcing its value in understanding oxidizing power in this protein family. In the discussion of methods for full pK a calculations, DsbA provides a valuable diversion. It is not well characterized in terms of general pK a measurements, but the large and functionally significant Cys 30 ⌬pK a provides a crucial test that higher ⑀ p models fail. Because the modified ⑀ p ϭ 4 model recovers a reasonable (moderately stabilizing) ⌬pK a profile as well as the Cys 30 ⌬pK a , potential clearly exists for detailed model development against proteins with well characterized pH-dependent properties.