20156
J. Phys. Chem. 1996, 100, 20156-20163
Including Side Chain Flexibility in Continuum Electrostatic Calculations of Protein Titration Paul Beroza and David A. Case* Department of Molecular Biology, The Scripps Research Institute, La Jolla, California 92037 ReceiVed: August 6, 1996; In Final Form: October 10, 1996X
We have extended Monte Carlo procedures for computing statistical averages over protonation states of a protein to include conformational states of the titrating amino acid side chains. This computational method couples side chain motion and protonation with changes in solution pH. Using a continuum electrostatic model for protein titration, we applied this sampling method to calculate titration curves for lysozyme, myoglobin, and hemoglobin. In addition to the X-ray conformation, each titrating site was allowed to reorient to a conformation with maximum solvent accessibility. For all proteins considered, inclusion of these additional conformations improved agreement with experimental measurements for both overall titration and individual pKas. The results suggest that well-solvated orientations of amino acid side chains are an important factor in determining proton binding characteristics of proteins.
1. Introduction The proton titration of an amino acid necessarily changes the amino acid’s net charge. Consequently, it is likely that the free energy of proton binding differs in different electrostatic environments within a protein, resulting in pKa values of titrating side chains that differ from isolated amino acids in solution. Thus, current models for protein titration concentrate on how the protein environment alters the electrostatic potential at titrating amino acids. The most widely adopted models use a continuum approximation, in which the protein and solvent are represented as volumes of continuous dielectric, and solve the PoissonBoltzmann equation1 to determine the electrostatic potential, from which pKa shifts are obtained. Early applications of the model (dating to the 1950s) used simplified geometries and concentrated on interactions between titrating side chains,2,3 while neglecting the polarity of bonds and important self-energy contributions that arise when a site is poorly solvated.4 Advances in numerical methods for solving the PoissonBoltzmann equation5 have permitted continuum models to deal with electrostatics on an atomic level of detail.6 It was found that local protein conformation had a large effect on the energetics of proton binding, and recent work has shown that systematic sampling of the torsional space of titrating side chains also improves the agreement between calculated and experimental pKas.7 Flexibility is clearly important, but it is difficult to include in titration calculations because of the large number of states that must be sampled. Even with the simplest models, in which only the X-ray conformation of each titrating side chain is considered, the number of possible protonation states is 2N (for N titrating sites). If Nc alternative conformations are considered for each titrating site, the number of states is given by (2Nc)N. Moreover, each combination of side chain conformations alters the protein geometry and therefore, in principle, requires a separate calculation of the energy. Thus, an efficient method to estimate electrostatic energies is required, along with a method to restrict the conformational sampling to those states that are most statistically important (i.e., have the lowest free energy). The purpose of this paper is twofold. First, we present a method whereby side chain conformations can be included in X
Abstract published in AdVance ACS Abstracts, December 1, 1996.
S0022-3654(96)02370-2 CCC: $12.00
the Monte Carlo sampling procedure described previously.8 Second, as an application of the method, we calculate pKas for amino acids using two conformations for each side chain: the X-ray coordinates and and an alternative “extended” conformation, based on maximum solvent accessibility. We apply this method to calculate pKas in three proteins: lysozyme, myoglobin, and hemoglobin. For all proteins, the addition of the extended conformation to the ensemble of states improves overall agreement with experiment. 2. Methods 2.1. Continuum Electrostatic Model for Protein Titration. Statistical models for protein titration are based on computing the protonation of an amino acid averaged over the accessible states of the protein:
∑ξ xi(ξ)e-G(ξ)/k T B
〈xi〉 )
Z
(1)
where 〈xi〉 is the average protonation state of titrating site i, xi(ξ) is the protonation state of site i in protein state ξ, G(ξ) is the free energy of this state, and the sum is over all such states. The distribution is normalized by the partition function Z ) ∑ξe-G(ξ)/kBT, where kB is Boltzmann’s constant and T is the temperature. Evaluation of (1) at a number of pH values yields a titration curve.6 It has often been assumed that the conformation of the protein is restricted to that given by the X-ray coordinates, so that the states of the protein are determined solely by the protonation state of the titrating sites.6,8-13 Given this restriction, the state of a protein can be defined by a vector x whose elements specify the number of titratable protons bound to each titrating site:
x ) (x1, x2, ..., xN)
(2)
where xi is the protonation state of site i. The energy of a state x is given by13 © 1996 American Chemical Society
Side Chain Flexibility in Protein Titration
J. Phys. Chem., Vol. 100, No. 51, 1996 20157
G(x) ) -kBT ln 10∑xi(pKint,i - pH) + i
1
(xi + qiu)(xj + qju)Wij ∑ 2 i,j
(3)
where xi is the protonation state of site i in protein state x, pKint,i is the intrinsic pKa of site i, qui is the charge on site i in its unprotonated state, and Wij is the interaction energy between two charged sites. The intrinsic pKas and site-site interactions are calculated from continuum electrostatic theory. Starting from the static atomic coordinates, as given by X-ray crystallography, regions internal to the molecular surface14 are assigned a low dielectric constant while external regions are assigned a high dielectric constant. Intrinsic pKas are defined by adding to the pKa of the amino acid in solution a perturbation determined by the additional electrostatic work required to change the charge distribution of the amino acid from its unprotonated to protonated state in the protein environment (assuming other titrating sites are in their neutral state).3 This perturbation, and the electrostatic interaction between charged titrating sites, are calculated by numerical methods.6,9-13 Once the intrinsic pKas and site-site interactions have been calculated, the titration curve is calculated by computing a Boltzmann average of the 2N possible protonation states, where N is the number of titrating sites. For relatively small values of N, this average can be calculated directly6,15 or from clustering methods.10,16 For larger problems, one can use Monte Carlo methods to sample the low-energy states.8 In this work, we extend the Monte Carlo procedure to include additional conformational states. 2.2. Conformational Flexibility. Perhaps the most significant assumption in the continuum model of protein titration is the absence of flexibility that would allow an amino acid side chain to reorient in response to changes in the electric field. As pH changes and titrating side chains bind or release protons, the electric fields within the protein change. Amino acids in the protein that are not constrained would respond to these changes by orienting to lower the total electrostatic energy of the system. This is particularly relevant for titrating residues on the protein surface, which have considerable freedom to reorient. Recently, You and Bashford included alternative conformations for the titrating sites in the continuum model.7 However, the number of unique combinations of side chain protonation and conformation grows exponentially with the number of alternative conformations allowed for each side chain. Moreover, each unique conformational state of the protein requires, in principle, a separate set of finite difference calculations to determine the intrinsic pKas and site-site interactions. This proves to be computationally impractical.7 Monte Carlo sampling, which was successful extending the number of sites that could be considered in the static model, can also be used to accommodate the additional conformations of titrating amino acids. However, it is necessary to modify the formalism of the energetics and Monte Carlo sampling from that described previously.8 The state vector, x, which consisted of one element for each titrating site, can be extended to include different conformations for each site:
x ) (s1, s2, ..., sN)
(4)
where si is a vector representing not only a protonation state i, but the conformation as well. The state of each titrating site is represented by a subset of elements, only one of which can be nonzero, indicating the occupied state (see Figure 1). Thus,
Figure 1. Modification of the state vector used to define the statistical state of the protein. Each site is represented by a subset s of the state vector x. Odd elements are conformations of the unprotonated side chain; even elements are conformations of the protonated side chain. Only one nonzero element, si, is allowed for each subset.
even when only one conformation is considered for each site, the dimension of the state vector is twice that given by eq 2. If we let si be the nonzero element of vector si, the expression for the energy of a state, given by eq 3 in the static model, can be recast:
G(x) ) -kBT ln 10∑([(si + 1) mod 2](pKmodel,i - pH) + i
N
∑∆G°(si) + 1/2 ∑ W(si, sj) (5) i
i,j)1 i*j
where si is the nonzero element of si, which represents both the conformation and protonation state of site i; pKmodel,i is the pKa for the model compound in solution. Protonated states are referenced to the external proton concentration (i.e., the [(si + 1) mod2] term is nonzero for only protonated states). The ∆G° term includes the Born energy terms and background interactions with nontitrating charges. However, unlike the intrinsic pKa, ∆G° does not include the interactions with neutral titrating side chains; these are treated explicitly in the W(si,sj) term. This energy is referenced to a “model compound energy” that is calculated separately (i.e., ∆G° ) G°(protein) - G°(model)). Following Bashford and Karplus,6 we use the blocked peptide for the model compound, with the conformation found in the protein. Thus, there is a unique model compound for each side chain conformation. The W(si,sj) term is the electrostatic energy between titrating sites i and j when they are in states si and sj (see Figure 2). Note that W(si,sj) is nonzero even if both sites are neutral, because neutral residues have dipoles and higher multipoles (represented by point charges) that interact. Thus, part of the intrinsic pKa has been incorporated into this term. We have assumed that ∆G°i(si) and W(si, sj) are independent of the state of other titrating sites. This is correct in the static model, but need not be if the side chains are flexible because the dielectric boundary depends on the conformation of the whole protein (see You and Bashford for discussion7). In addition, we have not included any conformationally dependent nonelectrostatic energy (aside from the hard sphere model used to determine the extended conformation), although this could be incorporated in the ∆G°i(si) term. This formalism also provides a straightforward means of including tautomeric states of side chains. Tautomers can be
20158 J. Phys. Chem., Vol. 100, No. 51, 1996
Figure 2. Modification of the interaction matrix between titrating sites. In the static model, the interaction between two sites was given by a scalar Wij. In the flexible model, the interaction depends on the protonation and conformational states of the two sites and, therefore, is represented by a matrix wij.
considered as different conformations of a residue, with the proviso that the redundant state not be oversampled. For example, if site i is a histidine residue, the δ and tautomers can be represented by a vector s (a subset of the entire state vector x) of dimension 4. s ) (1 0 0 0) represents the δ tautomer, while s ) (0 0 1 0) represents the tautomer. However, the states (0 1 0 0) and (0 0 0 1) would both refer to the protonated imidazolium ion, and hence are identical, and the Monte Carlo sampling must be modified to sample only one of them. This model is essentially identical to that used by Bashford et al.17 to study histidine tautomerism in myoglobin but provides a straightforward interface to extend this idea to other types of tautomers (such as carboxylic acids) and to alternative conformations. All results reported here included both δ and tautomers for every conformation of each histidine residue. 2.3. Application of the Method. As an example of the Monte Carlo sampling of multiple side chain conformations in the continuum model, we allow two conformations for each titrating amino acid: (1) the X-ray conformation and (2) an “extended” conformation generated by choosing the set of dihedral angles that maximizes the solvent accessibility of the titrating amino acid side chain. The extended conformation will maximize the solvent screening for both the self-energy and the site-site interactions with most other residues (the exception being when the most solvent exposed conformation of one residue increases its interaction with a nearby titrating site). This relatively simple choice of alternative conformation has four important features. (1) It is a computationally tractable and reproducible procedure that does not depend on the details of force fields, charge models, or minimization procedures. (2) It lowers the solvation penalty for residues that are not wellexposed to solvent because they were in their neutral state in the protein crystal from which the X-ray coordinates were derived. This will be especially important at pH values far from the conditions under which the protein was crystallized. (3) It lowers the interaction between sites that are not both charged in the X-ray structure. In effect, this is similar to increasing the internal dielectric constant of the protein (both will lower site-site interactions) but has a larger effect on surface residues. (4) Finally, it chooses what is likely to be a low-energy conformation for a residue that is poorly resolved in the X-ray structure and may be disordered. A set of self-energies and site-site interaction energies are
Beroza and Case calculated for this “extended” structure and included in the statistical averaging used to calculate the titration curve of the protein. Intermediate states, in which some titrating side chains are in their X-ray conformation while others are in the extended conformation, are allowed, even though energies were not explicitly calculated for such intermediate states. We assume that the self-energy of a titrating site in its X-ray (extended) conformation is that given by the calculation with all other titrating side chains in their X-ray (extended) conformation. For site-site interactions we assume that interactions between side chains in mixed conformations are given by the interactions calculated from the extended conformation. This formalism for including protein flexibility differs from those presented previously in that it handles conformational flexibility in a single calculation of statistical averages. Antosiewicz et al.18 computed pKas by taking averages of sets of pK1/2 values obtained from different NMR-derived structures. This neglects the Boltzmann probability weighting associated with the different conformations of each site. You and Bashford7 computed the interaction between two titrating sites from a Boltzmann average of the set of interaction energies obtained by fixing one residue and allowing the other to be flexible. Thus, the interaction energy between two sites was averaged and correlation between protonation and side chain motion was not possible, although the ensemble of conformations used was more complete that that reported here. The Monte Carlo method we have described gives a more correct statistical distribution for the different conformation and protonation states, within the constraints of the ensemble and energy terms used. 2.4. Computational Details. For the lysozyme and myoglobin calculations, we used the same coordinates and parameters that were used in the previous studies to which we compare.7,17 For lysozyme, we used the triclinic structure19 (Brookhaven databank20 entry 1lzt) as minimized by the ALLME0E4 protocol described by You and Bashford.7 As in that paper, we used the charges and radii from the CHARMM 22 molecular mechanics program.21 For myoglobin, we used the sperm whale carbonmonoxymyoglobin X-ray coordinates22 (Brookhaven databank entry 1mbc). Again, following the procedure in the original paper,17 we used charges from the AMBER potential function23 and van der Walls radii reported by Bondi.24 Further structural details (e.g., positioning of polar hydrogens) can be found in the original papers.7,17 The coordinates for the three hemoglobin structures were obtained from the Brookhaven protein databank: R state structure (databank entry 1hho) coordinates for oxyhemoglobin from Shanaan et al.;25 R2 state structure (databank entry 1bbb) coordinates for carbon monoxyhemoglobin from Silva et al.;26 and T state structure (databank entry 2hhb) coordinates for deoxyhemoglobin from Fermi et al.27 Proton positions were determined using the Protonate and pol_h programs from the Amber suite of molecular modeling utilities.28 Unfavorable steric overlaps were removed by 500 iterations of energy minimization (100 steps of steepest descent followed by 400 steps of conjugate gradient) in which atoms were constrained to their X-ray coordinates by a 10 kcal/(mol Å) force. In order not to bias the minimization toward a particular charge distribution, electrostatic interactions were reduced by using a dielectric constant of 1000. Partial charges for the amino acids were taken from Cornell et al.,29 and atomic “Born” radii were taken from the PARSE parameter set.30 The dielectric interface (between the protein and solvent) was assumed to be the molecular surface, defined by Lee and Richards31 as the surface of contact between the protein and a
Side Chain Flexibility in Protein Titration spherical probe of radius 1.4 Å. Regions inside the molecular surface were assigned a dielectric constant of 4, and regions beyond the molecular surface were assigned a dielectric constant of 80. An ionic strength of 100 mM monovalent salt was used for all calculations, and counterions were excluded within 2 Å of the molecular surface. Energy terms in eq 5 were determined using a Green’s function formalism and a finite difference solution to the Poisson-Boltzmann equation, as described previously.13 This approach minimizes the error associated with mapping charges onto a discrete grid, permits focusing to convergence in the self-energy, and was shown to yield terms in in eq 3 that were accurate to 15 11.53 3.95 3.56 6.76 3.12
a “X-ray”: calculations based on X-ray coordinates. “Extended”: calculations based on the X-ray and extended conformations of titrating side chains as described in the text.
were greater that 2 pH units. The upward shifts of pK1/2s for acidic side chains were pronounced, but many pK1/2 values remained outside the experimental pH range, so the effect is not observable in Figure 5. Predicted individual pK1/2 values for histidines are also closer to experiment for the X-ray + extended calculation (see Table 3). Although for some residues the improvement is dramatic (e.g., reduction in error by 2.4 pKa units for His 119), it is
Figure 6. Calculated (dotted and dashed lines) and experimental (solid line) titration curves for R state (oxy) hemoglobin.
typically under 0.5 pKa units. For three of the histidines the pKas the experimental pKa is between that calculated from the X-ray and X-ray + extended methods, which may indicate improved sampling of the accessible states of the protein. 3.3. Overall Titration for Oxyhemoglobin and Deoxyhemoglobin. Changes in hemoglobin’s binding of protons with changes in pH (the Bohr effect) are crucial to its oxygen transport properties. The ligand-bound (oxy and carbon monoxy) states have a lower affinity for protons, likely the result of pKa shifts that result from structural changes accompanying ligand binding. Individual pKas have been measured experimentally by site-directed mutatgenesis, chemical modification, and proton NMR (see ref 39 and references therein). Hemoglobin was the subject of early computational studies using simplified continuum models.40,41 With four subunits and 170 titratable side chains, hemoglobin is a challenging system for modern continuum models. Here we report results for the overall titration curves for the different conformational states using the more sophisticated model described above. Analysis of the individual pKas of amino acids in the different states of hemoglobin and an analysis of the Bohr effect will be the subject of a separate paper. The calculated and experimental proton titration curves for oxy-, carbon monoxy-, and deoxyhemoglobin are shown in Figures 6-8. Experimental values were from measurements in the changes in net protonation with pH for oxy- and deoxyhemoglobin.42 The additive constant needed to obtain the net charge was determined by setting the experimental protonation equal to 0 at the isoionic pH (7.2 for the R state and 7.4 for the T state). The calculated net charge was obtained from
20162 J. Phys. Chem., Vol. 100, No. 51, 1996
Figure 7. Calculated (dotted and dashed lines) and experimental (solid line) titration curves for R2 state (carbon monoxy) hemoglobin.
Figure 8. Calculated (dotted and dashed lines) and experimental (solid line) titration curves for T state (deoxy) hemoglobin.
the average protonation given by the Monte Carlo sampling of the protonation and conformation of the titrating amino acids. For the R and R2 states, including the extended conformations improves agreement with experiment. Not only is the net number of protons titrated in better agreement (as in the myoglobin calculation), but the deviation from the experimental curve lies within one proton of the experimental curve throughout the pH range. This corresponds to an accuracy of significantly better than one proton for each of the four monomers of the protein. For the T (deoxy) state, the calculation that includes alternative conformations lies closer to the experimental curve than the calculation based on the X-ray coordinates alone. However, both calculations lay significantly below the experimental curve, deviating by as much as 3-4 protons near pH 7. Nevertheless, the extended calculation has the correct number of residues titrating in the experimental pH range, but some pKas have apparently been underestimated. In contrast, the static model result remains three protons below the experimental curve at low pH, which indicates a more significant underestimate of individual pKas. 4. Conclusions We have presented a Monte Carlo method for calculating amino acid pKas that incorporates side chain motion. The
Beroza and Case method can, in principle, accommodate many different conformations for each titrating site, although the amount of work required to provide energies for large numbers of states is formidable. As an example of the usefulness of this method, we calculated titration curves for a protein using two structures; that provided by the X-ray coordinates and a structure generated by chosing the orientation of each titrating side chain that maximizes its solvent accessibility. This simple example rivals the much more exhaustive dihedral search followed by global minimization procedure presented by You and Bashford7 but requires far less effort. For the five protein structures considered, significant improvement in overall calculated titration curves was obtained by including an additional orientation for each side chain. The distribution of X-ray conformations and extended conformations is fairly even. For all proteins and all pH values, the number of sites in the X-ray conformation is between 30 and 70%. Thus, both conformations are statistically important, and although calculations based solely on extended structures give better agreement with experiment than those based solely on the X-ray conformation (data not shown), both orientations of titrating side chains are important for correct statistics. Other continuum electrostatic calculations of amino acid pKas in proteins suggest that solvation energetics are not properly taken into account when static X-ray structures are used. Agreement between calculated and experimental pKas was found to improve when the charge on the titrating site was localized on the most solvent accessible atom.9 In addition, use of a dielectric constant of 20 for the protein interior, much higher than the value of 4 that we and others have used,6,10 improves the agreement between calculated and experimental pKas.11 Our results suggest that when side chain flexibility is included in the continuum model, calculated pKas tend to shift toward their values in solution because the charged side chains can be better solvated. It may be that part of the success of these other models results from mimicking this solvation. This is in general accord with molecular dynamics studies of protein dielectric behavior that attribute much of the dielectric response of a protein to motion of charged side chains,43,44 which is taken into account explicitly in our continuum approach. The question of how to sample the conformational states of titrating side chains most efficiently remains open. We have only included one additional conformation based on maximizing the solvent accessibility of the titrating moiety. This improves agreement for both overall titration curves and for most side chains. However, there are certain cases (Glu 35 in lysozyme, for example) in which interactions with other charges are probably underestimated by choosing the conformation with the greatest solvent accessibility. A conformation which maximizes stabilizing interactions with nearby sites would be a logical choice for an additional conformation. Another alternative is to use ensembles of solution structures that are available from NMR studies,18 although this option limits the size of the protein that can be considered. Whatever the method for obtaining the various structures, the Monte Carlo methods presented here can be used to attack the coupled titration problem in a statistically consistent and accurate way. Finally, there is the question of the accuracy of the calculated energies of the individual conformational states sampled. It may be that continuum models are not detailed enough to properly account for the energies that determine proton binding affinities for individual side chains. However, as more conformational space is sampled, the results appear to systematically improve, which suggests that, rather than refining continuum models or
Side Chain Flexibility in Protein Titration adjusting parameters, improving conformational sampling may be a more successful approach. Acknowledgment. This work was supported by National Institutes of Health Grant HL-51030. We thank Don Bashford for many helpful discussions and for providing details of the calculations reported in refs 7 and 17. References and Notes (1) McQuarrie, D. Statistical Mechanics; Harper and Row: New York, 1976. (2) Linderstrom-Lang, K. C. R. TraV. Lab. Carlsberg 1924, 15, 1. (3) Tanford, C.; Kirkwood, J. G. J. Am. Chem. Soc. 1957, 79, 5333. (4) Warshel, A.; Russell, S. T.; Churg, A. K. Proc. Natl. Acad. Sci. U.S.A. 1984, 81, 4785. (5) Warwicker, J.; Watson, H. C. J. Mol. Biol. 1982, 157, 671. (6) Bashford, D.; Karplus, M. Biochemistry 1990, 29, 10219. (7) You, T. J.; Bashford, D. Biophys. J. 1995, 69, 1721. (8) Beroza, P.; Fredkin, D. R.; Okamura, M. Y.; Feher, G. Proc. Natl. Acad. Sci. U.S.A. 1991, 88, 5804. (9) Oberoi, H.; Allewell, N. M. Biophys. J. 1993, 65, 48. (10) Yang, A.-S.; Gunner, M. R.; Sampogna, R.; Sharp, K.; Honig, B. Proteins 1993, 15, 252. (11) Antosiewicz, J.; McCammon, J. A.; Gilson, M. K. J. Mol. Biol. 1994, 238, 415. (12) Beroza, P.; Fredkin, D. R.; Okamura, M. Y.; Feher, G. Biophys. J. 1995, 68, 2233. (13) Beroza, P.; Fredkin, D. R. J. Comput. Chem. 1996, 17, 1229. (14) Richards, F. M. Annu. ReV. Biophys. Bioeng. 1977, 6, 151. (15) Bashford, D.; Karplus, M. J. Phys. Chem. 1991, 95, 9556. (16) Gilson, M. K. Proteins 1993, 15, 266. (17) Bashford, D.; Case, D. A.; Dalvit, C.; Tennant, L.; Wright, P. E. Biochemistry 1993, 32, 8045. (18) Antosiewicz, J.; McCammon, J. A.; Gilson, M. K. Biochemistry 1996, 35, 7819. (19) Hodsdon, J.; Brown, G.; Sieker, L.; Jensen, L. Acta Crystallogr. B 1990, 46, 54. (20) Bernstein, F. C.; Koetzle, T. F.; Williams, G. J. B.; Meyer Jr., E. F.; Brice, M. D.; Rodgers, J. R.; Kennard, O.; Shimanouchi, T.; Tasumi, M. J. Mol. Biol. 1977, 112, 535.
J. Phys. Chem., Vol. 100, No. 51, 1996 20163 (21) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M. J. Comput. Chem. 1983, 4, 187. (22) Kuriyan, J.; Wilz, S.; Karplus, M.; Petsko, G. A. J. Mol. Biol. 1986, 192, 133. (23) Weiner, S. J.; Kollman, P. A.; Nguyen, D. T.; Case, D. A. J. Comput. Chem. 1986, 7, 230. (24) Bondi, A. J. Phys. Chem. 1964, 68, 441. (25) Shanaan, B. J. Mol. Biol. 1983, 171, 31. (26) Silva, M. M.; Rogers, P. H.; Arnone, A. J. Biol. Chem. 1992, 267, 17248. (27) Fermi, G.; Perutz, M. F.; Shaanan, B.; Fourme, R. J. Mol. Biol. 1984, 175, 159. (28) Pearlman, D. A.; Case, D. A.; Caldwell, J. C.; Ross, W. R.; Cheatham III, T. E.; DeBolt, S.; Ferguson, D.; Seibel, G.; Kollman, P. Comput. Phys. Commun. 1995, 91, 1. (29) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz Jr., K. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. J. Am. Chem. Soc. 1995, 117, 5179. (30) Sitkoff, D.; Sharp, K. A.; Honig, B. J. Phys. Chem. 1994, 98, 1978. (31) Lee, B.; Richards, F. M. J. Mol. Biol. 1971, 55, 379. (32) Havel, T. F. Prog. Biophys. Mol. Biol. 1991, 56, 43. (33) Guntert, P.; Braun, W.; Wuthrich, K. J. Mol. Biol. 1991, 217, 517. (34) Macke, T. J. NAB: A Language for Molecular Manipulation, PhD Thesis, The Scripps Research Institute, 1996. (35) Kuramitsu, S.; Hamaguchi, K. J. Biochem. 1980, 87, 1215. (36) Bartik, K.; Redfield, C.; Dobson, C. M. Biophys. J. 1994, 66, 1180. (37) Breslow, E.; Gurd, F. R. N. J. Biol. Chem. 1962, 237, 371. (38) Shire, S. J.; Hanania, G. I. H.; Gurd, F. R. N. Biochemistry 1974, 13, 2967. (39) Ho, C.; Russu, I. M. Biochemistry 1987, 26, 6299. (40) Matthew, J. B.; Hanania, G. I. H.; Gurd, F. R. N. Biochemistry 1979, 18, 1919. (41) Matthew, J. B.; Hanania, G. I. H.; Gurd, F. R. N. Biochemistry 1979, 18, 1928. (42) Rollema, H. S.; de Bruin, S. H.; Janssen, L. H. M.; van Os, G. A. J. J. Biol. Chem. 1975, 250, 1333. (43) Simonson, T.; Perahia, D. Proc. Natl. Acad. Sci. U.S.A. 1995, 92, 1082. (44) Simonson, T.; Brooks III, C. L. J. Am. Chem. Soc. 1996, 118, 8452.
