Search

The development of models of macromolecular electrostatics capable of delivering improved fidelity to quantum mechanical calculations is an active field of research in computational chemistry. Most molecular force field development takes place in the context of models with full Cartesian coordinate degrees of freedom. Nevertheless, a number of macromolecular modeling programs use a reduced set of conformational variables limited to rotatable bonds. Efficient algorithms for minimizing the energies of macromolecular systems with torsional degrees of freedom have been developed with the assumption that all atom-atom interaction potentials are isotropic. We describe novel modifications to address the anisotropy of higher order multipole terms while retaining the efficiency of these approaches. In addition, we present a treatment for obtaining derivatives of atom-centered tensors with respect to torsional degrees of freedom. We apply these results to enable minimization of the Amoeba multipole electrostatics potential in a system with torsional degrees of freedom, and validate the correctness of the gradients by comparison to finite difference approximations. In the interest of enabling a complete model of electrostatics with implicit treatment of solvent-mediated effects, we also derive expressions for the derivative of solvent accessible surface area with respect to torsional degrees of freedom.

Background: Sequence-specific binding by transcription factors (TFs) plays a significant role in the selection and regulation of target genes. At the protein:DNA interface, amino acid side-chains construct a diverse physicochemical network of specific and non-specific interactions, and seemingly subtle changes in amino acid identity at certain positions may dramatically impact TF:DNA binding. Variation of these specificity-determining residues (SDRs) is a major mechanism of functional divergence between TFs with strong structural or sequence homology., Methods: In this study, we employed a combination of high-throughput specificity profiling by SELEX and Spec-seq, structural modeling, and evolutionary analysis to probe the binding preferences of winged helix-turn-helix TFs belonging to the OmpR sub-family in Escherichia coli ., Results: We found that E. coli OmpR paralogs recognize tandem, variably spaced repeats composed of "GT-A" or "GCT"-containing half-sites. Some divergent sequence preferences observed within the "GT-A" mode correlate with amino acid similarity; conversely, "GCT"-based motifs were observed for a subset of paralogs with low sequence homology. Direct specificity profiling of a subset of OmpR homologues (CpxR, RstA, and OmpR) as well as predicted "SDR-swap" variants revealed that individual SDRs may impact sequence preferences locally through direct contact with DNA bases or distally via the DNA backbone., Conclusions: Overall, our work provides evidence for a common structural code for sequence-specific wHTH:DNA interactions, and demonstrates that surprisingly modest residue changes can enable recognition of highly divergent sequence motifs. Further examination of SDR predictions will likely reveal additional mechanisms controlling the evolutionary divergence of this important class of transcriptional regulators., Competing Interests: COMPLIANCE WITH ETHICS GUIDELINES The authors Adam P. Joyce and James J. Havranek declare that they have no conflict of interests. This article does not contain any studies with human or animal subjects performed by any of the authors.

Cooperative binding of transcription factors is known to be important in the regulation of gene expression programs conferring cellular identities. However, current methods to measure cooperativity parameters have been laborious and therefore limited to studying only a few sequence variants at a time. We developed Coop-seq (cooperativity by sequencing) that is capable of efficiently and accurately determining the cooperativity parameters for hundreds of different DNA sequences in a single experiment. We apply Coop-seq to 12 dimer pairs from the Sox and POU families of transcription factors using 324 unique sequences with changed half-site orientation, altered spacing and discrete randomization within the binding elements. The study reveals specific dimerization profiles of different Sox factors with Oct4. By contrast, Oct4 and the three neural class III POU factors Brn2, Brn4 and Oct6 assemble with Sox2 in a surprisingly indistinguishable manner. Two novel half-site configurations can support functional Sox/Oct dimerization in addition to known composite motifs. Moreover, Coop-seq uncovers a nucleotide switch within the POU half-site when spacing is altered, which is mirrored in genomic loci bound by Sox2/Oct4 complexes., (© The Author(s) 2016. Published by Oxford University Press on behalf of Nucleic Acids Research.)

Cre recombinase catalyzes the cleavage and religation of DNA at loxP sites. The enzyme is a homotetramer in its functional state, and the symmetry of the protein complex enforces a pseudo-palindromic symmetry upon the loxP sequence. The Cre-lox system is a powerful tool for many researchers. However, broader application of the system is limited by the fixed sequence preferences of Cre, which are determined by both the direct DNA contacts and the homotetrameric arrangement of the Cre monomers. As a first step toward achieving recombination at arbitrary asymmetric target sites, we have broken the symmetry of the Cre tetramer assembly. Using a combination of computational and rational protein design, we have engineered an alternative interface between Cre monomers that is functional yet incompatible with the wild-type interface. Wild-type and engineered interface halves can be mixed to create two distinct Cre mutants, neither of which are functional in isolation, but which can form an active heterotetramer when combined. When these distinct mutants possess different DNA specificities, control over complex assembly directly discourages recombination at unwanted half-site combinations, enhancing the specificity of asymmetric site recombination. The engineered Cre mutants exhibit this assembly pattern in a variety of contexts, including mammalian cells., (© The Author(s) 2015. Published by Oxford University Press on behalf of Nucleic Acids Research.)

Leber congenital amaurosis 9 (LCA9) is an autosomal recessive retinal degeneration condition caused by mutations in the NAD(+) biosynthetic enzyme NMNAT1. This condition leads to early blindness but no other consistent deficits have been reported in patients with NMNAT1 mutations despite its central role in metabolism and ubiquitous expression. To study how these mutations affect NMNAT1 function and ultimately lead to the retinal degeneration phenotype, we performed detailed analysis of LCA-associated NMNAT1 mutants, including the expression, nuclear localization, enzymatic activity, secondary structure, oligomerization, and promotion of axonal and cellular integrity in response to injury. In many assays, most mutants produced results similar to wild type NMNAT1. Indeed, NAD(+) synthetic activity is unlikely to be a primary mechanism underlying retinal degeneration as most LCA-associated NMNAT1 mutants had normal enzymatic activity. In contrast, the secondary structure of many NMNAT1 mutants was relatively less stable as they lost enzymatic activity after heat shock, whereas wild type NMNAT1 retains significant activity after this stress. These results suggest that LCA-associated NMNAT1 mutants are more vulnerable to stressful conditions that lead to protein unfolding, a potential contributor to the retinal degeneration observed in this syndrome., (© 2015 by The American Society for Biochemistry and Molecular Biology, Inc.)

Protein:DNA interactions are essential to a range of processes that maintain and express the information encoded in the genome. Structural modeling is an approach that aims to understand these interactions at the physicochemical level. It has been proposed that structural modeling can lead to deeper understanding of the mechanisms of protein:DNA interactions, and that progress in this field can not only help to rationalize the observed specificities of DNA-binding proteins but also to allow researchers to engineer novel DNA site specificities. In this review we discuss recent developments in the structural description of protein:DNA interactions and specificity, as well as the challenges facing the field in the future., (© The Author 2014. Published by Oxford University Press. All rights reserved. For permissions, please email: journals.permissions@oup.com.)

Computational protein design relies on several approximations, including the use of fixed backbones and rotamers, to reduce protein design to a computationally tractable problem. However, allowing backbone and off-rotamer flexibility leads to more accurate designs and greater conformational diversity. Exhaustive sampling of this additional conformational space is challenging, and often impossible. Here, we report a computational method that utilizes a preselected library of native interactions to direct backbone flexibility to accommodate placement of these functional contacts. Using these native interaction modules, termed motifs, improves the likelihood that the interaction can be realized, provided that suitable backbone perturbations can be identified. Furthermore, it allows a directed search of the conformational space, reducing the sampling needed to find low energy conformations. We implemented the motif-based design algorithm in Rosetta, and tested the efficacy of this method by redesigning the substrate specificity of methionine aminopeptidase. In summary, native enzymes have evolved to catalyze a wide range of chemical reactions with extraordinary specificity. Computational enzyme design seeks to generate novel chemical activities by altering the target substrates of these existing enzymes. We have implemented a novel approach to redesign the specificity of an enzyme and demonstrated its effectiveness on a model system., (© 2014 The Authors Protein Science published by Wiley Periodicals, Inc. on behalf of The Protein Society.)

The Rosetta molecular modeling software package provides experimentally tested and rapidly evolving tools for the 3D structure prediction and high-resolution design of proteins, nucleic acids, and a growing number of non-natural polymers. Despite its free availability to academic users and improving documentation, use of Rosetta has largely remained confined to developers and their immediate collaborators due to the code's difficulty of use, the requirement for large computational resources, and the unavailability of servers for most of the Rosetta applications. Here, we present a unified web framework for Rosetta applications called ROSIE (Rosetta Online Server that Includes Everyone). ROSIE provides (a) a common user interface for Rosetta protocols, (b) a stable application programming interface for developers to add additional protocols, (c) a flexible back-end to allow leveraging of computer cluster resources shared by RosettaCommons member institutions, and (d) centralized administration by the RosettaCommons to ensure continuous maintenance. This paper describes the ROSIE server infrastructure, a step-by-step 'serverification' protocol for use by Rosetta developers, and the deployment of the first nine ROSIE applications by six separate developer teams: Docking, RNA de novo, ERRASER, Antibody, Sequence Tolerance, Supercharge, Beta peptide design, NCBB design, and VIP redesign. As illustrated by the number and diversity of these applications, ROSIE offers a general and speedy paradigm for serverification of Rosetta applications that incurs negligible cost to developers and lowers barriers to Rosetta use for the broader biological community. ROSIE is available at http://rosie.rosettacommons.org.

Accurate energy functions are critical to macromolecular modeling and design. We describe new tools for identifying inaccuracies in energy functions and guiding their improvement, and illustrate the application of these tools to the improvement of the Rosetta energy function. The feature analysis tool identifies discrepancies between structures deposited in the PDB and low-energy structures generated by Rosetta; these likely arise from inaccuracies in the energy function. The optE tool optimizes the weights on the different components of the energy function by maximizing the recapitulation of a wide range of experimental observations. We use the tools to examine three proposed modifications to the Rosetta energy function: improving the unfolded state energy model (reference energies), using bicubic spline interpolation to generate knowledge-based torisonal potentials, and incorporating the recently developed Dunbrack 2010 rotamer library (Shapovalov & Dunbrack, 2011)., (Copyright © 2013 Elsevier Inc. All rights reserved.)

The ability to engineer novel protein folds, conformations, and enzymatic activities offers enormous potential for the development of new protein therapeutics and biocatalysts. However, many de novo and redesigned proteins exhibit poor hydrophobic packing in their predicted structures, leading to instability or insolubility. The general utility of rational, structure-based design would greatly benefit from an improved ability to generate well-packed conformations. Here we present an automated protocol within the RosettaDesign framework that can identify and improve poorly packed protein cores by selecting a series of stabilizing point mutations. We apply our method to previously characterized designed proteins that exhibited a decrease in stability after a full computational redesign. We further demonstrate the ability of our method to improve the thermostability of a well-behaved native protein. In each instance, biophysical characterization reveals that we were able to stabilize the original proteins against chemical and thermal denaturation. We believe our method will be a valuable tool for both improving upon designed proteins and conferring increased stability upon native proteins.

We have recently completed a full re-architecturing of the ROSETTA molecular modeling program, generalizing and expanding its existing functionality. The new architecture enables the rapid prototyping of novel protocols by providing easy-to-use interfaces to powerful tools for molecular modeling. The source code of this rearchitecturing has been released as ROSETTA3 and is freely available for academic use. At the time of its release, it contained 470,000 lines of code. Counting currently unpublished protocols at the time of this writing, the source includes 1,285,000 lines. Its rapid growth is a testament to its ease of use. This chapter describes the requirements for our new architecture, justifies the design decisions, sketches out central classes, and highlights a few of the common tasks that the new software can perform., (© 2011 Elsevier Inc. All rights reserved.)

Modeling the conformational changes that occur on binding of macromolecules is an unsolved challenge. In previous rounds of the Critical Assessment of PRediction of Interactions (CAPRI), it was demonstrated that the Rosetta approach to macromolecular modeling could capture side chain conformational changes on binding with high accuracy. In rounds 13-19 we tested the ability of various backbone remodeling strategies to capture the main-chain conformational changes observed during binding events. These approaches span a wide range of backbone motions, from limited refinement of loops to relieve clashes in homologous docking, through extensive remodeling of loop segments, to large-scale remodeling of RNA. Although the results are encouraging, major improvements in sampling and energy evaluation are clearly required for consistent high accuracy modeling. Analysis of our failures in the CAPRI challenges suggest that conformational sampling at the termini of exposed beta strands is a particularly pressing area for improvement., (© Wiley-Liss, Inc.)

A long-standing goal of computational protein design is to create proteins similar to those found in Nature. One motivation is to harness the exquisite functional capabilities of proteins for our own purposes. The extent of similarity between designed and natural proteins also reports on how faithfully our models represent the selective pressures that determine protein sequences. As the field of protein design shifts emphasis from reproducing native-like protein structure to function, it has become important that these models treat the notion of specificity in molecular interactions. Although specificity may, in some cases, be achieved by optimization of a desired protein in isolation, methods have been developed to address directly the desire for proteins that exhibit specific functions and interactions.

Site-specific homing endonucleases are capable of inducing gene conversion via homologous recombination. Reprogramming their cleavage specificities allows the targeting of specific biological sites for gene correction or conversion. We used computational protein design to alter the cleavage specificity of I-MsoI for three contiguous base pair substitutions, resulting in an endonuclease whose activity and specificity for its new site rival that of wild-type I-MsoI for the original site. Concerted design for all simultaneous substitutions was more successful than a modular approach against individual substitutions, highlighting the importance of context-dependent redesign and optimization of protein-DNA interactions. We then used computational design based on the crystal structure of the designed complex, which revealed significant unanticipated shifts in DNA conformation, to create an endonuclease that specifically cleaves a site with four contiguous base pair substitutions. Our results demonstrate that specificity switches for multiple concerted base pair substitutions can be computationally designed, and that iteration between design and structure determination provides a route to large scale reprogramming of specificity.

Enzymes use substrate-binding energy both to promote ground-state association and to stabilize the reaction transition state selectively. The monomeric homing endonuclease I-AniI cleaves with high sequence specificity in the centre of a 20-base-pair (bp) DNA target site, with the amino (N)-terminal domain of the enzyme making extensive binding interactions with the left (-) side of the target site and the similarly structured carboxy (C)-terminal domain interacting with the right (+) side. Here we show that, despite the approximate twofold symmetry of the enzyme-DNA complex, there is almost complete segregation of interactions responsible for substrate binding to the (-) side of the interface and interactions responsible for transition-state stabilization to the (+) side. Although single base-pair substitutions throughout the entire DNA target site reduce catalytic efficiency, mutations in the (-) DNA half-site almost exclusively increase the dissociation constant (K(D)) and the Michaelis constant under single-turnover conditions (K(M)*), and those in the (+) half-site primarily decrease the turnover number (k(cat)*). The reduction of activity produced by mutations on the (-) side, but not mutations on the (+) side, can be suppressed by tethering the substrate to the endonuclease displayed on the surface of yeast. This dramatic asymmetry in the use of enzyme-substrate binding energy for catalysis has direct relevance to the redesign of endonucleases to cleave genomic target sites for gene therapy and other applications. Computationally redesigned enzymes that achieve new specificities on the (-) side do so by modulating K(M)*, whereas redesigns with altered specificities on the (+) side modulate k(cat)*. Our results illustrate how classical enzymology and modern protein design can each inform the other.

Computational protein design relies on a number of approximations to efficiently search the huge sequence space available to proteins. The fixed backbone and rotamer approximations in particular are important for formulating protein design as a discrete combinatorial optimization problem. However, the resulting coarse-grained sampling of possible side-chain terminal positions is problematic for the design of protein function, which depends on precise positioning of side-chain atoms. Although backbone flexibility can greatly increase the conformation freedom of side-chain functional groups, it is not obvious which backbone movements will generate the critical constellation of atoms responsible for protein function. Here, we report an automated method for identifying protein backbone movements that can give rise to any specified set of desired side-chain atomic placements and interactions, using protein-DNA interfaces as a model system. We use a library of previously observed protein-DNA interactions (motifs) and a rotamer-based description of side-chain conformation freedom to identify placements for the protein backbone that can give rise to a favorable side-chain interaction with DNA. We describe a tree-search algorithm for identifying those combinations of interactions from the library that can be realized with minimal perturbation of the protein backbone. We compare the efficiency of this method with the alternative approach of building and screening alternate backbone conformations.

Protein tyrosine kinases are critical cell signaling enzymes. These enzymes have a highly conserved Arg residue in their catalytic loop which is present two residues or four residues downstream from an absolutely conserved Asp catalytic base. Prior studies on protein tyrosine kinases Csk and Src revealed the potential for chemical rescue of catalytically deficient mutant kinases (Arg to Ala mutations) by small diamino compounds, particularly imidazole; however, the potency and efficiency of rescue was greater for Src. This current study further examines the structural and kinetic basis of rescue for mutant Src as compared to mutant Abl tyrosine kinase. An X-ray crystal structure of R388A Src revealed the surprising finding that a histidine residue of the N-terminus of a symmetry-related kinase inserts into the active site of the adjacent Src and mimics the hydrogen-bonding pattern seen in wild-type protein tyrosine kinases. Abl R367A shows potent and efficient rescue more comparable to Src, even though its catalytic loop is more like that of Csk. Various enzyme redesigns of the active sites indicate that the degree and specificity of rescue are somewhat flexible, but the overall properties of the enzymes and rescue agents play an overarching role. The newly discovered rescue agent 2-aminoimidazole is about as efficient as imidazole in rescuing R/A Src and Abl. Rate vs pH studies with these imidazole analogues suggest that the protonated imidazolium is the preferred form for chemical rescue, consistent with structural models. The efficient rescue seen with mutant Abl points to the potential of this approach to be used effectively to analyze Abl phosphorylation pathways in cells.

We recently reported that certain mutations in the FK506-rapamycin binding (FRB) domain disrupt its stability in vitro and in vivo (Stankunas et al. Mol. Cell, 2003, 12, 1615). To determine the precise residues that cause instability, we calculated the folding free energy (Delta G) of a collection of FRB mutants by measuring their intrinsic tryptophan fluorescence during reversible chaotropic denaturation. Our results implicate the T2098L point mutation as a key determinant of instability. Further, we found that some of the mutants in this collection were destabilized by up to 6 kcal mol(-1) relative to the wild type. To investigate how these mutants behave in cells, we expressed firefly luciferase fused to FRB mutants in African green monkey kidney (COS) cell lines and mouse embryonic fibroblasts (MEFs). When unstable FRB mutants were used, we found that the protein levels and the luminescence intensities were low. However, addition of a chemical ligand for FRB, rapamycin, restored luciferase activity. Interestingly, we found a roughly linear relationship between the Delta G of the FRB mutants calculated in vitro and the relative chemical rescue in cells. Because rapamycin is capable of simultaneously binding both FRB and the chaperone, FK506-binding protein (FKBP), we next examined whether FKBP might contribute to the protection of FRB mutants. Using both in vitro experiments and a cell-based model, we found that FKBP stabilizes the mutants. These findings are consistent with recent models that suggest damage to intrinsic Delta G can be corrected by pharmacological chaperones. Further, these results provide a collection of conditionally stable fusion partners for use in controlling protein stability.

Recent efforts to design de novo or redesign the sequence and structure of proteins using computational techniques have met with significant success. Most, if not all, of these computational methodologies attempt to model atomic-level interactions, and hence high-resolution structural characterization of the designed proteins is critical for evaluating the atomic-level accuracy of the underlying design force-fields. We previously used our computational protein design protocol RosettaDesign to completely redesign the sequence of the activation domain of human procarboxypeptidase A2. With 68% of the wild-type sequence changed, the designed protein, AYEdesign, is over 10 kcal/mol more stable than the wild-type protein. Here, we describe the high-resolution crystal structure and solution NMR structure of AYEdesign, which show that the experimentally determined backbone and side-chains conformations are effectively superimposable with the computational model at atomic resolution. To isolate the origins of the remarkable stabilization, we have designed and characterized a new series of procarboxypeptidase mutants that gain significant thermodynamic stability with a minimal number of mutations; one mutant gains more than 5 kcal/mol of stability over the wild-type protein with only four amino acid changes. We explore the relationship between force-field smoothing and conformational sampling by comparing the experimentally determined free energies of the overall design and these focused subsets of mutations to those predicted using modified force-fields, and both fixed and flexible backbone sampling protocols.

The reprogramming of DNA-binding specificity is an important challenge for computational protein design that tests current understanding of protein-DNA recognition, and has considerable practical relevance for biotechnology and medicine. Here we describe the computational redesign of the cleavage specificity of the intron-encoded homing endonuclease I-MsoI using a physically realistic atomic-level forcefield. Using an in silico screen, we identified single base-pair substitutions predicted to disrupt binding by the wild-type enzyme, and then optimized the identities and conformations of clusters of amino acids around each of these unfavourable substitutions using Monte Carlo sampling. A redesigned enzyme that was predicted to display altered target site specificity, while maintaining wild-type binding affinity, was experimentally characterized. The redesigned enzyme binds and cleaves the redesigned recognition site approximately 10,000 times more effectively than does the wild-type enzyme, with a level of target discrimination comparable to the original endonuclease. Determination of the structure of the redesigned nuclease-recognition site complex by X-ray crystallography confirms the accuracy of the computationally predicted interface. These results suggest that computational protein design methods can have an important role in the creation of novel highly specific endonucleases for gene therapy and other applications.

Protein-DNA interactions play a central role in transcriptional regulation and other biological processes. Investigating the mechanism of binding affinity and specificity in protein-DNA complexes is thus an important goal. Here we develop a simple physical energy function, which uses electrostatics, solvation, hydrogen bonds and atom-packing terms to model direct readout and sequence-specific DNA conformational energy to model indirect readout of DNA sequence by the bound protein. The predictive capability of the model is tested against another model based only on the knowledge of the consensus sequence and the number of contacts between amino acids and DNA bases. Both models are used to carry out predictions of protein-DNA binding affinities which are then compared with experimental measurements. The nearly additive nature of protein-DNA interaction energies in our model allows us to construct position-specific weight matrices by computing base pair probabilities independently for each position in the binding site. Our approach is less data intensive than knowledge-based models of protein-DNA interactions, and is not limited to any specific family of transcription factors. However, native structures of protein-DNA complexes or their close homologs are required as input to the model. Use of homology modeling can significantly increase the extent of our approach, making it a useful tool for studying regulatory pathways in many organisms and cell types.

Protein-DNA interactions are crucial for many biological processes. Attempts to model these interactions have generally taken the form of amino acid-base recognition codes or purely sequence-based profile methods, which depend on the availability of extensive sequence and structural information for specific structural families, neglect side-chain conformational variability, and lack generality beyond the structural family used to train the model. Here, we take advantage of recent advances in rotamer-based protein design and the large number of structurally characterized protein-DNA complexes to develop and parameterize a simple physical model for protein-DNA interactions. The model shows considerable promise for redesigning amino acids at protein-DNA interfaces, as design calculations recover the amino acid residue identities and conformations at these interfaces with accuracies comparable to sequence recovery in globular proteins. The model shows promise also for predicting DNA-binding specificity for fixed protein sequences: native DNA sequences are selected correctly from pools of competing DNA substrates; however, incorporation of backbone movement will likely be required to improve performance in homology modeling applications. Interestingly, optimization of zinc finger protein amino acid sequences for high-affinity binding to specific DNA sequences results in proteins with little or no predicted specificity, suggesting that naturally occurring DNA-binding proteins are optimized for specificity rather than affinity. When combined with algorithms that optimize specificity directly, the simple computational model developed here should be useful for the engineering of proteins with novel DNA-binding specificities.

Specific protein-protein interactions are crucial in signaling networks and for the assembly of multi-protein complexes, and represent a challenging goal for protein design. Optimizing interaction specificity requires both positive design, the stabilization of a desired interaction, and negative design, the destabilization of undesired interactions. Currently, no automated protein-design algorithms use explicit negative design to guide a sequence search. We describe a multi-state framework for engineering specificity that selects sequences maximizing the transfer free energy of a protein from a target conformation to a set of undesired competitor conformations. To test the multi-state framework, we engineered coiled-coil interfaces that direct the formation of either homodimers or heterodimers. The algorithm identified three specificity motifs that have not been observed in naturally occurring coiled coils. In all cases, experimental results confirm the predicted specificities.

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 self-energies of the polar side-chains 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 root-mean-square 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 side-chain relaxation (including >2 x 10(23) rotamer conformations of the protein) provide pKa values of unprecedented accuracy.