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.

Toward Atomistic Modeling of Irreversible Covalent Inhibitor Binding Kinetics.

Covalent inhibitors have emerged as an important drug class in recent years, largely due to their many unique advantages as compared to noncovalent inhibitors, including longer duration of action, lower prolonged systemic exposure, higher potency, and selectivity. However, the potential off-target toxicity of covalent inhibitors, particularly of irreversible covalent inhibitors, represents a great challenge in covalent drug development. Therefore, accurate calculation of protein covalent inhibitor reaction kinetics to guide the design of selective inhibitors would greatly benefit covalent drug discovery efforts. In the present paper, we present a computational method to calculate the relative reaction kinetics between congeneric irreversible covalent inhibitors and their protein receptors. The method combines density functional theory calculations of the transition state barrier height of the rate-limiting step for reaction between the warhead of the inhibitor and a single protein residue, and molecular-mechanics-based free energy calculations to account for the interactions between the ligand in the transition state and the protein environment. The method was tested on four pharmaceutically interesting irreversible covalent binding systems involving 28 ligands; the mean unsigned error (MUE) of the relative reaction rate for all pairs of ligands between the predictions and experimental results for these tested systems is 0.79 log unit. This is to our knowledge the first time where the reaction kinetics of protein irreversible covalent inhibition have been directly calculated with physics-based free energy calculation methods and transition state theory. We anticipate the outstanding accuracy demonstrated here across a broad range of target classes will have a strong impact on the design of selective covalent inhibitors.

Accurate and reliable prediction of relative ligand binding potency in prospective drug discovery by way of a modern free-energy calculation protocol and force field.

Designing tight-binding ligands is a primary objective of small-molecule drug discovery. Over the past few decades, free-energy calculations have benefited from improved force fields and sampling algorithms, as well as the advent of low-cost parallel computing. However, it has proven to be challenging to reliably achieve the level of accuracy that would be needed to guide lead optimization (∼5× in binding affinity) for a wide range of ligands and protein targets. Not surprisingly, widespread commercial application of free-energy simulations has been limited due to the lack of large-scale validation coupled with the technical challenges traditionally associated with running these types of calculations. Here, we report an approach that achieves an unprecedented level of accuracy across a broad range of target classes and ligands, with retrospective results encompassing 200 ligands and a wide variety of chemical perturbations, many of which involve significant changes in ligand chemical structures. In addition, we have applied the method in prospective drug discovery projects and found a significant improvement in the quality of the compounds synthesized that have been predicted to be potent. Compounds predicted to be potent by this approach have a substantial reduction in false positives relative to compounds synthesized on the basis of other computational or medicinal chemistry approaches. Furthermore, the results are consistent with those obtained from our retrospective studies, demonstrating the robustness and broad range of applicability of this approach, which can be used to drive decisions in lead optimization.

Target-template relationships in protein structure prediction and their effect on the accuracy of thermostability calculations.

Improving protein thermostability has been a labor- and time-consuming process in industrial applications of protein engineering. Advances in computational approaches have facilitated the development of more efficient strategies to allow the prioritization of stabilizing mutants. Among these is FEP+, a free energy perturbation implementation that uses a thoroughly tested physics-based method to achieve unparalleled accuracy in predicting changes in protein thermostability. To gauge the applicability of FEP+ to situations where crystal structures are unavailable, here we have applied the FEP+ approach to homology models of 12 different proteins covering 316 mutations. By comparing predictions obtained with homology models to those obtained using crystal structures, we have identified that local rather than global sequence conservation between target and template sequence is a determining factor in the accuracy of predictions. By excluding mutation sites with low local sequence identity (<40%) to a template structure, we have obtained predictions with comparable performance to crystal structures (R2 of 0.67 and 0.63 and an RMSE of 1.20 and 1.16 kcal/mol for crystal structure and homology model predictions, respectively) for identifying stabilizing mutations when incorporating residue scanning into a cascade screening strategy. Additionally, we identify and discuss inherent limitations in sequence alignments and homology modeling protocols that translate into the poor FEP+ performance of a few select examples. Overall, our retrospective study provides detailed guidelines for the application of the FEP+ approach using homology models for protein thermostability predictions, which will greatly extend this approach to studies that were previously limited by structure availability.

Close intramolecular sulfur-oxygen contacts: modified force field parameters for improved conformation generation.

The Cambridge Structural Database (CSD) offers an excellent data source to study small molecule conformations and molecular interactions. We have analyzed 130 small molecules from the CSD containing an intramolecular sulfur-oxygen distance less than the sum of their van der Waals (vdW) radii. Close S···O distances are observed in several important medicinal chemistry motifs (e.g. a carbonyl oxygen connected by a carbon or nitrogen linker to a sulfur) and are not treated well with existing parameters in the MMFFs or OPLS_2005 force fields, resulting in suboptimal geometries and energetics. In this work, we develop modified parameters for the OPLS_2005 force field to better treat this specific interaction in order to generate conformations close to those found in the CSD structures. We use a combination of refitting a force field torsional parameter, adding a specific atom pair vdW term, and attenuating the electrostatic interactions to obtain an improvement in the accuracy of geometry minimizations and conformational searches for these molecules. Specifically, in a conformational search 58 % of the cases produced a conformation less than 0.25 Å from the CSD crystal conformation with the modified OPLS force field parameters developed in this work. In contrast, 25 and 37 % produced a conformation less than 0.25 Å with the MMFFs and OPLS_2005 force fields, respectively. As an application of the new parameters, we generated conformations for the tyrosine kinase inhibitor axitinib (trade name Inlyta) that could be correctly repacked into three observed polymorphic structures, which was not possible with conformations generated using MMFFs or OPLS_2005. The improved parameters can be mapped directly onto physical characteristics of the systems that are treated inadequately with the molecular mechanics force fields used in this study and potentially other force fields as well.

A molecular dynamics investigation of lipid bilayer perturbation by PIP2.

Phosphoinositides like phosphatidylinositol 4,5-bisphosphate (PIP(2)) are negatively charged lipids that play a pivotal role in membrane trafficking, signal transduction, and protein anchoring. We have designed a force field for the PIP(2) headgroup using quantum mechanical methods and characterized its properties inside a lipid bilayer using molecular dynamics simulations. Macroscopic properties such as area/headgroup, density profiles, and lipid order parameters calculated from these simulations agree well with the experimental values. However, microscopically, the PIP(2) introduces a local perturbation of the lipid bilayer. The average PIP(2) headgroup orientation of 45 degrees relative to the bilayer normal induces a unique, distance-dependent organization of the lipids that surround PIP(2). The headgroups of these lipids preferentially orient closer to the bilayer normal. This perturbation creates a PIP(2) lipid microdomain with the neighboring lipids. We propose that the PIP(2) lipid microdomain enables the PIP(2) to function as a membrane-bound anchoring molecule.

OPLS3: A Force Field Providing Broad Coverage of Drug-like Small Molecules and Proteins.

The parametrization and validation of the OPLS3 force field for small molecules and proteins are reported. Enhancements with respect to the previous version (OPLS2.1) include the addition of off-atom charge sites to represent halogen bonding and aryl nitrogen lone pairs as well as a complete refit of peptide dihedral parameters to better model the native structure of proteins. To adequately cover medicinal chemical space, OPLS3 employs over an order of magnitude more reference data and associated parameter types relative to other commonly used small molecule force fields (e.g., MMFF and OPLS_2005). As a consequence, OPLS3 achieves a high level of accuracy across performance benchmarks that assess small molecule conformational propensities and solvation. The newly fitted peptide dihedrals lead to significant improvements in the representation of secondary structure elements in simulated peptides and native structure stability over a number of proteins. Together, the improvements made to both the small molecule and protein force field lead to a high level of accuracy in predicting protein-ligand binding measured over a wide range of targets and ligands (less than 1 kcal/mol RMS error) representing a 30% improvement over earlier variants of the OPLS force field.

Core deformations in protein families: a physical perspective.

An analysis is presented on how structural cores change shape within protein families, and whether or not there is a relationship between these structural changes and the vibrational modes that proteins experiment due to topological constraints. A set of 13 representative and well-populated protein families are studied. The evolutionary directions of deformation are obtained by applying a new multiple structural alignment technique to superimpose the structures and extract a conserved core, together with Principal Components Analysis (PCA) to extract the main deformation modes. A low-resolution Normal Mode Analysis (NMA) technique is used in parallel to study the properties of the mechanical core plasticity of the same proteins. We find that the evolutionary deformations span a low dimensional space. A statistically significant correspondence exists between these principal deformations and the vibrational modes accessible to a particular topology. We conclude that, to a significant extent, the structures of evolving proteins seem to respond to sequence changes by collective deformations along combinations of low-frequency modes. The findings have implications in structure prediction by homology modeling.

A role for a specific cholesterol interaction in stabilizing the Apo configuration of the human A(2A) adenosine receptor.

The function of G-protein-coupled receptors is tightly modulated by the lipid environment. Long-timescale molecular dynamics simulations (totaling approximately 3 mus) of the A(2A) receptor in cholesterol-free bilayers, with and without the antagonist ZM241385 bound, demonstrate the instability of helix II in the apo receptor in cholesterol-poor membrane regions. We directly observe that the effect of cholesterol binding is to stabilize helix II against a buckling-type deformation, perhaps rationalizing the observation that the A(2A) receptor couples to G protein only in the presence of cholesterol (Zezula and Freissmuth, 2008). The results suggest a mechanism by which the A(2A) receptor may function as a coincidence detector, activating only in the presence of both cholesterol and agonist. We also observed a previously hypothesized conformation of the tryptophan "rotameric switch" on helix VI in which a phenylalanine on helix V positions the tryptophan out of the ligand binding pocket.

Effective mast cell degranulating peptide inhibitors of the IgE/Fc epsilonRI receptor interaction.

Previous studies with mast cell degranulating (MCD) peptide have shown that peptide [Ala(12)]MCD 8 was an inhibitor of IgE binding to mast cell receptors. In an attempt to produce increased inhibition, analogs were synthesized that maintained the alanine residue in position 12 in the MCD peptide sequence and were further modified at both termini. Analogs modified at the C-terminus were [Ala(12),desLys(21)]MCD 2 and [Ala(12),D-Lys(21)]MCD 4. N-terminus modifications were [desLys(6)-Arg(7)-His(8),Ala(12)]MCD 1, [Ala(6), Ala(12)]MCD 6, and [Val(6),Ala(12)]MCD 7. To assess the role of the Proline(12), analogs [D-Ala(12)]MCD 3 and [Meleu(12)]MCD 5 were also synthesized. The analogs were tested for binding to the IgE receptor in cultured mast cells. Inhibitory activity of IgE-caused degranulation was measured using a beta-hexosaminidase assay. Circular dichroism (CD) and molecular modeling of selected analogs were used to follow possible structural differences among these analogs. All analogs showed binding affinity to the IgE receptor and inhibition of IgE-induced mast cell degranulation at different levels. Differences in inhibition were most likely because of diverse interactions of the analogs with the receptor as inferred by the CD and modeling studies. Based on the results of the beta-hexosaminidase assay, analog [Val(6), Ala(12)]MCD 7 proved to be an excellent inhibitor of IgE-mediated mast cell degranulation.

Diverse Kir modulators act in close proximity to residues implicated in phosphoinositide binding.

Inwardly rectifying potassium (Kir) channels were the first shown to be directly activated by phosphoinositides in general and phosphatidylinositol bisphosphate (PIP(2)) in particular. Atomic resolution structures have been determined for several mammalian and bacterial Kir channels. Basic residues, identified through mutagenesis studies to contribute to the sensitivity of the channel to PIP(2), have been mapped onto the three dimensional channel structure and their localization has given rise to a plausible model that can explain channel activation by PIP(2). Moreover, mapping onto the three-dimensional channel structure sites involved in the modulation of Kir channel activity by a diverse group of regulatory molecules, revealed a striking proximity to residues implicated in phosphoinositide binding. These observations support the hypothesis that the observed dependence of diverse modulators on channel-PIP(2) interactions stems from their localization within distances that can affect PIP(2)-interacting residues.

Phosphoinositide-mediated gating of inwardly rectifying K(+) channels.

Phosphoinositides, such as phosphatidylinositol-bisphosphate (PIP(2)), control the activity of many ion channels in yet undefined ways. Inwardly, rectifying potassium (Kir) channels were the first shown to be dependent on direct interactions with phosphoinositides. Alterations in channel-PIP(2) interactions affect Kir single-channel gating behavior. Aberrations in channel-PIP(2) interactions can lead to human disease. As the activity of all Kir channels depends on their interactions with phosphoinositides, future research will aim to understand the molecular events that occur from phosphoinositide binding to channel gating. The determination of atomic resolution structures for several mammalian and bacterial Kir channels provides great promise towards this goal. We have mapped onto the three-dimensional channel structure the position of basic residues identified through mutagenesis studies that contribute to the sensitivity of a Kir channel to PIP(2). The localization of these putative PIP(2)-interacting residues relative to the channel's permeation pathway has given rise to a testable model, which could account for channel activation by PIP(2).

A new progressive-iterative algorithm for multiple structure alignment.

MOTIVATION: Multiple structure alignments are becoming important tools in many aspects of structural bioinformatics. The current explosion in the number of available protein structures demands multiple structural alignment algorithms with an adequate balance of accuracy and speed, for large scale applications in structural genomics, protein structure prediction and protein classification. RESULTS: A new multiple structural alignment program, MAMMOTH-mult, is described. It is demonstrated that the alignments obtained with the new method are an improvement over previous manual or automatic alignments available in several widely used databases at all structural levels. Detailed analysis of the structural alignments for a few representative cases indicates that MAMMOTH-mult delivers biologically meaningful trees and conservation at the sequence and structural levels of functional motifs in the alignments. An important improvement over previous methods is the reduction in computational cost. Typical alignments take only a median time of 5 CPU seconds in a single R12000 processor. MAMMOTH-mult is particularly useful for large scale applications. AVAILABILITY: http://ub.cbm.uam.es/mammoth/mult.

An analysis of core deformations in protein superfamilies.

An analysis is presented on how structural cores modify their shape across homologous proteins, and whether or not a relationship exists between these structural changes and the vibrational normal modes that proteins experience as a result of the topological constraints imposed by the fold. A set of 35 representative, well-populated protein families is studied. The evolutionary directions of deformation are obtained by using multiple structural alignments to superimpose the structures and extract a conserved core, together with principal components analysis to extract the main deformation modes from the three-dimensional superimposition. In parallel, a low-resolution normal mode analysis technique is employed to study the properties of the mechanical core plasticity of these same families. We show that the evolutionary deformations span a low dimensional space of 4-5 dimensions on average. A statistically significant correspondence exists between these principal deformations and the approximately 20 slowest vibrational modes accessible to a particular topology. We conclude that, to a significant extent, the structural response of a protein topology to sequence changes takes place by means of collective deformations along combinations of a small number of low-frequency modes. The findings have implications in structure prediction by homology modeling.