Free Energy Perturbation Research Articles

Computer-aided drug design approaches applied toscreen natural product's structural analogs targeting arginase in Leishmania spp.

Introduction: Leishmaniasis is a disease with high mortality rates and approximately 1.5 million new cases each year. Despite the new approaches and advances to fight the disease, there are no effective therapies. Methods: Hence, this study aims to screen for natural products' structural analogs as new drug candidates against leishmaniasis. We applied Computer-aided drug design (CADD) approaches, such as virtual screening, molecular docking, molecular dynamics simulation, molecular mechanics-generalized Born surface area (MM-GBSA) binding free estimation, and free energy perturbation (FEP) aiming to select structural analogs from natural products that have shown anti-leishmanial and anti-arginase activities and that could bind selectively against the Leishmania arginase enzyme. Results: The compounds 2H-1-benzopyran, 3,4-dihydro-2-(2-methylphenyl)-(9CI), echioidinin, and malvidin showed good results against arginase targets from three parasite species and negative results for potential toxicities. The echioidinin and malvidin ligands generated interactions in the active center at pH 2.0 conditions by MM-GBSA and FEP methods. Conclusions: This work suggests the potential anti-leishmanial activity of the compounds and thus can be further in vitro and in vivo experimentally validated.

Deciphering the Anomalous Acidic Tendency of Terminal Water at Rutile(110)–Water Interfaces

Understanding the mechanism of the oxygen evolution reaction (OER) is essential to improve the efficiency of photocatalysis for TiO2. Previous studies have highlighted the importance of terminal hydroxide radical (TiOH•) in the OER. Ab initio molecular dynamics simulations (AIMD) with hybrid functional have revealed that this radical readily loses its proton, creating the key intermediate, oxygen radical anion (Ti5cO•–). Herein, we combine machine-learning-accelerated molecular dynamics with density functional theory calculations to demonstrate that the Ti5cO•– can alternatively be generated through the trapping of a hole in a terminal oxygen anion (Ti5cO2–) at rutile(110)–water interfaces. Further examination reveals that the Ti5cO2– results from the deprotonation of Ti5cOH– and remains stable at the charge-neutral interfaces for a transient time period of ca. 100 ps. The AIMD-based free energy perturbation method predicts that the acidity constant of Ti5cOH– is even smaller than that of Ti5cOH2, thereby rationalizing the stability of Ti5cO2–. Structural analyses show that this anomalous acidic tendency of terminal water originates from the decrease of Ti–O bond length and the transition of Titanium’s coordination from octahedral to pyramidal in Ti5cO2–. Our findings provide valuable insights into the surface acid–base chemistry and a potential explanation for the pH-dependent behavior of photogenerated holes for TiO2.

Novel Allosteric Effectors Targeting Human Transcription Factor TEAD.

The Hippo pathway is an evolutionary conserved signaling network involved in several cellular regulatory processes. Dephosphorylation and overexpression of Yes-associated proteins (YAPs) in the Hippo-off state are common in several types of solid tumors. YAP overexpression results in its nuclear translocation and interaction with transcriptional enhanced associate domain 1-4 (TEAD1-4) transcription factors. Covalent and non-covalent inhibitors have been developed to target several interaction sites between TEAD and YAP. The most targeted and effective site for these developed inhibitors is the palmitate-binding pocket in the TEAD1-4 proteins. Screening of a DNA-encoded library against the TEAD central pocket was performed experimentally to identify six new allosteric inhibitors. Inspired by the structure of the TED-347 inhibitor, chemical modification was performed on the original inhibitors by replacing secondary methyl amide with a chloromethyl ketone moiety. Various computational tools, including molecular dynamics, free energy perturbation, and Markov state model analysis, were employed to study the effect of ligand binding on the protein conformational space. Four of the six modified ligands were associated with enhanced allosteric communication between the TEAD4 and YAP1 domains indicated by the relative free energy perturbation to original molecules. Phe229, Thr332, Ile374, and Ile395 residues were revealed to be essential for the effective binding of the inhibitors.

Alchemical Free Energy and Hamiltonian Replica Exchange Molecular Dynamics to Compute Hydrofluorocarbon Isotherms in Imidazolium-Based Ionic Liquids.

Ionic liquids (ILs) have shown promise for applications that leverage differential gas solubility in an IL solvent, e.g., gas separations. Although most available literature provides Henry's law constants, the ability to efficiently estimate full isotherms is important for engineering design calculations. Molecular simulation can be used as a tool to predict full isotherms of gas in ILs. However, particle insertions or deletions in a charge-dense IL medium and the sluggish conformational dynamics of ILs present two sampling challenges for these systems. We therefore devised a method that uses Hamiltonian replica exchange (HREX) molecular dynamics (MD) combined with alchemical free energy calculations to compute full solubility isotherms of two different hydrofluorocarbons (HFCs) in imidazolium-based IL binary mixtures. This workflow is significantly faster than the Gibbs ensemble Monte Carlo (GEMC) simulations which fail to deal with the slow conformational relaxation caused by the sluggish dynamics of ILs. Multiple free energy estimators, including thermodynamic integration, free energy perturbation, and multistate Bennett acceptance ratio method, provided consistent results. Overall, the simulated Henry's law constant, isotherm curvature, and solubility trends match experimental results reasonably well. We close by calculating the full solubility isotherms of two HFCs in IL mixtures that have not been reported in the literature, demonstrating the potential of this method to be used for solubility prediction and setting the stage for future computational screening studies that search for the "best" IL to separate azeotropic HFC mixtures.

Enhancing Hit Discovery in Virtual Screening through Absolute Protein-Ligand Binding Free-Energy Calculations.

In the hit identification stage of drug discovery, a diverse chemical space needs to be explored to identify initial hits. Contrary to empirical scoring functions, absolute protein-ligand binding free-energy perturbation (ABFEP) provides a theoretically more rigorous and accurate description of protein-ligand binding thermodynamics and could, in principle, greatly improve the hit rates in virtual screening. In this work, we describe an implementation of an accurate and reliable ABFEP method in FEP+. We validated the ABFEP method on eight congeneric compound series binding to eight protein receptors including both neutral and charged ligands. For ligands with net charges, the alchemical ion approach is adopted to avoid artifacts in electrostatic potential energy calculations. The calculated binding free energies correlate with experimental results with a weighted average of R2 = 0.55 for the entire dataset. We also observe an overall root-mean-square error (RMSE) of 1.1 kcal/mol after shifting the zero-point of the simulation data to match the average experimental values. Through ABFEP calculations using apo versus holo protein structures, we demonstrated that the protein conformational and protonation state changes between the apo and holo proteins are the main physical factors contributing to the protein reorganization free energy manifested by the overestimation of raw ABFEP calculated binding free energies using the holo structures of the proteins. Furthermore, we performed ABFEP calculations in three virtual screening applications for hit enrichment. ABFEP greatly improves the hit rates as compared to docking scores or other methods like metadynamics. The good performance of ABFEP in rank ordering compounds demonstrated in this work confirms it as a useful tool to improve the hit rates in virtual screening, thus facilitating hit discovery.

Fragment Linker Prediction Using the Deep Encoder-Decoder Network for PROTACs Drug Design.

A drug discovery and development pipeline is a prolonged and complex process that remains challenging for both computational methods and medicinal chemists and has not been able to be resolved using computational methods. Deep learning has been utilized in various fields and achieved tremendous success in designing novel molecules in the pharmaceutical industry. Herein, we use state-of-the-art techniques to propose a deep neural network, AIMLinker, to rapidly design and generate meaningful drug-like proteolysis targeting chimeras (PROTACs) analogs. The model extracts the structural information from the input fragments and generates linkers to incorporate them. We integrate filters in the model to exclude nondruggable structures guided via protein-protein complexes while retaining molecules with potent chemical properties. The novel PROTACs subsequently pass through molecular docking, taking root-mean-square deviation (RMSD), relative Gibbs free energy (ΔΔGbinding), molecular dynamics (MD) simulation, and free energy perturbation (FEP) calculations as the measurement criteria for testing the robustness and feasibility of the model. The generated novel PROTACs molecules possess similar structural information with superior binding affinity to the binding pockets compared to the existing CRBN-dBET6-BRD4 ternary complexes. We demonstrate the effectiveness of the methodology of leveraging AIMLinker to design novel compounds for PROTACs molecules exhibiting better chemical properties compared to the dBET6 crystal pose.

Topical selective TYK2 inhibition attenuates inflammation and clinical symptoms in a mouse model of human psoriasis

Abstract Psoriasis is a chronic skin condition, that affects ~8M patients in the United States. A majority (~7M) of patients suffer from mild to moderate form of disease that lacks safe & effective treatment options. Although antibodies targeting Th1/Th17 cytokines (eg IL12/IL23/IL17) are highly effective in moderate to severe psoriasis, they are not prescribed to mild to moderate patients due to high costs & safety concerns. There is a huge unmet need for developing safe and effective non-steroidal topical therapies. Small molecule inhibitors of Th1/Th17 pathways with optimal topical drug like properties can effectively address patient needs; making TYK2 inhibitors (TYK2i) attractive candidates for treating mild to moderate psoriasis. However, given the high degree of structural similarity among JAK family, designing topical selective TYK2i have remained a challenge. Using an innovative structure based approach, we have designed, synthesized and characterized TYK2i optimized for topical drug like properties and JAK family selectivity using free energy perturbation (FEP) methods and medchem SAR. We have identified TYK2i with pM to low nM potencies that are highly selective for TYK2 among 468 kinases. These analogs have enhanced topical bioavailability with suitable drug like properties. In a mouse model of psoriasis, these analogs dose-dependently suppressed PASI clinical scores, ear & back skin thickness and inflammation. Additionally, skin histology showed significant attenuation of epidermal thickening and pSTAT activation suggesting that topical TYK2 inhibition retains anti-inflammatory activity while eliminating potential for dose-limiting side effects observed with systemic non-selective TYK2 or other JAK inhibitors. This research has been fully funded by Alembic Pharmaceuticals Ltd., including salary support to all authors.

Computing chemical potentials of adsorbed or confined fluids.

The chemical potential of adsorbed or confined fluids provides insight into their unique thermodynamic properties and determines adsorption isotherms. However, it is often difficult to compute this quantity from atomistic simulations using existing statistical mechanical methods. We introduce a computational framework that utilizes static structure factors, thermodynamic integration, and free energy perturbation for calculating the absolute chemical potential of fluids. For demonstration, we apply the method to compute the adsorption isotherms of carbon dioxide in a metal-organic framework and water in carbon nanotubes.

SophosQM: Accurate Binding Affinity Prediction in Compound Optimization.

The optimization of compounds' binding affinity for a biological target is a crucial aspect of the drug development process. Being able to accurately predict binding energies in advance of synthesizing compounds would have a massive impact on the speed of the drug discovery process. The ideal binding affinity prediction method should combine accuracy, reliability, and speed. In this paper, we present SophosQM, a quantum mechanics (QM)-based approach, which can accurately predict the binding affinities of compounds to proteins. The binding affinity predictive models generated by SophosQM are based on the fragment molecular orbital (FMO) method to estimate the enthalpic component of the binding free energy, and a macroscopic descriptor, clog P, is used as an approximation of the entropic component. The affinity prediction is performed using multilinear regression, fitting the experimental values against the FMO-computed enthalpic term and clog P. The quality of the prediction can be assessed in terms of the correlation coefficient between experimental and predicted values. In this work, the method's reliability and accuracy are exemplified by applying SophosQM to 70 compounds binding to six different targets of pharmaceutical relevance. Overall, the results show a very satisfactory performance with a global correlation coefficient in the order of 0.9. Our predictions also show a satisfactory performance compared to data based on free energy perturbation. Finally, SophosQM can also be applied in high-throughput mode by using semiempirical QM methods to evaluate large portions of chemical space, while retaining a good level of accuracy, but decreasing the computing time to just a few seconds per compound.

Combining QM/MM Calculations with Classical Mining Minima to Predict Protein-Ligand Binding Free Energy.

We developed an effective binding free energy prediction protocol which incorporates quantum mechanical/molecular mechanical (QM/MM) calculations to substitute the specified atomic charges of force fields with quantum-mechanically recalculated ones at a proposed pose using a mining minima approach with the VeraChem mining minima engine. We tested this protocol using seven well-known targets with 147 different ligands and compared it with classical mining minima and the most popular binding free energy (BFE) methods using different metrics. Our new protocol, dubbed Qcharge-VM2, yielded an overall Pearson correlation of 0.86, which was better than all the methods examined. Qcharge-VM2 performed significantly better than implicit solvent-based methods, such as MM-GBSA and MM-PBSA, but not as good as explicit water-based free energy perturbation methods, such as FEP+, in terms of root-mean-square error, RMSE (1.75 kcal/mol) and mean unsigned error, MUE (1.39 kcal/mol) on a limited set of targets. However, our protocol is substantially less computationally demanding compared with FEP+. The combined accuracy and efficiency of our method can be valuable in drug discovery campaigns.

Enhanced Sampling for Free Energy Profiles with Post-Transition-State Bifurcations

We present a strategy to explore the free energy landscapes of chemical reactions with post-transition-state bifurcations using an enhanced sampling method based on well-tempered metadynamics. Obviating the need for computationally expensive density functional theory-level ab initio molecular dynamics simulations, we obtain accurate energetics by utilizing a free energy perturbation scheme and deep learning estimator for the single-point energies of substrate configurations. Using a pair of easily interpretable collective variables, we present a quantitative free energy surface that is compatible with harmonic transition state theory calculations and in which the bifurcations are clearly visible. We demonstrate our approach with the example of the SpnF-catalyzed Diels-Alder reaction, a cycloaddition reaction in which post-transition-state bifurcation leads to the [4+2] as well as the [6+4] cycloadduct. We obtain the free energy landscapes for different stereochemical reaction pathways and characterize the mechanistic continuum between relevant reaction channels without explicitly searching for the pertinent transition state structures.

Scaffold Hopping and Optimization of Small Molecule Soluble Adenyl Cyclase Inhibitors Led by Free Energy Perturbation.

Free energy perturbation is a computational technique that can be used to predict how small changes to an inhibitor structure will affect the binding free energy to its target. In this paper, we describe the utility of free energy perturbation with FEP+ in the hit-to-lead stage of a drug discovery project targeting soluble adenyl cyclase. The project was structurally enabled by X-ray crystallography throughout. We employed free energy perturbation to first scaffold hop to a preferable chemotype and then optimize the binding affinity to sub-nanomolar levels while retaining druglike properties. The results illustrate that effective use of free energy perturbation can enable a drug discovery campaign to progress rapidly from hit to lead, facilitating proof-of-concept studies that enable target validation.

Implementation of solvation free energy framework to predict the vapor-liquid equilibrium behaviors for the water–hydrazine and ethanol–water azeotropic systems

An azeotrope is a special vapor–liquid equilibrium in which the composition of two coexisting phases is equal. Due to its strong deviation from the Raoult’s law, it poses a great challenge to the simple distillation process and the understanding of fundamental molecular mechanisms. In this study, we calculated the coupling work of solute species based on the solvation free energy framework to predict the vapor–liquid equilibrium (VLE) behavior and azeotrope point of mixtures. Two prototypical binary mixtures: the negative water-hydrazine and positive ethanol–water azeotrope were investigated using molecular dynamic simulations and the OPLS-AA force field. All the simulations were performed with the TIP3P water model, because it has a good capability to represent enthalpy properties. The molecular models of hydrazine and ethanol have been re-parameterized to reproduce the density and solvation free energy properties of pure substance. Using free energy perturbation simulation, the coupling work of each substance was studied against liquid mixture composition. The coupling work is also translated to x-y or T-x-y phase diagrams for comparison with available experimental data. The tendency of two azeotropes is that the more volatile components experience the less pronounced coupling work changes in their local environments as the azeotropes form. The coupling work is also divided into van der Waals and Coulomb contributions to investigate the interactions between particles. Lithium Chloride was investigated as an entrainer for the separation of ethanol–water system. The coupling works of binary mixtures were calculated at different salt concentrations to study the salt-out effect for the violate species and salt-in effect for the less violate species. By adding LiCl, it was found that the azeotrope point was efficiently eliminated at a salinity 0.05 mol fraction at 1 bar. All the results show that the simulation methodology of coupling work are reliable for calculating the liquid–vapor equilibrium behavior. This method is also valuable for selecting and designing new extractants for distillation process.

Validation of the Alchemical Transfer Method for the Estimation of Relative Binding Affinities of Molecular Series.

The accurate prediction of protein-ligand binding affinities is crucial for drug discovery. Alchemical free energy calculations have become a popular tool for this purpose. However, the accuracy and reliability of these methods can vary depending on the methodology. In this study, we evaluate the performance of a relative binding free energy protocol based on the alchemical transfer method (ATM), a novel approach based on a coordinate transformation that swaps the positions of two ligands. The results show that ATM matches the performance of more complex free energy perturbation (FEP) methods in terms of Pearson correlation but with marginally higher mean absolute errors. This study shows that the ATM method is competitive compared to more traditional methods in speed and accuracy and offers the advantage of being applicable with any potential energy function.

QSAR via Multisite λ-Dynamics in the Orphaned TSSK1B Kinase.

Multisite λ-dynamics (MSλD) is a novel method for the calculation of relative free energies of binding for ligands to their targeted receptors. It can be readily used to examine a large number of molecules with multiple functional groups at multiple sites around a common core. This makes MSλD a powerful tool in structure-based drug design. In the present study, MSλD is applied to calculate the relative binding free energies of 1296 inhibitors to the testis specific serine kinase 1B (TSSK1B), a validated target for male contraception. For this system, MSλD requires significantly fewer computational resources compared to traditional free energy methods like free energy perturbation or thermodynamic integration. From MSλD simulations, we examined whether modifications of a ligand at two different sites are coupled or not. Based on our calculations, we established a quantitative structure-activity relationship (QSAR) for this set of molecules and identified a site in the ligand where further modification, such as adding more polar groups, may lead to increased binding affinity. This article is protected by copyright. All rights reserved.

Aptamer affinity to P53 DBD: A molecular dynamics study

Protein P53, as a tumor suppressor plays an important role in many biological processes. Its DNA binding domain (DBD) contains a Zn ion which has an important role in aggregation and binding. The presence of abnormal wild-type P53 and its mutations in high concentration can be indicative of a pathological cancer state. Because of the intrinsic disordering in its structure, finding structures to overcome disordering and to interact with this protein and its mutations is one of the challenges of scientists for using P53 as a biomarker to early diagnosis of cancers. Here, we used the molecular dynamics simulation to study the interaction of wild-type P53 DBD and its mutations with an aptamer to find one solution for the mentioned challenge. The analyses of RMSD, Electrostatic energy, Number of hydrogen binding are used to discuss the quality of interactions. Aptamer- protein binding free energy is calculated with free energy perturbation method and value of ΔG= -181 Kcal.mol−1 for the interaction of aptamer- R249S, allow us to introduce the aptamer as a material to capture and recognize of R249S.

Are Deep Learning Structural Models Sufficiently Accurate for Virtual Screening? Application of Docking Algorithms to AlphaFold2 Predicted Structures.

Machine learning-based protein structure prediction algorithms, such as RosettaFold and AlphaFold2, have greatly impacted the structural biology field, arousing a fair amount of discussion around their potential role in drug discovery. While there are few preliminary studies addressing the usage of these models in virtual screening, none of them focus on the prospect of hit-finding in a real-world virtual screen with a model based on low prior structural information. In order to address this, we have developed an AlphaFold2 version where we exclude all structural templates with more than 30% sequence identity from the model-building process. In a previous study, we used those models in conjunction with state-of-the-art free energy perturbation methods and demonstrated that it is possible to obtain quantitatively accurate results. In this work, we focus on using these structures in rigid receptor-ligand docking studies. Our results indicate that using out-of-the-box Alphafold2 models is not an ideal scenario for virtual screening campaigns; in fact, we strongly recommend to include some post-processing modeling to drive the binding site into a more realistic holo model.

To Design Scalable Free Energy Perturbation Networks, Optimal Is Not Enough.

Drug discovery is accelerated with computational methods such as alchemical simulations to estimate ligand affinities. In particular, relative binding free energy (RBFE) simulations are beneficial for lead optimization. To use RBFE simulations to compare prospective ligands in silico, researchers first plan the simulation experiment, using graphs where nodes represent ligands and graph edges represent alchemical transformations between ligands. Recent work demonstrated that optimizing the statistical architecture of these perturbation graphs improves the accuracy of the predicted changes in the free energy of ligand binding. Therefore, to improve the success rate of computational drug discovery, we present the open-source software package High Information Mapper (HiMap)─a new take on its predecessor, Lead Optimization Mapper (LOMAP). HiMap removes heuristics decisions from design selection and instead finds statistically optimal graphs over ligands clustered with machine learning. Beyond optimal design generation, we present theoretical insights for designing alchemical perturbation maps. Some of these results include that for n number of nodes, the precision of perturbation maps is stable at n·ln(n) edges. This result indicates that even an "optimal" graph can result in unexpectedly high errors if a plan includes too few alchemical transformations for the given number of ligands and edges. And, as a study compares more ligands, the performance of even optimal graphs will deteriorate with linear scaling of the edge count. In this sense, ensuring an A- or D-optimal topology is not enough to produce robust errors. We additionally find that optimal designs will converge more rapidly than radial and LOMAP designs. Moreover, we derive bounds for how clustering reduces cost for designs with a constant expected relative error per cluster, invariant of the size of the design. These results inform how to best design perturbation maps for computational drug discovery and have broader implications for experimental design.

The stepwise dissociation of the Zn(II)-bound Atox1 homodimer and its energetic asymmetry

The copper chaperone Atox1, which is featured by a peculiar CxxC metal-binding motif, is a protein involved in copper trafficking to the Golgi organelle, allowing the delivery of copper ions to specific target proteins in the secretory pathway. Various experimental studies suggest that this cysteine-containing protein may also be implicated in Zn(II) ion binding, functioning as a regulatory factor. However, despite the relevance of this interaction, the mechanism of zinc ion dissociation from the Atox1 homodimer has not yet been elucidated.We herein devised a simulation workflow based on the combination of parallel bias metadynamics and free energy perturbation theory, to predict the dissociation of the zinc ion from the Atox1 homodimer. The freeenergy simulations successfully predicted the dissociation free energy differences of the process, disclosing a mechanism based on a stepwise dissociation pathway of the zinc ion from the Atox1 homodimer.

Unbiased molecular dynamics simulation of a first-in-class small molecule inhibitor binds to oncostatin M

Small molecule inhibitors (SMIs) targeting oncostatin M (OSM) signaling pathway represent new therapeutics to combat cancer, inflammatory bowel disease (IBD) and CNS disease. Recently, the first-in-class SMI named SMI-10B that target OSM and block its interaction with receptor (OSMR) were reported. However, the binding pocket and interaction mode of the compound on OSM remain poorly understood, which hampering the rational design of SMIs that target OSM. Here, using SMI-10B as a probe, the multiple pockets on OSM for small molecules binding were extensively explored by unbiased molecular dynamics (MD) simulations. Then, the near-native structure of the complex was identified by molecular mechanics generalized Born surface area (MM/GBSA) binding energy funnel. Moreover, the binding stabilities of the protein-ligand complexes in near- and non-native conformations were verified by additional independent MD runs and absolute free energy perturbation (FEP) calculation. In summary, the unique feature of SMI-10B spontaneously binds to OSM characterized here not only provide detailed information for understanding the molecular mechanism of SMI-10B binding to OSM, but also will facilitate the rational design of novel and more potent SMIs to block OSM signaling.