Energy functions for protein design: adjustment with protein-protein complex affinities, models for the unfolded state, and negative design of solubility and specificity

The development of the EGAD program and energy function for protein design is described. In contrast to most protein design methods, which require several empirical parameters or heuristics such as patterning of residues or rotamers, EGAD has a minimalist philosophy; it uses very few empirical factors to account for inaccuracies resulting from the use of fixed backbones and discrete rotamers in protein design calculations, and describes the unfolded state, aggregates, and alternative conformers explicitly with physical models instead of fitted parameters. This approach unveils important issues in protein design that are often camouflaged by heuristic-emphasizing methods. Inter-atom energies are modeled with the OPLS-AA all-atom forcefield, electrostatics with the generalized Born continuum model, and the hydrophobic effect with a solvent-accessible surface area-dependent term. Experimental characterization of proteins designed with an unmodified version of the energy function revealed problems with under-packing, stability, aggregation, and structural specificity. Under-packing was addressed by modifying the van der Waals function. By optimizing only three parameters, the effects of >400 mutations on protein-protein complex formation were predicted to within 1.0 kcal mol(-1). As an independent test, this modified energy function was used to predict the stabilities of >1500 mutants to within 1.0 kcal mol(-1); this required a physical model of the unfolded state that includes more interactions than traditional tripeptide-based models. Solubility and structural specificity were addressed with simple physical approximations of aggregation and conformational equilibria. The complete energy function can design protein sequences that have high levels of identity with their natural counterparts, and have predicted structural properties more consistent with soluble and uniquely folded proteins than the initial designs.     

Chemically accurate protein structures: Validation of protein NMR structures by comparison of measured and predicted pK a values

A new method is presented for evaluating the quality of protein structures obtained by NMR. This method exploits the dependence between measurable chemical properties of a protein, namely pK _a values of acidic residues, and protein structure. The accurate and fast empirical computational method employed by the PROPKA program () allows the user to test the ability of a given structure to reproduce known pK _a values, which in turn can be used as a criterion for the selection of more accurate structures. We demonstrate the feasibility of this novel idea for a series of proteins for which both␣NMR and X-ray structures, as well as pK _a values of all ionizable residues, have been determined. For the 17 NMR ensembles used in this study, this criterion is shown effective in the elimination of a large number of NMR structure ensemble members.     

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.     

Electrostatic effects in unfolded staphylococcal nuclease

Structure-based calculations of pKa values and electrostatic free energies of proteins assume that electrostatic effects in the unfolded state are negligible. In light of experimental evidence showing that this assumption is invalid for many proteins, and with increasing awareness that the unfolded state is more structured and compact than previously thought, a detailed examination of electrostatic effects in unfolded proteins is warranted. Here we address this issue with structure-based calculations of electrostatic interactions in unfolded staphylococcal nuclease. The approach involves the generation of ensembles of structures representing the unfolded state, and calculation of Coulomb energies to Boltzmann weight the unfolded state ensembles. Four different structural models of the unfolded state were tested. Experimental proton binding data measured with a variant of nuclease that is unfolded under native conditions were used to establish the validity of the calculations. These calculations suggest that weak Coulomb interactions are an unavoidable property of unfolded proteins. At neutral pH, the interactions are too weak to organize the unfolded state; however, at extreme pH values, where the protein has a significant net charge, the combined action of a large number of weak repulsive interactions can lead to the expansion of the unfolded state. The calculated pKa values of ionizable groups in the unfolded state are similar but not identical to the values in small peptides in water. These studies suggest that the accuracy of structure-based calculations of electrostatic contributions to stability cannot be improved unless electrostatic effects in the unfolded state are calculated explicitly.     

Protonation state of Asp30 exerts crucial influence over surface loop rearrangements responsible for NO release in nitrophorin 4

pK(a) values of ionizable residues were calculated for the crystal structures describing the pH and NO binding dependant conformations of nitrophorin 4, a pH sensitive NO carrier heme protein. Comparison of resultant H-bonding patterns allowed the identification of the amino acids that take part in signaling pH change. We carried out MD simulations to show that the protonation state of Asp30, buried in the closed conformation, is crucial for maintaining the tight packed conformation of the closed form of the complex - presenting a model for the functional decrease of NO binding affinity of nitrophorins at physiological pH.     

Structural origins of high apparent dielectric constants experienced by ionizable groups in the hydrophobic core of a protein

The side chains of Lys66, Asp66, and Glu66 in staphylococcal nuclease are fully buried and surrounded mainly by hydrophobic matter, except for internal water molecules associated with carboxylic oxygen atoms. These ionizable side chains titrate with pK_a values of 5.7, 8.8, and 8.9, respectively. To reproduce these pK_a values with continuum electrostatics calculations, we treated the protein with high dielectric constants. We have examined the structural origins of these high apparent dielectric constants by using NMR spectroscopy to characterize the structural response to the ionization of these internal side chains. Substitution of Val66 with Lys66 and Asp66 led to increased conformational fluctuations of the microenvironments surrounding these groups, even under pH conditions where Lys66 and Asp66 are neutral. When Lys66, Asp66, and Glu66 are charged, the proteins remain almost fully folded, but resonances for a few backbone amides adjacent to the internal ionizable residues are broadened. This suggests that the ionization of the internal groups promotes a local increase in dynamics on the intermediate timescale, consistent with either partial unfolding or increased backbone fluctuations of helix 1 near residue 66, or, less likely, with increased fluctuations of the charged side chains at position 66. These experiments confirm that the high apparent dielectric constants reported by internal Lys66, Asp66, and Glu66 reflect localized changes in conformational fluctuations without incurring detectable global structural reorganization. To improve structure-based pK_a calculations in proteins, we will need to learn how to treat this coupling between ionization of internal groups and local changes in conformational fluctuations explicitly.     

One- and two-body decomposable Poisson-Boltzmann methods for protein design calculations

Successfully modeling electrostatic interactions is one of the key factors required for the computational design of proteins with desired physical, chemical, and biological properties. In this paper, we present formulations of the finite difference Poisson-Boltzmann (FDPB) model that are pairwise decomposable by side chain. These methods use reduced representations of the protein structure based on the backbone and one or two side chains in order to approximate the dielectric environment in and around the protein. For the desolvation of polar side chains, the two-body model has a 0.64 kcal/mol RMSD compared to FDPB calculations performed using the full representation of the protein structure. Screened Coulombic interaction energies between side chains are approximated with an RMSD of 0.13 kcal/mol. The methods presented here are compatible with the computational demands of protein design calculations and produce energies that are very similar to the results of traditional FDPB calculations.     

Are acidic and basic groups in buried proteins predicted to be ionized?

Ionizable residues play essential roles in proteins, modulating protein stability, fold and function. Asp, Glu, Arg, and Lys make up about a quarter of the residues in an average protein. Multi-conformation continuum electrostatic (MCCE) calculations were used to predict the ionization states of all acidic and basic residues in 490 proteins. Of all 36,192 ionizable residues, 93.5% were predicted to be ionized. Thirty-five percent have lost 4.08 kcal/mol solvation energy (DeltaDeltaG(rxn)) sufficient to shift a pK(a) by three pH units in the absence of other interactions and 17% have DeltaDeltaG(rxn) sufficient to shift pK(a) by five pH units. Overall 85% of these buried residues (DeltaDeltaG(rxn)>5DeltapK units) are ionized, including 92% of the Arg, 86% of the Asp, 77% of the Glu, and 75% of the Lys. Ion-pair interactions stabilize the ionization of both acids and bases. The backbone dipoles stabilize anions more than cations. The interactions with polar side-chains are also different for acids and bases. Asn and Gln stabilize all charges, Ser and Thr stabilize only acids while Tyr rarely stabilize Lys. Thus, hydroxyls are better hydrogen bond donors than acceptors. Buried ionized residues are more likely to be conserved than those on the surface. There are 3.95 residues buried per 100 residues in an average protein.     

Carbanions from decarboxylation of orotate analogs: stability and mechanistic implications

The pKa's of the 6-CH groups of 1,3-dimethyluracil, N-methyl-2-pyridone, and N-methyl-4-pyridone were determined through their reactions with bases derived from carbon acids with known pKa and the reactions of their corresponding carbanions with the carbon acids. No correlation between the stability of the carbanions and the rate of decarboxylation of corresponding carboxylic acids was found.     

Distance dependence and salt sensitivity of pairwise, coulombic interactions in a protein

Histidine pK(a) values were measured in charge-reversal (K78E, K97E, K127E, and K97E/K127E) and charge-neutralization (E10A, E101A, and R35A) mutants of staphylococcal nuclease (SNase) by (1)H-NMR spectroscopy. Energies of interaction between pairs of charges (DeltaG(ij)) were obtained from the shifts in pK(a) values relative to wild-type values. The data describe the distance dependence and salt sensitivity of pairwise coulombic interactions. Calculations with a continuum electrostatics method captured the experimental DeltaG(ij) when static structures were used and when the protein interior was treated empirically with a dielectric constant of 20. The DeltaG(ij) when r(ij) < or = 10 A were exaggerated slightly in the calculations. Coulomb's law with a dielectric constant near 80 and a Debye-Hückel term to account for screening by the ionic strength reproduced the salt sensitivity and distance dependence of DeltaG(ij) as well as the structure-based method. In their interactions with each other, surface charges behave as if immersed in water; the Debye length describes realistically the distance where interactions become negligible at a given ionic strength. On average, charges separated by distances (r(ij)) approximately 5 A interacted with DeltaG(ij) approximately 0.6 kcal/mole in 0.01 M KCl, but DeltaG(ij) decayed to < or =0.10 kcal/mole when r(ij) = 20 A. In 0.10 M KCl, DeltaG(ij) approximately 0.10 kcal/mole when r(ij) = 10 A. In 1.5 M KCl, only short-range interactions with r(ij) < or = 5 A persisted. Although at physiological ionic strengths the interactions between charges separated by more than 10 A are extremely weak, in situations where charge imbalance exists many weak interactions can cumulatively produce substantial effects.     

Structural and mechanistic insights into unusual thiol disulfide oxidoreductase

Cytoplasmic desulfothioredoxin (Dtrx) from the anaerobe Desulfovibrio vulgaris Hildenborough has been identified as a new member of the thiol disulfide oxidoreductase family. The active site of Dtrx contains a particular consensus sequence, CPHC, never seen in the cytoplasmic thioredoxins and generally found in periplasmic oxidases. Unlike canonical thioredoxins (Trx), Dtrx does not present any disulfide reductase activity, but it presents instead an unusual disulfide isomerase activity. We have used NMR spectroscopy to gain insights into the structure and the catalytic mechanism of this unusual Dtrx. The redox potential of Dtrx (-181 mV) is significantly less reducing than that of canonical Trx. A pH dependence study allowed the determination of the pK(a) of all protonable residues, including the cysteine and histidine residues. Thus, the pK(a) values for the thiol group of Cys(31) and Cys(34) are 4.8 and 11.3, respectively. The His(33) pK(a) value, experimentally determined for the first time, differs notably as a function of the redox states, 7.2 for the reduced state and 4.6 for the oxidized state. These data suggest an important role for His(33) in the molecular mechanism of Dtrx catalysis that is confirmed by the properties of mutant DtrxH33G protein. The NMR structure of Dtrx shows a different charge repartition compared with canonical Trx. The results presented are likely indicative of the involvement of this protein in the catalysis of substrates specific of the anaerobe cytoplasm of DvH. The study of Dtrx is an important step toward revealing the molecular details of the thiol-disulfide oxidoreductase catalytic mechanism.     

A conserved arginine plays a role in the catalytic cycle of the protein disulphide isomerases

The pK(a) values of the CXXC active-site cysteine residues play a critical role in determining the physiological function of the thioredoxin superfamily. To act as an efficient thiol-disulphide oxidant the thiolate state of the N-terminal cysteine must be stabilised and the thiolate state of the C-terminal cysteine residue destabilised. While increasing the pK(a) value of the C-terminal cysteine residue promotes oxidation of substrates, it has an inhibitory effect on the reoxidation of the enzyme, which is promoted by the formation of a thiolate at this position. Since reoxidation is essential to complete the catalytic cycle, the differential requirement for a high and a low pK(a) value for the C-terminal cysteine residue for different steps in the reaction presents us with a paradox. Here, we report the identification of a conserved arginine residue, located in the loop between beta5 and alpha4 of the catalytic domains of the human protein disulphide isomerase (PDI) family, which is critical for the catalytic function of PDI, ERp57, ERp72 and P5, specifically for reoxidation. An examination of the published NMR structure for the a domain of PDI combined with molecular dynamic studies suggest that the side-chain of this arginine residue moves into and out of the active-site locale and that this has a very marked effect on the pK(a) value of the active-site cysteine residues. This intra-domain motion resolves the apparent dichotomy of the pK(a) requirements for the C-terminal active-site cysteine.     

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.     

Roles of N-terminal pyroglutamate in maintaining structural integrity and pKa values of catalytic histidine residues in bullfrog ribonuclease 3

Many proteins and bioactive peptides contain an N-terminal pyroglutamate residue (Pyr1). This residue reduces the susceptibility of the protein to aminopeptidases and often has important functional roles. The antitumor ribonuclease RC-RNase 3 (RNase 3) from oocytes of Rana catesbeiana (bullfrog) is one such protein. We have produced recombinant RNase 3 containing the N-terminal Pyr1 (pRNase 3) and found it to be indistinguishable from the native RNase 3 by mass spectrometry and a variety of other biochemical and immunological criteria. We demonstrated by NMR analysis that the Pyr1 of pRNase 3 forms hydrogen bonds with Lys9 and Ile96 and stabilizes the N-terminal alpha-helix in a rigid conformation. In contrast, the N-terminal alpha-helix becomes flexible and the pKa values of the catalytic residues His10 and His97 altered when Pyr1 formation is blocked by an extra methionine at the N terminus in the recombinant mqRNase 3. Thus, our results provide a mechanistic explanation on the essential role of Pyr1 in maintaining the structural integrity, especially at the N-terminal alpha-helix, and in providing the proper environment for the ionization of His10 and His97 residues for catalysis and cytotoxicity against HeLa cells.     

Maturation efficiency, trypsin sensitivity, and optical properties of Arg96, Glu222, and Gly67 variants of green fluorescent protein

Spontaneous chromophore biosynthesis in green fluorescent protein (GFP) is initiated by a main-chain cyclization reaction catalyzed by the protein fold. To investigate the structural prerequisites for chromophore formation, we have substituted the conserved residues Arg96, Glu222, and Gly67. Upon purification, the variants can be ordered based on their decreasing extent of chromophore maturation according to the series EGFP, E222Q, R96K, G67A, and R96M. Arg96 and Glu222 appear to play catalytic roles, whereas Gly67 is likely important in interior packing to enforce correct hydrogen bonding to Arg96. The effect of Arg96 can be partially compensated for by a lysine, but not by a methionine residue, confirming its electrophilic role. Limited trypsinolysis data suggest that protein stability is largely unaffected by the presence of the chromophore, inconsistent with the mechanical compression hypothesis. Trends in optical properties may be related to the degree of chromophore charge delocalization, which is modulated by residue 96.     

Novel bisubstrate analog inhibitors of serotonin N-acetyltransferase: the importance of being neutral

Linker modified novel bisubstrate analog inhibitors 4-7 for serotonin N-acetyltransferase (arylalkylamine N-acetyltransferase, AANAT) have been designed and synthesized. Examination of these inhibitors with AANAT in vitro suggested that: (i) linker hydrogen bonding makes only modest contributions to the affinity of bisubstrate analog inhibitors studied; (ii) greater than or equal to four methylene groups between the indole and the coenzyme A (CoASH) moieties are required for a bisubstrate analog inhibitor to achieve strong AANAT inhibition; (iii) the AANAT active site appears not to accommodate positively charged linkers as well as neutral ones; and (iv) substrate amine pKa depression may constitute one strategy for AANAT substrate recognition and catalysis. The results reported here have enhanced our understanding of AANAT substrate recognition/catalysis, and are important for novel inhibitor design. Since AANAT belongs to the GCN5-related N-acetyltransferase (GNAT) superfamily, our experimental strategies should find applications for other acetyltransferases.     

Cholesterol-induced protein sorting: an analysis of energetic feasibility

The mechanism(s) underlying the sorting of integral membrane proteins between the Golgi complex and the plasma membrane remain uncertain because no specific Golgi retention signal has been found. Moreover one can alter a protein's eventual localization simply by altering the length of its transmembrane domain (TMD). M. S. Bretscher and S. Munro (SCIENCE: 261:1280-1281, 1993) therefore proposed a physical sorting mechanism based on the hydrophobic match between the proteins' TMD and the bilayer thickness, in which cholesterol would regulate protein sorting by increasing the lipid bilayer thickness. In this model, Golgi proteins with short TMDs would be excluded from cholesterol-enriched domains (lipid rafts) that are incorporated into transport vesicles destined for the plasma membrane. Although attractive, this model remains unproven. We therefore evaluated the energetic feasibility of a cholesterol-dependent sorting process using the theory of elastic liquid crystal deformations. We show that the distribution of proteins between cholesterol-enriched and cholesterol-poor bilayer domains can be regulated by cholesterol-induced changes in the bilayer physical properties. Changes in bilayer thickness per se, however, have only a modest effect on sorting; the major effect arises because cholesterol changes also the bilayer material properties, which augments the energetic penalty for incorporating short TMDs into cholesterol-enriched domains. We conclude that cholesterol-induced changes in the bilayer physical properties allow for effective and accurate sorting which will be important generally for protein partitioning between different membrane domains.     

Probing the sweet determinants of brazzein: wild-type brazzein and a tasteless variant, brazzein-ins(R18a-I18b), exhibit different pH-dependent NMR chemical shifts

Brazzein is a small, intensely sweet protein. As a probe of the functional properties of its solvent-exposed loop, two residues (Arg-Ile) were inserted between Leu18 and Ala19 of brazzein. Psychophysical testing demonstrated that this mutant is totally tasteless. NMR chemical shift mapping of differences between this mutant and brazzein indicated that residues affected by the insertion are localized to the mutated loop, the region of the single alpha-helix, and around the Cys16-Cys37 disulfide bond. Residues unaffected by this mutation included those near the C-terminus and in the loop connecting the alpha-helix and the second beta-strand. In particular, several residues of brazzein previously shown to be essential for its sweetness (His31, Arg33, Glu41, Arg43, Asp50, and Tyr54) exhibited negligible chemical shift changes. Moreover, the pH dependence of the chemical shifts of His31, Glu41, Asp50, and Tyr54 were unaltered by the insertion. The insertion led to large chemical shift and pKa perturbation of Glu36, a residue shown previously to be important for brazzein's sweetness. These results serve to refine the known sweetness determinants of brazzein and lend further support to the idea that the protein interacts with a sweet-taste receptor through a multi-site interaction mechanism, as has been postulated for brazzein and other sweet proteins (monellin and thaumatin).     

Disruption of an intermonomer salt bridge in the p53 tetramerization domain results in an increased propensity to form amyloid fibrils

We describe in molecular detail how disruption of an intermonomer salt bridge (Arg337-Asp352) leads to partial destabilization of the p53 tetramerization domain and a dramatically increased propensity to form amyloid fibrils. At pH 4.0 and 37 degrees C, a p53 tetramerization domain mutant (p53tet-R337H), associated with adrenocortical carcinoma in children, readily formed amyloid fibrils, while the wild-type (p53tet-wt) did not. We characterized these proteins by equilibrium denaturation, 13C(alpha) secondary chemical shifts, (1H)-15N heteronuclear NOEs, and H/D exchange. Although p53tet-R337H was thermodynamically less stable, NMR data indicated that the two proteins had similar secondary structure and molecular dynamics. NMR derived pK(a) values indicated that at low pH the R337H mutation partially disrupted an intermonomer salt bridge. Backbone H/D exchange results showed that for at least a small population of p53tet-R337H molecules disruption of this salt bridge resulted in partial destabilization of the protein. It is proposed that this decrease in p53tet-R337H stability resulted in an increased propensity to form amyloid fibrils.