Ligand Gaussian Accelerated Molecular Dynamics (LiGaMD): Characterization of Ligand Binding Thermodynamics and Kinetics

Computer-Aided Drug Discovery for Undruggable Targets

Undruggable targets are those of therapeutical significance but challenging for conventional drug design approaches. Such targets often exhibit unique features, including highly dynamic structures, a lack of well-defined ligand-binding pockets, the presence of highly conserved active sites, and functional modulation by protein-protein interactions. Recent advances in computational simulations and artificial intelligence have revolutionized the drug design landscape, giving rise to innovative strategies for overcoming these obstacles. In this review, we highlight the latest progress in computational approaches for drug design against undruggable targets, present several successful case studies, and discuss remaining challenges and future directions. Special emphasis is placed on four primary target categories: intrinsically disordered proteins, protein allosteric regulation, protein-protein interactions, and protein degradation, along with discussion of emerging target types. We also examine how AI-driven methodologies have transformed the field, from applications in protein-ligand complex structure prediction and virtual screening to de novo ligand generation for undruggable targets. Integration of computational methods with experimental techniques is expected to bring further breakthroughs to overcome the hurdles of undruggable targets. As the field continues to evolve, these advancements hold great promise to expand the druggable space, offering new therapeutic opportunities for previously untreatable diseases.

On the application of artificial intelligence in virtual screening

INTRODUCTION

Artificial intelligence (AI) has emerged as a transformative tool in drug discovery, particularly in virtual screening (VS), a crucial initial step in identifying potential drug candidates. This article highlights the significance of AI in revolutionizing both ligand-based virtual screening (LBVS) and structure-based virtual screening (SBVS) approaches, streamlining and enhancing the drug discovery process.

AREAS COVERED

The authors provide an overview of AI applications in drug discovery, with a focus on LBVS and SBVS approaches utilized in prospective cases where new bioactive molecules were identified and experimentally validated. Discussion includes the use of AI in quantitative structure-activity relationship (QSAR) modeling for LBVS, as well as its role in enhancing SBVS techniques such as molecular docking and molecular dynamics simulations. The article is based on literature searches on studies published up to March 2025.

EXPERT OPINION

AI is rapidly transforming VS in drug discovery, by leveraging increasing amounts of experimental data and expanding its scalability. These innovations promise to enhance efficiency and precision across both LBVS and SBVS approaches, yet challenges such as data curation, rigorous and prospective validation of new models, and efficient integration with experimental methods remain critical for realizing AI's full potential in drug discovery.

Design, optimization, and ADMET evaluation of S11a-0000168202: A promising LIMK1 inhibitor for gastric cancer treatment

This study focuses on the development and optimization of S11a-0000168202, a novel LIMK1 inhibitor with potential therapeutic applications in gastric cancer. Through scaffold hopping and structural modification of HIT100844099, S11a-0000168202 demonstrated enhanced binding stability and stronger interactions with key LIMK1 residues, including GLU-414, ILE-416, and HIS-464. Molecular dynamics simulations and MMGBSA analyses confirmed the compound's stability, while ADMET evaluation revealed favorable properties such as moderate lipophilicity, good human intestinal absorption, and low P-glycoprotein inhibition. Despite the promising computational results, the lack of experimental validation remains a limitation. Future studies should focus on in vitro and in vivo testing to confirm S11a-0000168202's efficacy, pharmacokinetics, and safety. This compound holds significant potential as a therapeutic agent for LIMK1-targeted gastric cancer treatment.

Identification and characterization of binding thermodynamics and kinetics of inhibitors targeting FGFR1 via molecular modelling and ligand Gaussian accelerated molecular dynamics simulations

Fibroblast growth factor receptor1 (FGFR1) kinase has a crucial role in cell proliferation, migration, and differentiation. Any imbalance in its level can cause cancer and many other illnesses. Despite the availability of numerous treatments, cytotoxicity, selectivity, and drug resistance issues demand the development of new FGFR1 inhibitors. Herein, we performed a high-throughput virtual screening of 54 624 compounds from NPASS and HMDB databases using the Schrodinger software suite. Compounds with a docking score cutoff of -11.0 kcal mol-1 were further screened for ADMET properties. Following the all-atom molecular dynamics simulation of selected molecules in replica, the binding free energy was calculated using the molecular mechanics Poisson Boltzmann surface area (MM-PBSA) scheme. We obtained two compounds, namely, bevantolol and 3-hydroxy glabrol, which exhibited higher binding affinities than the control drug ponatinib. Bevantolol was further optimized via virtual screening and simulation studies of its 100 structural analogues to obtain the best analogue. Subsequently, we investigated the binding thermodynamics and kinetics of the best analogue molecule via ligand Gaussian accelerated molecular dynamics (LiGaMD) simulation, performing independent replica runs of 4 μs each. The 1D and 2D potential of mean force and Kramer's rate theory determined the kinetic rate constants (kon/koff) associated with the FGFR1 complex. The binding constant was estimated to be 7.4 ± 0.27 nM, which was similar to the type II tyrosine kinase inhibitor ponatinib. Overall, this study highlights the dynamics of FGFR1-ligand interaction while proposing bevantolol and its analogue molecule ANLG-2 as promising drug candidates for FGFR1 therapeutic intervention.

Cryoelectron microscopy as a tool for illuminating activation mechanisms of human class A orphan G protein-coupled receptors

G protein-coupled receptors (GPCRs) are critically important medicinal targets, and the cryogenic electron microscopy (cryo-EM) revolution is providing novel high-resolution GPCR structures at a rapid pace. Orphan G protein-coupled receptors (oGPCRs) are a group of approximately 100 nonolfactory GPCRs for which endogenous ligands are unknown or not validated. The absence of modulating ligands adds difficulties to understanding the physiologic significance of oGPCRs and in the determination of high-resolution structures of isolated receptors that could facilitate drug discovery. Despite the challenges, cryo-EM structures of oGPCR-G protein complexes are emerging. This is being facilitated by numerous developments to stabilize GPCR-G protein complexes such as the use of dominant-negative G proteins, mini-G proteins, complex-stabilizing nanobodies or antibody fragments, and protein tethering methods. Moreover, many oGPCRs are constitutively active, which can facilitate complex formation in the absence of a known activating ligand. Consequently, in addition to providing templates for drug discovery, active oGPCR structures shed light on constitutive GPCR activation mechanisms. These comprise self-activation, whereby mobile extracellular portions of the receptor act as tethered agonists by occupying a canonical orthosteric-binding site in the transmembrane core, constitutive activity due to alterations to conserved molecular switches that stabilize inactive states of GPCRs, as well as receptors activated by cryptic ligands that are copurified with the receptor. Cryo-EM structures of oGPCRs are now being determined at a rapid pace and are expected to be invaluable tools for oGPCR drug discovery. SIGNIFICANCE STATEMENT: Orphan G protein-coupled receptors (GPCRs) provide large untapped potential for development of new medicines. Many of these receptors display constitutive activity, enabling structure determination and insights into observed GPCR constitutive activity including (1) self-activation by mobile receptor extracellular portions that function as tethered agonists, (2) modification of conserved motifs canonically involved in receptor quiescence and/or activation, and (3) activation by cryptic lipid ligands. Collectively, these studies advance fundamental understanding of GPCR function and provide opportunities for novel drug discovery.

Dissecting Large-Scale Structural Transitions in Membrane Transporters Using Advanced Simulation Technologies

Membrane transporters are integral membrane proteins that act as gatekeepers of the cell, controlling fundamental processes such as recruitment of nutrients and expulsion of waste material. At a basic level, transporters operate using the "alternating access model," in which transported substances are accessible from only one side of the membrane at a time. This model usually involves large-scale structural changes in the transporter, which often cannot be captured using unbiased, conventional molecular simulation techniques. In this article, we provide an overview of some of the major simulation techniques that have been applied to characterize the structural dynamics and energetics involved in the transition of membrane transporters between their functional states. After briefly introducing each technique, we discuss some of their advantages and limitations and provide some recent examples of their application to membrane transporters.

Diverse toxins exhibit a common binding mode to the nicotinic acetylcholine receptors

Nicotinic acetylcholine receptors (nAChRs) are critical ligand-gated ion channels in the human nervous system. They are targets for various neurotoxins produced by algae, plants, and animals. While many structures of nAChRs bound by neurotoxins have been published, the binding mechanism of toxins to the nAChRs remains unclear. In this work, we have performed extensive Gaussian accelerated molecular dynamics simulations on several Aplysia californica nAChRs in complex with α-conotoxins, strychnine, and pinnatoxins, as well as human nAChRs in complex with α-bungarotoxin and α-conotoxin, to determine the binding and dissociation pathways of the toxins to the nAChRs and the associated effects. We uncovered two common binding and dissociation pathways shared by toxins and nAChRs. In the first binding pathway, the toxins diffused from the bulk solvent to bind a region near the extracellular pore before moving downwards along the nAChRs to the nAChR orthosteric pocket. The second binding pathway involved a direct diffusion of the toxins from the bulk solvent into the nAChR orthosteric pocket. The dissociation pathways were the reverse of the observed binding pathways. Notably, we determined that the electrostatically bipolar interactions between the nAChR orthosteric pocket and toxins provided an explanation for the common binding mode shared by diverse toxins.

Effect of hexyloxy position on antagonistic properties of KH-5 (1-{2-[4-(hexyloxy)benzoyloxy]ethyl}-1-methyl-1,2,3,6-tetrahydropyridin-1-ium iodide) at muscarinic acetylcholine receptors

Antagonists with a long residence time at the receptors are desired for the possibility of reducing daily doses and side effects. KH-5 (1-{2-[4-(hexyloxy)benzoyloxy]ethyl}-1-methyl-1,2,3,6-tetrahydropyridin-1-ium iodide) is the long-acting M1-preferring bitopic muscarinic antagonist with a half-life at muscarinic receptors of up to five hours. The binding of 2-hexyloxy and 3-hexyloxy analogues of KH-5 was simulated in silico, compounds were synthesized and their binding and antagonistic properties were measured experimentally in CHO cells expressing individual subtypes of muscarinic acetylcholine receptors. The overall binding affinities of the new compounds were similar to their respective parent compounds. Shifting the hexyloxy chain to ortho and meta positions led to a decrease in potency at the M1 receptor but an increase in potency at the M2 receptor and abolition of long-term antagonism. Preservation of the para position of the hexyloxy chain is essential for the further development of M1-preferring antagonists. Modifications of the basic centre may be the way to improve the geometry of antagonists towards long residence times to obtain the desired long-acting muscarinic antagonists in the future. The additional challenge for further development is the low metabolic stability of compounds.

Mechanistic Perspective on C-N and C-S Bond Construction Catalyzed by Cytochrome P450 Enzymes

Cytochrome P450 enzymes catalyze a large number of oxidative transformations that are responsible for natural product synthesis. Recent studies have revealed their unique ability to catalyze the formation of C-N and C-S bonds, broadening their biosynthetic applications. However, the enzymatic mechanisms behind these reactions are still unclear. This review focuses on theoretical insights into the mechanisms of P450-catalyzed C-N and C-S bond formation. The key roles of the conformational dynamics of substrate radicals, influenced by the enzyme environment, in modulating selectivity and reactivity are highlighted. Understanding these reaction mechanisms offers valuable guidance for P450 enzyme engineering and the design of biosynthetic applications.

ModBind, a Rapid Simulation-Based Predictor of Ligand Binding and Off-Rates

In rational drug discovery, both free energy of binding and the binding half-life (koff) are important factors in determining the efficacy of drugs. Numerous computational methods have been developed to predict these important properties, many of which rely on molecular dynamics (MD) simulations. While binding free-energy methods (thermodynamic equilibrium predictions) have been well validated and have demonstrated the ability to drive daily synthesis decisions in a commercial drug discovery setting, the prediction of koff (kinetics predictions) has had limited validation, and predictive methods have largely not been deployed in drug discovery settings. We developed ModBind, a novel method for MD simulation-based koff predictions. ModBind demonstrated similar accuracy to current state-of-the-art free-energy prediction methods. Additionally, ModBind performs ∼100 times faster than most available MD simulation-based free-energy or koff methods, allowing for widespread use by the molecular modeling community. While most free-energy methods rely on relative free-energy changes and are primarily useful for optimization of a congeneric series, our method requires no structural similarity between ligands, making ModBind an absolute predictor of koff. ModBind is thus a tool that can be used in virtual screening of diverse ligands, making it distinct from relative free-energy methods. We also discuss conditions that enable approximate prediction of ligand efficacy using ModBind and the limitations of this approach.

Potency Prediction of Covalent Inhibitors against SARS-CoV-2 3CL-like Protease and Multiple Mutants by Multiscale Simulations

3-Chymotrypsin-like protease (3CLpro) is a prominent target against pathogenic coronaviruses. Expert knowledge of the cysteine-targeted covalent reaction mechanism is crucial to predict the inhibitory potency of approved inhibitors against 3CLpros of SARS-CoV-2 variants and perform structure-based drug design against newly emerging coronaviruses. We carried out an extensive array of classical and hybrid QM/MM molecular dynamics simulations to explore covalent inhibition mechanisms of five well-characterized inhibitors toward SARS-CoV-2 3CLpro and its mutants. The calculated binding affinity and reactivity of the inhibitors are highly consistent with experimental data, and the predicted inhibitory potency of the inhibitors against 3CLpro with L167F, E166V, or T21I/E166V mutant is in full agreement with IC50s determined by the accompanying enzymatic assays. The explored mechanisms unveil the impact of residue mutagenesis on structural dynamics that communicates to change not only noncovalent binding strength but also covalent reaction free energy. Such a change is inhibitor dependent, corresponding to varied levels of drug resistance of these 3CLpro mutants against nirmatrelvir and simnotrelvir and no resistance to the 11a compound. These results together suggest that the present simulations with a suitable protocol can efficiently evaluate the reactivity and potency of covalent inhibitors along with the elucidated molecular mechanisms of covalent inhibition.

Experimental and computational biophysics to identify vasodilator drugs targeted at TRPV2 using agonists based on the probenecid scaffold

TRP channels are important pharmacological targets in physiopathology. TRPV2 plays distinct roles in cardiac and neuromuscular function, immunity, and metabolism, and is associated with pathologies like muscular dystrophy and cancer. However, TRPV2 pharmacology is unspecific and scarce at best. Using in silico similarity-based chemoinformatics we obtained a set of 270 potential hits for TRPV2 categorized into families based on chemical nature and similarity. Docking the compounds on available rat TRPV2 structures allowed the clustering of drug families in specific ligand binding sites. Starting from a probenecid docking pose in the piperlongumine binding site and using a Gaussian accelerated molecular dynamics approach we have assigned a putative probenecid binding site. In parallel, we measured the EC50 of 7 probenecid derivatives on TRPV2 expressed in Pichia pastoris using a novel medium-throughput Ca2+ influx assay in yeast membranes together with an unbiased and unsupervised data analysis method. We found that 4-(piperidine-1-sulfonyl)-benzoic acid had a better EC50 than probenecid, which is one of the most specific TRPV2 agonists to date. Exploring the TRPV2-dependent anti-hypertensive potential in vivo, we found that 4-(piperidine-1-sulfonyl)-benzoic acid shows a sex-biased vasodilator effect producing larger vascular relaxations in female mice. Overall, this study expands the pharmacological toolbox for TRPV2, a widely expressed membrane protein and orphan drug target.

Advances and Challenges in Milestoning Simulations for Drug-Target Kinetics

Molecular dynamics simulations have become indispensable for exploring complex biological processes, yet their limitations in capturing rare events hinder our understanding of drug-target kinetics. In this Perspective, we investigate the domain of milestoning simulations to understand this challenge. The milestoning approach divides the phase space of the drug-target complex into discrete cells, offering extended time scale insights. This Perspective traces the history, applications, and future potential of milestoning simulations in the context of drug-target kinetics. It explores the fundamental principles of milestoning, highlighting the importance of probabilistic transitions and transition time independence. Markovian milestoning with Voronoi tessellations is revisited to address the traditional milestoning challenges. While observing the advancements in this field, this Perspective also addresses impending challenges in estimating drug-target unbinding rate constants through milestoning simulations, paving the way for more effective drug design strategies.

Mixture-of-Experts Based Dissociation Kinetic Model for De Novo Design of HSP90 Inhibitors with Prolonged Residence Time

The dissociation rate constant (koff) significantly impacts the drug potency and dosing frequency. This work proposes a powerful optimization-based framework for de novo drug design guided by koff. First, a comprehensive database containing 2,773 unique koff values is created. Based on the database, a novel generic dissociation kinetic model is developed with a mixture-of-experts architecture, enabling high-throughput predictions of koff with high accuracy. The developed model is then integrated with an optimization-based mathematical programming approach to design drug candidates with low koff. Finally, the τ-RAMD method is utilized to rigorously verify the designed potential drug candidates. In a case study, the framework successfully identified numerous new potential HSP90 inhibitor candidates, achieving a maximum 45.7% improvement in residence time (τ = 1/koff) compared to that of a known exceptional HSP90 inhibitor. These findings demonstrate the feasibility and effectiveness of the kinetics-guided optimization-based de novo drug design framework in designing drug candidates with prolonged τ.

Harnessing AlphaFold to reveal hERG channel conformational state secrets

To design safe, selective, and effective new therapies, there must be a deep understanding of the structure and function of the drug target. One of the most difficult problems to solve has been resolution of discrete conformational states of transmembrane ion channel proteins. An example is KV11.1 (hERG), comprising the primary cardiac repolarizing current, I Kr. hERG is a notorious drug anti-target against which all promising drugs are screened to determine potential for arrhythmia. Drug interactions with the hERG inactivated state are linked to elevated arrhythmia risk, and drugs may become trapped during channel closure. However, the structural details of multiple conformational states have remained elusive. Here, we guided AlphaFold2 to predict plausible hERG inactivated and closed conformations, obtaining results consistent with multiple available experimental data. Drug docking simulations demonstrated hERG state-specific drug interactions in good agreement with experimental results, revealing that most drugs bind more effectively in the inactivated state and are trapped in the closed state. Molecular dynamics simulations demonstrated ion conduction for an open but not AlphaFold2 predicted inactivated state that aligned with earlier studies. Finally, we identified key molecular determinants of state transitions by analyzing interaction networks across closed, open, and inactivated states in agreement with earlier mutagenesis studies. Here, we demonstrate a readily generalizable application of AlphaFold2 as an effective and robust method to predict discrete protein conformations, reconcile seemingly disparate data and identify novel linkages from structure to function.

Molecular Dynamics Unveils Multiple-Site Binding of Inhibitors with Reduced Activity on the Surface of Dihydrofolate Reductase

The sensitivity to protein inhibitors is altered by modifications or protein mutations, as represented by drug resistance. The mode of stable drug binding to the protein pocket has been experimentally clarified. However, the nature of the binding of inhibitors with reduced sensitivity remains unclear at the atomic level. In this study, we analyzed the thermodynamics and kinetics of inhibitor binding to the surface of wild-type and mutant dihydrofolate reductase (DHFR) using molecular dynamics simulations combined with Markov state modeling. A strong inhibitor of methotrexate (MTX) showed a preference for the active site of wild-type DHFR with minimal binding to unrelated (secondary) sites. Deletion of a side-chain fragment in MTX largely destabilized the active site-bound state, with clear evidence of binding to secondary sites. Similarly, the F31V mutation in DHFR diminished the specificity of MTX binding to the active site. These results reveal the presence of multiple-bound states whose stabilities are comparable to or higher than those of the unbound state, suggesting that a reduction in the binding affinity for the active site significantly elevates the fractions of these states. This study presents a theoretical model that more accurately interprets the altered drug sensitivity than the traditional two-state model.

Large-scale annotation of biochemically relevant pockets and tunnels in cognate enzyme-ligand complexes

Tunnels in enzymes with buried active sites are key structural features allowing the entry of substrates and the release of products, thus contributing to the catalytic efficiency. Targeting the bottlenecks of protein tunnels is also a powerful protein engineering strategy. However, the identification of functional tunnels in multiple protein structures is a non-trivial task that can only be addressed computationally. We present a pipeline integrating automated structural analysis with an in-house machine-learning predictor for the annotation of protein pockets, followed by the calculation of the energetics of ligand transport via biochemically relevant tunnels. A thorough validation using eight distinct molecular systems revealed that CaverDock analysis of ligand un/binding is on par with time-consuming molecular dynamics simulations, but much faster. The optimized and validated pipeline was applied to annotate more than 17,000 cognate enzyme-ligand complexes. Analysis of ligand un/binding energetics indicates that the top priority tunnel has the most favourable energies in 75% of cases. Moreover, energy profiles of cognate ligands revealed that a simple geometry analysis can correctly identify tunnel bottlenecks only in 50% of cases. Our study provides essential information for the interpretation of results from tunnel calculation and energy profiling in mechanistic enzymology and protein engineering. We formulated several simple rules allowing identification of biochemically relevant tunnels based on the binding pockets, tunnel geometry, and ligand transport energy profiles.Scientific contributionsThe pipeline introduced in this work allows for the detailed analysis of a large set of protein-ligand complexes, focusing on transport pathways. We are introducing a novel predictor for determining the relevance of binding pockets for tunnel calculation. For the first time in the field, we present a high-throughput energetic analysis of ligand binding and unbinding, showing that approximate methods for these simulations can identify additional mutagenesis hotspots in enzymes compared to purely geometrical methods. The predictor is included in the supplementary material and can also be accessed at https://github.com/Faranehhad/Large-Scale-Pocket-Tunnel-Annotation.git . The tunnel data calculated in this study has been made publicly available as part of the ChannelsDB 2.0 database, accessible at https://channelsdb2.biodata.ceitec.cz/ .

Self-organizing maps of unbiased ligand-target binding pathways and kinetics

The interpretation of ligand-target interactions at atomistic resolution is central to most efforts in computational drug discovery and optimization. However, the highly dynamic nature of protein targets, as well as possible induced fit effects, makes difficult to sample many interactions effectively with docking studies or even with large-scale molecular dynamics (MD) simulations. We propose a novel application of Self-Organizing Maps (SOMs) to address the sampling and dynamic mapping tasks, particularly in cases involving ligand flexibility and induced fit. The SOM approach offers a data-driven strategy to create a map of the interaction process and pathways based on unbiased MD. Furthermore, we show how the preliminary SOM mapping is complementary to kinetic analysis, with the employment of both network-based approaches and Markov state models. We demonstrate the method by comprehensively mapping a large dataset of 640 μs of unbiased trajectories sampling the recognition process between the phosphorylated YEEI peptide and its high-specificity target lck-SH2. The integration of SOM into unbiased simulation protocols significantly advances our understanding of the ligand binding mechanism. This approach serves as a potent tool for mapping intricate ligand-target interactions with unprecedented detail, thereby enhancing the characterization of kinetic properties crucial to drug design.

Structural Models of Human Norepinephrine Transporter Ensemble Reveal the Allosteric Sites and Ligand-Binding Mechanism

The norepinephrine transporter (NET) plays a pivotal role in recycling norepinephrine (NE) from the synaptic cleft. However, the structures referring to the conformational heterogeneity of NET during the transport cycle remain poorly understood. Here, three structural models of NE bound to the orthosteric site of NET in outward-open (OOholo), outward-occluded (OCholo), and inward-open (IOholo) conformations were first obtained using the multistate structures of serotonin transporter as templates and further characterized through Gaussian-accelerated molecular dynamics and free energy reweighting. Analysis of the structures revealed eight potential allosteric sites on the functional-specific states of NET. One of the pharmacologically relevant pockets located at the extracellular vestibule was further verified by simulating the binding behaviors of a clinical trial drug χ-MrIA that is allosterically regulating NET. These structural and energetic insights into NET advanced our understanding of NE reuptake and paved the way for discovering novel molecules targeting the allosteric sites.

Reinforcing Tunnel Network Exploration in Proteins Using Gaussian Accelerated Molecular Dynamics

Tunnels are structural conduits in biomolecules responsible for transporting chemical compounds and solvent molecules from the active site. They have been shown to be present in a wide variety of enzymes across all functional and structural classes. However, the study of such pathways is experimentally challenging, because they are typically transient. Computational methods, such as molecular dynamics (MD) simulations, have been successfully proposed to explore tunnels. Conventional MD (cMD) provides structural details to characterize tunnels but suffers from sampling limitations to capture rare tunnel openings on longer time scales. Therefore, in this study, we explored the potential of Gaussian accelerated MD (GaMD) simulations to improve the exploration of complex tunnel networks in enzymes. We used the haloalkane dehalogenase LinB and its two variants with engineered transport pathways, which are not only well-known for their application potential but have also been extensively studied experimentally and computationally regarding their tunnel networks and their importance in multistep catalytic reactions. Our study demonstrates that GaMD efficiently improves tunnel sampling and allows the identification of all known tunnels for LinB and its two mutants. Furthermore, the improved sampling provided insight into a previously unknown transient side tunnel (ST). The extensive conformational landscape explored by GaMD simulations allowed us to investigate in detail the mechanism of ST opening. We determined variant-specific dynamic properties of ST opening, which were previously inaccessible due to limited sampling of cMD. Our comprehensive analysis supports multiple indicators of the functional relevance of the ST, emphasizing its potential significance beyond structural considerations. In conclusion, our research proves that the GaMD method can overcome the sampling limitations of cMD for the effective study of tunnels in enzymes, providing further means for identifying rare tunnels in enzymes with the potential for drug development, precision medicine, and rational protein engineering.

Activating Point Mutations in the MET Kinase Domain Represent a Unique Molecular Subset of Lung Cancer and Other Malignancies Targetable with MET Inhibitors

Activating point mutations in the MET tyrosine kinase domain (TKD) are oncogenic in a subset of papillary renal cell carcinomas. Here, using comprehensive genomic profiling among >600,000 patients, we identify activating MET TKD point mutations as putative oncogenic driver across diverse cancers, with a frequency of ∼0.5%. The most common mutations in the MET TKD defined as oncogenic or likely oncogenic according to OncoKB resulted in amino acid substitutions at positions H1094, L1195, F1200, D1228, Y1230, M1250, and others. Preclinical modeling of these alterations confirmed their oncogenic potential and also demonstrated differential patterns of sensitivity to type I and type II MET inhibitors. Two patients with metastatic lung adenocarcinoma harboring MET TKD mutations (H1094Y, F1200I) and no other known oncogenic drivers achieved confirmed partial responses to a type I MET inhibitor. Activating MET TKD mutations occur in multiple malignancies and may confer clinical sensitivity to currently available MET inhibitors. Significance: The identification of targetable genomic subsets of cancer has revolutionized precision oncology and offers patients treatments with more selective and effective agents. Here, we demonstrate that activating, oncogenic MET tyrosine kinase domain mutations are found across a diversity of cancer types and are responsive to MET tyrosine kinase inhibitors.

Pathways and mechanism of MRTX1133 binding to KRAS G12D elucidated by molecular dynamics simulations and Markov state models

KRAS G12D is the most common oncogenic mutation identified in several types of cancer. Therefore, design of inhibitors targeting KRAS G12D represents a promising strategy for anticancer therapy. MRTX1133 is a highly potent inhibitor (approximate experiment Kd ≈ 0.0002 nM) of KRAS G12D and is currently in Phase 1/2 study, however, pathways of the compound binding to KRAS G12D has remained unknown, and the mechanism underlying the complicated dynamic process are challenging to capture experimentally, which hinder the structure-based anti-cancer drug design. Here, using MRTX1133 as a probe, unbiased molecular dynamics (MD) was used to simulate the process of MRTX1133 spontaneously binding to KRAS G12D. In six of 42 independent MD simulation (a total of 99 μs), MRTX1133 was observed to successfully associate with KRAS G12D. The kinetically metastable states refer to the potential pathways of MRTX1133 binding to KRAS G12D were revealed by Markov state models (MSM) analysis. Additionally, 8 key residues that are essential for MRTX1133 recognition and tight binding at the preferred low energy states were identified by MM/GBSA analysis. In sum, this study provides a new perspective on understanding the pathways and mechanism of MRTX1133 binding to KRAS G12D.

Increased throughput in methods for simulating protein ligand binding and unbinding

By incorporating full flexibility and enabling the quantification of crucial parameters such as binding free energies and residence times, methods for investigating protein-ligand binding and unbinding via molecular dynamics provide details on the involved mechanisms at the molecular level. While these advancements hold promise for impacting drug discovery, a notable drawback persists: their relatively time-consuming nature limits throughput. Herein, we survey recent implementations which, employing a blend of enhanced sampling techniques, a clever choice of collective variables, and often machine learning, strive to enhance the efficiency of new and previously reported methods without compromising accuracy. Particularly noteworthy is the validation of these methods that was often performed on systems mirroring real-world drug discovery scenarios.

Ligand Gaussian Accelerated Molecular Dynamics 3 (LiGaMD3): Improved Calculations of Binding Thermodynamics and Kinetics of Both Small Molecules and Flexible Peptides

Binding thermodynamics and kinetics play critical roles in drug design. However, it has proven challenging to efficiently predict ligand binding thermodynamics and kinetics of small molecules and flexible peptides using conventional molecular dynamics (cMD), due to limited simulation time scales. Based on our previously developed ligand Gaussian accelerated molecular dynamics (LiGaMD) method, we present a new approach, termed "LiGaMD3″, in which we introduce triple boosts into three individual energy terms that play important roles in small-molecule/peptide dissociation, rebinding, and system conformational changes to improve the sampling efficiency of small-molecule/peptide interactions with target proteins. To validate the performance of LiGaMD3, MDM2 bound by a small molecule (Nutlin 3) and two highly flexible peptides (PMI and P53) were chosen as the model systems. LiGaMD3 could efficiently capture repetitive small-molecule/peptide dissociation and binding events within 2 μs simulations. The predicted binding kinetic constant rates and free energies from LiGaMD3 were in agreement with the available experimental values and previous simulation results. Therefore, LiGaMD3 provides a more general and efficient approach to capture dissociation and binding of both small-molecule ligands and flexible peptides, allowing for accurate prediction of their binding thermodynamics and kinetics.

Molecular basis for the activation of PAF receptor by PAF

Platelet-activating factor (PAF) is a potent phospholipid mediator crucial in multiple inflammatory and immune responses through binding and activating the PAF receptor (PAFR). However, drug development targeting the PAFR has been limited, partly due to an incomplete understanding of its activation mechanism. Here, we present a 2.9-Å structure of the PAF-bound PAFR-Gi complex. Structural and mutagenesis analyses unveil a specific binding mode of PAF, with the choline head forming cation-π interactions within PAFR hydrophobic pocket, while the alkyl tail penetrates deeply into an aromatic cleft between TM4 and TM5. Binding of PAF modulates conformational changes in key motifs of PAFR, triggering the outward movement of TM6, TM7, and helix 8 for G protein coupling. Molecular dynamics simulation suggests a membrane-side pathway for PAF entry into PAFR via the TM4-TM5 cavity. By providing molecular insights into PAFR signaling, this work contributes a foundation for developing therapeutic interventions targeting PAF signal axis.

Insights into the Interaction Mechanisms of Peptide and Non-Peptide Inhibitors with MDM2 Using Gaussian-Accelerated Molecular Dynamics Simulations and Deep Learning

Inhibiting MDM2-p53 interaction is considered an efficient mode of cancer treatment. In our current study, Gaussian-accelerated molecular dynamics (GaMD), deep learning (DL), and binding free energy calculations were combined together to probe the binding mechanism of non-peptide inhibitors K23 and 0Y7 and peptide ones PDI6W and PDI to MDM2. The GaMD trajectory-based DL approach successfully identified significant functional domains, predominantly located at the helixes α2 and α2', as well as the β-strands and loops between α2 and α2'. The post-processing analysis of the GaMD simulations indicated that inhibitor binding highly influences the structural flexibility and collective motions of MDM2. Calculations of molecular mechanics-generalized Born surface area (MM-GBSA) and solvated interaction energy (SIE) not only suggest that the ranking of the calculated binding free energies is in agreement with that of the experimental results, but also verify that van der Walls interactions are the primary forces responsible for inhibitor-MDM2 binding. Our findings also indicate that peptide inhibitors yield more interaction contacts with MDM2 compared to non-peptide inhibitors. Principal component analysis (PCA) and free energy landscape (FEL) analysis indicated that the piperidinone inhibitor 0Y7 shows the most pronounced impact on the free energy profiles of MDM2, with the piperidinone inhibitor demonstrating higher fluctuation amplitudes along primary eigenvectors. The hot spots of MDM2 revealed by residue-based free energy estimation provide target sites for drug design toward MDM2. This study is expected to provide useful theoretical aid for the development of selective inhibitors of MDM2 family members.

Computational Insights into the Intramolecular Aromatic C-C Coupling Catalyzed by the Cytochrome P450 Enzyme CYP121 from Mycobacterium tuberculosis

CYP121 is a P450 enzyme that catalyzes the intramolecular C-C coupling of its native substrate, dicyclotyrosine (cYY). According to previous suggestions, when the cosubstrate peracetic acid was used to generate Cpd I, the substrate cYY was suggested to participate in the cleavage of the O-O bond; however, whether cYY is involved in the formation of Cpd I and how two distant aromatic carbon atoms are activated are still unclear. Here, we constructed computational models and performed QM/MM calculations to clarify the reaction mechanism. On the basis of our calculation results, cYY is not involved in the formation of Cpd I, and the C-C coupling reaction starts from hydrogen abstraction. In the second stage, the substrate should first undergo a complex conformational change, leading to two phenolic hydroxyls of cYY close to each other. In the subsequent reaction, the resultant Cpd II again abstracts a hydrogen atom from the proximal tyrosine to generate the diradical intermediate. In addition, the C-C coupling occurs in the active site, but the final aromatization may be a nonenzymatic reaction. In general, the intramolecular C-C coupling requires two basic conditions, including the active site having good flexibility and the substrate itself having a suitable and rotatable skeleton.

Learning Protein-Ligand Unbinding Pathways via Single-Parameter Community Detection

Understanding the dynamics of biomolecular complexes, e.g., of protein-ligand (un)binding, requires the comprehension of paths such systems take between metastable states. In MD simulations, paths are usually not observable per se, but they need to be inferred from simulation trajectories. Here, we present a novel approach to cluster trajectories based on a community detection algorithm that necessitates only the definition of a single parameter. The unbinding of the streptavidin-biotin complex is used as a benchmark system and the A2a adenosine receptor in complex with the inhibitor ZM241385 as an elaborate application. We demonstrate how such clusters of trajectories correspond to pathways and how the approach helps in the identification of reaction coordinates for a considered (un)binding process.

Machine Learning Deciphered Molecular Mechanistics with Accurate Kinetic and Thermodynamic Prediction

Time-lagged independent component analysis (tICA) and the Markov state model (MSM) have been extensively employed for extracting conformational dynamics and kinetic community networks from unbiased trajectory ensembles. However, these techniques may not be the optimal choice for elucidating transition mechanisms within low-dimensional representations, especially for intricate biosystems. Unraveling the association mechanism in such complex systems always necessitates permutations of several essential independent components or collective variables, a process that is inherently obscure and may require empirical knowledge for selection. To address these challenges, we have implemented an integrated unsupervised dimension reduction model: uniform manifold approximation and projection (UMAP) with hierarchy density-based spatial clustering of applications with noise (HDBSCAN). This approach effectively generates low-dimensional configurational embeddings. The hierarchical application of this architecture, in conjunction with MSM, reveals global kinetic connectivity while identifying local conformational states. Consequently, our methodology establishes a multiscale mechanistic elucidation framework. Leveraging the benefits of the uniform sample distribution and a denoising approach, our model demonstrates robustness in preserving global and local data structures compared to traditional dimension reduction methods in the field of MD analysis area. The interpretability of hyperparameter selection and compatibility with downstream tasks are cross-validated across various simulation data sets, utilizing both computational evaluation metrics and experimental kinetic observables. Furthermore, the predicted Mcl1-BH3 association kinetics (0.76 s-1) is in close agreement with surface plasmon resonance experiments (0.12 s-1), affirming the plausibility of the identified pathway composed of representative conformations. We anticipate that the devised workflow will serve as a foundational framework for studying recognition patterns in complex biological systems. Its contributions extend to the exploration of protein functional dynamics and rational drug design, offering a potent avenue for advancing research in these domains.

Molecular Mechanisms of the Impaired Heparin Pentasaccharide Interactions in 10 Antithrombin Heparin Binding Site Mutants Revealed by Enhanced Sampling Molecular Dynamics

Antithrombin (AT) is a critical regulator of the coagulation cascade by inhibiting multiple coagulation factors including thrombin and FXa. Binding of heparinoids to this serpin enhances the inhibition considerably. Mutations located in the heparin binding site of AT result in thrombophilia in affected individuals. Our aim was to study 10 antithrombin mutations known to affect their heparin binding in a heparin pentasaccharide bound state using two molecular dynamics (MD) based methods providing enhanced sampling, GaMD and LiGaMD2. The latter provides an additional boost to the ligand and the most important binding site residues. From our GaMD simulations we were able to identify four variants (three affecting amino acid Arg47 and one affecting Lys114) that have a particularly large effect on binding. The additional acceleration provided by LiGaMD2 allowed us to study the consequences of several other mutants including those affecting Arg13 and Arg129. We were able to identify several conformational types by cluster analysis. Analysis of the simulation trajectories revealed the causes of the impaired pentasaccharide binding including pentasaccharide subunit conformational changes and altered allosteric pathways in the AT protein. Our results provide insights into the effects of AT mutations interfering with heparin binding at an atomic level and can facilitate the design or interpretation of in vitro experiments.

Recent advances in computational and experimental protein-ligand affinity determination techniques

INTRODUCTION

Modern drug discovery revolves around designing ligands that target the chosen biomolecule, typically proteins. For this, the evaluation of affinities of putative ligands is crucial. This has given rise to a multitude of dedicated computational and experimental methods that are constantly being developed and improved.

AREAS COVERED

In this review, the authors reassess both the industry mainstays and the newest trends among the methods for protein - small-molecule affinity determination. They discuss both computational affinity predictions and experimental techniques, describing their basic principles, main limitations, and advantages. Together, this serves as initial guide to the currently most popular and cutting-edge ligand-binding assays employed in rational drug design.

EXPERT OPINION

The affinity determination methods continue to develop toward miniaturization, high-throughput, and in-cell application. Moreover, the availability of data analysis tools has been constantly increasing. Nevertheless, cross-verification of data using at least two different techniques and careful result interpretation remain of utmost importance.

Understanding the impact of binding free energy and kinetics calculations in modern drug discovery

INTRODUCTION

For rational drug design, it is crucial to understand the receptor-drug binding processes and mechanisms. A new era for the use of computer simulations in predicting drug-receptor interactions at an atomic level has begun with remarkable advances in supercomputing and methodological breakthroughs.

AREAS COVERED

End-point free energy calculation methods such as Molecular Mechanics/Poisson Boltzmann Surface Area (MM/PBSA) or Molecular-Mechanics/Generalized Born Surface Area (MM/GBSA), free energy perturbation (FEP), and thermodynamic integration (TI) are commonly used for binding free energy calculations in drug discovery. In addition, kinetic dissociation and association rate constants (k off and k on ) play critical roles in the function of drugs. Nowadays, Molecular Dynamics (MD) and enhanced sampling simulations are increasingly being used in drug discovery. Here, the authors provide a review of the computational techniques used in drug binding free energy and kinetics calculations.

EXPERT OPINION

The applications of computational methods in drug discovery and design are expanding, thanks to improved predictions of the binding free energy and kinetic rates of drug molecules. Recent microsecond-timescale enhanced sampling simulations have made it possible to accurately capture repetitive ligand binding and dissociation, facilitating more efficient and accurate calculations of ligand binding free energy and kinetics.

Dynamical Reweighting for Biased Rare Event Simulations

Dynamical reweighting techniques aim to recover the correct molecular dynamics from a simulation at a modified potential energy surface. They are important for unbiasing enhanced sampling simulations of molecular rare events. Here, we review the theoretical frameworks of dynamical reweighting for modified potentials. Based on an overview of kinetic models with increasing level of detail, we discuss techniques to reweight two-state dynamics, multistate dynamics, and path integrals. We explore the natural link to transition path sampling and how the effect of nonequilibrium forces can be reweighted. We end by providing an outlook on how dynamical reweighting integrates with techniques for optimizing collective variables and with modern potential energy surfaces.

Pharmacological modulation of septins restores calcium homeostasis and is neuroprotective in models of Alzheimer's disease

Abnormal calcium signaling is a central pathological component of Alzheimer's disease (AD). Here, we describe the identification of a class of compounds called ReS19-T, which are able to restore calcium homeostasis in cell-based models of tau pathology. Aberrant tau accumulation leads to uncontrolled activation of store-operated calcium channels (SOCCs) by remodeling septin filaments at the cell cortex. Binding of ReS19-T to septins restores filament assembly in the disease state and restrains calcium entry through SOCCs. In amyloid-β and tau-driven mouse models of disease, ReS19-T agents restored synaptic plasticity, normalized brain network activity, and attenuated the development of both amyloid-β and tau pathology. Our findings identify the septin cytoskeleton as a potential therapeutic target for the development of disease-modifying AD treatments.

Quantitative Prediction of Dissociation Rates of PYK2 Ligands Using Umbrella Sampling and Milestoning

We used umbrella sampling and the milestoning simulation method to study the dissociation of multiple ligands from protein kinase PYK2. The activation barriers obtained from the potential of mean force of the umbrella sampling simulations correlated well with the experimental dissociation rates. Using the zero-temperature string method, we obtained optimized paths along the free-energy surfaces for milestoning simulations of three ligands with a similar structure. The milestoning simulations gave an absolute dissociation rate within 2 orders of magnitude of the experimental value for two ligands but at least 3 orders of magnitude too high for the third. Despite the similarity in their structures, the ligands took different pathways to exit from the binding site of PYK2, making contact with different sets of residues. In addition, the protein experienced different conformational changes for dissociation of the three ligands.

Ligand Gaussian accelerated Molecular Dynamics 3 (LiGaMD3): Improved Calculations of Binding Thermodynamics and Kinetics of Both Small Molecules and Flexible Peptides

Binding thermodynamics and kinetics play critical roles in drug design. However, it has proven challenging to efficiently predict ligand binding thermodynamics and kinetics of small molecules and flexible peptides using conventional Molecular Dynamics (cMD), due to limited simulation timescales. Based on our previously developed Ligand Gaussian accelerated Molecular Dynamics (LiGaMD) method, we present a new approach, termed "LiGaMD3", in which we introduce triple boosts into three individual energy terms that play important roles in small-molecule/peptide dissociation, rebinding and system conformational changes to improve the sampling efficiency of small-molecule/peptide interactions with target proteins. To validate the performance of LiGaMD3, MDM2 bound by a small molecule (Nutlin 3) and two highly flexible peptides (PMI and P53) were chosen as model systems. LiGaMD3 could efficiently capture repetitive small-molecule/peptide dissociation and binding events within 2 microsecond simulations. The predicted binding kinetic constant rates and free energies from LiGaMD3 agreed with available experimental values and previous simulation results. Therefore, LiGaMD3 provides a more general and efficient approach to capture dissociation and binding of both small-molecule ligand and flexible peptides, allowing for accurate prediction of their binding thermodynamics and kinetics.

The rotamer of the second-sphere histidine in AA9 lytic polysaccharide monooxygenase is pH dependent

Lytic polysaccharide monooxygenases (LPMOs) catalyze a reaction that is crucial for the biological decomposition of various biopolymers and for the industrial conversion of plant biomass. Despite the importance of LPMOs, the exact molecular-level nature of the reaction mechanism is still debated today. Here, we investigated the pH-dependent conformation of a second-sphere histidine (His) that we call the stacking histidine, which is conserved in fungal AA9 LPMOs and is speculated to assist catalysis in several of the LPMO reaction pathways. Using constant-pH and accelerated molecular dynamics simulations, we monitored the dynamics of the stacking His in different protonation states for both the resting Cu(II) and active Cu(I) forms of two fungal LPMOs. Consistent with experimental crystallographic and neutron diffraction data, our calculations suggest that the side chain of the protonated and positively charged form is rotated out of the active site toward the solvent. Importantly, only one of the possible neutral states of histidine (HIE state) is observed in the stacking orientation at neutral pH or when bound to cellulose. Our data predict that, in solution, the stacking His may act as a stabilizer (via hydrogen bonding) of the Cu(II)-superoxo complex after the LPMO-Cu(I) has reacted with O2 in solution, which, in fine, leads to H2O2 formation. Also, our data indicate that the HIE-stacking His is a poor acid/base catalyst when bound to the substrate and, in agreement with the literature, may play an important stabilizing role (via hydrogen bonding) during the peroxygenase catalysis. Our study reveals the pH titration midpoint values of the pH-dependent orientation of the stacking His should be considered when modeling and interpreting LPMO reactions, whether it be for classical LPMO kinetics or in industry-oriented enzymatic cocktails, and for understanding LPMO behavior in slightly acidic natural processes such as fungal wood decay.

Computational insights into the binding modes, keto-enol tautomerization and stereo-electronically controlled decarboxylation of oxaloacetate in the active site of macrophomate synthase

Oxaloacetic acid (OAA) is a β-ketocarboxylic acid, which plays an important role as an intermediate in some metabolic pathways, including the tricarboxylic acid cycle, gluconeogenesis and fatty acid biosynthesis. Animal studies have indicated that supplementing oxaloacetic acid shows an increase of lifespan and other substantial health benefits including mitochondrial DNA protection, and protection of retinal, neural and pancreatic tissues. Most of the chemical transformations of OAA in the metabolic pathways have been extensively studied; however, the understanding of decarboxylation of OAA at the atomic level is relatively lacking. Here, we carried out MD simulations and combined quantum mechanical/molecular mechanical (QM/MM) calculations as an example to systematically elucidate the binding modes, keto-enol tautomerization and decarboxylation of OAA in the active site of macrophomate synthase (MPS), which is a Mg(II)-dependent bifunctional enzyme that catalyzes both the decarboxylation of OAA and [4+2] cycloaddition of 2-pyrone with the decarboxylated intermediate of OAA (pyruvate enolate). On the basis of our calculations, it was found that the Mg2+-coordinated oxaloacetate may exist in enol forms and keto forms. The four keto forms can be transformed into each other by simply rotating the C2-C3 single bond, nevertheless, the keto-enol tautomerization strictly requires the assistance of pocket water molecules. In addition, the decarboxylation is stereo-electronically controlled, i.e., it is the relative orientation of the terminal carboxyl anion that determines the rate of decarboxylation. As such, the chemistry of oxaloacetate in the active site of MPS is complex. On one hand, the most stable binding mode (K-I) may undergo enol-keto tautomerization to isomerize to the enol form, which may further react with the second substrate; on the other hand, K-I may isomerize to another binding mode K-II to proceed decarboxylation to generate pyruvate enolate and CO2. Starting from K-I, the enol-keto tautomerization corresponds to a barrier of 16.2 kcal mol-1, whereas the decarboxylation is associated with an overall barrier of 19.7 kcal mol-1. These findings may provide useful information for understanding the chemistry of OAA and the catalysis of related enzymes, and they are basically in agreement with the available experimental kinetic data.

Exploring Product Release from Yeast Cytosine Deaminase with Metadynamics

The yeast cytosine deaminase (yCD) enzyme/5-fluorocytosine prodrug system is a promising candidate for targeted chemotherapeutics. After conversion of the prodrug into the toxic chemotherapeutic drug, 5-fluorouracil (5-FU), the slow product release from the enzyme limits the overall catalytic efficiency of the enzyme/prodrug system. Here, we present a computational study of the product release of the anticancer drug, 5-FU, from yCD using metadynamics. We present a comparison of the 5-FU drug to the natural enzyme product, uracil. We use volume-based metadynamics to compute the free energy landscape for product release and show a modest binding affinity for the product to the enzyme, consistent with experiments. Next, we use infrequent metadynamics to estimate the unbiased release rate from Kramers time-dependent rate theory and find a favorable comparison to experiment with a slower rate of product release for the 5-FU system. Our work demonstrates how adaptive sampling methods can be used to study the protein-ligand unbinding process for engineering enzyme/prodrug systems and gives insights into the molecular mechanism of product release for the yCD/5-FU system.

Modeling of noncovalent inhibitors of the papain-like protease (PLpro) from SARS-CoV-2 considering the protein flexibility by using molecular dynamics and cross-docking

The papain-like protease (PLpro) found in coronaviruses that can be transmitted from animals to humans is a critical target in respiratory diseases linked to Severe Acute Respiratory Syndrome (SARS-CoV). Researchers have proposed designing PLpro inhibitors. In this study, a set of 89 compounds, including recently reported 2-phenylthiophenes with nanomolar inhibitory potency, were investigated as PLpro noncovalent inhibitors using advanced molecular modeling techniques. To develop the work with these inhibitors, multiple structures of the SARS-CoV-2 PLpro binding site were generated using a molecular sampling method. These structures were then clustered to select a group that represents the flexibility of the site. Subsequently, models of the protein-ligand complexes were created for the set of inhibitors within the chosen conformations. The quality of the complex models was assessed using LigRMSD software to verify similarities in the orientations of the congeneric series and interaction fingerprints to determine the recurrence of chemical interactions. With the multiple models constructed, a protocol was established to choose one per ligand, optimizing the correlation between the calculated docking energy values and the biological activities while incorporating the effect of the binding site's flexibility. A strong correlation (R2 = 0.922) was found when employing this flexible docking protocol.

Spectroscopy and molecular simulation on the interaction of Nano-Kaempferol prepared by oil-in-water with two carrier proteins: An investigation of protein-protein interaction

In this work, the interaction of human serum albumin (HSA) and human holo-transferrin (HTF) with the prepared Nano-Kaempferol (Nano-KMP) through oil-in-water procedure was investigated in the form of binary and ternary systems by the utilization of different spectroscopy techniques along with molecular simulation and cancer cell experiments. According to fluorescence spectroscopy outcomes, Nano-KMP is capable of quenching both proteins as binary systems by a static mechanism, while in the form of (HSA-HTF) Nano-KMP as the ternary system, an unlinear Stern-Volmer plot was elucidated with the occurrence of both dynamic and static fluorescence quenching mechanisms in the binding interaction. In addition, the two acquired Ksv values in the ternary system signified the existence of two sets of binding sites with two different interaction behaviors. The binding constant values of HSA-Nano KMP, HTF-Nano-KMP, and (HSA-HTF) Nano-KMP complexes formation were (2.54 ± 0.03) × 104, (2.15 ± 0.02) × 104 and (1.43 ± 0.04) × 104M-1at the first set of binding sites and (4.68 ± 0.05) × 104 M-1 at the second set of binding sites, respectively. The data of thermodynamic parameters confirmed the major roles of hydrogen binding and van der Waals forces in the formation of HSA-Nano KMP and HTF-Nano KMP complexes. The thermodynamic parameter values of (HSA-HTF) Nano KMP revealed the dominance of hydrogen binding and van der Waals forces in the first set of binding sites and hydrophobic forces for the second set of binding sites. Resonance light scattering (RLS) analysis displayed the existence of a different interaction behavior for HSA-HTF complex in the presence of Nano-KMP as the ternary system. Moreover, circular dichroism (CD) technique affirmed the conformational changes of the secondary structure of proteins as binary and ternary systems. Molecular docking and molecular dynamics simulations (for 100 ns) were performed to investigate the mechanism of KMP binding to HSA, HTF, and HSA-HTF. Next to observing a concentration and time-dependent cytotoxicity, the down regulation of PI3K/AkT/mTOR pathway resulted in cell cycle arrest in SW480 cells.

Computational Chemistry in Structure-Based Solute Carrier Transporter Drug Design: Recent Advances and Future Perspectives

Solute carrier transporters (SLCs) are a class of important transmembrane proteins that are involved in the transportation of diverse solute ions and small molecules into cells. There are approximately 450 SLCs within the human body, and more than a quarter of them are emerging as attractive therapeutic targets for multiple complex diseases, e.g., depression, cancer, and diabetes. However, only 44 unique transporters (∼9.8% of the SLC superfamily) with 3D structures and specific binding sites have been reported. To design innovative and effective drugs targeting diverse SLCs, there are a number of obstacles that need to be overcome. However, computational chemistry, including physics-based molecular modeling and machine learning- and deep learning-based artificial intelligence (AI), provides an alternative and complementary way to the classical drug discovery approach. Here, we present a comprehensive overview on recent advances and existing challenges of the computational techniques in structure-based drug design of SLCs from three main aspects: (i) characterizing multiple conformations of the proteins during the functional process of transportation, (ii) identifying druggability sites especially the cryptic allosteric ones on the transporters for substrates and drugs binding, and (iii) discovering diverse small molecules or synthetic protein binders targeting the binding sites. This work is expected to provide guidelines for a deep understanding of the structure and function of the SLC superfamily to facilitate rational design of novel modulators of the transporters with the aid of state-of-the-art computational chemistry technologies including artificial intelligence.

Molecular dynamics-based identification of binding pathways and two distinct high-affinity sites for succinate in succinate receptor 1/GPR91

SUCNR1 is an auto- and paracrine sensor of the metabolic stress signal succinate. Using unsupervised molecular dynamics (MD) simulations (170.400 ns) and mutagenesis across human, mouse, and rat SUCNR1, we characterize how a five-arginine motif around the extracellular pole of TM-VI determines the initial capture of succinate in the extracellular vestibule (ECV) to either stay or move down to the orthosteric site. Metadynamics demonstrate low-energy succinate binding in both sites, with an energy barrier corresponding to an intermediate stage during which succinate, with an associated water cluster, unlocks the hydrogen-bond-stabilized conformationally constrained extracellular loop (ECL)-2b. Importantly, simultaneous binding of two succinate molecules through either a "sequential" or "bypassing" mode is a frequent endpoint. The mono-carboxylate NF-56-EJ40 antagonist enters SUCNR1 between TM-I and -II and does not unlock ECL-2b. It is proposed that occupancy of both high-affinity sites is required for selective activation of SUCNR1 by high local succinate concentrations.

Multiscale Study on the Intramolecular C-S Bond Formation Catalyzed by P450 Monooxygenase CxnD Involved in the Biosynthesis of Chuangxinmycin: The Critical Roles of Noncrystal Water Molecule and Conformational Change

Cytochrome P450 monooxygenase CxnD catalyzes intramolecular C-S bond formation in the biosynthesis of chuangxinmycin, which is representative of the synthesis of sulfur-containing natural heterocyclic compounds. The intramolecular cyclization usually requires the activation of two reaction sites and a large conformational change; thus, illuminating its detailed reaction mechanism remains challengeable. Here, the reaction pathway of CxnD-catalyzed C-S bond formation was clarified by a series of calculations, including Gaussian accelerated molecular dynamics simulations and quantum mechanical-molecular mechanical calculations. Our results revealed that the C-S formation follows a diradical coupling mechanism. CxnD first employs Cpd I to abstract the hydrogen atom from the imino group of the indole ring, and then, the resulted Cpd II further extracts another hydrogen atom from the thiol group of the side chain to afford a diradical intermediate, in which a noncrystal water molecule entering into the active site after the formation of Cpd I was proved to play an indispensable role. Moreover, the diradical intermediate cannot directly perform the coupling reaction. It should first undergo a series of conformational changes leading to the proximity of two reaction sites. It is the flexibility of the active site of the enzyme and the side chain of the substrate that makes the diradical coupling to be successful.

Mechanistic insights into ligand dissociation from the SARS-CoV-2 spike glycoprotein

The COVID-19 pandemic, driven by the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), has spurred an urgent need for effective therapeutic interventions. The spike glycoprotein of the SARS-CoV-2 is crucial for infiltrating host cells, rendering it a key candidate for drug development. By interacting with the human angiotensin-converting enzyme 2 (ACE2) receptor, the spike initiates the infection of SARS-CoV-2. Linoleate is known to bind the spike glycoprotein, subsequently reducing its interaction with ACE2. However, the detailed mechanisms underlying the protein-ligand interaction remain unclear. In this study, we characterized the pathways of ligand dissociation and the conformational changes associated with the spike glycoprotein by using ligand Gaussian accelerated molecular dynamics (LiGaMD). Our simulations resulted in eight complete ligand dissociation trajectories, unveiling two distinct ligand unbinding pathways. The preference between these two pathways depends on the gate distance between two α-helices in the receptor binding domain (RBD) and the position of the N-linked glycan at N343. Our study also highlights the essential contributions of K417, N121 glycan, and N165 glycan in ligand unbinding, which are equally crucial in enhancing spike-ACE2 binding. We suggest that the presence of the ligand influences the motions of these residues and glycans, consequently reducing accessibility for spike-ACE2 binding. These findings enhance our understanding of ligand dissociation from the spike glycoprotein and offer significant implications for drug design strategies in the battle against COVID-19.

Enhanced Sampling of Buried Charges in Free Energy Calculations Using Replica Exchange with Charge Tempering

Buried ionizable groups in proteins often play important structural and functional roles. However, it is generally challenging to study the detailed molecular mechanisms solely based on experimental measurements. Free energy calculations using atomistic simulations, on the other hand, complement experimental studies and can provide high temporal and spatial resolution information that can lead to mechanistic insights. Nevertheless, it is also well recognized that sufficient sampling of such atomistic simulations can be challenging, considering that structural changes related to the buried charges may be very slow. In the present study, we describe a simple but effective enhanced sampling technique called replica exchange with charge tempering (REChgT) with a novel free energy method, multisite λ dynamics (MSλD), to study two systems containing buried charges, pKa prediction of a small molecule, orotate, in complex with the dihydroorotate dehydrogenase, and relative stability of a Glu-Lys pair buried in the hydrophobic core of two variants of Staphylococcal nuclease. Compared to the original MSλD simulations, the usage of REChgT dramatically increases sampling in both conformational and alchemical spaces, which directly translates into a significant reduction of wall time to converge the free energy calculations. This study highlights the importance of sufficient sampling toward developing improved free energy methods.

Mechanism of Ligand Binding to Theophylline RNA Aptamer

Studying RNA-ligand interactions and quantifying their binding thermodynamics and kinetics are of particular relevance in the field of drug discovery. Here, we combined biochemical binding assays and accelerated molecular simulations to investigate ligand binding and dissociation in RNA using the theophylline-binding RNA as a model system. All-atom simulations using a Ligand Gaussian accelerated Molecular Dynamics method (LiGaMD) have captured repetitive binding and dissociation of theophylline and caffeine to RNA. Theophylline's binding free energy and kinetic rate constants align with our experimental data, while caffeine's binding affinity is over 10,000 times weaker, and its kinetics could not be determined. LiGaMD simulations allowed us to identify distinct low-energy conformations and multiple ligand binding pathways to RNA. Simulations revealed a "conformational selection" mechanism for ligand binding to the flexible RNA aptamer, which provides important mechanistic insights into ligand binding to the theophylline-binding model. Our findings suggest that compound docking using a structural ensemble of representative RNA conformations would be necessary for structure-based drug design of flexible RNA.

Accelerated Molecular Dynamics and AlphaFold Uncover a Missing Conformational State of Transporter Protein OxlT

Transporter proteins change their conformations to carry their substrate across the cell membrane. The conformational dynamics is vital to understanding the transport function. We have studied the oxalate transporter (OxlT), an oxalate:formate antiporter from Oxalobacter formigenes, significant in avoiding kidney stone formation. The atomic structure of OxlT has been recently solved in the outward-open and occluded states. However, the inward-open conformation is still missing, hindering a complete understanding of the transporter. Here, we performed a Gaussian accelerated molecular dynamics simulation to sample the extensive conformational space of OxlT and successfully predicted the inward-open conformation where cytoplasmic substrate formate binding was preferred over oxalate binding. We also identified critical interactions for the inward-open conformation. The results were complemented by an AlphaFold2 structure prediction. Although AlphaFold2 solely predicted OxlT in the outward-open conformation, mutation of the identified critical residues made it partly predict the inward-open conformation, identifying possible state-shifting mutations.

Accelerating Molecular Dynamics Simulations for Drug Discovery

Accurate prediction of ligand binding thermodynamics and kinetics is crucial in drug design. However, it remains challenging for conventional molecular dynamics (MD) simulations due to sampling issues. Gaussian accelerated MD (GaMD) is an enhanced sampling method that adds a harmonic boost to overcome energy barriers, which has demonstrated significant benefits in exploring protein-ligand interactions. Especially, the ligand GaMD (LiGaMD) applies a selective boost potential to the ligand nonbonded potential energy, significantly improving sampling for ligand binding and dissociation. Furthermore, a selective boost potential is applied to the potential of both ligand and protein residues around binding pocket in LiGaMD2 to further increase the sampling of protein-ligand interaction. LiGaMD and LiGaMD2 simulations could capture repetitive ligand binding and unbinding events within microsecond simulations, allowing to simultaneously characterize ligand binding thermodynamics and kinetics, which is expected to greatly facilitate drug design. In this chapter, we provide a brief review of the status of LiGaMD in drug discovery and outline its usage.

Roles of Accelerated Molecular Dynamics Simulations in Predictions of Binding Kinetic Parameters

Rational predictions on binding kinetics parameters of drugs to targets play significant roles in future drug designs. Full conformational samplings of targets are requisite for accurate predictions of binding kinetic parameters. In this review, we mainly focus on the applications of enhanced sampling technologies in calculations of binding kinetics parameters and residence time of drugs. The methods involved in molecular dynamics simulations are applied to not only probe conformational changes of targets but also reveal calculations of residence time that is significant for drug efficiency. For this review, special attention are paid to accelerated molecular dynamics (aMD) and Gaussian aMD (GaMD) simulations that have been adopted to predict the association or disassociation rate constant. We also expect that this review can provide useful information for future drug design.