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.     

Use of MM-PB/SA in estimating the free energies of proteins: application to native, intermediates, and unfolded villin headpiece

We investigated the stability of three different ensembles of the 36-mer villin headpiece subdomain, the native, a compact folding intermediate, and the random coil. Structures were taken from a 1-micros molecular dynamics folding simulation and a 100-ns control simulation on the native structure. Our approach for each conformation is to first determine the solute internal energy from the molecular mechanics potential and then to add the change resulting from solvation (DeltaG(solv)). Explicit water was used to run the simulation, and a continuum model was used to estimate DeltaG(solv) with the finite difference Poisson-Boltzmann model accounting for the polarization part and a linearly surface area-dependent term for the non-polar part. We leave out the solute vibrational entropy from these values but demonstrate that there is no statistical difference among the native, folding intermediate, and random coil ensembles. We find the native ensemble to be approximately 26 kcal/mol more stable than the folding intermediate and approximately 39 kcal/mol more stable than the random coil ensemble. With an experimental estimate for the free energy of denaturation equal to 3 kcal/mol, we approximate the non-native degeneracy to lie between 10(16) and 10.(25) We also present a possible scheme for the mechanism of folding, first-order exponential decay of a putative transition state, with an estimate for the t(1/2) of folding of approximately 1 micros.     

Structural bases of lectin-carbohydrate affinities: comparison with protein-folding energetics.

We have made a comparative structure based analysis of the thermodynamics of lectin-carbohydrate (L-C) binding and protein folding. Examination of the total change in accessible surface area in those processes revealed a much larger decrease in free energy per unit of area buried in the case of L-C associations. According to our analysis, this larger stabilization of L-C interactions arises from a more favorable enthalpy of burying a unit of polar surface area, and from higher proportions of polar areas. Hydrogen bonds present at 14 L-C interfaces were identified, and their overall characteristics were compared to those reported before for hydrogen bonds in protein structures. Three major factors might explain why polar-polar interactions are stronger in L-C binding than in protein folding: (1) higher surface density of hydrogen bonds; (2) better hydrogen-bonding geometry; (3) larger proportion of hydrogen bonds involving charged groups. Theoretically, the binding entropy can be partitioned into three main contributions: entropy changes due to surface desolvation, entropy losses arising from freezing rotatable bonds, and entropic effects that result from restricting translation and overall rotation motions. These contributions were estimated from structural information and added up to give calculated binding entropies. Good correlation between experimental and calculated values was observed when solvation effects were treated according to a parametrization developed by other authors from protein folding studies. Finally, our structural parametrization gave calculated free energies that deviate from experimental values by 1.1 kcal/mol on the average; this amounts to an uncertainty of one order of magnitude in the binding constant.     

Recent advances in the development and application of implicit solvent models in biomolecule simulations

Advances have recently been made in the development of implicit solvent methodologies and their application to the modeling of biomolecules, particularly with regard to generalized Born approaches, dielectric screening function formulations and models based on solvent-accessible surface areas. Interesting new developments include more refined non-polar solvation energy estimators, and implicit methods for modeling low-dielectric and heterogeneous environments such as membrane systems. These have been successfully applied to molecular dynamics simulations, the scoring of protein conformations, and the calculation of binding affinities and folding free energy landscapes.     

Computational analysis of the binding affinities of the natural-product cyclopentapeptides argifin and argadin to chitinase B from Serratia marcescens

Molecular dynamics (MD) simulations and the molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) method were applied to study the interaction of the natural-product cyclopentapeptide chitinase inhibitors argifin and argadin with chitinase B (ChiB) from Serratia marcescens. Argadin inhibited ChiB with an inhibition constant (K(i)) value of 20 nM, which was three orders of magnitude greater than that of argifin (K(i)=33,000 nM). The MM-PBSA free-energy analysis provided absolute binding free energies of -6.98 and -11.16 kcal/mol for the argifin and argadin complexes, respectively. These estimates were in good agreement with the free energies derived from the experimental K(i) values (-6.36 and -10.92 kcal/mol for the argifin and argadin complexes, respectively). The energetic analysis revealed that the van der Waals and nonpolar solvation energies drove the binding of both argifin and argadin. We found that the binding of argadin gained approximately 12 kcal/mol more van der Waals energy than that of argifin, which was mainly responsible for the difference in binding free energy between argifin and argadin. In particular, W220 and W403 of ChiB were found to contribute to the more favorable van der Waals interaction with argadin. We also designed argifin derivatives with better binding affinity, in which a constituent amino-acid residue of argifin was mutated to one with a bulky side chain. The derivative in which D-Ala of argifin was replaced with D-Trp appeared to possess a binding affinity that was equally potent to that of argadin.     

Protein surface roughness accounts for binding free energy of Plasmepsin II‐ligand complexes

The calculation of absolute binding affinities for protein‐inhibitor complexes remains as one of the main challenges in computational structure‐based ligand design. The present work explored the calculations of surface fractal dimension (as a measure of surface roughness) and the relationship with experimental binding free energies of Plasmepsin II complexes. Plasmepsin II is an attractive target for novel therapeutic compounds to treat malaria. However, the structural flexibility of this enzyme is a drawback when searching for specific inhibitors. Concerning that, we performed separate explicitly solvated molecular dynamics simulations using the available high‐resolution crystal structures of different Plasmepsin II complexes. Molecular dynamics simulations allowed a better approximation to systems dynamics and, therefore, a more reliable estimation of surface roughness. This constitutes a novel approximation in order to obtain more realistic values of fractal dimension, because previous works considered only x‐ray structures. Binding site fractal dimension was calculated considering the ensemble of structures generated at different simulation times. A linear relationship between binding site fractal dimension and experimental binding free energies of the complexes was observed within 20 ns. Previous studies of the subject did not uncover this relationship. Regression model, coined FD model, was built to estimate binding free energies from binding site fractal dimension values. Leave‐one‐out cross‐validation showed that our model reproduced accurately the absolute binding free energies for our training set (R² = 0.76; <|error|> =0.55 kcal/mol; SD_error = 0.19 kcal/mol). The fact that such a simple model may be applied raises some questions that are addressed in the article.     

Estimation of the binding affinities of FKBP12 inhibitors using a linear response method.

A series of non-immunosuppressive inhibitors of FK506 binding protein (FKBP12) are investigated using Monte Carlo statistical mechanics simulations. These small molecules may serve as scaffolds for chemical inducers of protein dimerization, and have recently been found to have FKBP12-dependent neurotrophic activity. A linear response model was developed for estimation of absolute binding free energies based on changes in electrostatic and van der Waals energies and solvent-accessible surface areas, which are accumulated during simulations of bound and unbound ligands. With average errors of 0.5 kcal/mol, this method provides a relatively rapid way to screen the binding of ligands while retaining the structural information content of more rigorous free energy calculations.     

Tryptophan-lipid interactions in membrane protein folding probed by ultraviolet resonance Raman and fluorescence spectroscopy

Aromatic amino acids of membrane proteins are enriched at the lipid-water interface. The role of tryptophan on the folding and stability of an integral membrane protein is investigated with ultraviolet resonance Raman and fluorescence spectroscopy. We investigate a model system, the β-barrel outer membrane protein A (OmpA), and focus on interfacial tryptophan residues oriented toward the lipid bilayer (trp-7, trp-170, or trp-15) or the interior of the β-barrel pore (trp-102). OmpA mutants with a single tryptophan residue at a nonnative position 170 (Trp-170) or a native position 7 (Trp-7) exhibit the greatest stability, with Gibbs free energies of unfolding in the absence of denaturant of 9.4 and 6.7 kcal/mol, respectively. These mutants are more stable than the tryptophan-free OmpA mutant, which exhibits a free energy of unfolding of 2.6 kcal/mol. Ultraviolet resonance Raman spectra of Trp-170 and Trp-7 reveal evolution of a hydrogen bond in a nonpolar environment during the folding reaction, evidenced by systematic shifts in hydrophobicity and hydrogen bond markers. These observations suggest that the hydrogen bond acceptor is the lipid acyl carbonyl group, and this interaction contributes significantly to membrane protein stabilization. Other spectral changes are observed for a tryptophan residue at position 15, and these modifications are attributed to development of a tryptophan-lipid cation-π interaction that is more stabilizing than an intraprotein hydrogen bond by ∼2 kcal/mol. As expected, there is no evidence for lipid-protein interactions for the tryptophan residue oriented toward the interior of the β-barrel pore. These results highlight the significance of lipid-protein interactions, and indicate that the bilayer provides more than a hydrophobic environment for membrane protein folding. Instead, a paradigm of lipid-assisted membrane protein folding and stabilization must be adopted.     

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.     

Thermodynamics of folding and binding in an affibody:affibody complex

Affibody binding proteins are selected from phage-displayed libraries of variants of the 58 residue Z domain. Z(Taq) is an affibody originally selected as a binder to Taq DNA polymerase. The anti-Z(Taq) affibody was selected as a binder to Z(Taq) and the Z(Taq):anti-Z(Taq) complex is formed with a dissociation constant K(d)=0.1 microM. We have determined the structure of the Z(Taq):anti-Z(Taq) complex as well as the free state structures of Z(Taq) and anti-Z(Taq) using NMR. Here we complement the structural data with thermodynamic studies of Z(Taq) and anti-Z(Taq) folding and complex formation. Both affibody proteins show cooperative two-state thermal denaturation at melting temperatures T(M) approximately 56 degrees C. Z(Taq):anti-Z(Taq) complex formation at 25 degrees C in 50 mM NaCl and 20 mM phosphate buffer (pH 6.4) is enthalpy driven with DeltaH degrees (bind) = -9.0 (+/-0.1) kcal mol(-1)(.) The heat capacity change DeltaC(P) degrees (,bind)=-0.43 (+/-0.01) kcal mol(-1) K(-1) is in accordance with the predominantly non-polar character of the binding surface, as judged from calculations based on changes in accessible surface areas. A further dissection of the small binding entropy at 25 degrees C (-TDeltaS degrees (bind) = -0.6 (+/-0.1) kcal mol(-1)) suggests that a favourable desolvation of non-polar surface is almost completely balanced by unfavourable conformational entropy changes and loss of rotational and translational entropy. Such effects can therefore be limiting for strong binding also when interacting protein components are stable and homogeneously folded. The combined structure and thermodynamics data suggest that protein properties are not likely to be a serious limitation for the development of engineered binding proteins based on the Z domain.     

Thermodynamics of amide hydrogen bond formation in polar and apolar solvents

We present the initial findings of a theoretical study of hydrogen bond formation between two formamide molecules in water and in carbon tetrachloride. These systems were chosen as the simplest models for secondary structure formation in the polar environment near the protein surface and the apolar environment of the protein interior. We have employed thermodynamic simulation methods to obtain absolute binding free energies and free energy profiles for the formation of peptide hydrogen bonds in the two solvents. We find that the amide hydrogen bond is stable by 8.4 kcal/mol in CCl4, and by 0.3 kcal/mol in water. Our results indicate also that the hydrogen-bonded dimer is 2.2 kcal/mol more stable in water than it is in CCl4. We compare our results with those from experiment, and discuss their use in interpreting mechanisms of protein folding.     

Comparison of binding energies of SrcSH2-phosphotyrosyl peptides with structure-based prediction using surface area based empirical parameterization

The prediction of binding energies from the three-dimensional (3D) structure of a protein-ligand complex is an important goal of biophysics and structural biology. Here, we critically assess the use of empirical, solvent-accessible surface area-based calculations for the prediction of the binding of Src-SH2 domain with a series of tyrosyl phosphopeptides based on the high-affinity ligand from the hamster middle T antigen (hmT), where the residue in the pY+ 3 position has been changed. Two other peptides based on the C-terminal regulatory site of the Src protein and the platelet-derived growth factor receptor (PDGFR) are also investigated. Here, we take into account the effects of proton linkage on binding, and test five different surface area-based models that include different treatments for the contributions to conformational change and protein solvation. These differences relate to the treatment of conformational flexibility in the peptide ligand and the inclusion of proximal ordered solvent molecules in the surface area calculations. This allowed the calculation of a range of thermodynamic state functions (deltaCp, deltaS, deltaH, and deltaG) directly from structure. Comparison with the experimentally derived data shows little agreement for the interaction of SrcSH2 domain and the range of tyrosyl phosphopeptides. Furthermore, the adoption of the different models to treat conformational change and solvation has a dramatic effect on the calculated thermodynamic functions, making the predicted binding energies highly model dependent. While empirical, solvent-accessible surface area based calculations are becoming widely adopted to interpret thermodynamic data, this study highlights potential problems with application and interpretation of this type of approach. There is undoubtedly some agreement between predicted and experimentally determined thermodynamic parameters: however, the tolerance of this approach is not sufficient to make it ubiquitously applicable.     

Insights into protein-protein binding by binding free energy calculation and free energy decomposition for the Ras-Raf and Ras-RalGDS complexes

Absolute binding free energy calculations and free energy decompositions are presented for the protein-protein complexes H-Ras/C-Raf1 and H-Ras/RalGDS. Ras is a central switch in the regulation of cell proliferation and differentiation. In our study, we investigate the capability of the molecular mechanics (MM)-generalized Born surface area (GBSA) approach to estimate absolute binding free energies for the protein-protein complexes. Averaging gas-phase energies, solvation free energies, and entropic contributions over snapshots extracted from trajectories of the unbound proteins and the complexes, calculated binding free energies (Ras-Raf: -15.0(+/-6.3)kcal mol(-1); Ras-RalGDS: -19.5(+/-5.9)kcal mol(-1)) are in fair agreement with experimentally determined values (-9.6 kcal mol(-1); -8.4 kcal mol(-1)), if appropriate ionic strength is taken into account. Structural determinants of the binding affinity of Ras-Raf and Ras-RalGDS are identified by means of free energy decomposition. For the first time, computationally inexpensive generalized Born (GB) calculations are applied in this context to partition solvation free energies along with gas-phase energies between residues of both binding partners. For selected residues, in addition, entropic contributions are estimated by classical statistical mechanics. Comparison of the decomposition results with experimentally determined binding free energy differences for alanine mutants of interface residues yielded correlations with r(2)=0.55 and 0.46 for Ras-Raf and Ras-RalGDS, respectively. Extension of the decomposition reveals residues as far apart as 25A from the binding epitope that can contribute significantly to binding free energy. These "hotspots" are found to show large atomic fluctuations in the unbound proteins, indicating that they reside in structurally less stable regions. Furthermore, hotspot residues experience a significantly larger-than-average decrease in local fluctuations upon complex formation. Finally, by calculating a pair-wise decomposition of interactions, interaction pathways originating in the binding epitope of Raf are found that protrude through the protein structure towards the loop L1. This explains the finding of a conformational change in this region upon complex formation with Ras, and it may trigger a larger structural change in Raf, which is considered to be necessary for activation of the effector by Ras.     

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.     

Calculations of electrostatic energies in proteins. The energetics of ionized groups in bovine pancreatic trypsin inhibitor

The importance of including different energy contributions in calculations of electrostatic energies in proteins is examined by calculating the intrinsic pKa values of the acidic groups of bovine pancreatic trypsin inhibitor. It appears that such calculations provide a powerful and revealing test; the relevant solvation energies of the ionized acids are of the order of -70 kcal/mol (1 cal = 4.184 J), and microscopic calculations that do not attempt to simulate the complete protein dielectric effect (including the surrounding solvent) can underestimate the solvation energy by as much as 50 kcal/mol. Reproducing correctly, by the same set of parameters, the solvation energies of ionized acids in different sites of a protein cannot be accomplished by including only part of the key energy contributions. The problems associated with macroscopic calculations are also considered and illustrated by the specific case of bovine pancreatic trypsin inhibitor. A promising approach is shown to be provided by a refinement of the previously developed Protein Dipoles Langevin Dipoles model. This model seems to represent consistently the microscopic dielectric of the protein and the surrounding water molecules. The model overcomes the problems associated with the macroscopic models (by treating explicitly the solvent molecules) and avoids the convergence problems associated with all-atom solvent models (by treating the average solvent polarization rather than averaging the actual polarization energy). This paper describes in detail the actual implementation of the model and examines its performance in evaluating intrinsic pKa values. Preliminary microscopic considerations of charge-charge interactions are presented.     

Quantitative analysis of DNA secondary structure from solvent-accessible surfaces: the B- to Z-DNA transition as a model

Solvent structure and its interactions have been suggested to play a critical role in defining the conformation of polynucleotides and other macromolecules. In this work, we attempt to quantitate solvent effects on the well-studied conformational transition between right-handed B- and left-handed Z-DNA. The solvent-accessible surfaces of the hexamer sequences d(m5CG)3, d(CG)3, d(CA)3, and d(TA)3 were calculated in their B- and Z-DNA conformations. The difference in hydration free energies between the Z and the B conformations (delta delta GH(Z-B] was determined from these surfaces to be -0.494 kcal/mol for C-5 methylated d(CG), 0.228 kcal/mol for unmethylated d(CG), 0.756 kcal/mol for d(CA)-d(TG), and 0.896 kcal/mol for d(TA) dinucleotides. These delta delta GH(Z-B) values were compared to the experimental B- to Z-DNA transition energies of -0.56 kcal/mol that we measured for C-5 methylated d(CG), 0.69-1.30 kcal/mol reported for unmethylated d(CG), 1.32-1.48 kcal/mol reported for d(CA)-d(TG), and 2.3-2.4 kcal/mol for d(TA) dinucleotides. From this comparison, we found that the calculated delta delta GH(Z-B) of these dinucleotides could account for the previous observation that the dinucleotides were ordered as d(m5CG) greater than d(CG) greater than d(CA)-d(TG) greater than d(TA) in stability as Z-DNA. Furthermore, we predicted that one of the primary reasons for the inability of d(TA) sequences to form Z-DNA results from a decrease in exposed hydrophilic surfaces of adjacent base pairs due to the C-5 methyl group of thymine; thus, d(UA) dinucleotides should be more stable as Z-DNA than the analogous d(TA) dinucleotides.(ABSTRACT TRUNCATED AT 250 WORDS)     

Thermodynamic measurements and predictions of the adsorption of short-chain peptides on nanothin polymer films

This contribution describes experimental measurements of submolecular-level interaction energies involved in the process of peptide adsorption on polymer films. The objective of this study was to use surface plasmon resonance (SPR) spectroscopy to measure the Gibbs energy change on adsorption (DeltaG(ad)) for pairs of various homopeptides on highly uniform, nanothin polymer films and to use these data, along with the principle of additivity, to predict DeltaG(ad) for homologous homopeptides, as well as for a mixed-residue peptide. By using a graft polymerization methodology, a nanothin poly(2-vinylpyridine) film was prepared and adsorption energies were measured first for a homologous series of tyrosine (Y) homopeptides on this film to determine submolecular-level interaction energies. By using SPR, adsorption isotherms were measured for YY and YYY peptides; analysis of these isotherms provided DeltaG(ad) data for a midchain tyrosine unit and a set of chain-end tyrosine units; values were -0.75 +/- 0.07 kcal/mol and -2.12 +/- 0.04 kcal/mol, respectively. Combining the thermodynamic contributions for adsorption of individual tyrosine units allowed a predictive estimate of -5.12 +/- 0.32 kcal/mol for the adsorption energy for YYYYYY; this estimate deviated by only 2.3% from its measured value of -5.24 +/- 0.06 kcal/mol. Similarly, adsorption energies were found for phenylalanine, glycine, and tyrosine-leucine peptides. Combining the thermodynamic contributions for adsorption of individual residue units allowed a predictive estimate of -3.24 +/- 0.38 kcal/mol for a pentapeptide, leucine enkephalin; this estimate deviated by only 3.0% from its measured value of -3.34 +/-0.11 kcal/mol.