The relationship between the solvent-accessible surface area of a protein and the number of native contacts in its structure

The degree of protein folding is characterized either by the solvent-accessible surface area (S _ASA ) or the actual number of native contacts (N _cont). Evidently, these values should correlate with each other, as a decrease in S _ASA caused by the change in protein conformation during its folding must be accompanied by the corresponding increase in the number of native contacts. It is shown that this correlation does exist and is very strong (the correlation coefficient exceeds 99%), which can be used for an accurate and rapid estimation of the protein surface area from the number of native contacts. Among the methods commonly used for calculating the native contacts, the atom-atom approach gives the best fit if hydrogen atoms are taken into account and the cutoff value for the distance between the centers of atoms is taken to be 8 Å. The latter means that two layers of surface atoms are required to shield the protein core from the solvent.     

H-bonding in protein hydration revisited

H-bonding between protein surface polar/charged groups and water is one of the key factors of protein hydration. Here, we introduce an Accessible Surface Area (ASA) model for computationally efficient estimation of a free energy of water-protein H-bonding at any given protein conformation. The free energy of water-protein H-bonds is estimated using empirical formulas describing probabilities of hydrogen bond formation that were derived from molecular dynamics simulations of water molecules at the surface of a small protein, Crambin, from the Abyssinian cabbage (Crambe abyssinica) seed. The results suggest that atomic solvation parameters (ASP) widely used in continuum hydration models might be dependent on ASA for polar/charged atoms under consideration. The predictions of the model are found to be in qualitative agreement with the available experimental data on model compounds. This model combines the computational speed of ASA potential, with the high resolution of more sophisticated solvation methods.     

Energy functions for protein design I: efficient and accurate continuum electrostatics and solvation

Electrostatics and solvation energies are important for defining protein stability, structural specificity, and molecular recognition. Because these energies are difficult to compute quickly and accurately, they are often ignored or modeled very crudely in computational protein design. To address this problem, we have developed a simple, fast, and accurate approximation for calculating Born radii in the context of protein design calculations. When these approximate Born radii are used with the generalized Born continuum dielectric model, energies calculated by the 10(6)-fold slower finite difference Poisson-Boltzmann model are faithfully reproduced. A similar approach can be used for estimating solvent-accessible surface areas (SASAs). As an independent test, we show that these approximations can be used to accurately predict the experimentally determined pK(a)s of >200 ionizable groups from 15 proteins.     

CIRSE: a solvation energy estimator compatible with flexible protein docking and design applications

We present the Coordinate Internal Representation of Solvation Energy (CIRSE) for computing the solvation energy of protein configurations in terms of pairwise interactions between their atoms with analytic derivatives. Currently, CIRSE is trained to a Poisson/surface-area benchmark, but CIRSE is not meant to fit this benchmark exclusively. CIRSE predicts the overall solvation energy of protein structures from 331 NMR ensembles with 0.951+/-0.047 correlation and predicts relative solvation energy changes between members of individual ensembles with an accuracy of 15.8+/-9.6 kcal/mol. The energy of individual atoms in any of CIRSE's 17 types is predicted with at least 0.98 correlation. We apply the model in energy minimization, rotamer optimization, protein design, and protein docking applications. The CIRSE model shows some propensity to accumulate errors in energy minimization as well as rotamer optimization, but these errors are consistent enough that CIRSE correctly identifies the relative solvation energies of designed sequences as well as putative docked complexes. We analyze the errors accumulated by the CIRSE model during each type of simulation and suggest means of improving the model to be generally useful for all-atom simulations.     

Designed protein G core variants fold to native-like structures: sequence selection by ORBIT tolerates variation in backbone specification

The solution structures of two computationally designed core variants of the beta 1 domain of streptococcal protein G (G beta 1) were solved by (1)H NMR methods to assess the robustness of amino acid sequence selection by the ORBIT protein design package under changes in protein backbone specification. One variant has mutations at three of 10 core positions and corresponds to minimal perturbations of the native G beta 1 backbone. The other, with mutations at six of 10 positions, was calculated for a backbone in which the separation between G beta 1's alpha-helix and beta-sheet was increased by 15% relative to native G beta 1. Exchange broadening of some resonances and the complete absence of others in spectra of the sixfold mutant bespeak conformational heterogeneity in this protein. The NMR data were sufficiently abundant, however, to generate structures of similar, moderately high quality for both variants. Both proteins adopt backbone structures similar to their target folds. Moreover, the sequence selection algorithm successfully predicted all core chi(1) angles in both variants, five of six chi(2) angles in the threefold mutant and four of seven chi(2) angles in the sixfold mutant. We conclude that ORBIT calculates sequences that fold specifically to a geometry close to the template, even when the template is moderately perturbed relative to a naturally occurring structure. There are apparently limits to the size of acceptable perturbations: In this study, the larger perturbation led to undesired dynamic behavior.     

Backrub-like backbone simulation recapitulates natural protein conformational variability and improves mutant side-chain prediction

Incorporation of effective backbone sampling into protein simulation and design is an important step in increasing the accuracy of computational protein modeling. Recent analysis of high-resolution crystal structures has suggested a new model, termed backrub, to describe localized, hinge-like alternative backbone and side-chain conformations observed in the crystal lattice. The model involves internal backbone rotations about axes between C-alpha atoms. Based on this observation, we have implemented a backrub-inspired sampling method in the Rosetta structure prediction and design program. We evaluate this model of backbone flexibility using three different tests. First, we show that Rosetta backrub simulations recapitulate the correlation between backbone and side-chain conformations in the high-resolution crystal structures upon which the model was based. As a second test of backrub sampling, we show that backbone flexibility improves the accuracy of predicting point-mutant side-chain conformations over fixed backbone rotameric sampling alone. Finally, we show that backrub sampling of triosephosphate isomerase loop 6 can capture the millisecond/microsecond oscillation between the open and closed states observed in solution. Our results suggest that backrub sampling captures a sizable fraction of localized conformational changes that occur in natural proteins. Application of this simple model of backbone motions may significantly improve both protein design and atomistic simulations of localized protein flexibility.     

Exploring protein interfaces with a general photochemical reagent

Protein folding, natural conformational changes, or interaction between partners involved in recognition phenomena brings about differences in the solvent-accessible surface area (SASA) of the polypeptide chain. This primary event can be monitored by the differential chemical reactivity of functional groups along the protein sequence. Diazirine (DZN), a photoreactive gas similar in size to water, generates methylene carbene (:CH(2)). The extreme chemical reactivity of this species allows the almost instantaneous and indiscriminate modification of its immediate molecular cage. (3)H-DZN was successfully used in our laboratory for studying protein structure and folding. Here we address for the first time the usefulness of this probe to examine the area of interaction in protein-protein complexes. For this purpose we chose the complex formed between hen egg white lysozyme (HEWL) and the monoclonal antibody IgG(1) D1.3. :CH(2) labeling of free HEWL or complexed with IgG(1) D1.3 yields 2.76 and 2.32 mmol CH(2) per mole protein at 1 mM DZN concentration, respectively. This reduction (15%) becomes consistent with the expected decrement in the SASA of HEWL occurring upon complexation derived from crystallographic data (11%), in agreement with the known unspecific surface labeling reaction of :CH(2). Further comparative analysis at the level of tryptic peptides led to the identification of the sites involved in the interaction. Remarkably, those peptides implicated in the contact area show the highest differential labeling: H(15)GLDNYR(21), G(117)TDVQAWIR(125), andG(22)YSLGNWVCAAK(33). Thus, protein footprinting with DZN emerges as a feasible methodology useful for mapping contact regions of protein domains involved in macromolecular assemblies.     

Denaturant m values and heat capacity changes: relation to changes in accessible surface areas of protein unfolding.

Denaturant m values, the dependence of the free energy of unfolding on denaturant concentration, have been collected for a large set of proteins. The m value correlates very strongly with the amount of protein surface exposed to solvent upon unfolding, with linear correlation coefficients of R = 0.84 for urea and R = 0.87 for guanidine hydrochloride. These correlations improve to R = 0.90 when the effect of disulfide bonds on the accessible area of the unfolded protein is included. A similar dependence on accessible surface area has been found previously for the heat capacity change (delta Cp), which is confirmed here for our set of proteins. Denaturant m values and heat capacity changes also correlate well with each other. For proteins that undergo a simple two-state unfolding mechanism, the amount of surface exposed to solvent upon unfolding is a main structural determinant for both m values and delta Cp.     

Exploring the conformational space of protein side chains using dead-end elimination and the A* algorithm

We describe an algorithm which enables us to search the conformational space of the side chains of a protein to identify the global minimum energy combination of side chain conformations as well as all other conformations within a specified energy cutoff of the global energy minimum. The program is used to explore the side chain conformational energy surface of a number of proteins, to investigate how this surface varies with the energy model used to describe the interactions within the system and the rotamer library. Enumeration of the rotamer combinations enables us to directly evaluate the partition function, and thus calculate the side chain contribution to the conformational entropy of the folded protein. An investigation of these conformations and the relationships between them shows that most of the conformations near to the global energy minimum arise from changes in side chain conformations that are essentially independent; very few result from a concerted change in conformation of two or more residues. Some of the limitations of the approach are discussed. Proteins 33:227–239, 1998. © 1998 Wiley-Liss, Inc.     

Computational redesign of the SHV-1 beta-lactamase/beta-lactamase inhibitor protein interface

β-lactamases are enzymes that catalyze the hydrolysis of β-lactam antibiotics. β-lactamase/β-lactamase inhibitor protein (BLIP) complexes are emerging as a well characterized experimental model system for studying protein-protein interactions. BLIP is a 165 amino acid protein that inhibits several class A β-lactamases with a wide range of affinities: picomolar affinity for K1; nanomolar affinity for TEM-1, SME-1, and BlaI; but only micromolar affinity for SHV-1 β-lactamase. The large differences in affinity coupled with the availability of extensive mutagenesis data and high-resolution crystal structures for the TEM-1/BLIP and SHV-1/BLIP complexes make them attractive systems for the further development of computational design methodology. We used EGAD, a physics-based computational design program, to redesign BLIP in an attempt to increase affinity for SHV-1. Characterization of several of designs and point mutants revealed that in all cases, the mutations stabilize the interface by 10- to 1000-fold relative to wild type BLIP. The calculated changes in binding affinity for the mutants were within a mean absolute error of 0.87 kcal/mol from the experimental values, and comparison of the calculated and experimental values for a set of 30 SHV-1/BLIP complexes yielded a correlation coefficient of 0.77. Structures of the two complexes with the highest affinity, SHV-1/BLIP (E73M) and SHV-1/BLIP (E73M, S130K, S146M), are presented at 1.7 Å resolution. While the predicted structures have much in common with the experimentally determined structures, they do not coincide perfectly; in particular a salt bridge between SHV-1 D104 and BLIP K74 is observed in the experimental structures, but not in the predicted design conformations. This discrepancy highlights the difficulty of modeling salt bridge interactions with a protein design algorithm that approximates side chains as discrete rotamers. Nevertheless, while local structural features of the interface were sometimes miscalculated, EGAD is globally successful in designing complexes with increased affinity.     

Folding landscapes of the Alzheimer amyloid-beta(12-28) peptide

The energy landscape for folding of the 12-28 fragment of the Alzheimer amyloid beta (Abeta) peptide is characterized using replica-exchange molecular dynamics simulations with an all-atom peptide model and explicit solvent. At physiological temperatures, the peptide exists mostly as a collapsed random coil, populating a small fraction (less than 10%) of hairpins with a beta-turn at position V18F19, with another 10% of hairpin-like conformations possessing a bend rather than a turn in the central VFFA positions. A small fraction of the populated states, approximately 14%, adopt polyproline II (PPII) conformations. Folding of the structured hairpin states proceeds through the assembly of two locally stable segments, VFFAE and EDVGS. The interactions stabilizing these locally folded structural motifs are in conflict with those stabilizing the global fold of A12-28, a signature of underlying residual frustration in this peptide. At increased temperature, the population of both beta-strand and PPII conformations diminishes in favor of beta-turn and random-coil states. On the basis of the conformational preferences of Abeta 12-28 monomers, two models for the molecular structure of amyloid fibrils formed by this peptide are proposed.     

Protein fluorescence quenching by small molecules: protein penetration versus solvent exposure

Experiments were done to test the thesis that acrylamide and similar small molecules can penetrate into proteins on a nanosecond time scale. The approach taken was to measure the pattern of fluorescence quenching exhibited by quenching molecules differing in molecular character (size, polarity, charge) when these are directed against protein tryptophans that cover the whole range of tryptophan accessibility. If quenching involves protein penetration and internal quencher migration, one expects that larger quenchers and more polar quenchers should display lesser quenching. In fact, no significant dependence on quencher character was found. For proteins that display measurable quenching, the disparate quenchers studied display very similar quenching rate constants when directed against any particular protein tryptophan. For several proteins having tryptophans known to be buried, no quenching occurs. These results are not consistent with the view that the kinds of small molecules studied can quite generally penetrate into and diffuse about within proteins at near-diffusion-limited rates. Rather the results suggest that when quenching is observed, the pathway involves encounters with tryptophans that are partially exposed at the protein surface. Available crystallographic results support this conclusion.     

Physical scoring function based on AMBER force field and Poisson-Boltzmann implicit solvent for protein structure prediction

A well-behaved physics-based all-atom scoring function for protein structure prediction is analyzed with several widely used all-atom decoy sets. The scoring function, termed AMBER/Poisson-Boltzmann (PB), is based on a refined AMBER force field for intramolecular interactions and an efficient PB model for solvation interactions. Testing on the chosen decoy sets shows that the scoring function, which is designed to consider detailed chemical environments, is able to consistently discriminate all 62 native crystal structures after considering the heteroatom groups, disulfide bonds, and crystal packing effects that are not included in the decoy structures. When NMR structures are considered in the testing, the scoring function is able to discriminate 8 out of 10 targets. In the more challenging test of selecting near-native structures, the scoring function also performs very well: for the majority of the targets studied, the scoring function is able to select decoys that are close to the corresponding native structures as evaluated by ranking numbers and backbone Calpha root mean square deviations. Various important components of the scoring function are also studied to understand their discriminative contributions toward the rankings of native and near-native structures. It is found that neither the nonpolar solvation energy as modeled by the surface area model nor a higher protein dielectric constant improves its discriminative power. The terms remaining to be improved are related to 1-4 interactions. The most troublesome term is found to be the large and highly fluctuating 1-4 electrostatics term, not the dihedral-angle term. These data support ongoing efforts in the community to develop protein structure prediction methods with physics-based potentials that are competitive with knowledge-based potentials.     

Evaluation of an inverse molecular design algorithm in a model binding site

Computational molecular design is a useful tool in modern drug discovery. Virtual screening is an approach that docks and then scores individual members of compound libraries. In contrast to this forward approach, inverse approaches construct compounds from fragments, such that the computed affinity, or a combination of relevant properties, is optimized. We have recently developed a new inverse approach to drug design based on the dead-end elimination and A* algorithms employing a physical potential function. This approach has been applied to combinatorially constructed libraries of small-molecule ligands to design high-affinity HIV-1 protease inhibitors (Altman et al., J Am Chem Soc 2008;130:6099-6013). Here we have evaluated the new method using the well-studied W191G mutant of cytochrome c peroxidase. This mutant possesses a charged binding pocket and has been used to evaluate other design approaches. The results show that overall the new inverse approach does an excellent job of separating binders from nonbinders. For a few individual cases, scoring inaccuracies led to false positives. The majority of these involve erroneous solvation energy estimation for charged amines, anilinium ions, and phenols, which has been observed previously for a variety of scoring algorithms. Interestingly, although inverse approaches are generally expected to identify some but not all binders in a library, due to limited conformational searching, these results show excellent coverage of the known binders while still showing strong discrimination of the nonbinders.     

Fast and accurate side-chain topology and energy refinement (FASTER) as a new method for protein structure optimization

We have developed an original method for global optimization of protein side-chain conformations, called the Fast and Accurate Side-Chain Topology and Energy Refinement (FASTER) method. The method operates by systematically overcoming local minima of increasing order. Comparison of the FASTER results with those of the dead-end elimination (DEE) algorithm showed that both methods produce nearly identical results, but the FASTER algorithm is 100-1000 times faster than the DEE method and scales in a stable and favorable way as a function of protein size. We also show that low-order local minima may be almost as accurate as the global minimum when evaluated against experimentally determined structures. In addition, the new algorithm provides significant information about the conformational flexibility of individual side-chains. We observed that strictly rigid side-chains are concentrated mainly in the core of the protein, whereas highly flexible side-chains are found almost exclusively among solvent-oriented residues.     

Evaluation of a fast implicit solvent model for molecular dynamics simulations.

A solvation term based on the solvent accessible surface area (SASA) is combined with the CHARMM polar hydrogen force field for the efficient simulation of peptides and small proteins in aqueous solution. Only two atomic solvation parameters are used: one is negative for favoring the direct solvation of polar groups and the other positive for taking into account the hydrophobic effect on apolar groups. To approximate the water screening effects on the intrasolute electrostatic interactions, a distance-dependent dielectric function is used and ionic side chains are neutralized. The use of an analytical approximation of the SASA renders the model extremely efficient (i.e., only about 50% slower than in vacuo simulations). The limitations and range of applicability of the SASA model are assessed by simulations of proteins and structured peptides. For the latter, the present study and results reported elsewhere show that with the SASA model it is possible to sample a significant amount of folding/unfolding transitions, which permit the study of the thermodynamics and kinetics of folding at an atomic level of detail.     

Estimating protein-ligand binding free energy: atomic solvation parameters for partition coefficient and solvation free energy calculation

Solvation energy calculation is one of the main difficulties for the estimation of protein-ligand binding free energy and the correct scoring in docking studies. We have developed a new solvation energy estimation method for protein-ligand binding based on atomic solvation parameter (ASP), which has been shown to improve the power of protein-ligand binding free energy predictions. The ASP set, designed to handle both proteins and organic compounds and derived from experimental n-octanol/water partition coefficient (log P) data, contains 100 atom types (united model that treats hydrogen atoms implicitly) or 119 atom types (all-atom model that treats hydrogen atoms explicitly). By using this unified ASP set, an algorithm was developed for solvation energy calculation and was further integrated into a score function for predicting protein-ligand binding affinity. The score function reproduced the absolute binding free energies of a test set of 50 protein-ligand complexes with a standard error of 8.31 kJ/mol. As a byproduct, a conformation-dependent log P calculation algorithm named ASPLOGP was also implemented. The predictive results of ASPLOGP for a test set of 138 compounds were r = 0.968, s = 0.344 for the all-atom model and r = 0.962, s = 0.367 for the united model, which were better than previous conformation-dependent approaches and comparable to fragmental and atom-based methods. ASPLOGP also gave good predictive results for small peptides. The score function based on the ASP model can be applied widely in protein-ligand interaction studies and structure-based drug design.