Unbinding Kinetics of a p38 MAP Kinase Type II Inhibitor from Metadynamics Simulations

Good Rates From Bad Coordinates: The Exponential Average Time-dependent Rate Approach

Our ability to calculate rate constants of biochemical processes using molecular dynamics simulations is severely limited by the fact that the time scales for reactions, or changes in conformational state, scale exponentially with the relevant free-energy barrier heights. In this work, we improve upon a recently proposed rate estimator that allows us to predict transition times with molecular dynamics simulations biased to rapidly explore one or several collective variables (CVs). This approach relies on the idea that not all bias goes into promoting transitions, and along with the rate, it estimates a concomitant scale factor for the bias termed the "CV biasing efficiency" γ. First, we demonstrate mathematically that our new formulation allows us to derive the commonly used Infrequent Metadynamics (iMetaD) estimator when using a perfect CV, where γ = 1. After testing it on a model potential, we then study the unfolding behavior of a previously well characterized coarse-grained protein, which is sufficiently complex that we can choose many different CVs to bias, but which is sufficiently simple that we are able to compute the unbiased rate directly. For this system, we demonstrate that predictions from our new Exponential Average Time-Dependent Rate (EATR) estimator converge to the true rate constant more rapidly as a function of bias deposition time than does the previous iMetaD approach, even for bias deposition times that are short. We also show that the γ parameter can serve as a good metric for assessing the quality of the biasing coordinate. We demonstrate that these results hold when applying the methods to an atomistic protein folding example. Finally, we demonstrate that our approach works when combining multiple less-than-optimal bias coordinates, and adapt our method to the related "OPES flooding" approach. Overall, our time-dependent rate approach offers a powerful framework for predicting rate constants from biased simulations.

How exascale computing can shape drug design: A perspective from multiscale QM/MM molecular dynamics simulations and machine learning-aided enhanced sampling algorithms

Molecular simulations are an essential asset in the first steps of drug design campaigns. However, the requirement of high-throughput limits applications mainly to qualitative approaches with low computational cost, but also low accuracy. Unlocking the potential of more rigorous quantum mechanical/molecular mechanics (QM/MM) models combined with molecular dynamics-based free energy techniques could have a tremendous impact. Indeed, these two relatively old techniques are emerging as promising methods in the field. This has been favored by the exponential growth of computer power and the proliferation of powerful data-driven methods. Here, we briefly review recent advances and applications, and give our perspective on the impact that QM/MM and free-energy methods combined with machine learning-aided algorithms can have on drug design.

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.

Accurate Characterization of Binding Kinetics and Allosteric Mechanisms for the HSP90 Chaperone Inhibitors Using AI-Augmented Integrative Biophysical Studies

The binding kinetics of drugs to their targets are gradually being recognized as a crucial indicator of the efficacy of drugs in vivo, leading to the development of various computational methods for predicting the binding kinetics in recent years. However, compared with the prediction of binding affinity, the underlying structure and dynamic determinants of binding kinetics are more complicated. Efficient and accurate methods for predicting binding kinetics are still lacking. In this study, quantitative structure-kinetics relationship (QSKR) models were developed using 132 inhibitors targeting the ATP binding domain of heat shock protein 90α (HSP90α) to predict the dissociation rate constant (koff), enabling a direct assessment of the drug-target residence time. These models demonstrated good predictive performance, where hydrophobic and hydrogen bond interactions significantly influence the koff prediction. In subsequent applications, our models were used to assist in the discovery of new inhibitors for the N-terminal domain of HSP90α (N-HSP90α), demonstrating predictive capabilities on an experimental validation set with a new scaffold. In X-ray crystallography experiments, the loop-middle conformation of apo N-HSP90α was observed for the first time (previously, the loop-middle conformation had only been observed in holo-N-HSP90α structures). Interestingly, we observed different conformations of apo N-HSP90α simultaneously in an asymmetric unit, which was also observed in a holo-N-HSP90α structure, suggesting an equilibrium of conformations between different states in solution, which could be one of the determinants affecting the binding kinetics of the ligand. Different ligands can undergo conformational selection or alter the equilibrium of conformations, inducing conformational rearrangements and resulting in different effects on binding kinetics. We then used molecular dynamics simulations to describe conformational changes of apo N-HSP90α in different conformational states. In summary, the study of the binding kinetics and molecular mechanisms of N-HSP90α provides valuable information for the development of more targeted therapeutic approaches.

Exploring Kinase Asp-Phe-Gly (DFG) Loop Conformational Stability with AlphaFold2-RAVE

Kinases compose one of the largest fractions of the human proteome, and their misfunction is implicated in many diseases, in particular, cancers. The ubiquitousness and structural similarities of kinases make specific and effective drug design difficult. In particular, conformational variability due to the evolutionarily conserved Asp-Phe-Gly (DFG) motif adopting in and out conformations and the relative stabilities thereof are key in structure-based drug design for ATP competitive drugs. These relative conformational stabilities are extremely sensitive to small changes in sequence and provide an important problem for sampling method development. Since the invention of AlphaFold2, the world of structure-based drug design has noticeably changed. In spite of it being limited to crystal-like structure prediction, several methods have also leveraged its underlying architecture to improve dynamics and enhanced sampling of conformational ensembles, including AlphaFold2-RAVE. Here, we extend AlphaFold2-RAVE and apply it to a set of kinases: the wild type DDR1 sequence and three mutants with single point mutations that are known to behave drastically differently. We show that AlphaFold2-RAVE is able to efficiently recover the changes in relative stability using transferable learned order parameters and potentials, thereby supplementing AlphaFold2 as a tool for exploration of Boltzmann-weighted protein conformations (Meller, A.; Bhakat, S.; Solieva, S.; Bowman, G. R. Accelerating Cryptic Pocket Discovery Using AlphaFold. J. Chem. Theory Comput. 2023, 19, 4355-4363).

Applications of Molecular Dynamics Simulations in Drug Discovery

In the current drug development process, molecular dynamics (MD) simulations have proven to be very useful. This chapter provides an overview of the current applications of MD simulations in drug discovery, from detecting protein druggable sites and validating drug docking outcomes to exploring protein conformations and investigating the influence of mutations on its structure and functions. In addition, this chapter emphasizes various strategies to improve the conformational sampling efficiency in molecular dynamics simulations. With a growing computer power and developments in the production of force fields and MD techniques, the importance of MD simulations in helping the drug development process is projected to rise significantly in the future.

Investigation of Rare Protein Conformational Transitions via Dissipation-Corrected Targeted Molecular Dynamics

To sample rare events, dissipation-corrected targeted molecular dynamics (dcTMD) applies a constant velocity constraint along a one-dimensional reaction coordinate s, which drives an atomistic system from an initial state into a target state. Employing a cumulant approximation of Jarzynski's identity, the free energy ΔG(s) is calculated from the mean external work and dissipated work of the process. By calculating the friction coefficient Γ(s) from the dissipated work, in a second step, the equilibrium dynamics of the process can be studied by propagating a Langevin equation. While so far dcTMD has been mostly applied to study the unbinding of protein-ligand complexes, here its applicability to rare conformational transitions within a protein and the prediction of their kinetics are investigated. As this typically requires the introduction of multiple collective variables {xj} = x, a theoretical framework is outlined to calculate the associated free energy ΔG(x) and friction Γ(x) from dcTMD simulations along coordinate s. Adopting the α-β transition of alanine dipeptide as well as the open-closed transition of T4 lysozyme as representative examples, the virtues and shortcomings of dcTMD to predict protein conformational transitions and the related kinetics are studied.

Assessing Metadynamics and Docking for Absolute Binding Free Energy Calculations Using Severe Acute Respiratory Syndrome Coronavirus 2 Main Protease Inhibitors

Absolute binding free energy (ABFE) calculations can be an important part of the drug discovery process by identifying molecules that have the potential to be strong binders for a biomolecular target. Recent work has used free energy perturbation (FEP) theory for these calculations, focusing on a set of 16 inhibitors of the severe acute respiratory syndrome coronavirus 2 main protease (Mpro). Herein, the same data set is evaluated by metadynamics (MetaD), four different docking programs, and molecular mechanics with generalized Born and surface area solvation. MetaD yields a Kendall τ distance of 0.28 and Pearson r2 of 0.49, which reflect somewhat less accuracy than that from the ABFE FEP results. Notably, it is demonstrated that an ensemble docking protocol by which each ligand is docked into the 13 crystal structures in this data set provides improved performance, particularly when docking is carried out with Glide XP (Kendall τ distance = 0.20, Pearson r2 = 0.71), Glide SP (Kendall τ distance = 0.19, Pearson r2 = 0.66), or AutoDock 4 (Kendall τ distance = 0.21, Pearson r2 = 0.55). The best results are obtained with "superconsensus" docking by averaging the 52 results for each compound using the 4 docking protocols and all 13 crystal structures (Kendall τ distance = 0.18, Pearson r2 = 0.73).

Multifunctional Roles of Clathrate Hydrate Nanoreactors for CO2 Reduction

In this study, we explore a possible platform for the CO2 reduction (CO2 R) in one of water's solid phases, namely clathrate hydrates (CHs), by ab initio molecular dynamics and well-tempered metadynamics simulations with periodic boundary conditions. We found that the stacked H2 O nanocages in CHs help to initialize CO2 R by increasing the electron-binding ability of CO2 . The substantial CO2 R processes are further influenced by the hydrogen bond networks in CHs. The first intermediate CO2 - in this process can be stabilized through cage structure reorganization into the H-bonded [CO2 - ⋅⋅⋅H-OHcage ] complex. Further cooperative structural dynamics enables the complex to convert into a vital transient [CO2 2- ⋅⋅⋅H-OHcage ] intermediate in a low-barrier disproportionation-like process. Such a highly reactive intermediate spontaneously triggers subsequent double proton transfer along its tethering H-bonds, finally converting it into HCOOH. These hydrogen-bonded nanoreactors feature multiple functions in facilitating CO2 R such as confining, tethering, H-bond catalyzing and proton pumping. Our findings have a general interest and extend the knowledge of CO2 R into porous aqueous systems.

Free-Energy Landscape and Rate Estimation of the Aromatic Ring Flips in Basic Pancreatic Trypsin Inhibitors Using Metadynamics

Aromatic side chains (phenylalanine and tyrosine) of a protein flip by 180° around the Cβ-Cγ axis (χ2 dihedral of the side chain), producing two symmetry-equivalent states. The study of ring flip dynamics with nuclear magnetic resonance (NMR) experiments helps to understand local conformational fluctuations. Ring flips are categorized as slow (milliseconds and onward) or fast (nanoseconds to near milliseconds) based on timescales accessible to NMR experiments. In this study, we investigated the ability of the infrequent metadynamics approach to estimate the flip rate and discriminate between slow and fast ring flips for eight individual aromatic side chains (F4, Y10, Y21, F22, Y23, F33, Y35, and F45) of the basic pancreatic trypsin inhibitor. Well-tempered metadynamics simulations were performed to estimate the ring-flipping free-energy surfaces for all eight aromatic residues. The results indicate that χ2 as a standalone collective variable (CV) is not sufficient to obtain computationally consistent results. Inclusion of a complementary CV, such as χ1(Cα-Cβ), solved the problem for most residues and enabled us to classify fast and slow ring flips. This indicates the importance of librational motions in ring flips. Multiple pathways and mechanisms were observed for residues F4, Y10, and F22. Recrossing events were observed for residues F22 and F33, indicating a possible role of friction effects in ring flipping. The results demonstrate the successful application of infrequent metadynamics to estimate ring flip rates and identify certain limitations of the approach.

Intrinsic structural dynamics dictate enzymatic activity and inhibition

Enzymes are known to sample various conformations, many of which are critical for their biological function. However, structural characterizations of enzymes predominantly focus on the most populated conformation. As a result, single-point mutations often produce structures that are similar or essentially identical to those of the wild-type enzyme despite large changes in enzymatic activity. Here, we show for mutants of a histone deacetylase enzyme (HDAC8) that reduced enzymatic activities, reduced inhibitor affinities, and reduced residence times are all captured by the rate constants between intrinsically sampled conformations that, in turn, can be obtained independently by solution NMR spectroscopy. Thus, for the HDAC8 enzyme, the dynamic sampling of conformations dictates both enzymatic activity and inhibitor potency. Our analysis also dissects the functional role of the conformations sampled, where specific conformations distinct from those in available structures are responsible for substrate and inhibitor binding, catalysis, and product dissociation. Precise structures alone often do not adequately explain the effect of missense mutations on enzymatic activity and drug potency. Our findings not only assign functional roles to several conformational states of HDAC8 but they also underscore the paramount role of dynamics, which will have general implications for characterizing missense mutations and designing inhibitors.

Perspectives on Ligand/Protein Binding Kinetics Simulations: Force Fields, Machine Learning, Sampling, and User-Friendliness

Computational techniques applied to drug discovery have gained considerable popularity for their ability to filter potentially active drugs from inactive ones, reducing the time scale and costs of preclinical investigations. The main focus of these studies has historically been the search for compounds endowed with high affinity for a specific molecular target to ensure the formation of stable and long-lasting complexes. Recent evidence has also correlated the in vivo drug efficacy with its binding kinetics, thus opening new fascinating scenarios for ligand/protein binding kinetic simulations in drug discovery. The present article examines the state of the art in the field, providing a brief summary of the most popular and advanced ligand/protein binding kinetics techniques and evaluating their current limitations and the potential solutions to reach more accurate kinetic models. Particular emphasis is put on the need for a paradigm change in the present methodologies toward ligand and protein parametrization, the force field problem, characterization of the transition states, the sampling issue, and algorithms' performance, user-friendliness, and data openness.

Kinetics from Metadynamics: Principles, Applications, and Outlook

Metadynamics is a popular enhanced sampling algorithm for computing the free energy landscape of rare events by using molecular dynamics simulation. Ten years ago, Tiwary and Parrinello introduced the infrequent metadynamics approach for calculating the kinetics of transitions across free energy barriers. Since then, metadynamics-based methods for obtaining rate constants have attracted significant attention in computational molecular science. Such methods have been applied to study a wide range of problems, including protein-ligand binding, protein folding, conformational transitions, chemical reactions, catalysis, and nucleation. Here, we review the principles of elucidating kinetics from metadynamics-like approaches, subsequent methodological developments in this area, and successful applications on chemical, biological, and material systems. We also highlight the challenges of reconstructing accurate kinetics from enhanced sampling simulations and the scope of future developments.

Temperature Accelerated Sliced Sampling to Probe Ligand Dissociation from Protein

Modeling ligand unbinding in proteins to estimate the free energy of binding and probing the mechanism presents several challenges. They primarily pertain to the entropic bottlenecks resulting from protein and solvent conformations. While exploring the unbinding processes using enhanced sampling techniques, very long simulations are required to sample all of the conformational states as the system gets trapped in local free energy minima along transverse coordinates. Here, we demonstrate that temperature accelerated sliced sampling (TASS) is an ideal approach to overcome some of the difficulties faced by conventional sampling methods in studying ligand unbinding. Using TASS, we study the unbinding of avibactam inhibitor molecules from the Class C β-lactamase (CBL) active site. Extracting CBL-avibactam unbinding free energetics, unbinding pathways, and identifying critical interactions from the TASS simulations are demonstrated.

On the Dynamic Mechanism of Long-Flexible Fatty Acid Binding to Fatty Acid Binding Protein: Resolving the Long-Standing Debate

Fatty acids (FAs) are one of the essential energy sources for physiological processes, and they play a vital role in regulating immune and inflammatory responses, promoting cell differentiation and apoptosis, and inhibiting tumor growth. These functions are carried out by FA binding proteins (FABPs) that recognize and transport FAs. Although the crystal structure of the FA-FABPs complex has long been characterized, the mechanism behind FA binding and dissociation from FABP remains unclear. This study employed conventional MD simulations and enhanced sampling technologies to investigate the atomic-scale complexes of heart fatty acid binding proteins and stearic acid (SA). The results revealed two primary pathways for the binding or dissociation of the flexible long-chain ligand, with the orientation of the SA carboxyl head during dissociation determining the chosen path. Conformational changes in the portal region of FABP during the ligand binding/unbinding were found to be trivial, and the overturn of the ″cap″ or the unfolding of the α2 helix was not required. This study resolves the long-standing debate on the binding mechanism of SA with the long-flexible tail to FABP, which significantly improves the understanding of the transport mechanism of FABPs and the development of related therapeutic agents.

Dynamic structural analysis-based epitope prediction of Exendin-4 in aqueous solution

The study of epitopes has a broad range of applications in drug discovery, vaccine design, and immunotherapy. In this study, an epitope prediction method was developed based on the dynamic structure of protein antigens. Solvent accessible surface area, charge, and root mean square fluctuation were introduced as the key residue property parameters. The epitope prediction algorithm was established by constructing a three-parameter complex metrics of seven-peptide groups. The method was applied to predict the epitopes of Exendin-4, an effective antidiabetic drug. The epitopes of both the natural and C-terminal amidated forms of Exendin-4 were predicted and compared in their folded and intermediate states. In the folded state, the epitopes of natural Exendin-4 (His1-Phe6 and Asp9-Val19) were found to be nearly identical to the epitopes of C-terminal aminated Exendin-4 (His1-Thr7 and Asp9-Val19). In the intermediate state, however, the epitopes of natural Exendin-4 (His1-Gly4, Phe6 and Lys12-Arg20) covered fewer amino acids than the epitopes of C-terminal aminated Exendin-4 (His1-Gly4, Phe6, Asp9-Val19 and Trp25-Lys27). The comparison with the results from other prediction tools demonstrates the reliability of our predicted epitopes of Exendin-4.

15 Years of molecular simulation of drug-binding kinetics

INTRODUCTION

Drug-binding kinetics has been increasingly recognized as an important factor to be considered in drug discovery. Long residence time could prolong the action of some drugs while produce toxicity on others. Early evaluation of the binding kinetics of drug candidates could reduce attrition rate late in the drug discovery process. Computational prediction of drug-binding kinetics is useful as compounds can be evaluated even before they are made. However, simulation of drug-binding kinetics is a challenging problem because of the long-time scale involved. Nevertheless, significant progress has been made.

AREAS COVERED

This review illustrates the rapid evolution of qualitative to quantitative molecular dynamics-based methods that have been developed over the last 15 years.

EXPERT OPINION

The development of new methods based on molecular dynamics simulations now enables computation of absolute association/dissociation rate constants. Cheaper methods capable of identifying candidates with fast or slow binding kinetics, or rank-ordering rate constants are also available. Together, these methods have generated useful insights into the molecular mechanisms of drug-binding kinetics, and the design of drug candidates with therapeutically favorable kinetics. Although predicting absolute rate constants is still expensive and challenging, rapid improvement is expected in the coming years with the continuing refinement of current technologies, development of new methodologies, and the utilization of machine learning.

Structure-Kinetic Relationship for Drug Design Revealed by a PLS Model with Retrosynthesis-Based Pre-Trained Molecular Representation and Molecular Dynamics Simulation

Drug design based on kinetic properties is growing in application. Here, we applied retrosynthesis-based pre-trained molecular representation (RPM) in machine learning (ML) to train 501 inhibitors of 55 proteins and successfully predicted the dissociation rate constant (k_off) values of 38 inhibitors from an independent dataset for the N-terminal domain of heat shock protein 90α (N-HSP90). Our RPM molecular representation outperforms other pre-trained molecular representations such as GEM, MPG, and general molecular descriptors from RDKit. Furthermore, we optimized the accelerated molecular dynamics to calculate the relative retention time (RT) for the 128 inhibitors of N-HSP90 and obtained the protein-ligand interaction fingerprints (IFPs) on their dissociation pathways and their influencing weights on the k_off value. We observed a high correlation among the simulated, predicted, and experimental -log(k_off) values. Combining ML, molecular dynamics (MD) simulation, and IFPs derived from accelerated MD helps design a drug for specific kinetic properties and selectivity profiles to the target of interest. To further validate our k_off predictive ML model, we tested our model on two new N-HSP90 inhibitors, which have experimental k_off values and are not in our ML training dataset. The predicted k_off values are consistent with experimental data, and the mechanism of their kinetic properties can be explained by IFPs, which shed light on the nature of their selectivity against N-HSP90 protein. We believe that the ML model described here is transferable to predict k_off of other proteins and will enhance the kinetics-based drug design endeavor.

Metadynamics simulations for the investigation of drug loading on functionalized inorganic nanoparticles

Inorganic nanoparticles show promising properties that allow them to be efficiently used as drug carriers. The main limitation in this type of application is currently the drug loading capacity, which can be overcome with a proper functionalization of the nanoparticle surface. In this study, we present, for the first time, a computational approach based on metadynamics to estimate the binding free energy of the doxorubicin drug (DOX) to a functionalized TiO₂ nanoparticle under different pH conditions. On a thermodynamic basis, we demonstrate the robustness of our approach to capture the overall mechanism behind the pH-triggered release of DOX due to environmental pH changes. Notably, binding free energy estimations align well with what is expected for a pH-sensitive drug delivery system. Based on our results, we envision the use of metadynamics as a promising computational tool for the rational design and in silico optimization of organic ligands with improved drug carrier properties.

Unbinding Kinetics of Muscarinic M3 Receptor Antagonists Explained by Metadynamics Simulations

The residence time (RT), the time for which a drug remains bound to its biological target, is a critical parameter for drug design. The prediction of this key kinetic property has been proven to be challenging and computationally demanding in the framework of atomistic simulations. In the present work, we setup and applied two distinct metadynamics protocols to estimate the RTs of muscarinic M3 receptor antagonists. In the first method, derived from the conformational flooding approach, the kinetics of unbinding is retrieved from a physics-based parameter known as the acceleration factor α (i.e., the running average over time of the potential deposited in the bound state). Such an approach is expected to recover the absolute RT value for a compound of interest. In the second method, known as the t_META-D approach, a qualitative estimation of the RT is given by the time of simulation required to drive the ligand from the binding site to the solvent bulk. This approach has been developed to reproduce the change of experimental RTs for compounds targeting the same target. Our analysis shows that both computational protocols are able to rank compounds in agreement with their experimental RTs. Quantitative structure-kinetics relationship (SKR) models can be identified and employed to predict the impact of a chemical modification on the experimental RT once a calibration study has been performed.

Girsanov Reweighting Enhanced Sampling Technique (GREST): On-the-Fly Data-Driven Discovery of and Enhanced Sampling in Slow Collective Variables

Molecular dynamics simulations of microscopic phenomena are limited by the short integration time steps which are required for numerical stability but which limit the practically achievable simulation time scales. Collective variable (CV) enhanced sampling techniques apply biases to predefined collective coordinates to promote barrier crossing, phase space exploration, and sampling of rare events. The efficacy of these techniques is contingent on the selection of good CVs correlated with the molecular motions governing the long-time dynamical evolution of the system. In this work, we introduce Girsanov Reweighting Enhanced Sampling Technique (GREST) as an adaptive sampling scheme that interleaves rounds of data-driven slow CV discovery and enhanced sampling along these coordinates. Since slow CVs are inherently dynamical quantities, a key ingredient in our approach is the use of both thermodynamic and dynamical Girsanov reweighting corrections for rigorous estimation of slow CVs from biased simulation data. We demonstrate our approach on a toy 1D 4-well potential, a simple biomolecular system alanine dipeptide, and the Trp-Leu-Ala-Leu-Leu (WLALL) pentapeptide. In each case GREST learns appropriate slow CVs and drives sampling of all thermally accessible metastable states starting from zero prior knowledge of the system. We make GREST accessible to the community via a publicly available open source Python package.

Selectivity and Ranking of Tight-Binding JAK-STAT Inhibitors Using Markovian Milestoning with Voronoi Tessellations

Janus kinases (JAK), a group of proteins in the nonreceptor tyrosine kinase (NRTKs) family, play a crucial role in growth, survival, and angiogenesis. They are activated by cytokines through the Janus kinase-signal transducer and activator of a transcription (JAK-STAT) signaling pathway. JAK-STAT signaling pathways have significant roles in the regulation of cell division, apoptosis, and immunity. Identification of the V617F mutation in the Janus homology 2 (JH2) domain of JAK2 leading to myeloproliferative disorders has stimulated great interest in the drug discovery community to develop JAK2-specific inhibitors. However, such inhibitors should be selective toward JAK2 over other JAKs and display an extended residence time. Recently, novel JAK2/STAT5 axis inhibitors (N-(1H-pyrazol-3-yl)pyrimidin-2-amino derivatives) have displayed extended residence times (hours or longer) on target and adequate selectivity excluding JAK3. To facilitate a deeper understanding of the kinase-inhibitor interactions and advance the development of such inhibitors, we utilize a multiscale Markovian milestoning with Voronoi tessellations (MMVT) approach within the Simulation-Enabled Estimation of Kinetic Rates v.2 (SEEKR2) program to rank order these inhibitors based on their kinetic properties and further explain the selectivity of JAK2 inhibitors over JAK3. Our approach investigates the kinetic and thermodynamic properties of JAK-inhibitor complexes in a user-friendly, fast, efficient, and accurate manner compared to other brute force and hybrid-enhanced sampling approaches.

Predicting Biomolecular Binding Kinetics: A Review

Biomolecular binding kinetics including the association (k_on) and dissociation (k_off) rates are critical parameters for therapeutic design of small-molecule drugs, peptides, and antibodies. Notably, the drug molecule residence time or dissociation rate has been shown to correlate with their efficacies better than binding affinities. A wide range of modeling approaches including quantitative structure-kinetic relationship models, Molecular Dynamics simulations, enhanced sampling, and Machine Learning has been developed to explore biomolecular binding and dissociation mechanisms and predict binding kinetic rates. Here, we review recent advances in computational modeling of biomolecular binding kinetics, with an outlook for future improvements.

Structure-Kinetics Relationships of Opioids from Metadynamics and Machine Learning Analysis

The nation's opioid overdose deaths reached an all-time high in 2021. The majority of deaths are due to synthetic opioids represented by fentanyl. Naloxone, which is a FDA-approved reversal agent, antagonizes opioids through competitive binding at the μ-opioid receptor (mOR). Thus, knowledge of the opioid's residence time is important for assessing the effectiveness of naloxone. Here, we estimated the residence times (τ) of 15 fentanyl and 4 morphine analogs using metadynamics and compared them with the most recent measurement of the opioid kinetic, dissociation, and naloxone inhibitory constants (Mann et al. Clin. Pharmacol. Therapeut. 2022, 120, 1020-1232). Importantly, the microscopic simulations offered a glimpse at the common binding mechanism and molecular determinants of dissociation kinetics for fentanyl analogs. The insights inspired us to develop a machine learning approach to analyze the kinetic impact of fentanyl's substituents based on the interactions with mOR residues. This proof-of-concept approach is general; for example, it may be used to tune ligand residence times in computer-aided drug discovery.

Ligand Gaussian Accelerated Molecular Dynamics 2 (LiGaMD2): Improved Calculations of Ligand Binding Thermodynamics and Kinetics with Closed Protein Pocket

Ligand binding thermodynamics and kinetics are critical parameters for drug design. However, it has proven challenging to efficiently predict ligand binding thermodynamics and kinetics from molecular simulations due to limited simulation timescales. Protein dynamics, especially in the ligand binding pocket, often plays an important role in ligand binding. Based on our previously developed Ligand Gaussian accelerated molecular dynamics (LiGaMD), here we present LiGaMD2 in which a selective boost potential was applied to both the ligand and protein residues in the binding pocket to improve sampling of ligand binding and dissociation. To validate the performance of LiGaMD2, the T4 lysozyme (T4L) mutants with open and closed pockets bound by different ligands were chosen as model systems. LiGaMD2 could efficiently capture repetitive ligand dissociation and binding within microsecond simulations of all T4L systems. The obtained ligand binding kinetic rates and free energies agreed well with available experimental values and previous modeling results. Therefore, LiGaMD2 provides an improved approach to sample opening of closed protein pockets for ligand dissociation and binding, thereby allowing for efficient calculations of ligand binding thermodynamics and kinetics.

Cooperation of structural motifs controls drug selectivity in cyclin-dependent kinases: an advanced theoretical analysis

Understanding drug selectivity mechanism is a long-standing issue for helping design drugs with high specificity. Designing drugs targeting cyclin-dependent kinases (CDKs) with high selectivity is challenging because of their highly conserved binding pockets. To reveal the underlying general selectivity mechanism, we carried out comprehensive analyses from both the thermodynamics and kinetics points of view on a representative CDK12 inhibitor. To fully capture the binding features of the drug-target recognition process, we proposed to use kinetic residue energy analysis (KREA) in conjunction with the community network analysis (CNA) to reveal the underlying cooperation effect between individual residues/protein motifs to the binding/dissociating process of the ligand. The general mechanism of drug selectivity in CDKs can be summarized as that the difference of structural cooperation between the ligand and the protein motifs leads to the difference of the energetic contribution of the key residues to the ligand. The proposed mechanisms may be prevalent in drug selectivity issues, and the insights may help design new strategies to overcome/attenuate the drug selectivity associated problems.

Complete reconstruction of dasatinib unbinding pathway from c-Src kinase by supervised molecular dynamics simulation method; assessing efficiency and trustworthiness of the method

Over the past years, rational drug design has gained lots of attention since employing it gave the world targeted therapy and more effective treatment solutions. Structure-based drug design (SBDD) is an excellent tool in rational drug design that takes advantage of accurate methods such as unbiased molecular dynamics (UMD) simulation for designing and optimizing molecular entities by understanding the binding and unbinding pathways of the binders. Supervised molecular dynamics (SuMD) simulation is a branch of UMD in which long-duration simulations are turned into short simulations, called replica, and a specific parameter is monitored throughout the simulation. In this work, we utilized this strategy to reconstruct the unbinding pathway of the anticancer drug dasatinib from its target protein, the c-Src kinase. Several unbinding events with valuable details were achieved. Then, to assess the efficiency and trustworthiness of the SuMD method, the unbinding pathway was also reconstructed by conventional UMD simulation, which uncovered some of the limitations of this method, such as limited sampling of the active site and finding the metastable states in the unbinding pathway. Furthermore, in times like these, when the world is desperate to find treatments for the Covid-19 disease, we think these methods are of exceptional value.Communicated by Ramaswamy H. Sarma.

Estimation of Drug-Target Residence Time by Targeted Molecular Dynamics Simulations

Drug-target residence time has emerged as a key selection factor in drug discovery since the binding duration of a drug molecule to its protein target can significantly impact its in vivo efficacy. The challenge in studying the residence time, in early drug discovery stages, lies in how to cost-effectively determine the residence time for the systematic assessment of compounds. Currently, there is still a lack of computational protocols to quickly estimate such a measure, particularly for large and flexible protein targets and drugs. Here, we report an efficient computational protocol, based on targeted molecular dynamics, to rank drug candidates by their residence time and to obtain insights into ligand-target dissociation mechanisms. The method was assessed on a dataset of 10 arylpyrazole inhibitors of CDK8, a large, flexible, and clinically important target, for which the experimental residence time of the inhibitors ranges from minutes to hours. The compounds were correctly ranked according to their estimated residence time scores compared to their experimental values. The analysis of protein-ligand interactions along the dissociation trajectories highlighted the favorable contribution of hydrophobic contacts to residence time and revealed key residues that strongly affect compound residence time.

Finite temperature string method with umbrella sampling using path collective variables: application to secondary structure change in a protein

The transition of an α-helix to a β-sheet in proteins is among the most complex conformational changes seen in biomolecular systems. Due to long time scales involved in the transition, it is challenging to study such protein conformational changes using direct molecular dynamics simulations. This limitation is typically overcome using an indirect approach wherein one computes the free energy landscape associated with the transition. Computation of free energy landscapes, however, requires a suitable set of collective variables that describe the transition. In this work, we demonstrate the use of path collective variables [D. Branduardi, F. L. Gervasio and M. Parrinello, J. Chem. Phys., 2007, 126, 054103] and combine it with the finite temperature string (FTS) method [E. Weinan, W. Ren and E. Vanden-Eijnden, J. Phys. Chem. B, 2005, 109, 6688-6693] to determine the molecular mechanisms involved during the structural transition of the mini G-protein from an α-helix to a β-hairpin. The transition from the α-helix proceeds via unfolding of the terminal residues, giving rise to a β-turn unfolded intermediate to eventually form the β-hairpin. Our proposed algorithm uses umbrella sampling simulations to simulate images along the string and the weighted histogram analysis to compute the free energy along the computed transition path. This work demonstrates that the string method in combination with path collective variables can be exploited to study complex protein conformational changes such as a complete change in the secondary structure.

Rare Event Kinetics from Adaptive Bias Enhanced Sampling

We introduce a novel enhanced sampling approach named on-the-fly probability enhanced sampling (OPES) flooding for calculating the kinetics of rare events from atomistic molecular dynamics simulation. This method is derived from the OPES approach [Invernizzi and Parrinello, J. Phys. Chem. Lett. 2020, 11, 7, 2731-2736], which has been recently developed for calculating converged free energy surfaces for complex systems. In this paper, we describe the theoretical details of the OPES flooding technique and demonstrate the application on three systems of increasing complexity: barrier crossing in a two-dimensional double-well potential, conformational transition in the alanine dipeptide in the gas phase, and the folding and unfolding of the chignolin polypeptide in an aqueous environment. From extensive tests, we show that the calculation of accurate kinetics not only requires the transition state to be bias-free, but the amount of bias deposited should also not exceed the effective barrier height measured along the chosen collective variables. In this vein, the possibility of computing rates from biasing suboptimal order parameters has also been explored. Furthermore, we describe the choice of optimum parameter combinations for obtaining accurate results from limited computational effort.

JND4135, a New Type II TRK Inhibitor, Overcomes TRK xDFG and Other Mutation Resistance In Vitro and In Vivo

The tropomyosin receptor kinases (TRKs) have been validated as effective targets in anticancer drug discovery. Two first-generation TRK inhibitors have been approved into market and displayed an encouraging therapeutic response in cancer patients harboring TRK fusion proteins. However, acquired resistance mediated by secondary TRK mutations especially in the xDFG motif remains an unsolved challenge in the clinic. Herein, we report the preclinical pharmacological results of JND4135, a new type II pan-TRK inhibitor, in overcoming TRK mutant resistance, including the xDFG mutations in vitro and in vivo. At a low nanomolar level, JND4135 displays a strong activity against wild-type TRKA/B/C and secondary mutations involving xDFG motif substitutions in kinase assays and cellular models; occupies the TRK proteins for an extended time; and has a slower dissociation rate than other TRK inhibitors. Moreover, by intraperitoneal injection, JND4135 exhibits tumor growth inhibition (TGI) of 81.0% at a dose of 40 mg/kg in BaF3-CD74-TRKA-G667C mice xenograft model. Therefore, JND4135 can be considered as a lead compound for drug discovery overcoming the resistance of TRK inhibitor drugs mediated by xDFG mutations.

Transition Rates and Efficiency of Collective Variables from Time-Dependent Biased Simulations

Simulations with adaptive time-dependent bias enable an efficient exploration of the conformational space of a system. However, the dynamic information is altered by the bias. Infrequent metadynamics recovers the transition rate of crossing a barrier, if the collective variables are ideal and there is no bias deposition near the transition state. Unfortunately, these conditions are not always fulfilled. To overcome these limitations, and inspired by single-molecule force spectroscopy, we use Kramers' theory for calculating the barrier-crossing rate when a time-dependent bias is added to the system. We assess the efficiency of collective variables parameter by measuring how efficiently the bias accelerates the transitions. We present approximate analytical expressions of the survival probability, reproducing the barrier-crossing time statistics and enabling the extraction of the unbiased transition rate even for challenging cases. We explore the limits of our method and provide convergence criteria to assess its validity.

Protein Flexibility and Dissociation Pathway Differentiation Can Explain Onset of Resistance Mutations in Kinases

Understanding how mutations render a drug ineffective is a problem of immense relevance. Often the mechanism through which mutations cause drug resistance can be explained purely through thermodynamics. However, the more perplexing situation is when two proteins have the same drug binding affinities but different residence times. In this work, we demonstrate how all-atom molecular dynamics simulations using recent developments grounded in statistical mechanics can provide a detailed mechanistic rationale for such variances. We discover dissociation mechanisms for the anti-cancer drug Imatinib (Gleevec) against wild-type and the N368S mutant of Abl kinase. We show how this point mutation triggers far-reaching changes in the protein's flexibility and leads to a different, much faster, drug dissociation pathway. We believe that this work marks an efficient and scalable approach to obtain mechanistic insight into resistance mutations in biomolecular receptors that are hard to explain using a structural perspective.

Enhanced-Sampling Simulations for the Estimation of Ligand Binding Kinetics: Current Status and Perspective

The dissociation rate (k off) associated with ligand unbinding events from proteins is a parameter of fundamental importance in drug design. Here we review recent major advancements in molecular simulation methodologies for the prediction of k off. Next, we discuss the impact of the potential energy function models on the accuracy of calculated k off values. Finally, we provide a perspective from high-performance computing and machine learning which might help improve such predictions.

Determination of Multidirectional Pathways for Ligand Release from the Receptor: A New Approach Based on Differential Evolution

Steered molecular dynamics (SMD) simulation is a powerful method in computer-aided drug design as it can be used to access the relative binding affinity with high precision but with low computational cost. The success of SMD depends on the choice of the direction along which the ligand is pulled from the receptor-binding site. In most simulations, the unidirectional pathway was used, but in some cases, this choice resulted in the ligand colliding with the complex surface of the exit tunnel. To overcome this difficulty, several variants of SMD with multidirectional pulling have been proposed, but they are not completely devoid of disadvantages. Here, we have proposed to determine the direction of pulling with a simple scoring function that minimizes the receptor–ligand interaction, and an optimization algorithm called differential evolution is used for energy minimization. The effectiveness of our protocol was demonstrated by finding expulsion pathways of Huperzine A and camphor from the binding site of Torpedo California acetylcholinesterase and P450cam proteins, respectively, and comparing them with the previous results obtained using memetic sampling and random acceleration molecular dynamics. In addition, by applying this protocol to a set of ligands bound with LSD1 (lysine specific demethylase 1), we obtained a much higher correlation between the work of pulling force and experimental data on the inhibition constant IC50 compared to that obtained using the unidirectional approach based on minimal steric hindrance.

Extracting Atomic Contributions to Binding Free Energy Using Molecular Dynamics Simulations with Mixed Solvents (MDmix)

BACKGROUND

Mixed solvents MD (MDmix) simulations have proved to be a useful and increasingly accepted technique with several applications in structure-based drug discovery. One of the assumptions behind the methodology is the transferability of free energy values from the simulated cosolvent molecules to larger drug-like molecules. However, the binding free energy maps (ΔGbind) calculated for the different moieties of the cosolvent molecules (e.g. a hydroxyl map for the ethanol) are largely influenced by the rest of the solvent molecule and do not reflect the intrinsic affinity of the moiety in question. As such, they are hardly transferable to different molecules.

METHOD

To achieve transferable energies, we present here a method for decomposing the molecular binding free energy into accurate atomic contributions.

RESULT

We demonstrate with two qualitative visual examples how the corrected energy maps better match known binding hotspots and how they can reveal hidden hotspots with actual drug design potential.

CONCLUSION

Atomic decomposition of binding free energies derived from MDmix simulations provides transferable and quantitative binding free energy maps.

From Data to Knowledge: Systematic Review of Tools for Automatic Analysis of Molecular Dynamics Output

An increasing number of crystal structures available on one side, and the boost of computational power available for computer-aided drug design tasks on the other, have caused that the structure-based drug design tools are intensively used in the drug development pipelines. Docking and molecular dynamics simulations, key representatives of the structure-based approaches, provide detailed information about the potential interaction of a ligand with a target receptor. However, at the same time, they require a three-dimensional structure of a protein and a relatively high amount of computational resources. Nowadays, as both docking and molecular dynamics are much more extensively used, the amount of data output from these procedures is also growing. Therefore, there are also more and more approaches that facilitate the analysis and interpretation of the results of structure-based tools. In this review, we will comprehensively summarize approaches for handling molecular dynamics simulations output. It will cover both statistical and machine-learning-based tools, as well as various forms of depiction of molecular dynamics output.

Combined Free-Energy Calculation and Machine Learning Methods for Understanding Ligand Unbinding Kinetics

The determination of drug residence times, which define the time an inhibitor is in complex with its target, is a fundamental part of the drug discovery process. Synthesis and experimental measurements of kinetic rate constants are, however, expensive and time consuming. In this work, we aimed to obtain drug residence times computationally. Furthermore, we propose a novel algorithm to identify molecular design objectives based on ligand unbinding kinetics. We designed an enhanced sampling technique to accurately predict the free-energy profiles of the ligand unbinding process, focusing on the free-energy barrier for unbinding. Our method first identifies unbinding paths determining a corresponding set of internal coordinates (ICs) that form contacts between the protein and the ligand; it then iteratively updates these interactions during a series of biased molecular dynamics (MD) simulations to reveal the ICs that are important for the whole of the unbinding process. Subsequently, we performed finite-temperature string simulations to obtain the free-energy barrier for unbinding using the set of ICs as a complex reaction coordinate. Importantly, we also aimed to enable the further design of drugs focusing on improved residence times. To this end, we developed a supervised machine learning (ML) approach with inputs from unbiased “downhill” trajectories initiated near the transition state (TS) ensemble of the string unbinding path. We demonstrate that our ML method can identify key ligand–protein interactions driving the system through the TS. Some of the most important drugs for cancer treatment are kinase inhibitors. One of these kinase targets is cyclin-dependent kinase 2 (CDK2), an appealing target for anticancer drug development. Here, we tested our method using two different CDK2 inhibitors for the potential further development of these compounds. We compared the free-energy barriers obtained from our calculations with those observed in available experimental data. We highlighted important interactions at the distal ends of the ligands that can be targeted for improved residence times. Our method provides a new tool to determine unbinding rates and to identify key structural features of the inhibitors that can be used as starting points for novel design strategies in drug discovery.

Decisive role of water and protein dynamics in residence time of p38α MAP kinase inhibitors

Target residence time plays a crucial role in the pharmacological activity of small molecule inhibitors. Little is known, however, about the underlying causes of inhibitor residence time at the molecular level, which complicates drug optimization processes. Here, we employ all-atom molecular dynamics simulations (~400 μs in total) to gain insight into the binding modes of two structurally similar p38α MAPK inhibitors (type I and type I½) with short and long residence times that otherwise show nearly identical inhibitory activities in the low nanomolar IC₅₀ range. Our results highlight the importance of protein conformational stability and solvent exposure, buried surface area of the ligand and binding site resolvation energy for residence time. These findings are further confirmed by simulations with a structurally diverse short residence time inhibitor SB203580. In summary, our data provide guidance in compound design when aiming for inhibitors with improved target residence time.

The molecular determinants of the residence time of a small molecule inhibitor at its target protein are not well understood. Here, Pantsar et al. show that the target protein’s conformational stability and solvent exposure are key factors governing the target residence time of kinase inhibitors.

Gaussian Accelerated Molecular Dynamics Simulations Investigation on the Mechanism of Angiotensin-Converting Enzyme (ACE) C-Domain Inhibition by Dipeptides

Angiotensin-converting enzyme (ACE)-inhibitory peptides extracted from food proteins can lower blood pressure by inhibiting ACE activity. A recent study showed that the inhibitory activity of IY (Ile-Tyr, a dipeptide derived from soybean protein) against ACE was much higher than that of LL (Leu-Leu), although they had similar hydrophobic and predicted activity values. It was difficult to reveal the deep molecular mechanism underlying this phenomenon by traditional experimental methods. The Apo and two complex systems (i.e., ACE–LL and ACE–IY) were therefore subjected to 1 μs long Gaussian accelerated molecular dynamics (GaMD) simulations. The results showed that the binding of IY can cause obvious contraction of the active site of ACE, mainly manifested by a significant lateral shift of α13, α14, and α15. In addition, hinge 2 and hinge 3 were more stable in the ACE–IY system, while these phenomena were not present in the ACE–LL system. Moreover, the α10 of the IY-bound ACE kept an inward state during the simulation progress, which facilitated the ACE to remain closed. However, for the LL-bound ACE, the α10 switched between two outward states. To sum up, our study provides detailed insights into inhibitor-induced conformational changes in ACE that may help in the design of specific inhibitors targeting ACE for the treatment of hypertension.

Binding Analysis Using Accelerated Molecular Dynamics Simulations and Future Perspectives

Biomolecular recognition such as binding of small molecules, nucleic acids, peptides and proteins to their target receptors plays key roles in cellular function and has been targeted for therapeutic drug design. Molecular dynamics (MD) is a computational approach to analyze these binding processes at an atomistic level, which provides valuable understandings of the mechanisms of biomolecular recognition. However, the rather slow biomolecular binding events often present challenges for conventional MD (cMD), due to limited simulation timescales (typically over hundreds of nanoseconds to tens of microseconds). In this regard, enhanced sampling methods, particularly accelerated MD (aMD), have proven useful to bridge the gap and enable all-atom simulations of biomolecular binding events. Here, we will review the recent method developments of Gaussian aMD (GaMD), ligand GaMD (LiGaMD) and peptide GaMD (Pep-GaMD), which have greatly expanded our capabilities to simulate biomolecular binding processes. Spontaneous binding of various biomolecules to their receptors has been successfully simulated by GaMD. Microsecond LiGaMD and Pep-GaMD simulations have captured repetitive binding and dissociation of small-molecule ligands and highly flexible peptides, and thus enabled ligand/peptide binding thermodynamics and kinetics calculations. We will also present relevant application studies in simulations of important drug targets and future perspectives for rational computer-aided drug design.

Deciphering the mechanism of action of 5FDQD and the design of new neutral analogues for the FMN riboswitch: a well-tempered metadynamics simulation study

The FMN riboswitch is a novel drug target for the design of new antibiotics, and efforts have been made to design new charged and uncharged ligands. Uncharged ligands have shown advantages of not requiring any transporter for intracellular transport or proteins for their phosphorylation. 5FDQD (5-(3-(4-fluorophenyl)butyl)-7,8-dimethylpyrido(3,4-b)quinoxaline-1,3(2H,5H)-dione) is a recently reported neutral ligand for the FMN riboswitch active against Clostridium difficile infection in mice. However, the crystal structure of the 5FDQD bound FMN riboswitch is not available, and the mechanism of ligand binding and triggering the function of the riboswitch is not well understood. We have examined 5FDQD for its binding affinity with the FMN riboswitch using the well-tempered metadynamics (WT-MtD) simulation technique. The crystal structure of the FMN riboswitch shows that the FMN interacts with the J4/5 region through the phosphate group with G62; however, the uncharged ligands take advantage of π-π stacking interactions with the same residue of the riboswitch observed from the WT-MtD simulation results. The simulation results show that the presence of fluorine on the phenyl ring in 5FDQD is important to enhance the binding affinity of the neutral ligands with the FMN riboswitch. The WT-MtD results showed that the 1,2-difluoro substitution on the phenyl ring in 5FDQD (FMN-difluoro2) and the 1,3 positions in the phenyl ring (FMN-difluoro1) showed weaker binding energy with the FMN riboswitch compared to 5FDQD. The substitution of another fluorine atom at the 5-position of the phenyl ring (FMN-trifluoro) showed a comparable binding affinity (∼-31.4 kcal mol^-1) to 5FDQD. Electron-donating substitution on the phenyl ring such as the amino group also lowered the binding affinity (-28.8 kcal mol^-1) with the riboswitch compared to 5FDQD. The computed results suggest that the position and nature of substitution in the phenyl ring of the uncharged ligands affect the overall binding and such a delicate balance is important to achieve superior binding affinity with the FMN riboswitch.

Challenges and frontiers of computational modelling of biomolecular recognition

Biomolecular recognition including binding of small molecules, peptides and proteins to their target receptors plays a key role in cellular function and has been targeted for therapeutic drug design. However, the high flexibility of biomolecules and slow binding and dissociation processes have presented challenges for computational modelling. Here, we review the challenges and computational approaches developed to characterise biomolecular binding, including molecular docking, molecular dynamics simulations (especially enhanced sampling) and machine learning. Further improvements are still needed in order to accurately and efficiently characterise binding structures, mechanisms, thermodynamics and kinetics of biomolecules in the future.

Creating Maps of the Ligand Binding Landscape for Kinetics-Based Drug Discovery

Simulations of ligand-protein interactions can be very useful for drug design and to gain biological insight. Full pathways of ligand-protein binding can be used to get information about ligand binding transition states, which form the rate-limiting step of the binding and release processes. However, these simulations are typically limited by the presence of large energy barriers that separate stable poses of interest. Here we describe a simulation protocol for exploring and analyzing landscapes of ligand-protein interactions that makes use of molecular docking, enhanced molecular simulation with the weighted ensemble algorithm, and network analysis. It can be accomplished using a modest cluster of graphics processing units and freely accessible software. This protocol focuses on the construction and analysis of a network model of ligand binding poses and provides links to resources that describe the other steps in more detail. The end result of this protocol is a map of the ligand-protein binding landscape that identifies transition states of the ligand binding pathway, as well as alternative bound poses that could be stabilized with modifications to the ligand.

Free Energy Surfaces and Barriers for Vacancy Diffusion on Al(100), Al(110), Al(111) Reconstructed Surfaces

Metadynamics is a popular enhanced sampling method based on the recurrent application of a history-dependent adaptive bias potential that is a function of a selected number of appropriately chosen collective variables. In this work, using metadynamics simulations, we performed a computational study for the diffusion of vacancies on three different Al surfaces [reconstructed Al(100), Al(110), and Al(111) surfaces]. We explored the free energy landscape of diffusion and estimated the barriers associated with this process on each surface. It is found that the surfaces are unique regarding vacancy diffusion. More specically, the reconstructed Al(110) surface presents four metastable states on the free energy surface having sizable and connected passage-ways with an energy barrier of height 0.55 eV. On the other hand, the reconstructed Al(100)/Al(111) surfaces exhibit two/three metastable states, respectively, with an energy barrier of height 0.33 eV. The findings in this study can help to understand surface vacancy diffusion in technologically relevant Al surfaces.

Kinetics and Mechanism of Fentanyl Dissociation from the μ-Opioid Receptor

Driven by illicit fentanyl, opioid related deaths have reached the highest level in 2020. Currently, an opioid overdose is resuscitated by the use of naloxone, which competitively binds and antagonizes the μ-opioid receptor (mOR). Thus, knowledge of the residence times of opioids at mOR and the unbinding mechanisms is valuable for assessing the effectiveness of naloxone. In the present study, we calculate the fentanyl-mOR dissociation time and elucidate the mechanism by applying an enhanced sampling molecular dynamics (MD) technique. Two sets of metadynamics simulations with different initial structures were performed while accounting for the protonation state of the conserved H297^6.52, which has been suggested to modulate the ligand-mOR affinity and binding mode. Surprisingly, with the Nδ-protonated H297^6.52, fentanyl can descend as much as 10 Å below the level of the conserved D147^3.32 before escaping the receptor and has a calculated residence time τ of 38 s. In contrast, with the Nϵ- and doubly protonated H297^6.52, the calculated τ are 2.6 and 0.9 s, respectively. Analysis suggests that formation of the piperidine-Hid297 hydrogen bond strengthens the hydrophobic contacts with the transmembrane helix (TM) 6, allowing fentanyl to explore a deep pocket. Considering the experimental τ of ∼4 min for fentanyl and the role of TM6 in mOR activation, the deep insertion mechanism may be biologically relevant. The work paves the way for large-scale computational predictions of opioid dissociation rates to inform evaluation of strategies for opioid overdose reversal. The profound role of the histidine protonation state found here may shift the paradigm in computational studies of ligand-receptor kinetics.

Enhanced Sampling Approach to the Induced-Fit Docking Problem in Protein-Ligand Binding: The Case of Mono-ADP-Ribosylation Hydrolase Inhibitors

Enhanced sampling methods can predict free-energy landscapes associated with protein/ligand binding, characterizing the involved intermolecular interactions in a precise way. However, these in silico approaches can be challenged by induced-fit effects. Here, we present a variant of volume-based metadynamics tailored to tackle this problem in a general and efficient way. The validity of the approach is established by applying it to substrate/enzyme complexes of pharmacological relevance: mono-ADP-ribose (ADPr) in complex with mono-ADP-ribosylation hydrolases (MacroD1 and MacroD2), where induced-fit phenomena are known to be significant. The calculated binding free energies are consistent with experiments, with an absolute error smaller than 0.5 kcal/mol. Our simulations reveal that in all circumstances, the active loops, delimiting the boundaries of the binding site, undergo significant conformation rearrangements upon ligand binding. The calculations further provide, for the first time, the molecular basis of ADPr specificity and the relative changes in its experimental binding affinity on passing from MacroD1 to MacroD2 and all its mutants. Our study paves the way to the quantitative description of induced-fit events in molecular recognition.

Structural Basis for Developing Multitarget Compounds Acting on Cysteinyl Leukotriene Receptor 1 and G-Protein-Coupled Bile Acid Receptor 1

G-protein-coupled receptors (GPCRs) are the molecular target of 40% of marketed drugs and the most investigated structures to develop novel therapeutics. Different members of the GPCRs superfamily can modulate the same cellular process acting on diverse pathways, thus representing an attractive opportunity to achieve multitarget drugs with synergic pharmacological effects. Here, we present a series of compounds with dual activity toward cysteinyl leukotriene receptor 1 (CysLT₁R) and G-protein-coupled bile acid receptor 1 (GPBAR1). They are derivatives of REV5901─the first reported dual compound─with therapeutic potential in the treatment of colitis and other inflammatory processes. We report the binding mode of the most active compounds in the two GPCRs, revealing unprecedented structural basis for future drug design studies, including the presence of a polar group opportunely spaced from an aromatic ring in the ligand to interact with Arg79^2.60 of CysLT₁R and achieve dual activity.