Can the stability of protein mutants be predicted by free energy calculations?

Effect of alanine versus serine at position 88 of human transthyretin mutants on the protein stability

Human transthyretin (TTR) is a homo-tetrameric plasma protein associated with a high percentage of β-sheet forming amyloid fibrils. It accumulates in tissues or extracellular matrices to cause amyloid diseases. Free energy simulations with thermodynamic integration based on all-atom molecular dynamics simulations have been carried out to analyze the effects of the His88 → Ala and Ser mutations on the stability of human TTR. The calculated free energy change differences (ΔΔG) caused by the His88 → Ala and His88 → Ser mutations are -1.84 ± 0.86 and 7.56 ± 0.55 kcal/mol, respectively, which are in excellent agreement with prior reported experimental values. The simulation results show that the H88A mutant is more stable than the wild type, whereas the H88S mutant is less stable than the wild type. The free energy component analysis shows that the contribution to the free energy change difference (ΔΔG) for the His88 → Ala and His88 → Ser mutations mainly arise from electrostatic and van der Waals interactions, respectively. The electrostatic term stabilizes the H88A mutant more than the wild type, but the van der Waals interaction destabilizes the H88S mutant relative to the wild type. Individual residue contributions to the free energy change show neighboring residues exert stabilizing and destabilizing influence on the mutants. The implications of the simulation results for understanding the stabilizing and destabilizing effect and its contribution to protein stability are discussed.

Free energy simulations to study mutational effect of a conserved residue, Trp24, on stability of human serum retinol-binding protein

Human serum retinol-binding protein (RBP) is a plasma transport protein for vitamin A. RBP is a prime subclass of lipocalins, which bind nonpolar ligands within a β-barrel. To understand the role of Trp 24, one of the highly conserved residues in RBP, free energy simulations have been carried out to understand the effects of the mutations from Trp at position 24 to Leu, Phe, and Tyr in the apo-RBP on its thermal stability. We examine various unfolded systems to study the dependence of the free energy differences on the denatured structure. Our calculated free energy difference values for the three mutations are in excellent agreement with the experimental values when the initial coordinates of the seven-residue peptide segments truncated from the crystal structure are used for the denatured systems. Our free energy change differences for the Trp→Leu, Trp→Phe, and Trp→Tyr mutations are 2.50 ± 0.69, 2.58 ± 0.50, and 2.49 ± 0.48 kcal/mol, respectively, when the native-like seven-residue peptides are used as models for the denatured systems. The main contributions to the free energy change differences for the Trp24→Leu and Trp24→Phe mutations are mainly from van der Waals and covalent interactions, respectively. Electrostatic, van der Waals and covalent terms equally contribute to the free energy change difference for the Trp24→Tyr mutation. The free energy simulation helps understand the detailed microscopic mechanism of the stability of the RBP mutants relative to the wild type and the role of the highly conserved residue, Trp24, of the human RBP.Communicated by Ramaswamy H. Sarma.

Exploring the Effect of Enhanced Sampling on Protein Stability Prediction

Changes in protein stability due to side-chain mutations are evaluated by alchemical free-energy calculations based on classical molecular dynamics (MD) simulations in explicit solvent using the GROMOS force field. Three proteins of different complexity with a total number of 93 single-point mutations are analyzed, and the relative free-energy differences are discussed with respect to configurational sampling and (dis)agreement with experimental data. For the smallest protein studied, a 34-residue WW domain, the starting structure dependence of the alchemical free-energy changes, is discussed in detail. Deviations from previous simulations for the two other proteins are shown to result from insufficient sampling in the earlier studies. Hamiltonian replica exchange in combination with multiple starting structures and sufficient sampling time of more than 100 ns per intermediate alchemical state is required in some cases to reach convergence.

Free energy simulations to understand the effect of Met → Ala mutations at positions 205, 206 and 213 on stability of human prion protein

Prion diseases are a family of infectious amyloid diseases affecting human and animals. Prion propagation in transmissible spongiform encephalopathies is associated with the unfolding and conversion of normal cellular prion protein into its pathogenic scrapie form. Understanding the fundamentals of prion protein aggregation caused by mutations is crucial to unravel the pathology of prion diseases. To help understand the contributions of individual residues to the stability of the human prion protein, we have carried out free energy simulations based on atomistic molecular dynamics trajectories. We focus on Met → Ala mutations at positions 205, 206 and 213, which are mostly buried residues located on helix 3 of the protein. The simulations predicted that all three mutations destabilize the prion protein. Changes in unfolding free energy upon mutation, ∆∆G, are 3.10 ± 0.79, 2.00 ± 0.26 and 3.06 ± 0.66 kcal/mol for M205A, M206A and M213A, respectively, in excellent agreement with the corresponding experimental values of 3.09 ± 0.28, 1.50 ± 0.34 and 3.12 ± 0.27 kcal/mol [T. Hart et al. (2009) PNAS 106, 5651-5656]. Component analysis indicates that the major contributions to the loss of protein stability arise from van der Waals interactions for the M205A and M206A mutations, and from van der Waals and covalent energy terms for M213A. Interestingly, while free energy contributions from a majority of residues neighboring the mutation sites tend to stabilize the wild type, there are a few residues stabilizing the mutant side chains. Our results show that this approach to free energy calculation can be very useful for understanding the detailed mechanism of human prion protein stability.

Rapid, accurate, precise and reproducible ligand-protein binding free energy prediction

A central quantity of interest in molecular biology and medicine is the free energy of binding of a molecule to a target biomacromolecule. Until recently, the accurate prediction of binding affinity had been widely regarded as out of reach of theoretical methods owing to the lack of reproducibility of the available methods, not to mention their complexity, computational cost and time-consuming procedures. The lack of reproducibility stems primarily from the chaotic nature of classical molecular dynamics (MD) and the associated extreme sensitivity of trajectories to their initial conditions. Here, we review computational approaches for both relative and absolute binding free energy calculations, and illustrate their application to a diverse set of ligands bound to a range of proteins with immediate relevance in a number of medical domains. We focus on ensemble-based methods which are essential in order to compute statistically robust results, including two we have recently developed, namely thermodynamic integration with enhanced sampling and enhanced sampling of MD with an approximation of continuum solvent. Together, these form a set of rapid, accurate, precise and reproducible free energy methods. They can be used in real-world problems such as hit-to-lead and lead optimization stages in drug discovery, and in personalized medicine. These applications show that individual binding affinities equipped with uncertainty quantification may be computed in a few hours on a massive scale given access to suitable high-end computing resources and workflow automation. A high level of accuracy can be achieved using these approaches.

Free energy calculations on the stability of the 14-3-3ζ protein

Mutations of cysteine are often introduced to e.g. avoid formation of non-physiological inter-molecular disulfide bridges in in-vitro experiments, or to maintain specificity in labeling experiments. Alanine or serine is typically preferred, which usually do not alter the overall protein stability, when the original cysteine was surface exposed. However, selecting the optimal mutation for cysteines in the hydrophobic core of the protein is more challenging. In this work, the stability of selected Cys mutants of 14-3-3ζ was predicted by free-energy calculations and the obtained data were compared with experimentally determined stabilities. Both the computational predictions as well as the experimental validation point at a significant destabilization of mutants C94A and C94S. This destabilization could be attributed to the formation of hydrophobic cavities and a polar solvation of a hydrophilic side chain. A L12E, M78K double mutant was further studied in terms of its reduced dimerization propensity. In contrast to naïve expectations, this double mutant did not lead to the formation of strong salt bridges, which was rationalized in terms of a preferred solvation of the ionic species. Again, experiments agreed with the calculations by confirming the monomerization of the double mutants. Overall, the simulation data is in good agreement with experiments and offers additional insight into the stability and dimerization of this important family of regulatory proteins.

Combining structural modeling with ensemble machine learning to accurately predict protein fold stability and binding affinity effects upon mutation

Advances in sequencing have led to a rapid accumulation of mutations, some of which are associated with diseases. However, to draw mechanistic conclusions, a biochemical understanding of these mutations is necessary. For coding mutations, accurate prediction of significant changes in either the stability of proteins or their affinity to their binding partners is required. Traditional methods have used semi-empirical force fields, while newer methods employ machine learning of sequence and structural features. Here, we show how combining both of these approaches leads to a marked boost in accuracy. We introduce ELASPIC, a novel ensemble machine learning approach that is able to predict stability effects upon mutation in both, domain cores and domain-domain interfaces. We combine semi-empirical energy terms, sequence conservation, and a wide variety of molecular details with a Stochastic Gradient Boosting of Decision Trees (SGB-DT) algorithm. The accuracy of our predictions surpasses existing methods by a considerable margin, achieving correlation coefficients of 0.77 for stability, and 0.75 for affinity predictions. Notably, we integrated homology modeling to enable proteome-wide prediction and show that accurate prediction on modeled structures is possible. Lastly, ELASPIC showed significant differences between various types of disease-associated mutations, as well as between disease and common neutral mutations. Unlike pure sequence-based prediction methods that try to predict phenotypic effects of mutations, our predictions unravel the molecular details governing the protein instability, and help us better understand the molecular causes of diseases.

The key to predicting the stability of protein mutants lies in an accurate description and proper configurational sampling of the folded and denatured states

BACKGROUND

The contribution of particular hydrogen bonds to the stability of a protein fold can be investigated experimentally as well as computationally by the construction of protein mutants which lack particular hydrogen-bond donors or acceptors with a subsequent determination of their structural stability. However, the comparison of experimental data with computational results is not straightforward. One of the difficulties is related to the representation of the unfolded state conformation.

METHODS

A series of molecular dynamics simulations of the 34-residue WW domain of protein Pin1 and 20 amide-to-ester mutants started from the X-ray crystal structure and the NMR solution structure are analysed in terms of backbone-backbone hydrogen bonding and differences in free enthalpies of folding in order to provide a structural interpretation of the experimental data available.

RESULTS

The contribution of the different β-sheet hydrogen bonds to the relative stability of the mutants with respect to wild type cannot be directly inferred from experimental thermal denaturation temperatures or free enthalpies of chaotrope denaturation for the different mutants, because some β-sheet hydrogen bonds show sizeable variation in occurrence between the different mutants.

CONCLUSIONS

A proper representation of unfolded state conformations appears to be essential for an adequate description of relative stabilities of protein mutants.

GENERAL SIGNIFICANCE

The simulations may be used to link the structural Boltzmann ensembles to relative free enthalpies of folding between mutants and wild-type protein and show that unfolded conformations have to be treated with a sufficient level of detail in free energy calculations of protein stability. This article is part of a Special Issue entitled Recent developments of molecular dynamics.

SDM--a server for predicting effects of mutations on protein stability and malfunction

The sheer volume of non-synonymous single nucleotide polymorphisms that have been generated in recent years from projects such as the Human Genome Project, the HapMap Project and Genome-Wide Association Studies means that it is not possible to characterize all mutations experimentally on the gene products, i.e. elucidate the effects of mutations on protein structure and function. However, automatic methods that can predict the effects of mutations will allow a reduced set of mutations to be studied. Site Directed Mutator (SDM) is a statistical potential energy function that uses environment-specific amino-acid substitution frequencies within homologous protein families to calculate a stability score, which is analogous to the free energy difference between the wild-type and mutant protein. Here, we present a web server for SDM (http://www-cryst.bioc.cam.ac.uk/~sdm/sdm.php), which has obtained more than 10 000 submissions since being online in April 2008. To run SDM, users must upload a wild-type structure and the position and amino acid type of the mutation. The results returned include information about the local structural environment of the wild-type and mutant residues, a stability score prediction and prediction of disease association. Additionally, the wild-type and mutant structures are displayed in a Jmol applet with the relevant residues highlighted.

Alchemical free energy simulations for biological complexes: powerful but temperamental

Free energy simulations compare multiple ligand:receptor complexes by "alchemically" transforming one into another, yielding binding free energy differences. Since their introduction in the 1980s, many technical and theoretical obstacles were surmounted, and the method ("MDFE," since molecular dynamics are often used) has matured into a powerful tool. We describe its current status, its effectiveness, and the challenges it faces. MDFE has provided chemical accuracy for many systems but remains expensive, with significant human overhead costs. The bottlenecks have shifted, partly due to increased computer power. To study diverse sets of ligands, force field availability and accuracy can be a major difficulty. Another difficulty is the frequent need to consider multiple states, related to sidechain protonation or buried waters, for example. Sophisticated, automated methods to sample these states are maturing, such as constant pH simulations. Meanwhile, combinations of MDFE and simpler approaches, like continuum dielectric models, can be very effective. As illustrations, we show how, with careful force field parameterization, MDFE accurately predicts binding specificities between complex tetracycline ligands and their targets. We describe substrate binding to the aspartyl-tRNA synthetase enzyme, where many distinct electrostatic states play a role, and a histidine and a Mg(2+) ion act as coupled switches that help enforce a strict preference for the aspartate substrate, relative to several analogs. Overall, MDFE has achieved a predictive status, where novel ligands can be studied and molecular recognition elucidated in depth. It should play an increasing role in the analysis of complex cellular processes and biomolecular engineering.

Calculation of the Redox Potential of the Protein Azurin and Some Mutants

Azurin from Pseudomonas aeruginosa is a small 128-residue, copper-containing protein. Its redox potential can be modified by mutating the protein. Free-energy calculations based on classical molecular-dynamics simulations of the protein and from mutants in aqueous solution at different pH values were used to compute relative redox potentials. The precision of the free-energy calculations with the lambda coupling-parameter approach is evaluated as function of the number and sequence of lambda values, the sampling time and initial conditions. It is found that the precision is critically dependent on the relaxation of hydrogen-bonding networks when changing the atomic-charge distribution due to a change of redox state or pH value. The errors in the free energies range from 1 to 10 k(B)T, depending on the type of process. Only qualitative estimates of the change in redox potential by protein mutation can be obtained.

Protein sequence randomization: efficient estimation of protein stability using knowledge-based potentials

Modifications of the amino acid sequence generally affect protein stability. Here, we use knowledge-based potentials to estimate the stability of protein structures under sequence variation. Calculations on a variety of protein scaffolds result in a clear distinction of known mutable regions from arbitrarily chosen control patches. For example, randomly changing the sequence of an antibody paratope yields a significantly lower number of destabilized mutants as compared to the randomization of comparable regions on the protein surface. The technique is computationally efficient and can be used to screen protein structures for regions that are amenable to molecular tinkering by preserving the stability of the mutated proteins.

Simulating proteins at constant pH: An approach combining molecular dynamics and Monte Carlo simulation

For the structure and function of proteins, the pH of the solution is one of the determining parameters. Current molecular dynamics (MD) simulations account for the solution pH only in a limited way by keeping each titratable site in a chosen protonation state. We present an algorithm that generates trajectories at a Boltzmann distributed ensemble of protonation states by a combination of MD and Monte Carlo (MC) simulation. The algorithm is useful for pH-dependent structural studies and to investigate in detail the titration behavior of proteins. The method is tested on the acidic residues of the protein hen egg white lysozyme. It is shown that small structural changes may have a big effect on the pK(A) values of titratable residues.

Contribution of Thr29 to the thermodynamic stability of goat alpha-lactalbumin as determined by experimental and theoretical approaches

The Thr29 residue in the hydrophobic core of goat alpha-lactalbumin (alpha-LA) was substituted with Val (Thr29Val) and Ile (Thr29Ile) to investigate the contribution of Thr29 to the thermodynamic stability of the protein. We carried out protein stability measurements, X-ray crystallographic analyses, and free energy calculations based on molecular dynamics simulation. The equilibrium unfolding transitions induced by guanidine hydrochloride demonstrated that the Thr29Val and Thr29Ile mutants were, respectively, 1.9 and 3.2 kcal/mol more stable than the wild-type protein (WT). The overall structures of the mutants were almost identical to that of WT, in spite of the disruption of the hydrogen bonding between the side-chain O-H group of Thr29 and the main-chain C=O group of Glu25. To analyze the stabilization mechanism of the mutants, we performed free energy calculations. The calculated free energy differences were in good agreement with the experimental values. The stabilization of the mutants was mainly caused by solvation loss in the denatured state. Furthermore, the O-H group of Thr29 favorably interacts with the C=O group of Glu25 to form hydrogen bonds and, simultaneously, unfavorably interacts electrostatically with the main-chain C=O group of Thr29. The difference in the free energy profile of the unfolding path between WT and the Thr29Ile mutant is discussed in light of our experimental and theoretical results.

Computed free energy differences between point mutations in a collagen-like peptide

We studied the results of mutating alanine --> glycine at three positions of a collagen-like peptide in an effort to develop a computational method for predicting the energetic and structural effects of a single point genetic mutation in collagen, which is associated with the clinical diagnosis of Osteogenesis Imperfecta (OI). The differences in free energy of denaturation were calculated between the collagen-like peptides [(POG)(4)(POA)(POG)(4)](3) and [(POG)(10)](3) (POG: proline-hydroxyproline-glycine).* Our computational results, which suggest significant destabilization of the collagen-like triple-helix upon the glycine --> alanine mutations, correlate very well with the experimental free energies of denaturation. The robustness of our collagen-like peptide model is shown by its reproduction of experimental results with both different simulation paths and different lengths of the model peptide. The individual free energy for each alanine --> glycine mutation (and the reverse free energy, glycine --> alanine mutation) in the collagen-like peptide has been calculated. We find that the first alanine introduced into the triple helix causes a very large destabilization of the helix, but the last alanine introduced into the same position of an adjacent chain causes a very small change in the peptide stability. Thus, our results demonstrate that each mutation does not contribute equally to the free energy. We find that the sum of the calculated individual residues' free energy can accurately model the experimental free energy for the whole peptide.

Computational studies on mutant protein stability: The correlation between surface thermal expansion and protein stability

Thermal stability of mutant proteins has been investigated using temperature dependent molecular dynamics (MD) simulations in vacuo. The numerical modeling was aimed at mimicking protein expansion upon heating. After the conditions for an expanding protein accessible surface area were established for T4 lysozyme and barnase wild-type proteins, MD simulations were carried out under the same conditions using the crystal structures of several mutant proteins. The computed thermal expansion of the accessible surface area of mutant proteins was found to be strongly correlated with their experimentally measured stabilities. A similar, albeit weaker, correlation was observed for model mutant proteins. This opens the possibility of obtaining stability information directly from protein structure.

Assignment of side-chain conformation using adiabatic energy mapping, free energy perturbation, and molecular dynamic simulations

NMR spectroscopic analysis of the C-terminal Kunitz domain fragment (alpha3(VI)) from the human alpha3-chain of type VI collagen has revealed that the side chain of Trp21 exists in two unequally populated conformations. The major conformation (M) is identical to the conformation observed in the X-ray crystallographic structure, while the minor conformation (m) cannot structurally be resolved in detail by NMR due to insufficient NOE data. In the present study, we have applied: (1) rigid and adiabatic mapping, (2) free energy simulations, and (3) molecular dynamic simulations to elucidate the structure of the m conformer and to provide a possible pathway of the Trp21 side chain between the two conformers. Adiabatic energy mapping of conformations of the Trp21 side chain obtained by energy minimization identified two energy minima: One corresponding to the conformation of Trp21 observed in the X-ray crystallographic structure and solution structure of alpha3(VI) (the M conformation) and the second corresponding to the m conformation predicted by NMR spectroscopy. A transition pathway between the M and m conformation is suggested. The free-energy difference between the two conformers obtained by the thermodynamic integration method is calculated to 1.77+/-0.7 kcal/mol in favor of the M form, which is in good agreement with NMR results. Structural and dynamic properties of the major and minor conformers of the alpha3(VI) molecule were investigated by molecular dynamic. Essential dynamics analysis of the two resulting 800 ps trajectories reveals that when going from the M to the m conformation only small, localized changes in the protein structure are induced. However, notable differences are observed in the mobility of the binding loop (residues Thr13-Ile18), which is more flexible in the m conformation than in the M conformation. This suggests that the reorientation of Trp2 might influence the inhibitory activity against trypsin, despite the relative large distance between the binding loop and Trp21.

Ligand solvation in molecular docking

Solvation plays an important role in ligand-protein association and has a strong impact on comparisons of binding energies for dissimilar molecules. When databases of such molecules are screened for complementarity to receptors of known structure, as often occurs in structure-based inhibitor discovery, failure to consider ligand solvation often leads to putative ligands that are too highly charged or too large. To correct for the different charge states and sizes of the ligands, we calculated electrostatic and non-polar solvation free energies for molecules in a widely used molecular database, the Available Chemicals Directory (ACD). A modified Born equation treatment was used to calculate the electrostatic component of ligand solvation. The non-polar component of ligand solvation was calculated based on the surface area of the ligand and parameters derived from the hydration energies of apolar ligands. These solvation energies were subtracted from the ligand-receptor interaction energies. We tested the usefulness of these corrections by screening the ACD for molecules that complemented three proteins of known structure, using a molecular docking program. Correcting for ligand solvation improved the rankings of known ligands and discriminated against molecules with inappropriate charge states and sizes.

Dependence of protein stability on the structure of the denatured state: free energy calculations of I56V mutation in human lysozyme

Free energy calculations were carried out to understand the effect of the I56V mutation of human lysozyme on its thermal stability. In the simulation of the denatured state, a short peptide including the mutation site in the middle is employed. To study the dependence of the stability on the denatured-state structure, five different initial conformations, native-like, extended, and three random-coil-like conformations, were examined. We found that the calculated free energy difference, DeltaDeltaGcal, depends significantly on the structure of the denatured state. When native-like structure is employed, DeltaDeltaGcal is in good agreement with the experimental free energy difference, DeltaDeltaGexp, whereas in the other four models, DeltaDeltaGcal differs sharply from DeltaDeltaGexp. It is therefore strongly suggested that the structure around the mutation site takes a native-like conformation rather than an extended or random-coil conformation. From the free energy component analysis, it has been shown that free energy components originating from Lennard-Jones and covalent interactions dominantly determine DeltaDeltaGcal. The contribution of protein-protein interactions to the nonbonded component of DeltaDeltaGcal is about the same as that from protein-water interactions. The residues that are located in a hydrophobic core (F3, L8, Y38, N39, T40, and I89) contribute significantly to the nonbonded free energy component of DeltaDeltaGcal. We also propose a general computational strategy for the study of protein stability that is equally conscious of the denatured and native states.

Computational analysis of thermal stability: effect of Ile-->Val mutations in human lysozyme

BACKGROUND

Free energy calculations are carried out to study the change of thermal stability caused by Ile23-->Val, Ile56-->Val, Ile89-->Val and Ile106-->Val mutations in human lysozyme. In order to examine the dependence of the free energy difference, DeltaDeltaG, on the denatured-state structure, extended and native-like conformations are employed as initial conformations in the denatured-state simulations.

RESULTS

Calculated values of DeltaDeltaG for the mutations, Ile56-->Val, Ile89-->Val and Ile106-->Val, were in good agreement with experimental values when the native-like structure was employed in the respective denatured-state simulations. In the case of Ile23-->Val, a considerable difference between the calculated and experimental values of DeltaDeltaG was observed.

CONCLUSIONS

The physical nature of Ile56-->Val, Ile89-->Val and Ile106-->Val mutations was rationally characterized by a free energy component analysis. It is suggested that the alpha domain in which Ile23 is included is considerably structured even in the denatured state.

Conformation of the Ras-binding domain of Raf studied by molecular dynamics and free energy simulations

Recognition of Ras by its downstream target Raf is mediated by a Ras-recognition region in the Ras-binding domain (RBD) of Raf. Residues 78-89 in this region occupy two different conformations in the ensemble of NMR solution structures of the RBD: a fully alpha-helical one, and one where 87-90 form a type IV beta-turn. Molecular dynamics simulations of the RBD in solution were performed to explore the stability of these and other possible conformations of both the wild-type RBD and the R89K mutant, which does not bind Ras. The simulations sample a fully helical conformation for residues 78-89 similar to the NMR helical structures, a conformation where 85-89 form a 3(10)-helical turn, and a conformation where 87-90 form a type I beta-turn, whose free energies are all within 0.3 kcal/mol of each other. NOE patterns and H(alpha) chemical shifts from the simulations are in reasonable agreement with experiment. The NMR turn structure is calculated to be 3 kcal/mol higher than the three above conformations. In a simulation with the same implicit solvent model used in the NMR structure generation, the turn conformation relaxes into the fully helical conformation, illustrating possible structural artifacts introduced by the implicit solvent model. With the Raf R89K mutant, simulations sample a fully helical and a turn conformation, the turn being 0.9 kcal/mol more stable. Thus, the mutation affects the population of RBD conformations, and this is expected to affect Ras binding. For example, if the fully helical conformation of residues 78-89 is required for binding, its free energy increase in R89K will increase the binding free energy by about 0.6 kcal/mol.

The nature of trypsin-pancreatic trypsin inhibitor binding: free energy calculation of Tyr39-->Phe39 mutation in trypsin

The main goal of this work is the detailed study of the binding interactions in the trypsin-pancreatic trypsin inhibitor (PTI) complex and, here, we present how meaningful the Tyr39-Ile19 interaction is to the stability of that particular complex using free energy methods. This knowledge should be very important in the design of new inhibitors for trypsin and enzymes homologous to it. In particular, it could help to decide whether it is possible to produce selective inhibitors for these enzymes by appropriate mutations of residues in the contact region of PTI.

Empirical free energy calculations: a blind test and further improvements to the method

Empirical Gibbs functions estimate free energies of non-covalent reactions (deltaG) from atomic coordinates of reaction products (e.g. antibody-antigen complexes). The function previously developed by us has four terms that quantify the effects of hydrophobic, electrostatic and entropy changes (conformational, association) upon complexation. The function was used to calculate delta deltaG of ten lysozyme mutants affecting the stability of the HyHEL-10 antibody-lysozyme complex. The mutants were computer-modeled from the X-ray structure of the wild-type, and free energy calculations produced a correlation coefficient of 0.5 with the experimental delta deltaG data (average absolute error +/-3 kcal). The following changes were then introduced into the Gibbs function: (1) the hydrophobic force was made proportional to the molecular surface, as calculated by the GEPOL93 algorithm, with the scaling constant of 70 cal/mol/A2; (2) calculation of the electrostatics of binding was carried out by the finite difference Poisson-Boltzmann algorithm, which employed uniform grid charging, dielectric boundary smoothing and charge anti-aliasing; and (3) side-chain conformational entropy was estimated from the CONGEN sampling of torsional degrees of freedom. In the new calculations, correlation with experimental data improved to 0.6 or 0.8 if a single outlying mutant, K96M, was neglected. Analysis of the errors remaining in our calculations indicated that molecular mechanics-based modeling of the mutants, rather than the form of our amended Gibbs function, was the main factor limiting the accuracy of the free energy estimates.

Computational design of a substrate specificity mutant of a protein

The wild-type trp repressor of E. coli bound 5-methoxytryptophan, a Trp analogue, less tightly than Trp. A mutant repressor (Val58-->Ala) that should bind 5-methoxytryptophan preferentially to Trp was computationally designed by free-energy calculations accompanied by free-energy decomposition. The designed mutant was demonstrated by experiments to bind 5-methoxytryptophan more tightly than Trp, consistent with the computational prediction. This success indicates the usefulness of free energy decomposition in protein design.

Exploring the interaction between D-xylose isomerase and D-xylose by free energy calculation

The numerical quadrature thermodynamic integration method is used to investigate enzyme-substrate interaction of D-xylose isomerase. A screening function for the coulombic interaction is introduced into the simulation to correct the effect of finite cut-off radius for the non-bonded interaction. The binding free energy difference for D-xylose with D-xylose isomerase and its N184D mutant has been calculated, and the result 3.9 +/- 1.2 kJ/mol agrees well with experimental data of 4.38 kJ/mol. In addition, the structure and dynamics of enzyme-substrate complex were simulated for mutant and wild-type enzyme, respectively. Analysis of the structures and intramolecular interactions of the complexes were found to be valuable for understanding the reaction mechanism of the enzyme D-xylose isomerase.

Energetic decomposition of the alpha-helix-coil equilibrium of a dynamic model system

Using molecular dynamics simulations to calculate free energies of molecular transformation, we have computed helix-coil transition free energies for alanine oligomers up to 14 residues long. The simulations have been done on the model in vacuo with dielectric constant, epsilon = 1, 5, 25, and infinity and on the model in solution with explicit representation of water molecules and with partial charges on the oligomer set to zero. (The analogous simulations of the solvated model with full charges on the oligomer were reported elsewhere [L. Wang et al. (1995) Proceedings of the National Academy of Science USA 92, 10924-10928]). In vacuo, both entropic and electrostatic contributions oppose formation of a 3-residue helical nucleus in the helix initiation step. The entropy change opposing helix growth is found to be 3 e.u., van der Waals interactions favor helix growth by 1.9 kcal/mol, and electrostatic interactions favor helix growth by 3 kcal/mol (for epsilon = 1; all these values are per residue). In water, helix stability is slightly greater for the zero-charge model than for the full-charge model, i.e., the polypeptide's electrostatic interactions, which include hydrogen bonds, slightly destabilize the helix. The helix stabilizing contribution of the hydrophobic effect was found to be identical to that of the van der Waals interactions in vacuo (i.e., 1.9 kcal/mol per residue). The zero-charge model has nearly identical helix stability in vacuo and in water, the almost identical free energies of transfer of helix and coil state of the zero-charge oligomer from vacuum to water are found to be small. Thus, the results of this systematic variation of the force field afford a meaningful decomposition of the free energies for helix initiation and growth.

Thermodynamic stability of water molecules in the bacteriorhodopsin proton channel: a molecular dynamics free energy perturbation study

The proton transfer activity of the light-driven proton pump, bacteriorhodopsin (bR) in the photochemical cycle might imply internal water molecules. The free energy of inserting water molecules in specific sites along the bR transmembrane channel has been calculated using molecular dynamics simulations based on a microscopic model. The existence of internal hydration is related to the free energy change on transfer of a water molecule from bulk solvent into a specific binding site. Thermodynamic integration and perturbation methods were used to calculate free energies of hydration for each hydrated model from molecular dynamics simulations of the creation of water molecules into specific protein-binding sites. A rigorous statistical mechanical formulation allowing the calculation of the free energy of transfer of water molecules from the bulk to a protein cavity is used to estimate the probabilities of occupancy in the putative bR proton channel. The channel contains a region lined primarily by nonpolar side-chains. Nevertheless, the results indicate that the transfer of four water molecules from bulk water to this apparently hydrophobic region is thermodynamically permitted. The column forms a continuous hydrogen-bonded chain over 12 A between a proton donor, Asp 96, and the retinal Schiff base acceptor. The presence of two water molecules in direct hydrogen-bonding association with the Schiff base is found to be strongly favorable thermodynamically. The implications of these results for the mechanism of proton transfer in bR are discussed.

Computational screening of combinatorial libraries

We first review existing computational methods with an intrinsic combinatorial feature, then describe a new computational method for screening combinatorial libraries using a recently developed multicopy sampling technique. The new method differs from the existing ones in that it can be used to screen simultaneously an entire library of molecules, instead of the individual molecules in a library. As an example, we have applied the method to study site-directed amino acid substitutions in a protein. After two rounds of library screening, we identified the energetically most stable substitutions along with their optimal conformations from all natural amino acids. In principle, the method is generally applicable to study ligand-host systems.

Understanding the thermodynamic stability of an RNA hairpin and its mutant

The thermodynamic stability of RNA hairpin loops has been a subject of considerable interest in the recent past (Wimberly et al., 1991). There have been experimental reports indicating that the hairpins with a C(UUCG)G loop sequence are thermodynamically very stable (Wimberly et al., 1991). We used the solution structure of GGAC(UUCG)GUCC (Cheong et al., 1990; Varani et al., 1991) as the starting conformation in our attempt to understand its thermodynamic stability. We carried out molecular dynamics/free energy simulations to understand the basis for the destabilization of the C(UUCG)G loop by mutating cytosine (C7)-->uracil. Because of the limited length of simulation and the presence of kinetic barriers (solvent intervention) to the uracil-->cytosine mutation, all of our computed free energy differences are based on multiple forward simulations. Based on these calculations we find that the cytosine-->uracil mutation in the loop destabilizes it by approximately 1.5kcal/mol relative to that of the reference state, an A-form RNA but with cytosine (C7) looped out. This is the same sign and magnitude as that observed in the thermodynamic studies carried out by Varani et al.(1991). We have carried out free energy component analysis to understand the effect of mutating the cytosine residue to uracil on the thermodynamic stability of the C(UUCG)G hairpin loops. Our calculations show that the most significant contribution to the stability is from the phosphate group linking U5 and U6, which favors the cytosine residue over uracil by about 6.0 kcal/mol. The residues U5, U6, and G8 in the loop region also contribute significantly to the stability. The contributions from the salt and solvent compensate each other, indicating the dynamic nature of interactions of the environment with the nucleic acid system and the coupling between these two components.

New approaches in molecular structure prediction

In the past years, much effort has been put on the development of new methodologies and algorithms for the prediction of protein secondary and tertiary structures from (sequence) data; this is reviewed in detail. New approaches for these predictions such as neural network methods, genetic algorithms, machine learning, and graph theoretical methods are discussed. Secondary structure prediction algorithms were improved mostly by considering families of related proteins; however, for the reliable tertiary structure modeling of proteins, knowledge-based techniques are still preferred. Methods and examples with more or less successful results are described. Also, programs and parameterizations for energy minimisations, molecular dynamics, and electrostatic interactions have been improved, especially with respect to their former limits of applicability. Other topics discussed in this review include the use of traditional and on-line databases, the docking problem and surface properties of biomolecules, packing of protein cores, de novo design and protein engineering, prediction of membrane protein structures, the verification and reliability of model structures, and progress made with currently available software and computer hardware. In summary, the prediction of the structure, function, and other properties of a protein is still possible only within limits, but these limits continue to be moved.

Predicting the structure of protein complexes: a step in the right direction

In a blind test of protein-docking algorithms, six groups used different methods to predict the structure of a protein complex. All six predicted structures were close enough to the experimental complex to be useful; nevertheless, several important details of the experimental complex were missed or only partially predicted.

Accuracy and reliability of the scaling-relaxation method for loop closure: an evaluation based on extensive and multiple copy conformational samplings

The accuracy and reliability of the recently proposed scaling-relaxation method for loop closure were examined by using extensive conformational sampling. For each of the eight heptapeptides chosen to represent a variety of protein conformations, 1,000-2,000 conformations were sampled. Each segment contained 14 rotatable backbone dihedral angles. The average root mean square deviations (RMSDs) between the predicted and the native conformations were 0.7 angstrom for the backbone and 1.2 angstroms for the side chain atoms. These predictions were substantially more accurate than the previous predictions (1.1 angstroms for the backbone and 2.2 angstroms for the side chain atoms) of the same eight protein segments based on limited conformational sampling (100 conformations for each segment). Large prediction errors mostly occurred at polar and surface side chains that are unlikely to have any meaningful conformation. Moreover, the reliability of seven of the eight predictions was demonstrated with their energy-RMSD and stability-RMSD correlations of the low-energy conformations, where the conformational stability was estimated by using the multiple copy simultaneous sampling method.

Testing homology modeling on mutant proteins: predicting structural and thermodynamic effects in the Ala98-->Val mutants of T4 lysozyme

BACKGROUND

Current approaches to homology modeling predict how amino acid substitutions will alter a protein's structure, primarily by modeling sidechain conformations upon essentially immobile backbone frameworks. However, recent crystal structures of T4 lysozyme mutants reveal significant shifts of the mainchain and other potentially serious problems for sidechain rotamer-based modeling. This paper evaluates the accuracy of structural and thermodynamic predictions from two common sidechain modeling approaches to measure errors caused by the fixed-backbone approximation.

RESULTS

Tested on a series of T4 lysozyme mutants, this sidechain rotamer library approach did not handle mainchain shifts well, correctly predicting the sidechain conformations of only two of six mutants. By contrast, allowing sidechains to move more flexibly appeared to compensate for the rigidity of the mainchain and gave reasonably accurate coordinate predictions (rms errors of 0.5-1.0 A for each mutated sidechain), better on average than 90% of possible conformations. The calculated packing energies correlated well with experimental stabilities (r2 = 0.81) and correctly captured the cooperative interactions of several neighboring mutations.

CONCLUSIONS

Mutant modeling can be relatively accurate despite the fixed-backbone approximation. Mainchain shifts (0.2-0.5 A) cause increased sidechain coordinate errors of 0.1-0.8 A, torsional errors of 10-30 degrees, and exaggerated strain energy for overpacked mutants, compared with the same calculations performed with the correct mutant backbones.