Using AMBER18 for Relative Free Energy Calculations

Hit-to-Lead Optimization of Heterocyclic Carbonyloxycarboximidamides as Selective Antagonists at Human Adenosine A3 Receptor

Antagonism of the human adenosine A3 receptor (hA3R) has potential therapeutic application. Alchemical relative binding free energy calculations of K18 and K32 suggested that the combination of a 3-(2,6-dichlorophenyl)-isoxazolyl group with 2-pyridinyl at the ends of a carbonyloxycarboximidamide group should improve hA3R affinity. Of the 25 new analogues synthesized, 37 and 74 showed improved hA3R affinity compared to K18 (and K32). This was further improved through the addition of a bromine group to the 2-pyridinyl at the 5-position, generating compound 39. Alchemical relative binding free energy calculations, mutagenesis studies and MD simulations supported the compounds' binding pattern while suggesting that the bromine of 39 inserts deep into the hA3R orthosteric pocket, so highlighting the importance of rigidification of the carbonyloxycarboximidamide moiety. MD simulations highlighted the importance of rigidification of the carbonyloxycarboximidamide, while suggesting that the bromine of 39 inserts deep into the hA3R orthosteric pocket, which was supported through mutagenesis studies 39 also selectively antagonized endogenously expressed hA3R in nonsmall cell lung carcinoma cells, while pharmacokinetic studies indicated low toxicity enabling in vivo evaluation. We therefore suggest that 39 has potential for further development as a high-affinity hA3R antagonist.

Identification of novel LRRK2 inhibitors by structure-based virtual screening and alchemical free energy calculation

The Leucine-rich repeat kinase 2 (LRRK2) target has been identified as a promising drug target for Parkinson's disease (PD) treatment. This study focuses on optimizing the activity of LRRK2 inhibitors using alchemical relative binding free energy (RBFE) calculations. Initially, we assessed various free energy calculation methods across different LRRK2 kinase inhibitor scaffolds. The results indicate that alchemical free energy calculations are promising for prospective predictions on LRRK2 inhibitors, especially for the aminopyrimidine scaffold with an RMSE of 1.15 kcal mol-1 and Rp of 0.83. Following this, we optimized a potent LRRK2 kinase inhibitor identified from previous virtual screenings, featuring a novel scaffold. Guided by RBFE predictions using alchemical methods, this optimization led to the discovery of compound LY2023-001. This compound, with a [1,2,4]triazolo[5,6-b]indole scaffold, exhibited enhanced inhibitory activity against G2019S LRRK2 (IC50 = 12.9 nM). Molecular dynamics (MD) simulations revealed that LY2023-001 formed stable hydrogen bonds with Glu1948, and Ala1950 in the G2019S LRRK2 protein. Additionally, its phenyl substituents engage in strong electrostatic interactions with Lys1906 and van der Waals interactions with Leu1885, Phe1890, Val1893, Ile1933, Met1947, Leu1949, Leu2001, Ala2016, and Asp2017. Our findings underscore the potential of computational methods in the successful optimization of small molecules, offering important insights for the development of novel LRRK2 inhibitors.

Ligand Binding Affinity Prediction for Membrane Proteins with Alchemical Free Energy Calculation Methods

Alchemical relative binding free energy (ΔΔG) calculations have shown high accuracy in predicting ligand binding affinity and have been used as important tools in computer-aided drug discovery and design. However, there has been limited research on the application of ΔΔG methods to membrane proteins despite the fact that these proteins represent a significant proportion of drug targets, play crucial roles in biological processes, and are implicated in numerous diseases. In this study, to predict the binding affinity of ligands to G protein-coupled receptors (GPCRs), we employed two ΔΔG calculation methods: thermodynamic integration (TI) with AMBER and the alchemical transfer method (AToM) with OpenMM. We calculated ΔΔG values for 53 transformations involving four class A GPCRs and evaluated the performance of AMBER-TI and AToM-OpenMM. In addition, we conducted tests using different numbers of windows and varying simulation times to achieve reliable ΔΔG results and to optimize resource utilization. Overall, both AMBER-TI and AToM-OpenMM show good agreement with the experimental data. Our results validate the applicability of AMBER-TI and AToM-OpenMM for optimization of lead compounds targeting membrane proteins.

Enhancing SILCS-MC via GPU Acceleration and Ligand Conformational Optimization with Genetic and Parallel Tempering Algorithms

In the domain of computer-aided drug design, achieving precise and accurate estimates of ligand-protein binding is paramount in the context of screening extensive drug libraries and performing ligand optimization. A fundamental aspect of the SILCS (site identification by ligand competitive saturation) methodology lies in the generation of comprehensive 3D free-energy functional group affinity maps (FragMaps), encompassing the entirety of the target molecule structure. These FragMaps offer an intricate landscape of functional group affinities across the protein, bilayer, or RNA, acting as the basis for subsequent SILCS-Monte Carlo (MC) simulations wherein ligands are docked to the target molecule. To augment the efficiency and breadth of ligand sampling capabilities, we implemented an improved SILCS-MC methodology. By harnessing the parallel computing capability of GPUs, our approach facilitates concurrent calculations over multiple ligands and binding sites, markedly enhancing the computational efficiency. Moreover, the integration of a genetic algorithm (GA) with MC allows us to employ an evolutionary approach to perform ligand sampling, assuring enhanced convergence characteristics. In addition, the potential utility of parallel tempering (PT) to improve sampling was investigated. Implementation of SILCS-MC on GPU architecture is shown to accelerate the speed of SILCS-MC calculations by over 2-orders of magnitude. Use of GA and PT yield improvements over Markov-chain MC, increasing the precision of the resultant docked orientations and binding free energies, though the extent of improvements is relatively small. Accordingly, significant improvements in speed are obtained through the GPU implementation with minor improvements in the precision of the docking obtained via the tested GA and PT algorithms.

Caught between a ROCK and a hard place: current challenges in structure-based drug design

The discipline of structure-based drug design (SBDD) is several decades old and it is tempting to think that the proliferation of experimental structures for many drug targets might make computer-aided drug design (CADD) straightforward. However, this is far from true. In this review, we illustrate some of the challenges that CADD scientists face every day in their work, even now. We use Rho-associated protein kinase (ROCK), and public domain structures and data, as an example to illustrate some of the challenges we have experienced during our project targeting this protein. We hope that this will help to prevent unrealistic expectations of what CADD can accomplish and to educate non-CADD scientists regarding the challenges still facing their CADD colleagues.

Computational Investigation of BMAA and Its Carbamate Adducts as Potential GluR2 Modulators

Beta-N-methylamino-l-alanine (BMAA) is a potential neurotoxic nonprotein amino acid, which can reach the human body through the food chain. When BMAA interacts with bicarbonate in the human body, carbamate adducts are produced, which share a high structural similarity with the neurotransmitter glutamate. It is believed that BMAA and its l-carbamate adducts bind in the glutamate binding site of ionotropic glutamate receptor 2 (GluR2). Chronic exposure to BMAA and its adducts could cause neurological illness such as neurodegenerative diseases. However, the mechanism of BMAA action and its carbamate adducts bound to GluR2 has not yet been elucidated. Here, we investigate the binding modes and the affinity of BMAA and its carbamate adducts to GluR2 in comparison to the natural agonist, glutamate, to understand whether these can act as GluR2 modulators. Initially, we perform molecular dynamics simulations of BMAA and its carbamate adducts bound to GluR2 to examine the stability of the ligands in the S1/S2 ligand-binding core of the receptor. In addition, we utilize alchemical free energy calculations to compute the difference in the free energy of binding of the beta-carbamate adduct of BMAA to GluR2 compared to that of glutamate. Our findings indicate that carbamate adducts of BMAA and glutamate remain stable in the binding site of the GluR2 compared to BMAA. Additionally, alchemical free energy results reveal that glutamate and the beta-carbamate adduct of BMAA have comparable binding affinity to the GluR2. These results provide a rationale that BMAA carbamate adducts may be, in fact, the modulators of GluR2 and not BMAA itself.

Current State of Open Source Force Fields in Protein-Ligand Binding Affinity Predictions

In drug discovery, the in silico prediction of binding affinity is one of the major means to prioritize compounds for synthesis. Alchemical relative binding free energy (RBFE) calculations based on molecular dynamics (MD) simulations are nowadays a popular approach for the accurate affinity ranking of compounds. MD simulations rely on empirical force field parameters, which strongly influence the accuracy of the predicted affinities. Here, we evaluate the ability of six different small-molecule force fields to predict experimental protein-ligand binding affinities in RBFE calculations on a set of 598 ligands and 22 protein targets. The public force fields OpenFF Parsley and Sage, GAFF, and CGenFF show comparable accuracy, while OPLS3e is significantly more accurate. However, a consensus approach using Sage, GAFF, and CGenFF leads to accuracy comparable to OPLS3e. While Parsley and Sage are performing comparably based on aggregated statistics across the whole dataset, there are differences in terms of outliers. Analysis of the force field reveals that improved parameters lead to significant improvement in the accuracy of affinity predictions on subsets of the dataset involving those parameters. Lower accuracy can not only be attributed to the force field parameters but is also dependent on input preparation and sampling convergence of the calculations. Especially large perturbations and nonconverged simulations lead to less accurate predictions. The input structures, Gromacs force field files, as well as the analysis Python notebooks are available on GitHub.

Interface-driven structural evolution on diltiazem as novel uPAR inhibitors: from in silico design to in vitro evaluation

The urokinase-type plasminogen activator receptor (uPAR) emerges as a key target for anti-metastasis owing to its pivotal role in facilitating the invasive and migratory processes of cancer cells. Recently, we identified the uPAR-targeting anti-metastatic ability of diltiazem (22), a commonly used antihypertensive agent. Fine-tuning the chemical structures of known hits represents a vital branch of drug development. To develop novel anti-metastatic drugs, we performed an interface-driven structural evolution strategy on 22. The uPAR-targeting and anti-cancer abilities of this antihypertensive drug wereidentified by us recently. Based on in silico strategy, including extensive molecular dynamics (MD) simulations, hierarchical binding free energy predictions, and ADMET profilings, we designed, synthesized, and identified three new diltiazem derivatives (221-8, 221-57, and 221-68) as uPAR inhibitors. Indeed, all of these three derivatives exhibited uPAR-depending inhibitory activity against PC-3 cell line invasion at micromolar level. Particularly, derivatives 221-68 and 221-8 showed enhanced uPAR-dependent inhibitory activity against the tumor cell invasion compared to the original compound. Microsecond timesclae MD simulations demonstrated the optimized moiety of 221-68 and 221-8 forming more comprehensive interactions with the uPAR, highlighting the reasonability of our strategy. This work introduces three novel uPAR inhibitors, which not only pave the way for the development of effective anti-metastatic therapeutics, but also emphasize the efficacy and robustness of an in silico-based lead compound optimization strategy in drug design.

An artificial biocatalytic cascade for efficient synthesis of norepinephrine by combination of engineered L-threonine transaldolase with multi-enzyme expression fine-tuning

Norepinephrine, a kind of β-adrenergic receptor agonist, is commonly used for treating shocks and hypotension caused by a variety of symptoms. The development of a straightforward, efficient and environmentally friendly biocatalytic route for manufacturing norepinephrine remains a challenge. Here, we designed and realized an artificial biocatalytic cascade to access norepinephrine starting from 3, 4-dihydroxybenzaldehyde and L-threonine mediated by a tailored-made L-threonine transaldolase PsLTTA-Mu1 and a newly screened tyrosine decarboxylase ErTDC. To overcome the imbalance of multi-enzymes in a single cell, engineering of PsLTTA for improved activity and fine-tuning expression mode of multi-enzymes in single E.coli cells were combined, leading to a robust whole cell biocatalyst ES07 that could produce 100 mM norepinephrine with 99% conversion, delivering a highest time-space yield (3.38 g/L/h) ever reported. To summarized, the current study proposed an effective biocatalytic approach for the synthesis of norepinephrine from low-cost substrates, paving the way for industrial applications of enzymatic norepinephrine production.

Rational Design of Benzobisheterocycle Metallo-β-Lactamase Inhibitors: A Tricyclic Scaffold Enhances Potency against Target Enzymes

Antimicrobial resistance is a global public health threat. Metallo-β-lactamases (MBLs) inactivate β-lactam antibiotics, including carbapenems, are disseminating among Gram-negative bacteria, and lack clinically useful inhibitors. The evolving bisthiazolidine (BTZ) scaffold inhibits all three MBL subclasses (B1-B3). We report design, synthesis, and evaluation of BTZ analogues. Structure-activity relationships identified the BTZ thiol as essential, while carboxylate is replaceable, with its removal enhancing potency by facilitating hydrophobic interactions within the MBL active site. While the introduction of a flexible aromatic ring is neutral or detrimental for inhibition, a rigid (fused) ring generated nM benzobisheterocycle (BBH) inhibitors that potentiated carbapenems against MBL-producing strains. Crystallography of BBH:MBL complexes identified hydrophobic interactions as the basis of potency toward B1 MBLs. These data underscore BTZs as versatile, potent broad-spectrum MBL inhibitors (with activity extending to enzymes refractory to other inhibitors) and provide a rational approach to further improve the tricyclic BBH scaffold.

RNF43 Inactivation Enhances the B-RAF/MEK Signaling and Creates a Combinatory Therapeutic Target in Cancer Cells

RING finger 43 (RNF43), a RING-type E3 ubiquitin ligase, is a key regulator of WNT signaling and is mutated in 6-10% of pancreatic tumors. However, RNF43-mediated effects remain unclear, as only a few in vivo substrates of RNF43 are identified. Here, it is found that RNF43-mutated pancreatic cancer cells exhibit elevated B-RAF/MEK activity and are highly sensitive to MEK inhibitors. The depletion of RNF43 in normal pancreatic ductal cells also enhances MEK activation, suggesting that it is a physiologically regulated process. It is confirmed that RNF43 ubiquitinates B-RAF at K499 to promote proteasome-dependent degradation, resulting in reduced MEK activity and proliferative ability in cancer cells. In addition, phosphorylation of B-RAF at T491 suppresses B-RAF ubiquitination by decreasing the interaction between RNF43 and B-RAF. Mutations at K499 in B-RAF are identified in various cancer types. MEK and WNT inhibitors synergistically suppress the growth of RNF43-mutated pancreatic cancer cells in vitro and in vivo. Collectively, the research reveals a novel mechanism by which RNF43 inhibits B-RAF/MEK signaling to suppress tumor growth and provide a new strategy for the treatment of RNF43-inactivated pancreatic cancer.

Dissecting the Interaction Fingerprints and Binding Affinity of BYL719 Analogs Targeting PI3Kα

Phosphatidylinositol-3-kinase Alpha (PI3Kα) is a lipid kinase which regulates signaling pathways involved in cell proliferation. Dysregulation of these pathways promotes several human cancers, pushing for the development of anticancer drugs to target PI3Kα. One such medicinal chemistry campaign at Novartis led to the discovery of BYL719 (Piqray, Alpelicib), a PI3Kα inhibitor approved by the FDA in 2019 for treatment of HR+/HER2-advanced breast cancer with a PIK3CA mutation. Structure-based drug design played a key role in compound design and optimization throughout the discovery process. However, further characterization of potency drivers via structural dynamics and energetic analyses can be advantageous for ensuing PI3Kα programs. Here, our goal is to employ various in-silico techniques, including molecular simulations and machine learning, to characterize 14 ligands from the BYL719 analogs and predict their binding affinities. The structural insights from molecular simulations suggest that although the ligand-hinge interaction is the primary driver of ligand stability at the pocket, the R group positioning at C2 or C6 of pyridine/pyrimidine also plays a major role. Binding affinities predicted via thermodynamic integration (TI) are in good agreement with previously reported IC50s. Yet, computationally demanding techniques such as TI might not always be the most efficient approach for affinity prediction, as in our case study, fast high-throughput techniques were capable of classifying compounds as active or inactive, and one docking approach showed accuracy comparable to TI.

An Efficient Approach to the Accurate Prediction of Mutational Effects in Antigen Binding to the MHC1

The major histocompatibility complex (MHC) can recognize and bind to external peptides to generate effective immune responses by presenting the peptides to T cells. Therefore, understanding the binding modes of peptide-MHC complexes (pMHC) and predicting the binding affinity of pMHCs play a crucial role in the rational design of peptide vaccines. In this study, we employed molecular dynamics (MD) simulations and free energy calculations with an Alanine Scanning with Generalized Born and Interaction Entropy (ASGBIE) method to investigate the protein-peptide interaction between HLA-A*02:01 and the G9209 peptide derived from the melanoma antigen gp100. The energy contribution of individual residue was calculated using alanine scanning, and hotspots on both the MHC and the peptides were identified. Our study shows that the pMHC binding is dominated by the van der Waals interactions. Furthermore, we optimized the ASGBIE method, achieving a Pearson correlation coefficient of 0.91 between predicted and experimental binding affinity for mutated antigens. This represents a significant improvement over the conventional MM/GBSA method, which yields a Pearson correlation coefficient of 0.22. The computational protocol developed in this study can be applied to the computational screening of antigens for the MHC1 as well as other protein-peptide binding systems.

Machine Learning Methods as a Cost-Effective Alternative to Physics-Based Binding Free Energy Calculations

The rank ordering of ligands remains one of the most attractive challenges in drug discovery. While physics-based in silico binding affinity methods dominate the field, they still have problems, which largely revolve around forcefield accuracy and sampling. Recent advances in machine learning have gained traction for protein-ligand binding affinity predictions in early drug discovery programs. In this article, we perform retrospective binding free energy evaluations for 172 compounds from our internal collection spread over four different protein targets and five congeneric ligand series. We compared multiple state-of-the-art free energy methods ranging from physics-based methods with different levels of complexity and conformational sampling to state-of-the-art machine-learning-based methods that were available to us. Overall, we found that physics-based methods behaved particularly well when the ligand perturbations were made in the solvation region, and they did not perform as well when accounting for large conformational changes in protein active sites. On the other end, machine-learning-based methods offer a good cost-effective alternative for binding free energy calculations, but the accuracy of their predictions is highly dependent on the experimental data available for training the model.

Comprehensive, Open-Source, and Automated Workflow for Multisite λ-Dynamics in Lead Optimization

Multisite λ-dynamics (MSLD) is a highly efficient binding free energy calculation method that samples multiple ligands in a single round by assigning different λ values to the alchemical part of each ligand. This method holds great promise for lead optimization (LO) in drug discovery. However, the complex data preparation and simulation process limits its widespread application in diverse protein-ligand systems. To address this challenge, we developed a comprehensive, open-source, and automated workflow for MSLD calculations based on the BLaDE dynamics engine. This workflow incorporates the Ligand Internal and Cartesian coordinate reconstruction-based alignment algorithm (LIC-align) and an optimized maximum common substructure (MCS) search algorithm to accurately generate MSLD multiple topologies with ideal perturbation patterns. Furthermore, our workflow is highly modularized, allowing straightforward integration and extension of various simulation techniques, and is highly accessible to nonexperts. This workflow was validated by calculating the relative binding free energies of large-scale congeneric ligands, many of which have large perturbing groups. The agreement between the calculations and experiments was excellent, with an average unsigned error of 1.08 ± 0.47 kcal/mol. More than 57.1% of the ligands had an error of less than 1.0 kcal/mol, and the perturbations of 6 targets were fully connected via the calculations, while those of 2 targets were connected via both calculations and experimental data. The Pearson correlation coefficient reached 0.88, indicating that the MSLD workflow provides accurate predictions that can guide lead optimization in drug discovery. We also examined the impact of single-site versus multisite perturbations, ligand grouping by perturbing group size, and the position of the anchor atom on the MSLD performance. By integrating our proposed LIC-align and optimized MCS search algorithm along with the coping strategies to handle challenging molecular substructures, our workflow can handle many realistic scenarios more reasonably than all previously published methods. Moreover, we observed that our MSLD workflow achieved similar accuracy to free energy perturbation (FEP) while improving computational efficiency by over 1 order of magnitude in speedup. These findings provide valuable insights and strategies for further MSLD development, making MSLD a competitive tool for lead optimization.

Accurate calculation of affinity changes to the close state of influenza A M2 transmembrane domain in response to subtle structural changes of adamantyl amines using free energy perturbation methods in different lipid bilayers

Experimental binding free energies of 27 adamantyl amines against the influenza M2(22-46) WT tetramer, in its closed form at pH 8, were measured by ITC in DPC micelles. The measured Kd's range is ~44 while the antiviral potencies (IC50) range is ~750 with a good correlation between binding free energies computed with Kd and IC50 values (r = 0.76). We explored with MD simulations (ff19sb, CHARMM36m) the binding profile of complexes with strong, moderate and weak binders embedded in DMPC, DPPC, POPC or a viral mimetic membrane and using different experimental starting structures of M2. To predict accurately differences in binding free energy in response to subtle changes in the structure of the ligands, we performed 18 alchemical perturbative single topology FEP/MD NPT simulations (OPLS2005) using the BAR estimator (Desmond software) and 20 dual topology calculations TI/MD NVT simulations (ff19sb) using the MBAR estimator (Amber software) for adamantyl amines in complex with M2(22-46) WT in DMPC, DPPC, POPC. We observed that both methods with all lipids show a very good correlation between the experimental and calculated relative binding free energies (r = 0.77-0.87, mue = 0.36-0.92 kcal mol-1) with the highest performance achieved with TI/MBAR and lowest performance with FEP/BAR in DMPC bilayers. When antiviral potencies are used instead of the Kd values for computing the experimental binding free energies we obtained also good performance with both FEP/BAR (r = 0.83, mue = 0.75 kcal mol-1) and TI/MBAR (r = 0.69, mue = 0.77 kcal mol-1).

Computational Workflow for Refining AlphaFold Models in Drug Design Using Kinetic and Thermodynamic Binding Calculations: A Case Study for the Unresolved Inactive Human Adenosine A3 Receptor

A structure-based drug design pipeline that considers both thermodynamic and kinetic binding data of ligands against a receptor will enable the computational design of improved drug molecules. For unresolved GPCR-ligand complexes, a workflow that can apply both thermodynamic and kinetic binding data in combination with alpha-fold (AF)-derived or other homology models and experimentally resolved binding modes of relevant ligands in GPCR-homologs needs to be tested. Here, as test case, we studied a congeneric set of ligands that bind to a structurally unresolved G protein-coupled receptor (GPCR), the inactive human adenosine A3 receptor (hA3R). We tested three available homology models from which two have been generated from experimental structures of hA1R or hA2AR and one model was a multistate alphafold 2 (AF2)-derived model. We applied alchemical calculations with thermodynamic integration coupled with molecular dynamics (TI/MD) simulations to calculate the experimental relative binding free energies and residence time (τ)-random accelerated MD (τ-RAMD) simulations to calculate the relative residence times (RTs) for antagonists. While the TI/MD calculations produced, for the three homology models, good Pearson correlation coefficients, correspondingly, r = 0.74, 0.62, and 0.67 and mean unsigned error (mue) values of 0.94, 1.31, and 0.81 kcal mol-1, the τ-RAMD method showed r = 0.92 and 0.52 for the first two models but failed to produce accurate results for the multistate AF2-derived model. With subsequent optimization of the AF2-derived model by reorientation of the side chain of R1735.34 located in the extracellular loop 2 (EL2) that blocked ligand's unbinding, the computational model showed r = 0.84 for kinetic data and improved performance for thermodynamic data (r = 0.81, mue = 0.56 kcal mol-1). Overall, after refining the multistate AF2 model with physics-based tools, we were able to show a strong correlation between predicted and experimental ligand relative residence times and affinities, achieving a level of accuracy comparable to an experimental structure. The computational workflow used can be applied to other receptors, helping to rank candidate drugs in a congeneric series and enabling the prioritization of leads with stronger binding affinities and longer residence times.

Discovery of Novel Inhibitors of BRD4 for Treating Prostate Cancer: A Comprehensive Case Study for Considering Water Networks in Virtual Screening and Drug Design

Androgen receptor (AR) is the primary target for treating prostate cancer (PCa), which inevitably progresses due to drug-resistant mutations. Bromodomain-containing protein 4 (BRD4) has been a new potential drug target for PCa treatment. Herein, we report the rational design and discovery of novel BRD4 inhibitors through computer-aided drug design (CADD), and a hit compound SQ-1 (IC50 = 676 nM) was identified by structure-based virtual screening (SBVS) with the conserved water network. To optimize the structure of SQ-1, the free energy landscape was constructed, and the binding mechanism was explored by characterizing the water profile and the dissociation mechanism. Finally, the compound SQ-17 with improved inhibitory activity (IC50 < 100 nM) was discovered, which showed potent antiproliferative activity against LNCaP. These data highlighted a successful attempt to identify and optimize a small molecule by comprehensive CADD application and provided essential clues for developing novel therapeutics for PCa treatment.

Transferable Machine Learning Interatomic Potential for Bond Dissociation Energy Prediction of Drug-like Molecules

We present a transferable MACE interatomic potential that is applicable to open- and closed-shell drug-like molecules containing hydrogen, carbon, and oxygen atoms. Including an accurate description of radical species extends the scope of possible applications to bond dissociation energy (BDE) prediction, for example, in the context of cytochrome P450 (CYP) metabolism. The transferability of the MACE potential was validated on the COMP6 data set, containing only closed-shell molecules, where it reaches better accuracy than the readily available general ANI-2x potential. MACE achieves similar accuracy on two CYP metabolism-specific data sets, which include open- and closed-shell structures. This model enables us to calculate the aliphatic C-H BDE, which allows us to compare reaction energies of hydrogen abstraction, which is the rate-limiting step of the aliphatic hydroxylation reaction catalyzed by CYPs. On the "CYP 3A4" data set, MACE achieves a BDE RMSE of 1.37 kcal/mol and better prediction of BDE ranks than alternatives: the semiempirical AM1 and GFN2-xTB methods and the ALFABET model that directly predicts bond dissociation enthalpies. Finally, we highlight the smoothness of the MACE potential over paths of sp3C-H bond elongation and show that a minimal extension is enough for the MACE model to start finding reasonable minimum energy paths of methoxy radical-mediated hydrogen abstraction. Altogether, this work lays the ground for further extensions of scope in terms of chemical elements, (CYP-mediated) reaction classes and modeling the full reaction paths, not only BDEs.

Binding of Cholesterol to the N-Terminal Domain of the NPC1L1 Transporter: Analysis of the Epimerization-Related Binding Selectivity and Loop Mutations

Cholesterol is a fat-like substance with a pivotal physiological relevance in humans, and its homeostasis is tightly regulated by various cellular processes, including the import in the small intestine and the reabsorption in the biliary ducts by the Niemann-Pick C1 Like 1 (NPC1L1) importer. NPC1L1 can mediate the absorption of a variety of sterols but strikingly exhibits a large sensitivity to cholesterol epimerization. This study examines the molecular basis of the epimerization-related selective binding of cholesterol by combining extended unbiased molecular dynamics simulations of the apo and holo species of the N-terminal domain of wild-type NPC1L1, in conjunction with relative binding free energy, umbrella sampling, and well-tempered metadynamics calculations. The analysis of the results discloses the existence of two distinct binding modes for cholesterol and epi-cholesterol. The former binds deeper in the cavity, forming key hydrogen-bond interactions with Q95, S56, and a water molecule. In contrast, epi-cholesterol is shifted ca. 3 Å to the mouth of the cavity and the transition to the Q95 site is prevented by an energetic barrier of 4.1 kcal·mol-1. Thus, the configuration of the hydroxyl group of cholesterol, together with the presence of a structural water molecule, is a key feature for effective absorption. Finally, whereas these findings may seemingly be challenged by single-point mutations that impair cholesterol transport but have a mild impact on the binding of cholesterol to the Q95 binding site, our results reveal that they have a drastic influence on the conformational landscape of the α8/β7 loop in the apo species compared to the wild-type protein. Overall, the results give support to the functional role played by the α8/β7 loop in regulating the access of ligands to NPC1L1, and hence to interpreting the impact of these mutations on diseases related to disruption of sterol absorption, paving the way to understanding certain physiological dysfunctions.

Molecular Dynamics as a Tool for Virtual Ligand Screening

Rational drug design is essential for new drugs to emerge, especially when the structure of a target protein or nucleic acid is known. To that purpose, high-throughput virtual ligand screening campaigns aim at discovering computationally new binding molecules or fragments to modulate particular biomolecular interactions or biological activities, related to a disease process. The structure-based virtual ligand screening process primarily relies on docking methods which allow predicting the binding of a molecule to a biological target structure with a correct conformation and the best possible affinity. The docking method itself is not sufficient as it suffers from several and crucial limitations (lack of full protein flexibility information, no solvation and ion effects, poor scoring functions, and unreliable molecular affinity estimation).At the interface of computer techniques and drug discovery, molecular dynamics (MD) allows introducing protein flexibility before or after a docking protocol, refining the structure of protein-drug complexes in the presence of water, ions, and even in membrane-like environments, describing more precisely the temporal evolution of the biological complex and ranking these complexes with more accurate binding energy calculations. In this chapter, we describe the up-to-date MD, which plays the role of supporting tools in the virtual ligand screening (VS) process.Without a doubt, using docking in combination with MD is an attractive approach in structure-based drug discovery protocols nowadays. It has proved its efficiency through many examples in the literature and is a powerful method to significantly reduce the amount of required wet experimentations (Tarcsay et al, J Chem Inf Model 53:2990-2999, 2013; Barakat et al, PLoS One 7:e51329, 2012; De Vivo et al, J Med Chem 59:4035-4061, 2016; Durrant, McCammon, BMC Biol 9:71-79, 2011; Galeazzi, Curr Comput Aided Drug Des 5:225-240, 2009; Hospital et al, Adv Appl Bioinforma Chem 8:37-47, 2015; Jiang et al, Molecules 20:12769-12786, 2015; Kundu et al, J Mol Graph Model 61:160-174, 2015; Mirza et al, J Mol Graph Model 66:99-107, 2016; Moroy et al, Future Med Chem 7:2317-2331, 2015; Naresh et al, J Mol Graph Model 61:272-280, 2015; Nichols et al, J Chem Inf Model 51:1439-1446, 2011; Nichols et al, Methods Mol Biol 819:93-103, 2012; Okimoto et al, PLoS Comput Biol 5:e1000528, 2009; Rodriguez-Bussey et al, Biopolymers 105:35-42, 2016; Sliwoski et al, Pharmacol Rev 66:334-395, 2014).

Dealing with Induced Fit, Conformational Selection, and Secondary Poses in Molecular Dynamics Simulations for Reliable Free Energy Predictions

In this study, we have tested the performance of standard molecular dynamics (MD) simulations, replicates of shorter standard MD simulations, and Hamiltonian Replica Exchange (HREM) simulations for the sampling of two macrocyclic hosts for guest delivery, characterized by induced fit (phenyl-based host) and conformation selection (naphthyl-based host) and of the ODR-BRD4(I) drug-receptor system where the ligand can assume two main poses. For the optimization of the HREM simulation, we have proposed and tested an on-the-fly iterative scheme for equalizing the acceptance ratio along the replica progression at a constant replica number resulting in a moderate impact of the sampling efficiency. Concerning standard MD, we have found that, while splitting the total allocated simulation time in short MD replicates can reproduce the sampling efficiency of HREM in the phenyl-based host and in the ODR-BRD4(I) complex, in the naphthyl-based macrocycle, characterized by long-lived metastable states, enhanced sampling techniques are the only viable alternative for a reliable canonical sampling of the rugged conformational landscape.

Alchemical Calculation of Relative Free Energies for Charge-Changing Mutations at Protein-Protein Interfaces Considering Fixed and Variable Protonation States

The calculation of relative free energies (ΔΔG) for charge-changing mutations at protein-protein interfaces through alchemical methods remains challenging due to variations in the system's net charge during charging steps, the possibility of mutated and contacting ionizable residues occurring in various protonation states, and undersampling issues. In this study, we present a set of strategies, collectively termed TIRST/TIRST-H+, to address some of these challenges. Our approaches combine thermodynamic integration (TI) with the prediction of pKa shifts to calculate ΔΔG values. Moreover, special sets of restraints are employed to keep the alchemically transformed molecules separated. The accuracy of the devised approaches was assessed on a large and diverse data set comprising 164 point mutations of charged residues (Asp, Glu, Lys, and Arg) to Ala at the protein-protein interfaces of complexes with known three-dimensional structures. Mean absolute and root-mean-square errors ranging from 1.38 to 1.66 and 1.89 to 2.44 kcal/mol, respectively, and Pearson correlation coefficients of ∼0.6 were obtained when testing the approaches on the selected data set using the GPU-TI module of Amber18 suite and the ff14SB force field. Furthermore, the inclusion of variable protonation states for the mutated acid residues improved the accuracy of the predicted ΔΔG values. Therefore, our results validate the use of TIRST/TIRST-H+ in prospective studies aimed at evaluating the impact of charge-changing mutations to Ala on the stability of protein-protein complexes.

Understanding base and backbone contributions of phosphorothioate DNA for molecular recognition with SBD proteins

Bacterial DNA phosphorothioate (PT) modification provides a specific anchoring site for sulfur-binding proteins (SBDs). Besides, their recognition patterns include phosphate links and bases neighboring the PT-modified site, thereby bringing about genome sequence-dependent properties in PT-related epigenetics. Here, we analyze the contributions of the DNA backbone (phosphates and deoxyribose) and bases bound with two SBD proteins in Streptomyces pristinaespiralis and coelicolor (SBDSco and SBDSpr). The chalcogen-hydrophobic interactions remained constantly at the anchoring site while the adjacent bases formed conditional and distinctive non-covalent interactions. More importantly, SBD/PT-DNA interactions were not limited within the traditional "4-bp core" range from 5'-I to 3'-III but extended to upstream 5'-II and 5'-III bases and even 5''-I to 5''-III at the non-PT-modified complementary strand. From the epigenetic viewpoint, bases 3'-II, 5''-I, and 5''-III of SBDSpr and 3'-II, 5''-II, and 5''-III of SBDSco present remarkable differentiations in the molecular recognitions. From the protein viewpoint, H102 in SBDSpr and R191 in SBDSco contribute significantly while proline residues at the PT-bound site are strictly conserved for the PT-chalcogen bond. The mutual and make-up mutations are proposed to alter the SBD/PT-DNA recognition pattern, besides additional chiral phosphorothioate modifications on phosphates 5'-II, 5'-II, 3'-I, and 3'-II.

Effects of astragaloside IV on glucocorticoid-induced avascular necrosis of the femoral head via regulating Akt-related pathways

We investigated the role of astragaloside IV (AS-IV) in preventing glucocorticoid-induced avascular necrosis of the femoral head (ANFH) and the underlying molecular mechanisms. Network pharmacology was used to predict the molecular targets of AS-IV. Molecular dynamic simulations were performed to explore the binding mechanism and interaction mode between AS-IV and Akt. Rat models of glucocorticoid-induced ANFH with AS-IV intervention were established, and osteogenesis, angiogenesis, apoptosis and oxidative stress were evaluated before and after blocking the PI3K/Akt pathway with LY294002. The effects of glucocorticoid and AS-IV on bone marrow mesenchymal stem cells and human umbilical vein endothelial cells incubated with and without LY294002 were determined. Downregulated p-Akt expression could be detected in the femoral heads of glucocorticoid-induced ANFH patients and rats. AS-IV increased trabecular bone integrity and vessel density of the femoral head in the model rats. AS-IV increased Akt phosphorylation and upregulated osteogenesis-, angiogenesis-, apoptosis- and oxidative stress-related proteins and mRNA and downregulated Bax, cleaved caspase-3 and cytochrome c levels. AS-IV promoted human umbilical vein endothelial cell migration, proliferation and tube formation ability; bone marrow mesenchymal stem cell proliferation; and osteogenic differentiation under glucocorticoid influence. AS-IV inhibited apoptosis. LY294002 inhibited these effects. AS-IV prevented glucocorticoid-induced ANFH by promoting osteogenesis and angiogenesis via the Akt/Runx2 and Akt/HIF-1α/VEGF pathways, respectively, and suppressing apoptosis and oxidative stress via the Akt/Bad/Bcl-2 and Akt/Nrf2/HO-1 pathways, respectively.

Rutin from Begonia roxburghii modulates iNOS and Sep A activity in treatment of Shigella flexneri induced diarrhoea in rats: An in vitro, in vivo and computational analysis

In developing countries, diarrhoea is a major issue of concern, where consistent use of antibiotics has resulted in several side effects along with development of resistance among pathogens against these antibiotics. Since natural products are becoming the treatment of choice, therefore present investigation involves mechanistic evaluation of antidiarrhoeal potential of Begonia roxburghii and its marker rutin against Shigella flexneri (SF) induced diarrhoea in rats following in vitro, in vivo and in silico protocols. The roots of the plant are used as vegetable in the North East India and are also used traditionally in treating diarrhoea. Phytochemically standardized ethanolic extract of B. roxburghii (EBR) roots and its marker rutin were first subjected to in vitro antibacterial evaluation against SF. Diarrhoea was induced in rats using suspension of SF and various diarrhoeagenic parameters were examined after first, third and fifth day of treatment at 100, 200 and 300 mg/kg, p.o. with EBR and 50 mg/kg, p.o. with rutin respectively. Additionally, density of SF in stools, stool water content, haematological and biochemical parameters, cytokine profiling, ion concentration, histopathology and Na+/K+-ATPase activity were also performed. Molecular docking and dynamics simulation studies of ligand rutin was studied against secreted extracellular protein A (Sep A, PDB: 5J44) from SF and Inducible nitric oxide synthase (iNOS, PDB: 1DD7) followed by network pharmacology. EBR and rutin demonstrated a potent antibacterial activity against SF and also showed significant recovery from diarrhoea (EBR: 81.29 ± 0.91% and rutin: 75.27 ± 0.89%) in rats after five days of treatment. EBR and rutin also showed significant decline in SF density in stools, decreased cytokine expression, potential antioxidant activity, cellular proliferative nature and recovered ion loss due to enhanced Na+/K+-ATPase activity, which was also supported by histopathology. Rutin showed a very high docking score of -11.61 and -9.98 kcal/mol against iNOS and Sep A respectively and their stable complex was also confirmed through dynamics, while network pharmacology suggested that, rutin is quite capable of modulating the pathways of iNOS and Sep A. Thus, we may presume that rutin played a key role in the observed antidiarrhoeal activity of B. roxburghii against SF induced diarrhoea.

The maximal and current accuracy of rigorous protein-ligand binding free energy calculations

Computational techniques can speed up the identification of hits and accelerate the development of candidate molecules for drug discovery. Among techniques for predicting relative binding affinities, the most consistently accurate is free energy perturbation (FEP), a class of rigorous physics-based methods. However, uncertainty remains about how accurate FEP is and can ever be. Here, we present what we believe to be the largest publicly available dataset of proteins and congeneric series of small molecules, and assess the accuracy of the leading FEP workflow. To ascertain the limit of achievable accuracy, we also survey the reproducibility of experimental relative affinity measurements. We find a wide variability in experimental accuracy and a correspondence between binding and functional assays. When careful preparation of protein and ligand structures is undertaken, FEP can achieve accuracy comparable to experimental reproducibility. Throughout, we highlight reliable protocols that can help maximize the accuracy of FEP in prospective studies.

Virtual screening and biological activity evaluation of novel efflux pump inhibitors targeting AdeB

The AdeABC efflux pump is an important mechanism causing multidrug resistance in Acinetobacter baumannii, and its main component AdeB can recognize carbapenems, aminoglycosides, and other multi-class antibiotics and efflux them intracellularly, which is an ideal target for the development of anti-multidrug resistant bacteria drugs. Here, we combined multiple computer-aided drug design methods to target AdeB to identify promising novel structural inhibitors. Virtual screening was performed by molecular docking and molecular dynamics simulation (MD) and 12 potential compounds were identified from the databases. Meanwhile, their biological activities were validated by in vitro activity assays, and ChemDiv L676-2179 (γ-IFN), ChemDiv L676-1461, and Chembridge 53717615 were confirmed to suppress efflux effects and restore antibiotic susceptibility of resistant bacteria, which are expected to be developed as adjuvant drugs for the treatment of multi-drug resistant Acinetobacter baumannii clinical infections.

TLK1-mediated RAD54 phosphorylation spatio-temporally regulates Homologous Recombination Repair

Environmental agents like ionizing radiation (IR) and chemotherapeutic drugs can cause severe damage to the DNA, often in the form of double-strand breaks (DSBs). Remaining unrepaired, DSBs can lead to chromosomal rearrangements, and cell death. One major error-free pathway to repair DSBs is homologous recombination repair (HRR). Tousled-like kinase 1 (TLK1), a Ser/Thr kinase that regulates the DNA damage checkpoint, has been found to interact with RAD54, a central DNA translocase in HRR. To determine how TLK1 regulates RAD54, we inhibited or depleted TLK1 and tested how this impacts HRR in human cells using a ISce-I-GR-DsRed fused reporter endonuclease. Our results show that TLK1 phosphorylates RAD54 at three threonines (T41, T59 and T700), two of which are located within its N-terminal domain (NTD) and one is located within its C-terminal domain (CTD). Phosphorylation at both T41 and T59 supports HRR and protects cells from DNA DSB damage. In contrast, phosphorylation of T700 leads to impaired HRR and engenders no protection to cells from cytotoxicity and rather results in repair delay. Further, our work enlightens the effect of RAD54-T700 (RAD54-CTD) phosphorylation by TLK1 in mammalian system and reveals a new site of interaction with RAD51.

The Coming of Age of AI/ML in Drug Discovery, Development, Clinical Testing, and Manufacturing: The FDA Perspectives

Artificial intelligence (AI) and machine learning (ML) represent significant advancements in computing, building on technologies that humanity has developed over millions of years-from the abacus to quantum computers. These tools have reached a pivotal moment in their development. In 2021 alone, the U.S. Food and Drug Administration (FDA) received over 100 product registration submissions that heavily relied on AI/ML for applications such as monitoring and improving human performance in compiling dossiers. To ensure the safe and effective use of AI/ML in drug discovery and manufacturing, the FDA and numerous other U.S. federal agencies have issued continuously updated, stringent guidelines. Intriguingly, these guidelines are often generated or updated with the aid of AI/ML tools themselves. The overarching goal is to expedite drug discovery, enhance the safety profiles of existing drugs, introduce novel treatment modalities, and improve manufacturing compliance and robustness. Recent FDA publications offer an encouraging outlook on the potential of these tools, emphasizing the need for their careful deployment. This has expanded market opportunities for retraining personnel handling these technologies and enabled innovative applications in emerging therapies such as gene editing, CRISPR-Cas9, CAR-T cells, mRNA-based treatments, and personalized medicine. In summary, the maturation of AI/ML technologies is a testament to human ingenuity. Far from being autonomous entities, these are tools created by and for humans designed to solve complex problems now and in the future. This paper aims to present the status of these technologies, along with examples of their present and future applications.

Promiscuous recognition of MR1 drives self-reactive mucosal-associated invariant T cell responses

Mucosal-associated invariant T (MAIT) cells use canonical semi-invariant T cell receptors (TCR) to recognize microbial riboflavin precursors displayed by the antigen-presenting molecule MR1. The extent of MAIT TCR crossreactivity toward physiological, microbially unrelated antigens remains underexplored. We describe MAIT TCRs endowed with MR1-dependent reactivity to tumor and healthy cells in the absence of microbial metabolites. MAIT cells bearing TCRs crossreactive toward self are rare but commonly found within healthy donors and display T-helper-like functions in vitro. Experiments with MR1-tetramers loaded with distinct ligands revealed significant crossreactivity among MAIT TCRs both ex vivo and upon in vitro expansion. A canonical MAIT TCR was selected on the basis of extremely promiscuous MR1 recognition. Structural and molecular dynamic analyses associated promiscuity to unique TCRβ-chain features that were enriched within self-reactive MAIT cells of healthy individuals. Thus, self-reactive recognition of MR1 represents a functionally relevant indication of MAIT TCR crossreactivity, suggesting a potentially broader role of MAIT cells in immune homeostasis and diseases, beyond microbial immunosurveillance.

Analysing the impact of eriosematin E from Eriosema chinense Vogel. against different diarrhoeagenic pathovars of Escherichia coli using in silico and in vitro approach

Since diarrhoea is reportedly the third largest cause of fatality among kids, therefore it is considered to be one of the major areas of concerns among developing nations. The main causative agents of diarrhoea include Escherichia coli, Vibrio cholera, and Shigella spp where E. coli shares the maximum contribution. The roots of the plant Eriosema chinense Vogel. (Fabaceae) are traditionally used by the native tribes of Meghalaya, India to treat diarrhoea. From previous reports, the plant and its marker eriosematin E have been reported to have antidiarrhoeal potential against pathogenic and nonpathogenic diarrhoea. Therefore, the objective of the current investigation was to use in silico studies to determine the efficacy of eriosematin E against different diarrhoeagenic strains of E. coli. Six different pathovars of E. coli i.e. enteropathogenic E. coli (EPEC), enterotoxigenic E. coli (ETEC), enterohaemorrhagic E. coli (EHEC), enteroaggregative E. coli (EAEC), uropathogenic E. coli (UPEC) and enteroinvasive E. coli (EIEC) were subjected to docking simulation studies utilizing Glide module of Schrodinger Maestro 2018-1 MM Share Version. Based on the obtained binding energy and balance between H-bonding, hydrophobic, and salt bridge interactions eriosematin E was found to be most effective against EPEC followed by EAEC and ETEC, while UPEC and EHEC were moderately affected. The molecular dynamics studies suggested a higher affinity of eriosematin E towards heat-labile enterotoxin b-pentamer from ETEC. The in vitro antibacterial studies against the universal strain S. aureus 12981 and E. coli 10418 revealed the effectiveness of eriosematin E showing MIC values of ≥256 µg/mL.Communicated by Ramaswamy H. Sarma.

Sustainable Electropolymerization of Zingerone and Its C2 Symmetric Dimer for Amperometric Biosensor Films

Polymeric permselective films are frequently used for amperometric biosensors to prevent electroactive interference present in the target matrix. Phenylenediamines are the most commonly used for the deposition of shielding polymeric films against interfering species; however, even phenolic monomers have been utilized in the creation of these films for microsensors and biosensors. The purpose of this paper is to evaluate the performances of electrosynthesized polymers, layered by means of constant potential amperometry (CPA), of naturally occurring compound zingerone (ZING) and its dimer dehydrozingerone (ZING DIM), which was obtained by straight oxidative coupling reaction. The polymers showed interesting shielding characteristics against the main interfering species, such as ascorbic acid (AA): actually, polyZING exhibited an AA shielding aptitude comprised between 77.6 and 99.6%, comparable to that obtained with PPD. Moreover, a marked capability of increased monitoring of hydrogen peroxide (HP), when data were compared with bare metal results, was observed. In particular, polyZING showed increases ranging between 55.6 and 85.6%. In the present work, the molecular structures of the obtained polymers have been theorized and docking analyses were performed to understand their peculiar characteristics better. The structures were docked using the Lamarckian genetic algorithm (LGA). Glutamate biosensors based on those polymers were built, and their performances were compared with biosensors based on PPD, which is the most widespread polymer for the construction of amperometric biosensors.

NE-RDFE: A protocol and toolkit for computing relative dissociation free energies with GROMACS between dissimilar molecules using bidirectional nonequilibrium dual topology schemes

We describe a step-by-step protocol and toolkit for the computation of the relative dissociation free energy (RDFE) with the GROMACS molecular dynamics package, based on a novel bidirectional nonequilibrium alchemical approach. The proposed methodology does not require any intervention on the code and allows computing with good accuracy the RDFE between small molecules with arbitrary differences in volume, charge, and chemical topology. The procedure is illustrated for the challenging SAMPL9 batch of host-guest pairs. The article is supplemented by a detailed online tutorial, available at https://procacci.github.io/vdssb_gromacs/NE-RDFE and by a public Zenodo repository available at https://zenodo.org/record/6982932.

Study of SQ109 analogs binding to mycobacterium MmpL3 transporter using MD simulations and alchemical relative binding free energy calculations

N-geranyl-N΄-(2-adamantyl)ethane-1,2-diamine (SQ109) is a tuberculosis drug that has high potency against Mycobacterium tuberculosis (Mtb) and may function by blocking cell wall biosynthesis. After the crystal structure of MmpL3 from Mycobacterium smegmatis in complex with SQ109 became available, it was suggested that SQ109 inhibits Mmpl3 mycolic acid transporter. Here, we showed using molecular dynamics (MD) simulations that the binding profile of nine SQ109 analogs with inhibitory potency against Mtb and alkyl or aryl adducts at C-2 or C-1 adamantyl carbon to MmpL3 was consistent with the X-ray structure of MmpL3 - SQ109 complex. We showed that rotation of SQ109 around carbon-carbon bond in the monoprotonated ethylenediamine unit favors two gauche conformations as minima in water and lipophilic solvent using DFT calculations as well as inside the transporter's binding area using MD simulations. The binding assays in micelles suggested that the binding affinity of the SQ109 analogs was increased for the larger, more hydrophobic adducts, which was consistent with our results from MD simulations of the SQ109 analogues suggesting that sizeable C-2 adamantyl adducts of SQ109 can fill a lipophilic region between Y257, Y646, F260 and F649 in MmpL3. This was confirmed quantitatively by our calculations of the relative binding free energies using the thermodynamic integration coupled with MD simulations method with a mean assigned error of 0.74 kcal mol^-1 compared to the experimental values.

Automated relative binding free energy calculations from SMILES to ΔΔG

In drug discovery, computational methods are a key part of making informed design decisions and prioritising experiments. In particular, optimizing compound affinity is a central concern during the early stages of development. In the last 10 years, alchemical free energy (FE) calculations have transformed our ability to incorporate accurate in silico potency predictions in design decisions, and represent the 'gold standard' for augmenting experiment-driven drug discovery. However, relative FE calculations are complex to set up, require significant expert intervention to prepare the calculation and analyse the results or are provided only as closed-source software, not allowing for fine-grained control over the underlying settings. In this work, we introduce an end-to-end relative FE workflow based on the non-equilibrium switching approach that facilitates calculation of binding free energies starting from SMILES strings. The workflow is implemented using fully modular steps, allowing various components to be exchanged depending on licence availability. We further investigate the dependence of the calculated free energy accuracy on the initial ligand pose generated by various docking algorithms. We show that both commercial and open-source docking engines can be used to generate poses that lead to good correlation of free energies with experimental reference data.

MELD-Bracket Ranks Binding Affinities of Diverse Sets of Ligands

Affinity ranking of structurally diverse small-molecule ligands is a challenging problem with important applications in structure-based drug discovery. Absolute binding free energy methods can model diverse ligands, but the high computational cost of the current methods limits application to data sets with few ligands. We recently developed MELD-Bracket, a Molecular Dynamics method for efficient affinity ranking of ligands [ JCTC 2022, 18 (1), 374-379]. It utilizes a Bayesian framework to guide sampling to relevant regions of phase space, and it couples this with a bracket-like competition on a pool of ligands. Here we find that 6-competitor MELD-Bracket can rank dozens of diverse ligands that have low structural similarity and different net charges. We benchmark it on four protein systems─PTB1B, Tyk2, BACE, and JAK3─having varied modes of interactions. We also validated 8-competitor and 12-competitor protocols. The MELD-Bracket protocols presented here may have the appropriate balance of accuracy and computational efficiency to be suitable for ranking diverse ligands from typical drug discovery campaigns.

System truncation accelerates binding affinity calculations with the fragment molecular orbital method: A benchmark study

The fragment molecular orbital (FMO) method is a fast quantum-mechanical method that divides systems into pieces of fragments and performs ab initio calculations. The system truncation enables further speed improvement. In this article, we systematically study the effects of system truncations on binding affinity calculations obtained with FMO in combination with either the polarizable continuum model (FMO/PCM) or in combination with the Møller-Plesset method (FMO-MP2). We have used five protein complexes with ligands of several charged states. The calculated binding energies of the size variants of the truncated system, including only a restricted number of atoms around the ligand, are compared to the energy obtained from a full system. The result shows that the systems could be truncated to a radius of 8 Å from neutral ligands within an error of 0.7 kcal/mol, and 12 Å from charged ligands within an error of 1.1 kcal/mol for calculating the binding energy in solution.

Redox sensitive human mitochondrial aconitase and its interaction with frataxin: In vitro and in silico studies confirm that it takes two to tango

Mitochondrial aconitase (ACO2) has been postulated as a redox sensor in the tricarboxylic acid cycle. Its high sensitivity towards reactive oxygen and nitrogen species is due to its particularly labile [4Fe-4S]²⁺ prosthetic group which yields an inactive [3Fe-4S]⁺ cluster upon oxidation. Moreover, ACO2 was found as a main oxidant target during aging and in pathologies where mitochondrial dysfunction is implied. Herein, we report the expression and characterization of recombinant human ACO2 and its interaction with frataxin (FXN), a protein that participates in the de novo biosynthesis of Fe-S clusters. A high yield of pure ACO2 (≥99%, 22 ± 2 U/mg) was obtained and kinetic parameters for citrate, isocitrate, and cis-aconitate were determined. Superoxide, carbonate radical, peroxynitrite, and hydrogen peroxide reacted with ACO2 with second-order rate constants of 10⁸, 10⁸, 10⁵, and 10² M^-1 s^-1, respectively. Temperature-induced unfolding assessed by tryptophan fluorescence of ACO2 resulted in apparent melting temperatures of 51.1 ± 0.5 and 43.6 ± 0.2 °C for [4Fe-4S]²⁺ and [3Fe-4S]⁺ states of ACO2, sustaining lower thermal stability upon cluster oxidation. Differences in protein dynamics produced by the Fe-S cluster redox state were addressed by molecular dynamics simulations. Reactivation of [3Fe-4S]⁺-ACO2 by FXN was verified by activation assays and direct iron-dependent interaction was confirmed by protein-protein interaction ELISA and fluorescence spectroscopic assays. Multimer modeling and protein-protein docking predicted an ACO2-FXN complex where the metal ion binding region of FXN approaches the [3Fe-4S]⁺ cluster, supporting that FXN is a partner for reactivation of ACO2 upon oxidative cluster inactivation.

Development and test of highly accurate endpoint free energy methods. 1: Evaluation of ABCG2 charge model on solvation free energy prediction and optimization of atom radii suitable for more accurate solvation free energy prediction by the PBSA method

Accurate estimation of solvation free energy (SFE) lays the foundation for accurate prediction of binding free energy. The Poisson-Boltzmann (PB) or generalized Born (GB) combined with surface area (SA) continuum solvation method (PBSA and GBSA) have been widely used in SFE calculations because they can achieve good balance between accuracy and efficiency. However, the accuracy of these methods can be affected by several factors such as the charge models, polar and nonpolar SFE calculation methods and the atom radii used in the calculation. In this work, the performance of the ABCG2 (AM1-BCC-GAFF2) charge model as well as other two charge models, that is, RESP (Restrained Electrostatic Potential) and AM1-BCC (Austin Model 1-bond charge corrections), on the SFE prediction of 544 small molecules in water by PBSA/GBSA was evaluated. In order to improve the performance of the PBSA prediction based on the ABCG2 charge, we further explored the influence of atom radii on the prediction accuracy and yielded a set of atom radius parameters for more accurate SFE prediction using PBSA based on the ABCG2/GAFF2 by reproducing the thermodynamic integration (TI) calculation results. The PB radius parameters of carbon, oxygen, sulfur, phosphorus, chloride, bromide and iodine, were adjusted. New atom types, on, oi, hn1, hn2, hn3, were introduced to further improve the fitting performance. Then, we tuned the parameters in the nonpolar SFE model using the experimental SFE data and the PB calculation results. By adopting the new radius parameters and new nonpolar SFE model, the root mean square error (RMSE) of the SFE calculation for the 544 molecules decreased from 2.38 to 1.05 kcal/mol. Finally, the new radius parameters were applied in the prediction of protein-ligand binding free energies using the MM-PBSA method. For the eight systems tested, we could observe higher correlation between the experiment data and calculation results and smaller prediction errors for the absolute binding free energies, demonstrating that our new radius parameters can improve the free energy calculation using the MM-PBSA method.

Role of Dynamics and Mutations in Interactions of a Zinc Finger Antiviral Protein with CG-rich Viral RNA

Zinc finger antiviral protein (ZAP) is a host antiviral factor that selectively inhibits the replication of a variety of viruses. ZAP recognizes the CG-enriched RNA sequences and activates the viral RNA degradation machinery. In this work, we investigated the dynamics of a ZAP/RNA complex and computed the energetics of mutations in ZAP that affect its binding to the viral RNA. The crystal structure of a mouse-ZAP/RNA complex showed that RNA interacts with the zinc finger 2 (ZF2) and ZF3 domains. However, we found that due to the dynamic behavior of the single-stranded RNA, the terminal nucleotides C1 and G2 of RNA change their positions from the ZF3 to the ZF1 domain. Moreover, the electrostatic interactions between the zinc ions and the viral RNA provide further stability to the ZAP/RNA complex. We also provide structural and thermodynamic evidence for seven residue pairs (C1-Arg74, C1-Arg179, G2-Arg74, U3-Lys76, C4-Lys76, G5-Arg95, and U6-Glu204) that show favorable ZAP/RNA interactions, although these interactions were not observed in the ZAP/RNA crystal structure. Consistent with the observations from the mouse-ZAP/RNA crystal structure, we found that four residue pairs (C4-Lys89, C4-Leu90, C4-Tyr108, and G5-Lys107) maintained stable interactions in MD simulations. Based on experimental mutagenesis studies and our residue-level interaction analysis, we chose seven residues (Arg74, Lys76, Lys89, Arg95, Lys107, Tyr108, and Arg179) for individual alanine mutations. In addition, we studied mutations in those residues that are only observed in the crystal structures as interacting with RNA (Tyr98, Glu148, and Arg170). Out of these 10 mutations, we found that the Ala mutation in each of the five residues Arg74, Lys76, Lys89, Lys107, and Glu148 significantly reduced the binding affinity of ZAP to RNA.

TIES 2.0: A Dual-Topology Open Source Relative Binding Free Energy Builder with Web Portal

Relative binding free energy (RBFE) calculations are widely used to aid the process of drug discovery. TIES, Thermodynamic Integration with Enhanced Sampling, is a dual-topology approach to RBFE calculations with support for NAMD and OpenMM molecular dynamics engines. The software has been thoroughly validated on publicly available datasets. Here we describe the open source software along with a web portal (https://ccs-ties.org) that enables users to perform such calculations correctly and rapidly.

Active Learning Guided Drug Design Lead Optimization Based on Relative Binding Free Energy Modeling

In silico identification of potent protein inhibitors commonly requires prediction of a ligand binding free energy (BFE). Thermodynamics integration (TI) based on molecular dynamics (MD) simulations is a BFE calculation method capable of acquiring accurate BFE, but it is computationally expensive and time-consuming. In this work, we have developed an efficient automated workflow for identifying compounds with the lowest BFE among thousands of congeneric ligands, which requires only hundreds of TI calculations. Automated machine learning (AutoML) orchestrated by active learning (AL) in an AL-AutoML workflow allows unbiased and efficient search for a small set of best-performing molecules. We have applied this workflow to select inhibitors of the SARS-CoV-2 papain-like protease and were able to find 133 compounds with improved binding affinity, including 16 compounds with better than 100-fold binding affinity improvement. We obtained a hit rate that outperforms that expected of traditional expert medicinal chemist-guided campaigns. Thus, we demonstrate that the combination of AL and AutoML with free energy simulations provides at least 20× speedup relative to the naïve brute force approaches.

Prioritizing Small Sets of Molecules for Synthesis through in-silico Tools: A Comparison of Common Ranking Methods

Prioritizing molecules for synthesis is a key role of computational methods within medicinal chemistry. Multiple tools exist for ranking molecules, from the cheap and popular molecular docking methods to more computationally expensive molecular-dynamics (MD)-based methods. It is often questioned whether the accuracy of the more rigorous methods justifies the higher computational cost and associated calculation time. Here, we compared the performance on ranking the binding of small molecules for seven scoring functions from five docking programs, one end-point method (MM/GBSA), and two MD-based free energy methods (PMX, FEP+). We investigated 16 pharmaceutically relevant targets with a total of 423 known binders. The performance of docking methods for ligand ranking was strongly system dependent. We observed that MD-based methods predominantly outperformed docking algorithms and MM/GBSA calculations. Based on our results, we recommend the application of MD-based free energy methods for prioritization of molecules for synthesis in lead optimization, whenever feasible.

Protein-Ligand Binding Free-Energy Calculations with ARROW─A Purely First-Principles Parameterized Polarizable Force Field

Protein-ligand binding free-energy calculations using molecular dynamics (MD) simulations have emerged as a powerful tool for in silico drug design. Here, we present results obtained with the ARROW force field (FF)─a multipolar polarizable and physics-based model with all parameters fitted entirely to high-level ab initio quantum mechanical (QM) calculations. ARROW has already proven its ability to determine solvation free energy of arbitrary neutral compounds with unprecedented accuracy. The ARROW FF parameterization is now extended to include coverage of all amino acids including charged groups, allowing molecular simulations of a series of protein-ligand systems and prediction of their relative binding free energies. We ensure adequate sampling by applying a novel technique that is based on coupling the Hamiltonian Replica exchange (HREX) with a conformation reservoir generated via potential softening and nonequilibrium MD. ARROW provides predictions with near chemical accuracy (mean absolute error of ∼0.5 kcal/mol) for two of the three protein systems studied here (MCL1 and Thrombin). The third protein system (CDK2) reveals the difficulty in accurately describing dimer interaction energies involving polar and charged species. Overall, for all of the three protein systems studied here, ARROW FF predicts relative binding free energies of ligands with a similar accuracy level as leading nonpolarizable force fields.

Practical Guidance for Consensus Scoring and Force Field Selection in Protein-Ligand Binding Free Energy Simulations

The advances in ligand binding affinity prediction have been fostered by system generation tools and improved force fields (FFs). CHARMM-GUI Free Energy Calculator provides input and postprocessing scripts for AMBER-TI free energy calculations with various FFs. In this study, we used 12 different FF combinations (ff14SB and ff19SB for protein, GAFF2.2 and OpenFF for ligand, and TIP3P, TIP4PEW, and OPC for water) to calculate relative binding free energies (ΔΔG_bind) for 80 alchemical transformations (among the JACS benchmark set) with different numbers of λ windows (12 or 21) and simulation times (1, 5, or 10 ns). Our results show that 12 λ windows and 5 ns simulation time for each window are sufficient to obtain reliable ΔΔG_bind with 4 independent runs for the current benchmark set. While there is no statistically noticeable performance difference among 12 different FF combinations compared to the experimental values, a combination of ff14SB + GAFF2.2 + TIP3P FFs appears to be best with a mean unsigned error of 0.87 [0.69, 1.07] kcal/mol, a root-mean-square error of 1.22 [0.94, 1.50] kcal/mol, a Pearson's correlation of 0.64 [0.52, 0.76], a Spearman's correlation of 0.73 [0.58, 0.83], and a Kendell's correlation of 0.54 [0.42, 0.64]. This large-scale ΔΔG_bind calculation study provides useful information about ΔΔG_bind prediction with different AMBER FF combinations and presents valuable suggestions for FF selection in AMBER-TI ΔΔG_bind calculations.

AMBER Drug Discovery Boost Tools: Automated Workflow for Production Free-Energy Simulation Setup and Analysis (ProFESSA)

We report an automated workflow for production free-energy simulation setup and analysis (ProFESSA) using the GPU-accelerated AMBER free-energy engine with enhanced sampling features and analysis tools, part of the AMBER Drug Discovery Boost package that has been integrated into the AMBER22 release. The workflow establishes a flexible, end-to-end pipeline for performing alchemical free-energy simulations that brings to bear technologies, including new enhanced sampling features and analysis tools, to practical drug discovery problems. ProFESSA provides the user with top-level control of large sets of free-energy calculations and offers access to the following key functionalities: (1) automated setup of file infrastructure; (2) enhanced conformational and alchemical sampling with the ACES method; and (3) network-wide free-energy analysis with the optional imposition of cycle closure and experimental constraints. The workflow is applied to perform absolute and relative solvation free-energy and relative ligand-protein binding free-energy calculations using different atom-mapping procedures. Results demonstrate that the workflow is internally consistent and highly robust. Further, the application of a new network-wide Lagrange multiplier constraint analysis that imposes key experimental constraints substantially improves binding free-energy predictions.

Molecular docking and molecular dynamics study Lianhua Qingwen granules (LHQW) treats COVID-19 by inhibiting inflammatory response and regulating cell survival

Purpose

2019 Coronavirus disease (COVID-19) is endangering health of populations worldwide. Latest research has proved that Lianhua Qingwen granules (LHQW) can reduce tissue damage caused by inflammatory reactions and relieve patients' clinical symptoms. However, the mechanism of LHQW treats COVID-19 is currently lacking. Therefore, we employed computer simulations to investigate the mechanism of LHQW treats COVID-19 by modulating inflammatory response.

Methods

We employed bioinformatics to screen active ingredients in LHQW and intersection gene targets. PPI, GO and KEGG was used to analyze relationship of intersection gene targets. Molecular dynamics simulations validated the binding stability of active ingredients and target proteins. Binding free energy, radius of gyration and the solvent accessible surface area were analyzed by supercomputer platform.

Results

COVID-19 had 4628 gene targets, LHQW had 1409 gene targets, intersection gene targets were 415. Bioinformatics analysis showed that intersection targets were closely related to inflammation and immunomodulatory. Molecular docking suggested that active ingredients (including: licopyranocoumarin, Glycyrol and 3-3-Oxopropanoic acid) in LHQW played a role in treating COVID-19 by acting on CSF2, CXCL8, CCR5, NLRP3, IFNG and TNF. Molecular dynamics was used to prove the binding stability of active ingredients and protein targets.

Conclusion

The mechanism of active ingredients in LHQW treats COVID-19 was investigated by computer simulations. We found that active ingredients in LHQW not only reduce cell damage and tissue destruction by inhibiting the inflammatory response through CSF2, CXCL8, CCR5 and IFNG, but also regulate cell survival and growth through NLRP3 and TNF thereby reducing apoptosis.

Identification of Selective BRD9 Inhibitor via Integrated Computational Approach

Bromodomain-containing protein 9 (BRD9), a member of the bromodomain and extra terminal domain (BET) protein family, works as an epigenetic reader. BRD9 has been considered an essential drug target for cancer, inflammatory diseases, and metabolic disorders. Due to its high similarity among other isoforms, no effective treatment of BRD9-associated disorders is available. For the first time, we performed a detailed comparative analysis among BRD9, BRD7, and BRD4. The results indicate that residues His42, Gly43, Ala46, Ala54, Val105, and Leu109 can confer the BRD9 isoform selectivity. The predicted crucial residues were further studied. The pharmacophore model's features were precisely mapped with some key residues including, Gly43, Phe44, Phe45, Asn100, and Tyr106, all of which play a crucial role in BRD9 inhibition. Docking-based virtual screening was utilized with the consideration of the conserved water network in the binding cavity to identify the potential inhibitors of BRD9. In this workflow, 714 compounds were shortlisted. To attain selectivity, 109 compounds were re-docked to BRD7 for negative selection. Finally, four compounds were selected for molecular dynamics studies. Our studies pave the way for the identification of new compounds and their role in causing noticeable, functional differences in isoforms and between orthologues.