pmx Webserver: A User Friendly Interface for Alchemistry

Alchemical Free Energy Workflows for the Computation of Protein-Ligand Binding Affinities

Alchemical free energy methods can be used for the efficient computation of relative binding free energies during preclinical drug discovery stages. In recent years, this has been facilitated further by the implementation of workflows that enable non-experts to quickly and consistently set up the required simulations. Given the correct input structures, workflows handle the difficult aspects of setting up perturbations, including consistently defining the perturbable molecule, its atom mapping and topology generation, perturbation network generation, running of the simulations via different sampling methods, and analysis of the results. Different academic and commercial workflows are discussed, including FEW, FESetup, FEPrepare, CHARMM-GUI, Transformato, PMX, QLigFEP, TIES, ProFESSA, PyAutoFEP, BioSimSpace, FEP+, Flare, and Orion. These workflows differ in various aspects, such as mapping algorithms or enhanced sampling methods. Some workflows can accommodate more than one molecular dynamics (MD) engine and use external libraries for tasks. Differences between workflows can present advantages for different use cases, however a lack of interoperability of the workflows' components hinders systematic comparisons.

Titin domains with reduced core hydrophobicity cause dilated cardiomyopathy

The underlying genetic defect in most cases of dilated cardiomyopathy (DCM), a common inherited heart disease, remains unknown. Intriguingly, many patients carry single missense variants of uncertain pathogenicity targeting the giant protein titin, a fundamental sarcomere component. To explore the deleterious potential of these variants, we first solved the wild-type and mutant crystal structures of I21, the titin domain targeted by pathogenic variant p.C3575S. Although both structures are remarkably similar, the reduced hydrophobicity of deeply buried position 3575 strongly destabilizes the mutant domain, a scenario supported by molecular dynamics simulations and by biochemical assays that show no disulfide involving C3575. Prompted by these observations, we have found that thousands of similar hydrophobicity-reducing variants associate specifically with DCM. Hence, our results imply that titin domain destabilization causes DCM, a conceptual framework that not only informs pathogenicity assessment of gene variants but also points to therapeutic strategies counterbalancing protein destabilization.

Conformational Dynamics of the Most Efficient Carboxylase Contributes to Efficient CO2 Fixation

Crotonyl-CoA carboxylase/reductase (Ccr) is one of the fastest CO2 fixing enzymes and has become part of efficient artificial CO2-fixation pathways in vitro, paving the way for future applications. The underlying mechanism of its efficiency, however, is not yet completely understood. X-ray structures of different intermediates in the catalytic cycle reveal tetramers in a dimer of dimers configuration with two open and two closed active sites. Upon binding a substrate, this active site changes its conformation from the open state to the closed state. It is challenging to predict how these coupled conformational changes will alter the CO2 binding affinity to the reaction's active site. To determine whether the open or closed conformations of Ccr affect binding of CO2 to the active site, we performed all-atom molecular simulations of the various conformations of Ccr. The open conformation without a substrate showed the highest binding affinity. The CO2 binding sites are located near the catalytic relevant Asn81 and His365 residues and in an optimal position for CO2 fixation. Furthermore, they are unaffected by substrate binding, and CO2 molecules stay in these binding sites for a longer time. Longer times at these reactive binding sites facilitate CO2 fixation through the nucleophilic attack of the reactive enolate in the closed conformation. We previously demonstrated that the Asn81Leu variant cannot fix CO2. Simulations of the Asn81Leu variant explain the loss of activity through the removal of the Asn81 and His365 binding sites. Overall, our findings show that the conformational dynamics of the enzyme controls CO2 binding. Conformational changes in Ccr increase the level of CO2 in the open subunit before the substrate is bound, the active site closes, and the reaction starts. The full catalytic Ccr cycle alternates among CO2 addition, conformational change, and chemical reaction in the four subunits of the tetramer coordinated by communication between the two dimers.

Integration of Genomic Sequencing Drives Therapeutic Targeting of PDGFRA in T-Cell Acute Lymphoblastic Leukemia/Lymphoblastic Lymphoma

PURPOSE

Patients with relapsed or refractory T-cell acute lymphoblastic leukemia (T-ALL) or lymphoblastic lymphoma (T-LBL) have limited therapeutic options. Clinical use of genomic profiling provides an opportunity to identify targetable alterations to inform therapy.

EXPERIMENTAL DESIGN

We describe a cohort of 14 pediatric patients with relapsed or refractory T-ALL enrolled on the Leukemia Precision-based Therapy (LEAP) Consortium trial (NCT02670525) and a patient with T-LBL, discovering alterations in platelet-derived growth factor receptor-α (PDGFRA) in 3 of these patients. We identified a novel mutation in PDGFRA, p.D842N, and used an integrated structural modeling and molecular biology approach to characterize mutations at D842 to guide therapeutic targeting. We conducted a preclinical study of avapritinib in a mouse patient-derived xenograft (PDX) model of FIP1L1-PDGFRA and PDGFRA p.D842N leukemia.

RESULTS

Two patients with T-ALL in the LEAP cohort (14%) had targetable genomic alterations affecting PDGFRA, a FIP1-like 1 protein/PDGFRA (FIP1L1-PDGFRA) fusion and a novel mutation in PDGFRA, p.D842N. The D842N mutation resulted in PDGFRA activation and sensitivity to tested PDGFRA inhibitors. In a T-ALL PDX model, avapritinib treatment led to decreased leukemia burden, significantly prolonged survival, and even cured a subset of mice. Avapritinib treatment was well tolerated and yielded clinical benefit in a patient with refractory T-ALL.

CONCLUSIONS

Refractory T-ALL has not been fully characterized. Alterations in PDGFRA or other targetable kinases may inform therapy for patients with refractory T-ALL who otherwise have limited treatment options. Clinical genomic profiling, in real time, is needed for fully informed therapeutic decision making.

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.

Molecular Insights into the Covalent Binding of Zoxamide to the β-Tubulin of Botrytis cinerea

Botrytis cinerea is a fungal plant pathogen that causes significant economic losses in the agricultural industry worldwide. Fungicides that target microtubules, such as carbendazim (CBZ), diethofencarb (DEF), and zoxamide (ZOX), are widely used in crop protection against this pathogen. These groups of compounds exert their fungicidal activity by disrupting the microtubule assembly by binding to the β-tubulin subunit, provoking cell-cycle arrest and cell death. However, with the appearance of isolates resistant to these compounds, it is necessary to search for new alternatives to control this pathogenic fungus. In this work, we gained insight into the binding and stability of these fungicides in the benzimidazole binding site of B. cinerea β-tubulin through different computational approaches. Our molecular dynamics simulation replicas showed that R enantiomers of ZOX and its analog RH-4032 had better interaction profiles at the site compared to S enantiomers. The simulations also revealed that while the R-isomer fungicides formed H-bonds with the main chain carbonyl of V236 or the side chain residue of S314, only CBZ interacted with E198. Previous experimental data have identified key mutations in B. cinerea's β-tubulin gene that lead to the development of resistance or, on the contrary, increased sensitivity for treatment with these fungicide compounds. In agreement with experimental findings, alchemical free energy calculations showed that E198A and E198V mutations in B. cinerea β-tubulin have high sensitivity to (R)-ZOX, whereas the E198K mutation decreased its affinity. Similarly, the results obtained explain the resistance to CBZ of B. cinerea isolates with E198A/V/K mutations and the insensitivity of the wild-type organism to DEF. Our work provides a deeper insight into the molecular mechanism of action of these fungicides, highlighting the importance of understanding the interaction profiles to develop more effective antifungal agents.

Web tools support predicting protein-nucleic acid complexes stability with affinity changes

Numerous biological processes, such as transcription, replication, and translation, rely on protein-nucleic acid interactions (PNIs). Demonstrating the binding stability of protein-nucleic acid complexes is vital to deciphering the code for PNIs. Numerous web-based tools have been developed to attach importance to protein-nucleic acid stability, facilitating the prediction of PNIs characteristics rapidly. However, the data and tools are dispersed and lack comprehensive integration to understand the stability of PNIs better. In this review, we first summarize existing databases for evaluating the stability of protein-nucleic acid binding. Then, we compare and evaluate the pros and cons of web tools for forecasting the interaction energies of protein-nucleic acid complexes. Finally, we discuss the application of combining models and capabilities of PNIs. We may hope these web-based tools will facilitate the discovery of recognition mechanisms for protein-nucleic acid binding stability. This article is categorized under: RNA Interactions with Proteins and Other Molecules > Protein-RNA Recognition RNA Interactions with Proteins and Other Molecules > RNA-Protein Complexes RNA Interactions with Proteins and Other Molecules > Protein-RNA Interactions: Functional Implications.

Imidazole-4-N-acetamide Derivatives as a Novel Scaffold for Selective Targeting of Cyclin Dependent Kinases

The rational design of cyclin-dependent protein kinase (CDK) inhibitors presumes the development of approaches for accurate prediction of selectivity and the activity of small molecular weight anticancer drug candidates. Aiming at attenuation of general toxicity of low selectivity compounds, we herein explored the new chemotype of imidazole-4-N-acetamide substituted derivatives of the pan-CDK inhibitor PHA-793887. Newly synthesized compounds 1-4 containing an aliphatic methyl group or aromatic radicals at the periphery of the scaffold were analyzed for the prediction of relative free energies of binding to CDK1, -2, -5, and -9 using a protocol based on non-equilibrium (NEQ) thermodynamics. This methodology allows for the demonstration of a good correlation between the calculated parameters of interaction of 1-4 with individual targets and the values of inhibitory potencies in in vitro kinase assays. We provide evidence in support of NEQ thermodynamics as a time sparing, precise, and productive approach for generating chemical inhibitors of clinically relevant anticancer targets.

Toward the design and development of peptidomimetic inhibitors of the Ataxin-1 aggregation pathway

Spinocerebellar ataxia type 1 is a degenerative disorder caused by polyglutamine expansions and aggregation of Ataxin-1. The interaction between Capicua (CIC) and the AXH domain of Ataxin-1 protein has been suggested as a possible driver of aggregation for the expanded Ataxin-1 protein and the subsequent onset of spinocerebellar ataxia 1. Experimental studies have demonstrated that short constructs of CIC may prevent such aggregation and suggested this as a possible candidate to inspire the rational design of peptidomimetics. In this work, molecular modeling techniques, namely the alchemical mutation and force field-based molecular dynamics, have been employed to propose a pipeline for the rational design of a CIC-inspired inhibitor of the ataxin-1 aggregation pathway. In particular, this study has shown that the alchemical mutation can estimate the affinity between AXH and CIC with good correlation with experimental data, while molecular dynamics shed light on molecular mechanisms that occur for stabilization of the interaction between the CIC-inspired construct and the AXH domain of Ataxin-1. This work lays the foundation for a rational methodology for the in silico screening and design of peptidomimetics, which can expedite and streamline experimental studies to identify strategies for inhibiting the ataxin-1 aggregation pathway.

High-Throughput Molecular Dynamics-Based Alchemical Free Energy Calculations for Predicting the Binding Free Energy Change Associated with the Selected Omicron Mutations in the Spike Receptor-Binding Domain of SARS-CoV-2

The ongoing pandemic caused by SARS-CoV-2 has gone through various phases. Since the initial outbreak, the virus has mutated several times, with some lineages showing even stronger infectivity and faster spread than the original virus. Among all the variants, omicron is currently classified as a variant of concern (VOC) by the World Health Organization, as the previously circulating variants have been replaced by it. In this work, we have focused on the mutations observed in omicron sub lineages BA.1, BA.2, BA.4 and BA.5, particularly at the receptor-binding domain (RBD) of the spike protein that is responsible for the interactions with the host ACE2 receptor and binding of antibodies. Studying such mutations is particularly important for understanding the viral infectivity, spread of the disease and for tracking the escape routes of this virus from antibodies. Molecular dynamics (MD) based alchemical free energy calculations have been shown to be very accurate in predicting the free energy change, due to a mutation that could have a deleterious or a stabilizing effect on either the protein itself or its binding affinity to another protein. Here, we investigated the significance of five spike RBD mutations on the stability of the spike protein binding to ACE2 by free energy calculations using high throughput MD simulations. For comparison, we also used conventional MD simulations combined with a Molecular Mechanics-Generalized Born Surface Area (MM-GBSA) based approach, and compared our results with the available experimental data. Overall, the alchemical free energy calculations performed far better than the MM-GBSA approach in predicting the individual impact of the mutations. When considering the experimental variation, the alchemical free energy method was able to produce a relatively accurate prediction for N501Y, the mutant that has previously been reported to increase the binding affinity to hACE2. On the other hand, the other individual mutations seem not to have a significant effect on the spike RBD binding affinity towards hACE2.

Allosteric Inhibition of PTP1B by a Nonpolar Terpenoid

Protein tyrosine phosphatases (PTPs) are promising drug targets for treating a wide range of diseases such as diabetes, cancer, and neurological disorders, but their conserved active sites have complicated the design of selective therapeutics. This study examines the allosteric inhibition of PTP1B by amorphadiene (AD), a terpenoid hydrocarbon that is an unusually selective inhibitor. Molecular dynamics (MD) simulations carried out in this study suggest that AD can stably sample multiple neighboring sites on the allosterically influential C-terminus of the catalytic domain. Binding to these sites requires a disordered α7 helix, which stabilizes the PTP1B-AD complex and may contribute to the selectivity of AD for PTP1B over TCPTP. Intriguingly, the binding mode of AD differs from that of the most well-studied allosteric inhibitor of PTP1B. Indeed, biophysical measurements and MD simulations indicate that the two molecules can bind simultaneously. Upon binding, both inhibitors destabilize the α7 helix by disrupting interactions at the α3-α7 interface and prevent the formation of hydrogen bonds that facilitate closure of the catalytically essential WPD loop. These findings indicate that AD is a promising scaffold for building allosteric inhibitors of PTP1B and illustrate, more broadly, how unfunctionalized terpenoids can engage in specific interactions with protein surfaces.

Molecular basis of Toxoplasma gondii oryzalin resistance from a novel α-tubulin binding site model

Oryzalin (ORY) is a dinitroaniline derivative that inhibits the microtubule polymerization in plants and parasitic protozoa by selectively binding to the α-tubulin subunit. This herbicidal agent exhibits good antiprotozoal activity against major human parasites, such as Toxoplasma gondii (toxoplasmosis), Leishmania mexicana (leishmaniasis), and Plasmodium falciparum (malaria). Previous chemical mutagenesis assays on T. gondii α-tubulin (TgAT) have identified key mutations that lead to ORY resistance. Herein, we employed alchemical free energy methods and molecular dynamics simulations to determine if the ORY resistance mutations either decrease the TgAT's affinity of the compound or increase the protein stability. Our results here suggest that L136F and V202F mutations significantly decrease the affinity of ORY to TgAT, while T239I and V252L mutations diminish TgAT's flexibility. On the other hand, protein stability predictors determined that R243S mutation reduces TgAT stability due to the loss of its salt bridge interaction with E27. Interestingly, molecular dynamics simulations confirm that the loss of this key interaction leads to ORY binding site closure. Our study provides a better insight into the TgAT-ORY interaction, further supporting our recently proposed ORY-binding site.

Oncogenic Mutations in the DNA-Binding Domain of FOXO1 that Disrupt Folding: Quantitative Insights from Experiments and Molecular Simulations

FOXO1, a member of the family of winged-helix motif Forkhead box (FOX) transcription factors, is the most abundantly expressed FOXO member in mature B cells. Sequencing of diffuse large B-cell lymphoma (DLBCL) tumors and cell lines identified specific mutations in the forkhead domain linked to loss of function. Differential scanning calorimetry and thermal shift assays were used to characterize how eight of these mutations affect the stability of the FOX domain. Mutations L183P and L183R were found to be particularly destabilizing. Electrophoresis mobility shift assays show these same mutations also disrupt FOXO1 binding to their canonical DNA sequences, suggesting that the loss of function is due to destabilization of the folded structure. Computational modeling of the effect of mutations on FOXO1 folding was performed using alchemical free energy perturbation (FEP), and a Markov model of the entire folding reaction was constructed from massively parallel molecular simulations, which predicts folding pathways involving the late folding of helix α3. Although FEP can qualitatively predict the destabilization from L183 mutations, we find that a simple hydrophobic transfer model, combined with estimates of unfolded-state solvent-accessible surface areas from molecular simulations, is able to more accurately predict changes in folding free energies due to mutations. These results suggest that the atomic detail provided by simulations is important for the accurate prediction of mutational effects on folding stability. Corresponding disease-associated mutations in other FOX family members support further experimental and computational studies of the folding mechanism of FOX domains.

Venetoclax Decreases the Expression of the Spike Protein through Amino Acids Q493 and S494 in SARS-CoV-2

The new coronavirus disease 2019 (COVID-19) caused by the severe acute respiratory syndrome coronavirus (SARS-CoV-2) has been reported and spread globally. There is an urgent need to take urgent measures to treat and prevent further infection of this virus. Here, we use virtual drug screening to establish pharmacophore groups and analyze the ACE2 binding site of the spike protein with the ZINC drug database and DrugBank database by molecular docking and molecular dynamics simulations. Screening results showed that Venetoclax, a treatment drug for chronic lymphocytic leukemia, has a potential ability to bind to the spike protein of SARS-CoV-2. In addition, our in vitro study found that Venetoclax degraded the expression of the spike protein of SARS-CoV-2 through amino acids Q493 and S494 and blocked the interaction with the ACE2 receptor. Our results suggest that Venetoclax is a candidate for clinical prevention and treatment and deserves further research.

Pre-exascale HPC approaches for molecular dynamics simulations. Covid-19 research: A use case

Exascale computing has been a dream for ages and is close to becoming a reality that will impact how molecular simulations are being performed, as well as the quantity and quality of the information derived for them. We review how the biomolecular simulations field is anticipating these new architectures, making emphasis on recent work from groups in the BioExcel Center of Excellence for High Performance Computing. We exemplified the power of these simulation strategies with the work done by the HPC simulation community to fight Covid-19 pandemics. This article is categorized under:Data Science > Computer Algorithms and ProgrammingData Science > Databases and Expert SystemsMolecular and Statistical Mechanics > Molecular Dynamics and Monte-Carlo Methods.

How Peptides Bind to PSD-95/Discs-Large/ZO-1 Domains

PSD-95/discs-large/ZO-1 (PDZ) domains form a large family of adaptor proteins that bind to the C-terminal tails of their binding partner proteins. Via extensive molecular dynamics simulations and alchemical free energy calculations, we characterized the binding modi of phosphorylated and unphosphorylated EQVSAV peptides and of a EQVEAV phosphate mimic to the hPTP1E PDZ2 and MAGI1 PDZ1 domains. The simulations reproduced the well-known binding characteristics such as tight coordination of the peptidic carboxyl tail and pronounced hydrogen bonding between the peptide backbone and the backbone atoms of a β-sheet in PDZ. Overall, coordination by hPTP1E PDZ2 appeared tighter than by MAGI1 PDZ1. Simulations of wild-type PDZ and arginine mutants suggest that contacts with Arg79/85 in hPTP1E/MAGI1 are more important for the EQVEAV peptide than for EQVSAV. Alchemical free energy calculations and PaCS-MD simulations could well reproduce the difference in binding free energy between unphosphorylated EQVSAV and EQVEAV peptides and the absolute binding free energy of EQVSAV. However, likely due to small force field inaccuracies, the simulations erroneously favored binding of the phosphorylated peptide instead of its unphosphorylated counterpart, which is in contrast to the experiment.

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 general theoretical framework to design base editors with reduced bystander effects

Base editors (BEs) hold great potential for medical applications of gene therapy. However, high precision base editing requires BEs that can discriminate between the target base and multiple bystander bases within a narrow active window (4 - 10 nucleotides). Here, to assist in the design of these optimized editors, we propose a discrete-state stochastic approach to build an analytical model that explicitly evaluates the probabilities of editing the target base and bystanders. Combined with all-atom molecular dynamic simulations, our model reproduces the experimental data of A3A-BE3 and its variants for targeting the "TC" motif and bystander editing. Analyzing this approach, we propose several general principles that can guide the design of BEs with a reduced bystander effect. These principles are then applied to design a series of point mutations at T218 position of A3G-BEs to further reduce its bystander editing. We verify experimentally that the new mutations provide different levels of stringency on reducing the bystander editing at different genomic loci, which is consistent with our theoretical model. Thus, our study provides a computational-aided platform to assist in the scientifically-based design of BEs with reduced bystander effects.

Expanded Ensemble Methods Can be Used to Accurately Predict Protein-Ligand Relative Binding Free Energies

Alchemical free energy methods have become indispensable in computational drug discovery for their ability to calculate highly accurate estimates of protein-ligand affinities. Expanded ensemble (EE) methods, which involve single simulations visiting all of the alchemical intermediates, have some key advantages for alchemical free energy calculation. However, there have been relatively few examples published in the literature of using expanded ensemble simulations for free energies of protein-ligand binding. In this paper, as a test of expanded ensemble methods, we compute relative binding free energies using the Open Force Field Initiative force field (codename "Parsley") for 24 pairs of Tyk2 inhibitors derived from a congeneric series of 16 compounds. The EE predictions agree well with the experimental values (root-mean-square error (RMSE) of 0.94 ± 0.13 kcal mol^-1 and mean unsigned error (MUE) of 0.75 ± 0.12 kcal mol^-1). We find that while increasing the number of alchemical intermediates can improve the phase space overlap, faster convergence can be obtained with fewer intermediates, as long as acceptance rates are sufficient. We also find that convergence can be improved using more aggressive updating of biases, and that estimates can be improved by performing multiple independent EE calculations. This work demonstrates that EE is a viable option for alchemical free energy calculation. We discuss the implications of these findings for rational drug design, as well as future directions for improvement.

Perturbation Free-Energy Toolkit: An Automated Alchemical Topology Builder

Free-energy calculations play an important role in the application of computational chemistry to a range of fields, including protein biochemistry, rational drug design, or materials science. Importantly, the free-energy difference is directly related to experimentally measurable quantities such as partition and adsorption coefficients, water activity, and binding affinities. Among several techniques aimed at predicting free-energy differences, perturbation approaches, involving the alchemical transformation of one molecule into another through intermediate states, stand out as rigorous methods based on statistical mechanics. However, despite the importance of free-energy calculations, the applicability of the perturbation approaches is still largely impeded by a number of challenges, including the definition of the perturbation path, i.e., alchemical changes leading to the transformation of one molecule to the other. To address this, an automatic perturbation topology builder based on a graph-matching algorithm is developed, which can identify the maximum common substructure (MCS) of two or multiple molecules and provide the perturbation topologies suitable for free-energy calculations using the GROMOS and the GROMACS simulation packages. Various MCS search options are presented leading to alternative definitions of the perturbation pathway. Moreover, perturbation topologies generated using the default multistate MCS search are used to calculate the changes in free energy between lysine and its two post-translational modifications, 3-methyllysine and acetyllysine. The pairwise free-energy calculations performed on this test system led to a cycle closure of 0.5 ± 0.3 and 0.2 ± 0.2 kJ mol-1, with GROMOS and GROMACS simulation packages, respectively. The same relative free energies between the three states are obtained by employing the enveloping distribution sampling (EDS) approach when compared to the pairwise perturbations. Importantly, this toolkit is made available online as an open-source Python package (https://github.com/drazen-petrov/SMArt).

How phosphorylation of peptides affects their interaction with 14-3-3η domains

Members of the 14-3-3 domain family have important functions as adapter domains. Via an amphipathic groove on their protein surface they typically bind to disordered C-terminals of other proteins. Importantly, binding partners of 14-3-3 domains usually contain a phosphorylated serine or threonine residue at their binding interface and possess one of three different sequence motifs. Binding of the respective unphosphorylated versions of the peptides is typically strongly disfavored. There is a wealth of structural and thermodynamic data available for the phosphorylated forms but not for the unphosphorylated forms as the binding affinities seem to be too weak to be measurable experimentally. Here, we characterized the mechanistic details that govern the preference for the binding of phosphorylated peptides to 14-3-3η domains by means of molecular dynamics (MD) simulations. We found that the phosphate group is ideally coordinated in the binding pocket whereas the respective unphosphorylated side-chain counterpart is not. Thus, the binding preference results from the tight coordination of the phosphorylated residue at the center of the binding interface. Furthermore, MD simulations of 14-3-3η dimers showed a preference for the simultaneous binding of two phosphorylated peptides in agreement with their experimentally observed cooperativity.

Controlling Protein Crystallization by Free Energy Guided Design of Interactions at Crystal Contacts

Protein crystallization can function as an effective method for protein purification or formulation. Such an application requires a comprehensive understanding of the intermolecular protein–protein interactions that drive and stabilize protein crystal formation to ensure a reproducible process. Using alcohol dehydrogenase from Lactobacillus brevis (LbADH) as a model system, we probed in our combined experimental and computational study the effect of residue substitutions at the protein crystal contacts on the crystallizability and the contact stability. Increased or decreased contact stability was calculated using molecular dynamics (MD) free energy simulations and showed excellent qualitative correlation with experimentally determined increased or decreased crystallizability. The MD simulations allowed us to trace back the changes to their physical origins at the atomic level. Engineered charge–charge interactions as well as engineered hydrophobic effects could be characterized and were found to improve crystallizability. For example, the simulations revealed a redesigning of a water mediated electrostatic interaction (“wet contact”) into a water depleted hydrophobic effect (“dry contact”) and the optimization of a weak hydrogen bonding contact towards a strong one. These findings explained the experimentally found improved crystallizability. Our study emphasizes that it is difficult to derive simple rules for engineering crystallizability but that free energy simulations could be a very useful tool for understanding the contribution of crystal contacts for stability and furthermore could help guide protein engineering strategies to enhance crystallization for technical purposes.

Implementing and Assessing an Alchemical Method for Calculating Protein-Protein Binding Free Energy

Protein-protein binding is fundamental to most biological processes. It is important to be able to use computation to accurately estimate the change in protein-protein binding free energy due to mutations in order to answer biological questions that would be experimentally challenging, laborious, or time-consuming. Although nonrigorous free-energy methods are faster, rigorous alchemical molecular dynamics-based methods are considerably more accurate and are becoming more feasible with the advancement of computer hardware and molecular simulation software. Even with sufficient computational resources, there are still major challenges to using alchemical free-energy methods for protein-protein complexes, such as generating hybrid structures and topologies, maintaining a neutral net charge of the system when there is a charge-changing mutation, and setting up the simulation. In the current study, we have used the pmx package to generate hybrid structures and topologies, and a double-system/single-box approach to maintain the net charge of the system. To test the approach, we predicted relative binding affinities for two protein-protein complexes using a nonequilibrium alchemical method based on the Crooks fluctuation theorem and compared the results with experimental values. The method correctly identified stabilizing from destabilizing mutations for a small protein-protein complex, and a larger, more challenging antibody complex. Strong correlations were obtained between predicted and experimental relative binding affinities for both protein-protein systems.

SAMPL7 Host-Guest Challenge Overview: assessing the reliability of polarizable and non-polarizable methods for binding free energy calculations

The SAMPL challenges focus on testing and driving progress of computational methods to help guide pharmaceutical drug discovery. However, assessment of methods for predicting binding affinities is often hampered by computational challenges such as conformational sampling, protonation state uncertainties, variation in test sets selected, and even lack of high quality experimental data. SAMPL blind challenges have thus frequently included a component focusing on host-guest binding, which removes some of these challenges while still focusing on molecular recognition. Here, we report on the results of the SAMPL7 blind prediction challenge for host-guest affinity prediction. In this study, we focused on three different host-guest categories-a familiar deep cavity cavitand series which has been featured in several prior challenges (where we examine binding of a series of guests to two hosts), a new series of cyclodextrin derivatives which are monofunctionalized around the rim to add amino acid-like functionality (where we examine binding of two guests to a series of hosts), and binding of a series of guests to a new acyclic TrimerTrip host which is related to previous cucurbituril hosts. Many predictions used methods based on molecular simulations, and overall success was mixed, though several methods stood out. As in SAMPL6, we find that one strategy for achieving reasonable accuracy here was to make empirical corrections to binding predictions based on previous data for host categories which have been studied well before, though this can be of limited value when new systems are included. Additionally, we found that alchemical free energy methods using the AMOEBA polarizable force field had considerable success for the two host categories in which they participated. The new TrimerTrip system was also found to introduce some sampling problems, because multiple conformations may be relevant to binding and interconvert only slowly. Overall, results in this challenge tentatively suggest that further investigation of polarizable force fields for these challenges may be warranted.

Methodological uncertainties in drug-receptor binding free energy predictions based on classical molecular dynamics

Computational approaches are becoming an essential tool in modern drug design and discovery, with fast compound triaging using a combination of machine learning and docking techniques followed by molecular dynamics binding free energies assessment using alchemical techniques. The traditional MD-based alchemical free energy perturbation (FEP) method faces severe sampling issues that may limits its reliability in automated workflows. Here we review the major sources of uncertainty in FEP protocols for drug discovery, showing how the sampling problem can be effectively tackled by switching to nonequilibrium alchemical techniques.

Evaluating inositol phospholipid interactions with inward rectifier potassium channels and characterising their role in disease

Membrane proteins are frequently modulated by specific protein-lipid interactions. The activation of human inward rectifying potassium (hKir) channels by phosphoinositides (PI) has been well characterised. Here, we apply a coarse-grained molecular dynamics free-energy perturbation (CG-FEP) protocol to capture the energetics of binding of PI lipids to hKir channels. By using either a single- or multi-step approach, we establish a consistent value for the binding of PIP2 to hKir channels, relative to the binding of the bulk phosphatidylcholine phospholipid. Furthermore, by perturbing amino acid side chains on hKir6.2, we show that the neonatal diabetes mutation E179K increases PIP2 affinity, while the congenital hyperinsulinism mutation K67N results in a reduced affinity. We show good agreement with electrophysiological data where E179K exhibits a reduction in neomycin sensitivity, implying that PIP2 binds more tightly E179K channels. This illustrates the application of CG-FEP to compare affinities between lipid species, and for annotating amino acid residues.

KRAS K104 modification affects the KRASG12D-GEF interaction and mediates cell growth and motility

Mutant RAS genes play an important role in regulating tumors through lysine residue 104 to impair GEF-induced nucleotide exchange, but the regulatory role of KRAS K104 modification on the KRAS^G12D mutant remains unclear. Therefore, we simulated the acetylation site on the KRAS^G12D three-dimensional protein structure, including KRAS^G12D, KRAS^G12D/K104A and KRAS^G12D/K104Q, and determined their trajectories and binding free energy with GEF. KRAS^G12D/K104Q induced structural changes in the α2- and α3-helices, promoted KRAS instability and hampered GEF binding (ΔΔG = 6.14 kJ/mol). We found decreased binding to the Raf1 RBD by KRAS^G12D/K104Q and reduced cell growth, invasion and migration. Based on whole-genome cDNA microarray analysis, KRAS^G12D/K104Q decreased expression of NPIPA2, DUSP1 and IL6 in lung and ovarian cancer cells. This study reports computational and experimental analyses of Lys104 of KRAS^G12D and GEF, and the findings provide a target for exploration for future treatment.

Effects of Mutating Trp42 Residue on γD-Crystallin Stability

Oligomerization and aggregation of γD-crystallins (HγDC) in the eye lens is one of the main causes of cataract development. To date, several congenital mutations related to this protein are known to promote the formation of aggregates. Previous studies have demonstrated that mutations in W42 residue of HγDC lead to the generation of partially unfolded intermediates that are more prone to aggregate. To understand the role of W42 in the stability of HγDC, we performed alchemical free-energy calculations and all-atom molecular dynamics simulations of different W42 mutant models. Our results suggest that substitution of W42 by small size and/or polar residues promotes HγDC denaturation due to the entry of water molecules into the hydrophobic core of the N-terminal domain. Similar behavior was observed in the C-terminal domain of HγDC when mutating the W130 residue located in a homologous position. Moreover, the exposure of the hydrophobic core residues could lead to the formation of aggregation-prone partially unfolded species. Overall, this study takes a step toward understanding the role of HγDC in cataract development.

Advances in the calculation of binding free energies

In recent years, calculations of binding affinities from molecular simulations seem to have matured significantly. While the number of applications of such methods in drug design and biotechnology increases, the number of truly new methodological developments decreases. This review provides an overview of the current status of the field as reflected in recent publications. The focus is on the challenges that remain when using endstate, alchemical and pathway methods. For endstate methods this is the calculation of entropic contributions. For alchemical methods there are unsolved problems associated with the solvation of the active site, sampling slow degrees of freedom and when modifying the net charge. For pathway methods achieving sufficient sampling remains challenging. New trends are also highlighted, including the use of pathway methods for the quantification of protein-protein interactions.

Telomere uncapping by common oxidative guanine lesions: Insights from atomistic models

Oxidative damage to DNA is widely known to contribute to aging and disease. This relationship has been extensively studied for telomeres - structures that cap chromosome ends - due to their role in cell proliferation and senescence, and exceptional susceptibility to oxidation. Indeed, the repetitive telomeric DNA sequence contains the 5'-GGG-3' motif that has the lowest ionization potential of all trinucleotides. Accordingly, experiments consistently show that telomeric oxidative lesions are more abundant and persistent than elsewhere in the genome. This led to a hypothesis that telomeres act as sensors of prolonged oxidative stress and prevent carcinogenesis, as disruption of telomeric integrity triggers senescence or apoptosis. Here, we use atomistic alchemical Molecular Dynamics simulations to perform a combinatorial assessment of changes in DNA binding affinity of telomeric proteins induced by oxidative guanine lesions. We rank lesions by their effect on telomere integrity, as well as telomeric proteins by their sensitivity to DNA oxidation. While the binding of most proteins is abolished by DNA oxidation, HOT1 emerges as a notable exception, suggesting its potential role in sensing of oxidative damage. Through statistical analysis and free energy decomposition, we also identify common trends in structural responses of protein-DNA complexes that contribute to decreased binding affinity.

Large scale relative protein ligand binding affinities using non-equilibrium alchemy

Ligand binding affinity calculations based on molecular dynamics (MD) simulations and non-physical (alchemical) thermodynamic cycles have shown great promise for structure-based drug design. However, their broad uptake and impact is held back by the notoriously complex setup of the calculations. Only a few tools other than the free energy perturbation approach by Schrödinger Inc. (referred to as FEP+) currently enable end-to-end application. Here, we present for the first time an approach based on the open-source software pmx that allows to easily set up and run alchemical calculations for diverse sets of small molecules using the GROMACS MD engine. The method relies on theoretically rigorous non-equilibrium thermodynamic integration (TI) foundations, and its flexibility allows calculations with multiple force fields. In this study, results from the Amber and Charmm force fields were combined to yield a consensus outcome performing on par with the commercial FEP+ approach. A large dataset of 482 perturbations from 13 different protein-ligand datasets led to an average unsigned error (AUE) of 3.64 ± 0.14 kJ mol-1, equivalent to Schrödinger's FEP+ AUE of 3.66 ± 0.14 kJ mol-1. For the first time, a setup is presented for overall high precision and high accuracy relative protein-ligand alchemical free energy calculations based on open-source software.

Thermodynamic Stability of Human γD-Crystallin Mutants Using Alchemical Free-Energy Calculations

γD-Crystallin (HγDC) is a key structural protein in the human lens, whose aggregation has been associated with the development of cataracts. Single-point mutations and post-translational modifications destabilize HγDC interactions, forming partially folded intermediates, where hydrophobic residues are exposed and thus triggering its aggregation. In this work, we used alchemical free-energy calculations to predict changes in thermodynamic stability (ΔΔG) of 10 alanine-scanning variants and 12 HγDC mutations associated with the development of congenital cataract. Our results show that W42R is the most destabilizing mutation in HγDC. This has been corroborated through experimental determination of ΔΔG employing differential scanning calorimetry. Calculations of hydration free energies from the HγDC WT and the W42R mutant suggested that the mutant has a higher aggregation propensity. Our combined theoretical and experimental results contribute to understand HγDC destabilization and aggregation mechanisms in age-onset cataracts.

Accurate Calculation of Free Energy Changes upon Amino Acid Mutation

Molecular dynamics based free energy calculations allow for a robust and accurate evaluation of free energy changes upon amino acid mutation in proteins. In this chapter we cover the basic theoretical concepts important for the use of calculations utilizing the non-equilibrium alchemical switching methodology. We further provide a detailed step-by-step protocol for estimating the effect of a single amino acid mutation on protein thermostability. In addition, the potential caveats and solutions to some frequently encountered issues concerning the non-equilibrium alchemical free energy calculations are discussed. The protocol comprises details for the hybrid structure/topology generation required for alchemical transitions, equilibrium simulation setup, and description of the fast non-equilibrium switching. Subsequently, the analysis of the obtained results is described. The steps in the protocol are complemented with an illustrative practical application: a destabilizing mutation in the Trp cage mini protein. The concepts that are described are generally applicable. The shown example makes use of the pmx software package for the free energy calculations using Gromacs as a molecular dynamics engine. Finally, we discuss how the current protocol can readily be adapted to carry out charge-changing or multiple mutations at once, as well as large-scale mutational scans.

Making sense of the past: hyperstability of ancestral thioredoxins explained by free energy simulations

Thioredoxin (Trx), a small and globular protein, is present in all kinds of organisms, from Archea to higher mammals. Throughout evolution, the Trx sequence has undergone subtle modifications to adapt to varying environmental conditions. The high degree of sequence conservation makes Trx very amenable to ancestral protein reconstruction techniques. In this work, we address the study of the structural and energetic determinants of thermostability in E. coli Trx using a dataset of mutations inspired by ancestral reconstruction. We compute, from first principles, the expected contribution of 19 different amino acid substitutions to the stability (ΔΔG) and the melting temperature (ΔT_m) of the protein. We also describe the specific changes in structure and protein dynamics responsible for the stabilizing or destabilizing effects of these mutations. Our results point to local and independent changes for most of the variants. Our predictions are accurate enough to substantiate the proposal of new hypotheses regarding evolutionary relationships between mutations, as in the case of T89R, P68A and G74S or K90L and F102A, and reach beyond the initial set to suggest improved variants, such as K90I or K90Y.

Computational Analysis of the Molecular Mechanism of RamR Mutations Contributing to Antimicrobial Resistance in Salmonella enterica

Antimicrobial resistance (AMR) in pathogenic microorganisms with multidrug resistance (MDR) constitutes a severe threat to human health. A major causative mechanism of AMR is mediated through the multidrug efflux pump (MEP). The resistance-nodulation-division superfamily (RND family) of Gram-negative bacteria is usually the major cause of MDR in clinical studies. In Salmonella enterica, the RND pump is translated from the acrAB gene, which is regulated by the activator RamA. Many MEP-caused AMR strains have high ramA gene expression due to mutations in RamR, which has a homodimeric structure comprising the dimerization domain and DNA-binding domain (DBD). Three mutations on the dimerization domain, namely Y59H, M84I, and E160D, are far from the DBD; the molecular mechanism through which they influence RamR’s binding affinity to the ramA gene promoter and consequently disrupt RamA remains unclear. The present study conducted molecular dynamics simulations, binding free energy calculations, and normal mode analysis to investigate the mechanism through which Y59H, M84I, and E160D mutations on the dimerization domain influence the binding affinity of RamR to the ramA promoter. The present results suggest that the three mutations alter the RamR structure, resulting in decreased DNA-binding affinity.