Kinetics of Huperzine A Dissociation from Acetylcholinesterase via Multiple Unbinding Pathways

How Robust Is the Ligand Binding Transition State?

For many drug targets, it has been shown that the kinetics of drug binding (e.g., on rate and off rate) is more predictive of drug efficacy than thermodynamic quantities alone. This motivates the development of predictive computational models that can be used to optimize compounds on the basis of their kinetics. The structural details underpinning these computational models are found not only in the bound state but also in the short-lived ligand binding transition states. Although transition states cannot be directly observed experimentally due to their extremely short lifetimes, recent successes have demonstrated that modeling the ligand binding transition state is possible with the help of enhanced sampling molecular dynamics methods. Previously, we generated unbinding paths for an inhibitor of soluble epoxide hydrolase (sEH) with a residence time of 11 min. Here, we computationally modeled unbinding events with the weighted ensemble method REVO (resampling of ensembles by variation optimization) for five additional inhibitors of sEH with residence times ranging from 14.25 to 31.75 min, with average prediction accuracy within an order of magnitude. The unbinding ensembles are analyzed in detail, focusing on features of the ligand binding transition state ensembles (TSEs). We find that ligands with similar bound poses can show significant differences in their ligand binding TSEs, in terms of their spatial distribution and protein-ligand interactions. However, we also find similarities across the TSEs when examining more general features such as ligand degrees of freedom. Together these findings show significant challenges for rational, kinetics-based drug design.

Toward Overcoming Pyrethroid Resistance in Mosquito Control: The Role of Sodium Channel Blocker Insecticides

Diseases spread by mosquitoes lead to the death of 700,000 people each year. The main way to reduce transmission is vector control by biting prevention with chemicals. However, the most commonly used insecticides lose efficacy due to the growing resistance. Voltage-gated sodium channels (VGSCs), membrane proteins responsible for the depolarizing phase of an action potential, are targeted by a broad range of neurotoxins, including pyrethroids and sodium channel blocker insecticides (SCBIs). Reduced sensitivity of the target protein due to the point mutations threatened malaria control with pyrethroids. Although SCBIs-indoxacarb (a pre-insecticide bioactivated to DCJW in insects) and metaflumizone-are used in agriculture only, they emerge as promising candidates in mosquito control. Therefore, a thorough understanding of molecular mechanisms of SCBIs action is urgently needed to break the resistance and stop disease transmission. In this study, by performing an extensive combination of equilibrium and enhanced sampling molecular dynamics simulations (3.2 μs in total), we found the DIII-DIV fenestration to be the most probable entry route of DCJW to the central cavity of mosquito VGSC. Our study revealed that F1852 is crucial in limiting SCBI access to their binding site. Our results explain the role of the F1852T mutation found in resistant insects and the increased toxicity of DCJW compared to its bulkier parent compound, indoxacarb. We also delineated residues that contribute to both SCBIs and non-ester pyrethroid etofenprox binding and thus could be involved in the target site cross-resistance.

Selecting High-Dimensional Representations of Physical Systems by Reweighted Diffusion Maps

Constructing reduced representations of high-dimensional systems is a fundamental problem in physical chemistry. Many unsupervised machine learning methods can automatically find such low-dimensional representations. However, an often overlooked problem is what high-dimensional representation should be used to describe systems before dimensionality reduction. Here, we address this issue using a recently developed method called the reweighted diffusion map [J. Chem. Theory Comput. 2022, 18, 7179-7192]. We show how high-dimensional representations can be quantitatively selected by exploring the spectral decomposition of Markov transition matrices built from data obtained from standard or enhanced sampling atomistic simulations. We demonstrate the performance of the method in several high-dimensional examples.

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.

Reweighted Manifold Learning of Collective Variables from Enhanced Sampling Simulations

Enhanced sampling methods are indispensable in computational chemistry and physics, where atomistic simulations cannot exhaustively sample the high-dimensional configuration space of dynamical systems due to the sampling problem. A class of such enhanced sampling methods works by identifying a few slow degrees of freedom, termed collective variables (CVs), and enhancing the sampling along these CVs. Selecting CVs to analyze and drive the sampling is not trivial and often relies on chemical intuition. Despite routinely circumventing this issue using manifold learning to estimate CVs directly from standard simulations, such methods cannot provide mappings to a low-dimensional manifold from enhanced sampling simulations, as the geometry and density of the learned manifold are biased. Here, we address this crucial issue and provide a general reweighting framework based on anisotropic diffusion maps for manifold learning that takes into account that the learning data set is sampled from a biased probability distribution. We consider manifold learning methods based on constructing a Markov chain describing transition probabilities between high-dimensional samples. We show that our framework reverts the biasing effect, yielding CVs that correctly describe the equilibrium density. This advancement enables the construction of low-dimensional CVs using manifold learning directly from the data generated by enhanced sampling simulations. We call our framework reweighted manifold learning. We show that it can be used in many manifold learning techniques on data from both standard and enhanced sampling simulations.

Computationally Elucidating the Binding Kinetics for Different AChE Inhibitors to Access the Rationale for Improving the Drug Efficacy

Traditional drug discovery is based on a binding affinity (thermodynamics)-driven paradigm. Numerous examples, however, demonstrated that drug efficacy does not always depend only on binding affinity but positively correlates with binding kinetics, that is, the dissociation rate constant (k_off). Binding free energy landscape (BFEL) constructor is a computational binding kinetics prediction method, previously developed by us, that estimates the binding kinetics for ligand-protein based on their constructed binding free energy landscape, but it also reveals the detailed molecular mechanism of the binding event, hence, providing the position of transition states at the molecular level to modify/improve the binding kinetics. Acetylcholinesterase (AChE) is a well-known Alzheimer's disease (AD) target for which there is still not an ideal drug on the market. Therefore, to improve the drug design strategy for AD, the binding kinetics and binding molecular mechanisms of the four inhibitors of AChE, that is, E2020 (Aricept), HupA, Rivastigmine, and Galantamine, were studied. Also, the differentiation of the binding kinetics between mAChE and TcAChE was studied to evaluate the sensitiveness of BFEL constructor. The flexibility of molecules has a noticeable effect on the nature of BFEL. To the same target, flexible molecules (i.e., E2020 and Rivastigmine) which contain more rotatable bonds tend to have more complicated BFELs reflecting more complicated molecular action mechanisms than the rigid ones (i.e., HupA and Galantamine), which therefore could be more challenging to be optimized. The binding kinetics is highly dependent on the structure of the molecules, such as the length and the functional groups. Therefore, E2020 presents better binding kinetic and thermodynamic properties with either TcAChE or mAChE. Therefore, it is the most promising lead drug for binding kinetics-based drug design. In addition, the binding kinetics of a drug may present different values in the proteins of different organisms because the residue compositions of the binding gorges of the targets are variant, that is, E2020 shows lower binding affinity and association energy barrier in binding with mAChE than TcAChE. However, HupA presents a better binding property with TcAChE than mAChE.

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.

Enhancing the Inhomogeneous Photodynamics of Canonical Bacteriophytochrome

The ability of phytochromes to act as photoswitches in plants and microorganisms depends on interactions between a bilin-like chromophore and a host protein. The interconversion occurs between the spectrally distinct red (Pr) and far-red (Pfr) conformers. This conformational change is triggered by the photoisomerization of the chromophore D-ring pyrrole. In this study, as a representative example of a phytochrome-bilin system, we consider biliverdin IXα (BV) bound to bacteriophytochrome (BphP) from Deinococcus radiodurans. In the absence of light, we use an enhanced sampling molecular dynamics (MD) method to overcome the photoisomerization energy barrier. We find that the calculated free energy (FE) barriers between essential metastable states agree with spectroscopic results. We show that the enhanced dynamics of the BV chromophore in BphP contributes to triggering nanometer-scale conformational movements that propagate by two experimentally determined signal transduction pathways. Most importantly, we describe how the metastable states enable a thermal transition known as the dark reversion between Pfr and Pr, through a previously unknown intermediate state of Pfr. We present the heterogeneity of temperature-dependent Pfr states at the atomistic level. This work paves a way toward understanding the complete mechanism of the photoisomerization of a bilin-like chromophore in phytochromes.

Comparative analysis of the unbinding pathways of antiviral drug Indinavir from HIV and HTLV1 proteases by supervised molecular dynamics simulation

Determining the unbinding pathways of potential small molecule compounds from their target proteins is of great significance for designing efficacious treatment solutions. One of these potential compounds is the approved HIV-1 protease inhibitor, Indinavir, which has a weak effect on the HTLV-1 protease. In this work, by employing the SuMD method, we reconstructed the unbinding pathways of Indinavir from HIV and HTLV-1 proteases to compare and understand the mechanism of the unbinding and to discover the reasons for the lack of inhibitory activity of Indinavir against the HTLV-1 protease. We achieved multiple unbinding events from both HIV and HTLV-1 proteases in which the RMSD values of Indinavir reached over 40 Å. Also, we found that the mobility and fluctuations of the flap region are higher in the HTLV-1 protease, making the drug less stable. We realized that critically positioned aromatic residues such as Trp98/Trp98' and Phe67/Phe67' in the HTLV-1 protease could make strong π-Stacking interactions with Indinavir in the unbinding pathway, which are unfavorable for the stability of Indinavir in the active site. The details found in this study can make a reasonable explanation for the lack of inhibitory activity of this drug against HTLV-1 protease. We believe the details discovered in this work can help design more effective and selective inhibitors for the HTLV-1 protease.

A Comprehensive Review of Cholinesterase Modeling and Simulation

The present article reviews published efforts to study acetylcholinesterase and butyrylcholinesterase structure and function using computer-based modeling and simulation techniques. Structures and models of both enzymes from various organisms, including rays, mice, and humans, are discussed to highlight key structural similarities in the active site gorges of the two enzymes, such as flexibility, binding site location, and function, as well as differences, such as gorge volume and binding site residue composition. Catalytic studies are also described, with an emphasis on the mechanism of acetylcholine hydrolysis by each enzyme and novel mutants that increase catalytic efficiency. The inhibitory activities of myriad compounds have been computationally assessed, primarily through Monte Carlo-based docking calculations and molecular dynamics simulations. Pharmaceutical compounds examined herein include FDA-approved therapeutics and their derivatives, as well as several other prescription drug derivatives. Cholinesterase interactions with both narcotics and organophosphate compounds are discussed, with the latter focusing primarily on molecular recognition studies of potential therapeutic value and on improving our understanding of the reactivation of cholinesterases that are bound to toxins. This review also explores the inhibitory properties of several other organic and biological moieties, as well as advancements in virtual screening methodologies with respect to these enzymes.

Data-guided Multi-Map variables for ensemble refinement of molecular movies

Driving molecular dynamics simulations with data-guided collective variables offer a promising strategy to recover thermodynamic information from structure-centric experiments. Here, the three-dimensional electron density of a protein, as it would be determined by cryo-EM or x-ray crystallography, is used to achieve simultaneously free-energy costs of conformational transitions and refined atomic structures. Unlike previous density-driven molecular dynamics methodologies that determine only the best map-model fits, our work employs the recently developed Multi-Map methodology to monitor concerted movements within equilibrium, non-equilibrium, and enhanced sampling simulations. Construction of all-atom ensembles along the chosen values of the Multi-Map variable enables simultaneous estimation of average properties, as well as real-space refinement of the structures contributing to such averages. Using three proteins of increasing size, we demonstrate that biased simulation along the reaction coordinates derived from electron densities can capture conformational transitions between known intermediates. The simulated pathways appear reversible with minimal hysteresis and require only low-resolution density information to guide the transition. The induced transitions also produce estimates for free energy differences that can be directly compared to experimental observables and population distributions. The refined model quality is superior compared to those found in the Protein Data Bank. We find that the best quantitative agreement with experimental free-energy differences is obtained using medium resolution density information coupled to comparatively large structural transitions. Practical considerations for probing the transitions between multiple intermediate density states are also discussed.

Slow-binding inhibitors of acetylcholinesterase of medical interest

Certain ligands slowly bind to acetylcholinesterase. As a result, there is a slow establishment of enzyme-inhibitor equilibrium characterized by a slow onset of inhibition prior reaching steady state. Three mechanisms account for slow-binding inhibition: a) slow binding rate constant k_on, b) slow ligand induced-fit following a fast binding step, c) slow conformational selection of an enzyme form. The slow equilibrium may be followed by a chemical step. This later that can be irreversible has been observed with certain alkylating agents and substrate transition state analogs. Slow-binding inhibitors present long residence times on target. This results in prolonged pharmacological or toxicological action. Through several well-known molecules (e.g. huperzine) and new examples (tocopherol, trifluoroacetophenone and a 6-methyluracil alkylammonium derivative), we show that slow-binding inhibitors of acetylcholinesterase are promising drugs for treatment of neurological diseases such as Alzheimer disease and myasthenia gravis. Moreover, they may be of interest for neuroprotection (prophylaxis) against organophosphorus poisoning. This article is part of the special issue entitled 'Acetylcholinesterase Inhibitors: From Bench to Bedside to Battlefield'.

Transient States and Barriers from Molecular Simulations and the Milestoning Theory: Kinetics in Ligand-Protein Recognition and Compound Design

This study presents a novel computational approach to study molecular recognition and binding kinetics for drug-like compounds dissociating from a flexible protein system. The intermediates and their free energy profile during ligand association and dissociation processes control ligand-protein binding kinetics and bring a more complete picture of ligand-protein binding. The method applied the milestoning theory to extract kinetics and thermodynamics information from running short classical molecular dynamics (MD) simulations for frames from a given dissociation path. High-dimensional ligand-protein motions (3N-6 degrees of freedom) during ligand dissociation were reduced by use of principal component modes for assigning more than 100 milestones, and classical MD runs were allowed to travel multiple milestones to efficiently obtain ensemble distribution of initial structures for MD simulations and estimate the transition time and rate during ligand traveling between milestones. We used five pyrazolourea ligands and cyclin-dependent kinase 8 with cyclin C (CDK8/CycC) as our model system as well as metadynamics and a pathway search method to sample dissociation pathways. With our strategy, we constructed the free energy profile for highly mobile biomolecular systems. The computed binding free energy and residence time correctly ranked the pyrazolourea ligand series, in agreement with experimental data. Guided by a barrier of a ligand passing an αC helix and activation loop, we introduced one hydroxyl group to parent compounds to design our ligands with increased residence time and validated our prediction by experiments. This work provides a novel and robust approach to investigate dissociation kinetics of large and flexible systems for understanding unbinding mechanisms and designing new small-molecule drugs with desired binding kinetics.

Trans-toxin ion-sensitivity of charybdotoxin-blocked potassium-channels reveals unbinding transitional states

In silico and in vitro studies have made progress in understanding protein–protein complex formation; however, the molecular mechanisms for their dissociation are unclear. Protein–protein complexes, lasting from microseconds to years, often involve induced-fit, challenging computational or kinetic analysis. Charybdotoxin (CTX), a peptide from the Leiurus scorpion venom, blocks voltage-gated K⁺-channels in a unique example of binding/unbinding simplicity. CTX plugs the external mouth of K⁺-channels pore, stopping K⁺-ion conduction, without inducing conformational changes. Conflicting with a tight binding, we show that external permeant ions enhance CTX-dissociation, implying a path connecting the pore, in the toxin-bound channel, with the external solution. This sensitivity is explained if CTX wobbles between several bound conformations, producing transient events that restore the electrical and ionic trans-pore gradients. Wobbling may originate from a network of contacts in the interaction interface that are in dynamic stochastic equilibria. These partially-bound intermediates could lead to distinct, and potentially manipulable, dissociation pathways.

Dissociation Process of a MDM2/p53 Complex Investigated by Parallel Cascade Selection Molecular Dynamics and the Markov State Model

Recently, we efficiently generated dissociation pathways of a protein-ligand complex without applying force bias with parallel cascade selection molecular dynamics (PaCS-MD) and showed that PaCS-MD in combination with the Markov state model (MSM) yielded a binding free energy comparable to experimental values. In this work, we applied the same procedure to a complex of MDM2 protein and the transactivation domain of p53 protein (TAD-p53), the latter of which is known to be very flexible in the unbound state. Using 30 independent MD simulations in PaCS-MD, we successfully generated 25 dissociation pathways of the complex, which showed complete or partial unfolding of the helical region of TAD-p53 during the dissociation process within an average simulation time of 154.8 ± 46.4 ns. The standard binding free energy obtained in combination with one-dimensional-, three-dimensional (3D)- or C_α-MSM was in good agreement with those determined experimentally. Using 3D-MSM based on the center of mass position of TAD-p53 relative to MDM2, the dissociation rate constant was calculated, which was comparable to those measured experimentally. C_α-MSM based on all C_α coordinates of TAD-p53 reproduced the experimentally measured standard binding free energy, and dissociation and association rate constants. We conclude that the combination of PaCS-MD and MSM offers an efficient computational procedure to calculate binding free energies and kinetic rates.