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

Increased throughput in methods for simulating protein ligand binding and unbinding

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

Measuring the Adsorption Cross Section of YOYO-1 to Immobilized DNA Molecules

Many interactions between small molecules and particles occur in solutions. They are surrounded by other molecules that do not react, for example, biological processes in water, chemical reactions in gas or liquid solutions, and environmental reactions in air and water. However, predicting the rate of such diffusive interactions remains challenging, due to the random motion of molecules in solutions, as exampled by the famous Brownian motion of pollen particles. In this report, we experimentally confirmed that a disruptive rate equation we have published before can predict the association rate of typical adsorption at interfaces, which has a surprising fraction order of 4/3 that has not been considered before. This could be an important step toward a generalized method to predict the adsorption rate of many reactions.

Incorporating Prior Knowledge in the Seeds of Adaptive Sampling Molecular Dynamics Simulations of Ligand Transport in Enzymes with Buried Active Sites

Because most proteins have buried active sites, protein tunnels or channels play a crucial role in the transport of small molecules into buried cavities for enzymatic catalysis. Tunnels can critically modulate the biological process of protein-ligand recognition. Various molecular dynamics methods have been developed for exploring and exploiting the protein-ligand conformational space to extract high-resolution details of the binding processes, a recent example being energetically unbiased high-throughput adaptive sampling simulations. The current study systematically contrasted the role of integrating prior knowledge while generating useful initial protein-ligand configurations, called seeds, for these simulations. Using a nontrivial system of a haloalkane dehalogenase mutant with multiple transport tunnels leading to a deeply buried active site, simulations were employed to derive kinetic models describing the process of association and dissociation of the substrate molecule. The most knowledge-based seed generation enabled high-throughput simulations that could more consistently capture the entire transport process, explore the complex network of transport tunnels, and predict equilibrium dissociation constants, koff/kon, on the same order of magnitude as experimental measurements. Overall, the infusion of more knowledge into the initial seeds of adaptive sampling simulations could render analyses of transport mechanisms in enzymes more consistent even for very complex biomolecular systems, thereby promoting drug development efforts and the rational design of enzymes with buried active sites.

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

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

Characterization of the Bottlenecks and Pathways for Inhibitor Dissociation from [NiFe] Hydrogenase

[NiFe] hydrogenases can act as efficient catalysts for hydrogen oxidation and biofuel production. However, some [NiFe] hydrogenases are inhibited by gas molecules present in the environment, such as O2 and CO. One strategy to engineer [NiFe] hydrogenases and achieve O2- and CO-tolerant enzymes is by introducing point mutations to block the access of inhibitors to the catalytic site. In this work, we characterized the unbinding pathways of CO in the complex with the wild-type and 10 different mutants of [NiFe] hydrogenase from Desulfovibrio fructosovorans using τ-random accelerated molecular dynamics (τRAMD) to enhance the sampling of unbinding events. The ranking provided by the relative residence times computed with τRAMD is in agreement with experiments. Extensive data analysis of the simulations revealed that from the two bottlenecks proposed in previous studies for the transit of gas molecules (residues 74 and 122 and residues 74 and 476), only one of them (residues 74 and 122) effectively modulates diffusion and residence times for CO. We also computed pathway probabilities for the unbinding of CO, O2, and H2 from the wild-type [NiFe] hydrogenase, and we observed that while the most probable pathways are the same, the secondary pathways are different. We propose that introducing mutations to block the most probable paths, in combination with mutations to open the main secondary path used by H2, can be a feasible strategy to achieve CO and O2 resistance in the [NiFe] hydrogenase from Desulfovibrio fructosovorans.

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

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

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.

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.

Integrating network pharmacology with molecular docking for elucidation of molecular biological mechanisms of Jiedu Qingjin formula for non-small cell lung cancer

Traditional Chinese medicine is an important part of complementary alternative medicine. Jiedu Qingjin formula (JDQJF) is an effective national invention patent for the treatment of non-small cell lung cancer (NSCLC). We investigated the molecular biological mechanisms based on network pharmacology, molecular docking, and molecular dynamics simulations. Compounds of JDQJF were screened through the TCMSP, ETCM, and literature. Targets were searched by DrugBank and predicted by SwissTargetPrediction. GEO database was applied for screening differentially expressed genes between cancerous tissues and healthy tissues of NSCLC. Subsequently, the protein-protein interaction between JDQJF and NSCLC were obtained by Cytoscape. Visual analyses were carried out to extract candidate genes, then subjected to Metascape for enrichment analyses. Finally, molecular docking was performed by AutoDock, and the best complexes were subjected to molecular dynamics simulation and binding energy calculations by MMPBSA. A total of 273 compounds, 390 targets, 3146 GO terms, and 174 KEGG pathways were obtained. Five potential compounds (quercetin, adenosine, apigenin, heptadecanoic acid, and luteolin) were notably modulated by key targets AKT1, MAPK3, and RAF1. Enrichment results included cell cycle process, growth transduction factor, immune response-activating transduction, and involved PI3K/AKT, MAPK, NF-κB and VEGF pathway. RAF1-quercetin showed the highest binding affinity (-9.1 kcal/mol), revealed stable interactions during the simulation, and the highest estimated relative binding energy of the RAF1-Heptadecanoic was -184.277 kcal/mol. This study suggested that EMT-related, inflammation-related, immune-related, and angiogenesis-related pathways may be associated with JDQJF, and involved in the advancement of NSCLC, which points out the research direction for subsequent utility mechanism validation.Communicated by Ramaswamy H. Sarma.

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.

Application of Molecular Dynamics Simulations to Determine Interactions between Canary Seed (Phalaris canariensis L.) Bioactive Peptides and Skin-Aging Enzymes

Food bioactive peptides are well recognized for their health benefits such as antimicrobial, antioxidant, and antihypertensive benefits, among others. Their drug-like behavior has led to their potential use in targeting skin-related aging factors like the inhibition of enzymes related with the skin-aging process. In this study, canary seed peptides (CSP) after simulated gastrointestinal digestion (<3 kDa) were fractioned by RP-HPLC and their enzyme-inhibition activity towards elastase and tyrosinase was evaluated in vitro. CSP inhibited elastase (IC50 = 6.2 mg/mL) and tyrosinase (IC50 = 6.1 mg/mL), while the hydrophobic fraction-VI (0.2 mg/mL) showed the highest inhibition towards elastase (93%) and tyrosinase (67%). The peptide fraction with the highest inhibition was further characterized by a multilevel in silico workflow, including physicochemical descriptor calculations, antioxidant activity predictions, and molecular dynamics-ensemble docking towards elastase and tyrosinase. To gain insights into the skin permeation process during molecular dynamics simulations, based on their docking scores, five peptides (GGWH, VPPH, EGLEPNHRVE, FLPH, and RPVNKYTPPQ) were identified to have favorable intermolecular interactions, such as hydrogen bonding of polar residues (W, H, and K) to lipid polar groups and 2-3 Å van der Waals close contact of hydrophobic aliphatic residues (P, V, and L). These interactions can play a critical role for the passive insertion of peptides into stratum corneum model skin-membranes, suggesting a promising application of CSP for skin-aging treatments.