Free Energy Perturbation Calculations Research Articles

The Discovery of MORF-627, a Highly Selective Conformationally-Biased Zwitterionic Integrin αvβ6 Inhibitor for Fibrosis.

Inhibition of integrin αvβ6 is a promising approach to the treatment of fibrotic disease such as idiopathic pulmonary fibrosis. Screening a small library combining head groups that stabilize the bent-closed conformation of integrin αIIbβ3 with αv integrin binding motifs resulted in the identification of hit compounds that bind the bent-closed conformation of αvβ6. Crystal structures of these compounds bound to αvβ6 and related integrins revealed opportunities to increase potency and selectivity, and these efforts were accelerated using accurate free energy perturbation (FEP+) calculations. Optimization of PK parameters including permeability, bioavailability, clearance, and half-life resulted in the discovery of development candidate MORF-627, a highly selective inhibitor of αvβ6 that stabilizes the bent-closed conformation and has good oral PK. Unfortunately, the compound showed toxicity in a 28-day NHP safety study, precluding further development. Nevertheless, MORF-627 is a useful tool compound for studying the biology of integrin αvβ6.

Coarse-Graining the Recognition of a Glycolipid by the C-Type Lectin Mincle Receptor.

Macrophage inducible Ca2+-dependent lectin (Mincle) receptor recognizes Mycobacterium tuberculosis glycolipids to trigger an immune response. This host membrane receptor is thus a key player in the modulation of the immune response to infection by M. tuberculosis and has emerged as a promising target for the development of new vaccines against tuberculosis. The recent development of the Martini 3 force field for coarse-grained (CG) molecular modeling allows the study of interactions of soluble proteins with small ligands which was not typically modeled well with the previous Martini 2 model. Here, we present a refined approach detailing a protocol for modeling interactions between a glycolipid and its receptor at a CG level using the Martini 3 force field. Using this approach, we studied Mincle and identified critical parameters governing ligand recognition, such as loop flexibility and the regulation of hydrophobic groove formation by calcium ions. In addition, we assessed ligand affinity using free energy perturbation calculations. Our results offer mechanistic insight into the interactions between Mincle and glycolipids, providing a basis for the rational design of molecules targeting this type of membrane receptors.

Alchemical Transformations and Beyond: Recent Advances and Real-World Applications of Free Energy Calculations in Drug Discovery.

Computational methods constitute efficient strategies for screening and optimizing potential drug molecules. A critical factor in this process is the binding affinity between candidate molecules and targets, quantified as binding free energy. Among various estimation methods, alchemical transformation methods stand out for their theoretical rigor. Despite challenges in force field accuracy and sampling efficiency, advancements in algorithms, software, and hardware have increased the application of free energy perturbation (FEP) calculations in the pharmaceutical industry. Here, we review the practical applications of FEP in drug discovery projects since 2018, covering both ligand-centric and residue-centric transformations. We show that relative binding free energy calculations have steadily achieved chemical accuracy in real-world applications. In addition, we discuss alternative physics-based simulation methods and the incorporation of deep learning into free energy calculations.

Effects of Normal and Lateral Electric Fields on Membrane Mechanical Properties.

As a core component of biological and synthetic membranes, lipid bilayers are key to compartmentalizing chemical processes. Bilayer morphology and mechanical properties are heavily influenced by electric fields, such as those caused by biological ion concentration gradients. We present atomistic simulations exploring the effects of electric fields applied normally and laterally to lipid bilayers. We find that normal fields decrease membrane tension, while lateral fields increase it. Free energy perturbation calculations indicate the importance of dipole-dipole interactions to these tension changes, especially for lateral fields. We additionally show that membrane area compressibilities can be related to their cohesive energies, allowing us to estimate changes in membrane bending rigidity under applied fields. We find that normal and lateral fields decrease and increase bending rigidity, respectively. These results point to the use of directed electric fields to locally control membrane stiffness, thereby modulating associated cellular processes.

Computational discovery of tripeptide inhibitors targeting monkeypox virus A42R profilin-like protein

Monkeypox is an infectious disease caused by the monkeypox virus (MPXV), a member of the Orthopoxvirus genus closely related to smallpox. The structure of the A42R profilin-like protein is the first and only available structure among MPXV proteins. Biochemical studies of A42R were conducted in the 1990s and later work also analyzed the protein's function in viral replication in cells. This study aims to screen tripeptides for their potential inhibition of the A42R profilin-like protein using computational methods, with implications for MPXV therapy. A total of 8000 tripeptides underwent molecular docking simulations, resulting in the identification of 20 compounds exhibiting strong binding affinity to A42R. To validate the docking results, molecular dynamics simulations and free energy perturbation calculations were performed. These analyses revealed two tripeptides with sequences TRP-THR-TRP and TRP-TRP-TRP, which displayed robust binding affinity to A42R. Markedly, electrostatic interactions predominated over van der Waals interactions in the binding process between tripeptides and A42R. Three A42R residues, namely Glu9, Ser12, and Arg38, appear to be pivotal in mediating the interaction between A42R and the tripeptide ligands. Notably, tripeptides containing two or three tryptophan residues demonstrate a pronounced binding affinity, with the tripeptide comprising three tryptophan amino acids showing the highest level of affinity. These findings offer valuable insights for the selection of compounds sharing a similar structure and possessing a high affinity for A42R, potentially capable of inhibiting its enzyme activity. The study highlights a structural advantage and paves the way for the development of targeted therapies against MPXV infections.

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.

Exploring ligands that target von Willebrand factor selectively under oxidizing conditions through docking and molecular dynamics simulations.

The blood protein von Willebrand factor (VWF) is a large multimeric protein that, when activated, binds to blood platelets, tethering them to the site of vascular injury and initiating blood coagulation. This process is critical for the normal hemostatic response, but especially under inflammatory conditions, it is thought to be a major player in pathological thrombus formation. For this reason, VWF has been the target for the development of anti-thrombotic therapeutics. However, it is challenging to prevent pathological thrombus formation while still allowing normal physiological blood coagulation, as currently available anti-thrombotic therapeutics are known to cause unwanted bleeding, in particular intracranial hemorrhage. This work explores the possibility of inhibiting VWF selectively under the inflammatory conditions present during pathological thrombus formation. In particular, the A2 domain of VWF is known to inhibit the neighboring A1 domain from binding to the platelet surface receptor GpIbα, and this auto-inhibitory mechanism has been shown to be removed by oxidizing agents released during inflammation. Hence, finding drug molecules that bind at the interface between A1 and A2 only under oxidizing conditions could restore such an auto-inhibitory mechanism. Here, by using a combination of computational docking, molecular dynamics simulations, and free energy perturbation calculations, a ligand from the ZINC15 database was identified that binds at the A1A2 interface, with the interaction being stronger under oxidizing conditions. The results provide a framework for the discovery of drug molecules that bind to a protein selectively in the presence of inflammatory conditions.

Distinct roles of the major binding residues in the cation-binding pocket of the melibiose transporter MelB

Salmonella enterica serovar Typhimurium melibiose permease (MelBSt) is a prototype of the major facilitator superfamily (MFS) transporters, which play important roles in human health and diseases. MelBSt catalyzed the symport of galactosides with either Na+, Li+, or H+, but prefers the coupling with Na+. Previously, we determined the structures of the inward- and outward-facing conformation of MelBSt and the molecular recognition for galactoside and Na+. However, the molecular mechanisms for H+- and Na+-coupled symport still remain poorly understood. In this study, we solved two x-ray crystal structures of MelBSt, the cation-binding site mutants D59C at an unliganded apo-state and D55C at a ligand-bound state, and both structures display the outward-facing conformations virtually identical as published. We determined the energetic contributions of three major Na+-binding residues for the selection of Na+ and H+ by free energy simulations. Transport assays showed that the D55C mutant converted MelBSt to a solely H+-coupled symporter, and together with the free-energy perturbation calculation, Asp59 is affirmed to be the sole protonation site of MelBSt. Unexpectedly, the H+-coupled melibiose transport exhibited poor activities at greater bulky ΔpH and better activities at reversal ΔpH, supporting the novel theory of transmembrane-electrostatically localized protons and the associated membrane potential as the primary driving force for the H+-coupled symport mediated by MelBSt. This integrated study of crystal structure, bioenergetics, and free energy simulations, demonstrated the distinct roles of the major binding residues in the cation-binding pocket of MelBSt.

Ionic hydration at ambient and higher pressures: Computed chemical potentials from simulations and finite-size effects

For the chemical potential of a hydrated mono-atomic ion, finite- size effects on simulation results obtained using molecular potential truncation are investigated. Free energy perturbation (FEP) calculations were carried out by scaling in two processes the Lennard-Jones (LJ) ion-water parameters and the ion-charge (q) for Br−, K+ and Ca2+ interacting with TIP4P water. Corrections which scale with q2 enable us to reduce finite-size effects. However, at ambient conditions, discrepancies which depend on q are shown by the corrected values when comparison is made with the experimental data of the Marcus compilation. Similar behavior was observed by extrapolating the original FEP results to an infinitely large system. Hence, these errors were assumed to depend on water density and corrected at high pressures. Consistency, within statistical uncertainties, is shown when comparing with results derived from computed volumetric quantities. Results are also compared with those derived from experimental values of excess volumes at ambient conditions.

Origin of Enantioselectivity in Engineered Cytochrome c-Catalyzed Carbon-Radical FePP Hydrolysis Revealed Using QM/MM (ABEEM Polarizable Force Field) and MD Simulations.

The origin of highly efficient asymmetric aminohydroxylation of styrene catalyzed by engineered cytochrome c is investigated by the developed Atom-Bond Electronegativity Equalization Method polarizable force field (ABEEM PFF), which is a combined outcome of electronic and steric effects. Model molecules were used to establish the charge parameters of the ABEEM PFF, for which the bond-stretching and angle-bending parameters were obtained by using a combination of modified Seminario and scan methods. The interactions between carbon-radical Fe-porphyrin (FePP) and waters are simulated by molecular dynamics, which shows a clear preference for the pre-R over the pre-S. This preference is attributed to the hydrogen-bond between the mutated 100S and 101P residues as well as van der Waals interactions, enforcing a specific conformation of the carbon-radical FePP complex within the binding pocket. Meanwhile, the hydrogen-bond between water and the nitrogen atom in the active intermediate dictates the stereochemical outcome. Quantum mechanics/molecular mechanics (QM/MM (ABEEM PFF)) and free-energy perturbation calculations elucidate that the 3RTS is characterized by sandwich-like structure among adjacent amino acid residues, which exhibits greater stability than crowed arrangement in 3STS and enables the R enantiomer to form more favorably. Thus, this study provides mechanistic insight into the catalytic reaction of hemoproteins.

Fastlomap: faster lead optimization mapper algorithm for large-scale relative free energy perturbation

In recent years, free energy perturbation calculations have garnered increasing attention as tools to support drug discovery. The lead optimization mapper (Lomap) was proposed as an algorithm to calculate the relative free energy between ligands efficiently. However, Lomap requires checking whether each edge in the FEP graph is removable, which necessitates checking the constraints for all edges. Consequently, conventional Lomap requires significant computation time, at least several hours for cases involving hundreds of compounds, and is impractical for cases with more than tens of thousands of edges. In this study, we aimed to reduce the computational cost of Lomap to enable the construction of FEP graphs for hundreds of compounds. We can reduce the overall number of constraint checks required from an amount dependent on the number of edges to one dependent on the number of nodes by using the chunk check process to check the constraints for as many edges as possible simultaneously. Based on the analysis of the execution profiles, we also improved the speed of cycle constraint and diameter constraint checks. Moreover, the output graph is the same as that obtained using the conventional Lomap, enabling direct replacement of the original one with our method. With our improvement, the execution was hundreds of times faster than that of the original Lomap.

Structure-aware dual-target drug design through collaborative learning of pharmacophore combination and molecular simulation.

Dual-target drug design has gained significant attention in the treatment of complex diseases, such as cancers and autoimmune disorders. A widely employed design strategy is combining pharmacophores to leverage the knowledge of structure-activity relationships of both targets. Unfortunately, pharmacophore combination often struggles with long and expensive trial and error, because the protein pockets of the two targets impose complex structural constraints. In this study, we propose AIxFuse, a structure-aware dual-target drug design method that learns pharmacophore fusion patterns to satisfy the dual-target structural constraints simulated by molecular docking. AIxFuse employs two self-play reinforcement learning (RL) agents to learn pharmacophore selection and fusion by comprehensive feedback including dual-target molecular docking scores. Collaboratively, the molecular docking scores are learned by active learning (AL). Through collaborative RL and AL, AIxFuse learns to generate molecules with multiple desired properties. AIxFuse is shown to outperform state-of-the-art methods in generating dual-target drugs against glycogen synthase kinase-3 beta (GSK3β) and c-Jun N-terminal kinase 3 (JNK3). When applied to another task against retinoic acid receptor-related orphan receptor γ-t (RORγt) and dihydroorotate dehydrogenase (DHODH), AIxFuse exhibits consistent performance while compared methods suffer from performance drops, leading to a 5 times higher performance in success rate. Docking studies demonstrate that AIxFuse can generate molecules concurrently satisfying the binding mode required by both targets. Further free energy perturbation calculation indicates that the generated candidates have promising binding free energies against both targets.

Development and Comprehensive Benchmark of a High-Quality AMBER-Consistent Small Molecule Force Field with Broad Chemical Space Coverage for Molecular Modeling and Free Energy Calculation.

Biomolecular simulations have become an essential tool in contemporary drug discovery, and molecular mechanics force fields (FFs) constitute its cornerstone. Developing a high quality and broad coverage general FF is a significant undertaking that requires substantial expert knowledge and computing resources, which is beyond the scope of general practitioners. Existing FFs originate from only a limited number of groups and organizations, and they either suffer from limited numbers of training sets, lower than desired quality because of oversimplified representations, or are costly for the molecular modeling community to access. To address these issues, in this work, we developed an AMBER-consistent small molecule FF with extensive chemical space coverage, and we provide Open Access parameters for the entire modeling community. To validate our FF, we carried out benchmarks of quantum mechanics (QM)/molecular mechanics conformer comparison and free energy perturbation calculations on several benchmark data sets. Our FF achieves a higher level of performance at reproducing QM energies and geometries than two popular open-source FFs, OpenFF2 and GAFF2. In relative binding free energy calculations for 31 protein-ligand data sets, comprising 1079 pairs of ligands, the new FF achieves an overall root-mean-square error of 1.19 kcal/mol for ΔΔG and 0.92 kcal/mol for ΔG on a subset of 463 ligands without bespoke fitting to the data sets. The results are on par with those of the leading commercial series of OPLS FFs.

Adaptive lambda schemes for efficient relative binding free energy calculation.

The relative free energy perturbation (RFEP) calculation is one of the most theoretically sound computational chemistry approaches for the binding affinity prediction. However, its application is often hindered by the complexity of the calculation choices and the requirement of expertise in the field. Improper lambda scheme of RFEP may result in deviations from an accurate description of the perturbation process and is prone to erroneous affinity predictions. To address such challenges, an automated adaptive lambda method is proposed where the adaptive lambda schemes are obtained through a split-and-merge algorithm based on the pilot runs. The newly established workflow along with a series of improvements to the perturbation settings increases the consistency of the RFEP calculation results. Comparing the pilot and adaptive lambda schemes, the latter demonstrated improvements in convergence and reproducibility and lowered the mean unsigned error and the root-mean-square error. Overall, the adaptive lambda method is a reliable and robust choice to predict small molecule relative binding free energy and can be capitalized to benefit routine RFEP calculations for drug discovery projects.

OPLS/2020 Force Field for Unsaturated Hydrocarbons, Alcohols, and Ethers.

The OPLS all-atom force field was updated and applied to modeling unsaturated hydrocarbons, alcohols, and ethers. Testing has included gas-phase conformational energetics, properties of pure liquids, and free energies of hydration. Monte Carlo statistical mechanics (MC) calculations were used to model 60 liquids. In addition, a robust, automated procedure was devised to compute the free energies of hydration with high precision via free-energy perturbation (FEP) calculations using double annihilation. Testing has included larger molecules than in the past, and parameters are reported for the first time for some less common groups including alkynes, allenes, dienes, and acetals. The average errors in comparison with experimental data for the computed properties of the pure liquids were improved with the modified force field (OPLS/2020). For liquid densities and heats of vaporization, the average unsigned errors are 0.01 g/cm3 and 0.2 kcal/mol. The average error and signed error for free energies of hydration are both 1.2 kcal/mol. As noted before, this reflects a systematic overestimate of the hydrophobicity of organic molecules when the parametrization is done to minimize the errors for properties of pure liquids. Implications for the modeling of biomolecular systems with standard force fields are considered.

Accurate Computation of Thermodynamic Activation Parameters in the Chorismate Mutase Reaction from Empirical Valence Bond Simulations.

Chorismate mutase (CM) enzymes have long served as model systems for benchmarking new methods and tools in computational chemistry. Despite the enzymes' prominence in the literature, the extent of the roles that activation enthalpy and entropy play in catalyzing the conversion of chorismate to prephenate is still subject to debate. Knowledge of these parameters is a key piece in fully understanding the mechanism of chorismate mutases. Within this study, we utilize EVB/MD free energy perturbation calculations at a range of temperatures, allowing us to extract activation enthalpies and entropies from an Arrhenius plot of activation free energies of the reaction catalyzed by a monofunctional Bacillus subtilis CM and the promiscuous enzyme isochorismate pyruvate lyase of Pseudomonas aeruginosa. In comparison to the uncatalyzed reaction, our results show that both enzyme-catalyzed reactions exhibit a substantial reduction in activation enthalpy, while the effect on activation entropy is relatively minor, demonstrating that enzyme-catalyzed CM reactions are enthalpically driven. Furthermore, we observe that the monofunctional CM from B. subtilis more efficiently catalyzes this reaction than its promiscuous counterpart. This is supported by a structural analysis of the reaction pathway at the transition state, from which we identified key residues explaining the enthalpically driven nature of the reactions and also the difference in efficiencies between the two enzymes.

Discovery of (R)-4-(8-methoxy-2-methyl-1-(1-phenylethy)-1H-imidazo[4,5-c]quinnolin-7-yl)-3,5-dimethylisoxazole as a potent and selective BET inhibitor for treatment of acute myeloid leukemia (AML) guided by FEP calculation

The functions of the bromodomain and extra terminal (BET) family of proteins have been proved to be involved in various diseases, particularly the acute myeloid leukemia (AML). In this work, guided by free energy perturbation (FEP) calculation, a methyl group was selected to be attached to the 1H-imidazo[4,5-c]quinoline skeleton, and a series of congeneric compounds were synthesized. Among them, compound 10 demonstrated outstanding activity against BRD4 BD1 with an IC50 value of 1.9 nM and exhibited remarkable antiproliferative effects against MV4-11 cells. The X-ray cocrystal structure proved that 10 occupied the acetylated lysine (KAc) binding cavity and the WPF shelf of BRD4 BD1. Additionally, 10 displayed high selectivity towards BET family members, effectively inhibiting the growth of AML cells, promoting apoptosis, and arresting the cell cycle at the G0/G1 phase. Further mechanistic studies demonstrated that compound 10 could suppress the expression of c-Myc and CDK6 while enhancing the expression of P21, PARP, and cleaved PARP. Moreover, 10 exhibited remarkable pharmacokinetic properties and significant antitumor efficacy in vivo. Therefore, compound 10 may represent a new, potent and selective BET bromodomain inhibitor for the development of therapeutics to treat AML.

Loops Mediate Agonist-Induced Activation of the Stimulator of Interferon Genes Protein.

The stimulator of interferon genes (STING) is an important therapeutic target for cancer diseases. The activated STING recruits downstream tank-binding kinase 1 (TBK1) to trigger several important immune responses. However, the molecular mechanism of how agonist molecules mediate the STING-TBK1 interactions remains elusive. Here, we performed molecular dynamics simulations to capture the conformational changes of STING and TBK1 upon agonist binding. Our simulations revealed that multiple helices (α5-α7) and especially three loops (loop 6, loop 8, and C-terminal tail) of STING participated in the allosteric mediation of the STING-TBK1 interactions. Consistent results were also observed in the simulations of the constitutive activating mutant of STING (R284S). We further identified α5 as a key region in this agonist-induced activation mechanism of STING. Free-energy perturbation calculations of multiple STING agonists demonstrated that an alkynyl group targeting α5 is a determinant for agonist activities. These results not only offer deeper insights into the agonist-induced allosteric mediation of STING-TKB1 interactions but also provide a guidance for future drug development of this important therapeutic target.

Dynamic extracellular vestibule of human SERT: Unveiling druggable potential with high-affinity allosteric inhibitors

The serotonin transporter (SERT) tightly regulates synaptic serotonin levels and has been the primary target of antidepressants. Binding of inhibitors to the allosteric site of human SERT (hSERT) impedes the dissociation of antidepressants bound at the central site and may enhance the efficacy of such antidepressants to potentially reduce their dosage and side effects. Here, we report the identification of a series of high-affinity allosteric inhibitors of hSERT in a unique scaffold, with the lead compound, Lu AF88273 (3-(1-(2-(1H-indol-3-yl)ethyl)piperidin-4-yl)-6-chloro-1H-indole), having 2.1 nM allosteric potency in inhibiting imipramine dissociation. In addition, we find that Lu AF88273 also inhibits serotonin transport in a noncompetitive manner. The binding pose of Lu AF88273 in the allosteric site of hSERT is determined with extensive molecular dynamics simulations and rigorous absolute binding free energy perturbation (FEP) calculations, which show that a part of the compound occupies a dynamically formed small cavity. The predicted binding location and pose are validated by site-directed mutagenesis and can explain much of the structure-activity relationship of these inhibitors using the relative binding FEP calculations. Together, our findings provide a promising lead compound and the structural basis for the development of allosteric drugs targeting hSERT. Further, they demonstrate that the divergent allosteric sites of neurotransmitter transporters can be selectively targeted.

Structure-Based Ligand Discovery Targeting the Transmembrane Domain of Frizzled Receptor FZD7.

Despite the essential roles of Frizzled receptors (FZDs) in mediating Wnt signaling in embryonic development and tissue homeostasis, ligands targeting FZDs are rare. A few antibodies and peptide modulators have been developed that mainly bind to the family-conserved extracellular cysteine-rich domain of FZDs, while the canonical binding sites in the transmembrane domain (TMD) are far from sufficiently addressed. Based on the recent structures of FZDs, we explored small-molecule ligand discovery by targeting TMD. From the ChemDiv library with ∼1.6 million compounds, we identified compound F7H as an antagonist of FZD7 with an IC50 at 1.25 ± 0.38 μM. Focusing on this hit, the structural dissection study, together with computing studies such as molecular docking, molecular dynamics simulation, and free energy perturbation calculations, defined the binding pocket with key residue recognition. Our results revealed the structural basis of ligand recognition and demonstrated the feasibility of structure-guided ligand discovery for FZD7-TMD.