New Research In
Physical Sciences
Social Sciences
Featured Portals
Articles by Topic
Biological Sciences
Featured Portals
Articles by Topic
 Agricultural Sciences
 Anthropology
 Applied Biological Sciences
 Biochemistry
 Biophysics and Computational Biology
 Cell Biology
 Developmental Biology
 Ecology
 Environmental Sciences
 Evolution
 Genetics
 Immunology and Inflammation
 Medical Sciences
 Microbiology
 Neuroscience
 Pharmacology
 Physiology
 Plant Biology
 Population Biology
 Psychological and Cognitive Sciences
 Sustainability Science
 Systems Biology
Tanford–Kirkwood electrostatics for protein modeling

Edited by Peter S. Kim, Massachusetts Institute of Technology, Cambridge, MA, and approved July 20, 1999 (received for review October 7, 1998)
Abstract
Solvent plays a significant role in determining the electrostatic potential energy of proteins, most notably through its favorable interactions with charged residues and its screening of electrostatic interactions. These energetic contributions are frequently ignored in computational protein design and protein modeling methodologies because they are difficult to evaluate rapidly and accurately. To address this deficiency, we report a revised form of the original Tanford–Kirkwood continuum electrostatic model [Tanford, C. & Kirkwood, J. G. (1957) J. Am. Chem. Soc. 79, 5333–5339], which accounts for the effects of solvent polarization on charged atoms in proteins. The Tanford–Kirkwood model was modified to increase its speed and to improve its sensitivity to the details of protein structure. For the 37 electrostatic selfenergies of the polar sidechains in bovine pancreatic trypsin inhibitor, and their 666 interaction energies, the modified Tanford–Kirkwood potential of mean force differs from a computationally intensive numerical potential (DelPhi) by rootmeansquare errors of 0.6 kcal/mol and 0.08 kcal/mol, respectively. The Tanford–Kirkwood approach makes possible a realistic treatment of electrostatics in computationally demanding protein modeling calculations. For example, pH titration calculations for ovomucoid third domain that model polar sidechain relaxation (including >2 × 10^{23} rotamer conformations of the protein) provide pKa values of unprecedented accuracy.
Substantial advances in protein design (1–3) and homology modeling (4–7) have resulted from the introduction of multicopy sampling algorithms (8), which simultaneously evaluate the fitness of multiple amino acids and sidechain rotamers at each position in a protein structure. Electrostatic energies in current studies have either been disregarded (1, 4–7, 9) or have been approximated by a highly damped Coulomb potential and empirical terms proportional to atomic solvent accessible surface areas (10). These methods ignore the substantial interactions of buried charges in proteins with solvent and incorrectly represent the screening of charge–charge interactions. Electrostatic energies in proteins are not small (11), and inaccurate treatments of the electrostatic potential would be expected to introduce large energetic errors into modeling calculations.
One reason for the omission of electrostatic energies is the lack of methods for their evaluation that are both rapid and accurate. Freeenergy perturbation calculations with explicit solvent molecules (12), as well as finitedifference and boundaryelement methods based on continuum solvent models (13), are thought to well approximate electrostatic free energies, but the extensive computation times they require prohibit their use in many protein modeling calculations. Fast analytical expressions for the electrostatic potential based on the Coulomb field approximation and the generalized Born equation (14) provide remarkably accurate smallmolecule solvation energies but give rise to large errors when applied to bulky solutes such as proteins.
To address these problems, we have revisited the classic Tanford–Kirkwood (TK) electrostatic model (15). The basis for the TK model is conceptually simple (Fig. 1). Regions of space occupied by solvent and solute are treated as dielectric continua. The solventexcluded volume of the protein is equated to a sphere of lowdielectric material, and the solvent is equated to a surrounding highdielectric material. Single charges and charge pairs are mapped from the protein into the sphere, and their interactions are evaluated in the spherical geometry by using Kirkwood’s analytical series solution to the Poisson–Boltzmann equation (16).
We have made three alterations to the TK model to facilitate its use in protein modeling. First, we replace the original TK pointcharge representation by shell charges, so that the electrostatic self energies of buried atoms in proteins can be calculated in a well defined manner (11). Second, we map charges between the protein and sphere geometries on the basis of a Coulomb field integral (17, 18), making the TK model sensitive to the structural details of the protein in the vicinity of a charge. Finally, we substitute Kirkwood’s equations with a complete and precise imagecharge solution to the potential of a charge in a spherical dielectric cavity, which allows fast evaluation of electrostatic self energies and interaction energies (19–21). The modified TanfordKirkwood (MTK) model matches the speed of analytical smallmolecule electrostatic methods but accurately reproduces in proteins the electrostatic free energies of computationally intensive numerical methods.
MATERIALS AND METHODS
Definitions.
The symbols R, ɛ_{s}, and ɛ_{p} denote, respectively, the radius of the lowdielectric sphere (LDS) used in the TK model, the dielectric constant assigned to the solvent, and the dielectric constant assigned to the protein interior. The symbols q_{i}, b_{i}, and d_{i} denote the charge, Born radius, and depth (distance to the surface of the LDS; negative when an atom is outside the sphere) of atom i. The symbol r_{ij} denotes the distance between atoms i and j. The symbols S_{k} and dS_{k} denote the geometricmean radius and thickness of a probe shell indexed k. The symbols ℱ_{k}^{i, pro} and ℱ_{k}^{i, sph} denote the fractional surface area of probe shell k concentric with atom i that falls inside the volume of low dielectric material in the protein and sphere geometries respectively. D⃗(x⃗) and Φ(x⃗) denote the dielectric displacement field and electrostatic potential as functions of spatial position (we assume a linear, osotropic dielectric medium). The self energy of atom i, W_{i}, is taken to be the electrostatic free energy with atom i charged and all other atoms of the solute neutral. The interaction potential, I_{ij}, of atoms i and j is defined I_{ij} = W_{ij} − W_{i} − W_{j}.
Overview.
Solute atoms are modeled as shells of charge (total charge q_{i}) characterized by a Born radius (b_{i}). The electrostatic free energy, W^{Elec}, is expressed as a sum over interactions between shells: W^{Elec} = 1/2 Σ_{i,j} q_{i}Φ_{j}. To determine q_{t}Φ_{j}, charges i and j are mapped from the solute into the LDS, where Φ_{j} is approximated by the potentials of a set of suitably chosen image charge shells (of charge q_{j}^{im} and radii b_{j}^{im}). If neither shell i nor j crosses the LDS boundary, then q_{i}Φ_{j} = Σ_{shells}^{image} q_{i}q_{j}^{im}/α_{ij}^{im}, where 1 F_{in}/F_{out} denotes the surface fraction of image shell j that lies inside/outside of shell i.
Shell charges crossing the LDS boundary are treated as transparent to the two dielectric media (22) and are divided by the boundary into two independent charge caps. Each cap is replaced by a set of image caps, whose electrostatic potentials are evaluated according to Jackson (page 128 of ref. 19; available as supplemental data on the PNAS web site, www.pnas.org). If shells i and j both cross the dielectric boundary with r_{ij} < (b_{i} + b_{j} + 2 Å), the potential of one shell is sampled over the surface of the other shell. To accelerate calculations, all interactions that occur at r_{ij} ≥ (b_{i} + b_{j} + 2 Å) are treated as though between point charges, regardless of the proximity of either i or j to the dielectric interface.
Charge Mapping.
For an atom i in a solute molecule, we seek a position in the LDS that preserves electrostatic self energy expressed as a spatial integral over energy density: 2 We replace the electric displacement with the displacement field of a point charge in a uniform dielectric medium appropriate to the domain of integration. This simplification is known as the Coulomb field approximation; implications of its use are discussed in refs. 17 and 18.
The integrals in Eq. 2 are approximated by discretely sampling points on a series of shells external to and concentric with atom i (Fig. 1) (supplementary material to ref. 15). Each point may be defined to be located in solvent or solute by using the prescription of ref. 23. Because of the spherical symmetry of the Coulomb field, knowledge of the fraction of each shell that is within the solute is sufficient to determine the self energy: 3 This approximate W_{i} is used only for mapping; a more accurate selfenergy is derived by using the image charge solution described below.
We define two methods for mapping atoms in proteins. Fine mapping is used for interactions within aminoacid residues whereas coarse mapping is used for interactions between residues.
Fine mapping.
To map a single atom into the LDS, we find the d_{i} and R, which minimize the quantity 4 by using a golden section search (24). This is equivalent to demanding that the Coulomb field integral be preserved between geometries. Pairs of charges i j are mapped into the LDS by minimizing at fixed r_{ij} the quantity (ℛ(i) + ℛ(j)) with respect to d_{i}, d_{j}, and R. We accomplish this via a Powell search (24). For pairs of charges, ℱ _{k}^{i, sph} is defined to be the fraction of shell k inside either the LDS or atom j (see supplemental data at www.pnas.org). Equations for computing this threebody ℱ_{k}^{i, sph} have been described by Richmond (25).
Coarse mapping.
We fix R such that 4/3πR^{3} gives the Connolly volume of the protein (26). For each atom, we minimize the quantity 5 yielding a complete set of d_{i}s. The mapping of pairs of charges i j is geometrically specified by d_{i}, d_{j}, R, and r_{ij}. There exist rare instances in which r_{ij} is incompatible with d_{i}, d_{j}, and R. On these occasions, when r_{ij} < d_{t} − d_{j} or r_{ij} > (2R − d_{i} − d_{j}), we adjust d_{i} by D⋅d_{i}/(d_{i} + d_{j}) and d_{j} by D⋅d_{j}/(d_{i}+d_{j}), where D = d_{i} −d_{j} − r_{ij} in the former case and D = r_{ij} − 2R + d_{i} + d_{j} in the latter. This method of partitioning D preferentially modifies the depths of atoms furthest from the dielectric boundary.
ImageCharge Solution.
Let x_{q}, x_{obs} denote the radial coordinates of a real charge and an observation point relative to the LDS center. The potential inside the LDS (x_{obs} < R) due to a charge inside the dielectric interface (x_{q} < R) is expressed exactly as a series of Legendre polynomials, which may be decomposed into Coulombic and solvent reaction contributions. It has been shown that the reaction potential can be approximated by an image point charge (20), either alone or in conjunction with the image shell charge at the dielectric boundary (21). Following Abagyan and Totrov (21), we approximate the electric potential at x_{obs} < R due to a charge at x_{q} < R with two point charges and a shell charge at the dielectric boundary: 6a where x_{inv} denotes the inverse point, colinear with x_{q} and the LDS center, at radius R^{2}/x_{q}. If either x_{q} > R or x_{obs} > R, the series solution for the potential takes on a different form. Analogous manipulations yield the following image charge sets, which define over all space the potential of arbitrarily placed charges: 6b 6c 6d Application of Eqs. 6a–6d to charge shells and caps gives rise to image shells and caps (see supplementary material). The electrostatic interaction of atoms i and j is obtained by summing in vacuo the interactions of atom i with the appropriate charge set representing the potential generated by atom j.
The image charge method does not account for salt effects. In principle, these effects could be incorporated in the MTK model by recourse to the original polynomial expressions for the electrostatic free energy derived by Kirkwood (16).
Small Molecules and Bovine Pancreatic Trypsin Inhibitor (BPTI).
All smallmolecule coordinates were kindly provided by D. Sitkoff (BristolMyers Squibb). Free energies of cavity formation in water for the small molecule were taken directly from ref. 27. For BPTI studies, the 4PTI structure (28) with no bound water molecules was used. Sidechain self energies were computed as the electrostatic free energy of the charged side chain in an otherwise neutral protein. Sidechain interaction energies were computed as the difference between the electrostatic free energy with both side chains charged and the sum of two sidechain self energies. DelPhi calculations were carried out on a 0.34Å grid at 0 ionic strength with full Coulombic boundary conditions and a 1.4Å probe sphere. The electrostatic free energy was evaluated as = 1/2Σ _{i,j}q_{i}q_{j}/ɛα_{ij}and α_{ij}defined by Eq. 1. DelPhi reactionfield energies were used for the testcharge calculations. Grid focusing and the DelPhi grid energies (29) were used for calculation of sidechain self energies and interaction energies in BPTI.
Timings.
The time (T_{total}) required to obtain electrostatic energies for a modeling calculation can be written in two terms that grow linearly and quadratically with the total number of rotamers n: 7 At a minimum, L and Q represent the average times required to calculate rotamer self and interaction energies respectively. L also includes the time spent obtaining depths for the MTK model. From the BPTI calculations, we have for the MTK model L = 11 sec/residue and Q = 0.13 sec/pair. For DelPhi, we find L = 111 sec/residue and Q = 106 sec/pair. The quadratic term dominates in most multicopy sampling problems (n > 170). In this regime, the MTK potential is evaluated nearly three orders of magnitude faster than the DelPhi potential.
pKa Values.
For pKa calculations, the protein was described by a conformational ensemble consisting of a fixed backbone and a probabilityweighted set of allowable rotamers at each position. Backbone coordinates were taken from 1TUR (30) or 1PPF (31), as were sidechain coordinates in calculations using fixed structures. The rotamer library of Tuffery et al. (32), as modified by Koehl and Delarue (6), was used to represent sidechain conformational freedom in calculations incorporating sidechain relaxation. Sidechain flexibility was modeled for all amino acids except alanine, glycine, isoleucine, leucine, proline, and valine. Side chain coordinates were built by using charmm19 geometric parameters (33). Protons for acidic side chains were placed trans to the preceding methylene group. One deprotonated rotamer was included for each protonated rotamer to prevent entropic bias for the protonated state (34).
The energy function for all calculations consisted of the MTK electrostatic potential, W ^{Elec}, the charmm19 Lennard–Jones potential (33), W ^{LJ}, and a protonation potential, W ^{H,p}. 8 The protonation potential is derived from a thermodynamic cycle connecting a titrating site in a protein to an equivalent site in an isolated Nformyl, Nmethyl amino amide (FMAA) of known pKa (35, 36) (see supplemental data). The protonation potential consists of the transfer free energy of the deprotonated FMAA from the protein into aqueous solution, W _{deprotonated}^{p→aq}, the free energy of protonation of the FMAA in aqueous solution, W ^{H,aq}, and the transfer free energy of the protonated FMAA back into the protein, W _{protonated}^{aq→p}: 9 W ^{H,aq} is evaluated as 2.3⋅RT(pH − IpK_{a}), where IpK_{a} is the intrinsic pKa of the appropriate FMAA. Transfer free energies were calculated as the difference in the potential energy of the FMAA in the protein and free in solution. In all cases, the φ and ψ dihedral angles of the isolated FMAA were set to the values in the protein, and the sidechain dihedral angles were set to the values of the most commonly occurring rotamer. Intrinsic pKa values for the formulated amino amides were assigned as Glu, 4.4; Asp, 4.0; His, 6.3; peptide C terminus, 3.8 (37).
A meanfield method (6) was used to refine the rotamer probabilities such that the free energy of the ensemble was minimized. Because of the inclusion of a protonation potential, the rotamer distributions that minimized the free energy varied with pH. Each titrating site’s pKa was assigned as the pH at which protonated and deprotonated states were equally populated. These values represent averages over many rotamer states of the protein, each with a likelihood equal to the product of the individual rotamer probabilities (see ref. 6).
For all calculations, experimentally determined atomic coordinates were used to define the shape of the lowdielectric region. Consequently, the effects of rotamer changes on the dielectric boundary were not treated, and the pairwise factorization of the electrostatic potential was preserved.
RESULTS
Tests of the Modified TK Potential.
We have introduced modifications into the TK electrostatic model that improve charge representation, charge mapping, and the speed with which the electric potential can be evaluated. Our first test of the MTK potential was to calculate small molecule vacuumtowater transfer free energies (Table 1). Using the parameters for solvation energy (PARSE) set of Born radii and partial charges (27), the calculated transfer free energies for a set of 53 canonical compounds (the training set of ref. 27) ranged over 80 kcal/mol with rms and (maximum) errors from experimentally measured values of 1.05 (3.54) kcal/mol.
We next tested the accuracy of the MTK model in the protein BPTI. All energies calculated by the MTK method were compared with values derived from the program DelPhi (29), which is based on a finitedifference solution to the Poisson–Boltzmann equation (FDPB). The protein dielectric constant was set to 4 and the solvent dielectric constant to 80. In Fig. 2B, the selfenergy of a 1Å shell charge passing along a line through the center of mass of BPTI is plotted against position. The MTK potential of mean force exhibits close agreement with the FDPB potential, even in the rapidly varying region near the molecular surface.
As a further test, we superimposed a 4Å grid of test charges onto BPTI originating at the center of mass and aligned along the principle axes of inertia of the protein (Fig. 2A). A 1Å shell with one electron unit of charge was placed at each grid point. The MTK self energies of the 441 test charges give an Rfactor with respect to the DelPhi self energies of 8% (Rfactor = Σ_{i}W_{i}^{MTK} − W_{i}^{DelPhi}/Σ_{i}W_{i}^{DelPhi}). Thus, the summed absolute values of selfenergy errors represented 8% of the cumulative DelPhi selfenergy. By comparison, the Coulomb field integral as implemented in the analytical algorithm ace (18) gives a selfenergy Rfactor with respect to DelPhi of 25%. As illustrated in Fig. 2B, the Coulomb field integral tends to underestimate selfenergies for buried charges and overestimate them for surface charges.
The interaction energies of each test charge with a test charge at the center of mass also were evaluated. The MTK interaction energies give an Rfactor with respect to DelPhi of 15% (i.e., the cumulative unsigned error is 15% of the summed unsigned interaction energies). By comparison, the generalized Born equation (14) (based on MTK selfenergies) gives an Rfactor of 42%. The generalized Born equation reliably describes the interaction of charges separated by long distances in proteins, but it incurs large errors for nearby charges.
Finally, we measured the electrostatic self and interaction energies of the entire set of 37 polar side chains in BPTI. These values represent the electrostatic inputs required for multicopy sampling algorithms. Relative to DelPhi, the MTK potential of mean force gives rms and (maximum) errors of 0.6 (1.27) kcal/mol for the 37 self energies (Fig. 3A) and rms and (maximum) errors of 0.08 (1.21) kcal/mol for the 666 sidechain interaction energies (Fig. 3B).
pKa Calculations.
As a functional application of the MTK model, we calculated pKa values for titratable groups in turkey ovomucoid third domain, which have been measured experimentally at low ionic strength (38, 39). Protonated and deprotonated forms of each titratable sidechain were treated as rotamers of the same amino acid. The free energy for sidechain protonation in the protein was calculated relative to the free energy for protonation of the corresponding Nformyl Nmethyl amide free in solution (40) (see Materials and Methods). At pH values between 1 and 10, the protonation states of the interacting titratable groups were relaxed by a conventional meanfield optimization algorithm (6). The pH at which each titratable group became 50% deprotonated was designated as its pKa. Although the meanfield approximation is known to give imperfect titration curves for strongly interacting sites of degenerate pKa (36, 41), we wished to evaluate the electrostatic potential under conditions used for protein sidechain repacking studies.
Following the convention of ref. 37, the protein conformation was represented by fixed coordinates, the protein dielectric constant was evaluated as 20, and atoms were assigned Born radii and partial charges according to the PARSE parameter set (27). The pKa values for ovomucoid third domain calculated with the modified TK model give an rms deviation from the experimental values of 0.62 pH units and a maximum error of 1.0 pH unit. These values compare favorably (Table 2) with the best results obtained by conventional FDPB methods [rms error 0.59 pH units (37)]. The MTK pKa values are significantly better than the null model, which assumes no pKa shifts from intrinsic values (rms error 1.1 pH units).
SideChain Conformational Freedom.
By coupling the speed of the MTK methodology with a rotamer repacking algorithm, we have examined how global sidechain motions perturb the calculated pKa values of titratable sites in ovomucoid third domain (Table 2). All polar side chains were allowed full rotamer freedom, giving rise to >2 × 10^{23} possible rotamer conformations of the protein. At pH values between 1 and 10, the rotamer distribution was refined to the freeenergy minimum defined by the MTK potential, the Lennard–Jones potential, and a protonation potential (see Materials and Methods). The pH at which each titratable group became 50% deprotonated was designated as its pKa. With a protein dielectric constant of 20, sidechain conformational sampling produced modest improvement in the agreement between calculated and observed pKa values (Table 2) (Δrms = 0.09 pH units). However, when the protein interior was assigned a dielectric constant of 4, the pKa values calculated with sidechain flexibility were dramatically improved relative to pKa’s calculated with a static structure (Δrms = 1.42 pH units), yielding computed values of unprecedented accuracy (rms error 0.48 pH units). In conjunction with explicitly modeled sidechain motion, the dielectric response of the protein is better represented by a constant of 4 than by a constant of 20.
DISCUSSION
To fully represent the electrostatic properties of proteins, it is necessary to model protein structural plasticity. For example, although pKa calculations on proteins have generally been carried out with static structures, charge changes at titratable sites are known to induce substantial rearrangements in protein sidechain conformation (42). These structural adjustments presumably diminish the energetic penalties associated with protonation and deprotonation. Unfortunately, the computation times required by conventional electrostatic potentials prohibit modeling of such global sidechain motions. Consequently, fixedstructure calculations have been performed. The dielectric response that would result from sidechain rearrangements has been represented implicitly by assigning an artificially large dielectric constant to the protein. Whereas the bulk dielectric constant of crystalline acetamide (a molecular analogue of the protein backbone) is 4, a protein dielectric constant of 20 is conventionally used for fixedstructure pKa calculations (37). Because it subsumes the effects of structural plasticity, the protein dielectric constant has come to be viewed as an adjustable parameter that depends on spatial position within the protein and on the electrostatic property being studied (43). Efforts to incorporate limited sidechain flexibility into electrostatic calculations have been reported [including allowance for hydroxyl rotamer relaxation (34), inclusion of multiple dihedral conformers at titrating sites (35, 44), and sampling of multiple protein conformations along a molecular dynamics trajectory (45–49)], but no methods exist to model the global rearrangements in protein structure that accompany pH changes.
The MTK potential was developed specifically to allow an accurate treatment of electrostatics in multicopy sampling calculations that treat global variation in sidechain chemistry and conformation. Using the MTK potential in conjunction with a sidechain repacking algorithm, we demonstrate that changes in the rotamer populations of polar residues contribute substantially to the dielectric response of a protein. The accuracy of calculated pKa values for ovomucoid third domain are improved by modeling these rearrangements. Moreover, when structural reorganization is treated microscopically, the remaining dielectric response of the protein is well represented by the dielectric constant of bulk crystalline acetamide (ɛ_{p} = 4). Macroscopic and microscopic estimates of the protein dielectric constant are thus reconciled. When the effects of protein structural reorganization are factored out of the electrostatic model, the protein dielectric constant can be treated as a uniform transferable property of the protein solute, as envisioned in the classical theory of dielectrics. We expect that the sidechain repacking approach presented here will prove generally useful for the computational analysis of sequence mutations in proteins. Indeed, fixedstructure calculations predict the destabilization resulting from single chargetoneutral residue substitutions to be much larger than is measured experimentally (50).
Current protein design methods sidestep the consideration of large electrostatic energies by restricting polar residues to the surface of designed proteins (9, 10, 51). This simplification would appear to be justified by the argument that buried polar interactions universally destabilize proteins (52). In at least two situations, however, buried charges and hence accurate electrostatic potentials are required. First, buried polar residues are essential for enforcing structural specificity in some protein folds, even if they are also destabilizing. For example, core asparagine residues have been shown to control the strand number (53), helix orientation (54, 55), and specificity of helix association (56) in dimeric coiledcoils. Furthermore, buried polar residues play a central role in catalysis by protein enzymes (57). As the breadth of computational protein design grows to encompass protein function and the specificity of proteinprotein interactions, the electrostatic potential will lie at the heart of the problem.
Acknowledgments
The authors thank P. Koehl and R. L. Baldwin for stimulating discussion and criticism throughout the course of this work and D. Sitkoff for generously providing smallmolecule coordinates. This research was supported by a junior faculty award from the Howard Hughes Medical Institute and a grant from the Chicago Community Trust to P.B.H. P.B.H. is a Searle scholar and a Terman fellow.
ABBREVIATIONS
 BPTI,
 bovine pancreatic trypsin inhibitor;
 TK,
 Tanford–Kirkwood;
 MTK,
 modified TK;
 LDS,
 lowdielectric sphere;
 FDPB,
 finitedifference solution to the Poisson–Boltzmann equation;
 FMAA,
 Nformyl, Nmethyl amino amide
 Received October 7, 1998.
 Copyright © 1999, The National Academy of Sciences
References
 ↵

 Dahiyat B I,
 Mayo S L
 ↵
 Saven J G,
 Wolynes P G
 ↵
 ↵
 ↵
 ↵
 ↵
 Harbury P B,
 Plecs J J,
 Tidor B,
 Alber T,
 Kim P S
 ↵
 ↵
 ↵
 Allen M P,
 Tildesley D J
 ↵
 Honig B,
 Nicholls A
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 Jackson J D
 ↵
 Friedman H L
 ↵
 ↵
 ↵
 ↵
 Press W H,
 Teukolsky S A,
 Vetterling W T,
 Flannery B P
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵

 Sham Y Y,
 Chu Z T,
 Warshel A
 ↵
 ↵
 ↵
 ↵
 ↵
 Harbury P B,
 Zhang T,
 Kim P S,
 Alber T
 ↵
 ↵
 ↵
 Zeng X,
 Herndon A M,
 Hu J C
 ↵
Shan, S. O. & Herschlag, D. (1999) Methods Enzymol, in press.