Simulation analysis of the stability mutant R96H of T4 lysozyme

Alchemical Free Energy Calculations on Membrane-Associated Proteins

Membrane proteins have diverse functions within cells and are well-established drug targets. The advances in membrane protein structural biology have revealed drug and lipid binding sites on membrane proteins, while computational methods such as molecular simulations can resolve the thermodynamic basis of these interactions. Particularly, alchemical free energy calculations have shown promise in the calculation of reliable and reproducible binding free energies of protein-ligand and protein-lipid complexes in membrane-associated systems. In this review, we present an overview of representative alchemical free energy studies on G-protein-coupled receptors, ion channels, transporters as well as protein-lipid interactions, with emphasis on best practices and critical aspects of running these simulations. Additionally, we analyze challenges and successes when running alchemical free energy calculations on membrane-associated proteins. Finally, we highlight the value of alchemical free energy calculations calculations in drug discovery and their applicability in the pharmaceutical industry.

Metal-Promoted Higher-Order Assembly of Disulfide-Stapled Helical Barrels

Peptide-based helical barrels are a noteworthy building block for hierarchical assembly, with a hydrophobic cavity that can serve as a host for cargo. In this study, disulfide-stapled helical barrels were synthesized containing ligands for metal ions on the hydrophilic face of each amphiphilic peptide helix. The major product of the disulfide-stapling reaction was found to be composed of five amphiphilic peptides, thereby going from a 16-amino-acid peptide to a stapled 80-residue protein in one step. The structure of this pentamer, 5HB1, was optimized in silico, indicating a significant hydrophobic cavity of ~6 Å within a helical barrel. Metal-ion-promoted assembly of the helical barrel building blocks generated higher order assemblies with a three-dimensional (3D) matrix morphology. The matrix was decorated with hydrophobic dyes and His-tagged proteins both before and after assembly, taking advantage of the hydrophobic pocket within the helical barrels and coordination sites within the metal ion-peptide framework. As such, this peptide-based biomaterial has potential for a number of biotechnology applications, including supplying small molecule and protein growth factors during cell and tissue growth within the matrix.

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.

Repulsive Soft-Core Potentials for Efficient Alchemical Free Energy Calculations

In alchemical free energy (FE) simulations, annihilation and creation of atoms are generally achieved with the soft-core potential that shifts the interparticle separations. While this soft-core potential eliminates the numerical instability occurring near the two end states of the transformation, it makes the hybrid Hamiltonian vary nonlinearly with respect to the parameter λ, which interpolates between the Hamiltonians representing the two end states. This complicates FE estimation by Bennett acceptance ratio (BAR), free energy perturbation (FEP), and thermodynamic integration (TI) methods, thus reducing their calculation efficiency. In this work, we develop a new type of repulsive soft-core potential, called Gaussian soft-core (GSC) potential, with two parameters controlling its maximum and width. The main advantage of this potential is the linearity of the hybrid Hamiltonian with respect to λ, thus permitting the direct application of BAR, FEP, TI, and other variant FE methods. The accuracy and efficiency of the GSC potential are demonstrated by comparing the free energies of annihilation determined for 13 small molecules and an alchemical mutation of a protein side chain. In addition, in combination with a TI integrand (∂H/∂λ) estimation strategy, we show that GSC can considerably reduce the number of λ simulations compared to the commonly used separation-shifted soft-core potential.

Improving the Accuracy of Protein Thermostability Predictions for Single Point Mutations

Accurately predicting the protein thermostability changes upon single point mutations in silico is a challenge that has implications for understanding diseases as well as industrial applications of protein engineering. Free energy perturbation (FEP) has been applied to predict the effect of single point mutations on protein stability for over 40 years and emerged as a potentially reliable prediction method with reasonable throughput. However, applications of FEP in protein stability calculations in industrial settings have been hindered by a number of limitations, including the inability to model mutations to and from prolines in which the bonded topology of the backbone is modified and the complexity in modeling charge-changing mutations. In this study, we have extended the FEP+ protocol to enable the accurate modeling of the effects on protein stability from proline mutations and from charge-changing mutations. We also evaluated the influence of the unfolded model in the stability calculations using increasingly longer peptides with native sequence and conformations. With the abovementioned improvements, the accuracy of FEP predictions of protein stability over a data set of 87 mutations on five different proteins has drastically improved compared with previous studies, with a mean unsigned error of 0.86 kcal/mol and root mean square error of 1.11 kcal/mol, comparable with the accuracy of previously published state-of-the-art small-molecule relative binding affinity calculations, which have been shown to be capable of driving discovery projects.

Learning the changes of barnase mutants thermostability from structural fluctuations obtained using anisotropic network modeling

In biotechnology applications, rational design of new proteins with improved physico-chemical properties includes a number of important tasks. One of the greatest practical and fundamental challenges is the design of highly thermostable protein enzymes that maintain catalytic activity at high temperatures. This problem may be solved by introducing mutations into the wild-type enzyme protein. In this work, to predict the impact of such mutations in barnase protein we applied the anisotropic network modeling approach, revealing atomic fluctuations in structural regions that are changed in mutants compared to the wild-type protein. A regression model was constructed based on these structural features that can allow one to predict the thermal stability of new barnase mutants. Moreover, the analysis of regression model provides a mechanistic explanation of how the structural features can contribute to the thermal stability of barnase mutants.

Computational chemistry methods to investigate the effects caused by DNA variants linked with disease

Computational chemistry offers variety of tools to study properties of biological macromolecules. These tools vary in terms of levels of details from quantum mechanical treatment to numerous macroscopic approaches. Here, we provide a review of computational chemistry algorithms and tools for modeling the effects of genetic variations and their association with diseases. Particular emphasis is given on modeling the effects of missense mutations on stability, conformational dynamics, binding, hydrogen bond network, salt bridges, and pH-dependent properties of the corresponding macromolecules. It is outlined that the disease may be caused by alteration of one or several of above-mentioned biophysical characteristics, and a successful prediction of pathogenicity requires detailed analysis of how the alterations affect the function of involved macromolecules. The review provides a short list of most commonly used algorithms to predict the molecular effects of mutations as well.

Sequence Tolerance of a Single-Domain Antibody with a High Thermal Stability: Comparison of Computational and Experimental Fitness Profiles

The sequence fitness of a llama single-domain antibody with an unusually high thermal stability is explored by a combined computational and experimental study. Starting with the X-ray crystallographic structure, RosettaBackrub simulations were applied to model sequence-structure tolerance profiles and identify key substitution sites. From the model calculations, an experimental site-directed mutagenesis was used to produce a panel of mutants, and their melting temperatures were determined by thermal denaturation. The results reveal a sequence fitness of an excess stability of approximately 12 °C, a value taken from a decrease in the melting temperature of an electrostatic charge-reversal substitution in the CRD3 without a deleterious effect on the binding affinity to the antigen. The tolerance for the disruption of antigen recognition without loss in the thermal stability was demonstrated by the introduction of a proline in place of a tyrosine in the CDR2, producing a mutant that eliminated binding. To further assist the sequence design and the selection of engineered single-domain antibodies, an assessment of different computational strategies is provided of their accuracy in the detection of substitution "hot spots" in the sequence tolerance landscape.

Spotting the difference in molecular dynamics simulations of biomolecules

Comparing two trajectories from molecular simulations conducted under different conditions is not a trivial task. In this study, we apply a method called Linear Discriminant Analysis with ITERative procedure (LDA-ITER) to compare two molecular simulation results by finding the appropriate projection vectors. Because LDA-ITER attempts to determine a projection such that the projections of the two trajectories do not overlap, the comparison does not suffer from a strong anisotropy, which is an issue in protein dynamics. LDA-ITER is applied to two test cases: the T4 lysozyme protein simulation with or without a point mutation and the allosteric protein PDZ2 domain of hPTP1E with or without a ligand. The projection determined by the method agrees with the experimental data and previous simulations. The proposed procedure, which complements existing methods, is a versatile analytical method that is specialized to find the "difference" between two trajectories.

A pH-dependent computational approach to the effect of mutations on protein stability

This article describes a novel software implementation for high-throughput scanning mutagenesis with a focus on protein stability. The approach combines molecular mechanics calculations with calculations of protein ionization and a Gaussian-chain model of electrostatic interactions in unfolded state. Comprehensive testing demonstrates a state-of-the-art accuracy for predicted free energy differences on single, double, and triple mutations with a correlation coefficient R above 0.7, which takes about 1.5 min per mutation on a single CPU. Unlike most of existing in silico methods for fast mutagenesis, the stability changes are reported as a continuous function of solution pH for wide pH intervals. We also propose a novel in silico strategy for searching stabilized protein variants that is based on combinatorial scanning mutagenesis using representative amino acid types. Our in silico predictions are in excellent agreement with the hyper-stabilized variants of mesophilic cold shock protein found using the Proside method of direct evolution. © 2016 Wiley Periodicals, Inc.

Extended Thermodynamic Integration: Efficient Prediction of Lambda Derivatives at Nonsimulated Points

Thermodynamic integration (TI) is one of the most commonly used free-energy calculation methods. The derivative of the Hamiltonian with respect to lambda, ⟨∂H/∂λ⟩, is determined at multiple λ-points. Because a numerical integration step is necessary, high curvature regions require simulations at densely spaced λ-points. Here, the principle of extended TI is introduced, where ⟨∂H/∂λ⟩ values are predicted at nonsimulated λ-points. On the basis of three model systems, it is shown that extended TI requires significantly fewer λ-points than regular TI to obtain similar accuracy.

Thermodynamic Cycle Without Turning Off Self-Interactions: Formal Discussion and a Numerical Example

The efficiency and accuracy of thermodynamic cycle calculations are considered. It is rigorously shown that the energy of the mutated part (MP) need not be scaled in a thermodynamic cycle computed with dual topology. Hence, there is no need to scale to zero any of the self-interactions (i.e. the interactions involving only particles of the same MP) regardless of whether the MP is bound or not to the main system. This observation carries a promise to lower computational resources and increase accuracy. A numerical test of a complete thermodynamic cycle illustrates cost and accuracy considerations.

Analyzing effects of naturally occurring missense mutations

Single-point mutation in genome, for example, single-nucleotide polymorphism (SNP) or rare genetic mutation, is the change of a single nucleotide for another in the genome sequence. Some of them will produce an amino acid substitution in the corresponding protein sequence (missense mutations); others will not. This paper focuses on genetic mutations resulting in a change in the amino acid sequence of the corresponding protein and how to assess their effects on protein wild-type characteristics. The existing methods and approaches for predicting the effects of mutation on protein stability, structure, and dynamics are outlined and discussed with respect to their underlying principles. Available resources, either as stand-alone applications or webservers, are pointed out as well. It is emphasized that understanding the molecular mechanisms behind these effects due to these missense mutations is of critical importance for detecting disease-causing mutations. The paper provides several examples of the application of 3D structure-based methods to model the effects of protein stability and protein-protein interactions caused by missense mutations as well.

Predicting folding free energy changes upon single point mutations

MOTIVATION

The folding free energy is an important characteristic of proteins stability and is directly related to protein's wild-type function. The changes of protein's stability due to naturally occurring mutations, missense mutations, are typically causing diseases. Single point mutations made in vitro are frequently used to assess the contribution of given amino acid to the stability of the protein. In both cases, it is desirable to predict the change of the folding free energy upon single point mutations in order to either provide insights of the molecular mechanism of the change or to design new experimental studies.

RESULTS

We report an approach that predicts the free energy change upon single point mutation by utilizing the 3D structure of the wild-type protein. It is based on variation of the molecular mechanics Generalized Born (MMGB) method, scaled with optimized parameters (sMMGB) and utilizing specific model of unfolded state. The corresponding mutations are built in silico and the predictions are tested against large dataset of 1109 mutations with experimentally measured changes of the folding free energy. Benchmarking resulted in root mean square deviation = 1.78 kcal/mol and slope of the linear regression fit between the experimental data and the calculations was 1.04. The sMMGB is compared with other leading methods of predicting folding free energy changes upon single mutations and results discussed with respect to various parameters.

AVAILABILITY

All the pdb files we used in this article can be downloaded from http://compbio.clemson.edu/downloadDir/mentaldisorders/sMMGB_pdb.rar.

SUPPLEMENTARY INFORMATION

Supplementary data are available at Bioinformatics online.

The linear interaction energy method for the prediction of protein stability changes upon mutation

The coupling of protein energetics and sequence changes is a critical aspect of computational protein design, as well as for the understanding of protein evolution, human disease, and drug resistance. To study the molecular basis for this coupling, computational tools must be sufficiently accurate and computationally inexpensive enough to handle large amounts of sequence data. We have developed a computational approach based on the linear interaction energy (LIE) approximation to predict the changes in the free-energy of the native state induced by a single mutation. This approach was applied to a set of 822 mutations in 10 proteins which resulted in an average unsigned error of 0.82 kcal/mol and a correlation coefficient of 0.72 between the calculated and experimental ΔΔG values. The method is able to accurately identify destabilizing hot spot mutations; however, it has difficulty in distinguishing between stabilizing and destabilizing mutations because of the distribution of stability changes for the set of mutations used to parameterize the model. In addition, the model also performs quite well in initial tests on a small set of double mutations. On the basis of these promising results, we can begin to examine the relationship between protein stability and fitness, correlated mutations, and drug resistance.

Evaluation at atomic resolution of the role of strain in destabilizing the temperature-sensitive T4 lysozyme mutant Arg 96 --> His

Mutant R96H is a classic temperature-sensitive mutant of bacteriophage T4 lysozyme. It was in fact the first variant of the protein to be characterized structurally. Subsequently, it has been studied extensively by a variety of experimental and computational techniques, but the reasons for the loss of stability of the mutant protein remain controversial. In the crystallographic refinement of the mutant structure at 1.9 A resolution one of the bond angles at the site of substitution appeared to be distorted by about 11( degrees ), and it was suggested that this steric strain was one of the major factors in destabilizing the mutant. Different computationally-derived models of the mutant structure, however, did not show such distortion. To determine the geometry at the site of mutation more reliably, we have extended the resolution of the data and refined the wildtype (WT) and mutant structures to be better than 1.1 A resolution. The high-resolution refinement of the structure of R96H does not support the bond angle distortion seen in the 1.9 A structure determination. At the same time, it does confirm other manifestations of strain seen previously including an unusual rotameric state for His96 with distorted hydrogen bonding. The rotamer strain has been estimated as about 0.8 kcal/mol, which is about 25% of the overall reduction in stability of the mutant. Because of concern that contacts from a neighboring molecule in the crystal might influence the geometry at the site of mutation we also constructed and analyzed supplemental mutant structures in which this crystal contact was eliminated. High-resolution refinement shows that the crystal contacts have essentially no effect on the conformation of Arg96 in WT or on His96 in the R96H mutant.

Contributions of all 20 amino acids at site 96 to the stability and structure of T4 lysozyme

To try to resolve the loss of stability in the temperature-sensitive mutant of T4 lysozyme, Arg 96 --> His, all of the remaining 18 naturally occurring amino acids were substituted at site 96. Also, in response to suggestions that the charged residues Lys85 and Asp89, which are 5-8 A away, may have important effects, each of these amino acids was replaced with alanine. Crystal structures were determined for many of the variants. With the exception of the tryptophan and valine mutants R96W and R96V, the crystallographic analysis shows that the substituted side chain following the path of Arg96 in wildtype (WT). The melting temperatures of the variants decrease by up to approximately 16 degrees C with WT being most stable. There are two site 96 replacements, with lysine or glutamine, that leave the stability close to that of WT. The only element that the side chains of these residues have in common with the WT arginine is the set of three carbon atoms at the C(alpha), C(beta), and C(gamma) positions. Although each side chain is long and flexible with a polar group at the distal position, the details of the hydrogen bonding to the rest of the protein differ in each case. Also, the glutamine replacement lacks a positive charge. This shows that there is some adaptability in achieving full stabilization at this site. At the other extreme, to be maximally destabilizing a mutation at site 96 must not only eliminate favorable interactions but also introduce an unfavorable element such as steric strain or a hydrogen-bonding group that remains unsatisfied. Overall, the study highlights the essential need for atomic resolution site-specific structural information to understand and to predict the stability of mutant proteins. It can be very misleading to simply assume that conservative amino acid substitutions cause small changes in stability, whereas large stability changes are associated with nonconservative replacements.

Fast and accurate predictions of protein stability changes upon mutations using statistical potentials and neural networks: PoPMuSiC-2.0

MOTIVATION

The rational design of proteins with modified properties, through amino acid substitutions, is of crucial importance in a large variety of applications. Given the huge number of possible substitutions, every protein engineering project would benefit strongly from the guidance of in silico methods able to predict rapidly, and with reasonable accuracy, the stability changes resulting from all possible mutations in a protein.

RESULTS

We exploit newly developed statistical potentials, based on a formalism that highlights the coupling between four protein sequence and structure descriptors, and take into account the amino acid volume variation upon mutation. The stability change is expressed as a linear combination of these energy functions, whose proportionality coefficients vary with the solvent accessibility of the mutated residue and are identified with the help of a neural network. A correlation coefficient of R = 0.63 and a root mean square error of sigma(c) = 1.15 kcal/mol between measured and predicted stability changes are obtained upon cross-validation. These scores reach R = 0.79, and sigma(c) = 0.86 kcal/mol after exclusion of 10% outliers. The predictive power of our method is shown to be significantly higher than that of other programs described in the literature.

AVAILABILITY

http://babylone.ulb.ac.be/popmusic

CHARMM: the biomolecular simulation program

CHARMM (Chemistry at HARvard Molecular Mechanics) is a highly versatile and widely used molecular simulation program. It has been developed over the last three decades with a primary focus on molecules of biological interest, including proteins, peptides, lipids, nucleic acids, carbohydrates, and small molecule ligands, as they occur in solution, crystals, and membrane environments. For the study of such systems, the program provides a large suite of computational tools that include numerous conformational and path sampling methods, free energy estimators, molecular minimization, dynamics, and analysis techniques, and model-building capabilities. The CHARMM program is applicable to problems involving a much broader class of many-particle systems. Calculations with CHARMM can be performed using a number of different energy functions and models, from mixed quantum mechanical-molecular mechanical force fields, to all-atom classical potential energy functions with explicit solvent and various boundary conditions, to implicit solvent and membrane models. The program has been ported to numerous platforms in both serial and parallel architectures. This article provides an overview of the program as it exists today with an emphasis on developments since the publication of the original CHARMM article in 1983.

Structural and functional implications of p53 missense cancer mutations

Most human cancers contain mutations in the transcription factor p53 and majority of these are missense and located in the DNA binding core domain. In this study, the stabilities of all core domain missense mutations are predicted and are used to infer their likely inactivation mechanisms. Overall, 47.0% non-PRO/GLY mutants are stable (DeltaDeltaG < 1.0 kT) and 36.3% mutants are unstable (DeltaDeltaG > 3.0 kT), 12.2% mutants are with 1.0 kT < DeltaDeltaG < 3.0 kT. Only 4.5% mutants are with no conclusive predictions. Certain types of either stable or unstable mutations are found not to depend on their local structures. Y, I, C, V, F and W (W, R and F) are the most common residues before (after) mutation in unstable mutants. Q, N, K, D, A, S and T (I, T, L and V) are the most common residues before (after) mutation in stable mutants. The stability correlations with sequence, structure, and molecular contacts are also analyzed. No direct correlation between secondary structure and stability is apparent, but a strong correlation between solvent exposure and stability is noticeable. Our correlation analysis shows that loss of protein-protein contacts may be an alternative cause for p53 inactivation. Correlation with clinical data shows that loss of stability and loss of DNA contacts are the two main inactivation mechanisms. Finally, correlation with functional data shows that most mutations which retain functions are stable, and most mutations that gain functions are unstable, indicating destabilized and deformed p53 proteins are more likely to find new binding partners.PACS codes: 87.14.E-

Reliable prediction of protein thermostability change upon double mutation from amino acid sequence

SUMMARY

The accurate prediction of protein stability change upon mutation is one of the important issues for protein design. In this work, we have focused on the stability change of double mutations and systematically analyzed the wild-type and mutant residues, patterns in amino acid sequence and locations of mutants. Based on the sequence information of wild-type, mutant and three neighboring residues, we have presented a weighted decision table method (WET) for predicting the stability changes of 180 double mutants obtained from thermal (DeltaDeltaG) denaturation. Using 10-fold cross-validation test, our method showed a correlation of 0.75 between experimental and predicted values of stability changes, and an accuracy of 82.2% for discriminating the stabilizing and destabilizing mutants.

Correlated motions and interactions at the onset of the DNA-induced partial unfolding of Ets-1

The binding of the Ets-1 transcription factor to its target DNA sequence is characterized by a highly unusual conformational change consisting of the unfolding of inhibitory helix 1 (HI-1). To probe the interactions that lead to this unfolding, we performed molecular dynamics simulations of the folded states of apo-Ets-1 and the Ets-1-DNA complex. The simulations showed large differences in correlated motions between helix 4 (H4) and HI-1. In apo-Ets-1, H4 and HI-1 moved in-phase and stabilized each other by hydrogen bonding and macrodipolar interactions, whereas in the DNA-bound state, the motion was out-of-phase, with a disruption of the stabilizing interactions. This change in motion was due to hydrogen-bonding interactions between helix 1 (H1) and the DNA. The dipolar energy between H1 and H4 was modulated by hydrogen bonds between H1 and DNA, and, in accordance with experiments, elimination of the hydrogen bonds increased the stability of HI-1. The simulations confirm that the hydrogen bonds between H1 and DNA act as a conformational switch and show that the presence of DNA is communicated from H1 to H4, destabilizing HI-1. The calculations reveal a critical role for correlated motions at the onset of the DNA-induced unfolding.

The PT1-Ca2+ Gla domain binds to a membrane through two dipalmitoylphosphatidylserines. A computational study

Binding of vitamin K-dependent proteins to cell membranes containing phosphatidylserine (PS) via gamma-carboxyglutamic acid (Gla) domains is one of the essential steps in the blood coagulation pathway. During activation of the coagulation cascade, prothrombin is converted to thrombin by prothrombinase, a complex consisting of serine protease FXa and cofactor FVa, anchored to anionic phospholipids on the surface of activated platelets in the presence of calcium ions. To investigate the binding of the Gla domain of prothrombin fragment 1 (PT1) to anionic lipids in the presence of Ca2+, we have conducted MD simulations of the protein with one and two dipalmitoylphosphatidylserines (DPPS) in a dipalmitoylphosphatidylcholine (DPPC) bilayer membrane. The results show a well-defined phosphatidylserine binding site, which agrees generally with crystallographic studies [Huang, M., et al. (2003) Nat. Struct. Biol. 10, 751-756]. However, in the presence of the lipid membrane, some of the interactions observed in the crystal structure adjust during the simulations possibly because in our system the PT1-Ca2+ complex is embedded in a DPPC lipid membrane. Our simulations confirm the existence of a second phospholipid headgroup binding site on the opposite face of the PT1-Ca2+ complex as suggested by MacDonald et al. [(1997) Biochemistry 36, 5120-5127]. The serine headgroup in the second site binds through a Gla domain-bound calcium ion Ca1, Gla30, and Lys11. On the basis of free energy simulations, we estimate the energy of binding of the PT1-Ca2+ complex to a single DPPS to be around -11.5 kcal/mol. The estimated free energy of binding of a DPPS lipid to the second binding site is around -8.8 kcal/mol and is in part caused by the nature of the second site and in part by entropic effects.

Simulation analysis of the stability mutants R96H of bacteriophage T4 lysozyme and I96A of barnase

Free energy simulation methods are used to analyse the effects of the mutation Arg-96----His on the stability of bacteriophage T4 lysozyme and of Ile-96----Ala on the stability of barnase. By use of thermodynamic integration, the contributions of specific interactions to the free energy change are evaluated. It is shown that a number of contributions that stabilize the wild-type or the mutant partially cancel in the overall free energy difference; some of these involve the unfolded state. Comparison of the results with conclusions based on structural and thermodynamic data leads to new insights into the origin of the stability difference between wild-type and mutant proteins. For the charged-to-charged amino acid mutation in T4 lysozyme, the importance of the contributions of more distant residues, solvent water and the covalent linkage involving the mutated amino acid are of particular interest. Also, the analysis of the Arg-96 to His mutation with respect to the interactions with the C-terminal end of a helix (residues 82-90) indicates that the nearby carbonyl groups (Tyr-88 and Asp-89) make the dominant contribution, that the amide groups do not contribute significantly and that the helix dipole model is inappropriate for this case. For the non-polar-to-non-polar amino acid mutation in barnase, the solvent contribution is unimportant, and covalent terms are shown to be significant because they do not cancel between the folded and unfolded state.

Ligand design by a combinatorial approach based on modeling and experiment: application to HLA-DR4

Combinatorial synthesis and large scale screening methods are being used increasingly in drug discovery, particularly for finding novel lead compounds. Although these "random" methods sample larger areas of chemical space than traditional synthetic approaches, only a relatively small percentage of all possible compounds are practically accessible. It is therefore helpful to select regions of chemical space that have greater likelihood of yielding useful leads. When three-dimensional structural data are available for the target molecule this can be achieved by applying structure-based computational design methods to focus the combinatorial library. This is advantageous over the standard usage of computational methods to design a small number of specific novel ligands, because here computation is employed as part of the combinatorial design process and so is required only to determine a propensity for binding of certain chemical moieties in regions of the target molecule. This paper describes the application of the Multiple Copy Simultaneous Search (MCSS) method, an active site mapping and de novo structure-based design tool, to design a focused combinatorial library for the class II MHC protein HLA-DR4. Methods for the synthesizing and screening the computationally designed library are presented; evidence is provided to show that binding was achieved. Although the structure of the protein-ligand complex could not be determined, experimental results including cross-exclusion of a known HLA-DR4 peptide ligand (HA) by a compound from the library. Computational model building suggest that at least one of the ligands designed and identified by the methods described binds in a mode similar to that of native peptides.

Association free energy of dipalmitoylphosphatidylserines in a mixed dipalmitoylphosphatidylcholine membrane

Blood coagulation is strongly dependent on the binding of vitamin K-dependent proteins to cell membranes containing phosphatidylserine (PS) via gamma-carboxyglutamic acid (Gla) domains. The process depends on calcium, which can induce nonideal behavior in membranes through domain formation. Such domain separation mediated by Ca(2+) ions or proteins can have an important contribution to the thermodynamics of the interaction between charged peripheral proteins and oppositely charged membranes. To characterize the properties of lipid-lipid interactions, molecular dynamics, and free energy simulations in a mixed bilayer membrane containing dipalmitoylphosphatidylcholine and dipalmitoylphosphatidylserine were carried out. The free energy of association between dipalmitoylphosphatidylserines in the environment of dipalmitoylphosphatidylcholines has been calculated by using a novel approach to the dual topology technique of the PS-PC hybrid. Two different methods, free energy perturbation and thermodynamic integration, were used to calculate the free energy difference. In thermodynamic integration runs three schemes were applied to evaluate the integral at the limits of lambda --> 0 or lambda --> 1. Our studies show that the association of two PSs in the environment of PCs is repulsive in the absence of Ca(2+) and becomes favorable in their presence. We also show that the mixed component membrane should exhibit nonideal behavior that will lead to PS clustering.

Calculation of the aqueous solvation energy and entropy, as well as free energy, of simple polar solutes

With the advent of more powerful computers, the question of calculating thermodynamic quantities, such as the energy and the entropy, in solute-solvent systems is revisited. The calculation of these thermodynamic quantitites was limited in the past by their slow convergence relative to the free energy. Using molecular dynamics simulations, the energy, entropy, and free energy of solvation of NMA and CH(3)NH(2), as well as their relative values, have been determined. Three different methods (the thermodynamic perturbation method, the thermodynamic integration method, and a finite-difference method) are compared. The finite difference method gives the best results and accurate values for the entropy and energy were obtained using a reasonable amount to computer time. The results suggest that a meaningful thermodynamic description of biomolecular processes can be realized with present methods and the available computer time.

Average assignment method for predicting the stability of protein mutants

Prediction of protein stability upon amino acid substitutions is an important problem in molecular biology and it will be helpful for designing stable mutants. In this work, we have analyzed the stability of protein mutants using three different data sets of 1791, 1396, and 2204 mutants, respectively, for thermal stability (DeltaTm), free energy change due to thermal (DeltaDeltaG), and denaturant denaturations (DeltaDeltaGH2O), obtained from the ProTherm database. We have classified the mutants into 380 possible substitutions and assigned the stability of each mutant using the information obtained with similar type of mutations. We observed that this assignment could distinguish the stabilizing and destabilizing mutants to an accuracy of 70-80% at different measures of stability. Further, we have classified the mutants based on secondary structure and solvent accessibility (ASA) and observed that the classification significantly improved the accuracy of prediction. The classification of mutants based on helix, strand, and coil distinguished the stabilizing/destabilizing mutants at an average accuracy of 82% and the correlation is 0.56; information about the location of residues at the interior, partially buried, and surface regions of a protein correctly identified the stabilizing/destabilizing residues at an average accuracy of 81% and the correlation is 0.59. The nine subclassifications based on three secondary structures and solvent accessibilities improved the accuracy of assigning stabilizing/destabilizing mutants to an accuracy of 84-89% for the three data sets. Further, the present method is able to predict the free energy change (DeltaDeltaG) upon mutations within a deviation of 0.64 kcal/mol. We suggest that this method could be used for predicting the stability of protein mutants.

Prediction of protein stability changes for single-site mutations using support vector machines

Accurate prediction of protein stability changes resulting from single amino acid mutations is important for understanding protein structures and designing new proteins. We use support vector machines to predict protein stability changes for single amino acid mutations leveraging both sequence and structural information. We evaluate our approach using cross-validation methods on a large dataset of single amino acid mutations. When only the sign of the stability changes is considered, the predictive method achieves 84% accuracy-a significant improvement over previously published results. Moreover, the experimental results show that the prediction accuracy obtained using sequence alone is close to the accuracy obtained using tertiary structure information. Because our method can accurately predict protein stability changes using primary sequence information only, it is applicable to many situations where the tertiary structure is unknown, overcoming a major limitation of previous methods which require tertiary information. The web server for predictions of protein stability changes upon mutations (MUpro), software, and datasets are available at http://www.igb.uci.edu/servers/servers.html.

Glycogen phosphorylase inhibitors: a free energy perturbation analysis of glucopyranose spirohydantoin analogues

GP catalyzes the phosphorylation of glycogen to Glc-1-P. Because of its fundamental role in the metabolism of glycogen, GP has been the target for a systematic structure-assisted design of inhibitory compounds, which could be of value in the therapeutic treatment of type 2 diabetes mellitus. The most potent catalytic-site inhibitor of GP identified to date is spirohydantoin of glucopyranose (hydan). In this work, we employ MD free energy simulations to calculate the relative binding affinities for GP of hydan and two spirohydantoin analogues, methyl-hydan and n-hydan, in which a hydrogen atom is replaced by a methyl- or amino group, respectively. The results are compared with the experimental relative affinities of these ligands, estimated by kinetic measurements of the ligand inhibition constants. The calculated binding affinity for methyl-hydan (relative to hydan) is 3.75 +/- 1.4 kcal/mol, in excellent agreement with the experimental value (3.6 +/- 0.2 kcal/mol). For n-hydan, the calculated value is 1.0 +/- 1.1 kcal/mol, somewhat smaller than the experimental result (2.3 +/- 0.1 kcal/mol). A free energy decomposition analysis shows that hydan makes optimum interactions with protein residues and specific water molecules in the catalytic site. In the other two ligands, structural perturbations of the active site by the additional methyl- or amino group reduce the corresponding binding affinities. The computed binding free energies are sensitive to the preference of a specific water molecule for two well-defined positions in the catalytic site. The behavior of this water is analyzed in detail, and the free energy profile for the translocation of the water between the two positions is evaluated. The results provide insights into the role of water molecules in modulating ligand binding affinities. A comparison of the interactions between a set of ligands and their surrounding groups in X-ray structures is often used in the interpretation of binding free energy differences and in guiding the design of new ligands. For the systems in this work, such an approach fails to estimate the order of relative binding strengths, in contrast to the rigorous free energy treatment.

The Allosteric Mechanism of Yeast Chorismate Mutase: A Dynamic Analysis

The effector-regulated allosteric mechanism of yeast chorismate mutase (YCM) was studied by normal mode analysis and targeted molecular dynamics. The normal mode analysis shows that the conformational change between YCM in the R state and in the T state can be represented by a relatively small number of low-frequency modes. This suggests that the transition is coded in the structure and is likely to have a low energetic barrier. Quantitative comparisons (i.e. frequencies) between the low-frequency modes of YCM with and without effectors (modeled structures) reveal that the binding of Trp increases the global flexibility, whereas Tyr decreases global flexibility. The targeted molecular dynamics simulation of substrate analog release from the YCM active site suggests that a series of residues are critical for orienting and "recruiting" the substrate. The simulation led to the switching of a series of substrate-release-coupled salt-bridge partners in the ligand-binding domain; similar changes occur in the transition between YCM R-state and T-state crystal structures. Thus, the normal mode analysis and targeted molecular dynamics results provide evidence that the effectors regulate YCM activity by influencing the global flexibility. The change in flexibility is coupled to the binding of substrate to the T state and release of the product from the R state, respectively.

Molecular basis of peptide recognition by the TCR: affinity differences calculated using large scale computing

Free energy calculations of the wild-type and the variant human T cell lymphotropic virus type 1 Tax peptide presented by the MHC to the TCR have been performed using large scale massively parallel molecular dynamics simulations. The computed free energy difference (-1.86 +/- 0.44 kcal/mol) using alchemical mutation-based thermodynamic integration agrees well with experimental data (-2.9 +/- 0.2 kcal/mol). Our simulations exploit state-of-the-art hardware and codes whose algorithms have been optimized for supercomputing platforms. This enables us to simulate larger, more realistic biological systems for longer durations without the imposition of artificial constraints.

Peptide recognition by the T cell receptor: comparison of binding free energies from thermodynamic integration, Poisson-Boltzmann and linear interaction energy approximations

The binding to the T cell receptor of wild-type and variant HTLV-1 Tax peptide complexed to the major histocompatibility complex has been investigated by means of molecular dynamics simulations. The binding free energy difference is calculated using the molecular mechanics Poisson-Boltzmann surface area and linear interaction energy methods. These methods extract useful information on the binding energetics from simulations of the physical states of the ligands, which are more computationally expedient than the commonly used thermodynamic integration method. The successful reproduction of the relative binding free energies shows that these methods can be useful for free energy calculations and the rational design of drugs and vaccines.

Structural stability of wild type and mutated alpha-keratin fragments: molecular dynamics and free energy calculations

The objective of this study is to investigate the influence of point mutations on the structural stability of coiled coil fragments of the human hair intermediate filament by molecular dynamics simulations and free energy calculations. Mutations in the helix termination motif of human hair keratin gene hHb6 seem to be connected to the hereditary hair dystrophy Monilethrix. The most common mutations reported are Glu413Lys and Glu413Asp, located at the C-terminal end of the coiled coil 2B rod domain of the IF. According to our simulations, significant conformational changes of the side chains at the mutation and neighboring sites occur due to the Glu413Lys mutation. Furthermore, the differences in electrostatic interactions cause a large change in free energy during transformation of Glu413 to Lys calculated by the thermodynamic integration approach. It is speculated that the structural rearrangement necessary to adapt the interactions in the mutated coiled coil leads to changes in the IF assembly or its stability. The second mutation, Glu413Asp, only leads to a small value of the calculated free energy difference that is within the error limits of the simulations. Thus, it has to be concluded that this mutation does not affect the coiled coil stability.

Relative importance of secondary structure and solvent accessibility to the stability of protein mutants

Understanding the factors influencing the stability of protein mutants is an important task in molecular and computational biology. In this work, we have approached this problem by examining the relative importance of secondary structure and solvent accessibility of the mutant residue for understanding/predicting the stability of protein mutants. We have used hydrophobic, electrostatic and hydrogen bond free energy terms and nine unique physicochemical, energetic and conformational properties of amino acids in the present study and these parameters have been related with changes in thermal stability (DeltaTm) of all the single mutants of lysozymes based on single and multiple correlation coefficients. As expected the properties reflecting hydrophobicity and hydrophobic free energy play a major role to distinguish stabilizing and destabilizing mutants. The hydrophobic free energy due to carbon and nitrogen atoms distinguish the stability of coil and strand mutations to the accuracy of 100 and 90%, respectively. In agreement with previous results, the subgroup classification based on secondary structure and the information about its location in the structure yielded good relationship with the experimental DeltaTm. We revealed that the secondary structure information is equally or more important than solvent accessibility for understanding the stability of protein mutants. The comparison of amino acid properties with free-energy terms indicate that the energetic contribution explains the mutant stability better in coil region whereas the amino acid properties do better in strand region. Further, the combination of free energies with amino acid properties increased the correlation significantly. The present study demonstrates the importance of classifying the mutants based on secondary structure to the stability of proteins upon mutations.

X-ray structural and simulation analysis of a protein mutant: the value of a combined approach

The effect of the mutation Arg 96 to His on the stability of bacteriophage T4 lysozyme has been previously studied by calorimetric experiments, X-ray crystallography, and free energy simulation techniques. The experimental and calculated values for the difference between the free energy of denaturation of the mutant and the wild type are in reasonable agreement. However, the two approaches led to different explanations for the loss in stability. To analyze the differences, a series of refinements based on the crystallographic data were performed, a number of aspects of the simulations were reexamined, and continuum electrostatic calculations were done to complement the latter. The results of those comparisons provide a better understanding of the origin of the free energy difference in this mutant. Furthermore, they show the importance of the combined use of simulations and crystallography for interpreting the effects of mutations on the energetics of the system.

Large-scale prediction of protein geometry and stability changes for arbitrary single point mutations

We have developed a method to both predict the geometry and the relative stability of point mutants that may be used for arbitrary mutations. The geometry optimization procedure was first tested on a new benchmark of 2141 ordered pairs of X-ray crystal structures of proteins that differ by a single point mutation, the largest data set to date. An empirical energy function, which includes terms representing the energy contributions of the folded and denatured proteins and uses the predicted mutant side chain conformation, was fit to a training set consisting of half of a diverse set of 1816 experimental stability values for single point mutations in 81 different proteins. The data included a substantial number of small to large residue mutations not considered by previous prediction studies. After removing 22 (approximately 2%) outliers, the stability calculation gave a standard deviation of 1.08 kcal/mol with a correlation coefficient of 0.82. The prediction method was then tested on the remaining half of the experimental data, giving a standard deviation of 1.10 kcal/mol and covariance of 0.66 for 97% of the test set. A regression fit of the energy function to a subset of 137 mutants, for which both native and mutant structures were available, gave a prediction error comparable to that for the complete training set with predicted side chain conformations. We found that about half of the variation is due to conformation-independent residue contributions. Finally, a fit to the experimental stability data using these residue parameters exclusively suggests guidelines for improving protein stability in the absence of detailed structure information.

The missing link between thermodynamics and structure in F1-ATPase

F(1)F(o)-ATP synthase is the enzyme responsible for most of the ATP synthesis in living systems. The catalytic domain F(1) of the F(1)F(o) complex, F(1)-ATPase, has the ability to hydrolyze ATP. A fundamental problem in the development of a detailed mechanism for this enzyme is that it has not been possible to determine experimentally the relation between the ligand binding affinities measured in solution and the different conformations of the catalytic beta subunits (beta(TP), beta(DP), beta(E)) observed in the crystal structures of the mitochondrial enzyme, MF(1). Using free energy difference simulations for the hydrolysis reaction ATP+H(2)O --> ADP+P(i) in the beta(TP) and beta(DP) sites and unisite hydrolysis data, we are able to identify beta(TP) as the "tight" (K(D) = 10(-12) M, MF(1)) binding site for ATP and beta(DP) as the "loose" site. An energy decomposition analysis demonstrates how certain residues, some of which have been shown to be important in catalysis, modulate the free energy of the hydrolysis reaction in the beta(TP) and beta(DP) sites, even though their structures are very similar. Combined with the recently published simulations of the rotation cycle of F(1)-ATPase, the present results make possible a consistent description of the binding change mechanism of F(1)-ATPase at an atomic level of detail.

Binding free energy differences in a TCR-peptide-MHC complex induced by a peptide mutation: a simulation analysis

Recognition by the T-cell receptor (TCR) of immunogenic peptides presented by class I major histocompatibility complexes (MHCs) is the determining event in the specific cellular immune response against virus-infected cells or tumor cells. It is of great interest, therefore, to elucidate the molecular principles upon which the selectivity of a TCR is based. These principles can in turn be used to design therapeutic approaches, such as peptide-based immunotherapies of cancer. In this study, free energy simulation methods are used to analyze the binding free energy difference of a particular TCR (A6) for a wild-type peptide (Tax) and a mutant peptide (Tax P6A), both presented in HLA A2. The computed free energy difference is 2.9 kcal/mol, in good agreement with the experimental value. This makes possible the use of the simulation results for obtaining an understanding of the origin of the free energy difference which was not available from the experimental results. A free energy component analysis makes possible the decomposition of the free energy difference between the binding of the wild-type and mutant peptide into its components. Of particular interest is the fact that better solvation of the mutant peptide when bound to the MHC molecule is an important contribution to the greater affinity of the TCR for the latter. The results make possible identification of the residues of the TCR which are important for the selectivity. This provides an understanding of the molecular principles that govern the recognition. The possibility of using free energy simulations in designing peptide derivatives for cancer immunotherapy is briefly discussed.

Free-energy simulations of the oxidation of c-terminal methionines in calmodulin

In the course of aging or under conditions of oxidative stress, methionine residues of calmodulin undergo oxidation, leading to loss of biological activity of the protein. We have performed free-energy simulations of the effects of C-terminal methionine side-chain oxidation on the properties of calmodulin. The simulation results indicate that oxidation should have a destabilizing effect on all three protein functional states: calcium free, calcium loaded, and with both calcium and target peptide bound. Because the different states are destabilized by different amounts, this leads to a more complex pattern in the observable effects on protein thermal stability, calcium affinity, and binding of a target peptide. The influence of oxidation on the free energy of CaM unfolding is estimated by comparing the free-energy cost of oxidizing a Met residue in a Gly-Met-Gly peptide and in the protein. The protein thermal stability of the oxidized forms is lowered by a moderate amount 1-3 kcal/mol, in qualitative agreement with experimental results of 0.3 kcal/mol. The calculated changes in affinity for calcium and for the target peptide show opposing trends. Oxidation at position 144 is predicted to enhance peptide binding and weaken calcium binding, whereas oxidation at 145 weakens peptide binding and enhances affinity for calcium. The lower affinity of Met 145-oxidized calmodulin toward the target peptide correlates with experimentally observed lowering of calmodulin-activated Ca-ATPase activity when oxidized calmodulin from aged rat brains is used. Thus, our simulations suggest that Met 145 is the oxidation site in the C-terminal fragment of calmodulin. The microscopic mechanism behind the calculated free energy changes appears to be a greater affinity for water of the oxidized Met side-chain relative to normal Met. Structures with Met exposed to solvent had consistently lower free energies than those with buried Met sidechains.

Free energy simulations come of age: protein-ligand recognition

In recent years, molecular dynamics simulations of biomolecular free energy differences have benefited from significant methodological advances and increased computer power. Applications to molecular recognition provide an understanding of the interactions involved that goes beyond, and is an important complement to, experimental studies. Poisson-Boltzmann electrostatic models provide a faster and simpler free energy method in cases where electrostatic interactions are important. We illustrate both molecular dynamics and Poisson-Boltzmann methods with a detailed study of amino acid recognition by aspartyl-tRNA synthetase, whose specificity is important for maintaining the integrity of the genetic code.

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.