Accuracy of Molecular Simulation-Based Predictions of koff Values: A Metadynamics Study

MiMiC: A high-performance framework for multiscale molecular dynamics simulations

MiMiC is a framework for performing multiscale simulations in which loosely coupled external programs describe individual subsystems at different resolutions and levels of theory. To make it highly efficient and flexible, we adopt an interoperable approach based on a multiple-program multiple-data (MPMD) paradigm, serving as an intermediary responsible for fast data exchange and interactions between the subsystems. The main goal of MiMiC is to avoid interfering with the underlying parallelization of the external programs, including the operability on hybrid architectures (e.g., CPU/GPU), and keep their setup and execution as close as possible to the original. At the moment, MiMiC offers an efficient implementation of electrostatic embedding quantum mechanics/molecular mechanics (QM/MM) that has demonstrated unprecedented parallel scaling in simulations of large biomolecules using CPMD and GROMACS as QM and MM engines, respectively. However, as it is designed for high flexibility with general multiscale models in mind, it can be straightforwardly extended beyond QM/MM. In this article, we illustrate the software design and the features of the framework, which make it a compelling choice for multiscale simulations in the upcoming era of exascale high-performance computing.

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.

Multiscale biomolecular simulations in the exascale era

The complexity of biological systems and processes, spanning molecular to macroscopic scales, necessitates the use of multiscale simulations to get a comprehensive understanding. Quantum mechanics/molecular mechanics (QM/MM) molecular dynamics (MD) simulations are crucial for capturing processes beyond the reach of classical MD simulations. The advent of exascale computing offers unprecedented opportunities for scientific exploration, not least within life sciences, where simulations are essential to unravel intricate molecular mechanisms underlying biological processes. However, leveraging the immense computational power of exascale computing requires innovative algorithms and software designs. In this context, we discuss the current status and future prospects of multiscale biomolecular simulations on exascale supercomputers with a focus on QM/MM MD. We highlight our own efforts in developing a versatile and high-performance multiscale simulation framework with the aim of efficient utilization of state-of-the-art supercomputers. We showcase its application in uncovering complex biological mechanisms and its potential for leveraging exascale computing.

Calcium-Driven In Silico Inactivation of a Human Olfactory Receptor

Conformational changes as well as molecular determinants related to the activation and inactivation of olfactory receptors are still poorly understood due to the intrinsic difficulties in the structural determination of this GPCR family. Here, we perform, for the first time, the in silico inactivation of human olfactory receptor OR51E2, highlighting the possible role of calcium in this receptor state transition. Using molecular dynamics simulations, we show that a divalent ion in the ion binding site, coordinated by two acidic residues at positions 2.50 and 3.39 conserved across most ORs, stabilizes the receptor in its inactive state. In contrast, protonation of the same two acidic residues is not sufficient to drive inactivation within the microsecond timescale of our simulations. Our findings suggest a novel molecular mechanism for OR inactivation, potentially guiding experimental validation and offering insights into the possible broader role of divalent ions in GPCR signaling.

Kinetic and structural details of urease inactivation by thiuram disulphides

This paper reports on the molecular details of the reactivity of urease, a nickel-dependent enzyme that catalyses the last step of organic nitrogen mineralization, with thiuram disulphides, a class of molecules known to inactivate the enzyme with high efficacy but for which the mechanism of action had not been yet established. IC50 values of tetramethylthiuram disulphide (TMTD or Thiram) and tetraethylthiuram disulphide (TETD or Disulfiram) in the low micromolar range were determined for plant and bacterial ureases. The X-ray crystal structure of Sporosarcina pasteurii urease inactivated by Thiram, determined at 1.68 Å resolution, revealed the presence of a covalent modification of the catalytically essential cysteine residue. This is located on the flexible flap that modulates the size of the active site channel and cavity. Formation of a Cys-S-S-C(S)-N(CH3)2 functionality responsible for enzyme inactivation was observed. Quantum-mechanical calculations carried out to rationalise the large reactivity of the active site cysteine support the view that a conserved histidine residue, adjacent to the cysteine in the active site flap, modulates the charge and electron density along the thiol SH bond by shifting electrons towards the sulphur atom and rendering the thiol proton more reactive. We speculate that this proton could be transferred to the nickel-coordinated urea amide group to yield a molecule of ammonia from the generated Curea-NH3+ functionality during catalysis.

QMrebind: incorporating quantum mechanical force field reparameterization at the ligand binding site for improved drug-target kinetics through milestoning simulations

Understanding the interaction of ligands with biomolecules is an integral component of drug discovery and development. Challenges for computing thermodynamic and kinetic quantities for pharmaceutically relevant receptor-ligand complexes include the size and flexibility of the ligands, large-scale conformational rearrangements of the receptor, accurate force field parameters, simulation efficiency, and sufficient sampling associated with rare events. Our recently developed multiscale milestoning simulation approach, SEEKR2 (Simulation Enabled Estimation of Kinetic Rates v.2), has demonstrated success in predicting unbinding (koff) kinetics by employing molecular dynamics (MD) simulations in regions closer to the binding site. The MD region is further subdivided into smaller Voronoi tessellations to improve the simulation efficiency and parallelization. To date, all MD simulations are run using general molecular mechanics (MM) force fields. The accuracy of calculations can be further improved by incorporating quantum mechanical (QM) methods into generating system-specific force fields through reparameterizing ligand partial charges in the bound state. The force field reparameterization process modifies the potential energy landscape of the bimolecular complex, enabling a more accurate representation of the intermolecular interactions and polarization effects at the bound state. We present QMrebind (Quantum Mechanical force field reparameterization at the receptor-ligand binding site), an ORCA-based software that facilitates reparameterizing the potential energy function within the phase space representing the bound state in a receptor-ligand complex. With SEEKR2 koff estimates and experimentally determined kinetic rates, we compare and interpret the receptor-ligand unbinding kinetics obtained using the newly reparameterized force fields for model host-guest systems and HSP90-inhibitor complexes. This method provides an opportunity to achieve higher accuracy in predicting receptor-ligand koff rate constants.

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.

Investigating the Unbinding of Muscarinic Antagonists from the Muscarinic 3 Receptor

Patient symptom relief is often heavily influenced by the residence time of the inhibitor-target complex. For the human muscarinic receptor 3 (hMR3), tiotropium is a long-acting bronchodilator used in conditions such as asthma or chronic obstructive pulmonary disease (COPD). The mechanistic insights into this inhibitor remain unclear; specifically, the elucidation of the main factors determining the unbinding rates could help develop the next generation of antimuscarinic agents. Using our novel unbinding algorithm, we were able to investigate ligand dissociation from hMR3. The unbinding paths of tiotropium and two of its analogues, N-methylscopolamin and homatropine methylbromide, show a consistent qualitative mechanism and allow us to identify the structural bottleneck of the process. Furthermore, our machine learning-based analysis identified key roles of the ECL2/TM5 junction involved in the transition state. Additionally, our results point to relevant changes at the intracellular end of the TM6 helix leading to the ICL3 kinase domain, highlighting the closest residue L482. This residue is located right between two main protein binding sites involved in signal transduction for hMR3's activation and regulation. We also highlight key pharmacophores of tiotropium that play determining roles in the unbinding kinetics and could aid toward drug design and lead optimization.

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.

Drug Design in the Exascale Era: A Perspective from Massively Parallel QM/MM Simulations

The initial phases of drug discovery - in silico drug design - could benefit from first principle Quantum Mechanics/Molecular Mechanics (QM/MM) molecular dynamics (MD) simulations in explicit solvent, yet many applications are currently limited by the short time scales that this approach can cover. Developing scalable first principle QM/MM MD interfaces fully exploiting current exascale machines - so far an unmet and crucial goal - will help overcome this problem, opening the way to the study of the thermodynamics and kinetics of ligand binding to protein with first principle accuracy. Here, taking two relevant case studies involving the interactions of ligands with rather large enzymes, we showcase the use of our recently developed massively scalable Multiscale Modeling in Computational Chemistry (MiMiC) QM/MM framework (currently using DFT to describe the QM region) to investigate reactions and ligand binding in enzymes of pharmacological relevance. We also demonstrate for the first time strong scaling of MiMiC-QM/MM MD simulations with parallel efficiency of ∼70% up to >80,000 cores. Thus, among many others, the MiMiC interface represents a promising candidate toward exascale applications by combining machine learning with statistical mechanics based algorithms tailored for exascale supercomputers.

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.

Predicting Protein-Ligand Binding and Unbinding Kinetics with Biased MD Simulations and Coarse-Graining of Dynamics: Current State and Challenges

The prediction of drug-target binding and unbinding kinetics that occur on time scales between milliseconds and several hours is a prime challenge for biased molecular dynamics simulation approaches. This Perspective gives a concise summary of the theory and the current state-of-the-art of such predictions via biased simulations, of insights into the molecular mechanisms defining binding and unbinding kinetics as well as of the extraordinary challenges predictions of ligand kinetics pose in comparison to binding free energy predictions.

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.

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.

Mechanism of Calcium Permeation in a Glutamate Receptor Ion Channel

The α-amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid receptors (AMPARs) are neurotransmitter-activated cation channels ubiquitously expressed in vertebrate brains. The regulation of calcium flux through the channel pore by RNA-editing is linked to synaptic plasticity while excessive calcium influx poses a risk for neurodegeneration. Unfortunately, the molecular mechanisms underlying this key process are mostly unknown. Here, we investigated calcium conduction in calcium-permeable AMPAR using Molecular Dynamics (MD) simulations with recently introduced multisite force-field parameters for Ca²⁺. Our calculations are consistent with experiment and explain the distinct calcium permeability in different RNA-edited forms of GluA2. For one of the identified metal binding sites, multiscale Quantum Mechanics/Molecular Mechanics (QM/MM) simulations further validated the results from MD and revealed small but reproducible charge transfer between the metal ion and its first solvation shell. In addition, the ion occupancy derived from MD simulations independently reproduced the Ca²⁺ binding profile in an X-ray structure of an NaK channel mimicking the AMPAR selectivity filter. This integrated study comprising X-ray crystallography, multisite MD, and multiscale QM/MM simulations provides unprecedented insights into Ca²⁺ permeation mechanisms in AMPARs, and paves the way for studying other biological processes in which Ca²⁺ plays a pivotal role.

Molecular dynamics simulations of the calmodulin-induced α-helix in the SK2 calcium-gated potassium ion channel

The family of small-conductance Ca²⁺-activated potassium ion channels (SK channels) is composed of four members (SK1, SK2, SK3, and SK4) involved in neuron-firing regulation. The gating of these channels depends on the intracellular Ca²⁺ concentration, and their sensitivity to this ion is provided by calmodulin (CaM). This protein binds to a specific region in SK channels known as the calmodulin-binding domain (CaMBD), an event which is essential for their gating. While CaMBDs are typically disordered in the absence of CaM, the SK2 channel subtype displays a small prefolded α-helical region in its CaMBD even if CaM is not present. This small helix is known to turn into a full α-helix upon CaM binding, although the molecular-level details for this conversion are not fully understood yet. In this work, we offer new insights on this physiologically relevant process by means of enhanced sampling, atomistic Hamiltonian replica exchange molecular dynamics simulations, providing a more detailed understanding of CaM binding to this target. Our results show that CaM is necessary for inducing a full α-helix along the SK2 CaMBD through hydrophobic interactions with V426 and L427. However, it is also necessary that W431 does not compete for these interactions; the role of the small prefolded α-helix in the SK2 CaMBD would be to stabilize W431 so that this is the case. In conclusion, our findings provide further insight into a key interaction between CaM and SK channels that is important for channel sensitivity to Ca²⁺.

Hierarchical Materials from High Information Content Macromolecular Building Blocks: Construction, Dynamic Interventions, and Prediction

Hierarchical materials that exhibit order over multiple length scales are ubiquitous in nature. Because hierarchy gives rise to unique properties and functions, many have sought inspiration from nature when designing and fabricating hierarchical matter. More and more, however, nature's own high-information content building blocks, proteins, peptides, and peptidomimetics, are being coopted to build hierarchy because the information that determines structure, function, and interfacial interactions can be readily encoded in these versatile macromolecules. Here, we take stock of recent progress in the rational design and characterization of hierarchical materials produced from high-information content blocks with a focus on stimuli-responsive and "smart" architectures. We also review advances in the use of computational simulations and data-driven predictions to shed light on how the side chain chemistry and conformational flexibility of macromolecular blocks drive the emergence of order and the acquisition of hierarchy and also on how ionic, solvent, and surface effects influence the outcomes of assembly. Continued progress in the above areas will ultimately usher in an era where an understanding of designed interactions, surface effects, and solution conditions can be harnessed to achieve predictive materials synthesis across scale and drive emergent phenomena in the self-assembly and reconfiguration of high-information content building blocks.

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.

Advances in computational methods for ligand binding kinetics

Binding kinetic parameters can be correlated with drug efficacy, which in recent years led to the development of various computational methods for predicting binding kinetic rates and gaining insight into protein-drug binding paths and mechanisms. In this review, we introduce and compare computational methods recently developed and applied to two systems, trypsin-benzamidine and kinase-inhibitor complexes. Methods involving enhanced sampling in molecular dynamics simulations or machine learning can be used not only to predict kinetic rates, but also to reveal factors modulating the duration of residence times, selectivity, and drug resistance to mutations. Methods which require less computational time to make predictions are highlighted, and suggestions to reduce the error of computed kinetic rates are presented.

The role of ATP in solubilizing RNA-binding protein fused in sarcoma

Intrinsically disordered protein (IDP) plays an important role in liquid-liquid phase separation (LLPS). RNA-binding protein fused in sarcoma (FUS) is a well-studied IDP that induces LLPS since its low-complexity core region (FUS-LC-core) is essential for droplet formation through contacts between FUS-LC-cores. Several experimental studies have reported that adenosine triphosphate (ATP) concentrations modulate LLPS-driven droplet formation through the dissolution of FUS. To elucidate the role of ATP in this dissolution, microsecond-order all-atom molecular dynamics (MD) simulations were performed for a crowded system of FUS-LC-cores in the presence of multiple ATP molecules. Our analysis revealed that the adenine group of ATP frequently contacted the FUS-LC-core, and the phosphoric acid group of ATP was exposed to the external solvent, which promoted both hydration and solubilization of FUS.

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.

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.

Nonlinear Reaction Coordinate of an Enzyme Catalyzed Proton Transfer Reaction

We present an in-depth study on the theoretical calculation of an optimum reaction coordinate as a linear or nonlinear combination of important collective variables (CVs) sampled from an ensemble of reactive transition paths for an intramolecular proton transfer reaction catalyzed by the enzyme human carbonic anhydrase (HCA) II. The linear models are optimized by likelihood maximization for a given number of CVs. The nonlinear models are based on an artificial neural network with the same number of CVs and optimized by minimizing the root-mean-square error in comparison to a training set of committor estimators generated for the given transition. The nonlinear reaction coordinate thus obtained yields the free energy of activation and rate constant as 9.46 kcal mol^-1 and 1.25 × 10⁶ s^-1, respectively. These estimates are found to be in quantitative agreement with the known experimental results. We have also used an extended autoencoder model to show that a similar analysis can be carried out using a single CV only. The resultant free energies and kinetics of the reaction slightly overestimate the experimental data. The implications of these results are discussed using a detailed microkinetic scheme of the proton transfer reaction catalyzed by HCA II.

Long-Timescale Simulations Revealed Critical Non-Conserved Residues of Phosphodiesterases Affecting Selectivity of BAY60-7550

A major obstacle of the selective inhibitor design for specific human phosphodiesterase (PDE) is that highly conserved catalytic pockets are difficult to be distinguished by inhibitor molecules. To overcome this, a feasible path is to understand the molecular determinants underlying the selectivity of current inhibitors. BAY60-7550 (BAY for short; IC₅₀ = 4.7 nM) is a highly selective inhibitor targeting PDE2A which is a dual-specificity PDE and an attractive target for therapeutic intervention of the central nervous system (CNS) disorders. Recent studies suggest that molecular determinants may be in binding processes of BAY. However, a detailed understanding of these processes are still lacking. To explore these processes, High-Throughput Molecular Dynamics (HTMD) simulations were performed to reproduce the spontaneous association of BAY with catalytic pockets of 4 PDE isoforms; Ligand Gaussian Accelerated Molecular Dynamics (LiGaMD) simulations were performed to reproduce the unbinding-rebinding processes of FKG and 10.13039/100016266MC2, two pyrazolopyrimidinone PDE2A selective inhibitors, in the PDE2A system. The produced molecular trajectories were analyzed by the Markov state model (MSM) and the molecular mechanics/generalized Born surface area (MM/GBSA). The results showed that the non-covalent interactions between the non-conserved residues and BAY, especially the hydrogen bonds, determined the unique binding pathways of BAY on the surface of PDE2A. These pathways were different from those of BAY on the surface of the other three PDE isoforms and the binding pathways of the other two PDE2A inhibitors in PDE2A systems. These differences were ultimately reflected in the high selectivity of this inhibitor for PDE2A. As a result, this study demonstrates the critical role of the binding processes in the selectivity of BAY, and also identifies the key non-conserved residues affecting the binding processes of BAY. Thus, this study provides a new perspective and data support for the further development of BAY-derived inhibitors targeting PDE2A.

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.

Predicting Ion Channel Conductance via Dissipation-Corrected Targeted Molecular Dynamics and Langevin Equation Simulations

Ion channels are important proteins for physiological information transfer and functional control. To predict the microscopic origins of their voltage-conductance characteristics, here we applied dissipation-corrected targeted molecular dynamics in combination with Langevin equation simulations to potassium diffusion through the gramicidin A channel as a test system. Performing a nonequilibrium principal component analysis on backbone dihedral angles, we find coupled protein-ion dynamics to occur during ion transfer. The dissipation-corrected free energy profiles correspond well to predictions from other biased simulation methods. The incorporation of an external electric field in Langevin simulations enables the prediction of macroscopic observables in the form of I-V characteristics.

Protein-Ligand Dissociation Rate Constant from All-Atom Simulation

Dissociation of a ligand isoniazid from a protein catalase was investigated using all-atom molecular dynamics (MD) simulations. Random acceleration MD (τ-RAMD) was used, in which a random artificial force applied to the ligand facilitates its dissociation. We have suggested a novel approach to extrapolate such obtained dissociation times to the zero-force limit assuming never before attempted universal exponential dependence of the bond strength on the applied force, allowing direct comparison with experimentally measured values. We have found that our calculated dissociation time was equal to 36.1 s with statistically significant values distributed in the interval of 0.2-72.0 s, which quantitatively matches the experimental value of 50 ± 8 s despite the extrapolation over 9 orders of magnitude in time.

G Protein-Coupled Receptor-Ligand Dissociation Rates and Mechanisms from τRAMD Simulations

There is a growing appreciation of the importance of drug-target binding kinetics for lead optimization. For G protein-coupled receptors (GPCRs), which mediate signaling over a wide range of time scales, the drug dissociation rate is often a better predictor of in vivo efficacy than binding affinity, although it is more challenging to compute. Here, we assess the ability of the τ-Random Acceleration Molecular Dynamics (τRAMD) approach to reproduce relative residence times and reveal dissociation mechanisms and the effects of allosteric modulation for two important membrane-embedded drug targets: the β2-adrenergic receptor and the muscarinic acetylcholine receptor M2. The dissociation mechanisms observed in the relatively short RAMD simulations (in which molecular dynamics (MD) simulations are performed using an additional force with an adaptively assigned random orientation applied to the ligand) are in general agreement with much more computationally intensive conventional MD and metadynamics simulations. Remarkably, although decreasing the magnitude of the random force generally reduces the number of egress routes observed, the ranking of ligands by dissociation rate is hardly affected and agrees well with experiment. The simulations also reproduce changes in residence time due to allosteric modulation and reveal associated changes in ligand dissociation pathways.

Targeted Free Energy Perturbation Revisited: Accurate Free Energies from Mapped Reference Potentials

We present an approach that extends the theory of targeted free energy perturbation (TFEP) to calculate free energy differences and free energy surfaces at an accurate quantum mechanical level of theory from a cheaper reference potential. The convergence is accelerated by a mapping function that increases the overlap between the target and the reference distributions. Building on recent work, we show that this map can be learned with a normalizing flow neural network, without requiring simulations with the expensive target potential but only a small number of single-point calculations, and, crucially, avoiding the systematic error that was found previously. We validate the method by numerically evaluating the free energy difference in a system with a double-well potential and by describing the free energy landscape of a simple chemical reaction in the gas phase.

Probing the Suitability of Different Ca2+ Parameters for Long Simulations of Diisopropyl Fluorophosphatase

Organophosphate hydrolases are promising as potential biotherapeutic agents to treat poisoning with pesticides or nerve gases. However, these enzymes often need to be further engineered in order to become useful in practice. One example of such enhancement is the alteration of enantioselectivity of diisopropyl fluorophosphatase (DFPase). Molecular modeling techniques offer a unique opportunity to address this task rationally by providing a physical description of the substrate-binding process. However, DFPase is a metalloenzyme, and correct modeling of metal cations is a challenging task generally coming with a tradeoff between simulation speed and accuracy. Here, we probe several molecular mechanical parameter combinations for their ability to empower long simulations needed to achieve a quantitative description of substrate binding. We demonstrate that a combination of the Amber19sb force field with the recently developed 12-6 Ca2+ models allows us to both correctly model DFPase and obtain new insights into the DFP binding process.

The two redox states of the human NEET proteins' [2Fe-2S] clusters

The NEET proteins constitute a unique class of [2Fe–2S] proteins. The metal ions bind to three cysteines and one histidine. The proteins’ clusters exist in two redox states; the oxidized protein (containing two Fe^III ions) can transfer the cluster to apo-acceptor protein(s), while the reduced form (containing one ferrous ion) remains bound to the protein frame. Here, we perform in silico and in vitro studies on human NEET proteins in both reduced and oxidized forms. Quantum chemical calculations on all available human NEET proteins structures suggest that reducing the cluster weakens the Fe–N^His and Fe–S^Cys bonds, similar to what is seen in other Fe–S proteins (e.g., ferredoxin and Rieske protein). We further show that the extra electron in the [2Fe–2S]⁺ clusters of one of the NEET proteins (mNT) is localized on the His-bound iron ion, consistently with our previous spectroscopic studies. Kinetic measurements demonstrate that the mNT [2Fe–2S]⁺ is released only by an increase in temperature. Thus, the reduced state of human NEET proteins [2Fe–2S] cluster is kinetically inert. This previously unrecognized kinetic inertness of the reduced state, along with the reactivity of the oxidized state, is unique across all [2Fe–2S] proteins. Finally, using a coevolutionary analysis, along with molecular dynamics simulations, we provide insight on the observed allostery between the loop L2 and the cluster region. Specifically, we show that W75, R76, K78, K79, F82 and G85 in the latter region share similar allosteric characteristics in both redox states.

An Appraisal of Computer Simulation Approaches in Elucidating Biomolecular Recognition Pathways

Computer simulation approaches in biomolecular recognition processes have come a long way. In this Perspective, we highlight a series of recent success stories in which computer simulations have played a remarkable role in elucidating the atomic resolution mechanism of kinetic processes of protein-ligand binding in a quantitative fashion. In particular, we show that a robust combination of unbiased simulation, harnessed by a high-fidelity computing environment, and Markov state modeling approaches has been instrumental in revealing novel protein-ligand recognition pathways in multiple systems. We also elucidate the role of recent developments in enhanced sampling approaches in providing the much-needed impetus in accelerating simulation of the ligand recognition process. We identify multiple key issues, including force fields and the sampling bottleneck, which are currently preventing the field from achieving quantitative reconstruction of experimental measurements. Finally, we suggest a possible way forward via adoption of multiscale approaches and coarse-grained simulations as next steps toward efficient elucidation of ligand binding kinetics.

Membrane-Mediated Ligand Unbinding of the PK-11195 Ligand from TSPO

The translocator protein (TSPO), previously known as the peripheral benzodiazepine receptor, is of longstanding medical interest as both a biomarker for neuroinjury and a potential drug target for neuroinflammation and other disorders. Recently, it was shown that ligand residence time is a key factor determining steroidogenic efficacy of TSPO-binding compounds. This spurs interest in simulations of (un)binding pathways of TSPO ligands, which could reveal the molecular interactions governing ligand residence time. In this study, we use a weighted ensemble algorithm to determine the unbinding pathway for different poses of PK-11195, a TSPO ligand used in neuroimaging. In contrast with previous studies, our results show that PK-11195 does not dissociate directly into the solvent but instead dissociates via the lipid membrane by going between the transmembrane helices. We analyze this path ensemble in detail, constructing descriptors that can facilitate a general understanding of membrane-mediated ligand binding. We construct a set of Markov state models augmented with additional straightforward simulations to determine pose-specific ligand residence times. Together, we combine over 40 μs of trajectory data to form a coherent picture of the ligand binding landscape. We find that multiple starting poses yield residence times that roughly agree with the experimental quantity. The ligand binding transition states predicted by these Markov state models occur when PK-11195 is already in the membrane and involves only minimal ligand-protein interactions. This has implications for the design of new long-residence-time TSPO ligands.

A combination of machine learning and infrequent metadynamics to efficiently predict kinetic rates, transition states, and molecular determinants of drug dissociation from G protein-coupled receptors

Determining the drug-target residence time (RT) is of major interest in drug discovery given that this kinetic parameter often represents a better indicator of in vivo drug efficacy than binding affinity. However, obtaining drug-target unbinding rates poses significant challenges, both computationally and experimentally. This is particularly palpable for complex systems like G Protein-Coupled Receptors (GPCRs) whose ligand unbinding typically requires very long timescales oftentimes inaccessible by standard molecular dynamics simulations. Enhanced sampling methods offer a useful alternative, and their efficiency can be further improved by using machine learning tools to identify optimal reaction coordinates. Here, we test the combination of two machine learning techniques, automatic mutual information noise omission and reweighted autoencoded variational Bayes for enhanced sampling, with infrequent metadynamics to efficiently study the unbinding kinetics of two classical drugs with different RTs in a prototypic GPCR, the μ-opioid receptor. Dissociation rates derived from these computations are within one order of magnitude from experimental values. We also use the simulation data to uncover the dissociation mechanisms of these drugs, shedding light on the structures of rate-limiting transition states, which, alongside metastable poses, are difficult to obtain experimentally but important to visualize when designing drugs with a desired kinetic profile.