Improving the Accuracy of Protein Thermostability Predictions for Single Point Mutations

Robust prediction of relative binding energies for protein-protein complex mutations using free energy perturbation calculations

Computational free energy-based methods have the potential to significantly improve throughput and decrease costs of protein design efforts. Such methods must reach a high level of reliability, accuracy, and automation to be effectively deployed in practical industrial settings in a way that impacts protein design projects. Here, we present a benchmark study for the calculation of relative changes in protein-protein binding affinity for single point mutations across a variety of systems from the literature, using free energy perturbation (FEP+) calculations. We describe a method for robust treatment of alternate protonation states for titratable amino acids, which yields improved correlation with and reduced error compared to experimental binding free energies. Following careful analysis of the largest outlier cases in our dataset, we assess limitations of the default FEP+ protocols and introduce an automated script which identifies probable outlier cases that may require additional scrutiny and calculates an empirical correction for a subset of charge-related outliers. Through a series of three additional case study systems, we discuss how protein FEP+ can be applied to real-world protein design projects, and suggest areas of further study.

Potts Hamiltonian Models and Molecular Dynamics Free Energy Simulations for Predicting the Impact of Mutations on Protein Kinase Stability

Single-point mutations in kinase proteins can affect their stability and fitness, and computational analysis of these effects can provide insights into the relationships among protein sequence, structure, and function for this enzyme family. To assess the impact of mutations on protein stability, we used a sequence-based Potts Hamiltonian model trained on a kinase family multiple-sequence alignment (MSA) to calculate the statistical energy (fitness) effects of mutations and compared these against relative folding free energies (ΔΔGs) calculated from all-atom molecular dynamics free energy perturbation (FEP) simulations in explicit solvent. The fitness effects of mutations in the Potts model (ΔEs) showed good agreement with experimental thermostability data (Pearson r = 0.68), similar to the correlation we observed with ΔΔGs predicted from structure-based relative FEP simulations. Recognizing the possible advantages of using Potts models to rapidly estimate protein stability effects of kinase mutations seen in cancer genomics data, we used the Potts statistical energy model to estimate the stability effects of 65 conservative and nonconservative mutations across three distinct kinases (Wee1, Abl1, and Cdc7) with somatic mutations reported in the Genomic Data Commons (GDC) database. The ΔEs of these mutations calculated from the Potts model are consistent with the corresponding ΔΔGs from FEP simulations (Pearson ratio of 0.72). The agreement between these methods suggests that the Potts model may be used as a sequence-based tool for high-throughput screening of mutational effects as part of a computational pipeline for predicting the stability effects of mutations. We also demonstrate how the scalability of the fitness-based Potts model calculations permits analyses that are not easily accessed using FEP simulations. To this end, we employed site-saturation mutagenesis in the Potts model in order to investigate the relative stability effects of mutations seen in different cancer evolutionary scenarios. We used this approach to analyze the effects of drug pressure in Abl kinase by contrasting the relative fitness penalties of somatic mutations seen in miscellaneous cancer types with those calculated for mutations associated with cancer drug resistance. We observed that, in contrast to somatic mutations of Abl seen in various tumors that appear to have evolved neutrally, cancer mutations that evolved under drug pressure in Abl-targeted therapies tend to preserve enzyme stability.

Selection pressures on evolution of ribonuclease H explored with rigorous free-energy-based design

Understanding natural protein evolution and designing novel proteins are motivating interest in development of high-throughput methods to explore large sequence spaces. In this work, we demonstrate the application of multisite λ dynamics (MSλD), a rigorous free energy simulation method, and chemical denaturation experiments to quantify evolutionary selection pressure from sequence-stability relationships and to address questions of design. This study examines a mesophilic phylogenetic clade of ribonuclease H (RNase H), furthering its extensive characterization in earlier studies, focusing on E. coli RNase H (ecRNH) and a more stable consensus sequence (AncCcons) differing at 15 positions. The stabilities of 32,768 chimeras between these two sequences were computed using the MSλD framework. The most stable and least stable chimeras were predicted and tested along with several other sequences, revealing a designed chimera with approximately the same stability increase as AncCcons, but requiring only half the mutations. Comparing the computed stabilities with experiment for 12 sequences reveals a Pearson correlation of 0.86 and root mean squared error of 1.18 kcal/mol, an unprecedented level of accuracy well beyond less rigorous computational design methods. We then quantified selection pressure using a simple evolutionary model in which sequences are selected according to the Boltzmann factor of their stability. Selection temperatures from 110 to 168 K are estimated in three ways by comparing experimental and computational results to evolutionary models. These estimates indicate selection pressure is high, which has implications for evolutionary dynamics and for the accuracy required for design, and suggests accurate high-throughput computational methods like MSλD may enable more effective protein design.

Structural-guided design to improve the catalytic performance of aldo-keto reductase KdAKR

Aldo-keto reductases (AKRs) are important biocatalysts that can be used to synthesize chiral pharmaceutical alcohols. In this study, the catalytic activity and stereoselectivity of a NADPH-dependent AKR from Kluyveromyces dobzhanskii (KdAKR) toward t-butyl 6-chloro (5S)-hydroxy-3-oxohexanoate ((5S)-CHOH) were improved by mutating its residues in the loop regions around the substrate-binding pocket. And the thermostability of KdAKR was improved by a consensus sequence method targeted on the flexible regions. The best mutant M6 (Y28A/L58I/I63L/G223P/Y296W/W297H) exhibited a 67-fold higher catalytic efficiency compared to the wild-type (WT) KdAKR, and improved R-selectivity toward (5S)-CHOH (dep value from 47.6% to >99.5%). Moreover, M6 exhibited a 6.3-fold increase in half-life (t1/2 ) at 40°C compared to WT. Under the optimal conditions, M6 completely converted 200 g/L (5S)-CHOH to diastereomeric pure t-butyl 6-chloro-(3R, 5S)-dihydroxyhexanoate ((3R, 5S)-CDHH) within 8.0 h, with a space-time yield of 300.7 g/L/day. Our results deepen the understandings of the structure-function relationship of AKRs, providing a certain guidance for the modification of other AKRs.

Comparing Supervised Learning and Rigorous Approach for Predicting Protein Stability upon Point Mutations in Difficult Targets

Accurate prediction of protein stability upon a point mutation has important applications in drug discovery and personalized medicine. It remains a challenging issue in computational biology. Existing computational prediction methods, which range from mechanistic to supervised learning approaches, have experienced limited progress over the last few decades. This stagnation is largely due to their heavy reliance on both the quantity and quality of the training data. This is evident in recent state-of-the-art methods that continue to yield substantial errors on two challenging blind test sets: frataxin and p53, with average root-mean-square errors exceeding 3 and 1.5 kcal/mol, respectively, which is still above the theoretical 1 kcal/mol prediction barrier. Rigorous approaches, on the other hand, offer greater potential for accuracy without relying on training data but are computationally demanding and require both wild-type and mutant structure information. Although they showed high accuracy for conserving mutations, their performance is still limited for charge-changing mutation cases. This might be due to the lack of an available mutant structure, often represented by a simplified capped peptide. The recent advances in protein structure prediction methods now make it possible to obtain structures comparable to experimental ones, including complete mutant structure information. In this work, we compare the performance of supervised learning-based methods and rigorous approaches for predicting protein stability on point mutations in difficult targets: frataxin and p53. The rigorous alchemical method significantly surpasses state-of-the-art techniques in terms of both the root-mean-squared error and Pearson correlation coefficient in these two challenging blind test sets. Additionally, we propose an improved alchemical method that employs the pmx double-system/single-box approach to accurately predict the folding free energy change upon both conserving and charge-changing mutations. The enhanced protocol can accurately predict both types of mutations, thereby outperforming existing state-of-the-art methods in overall performance.

Categorization of hotspots into three types - weak, moderate and strong to distinguish protein-protein versus protein-peptide interactions

Protein-protein and protein-peptide interactions (PPI and PPepI) belong to a similar category of interactions, yet seemingly subtle differences exist among them. To characterize differences between protein-protein (PP) and protein-peptide (PPep) interactions, we have focussed on two important classes of residues-hotspot and anchor residues. Using implicit solvation-based free energy calculations, a very large-scale alanine scanning has been performed on benchmark datasets, consisting of over 5700 interface residues. The differences in the two categories are more pronounced, if the data were divided into three distinct types, namely - weak hotspots (having binding free energy loss upon Ala mutation, ΔΔG, ∼2-10 kcal/mol), moderate hotspots (ΔΔG, ∼10-20 kcal/mol) and strong hotspots (ΔΔG ≥ ∼20 kcal/mol). The analysis suggests that for PPI, weak hotspots are predominantly populated by polar and hydrophobic residues. The distribution shifts towards charged and polar residues for moderate hotspot and charged residues (principally Arg) are overwhelmingly present in the strong hotspot. On the other hand, in the PPepI dataset, the distribution shifts from predominantly hydrophobic and polar (in the weak type) to almost similar preference for polar, hydrophobic and charged residues (in moderate type) and finally the charged residue (Arg) and Trp are mostly occupied in the strong type. The preferred anchor residues in both categories are Arg, Tyr and Leu, possessing bulky side chain and which also strike a delicate balance between side chain flexibility and rigidity. The present knowledge should aid in effective design of biologics, by augmentation or disruption of PPIs with peptides or peptidomimetics.Communicated by Ramaswamy H. Sarma.

Assessing variants of uncertain significance implicated in hearing loss using a comprehensive deafness proteome

Hearing loss is the leading sensory deficit, affecting ~ 5% of the population. It exhibits remarkable heterogeneity across 223 genes with 6328 pathogenic missense variants, making deafness-specific expertise a prerequisite for ascribing phenotypic consequences to genetic variants. Deafness-implicated variants are curated in the Deafness Variation Database (DVD) after classification by a genetic hearing loss expert panel and thorough informatics pipeline. However, seventy percent of the 128,167 missense variants in the DVD are "variants of uncertain significance" (VUS) due to insufficient evidence for classification. Here, we use the deep learning protein prediction algorithm, AlphaFold2, to curate structures for all DVD genes. We refine these structures with global optimization and the AMOEBA force field and use DDGun3D to predict folding free energy differences (∆∆GFold) for all DVD missense variants. We find that 5772 VUSs have a large, destabilizing ∆∆GFold that is consistent with pathogenic variants. When also filtered for CADD scores (> 25.7), we determine 3456 VUSs are likely pathogenic at a probability of 99.0%. Of the 224 genes in the DVD, 166 genes (74%) exhibit one or more missense variants predicted to cause a pathogenic change in protein folding stability. The VUSs prioritized here affect 119 patients (~ 3% of cases) sequenced by the OtoSCOPE targeted panel. Approximately half of these patients previously received an inconclusive report, and reclassification of these VUSs as pathogenic provides a new genetic diagnosis for six patients.

Mutation bias and the predictability of evolution

Predicting evolutionary outcomes is an important research goal in a diversity of contexts. The focus of evolutionary forecasting is usually on adaptive processes, and efforts to improve prediction typically focus on selection. However, adaptive processes often rely on new mutations, which can be strongly influenced by predictable biases in mutation. Here, we provide an overview of existing theory and evidence for such mutation-biased adaptation and consider the implications of these results for the problem of prediction, in regard to topics such as the evolution of infectious diseases, resistance to biochemical agents, as well as cancer and other kinds of somatic evolution. We argue that empirical knowledge of mutational biases is likely to improve in the near future, and that this knowledge is readily applicable to the challenges of short-term prediction. This article is part of the theme issue 'Interdisciplinary approaches to predicting evolutionary biology'.

Assessing Variants of Uncertain Significance Implicated in Hearing Loss Using a Comprehensive Deafness Proteome

Hearing loss is the leading sensory deficit, affecting ~ 5% of the population. It exhibits remarkable heterogeneity across 223 genes with 6,328 pathogenic missense variants, making deafness-specific expertise a prerequisite for ascribing phenotypic consequences to genetic variants. Deafness-implicated variants are curated in the Deafness Variation Database (DVD) after classification by a genetic hearing loss expert panel and thorough informatics pipeline. However, seventy percent of the 128,167 missense variants in the DVD are "variants of uncertain significance" (VUS) due to insufficient evidence for classification. Here, we use the deep learning protein prediction algorithm, AlphaFold2, to curate structures for all DVD genes. We refine these structures with global optimization and the AMOEBA force field and use DDGun3D to predict folding free energy differences (∆∆G Fold ) for all DVD missense variants. We find that 5,772 VUSs have a large, destabilizing ∆∆G Fold that is consistent with pathogenic variants. When also filtered for CADD scores (> 25.7), we determine 3,456 VUSs are likely pathogenic at a probability of 99.0%. These VUSs affect 119 patients (~ 3% of cases) sequenced by the OtoSCOPE targeted panel. Approximately half of these patients previously received an inconclusive report, and reclassification of these VUSs as pathogenic provides a new genetic diagnosis for six patients.

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.

Thermodynamic analysis of an entropically driven, high-affinity nanobody-HIV p24 interaction

Protein-protein interactions are fundamental to life processes. Complementary computational, structural, and biophysical studies of these interactions enable the forces behind their specificity and strength to be understood. Antibody fragments such as single-chain antibodies have the specificity and affinity of full antibodies but a fraction of their size, expediting whole molecule studies and distal effects without exceeding the computational capacity of modeling systems. We previously reported the crystal structure of a high-affinity nanobody 59H10 bound to HIV-1 capsid protein p24 and deduced key interactions using all-atom molecular dynamics simulations. We studied the properties of closely related medium (37E7) and low (48G11) affinity nanobodies, to understand how changes of three (37E7) or one (48G11) amino acids impacted these interactions; however, the contributions of enthalpy and entropy were not quantified. Here, we report the use of qualitative and quantitative experimental and in silico approaches to separate the contributions of enthalpy and entropy. We used complementary circular dichroism spectroscopy and molecular dynamics simulations to qualitatively delineate changes between nanobodies in isolation and complexed with p24. Using quantitative techniques such as isothermal titration calorimetry alongside WaterMap and Free Energy Perturbation protocols, we found the difference between high (59H10) and medium (37E7) affinity nanobodies on binding to HIV-1 p24 is entropically driven, accounted for by the release of unstable waters from the hydrophobic surface of 59H10. Our results provide an exemplar of the utility of parallel in vitro and in silico studies and highlight that differences in entropic interactions between amino acids and water molecules are sufficient to drive orders of magnitude differences in affinity.

Conformational Stability and Denaturation Processes of Proteins Investigated by Electrophoresis under Extreme Conditions

The functional structure of proteins results from marginally stable folded conformations. Reversible unfolding, irreversible denaturation, and deterioration can be caused by chemical and physical agents due to changes in the physicochemical conditions of pH, ionic strength, temperature, pressure, and electric field or due to the presence of a cosolvent that perturbs the delicate balance between stabilizing and destabilizing interactions and eventually induces chemical modifications. For most proteins, denaturation is a complex process involving transient intermediates in several reversible and eventually irreversible steps. Knowledge of protein stability and denaturation processes is mandatory for the development of enzymes as industrial catalysts, biopharmaceuticals, analytical and medical bioreagents, and safe industrial food. Electrophoresis techniques operating under extreme conditions are convenient tools for analyzing unfolding transitions, trapping transient intermediates, and gaining insight into the mechanisms of denaturation processes. Moreover, quantitative analysis of electrophoretic mobility transition curves allows the estimation of the conformational stability of proteins. These approaches include polyacrylamide gel electrophoresis and capillary zone electrophoresis under cold, heat, and hydrostatic pressure and in the presence of non-ionic denaturing agents or stabilizers such as polyols and heavy water. Lastly, after exposure to extremes of physical conditions, electrophoresis under standard conditions provides information on irreversible processes, slow conformational drifts, and slow renaturation processes. The impressive developments of enzyme technology with multiple applications in fine chemistry, biopharmaceutics, and nanomedicine prompted us to revisit the potentialities of these electrophoretic approaches. This feature review is illustrated with published and unpublished results obtained by the authors on cholinesterases and paraoxonase, two physiologically and toxicologically important enzymes.

Assessment of a Computational Approach to Predict Drug Resistance Mutations for HIV, HBV and SARS-CoV-2

Viral resistance is a worldwide problem mitigating the effectiveness of antiviral drugs. Mutations in the drug-targeting proteins are the primary mechanism for the emergence of drug resistance. It is essential to identify the drug resistance mutations to elucidate the mechanism of resistance and to suggest promising treatment strategies to counter the drug resistance. However, experimental identification of drug resistance mutations is challenging, laborious and time-consuming. Hence, effective and time-saving computational structure-based approaches for predicting drug resistance mutations are essential and are of high interest in drug discovery research. However, these approaches are dependent on accurate estimation of binding free energies which indirectly correlate to the computational cost. Towards this goal, we developed a computational workflow to predict drug resistance mutations for any viral proteins where the structure is known. This approach can qualitatively predict the change in binding free energies due to mutations through residue scanning and Prime MM-GBSA calculations. To test the approach, we predicted resistance mutations in HIV-RT selected by (-)-FTC and demonstrated accurate identification of the clinical mutations. Furthermore, we predicted resistance mutations in HBV core protein for GLP-26 and in SARS-CoV-2 3CLpro for nirmatrelvir. Mutagenesis experiments were performed on two predicted resistance and three predicted sensitivity mutations in HBV core protein for GLP-26, corroborating the accuracy of the predictions.

Large-scale application of free energy perturbation calculations for antibody design

Alchemical free energy perturbation (FEP) is a rigorous and powerful technique to calculate the free energy difference between distinct chemical systems. Here we report our implementation of automated large-scale FEP calculations, using the Amber software package, to facilitate antibody design and evaluation. In combination with Hamiltonian replica exchange, our FEP simulations aim to predict the effect of mutations on both the binding affinity and the structural stability. Importantly, we incorporate multiple strategies to faithfully estimate the statistical uncertainties in the FEP results. As a case study, we apply our protocols to systematically evaluate variants of the m396 antibody for their conformational stability and their binding affinity to the spike proteins of SARS-CoV-1 and SARS-CoV-2. By properly adjusting relevant parameters, the particle collapse problems in the FEP simulations are avoided. Furthermore, large statistical errors in a small fraction of the FEP calculations are effectively reduced by extending the sampling, such that acceptable statistical uncertainties are achieved for the vast majority of the cases with a modest total computational cost. Finally, our predicted conformational stability for the m396 variants is qualitatively consistent with the experimentally measured melting temperatures. Our work thus demonstrates the applicability of FEP in computational antibody design.

Whole-genome sequencing reveals de-novo mutations associated with nonsyndromic cleft lip/palate

The majority (85%) of nonsyndromic cleft lip with or without cleft palate (nsCL/P) cases occur sporadically, suggesting a role for de novo mutations (DNMs) in the etiology of nsCL/P. To identify high impact protein-altering DNMs that contribute to the risk of nsCL/P, we conducted whole-genome sequencing (WGS) analyses in 130 African case-parent trios (affected probands and unaffected parents). We identified 162 high confidence protein-altering DNMs some of which are based on available evidence, contribute to the risk of nsCL/P. These include novel protein-truncating DNMs in the ACTL6A, ARHGAP10, MINK1, TMEM5 and TTN genes; as well as missense variants in ACAN, DHRS3, DLX6, EPHB2, FKBP10, KMT2D, RECQL4, SEMA3C, SEMA4D, SHH, TP63, and TULP4. Many of these protein-altering DNMs were predicted to be pathogenic. Analysis using mouse transcriptomics data showed that some of these genes are expressed during the development of primary and secondary palate. Gene-set enrichment analysis of the protein-altering DNMs identified palatal development and neural crest migration among the few processes that were significantly enriched. These processes are directly involved in the etiopathogenesis of clefting. The analysis of the coding sequence in the WGS data provides more evidence of the opportunity for novel findings in the African genome.

Thermostability engineering of industrial enzymes through structure modification

Thermostability is an essential requirement of enzymes in the industrial processes to catalyze the reactions at high temperatures; thus, enzyme engineering through directed evolution, semi-rational design and rational design are commonly employed to construct desired thermostable mutants. Several strategies are implemented to fulfill enzymes' thermostability demand including decreasing the entropy of the unfolded state through substitutions Gly → Xxx or Xxx → Pro, hydrogen bond, salt bridge, introducing two different simultaneous interactions through single mutant, hydrophobic interaction, filling the hydrophobic cavity core, decreasing surface hydrophobicity, truncating loop, aromatic-aromatic interaction and introducing positively charged residues to enzyme surface. In the current review, horizons about compatibility between secondary structures and substitutions at preferable structural positions to generate the most desirable thermostability in industrial enzymes are broadened. KEY POINTS: • Protein engineering is a powerful tool for generating thermostable industrial enzymes. • Directed evolution and rational design are practical approaches in enzyme engineering. • Substitutions in preferable structural positions can increase thermostability.

The Use of Molecular Dynamics Simulation Method to Quantitatively Evaluate the Affinity between HBV Antigen T Cell Epitope Peptides and HLA-A Molecules

Chronic hepatitis B virus (HBV), a potentially life-threatening liver disease, makes people vulnerable to serious diseases such as cancer. T lymphocytes play a crucial role in clearing HBV virus, while the pathway depends on the strong binding of T cell epitope peptide and HLA. However, the experimental identification of HLA-restricted HBV antigenic peptides is extremely time-consuming. In this study, we provide a novel prediction strategy based on structure to assess the affinity between the HBV antigenic peptide and HLA molecule. We used residue scanning, peptide docking and molecular dynamics methods to obtain the molecular docking model of HBV peptide and HLA, and then adopted the MM-GBSA method to calculate the binding affinity of the HBV peptide–HLA complex. Overall, we collected 59 structures of HLA-A from Protein Data Bank, and finally obtained 352 numerical affinity results to figure out the optimal bind choice between the HLA-A molecules and 45 HBV T cell epitope peptides. The results were highly consistent with the qualitative affinity level determined by the competitive peptide binding assay, which confirmed that our affinity prediction process based on an HLA structure is accurate and also proved that the homologous modeling strategy for HLA-A molecules in this study was reliable. Hence, our work highlights an effective way by which to predict and screen for HLA-peptide binding that would improve the treatment of HBV infection.

Computer-Based Immunoinformatic Analysis to Predict Candidate T-Cell Epitopes for SARS-CoV-2 Vaccine Design

Since the first outbreak of coronavirus disease 2019 (COVID-19), caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) in 2019, its high infectivity led to its prevalence around the world in an exceptionally short time. Efforts have been made to control the ongoing outbreak, and among them, vaccine developments are going on high priority. New clinical trials add to growing evidence that vaccines from many countries were highly effective at preventing SARS-CoV-2 virus infection. One of them is B cell-based vaccines, which were common during a pandemic. However, neutralizing antibody therapy becomes less effective when viruses mutate. In order to tackle the problem, we focused on T-cell immune mechanism. In this study, the mutated strains of the virus were selected globally from India (B.1.617.1 and B.1.617.2), United Kingdom (B.1.1.7), South Africa (B.1.351), and Brazil (P.1), and the overlapping peptides were collected based on mutation sites of S-protein. After that, residue scanning was used to predict the affinity between overlapping peptide and HLA-A*11:01, the most frequent human leukocyte antigen (HLA) allele among the Chinese population. Then, the binding free energy was evaluated with molecular docking to further verify the affinity changes after the mutations happen in the virus genomes. The affinity test results of three epitopes on spike protein from experimental validation were consistent with our predicted results, thereby supporting the inclusion of the epitope ₃₇₄FSTFKCYGL₃₈₂ in future vaccine design and providing a useful reference route to improve vaccine development.

Addressing Intersite Coupling Unlocks Large Combinatorial Chemical Spaces for Alchemical Free Energy Methods

Alchemical free energy methods are playing a growing role in molecular design, both for computer-aided drug design of small molecules and for computational protein design. Multisite λ dynamics (MSλD) is a uniquely scalable alchemical free energy method that enables more efficient exploration of combinatorial alchemical spaces encountered in molecular design, but simulations have typically been limited to a few hundred ligands or sequences. Here, we focus on coupling between sites to enable scaling to larger alchemical spaces. We first discuss updates to the biasing potentials that facilitate MSλD sampling to include coupling terms and show that this can provide more thorough sampling of alchemical states. We then harness coupling between sites by developing a new free energy estimator based on the Potts models underlying direct coupling analysis, a method for predicting contacts from sequence coevolution, and find it yields more accurate free energies than previous estimators. The sampling requirements of the Potts model estimator scale with the square of the number of sites, a substantial improvement over the exponential scaling of the standard estimator. This opens up exploration of much larger alchemical spaces with MSλD for molecular design.

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.

A highly accurate metadynamics-based Dissociation Free Energy method to calculate protein-protein and protein-ligand binding potencies

Although seeking to develop a general and accurate binding free energy calculation method for protein-protein and protein-ligand interactions has been a continuous effort for decades, only limited successes have been obtained so far. Here, we report the development of a metadynamics-based procedure that calculates Dissociation Free Energy (DFE) and its application to 19 non-congeneric protein-protein complexes and hundreds of protein-ligand complexes covering eight targets. We achieved very high correlations in comparison to experimental binding free energies for these diverse sets of systems, demonstrating the generality and accuracy of the method. Since structures of most proteins are available owing to the recent success of prediction by artificial intelligence, a general free energy method such as DFE, combined with other methods, can make structure-based drug design a widely viable and reliable solution to develop both traditional small molecule drugs and biologic drugs as well as PROTACS.

Computational design of a cutinase for plastic biodegradation by mining molecular dynamics simulations trajectories

Polyethylene terephthalate (PET) has caused serious environmental concerns but could be degraded at high temperature. Previous studies show that cutinase from Thermobifida fusca KW3 (TfCut2) is capable of degrading and upcycling PET but is limited by its thermal stability. Nowadays, Popular protein stability modification methods rely mostly on the crystal structures, but ignore the fact that the actual conformation of protein is complex and constantly changing. To solve these problems, we developed a computational approach to design variants with enhanced protein thermal stability by mining Molecular Dynamics simulation trajectories using Machine Learning methods (MDL). The optimal classification accuracy and the optimal Pearson correlation coefficient of MDL model were 0.780 and 0.716, respectively. And we successfully designed variants with high ΔT_m values using MDL method. The optimal variant S121P/D174S/D204P had the highest ΔT_m value of 9.3 °C, and the PET degradation ratio increased by 46.42-fold at 70℃, compared with that of wild type TfCut2. These results deepen our understanding on the complex conformations of proteins and may enhance the plastic recycling and sustainability at glass transition temperature.

The reorganization energy of compounds upon binding to proteins, from dynamic and solvated bound and unbound states

The intramolecular reorganization energy (ΔE_Reorg) of compounds upon binding to proteins is a component of the binding free energy, which has long received particular attention, for fundamental and practical reasons. Understanding ΔE_Reorg would benefit the science of molecular recognition and drug design. For instance, the tolerable strain energy of compounds upon binding has been elusive. Prior studies found some large ΔE_Reorg values (e.g. > 10 kcal/mol), received with skepticism since they imply excessive opposition to binding. Indeed, estimating ΔE_Reorg is technically difficult. Typically, ΔE_Reorg has been approached by taking two energy-minimized conformers representing the bound and unbound states, and subtracting their conformational energy. This is a drastic oversimplification, liable to conformational collapse of the unbound conformer. Instead, the present work applies extensive molecular dynamics (MD) and the modern OPLS3 force-field to simulate compounds bound and unbound states, in explicit solvent under physically relevant conditions. The thermalized unbound compounds populate multiple conformations, not reducible to one or a few energy-minimized conformers. The intramolecular energies in the bound and unbound states were averaged over pertinent conformational ensembles, and the reorganization enthalpy upon binding (ΔH_Reorg) deduced by subtraction. This was applied to 76 systems, including 43 approved drugs, carefully selected for i) the quality of the bioactive X-ray structures and ii) the diversity of the chemotypes, their properties and protein targets. It yielded comparatively low ΔH_Reorg values (median = 1.4 kcal/mol, mean = 3.0 kcal/mol). A new finding is the observation of negative ΔH_Reorg values. Indeed, reorganization energies do not have to oppose binding, e.g. when intramolecular interactions stabilize preferentially the bound state. Conversely, even with competing water molecules, intramolecular interactions can occur predominantly in the unbound compound, and be replaced by intermolecular counterparts upon protein binding. Such disruption of intramolecular interactions upon binding gives rise to occasional larger ΔH_Reorg values. Such counterintuitive larger ΔH_Reorg values may be rationalized as a redistribution of interactions upon binding, qualitatively compatible with binding.

Accurate Prediction of Protein Thermodynamic Stability Changes upon Residue Mutation using Free Energy Perturbation

This work describes the application of a physics-based computational approach to predict the relative thermodynamic stability of protein variants, and evaluates the quantitative accuracy of those predictions compared to experimental data obtained from a diverse set of protein systems assayed at variable pH conditions. Physical stability is a key determinant of the clinical and commercial success of biological therapeutics, vaccines, diagnostics, enzymes and other protein-based products. Although experimental techniques for measuring the impact of amino acid residue mutation on the stability of proteins exist, they tend to be time consuming and costly, hence the need for accurate prediction methods. In contrast to many of the commonly available computational methods for stability prediction, the Free Energy Perturbation approach applied in this paper explicitly accounts for solvent effects and samples conformational dynamics using a rigorous molecular dynamics simulation process. On the entire validation dataset, consisting of 328 single point mutations spread across 14 distinct protein structures, our results show good overall correlation with experiment with an R² of 0.65 and a low mean unsigned error of 0.95 kcal/mol. Application of the FEP approach in conjunction with experimental assessment techniques offers opportunities to lower the time and expense of product development and reduce the risk of costly late-stage failures.

BLaDE: A Basic Lambda Dynamics Engine for GPU-Accelerated Molecular Dynamics Free Energy Calculations

There is an accelerating interest in practical applications of alchemical free energy methods to problems in protein design, constant pH simulations, and especially computer-aided drug design. In the present paper, we describe a basic lambda dynamics engine (BLaDE) that enables alchemical free energy simulations, including multisite λ dynamics (MSλD) simulations, on graphical processor units (GPUs). We find that BLaDE is 5 to 8 times faster than the current GPU implementation of MSλD-based free energy calculations in CHARMM. We also demonstrate that BLaDE running standard molecular dynamics attains a performance competitive with and sometimes exceeding that of the highly optimized OpenMM GPU code. BLaDE is available as a standalone program and through an API in CHARMM.

The future of biomolecular simulation in the pharmaceutical industry: what we can learn from aerodynamics modelling and weather prediction. Part 1. understanding the physical and computational complexity of in silico drug design

The predictive power of simulation has become embedded in the infrastructure of modern economies. Computer-aided design is ubiquitous throughout industry. In aeronautical engineering, built infrastructure and materials manufacturing, simulations are routinely used to compute the performance of potential designs before construction. The ability to predict the behaviour of products is a driver of innovation by reducing the cost barrier to new designs, but also because radically novel ideas can be piloted with relatively little risk. Accurate weather forecasting is essential to guide domestic and military flight paths, and therefore the underpinning simulations are critical enough to have implications for national security. However, in the pharmaceutical and biotechnological industries, the application of computer simulations remains limited by the capabilities of the technology with respect to the complexity of molecular biology and human physiology. Over the last 30 years, molecular-modelling tools have gradually gained a degree of acceptance in the pharmaceutical industry. Drug discovery has begun to benefit from physics-based simulations. While such simulations have great potential for improved molecular design, much scepticism remains about their value. The motivations for such reservations in industry and areas where simulations show promise for efficiency gains in preclinical research are discussed. In this, the first of two complementary papers, the scientific and technical progress that needs to be made to improve the predictive power of biomolecular simulations, and how this might be achieved, is firstly discussed (Part 1). In Part 2, the status of computer simulations in pharma is contrasted with aerodynamics modelling and weather forecasting, and comments are made on the cultural changes needed for equivalent computational technologies to become integrated into life-science industries.

Surface salt bridges contribute to the extreme thermal stability of an FN3-like domain from a thermophilic bacterium

This study uses differential scanning calorimetry, X-ray crystallography, and molecular dynamics simulations to investigate the structural basis for the high thermal stability (melting temperature 97.5°C) of a FN3-like protein domain from thermophilic bacteria Thermoanaerobacter tengcongensis (FN3tt). FN3tt adopts a typical FN3 fold with a three-stranded beta sheet packing against a four-stranded beta sheet. We identified three solvent exposed arginine residues (R23, R25, and R72), which stabilize the protein through salt bridge interactions with glutamic acid residues on adjacent strands. Alanine mutation of the three arginine residues reduced melting temperature by up to 22°C. Crystal structures of the wild type (WT) and a thermally destabilized (∆Tm -19.7°C) triple mutant (R23L/R25T/R72I) were found to be nearly identical, suggesting that the destabilization is due to interactions of the arginine residues. Molecular dynamics simulations showed that the salt bridge interactions in the WT were stable and provided a dynamical explanation for the cooperativity observed between R23 and R25 based on calorimetry measurements. In addition, folding free energy changes computed using free energy perturbation molecular dynamics simulations showed high correlation with melting temperature changes. This work is another example of surface salt bridges contributing to the enhanced thermal stability of thermophilic proteins. The molecular dynamics simulation methods employed in this study may be broadly useful for in silico surface charge engineering of proteins.

Scaffold Hopping Transformations Using Auxiliary Restraints for Calculating Accurate Relative Binding Free Energies

In silico screening of drug-target interactions is a key part of the drug discovery process. Changes in the drug scaffold via contraction or expansion of rings, the breaking of rings, and the introduction of cyclic structures from acyclic structures are commonly applied by medicinal chemists to improve binding affinity and enhance favorable properties of candidate compounds. These processes, commonly referred to as scaffold hopping, are challenging to model computationally. Although relative binding free energy (RBFE) calculations have shown success in predicting binding affinity changes caused by perturbing R-groups attached to a common scaffold, applications of RBFE calculations to modeling scaffold hopping are relatively limited. Scaffold hopping inevitably involves breaking and forming bond interactions of quadratic functional forms, which is highly challenging. A novel method for handling ring opening/closure/contraction/expansion and linker contraction/expansion is presented here. To the best of our knowledge, RBFE calculations on linker contraction/expansion have not been previously reported. The method uses auxiliary restraints to hold the atoms at the ends of a bond in place during the breaking and forming of the bonds. The broad applicability of the method was demonstrated by examining perturbations involving small-molecule macrocycles and mutations of proline in proteins. High accuracy was obtained using the method for most of the perturbations studied. The rigor of the method was isolated from the force field by validating the method using relative and absolute hydration free energy calculations compared to standard simulation results. Unlike other methods that rely on λ-dependent functional forms for bond interactions, the method presented here can be employed using modern molecular dynamics software without modification of codes or force field functions.

A strategy for proline and glycine mutations to proteins with alchemical free energy calculations

Computation of the thermodynamic consequences of protein mutations holds great promise in protein biophysics and design. Alchemical free energy methods can give improved estimates of mutational free energies, and are already widely used in calculations of relative and absolute binding free energies in small molecule design problems. In principle, alchemical methods can address any amino acid mutation with an appropriate alchemical pathway, but identifying a strategy that produces such a path for proline and glycine mutations is an ongoing challenge. Most current strategies perturb only side chain atoms, while proline and glycine mutations also alter the backbone parameters and backbone ring topology. Some strategies also perturb backbone parameters and enable glycine mutations. This work presents a strategy that enables both proline and glycine mutations and comprises two key elements: a dual backbone with restraints and scaling of bonded terms, facilitating backbone parameter changes, and a soft bond in the proline ring, enabling ring topology changes in proline mutations. These elements also have utility for core hopping and macrocycle studies in computer-aided drug design. This new strategy shows slight improvements over an alternative side chain perturbation strategy for a set T4 lysozyme mutations lacking proline and glycine, and yields good agreement with experiment for a set of T4 lysozyme proline and glycine mutations not previously studied. To our knowledge this is the first report comparing alchemical predictions of proline mutations with experiment. With this strategy in hand, alchemical methods now have access to the full palette of amino acid mutations.

A Cloud Computing Platform for Scalable Relative and Absolute Binding Free Energy Predictions: New Opportunities and Challenges for Drug Discovery

Free energy perturbation (FEP) has become widely used in drug discovery programs for binding affinity prediction between candidate compounds and their biological targets. However, limitations of FEP applications also exist, including, but not limited to, high cost, long waiting time, limited scalability, and breadth of application scenarios. To overcome these problems, we have developed XFEP, a scalable cloud computing platform for both relative and absolute free energy predictions using optimized simulation protocols. XFEP enables large-scale FEP calculations in a more efficient, scalable, and affordable way, for example, the evaluation of 5000 compounds can be performed in 1 week using 50-100 GPUs with a computing cost roughly equivalent to the cost for the synthesis of only one new compound. By combining these capabilities with artificial intelligence techniques for goal-directed molecule generation and evaluation, new opportunities can be explored for FEP applications in the drug discovery stages of hit identification, hit-to-lead, and lead optimization based not only on structure exploitation within the given chemical series but also including evaluation and comparison of completely unrelated molecules during structure exploration in a larger chemical space. XFEP provides the basis for scalable FEP applications to become more widely used in drug discovery projects and to speed up the drug discovery process from hit identification to preclinical candidate compound nomination.

HARP: a database of structural impacts of systematic missense mutations in drug targets of Mycobacterium leprae

Computational Saturation Mutagenesis is an in-silico approach that employs systematic mutagenesis of each amino acid residue in the protein to all other amino acid types, and predicts changes in thermodynamic stability and affinity to the other subunits/protein counterparts, ligands and nucleic acid molecules. The data thus generated are useful in understanding the functional consequences of mutations in antimicrobial resistance phenotypes. In this study, we applied computational saturation mutagenesis to three important drug-targets in Mycobacterium leprae (M. leprae) for the drugs dapsone, rifampin and ofloxacin namely Dihydropteroate Synthase (DHPS), RNA Polymerase (RNAP) and DNA Gyrase (GYR), respectively. M. leprae causes leprosy and is an obligate intracellular bacillus with limited protein structural information associating mutations with phenotypic resistance outcomes in leprosy. Experimentally solved structures of DHPS, RNAP and GYR of M. leprae are not available in the Protein Data Bank, therefore, we modelled the structures of these proteins using template-based comparative modelling and introduced systematic mutations in each model generating 80,902 mutations and mutant structures for all the three proteins. Impacts of mutations on stability and protein-subunit, protein-ligand and protein-nucleic acid affinities were computed using various in-house developed and other published protein stability and affinity prediction software. A consensus impact was estimated for each mutation using qualitative scoring metrics for physicochemical properties and by a categorical grouping of stability and affinity predictions. We developed a web database named HARP (a database of Hansen's Disease Antimicrobial Resistance Profiles), which is accessible at the URL - https://harp-leprosy.org and provides the details to each of these predictions.