Exhaustive Search of Ligand Binding Pathways via Volume-Based Metadynamics

Reinforcing Tunnel Network Exploration in Proteins Using Gaussian Accelerated Molecular Dynamics

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

Structural basis of the obligatory exchange mode of human neutral amino acid transporter ASCT2

ASCT2 is an obligate exchanger of neutral amino acids, contributing to cellular amino acid homeostasis. ASCT2 belongs to the same family (SLC1) as Excitatory Amino Acid Transporters (EAATs) that concentrate glutamate in the cytosol. The mechanism that makes ASCT2 an exchanger rather than a concentrator remains enigmatic. Here, we employ cryo-electron microscopy and molecular dynamics simulations to elucidate the structural basis of the exchange mechanism of ASCT2. We establish that ASCT2 binds three Na+ ions per transported substrate and visits a state that likely acts as checkpoint in preventing Na+ ion leakage, both features shared with EAATs. However, in contrast to EAATs, ASCT2 retains one Na+ ion even under Na+-depleted conditions. We demonstrate that ASCT2 cannot undergo the structural transition in TM7 that is essential for the concentrative transport cycle of EAATs. This structural rigidity and the high-affinity Na+ binding site effectively confine ASCT2 to an exchange mode.

Mechanism of phosphate release from actin filaments

After ATP-actin monomers assemble filaments, the ATP's [Formula: see text]-phosphate is hydrolyzedwithin seconds and dissociates over minutes. We used all-atom molecular dynamics simulations to sample the release of phosphate from filaments and study residues that gate release. Dissociation of phosphate from Mg2+ is rate limiting and associated with an energy barrier of 20 kcal/mol, consistent with experimental rates of phosphate release. Phosphate then diffuses within an internal cavity toward a gate formed by R177, as suggested in prior computational studies and cryo-EM structures. The gate is closed when R177 hydrogen bonds with N111 and is open when R177 forms a salt bridge with D179. Most of the time, interactions of R177 with other residues occlude the phosphate release pathway. Machine learning analysis reveals that the occluding interactions fluctuate rapidly, underscoring the secondary role of backdoor gate opening in Pi release, in contrast with the previous hypothesis that gate opening is the primary event.

Machine Learning Molecular Dynamics Shows Anomalous Entropic Effect on Catalysis through Surface Pre-melting of Nanoclusters

Due to the superior catalytic activity and efficient utilization of noble metals, nanocatalysts are extensively used in the modern industrial production of chemicals. The surface structures of these materials are significantly influenced by reactive adsorbates, leading to dynamic behavior under experimental conditions. The dynamic nature poses significant challenges in studying the structure-activity relations of catalysts. Herein, we unveil an anomalous entropic effect on catalysis via surface pre-melting of nanoclusters through machine learning accelerated molecular dynamics and free energy calculation. We find that due to the pre-melting of shell atoms, there exists a non-linear variation in the catalytic activity of the nanoclusters with temperature. Consequently, two notable changes in catalyst activity occur at the respective temperatures of melting for the shell and core atoms. We further study the nanoclusters with surface point defects, i.e. vacancy and ad-atom, and observe significant decrease in the surface melting temperatures of the nanoclusters, enabling the reaction to take place under more favorable and milder conditions. These findings not only provide novel insights into dynamic catalysis of nanoclusters but also offer new understanding of the role of point defects in catalytic processes.

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.

Understanding the Differences of Danusertib's Residence Time in Aurora Kinases A/B: Dissociation Paths and Key Residues Identified using Conventional and Enhanced Molecular Dynamics Simulations

The accurate experimental estimation of protein-ligand systems' residence time (τ) has become very relevant in drug design projects due to its importance in the last stages of refinement of the drug's pharmacodynamics and pharmacokinetics. It is now well-known that it is not sufficient to estimate the affinity of a protein-drug complex in the thermodynamic equilibrium process in in vitro experiments (closed systems), where the concentrations of the drug and protein remain constant. On the contrary, it is mandatory to consider the conformational dynamics of the system in terms of the binding and unbinding processes between protein and drugs in in vivo experiments (open systems), where their concentrations are in constant flux. This last model has been proven to dictate much of several drugs' pharmacological activities in vivo. At the atomistic level, molecular dynamics simulations can explain why some drugs are more effective than others or unveil the molecular aspects that make some drugs work better in one molecular target. Here, the protein kinases Aurora A/B, complexed with its inhibitor Danusertib, were studied using conventional and enhanced molecular dynamics (MD) simulations to estimate the dissociation paths and, therefore, the computational τ values and their comparison with experimental ones. Using classical molecular dynamics (cMD), three differential residues within the Aurora A/B active site, which seems to play an essential role in the observed experimental Danusertib's residence time against these kinases, were characterized. Then, using WT-MetaD, the relative Danusertib's residence times against Aurora A/B kinases were measured in a nanosecond time scale and were compared to those τ values observed experimentally. In addition, the potential dissociation paths of Danusertib in Aurora A and B were characterized, and differences that might be explained by the differential residues in the enzyme's active sites were found. In perspective, it is expected that this computational protocol can be applied to other protein-ligand complexes to understand, at the molecular level, the differences in residence times and amino acids that may contribute to it.

Mechanism of Phosphate Release from Actin Filaments

After ATP-actin monomers assemble filaments, the ATP's γ-phosphate is hydrolyzed within seconds and dissociates over minutes. We used all-atom molecular dynamics simulations to sample the release of phosphate from filaments and study residues that gate release. Dissociation of phosphate from Mg2+ is rate limiting and associated with an energy barrier of 20 kcal/mol, consistent with experimental rates of phosphate release. Phosphate then diffuses in an internal cavity toward a gate formed by R177 suggested in prior computational studies and cryo-EM structures. The gate is closed when R177 hydrogen bonds with N111 and is open when R177 forms a salt bridge with D179. Most of the time interactions of R177 with other residues occludes the phosphate release pathway. Machine learning analysis reveals that the occluding interactions fluctuate rapidly, underscoring the secondary role of backdoor gate opening in Pi release, in contrast with the previous hypothesis that gate opening is the primary event.

Exploring Product Release from Yeast Cytosine Deaminase with Metadynamics

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

Data-driven classification of ligand unbinding pathways

Studying the pathways of ligand-receptor binding is essential to understand the mechanism of target recognition by small molecules. The binding free energy and kinetics of protein-ligand complexes can be computed using molecular dynamics (MD) simulations, often in quantitative agreement with experiments. However, only a qualitative picture of the ligand binding/unbinding paths can be obtained through a conventional analysis of the MD trajectories. Besides, the higher degree of manual effort involved in analyzing pathways limits its applicability in large-scale drug discovery. Here, we address this limitation by introducing an automated approach for analyzing molecular transition paths with a particular focus on protein-ligand dissociation. Our method is based on the dynamic time-warping algorithm, originally designed for speech recognition. We accurately classified molecular trajectories using a very generic descriptor set of contacts or distances. Our approach outperforms manual classification by distinguishing between parallel dissociation channels, within the pathways identified by visual inspection. Most notably, we could compute exit-path-specific ligand-dissociation kinetics. The unbinding timescale along the fastest path agrees with the experimental residence time, providing a physical interpretation to our entirely data-driven protocol. In combination with appropriate enhanced sampling algorithms, this technique can be used for the initial exploration of ligand-dissociation pathways as well as for calculating path-specific thermodynamic and kinetic properties.

Assessing the Binding Mode of a Splicing Modulator Stimulating Pre-mRNA Binding to the Plastic U2AF2 Splicing Factor

RNA recognition motifs (RRMs) play a pivotal role in RNA metabolism and the regulation of gene expression. Owing to their plasticity and fuzziness, targeting RRM/RNA interfaces with small molecules is a daunting challenge for drug discovery campaigns. The U2AF2 splicing factor, which recognizes the polypyrimidine (polyPy) sequence of premature messenger (pre-m)RNA, exhibits a dynamic architecture consisting of two RRMs joined by a disordered linker. An inhibitor, NSC-194308, was shown to enhance the binding of pre-mRNA to U2AF2, selectively triggering cell death in leukemia cell lines containing spliceosome mutations. The NSC-194308 binding mode remains elusive; yet, unraveling its knowledge may offer intriguing insights for effectively targeting U2AF2 and other flexible protein/protein/RNA interfaces with small molecules. To infer plausible NSC-194308 binding poses to U2AF2, here, we applied and benchmarked the performance of static and dynamic docking approaches, elucidating the molecular basis of the NSC-194308-induced pre-mRNA stabilization on U2AF2. We demonstrate that introducing dynamic effects is mandatory to assess the binding mode of the inhibitors when they target plastic and modular architectures, such as those formed by interacting RRMs. The latter are widespread across RNA binding proteins; therefore, this mechanism may be broadly applicable to discover new therapeutics aimed at selectively modulating the RNA function by targeting protein/protein/RNA interfaces.

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).

Differences in Conformational Sampling and Intrinsic Electric Fields Drive Ion Binding in Telomeric and TERRA G-Quadruplexes

The formation of G-quadruplexes (GQs) occurs in guanine-rich sequences of DNA and RNA, producing highly stable and structurally diverse noncanonical nucleic acid structures. GQs play crucial roles in regulating transcription, translation, and replication and maintaining the genome, among others; thus, changes to their structures can lead to diseases such as cancer. Previous studies using polarizable molecular dynamics simulations have shown differences in ion binding properties between telomeric and telomeric repeat-containing RNA GQs despite architectural similarities. Here, we used volume-based metadynamics and repulsive potential simulations in conjunction with polarizable force fields to quantify the impact of ion binding on the GQ dynamics and ion binding free energies. Furthermore, we describe how GQs exert electric fields on their surroundings to link dynamics with variations in the electronic structure. Our findings provide new insights into the energetic, physical, and conformational properties of GQs and expose subtle but important differences between DNA and RNA GQs with the same fold.

Free energy along drug-protein binding pathways interactively sampled in virtual reality

We describe a two-step approach for combining interactive molecular dynamics in virtual reality (iMD-VR) with free energy (FE) calculation to explore the dynamics of biological processes at the molecular level. We refer to this combined approach as iMD-VR-FE. Stage one involves using a state-of-the-art 'human-in-the-loop' iMD-VR framework to generate a diverse range of protein-ligand unbinding pathways, benefitting from the sophistication of human spatial and chemical intuition. Stage two involves using the iMD-VR-sampled pathways as initial guesses for defining a path-based reaction coordinate from which we can obtain a corresponding free energy profile using FE methods. To investigate the performance of the method, we apply iMD-VR-FE to investigate the unbinding of a benzamidine ligand from a trypsin protein. The binding free energy calculated using iMD-VR-FE is similar for each pathway, indicating internal consistency. Moreover, the resulting free energy profiles can distinguish energetic differences between pathways corresponding to various protein-ligand conformations (e.g., helping to identify pathways that are more favourable) and enable identification of metastable states along the pathways. The two-step iMD-VR-FE approach offers an intuitive way for researchers to test hypotheses for candidate pathways in biomolecular systems, quickly obtaining both qualitative and quantitative insight.

Differences in Conformational Sampling and Intrinsic Electric Fields Drive Ion Binding in Telomeric and TERRA G-Quadruplexes

The formation of G-quadruplexes (GQs) occurs in guanine-rich sequences of DNA and RNA, producing highly stable and structurally diverse noncanonical nucleic acid structures. GQs play crucial roles in regulating transcription, translation, and replication; and maintaining the genome, among others, thus changes to their structures can lead to diseases such as cancer. Previous studies using polarizable molecular dynamics simulations have shown differences in ion binding properties between telomeric and TERRA GQs despite architectural similarities. Here, we used volume-based metadynamics and repulsive potential simulations in conjunction with polarizable force fields to quantify the impact of ion binding on GQ dynamics and ion binding free energies. Furthermore, we describe how GQs exert electric fields on their surroundings to link dynamics with variations in electronic structure. Our findings provide new insights into the energetic, physical, and conformational properties of GQs and expose subtle, but important, differences between DNA and RNA GQs with the same fold.

A deep encoder-decoder framework for identifying distinct ligand binding pathways

The pathway(s) that a ligand would adopt en route to its trajectory to the native pocket of the receptor protein act as a key determinant of its biological activity. While Molecular Dynamics (MD) simulations have emerged as the method of choice for modeling protein-ligand binding events, the high dimensional nature of the MD-derived trajectories often remains a barrier in the statistical elucidation of distinct ligand binding pathways due to the stochasticity inherent in the ligand's fluctuation in the solution and around the receptor. Here, we demonstrate that an autoencoder based deep neural network, trained using an objective input feature of a large matrix of residue-ligand distances, can efficiently produce an optimal low-dimensional latent space that stores necessary information on the ligand-binding event. In particular, for a system of L99A mutant of T4 lysozyme interacting with its native ligand, benzene, this deep encoder-decoder framework automatically identifies multiple distinct recognition pathways, without requiring user intervention. The intermediates involve the spatially discrete location of the ligand in different helices of the protein before its eventual recognition of native pose. The compressed subspace derived from the autoencoder provides a quantitatively accurate measure of the free energy and kinetics of ligand binding to the native pocket. The investigation also recommends that while a linear dimensional reduction technique, such as time-structured independent component analysis, can do a decent job of state-space decomposition in cases where the intermediates are long-lived, autoencoder is the method of choice in systems where transient, low-populated intermediates can lead to multiple ligand-binding pathways.

Metadynamics simulations of ligands binding to protein surfaces: a novel tool for rational drug design

Structure-based drug design protocols may encounter difficulties to investigate poses when the biomolecular targets do not exhibit typical binding pockets. In this study, by providing two concrete examples from our labs, we suggest that the combination of metadynamics free energy methods (validated against affinity measurements), along with experimental structural information (by X-ray crystallography and NMR), can help to identify the poses of ligands on protein surfaces. The simulation workflow proposed here was implemented in a widely used code, namely GROMACS, and it could straightforwardly be applied to various drug-design campaigns targeting ligands' binding to protein surfaces.

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.

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.

Resolving Protein Conformational Plasticity and Substrate Binding via Machine Learning

A long-standing target in elucidating the biomolecular recognition process is the identification of binding-competent conformations of the receptor protein. However, protein conformational plasticity and the stochastic nature of the recognition processes often preclude the assignment of a specific protein conformation to an individual ligand-bound pose. Here, we demonstrate that a computational framework coined as RF-TICA-MD, which integrates an ensemble decision-tree-based Random Forest (RF) machine learning (ML) technique with an unsupervised dimension reduction approach time-structured independent component analysis (TICA), provides an efficient and unambiguous solution toward resolving protein conformational plasticity and the substrate binding process. In particular, we consider multimicrosecond-long molecular dynamics (MD) simulation trajectories of a ligand recognition process in solvent-inaccessible cavities of archetypal proteins T4 lysozyme and cytochrome P450cam. We show that in a scenario in which clear correspondence between protein conformation and binding-competent macrostates could not be obtained via an unsupervised dimension reduction approach, an a priori decision-tree-based supervised classification of the simulated recognition trajectories via RF would help characterize key amino acid residue pairs of the protein that are deemed sensitive for ligand binding. A subsequent unsupervised dimensional reduction of the selected residue pairs via TICA would then delineate a conformational landscape of protein which is able to demarcate ligand-bound poses from unbound ones. The proposed RF-TICA-MD approach is shown to be data agnostic and found to be robust when using other ML-based classification methods such as XGBoost. As a promising spinoff of the protocol, the framework is found to be capable of identifying distal protein locations which would be allosterically important for ligand binding and would characterize their roles in recognition pathways. A Python implementation of a proposed ML workflow is available in GitHub https://github.com/navjeet0211/rf-tica-md.

Path separation of dissipation-corrected targeted molecular dynamics simulations of protein-ligand unbinding

Protein-ligand (un)binding simulations are a recent focus of biased molecular dynamics simulations. Such binding and unbinding can occur via different pathways in and out of a binding site. Here, we present a theoretical framework on how to compute kinetics along separate paths and on how to combine the path-specific rates into global binding and unbinding rates for comparison with experimental results. Using dissipation-corrected targeted molecular dynamics in combination with temperature-boosted Langevin equation simulations [S. Wolf et al., Nat. Commun. 11, 2918 (2020)] applied to a two-dimensional model and the trypsin-benzamidine complex as test systems, we assess the robustness of the procedure and discuss the aspects of its practical applicability to predict multisecond kinetics of complex biomolecular systems.

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.

Predictions of the Poses and Affinity of a Ligand over the Entire Surface of a NEET Protein: The Case of Human MitoNEET

Human NEET proteins contain two [2Fe-2S] iron-sulfur clusters, bound to three Cys residues and one His residue. They exist in two redox states. Recently, these proteins have revealed themselves as attractive drug targets for mitochondrial dysfunction-related diseases, such as type 2 diabetes, Wolfram syndrome 2, and cancers. Unfortunately, the lack of information and mechanistic understanding of ligands binding to the whole functional, cytoplasmatic domain has limited rational drug design approaches. Here, we use an enhanced sampling technique, volume-based metadynamics, recently developed by a team involving some of us, to predict the poses and affinity of the 2-benzamido-4-(1,2,3,4-tetrahydronaphthalen-2-yl)-thiophene-3-carboxylate ligand to the entire surface of the cytoplasmatic domain of the human NEET protein mitoNEET (mNT) in an aqueous solution. The calculations, based on the recently published X-ray structure of the complex, are consistent with the measured affinity. The calculated free energy landscape revealed that the ligand can bind in multiple sites and with poses other than the one found in the X-ray. This difference is likely to be caused by crystal packing effects that allow the ligand to interact with multiple adjacent NEET protein copies. Such extra contacts are of course absent in the solution; therefore, the X-ray pose is only transient in our calculations, where the binding free energy correlates with the number of contacts. We further evaluated how the reduction and protonation of the Fe-bound histidine, as well as temperature, can affect ligand binding. Both such modifications introduce the possibility for the ligand to bind in an area of the protein other than the one observed in the X-ray, with no or little impact on affinity. Overall, our study can provide insights on the molecular recognition mechanisms of ligand binding to mNT in different oxidative conditions, possibly helping rational drug design of NEET ligands.

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.

Role of Monovalent Ions in the NKCC1 Inhibition Mechanism Revealed through Molecular Simulations

The secondary active Na-K-Cl cotransporter 1 (NKCC1) promotes electroneutral uptake of two chloride ions, one sodium ion and one potassium ion. NKCC1 regulates Cl^- homeostasis, thus being implicated in transepithelial water transport and in neuronal excitability. Aberrant NKCC1 transport is linked to a variety of human diseases. The loop diuretic drugs bumetanide, furosemide, azosemide and ethacrynic acid target NKCC1, but are characterized by poor selectivity leading to severe side effects. Despite its therapeutic importance, the molecular details of the NKCC1 inhibition mechanism remain unclear. Using all-atom simulations, we predict a putative binding mode of these drugs to the zebrafish (z) and human (h) NKCC1 orthologs. Although differing in their specific interactions with NKCC1 and/or monovalent ions, all drugs can fit within the same cavity and engage in hydrophobic interactions with M304/M382 in z/hNKCC1, a proposed ion gating residue demonstrated to be key for bumetanide binding. Consistent with experimental evidence, all drugs take advantage of the K⁺/Na⁺ ions, which plastically respond to their binding. This study not only provides atomic-level insights useful for drug discovery campaigns of more selective/potent NKCC1 inhibitors aimed to tackle diseases related to deregulated Cl^- homeostasis, but it also supplies a paradigmatic example of the key importance of dynamical effects when drug binding is mediated by monovalent ions.

Ligand Unbinding Pathway and Mechanism Analysis Assisted by Machine Learning and Graph Methods

We present two methods to reveal protein-ligand unbinding mechanisms in biased unbinding simulations by clustering trajectories into ensembles representing unbinding paths. The first approach is based on a contact principal component analysis for reducing the dimensionality of the input data, followed by identification of unbinding paths and training a machine learning model for trajectory clustering. The second approach clusters trajectories according to their pairwise mean Euclidean distance employing the neighbor-net algorithm, which takes into account input data bias in the distances set and is superior to dendrogram construction. Finally, we describe a more complex case where the reaction coordinate relevant for path identification is a single intraligand hydrogen bond, highlighting the challenges involved in unbinding path reaction coordinate detection.

Interrogating RNA-Small Molecule Interactions with Structure Probing and Artificial Intelligence-Augmented Molecular Simulations

While there is increasing interest in the study of RNA as a therapeutic target, efforts to understand RNA-ligand recognition at the molecular level lag far behind our understanding of protein-ligand recognition. This problem is complicated due to the more than 10 orders of magnitude in time scales involved in RNA dynamics and ligand binding events, making it not straightforward to design experiments or simulations. Here, we make use of artificial intelligence (AI)-augmented molecular dynamics simulations to directly observe ligand dissociation for cognate and synthetic ligands from a riboswitch system. The site-specific flexibility profiles from our simulations are compared with in vitro measurements of flexibility using selective 2' hydroxyl acylation analyzed by primer extension and mutational profiling (SHAPE-MaP). Our simulations reproduce known relative binding affinity profiles for the cognate and synthetic ligands, and pinpoint how both ligands make use of different aspects of riboswitch flexibility. On the basis of our dissociation trajectories, we also make and validate predictions of pairs of mutations for both the ligand systems that would show differing binding affinities. These mutations are distal to the binding site and could not have been predicted solely on the basis of structure. The methodology demonstrated here shows how molecular dynamics simulations with all-atom force-fields have now come of age in making predictions that complement existing experimental techniques and illuminate aspects of systems otherwise not trivial to understand.

Exploration vs Convergence Speed in Adaptive-Bias Enhanced Sampling

In adaptive-bias enhanced sampling methods, a bias potential is added to the system to drive transitions between metastable states. The bias potential is a function of a few collective variables and is gradually modified according to the underlying free energy surface. We show that when the collective variables are suboptimal, there is an exploration–convergence tradeoff, and one must choose between a quickly converging bias that will lead to fewer transitions or a slower to converge bias that can explore the phase space more efficiently but might require a much longer time to produce an accurate free energy estimate. The recently proposed on-the-fly probability enhanced sampling (OPES) method focuses on fast convergence, but there are cases where fast exploration is preferred instead. For this reason, we introduce a new variant of the OPES method that focuses on quickly escaping metastable states at the expense of convergence speed. We illustrate the benefits of this approach in prototypical systems and show that it outperforms the popular metadynamics method.

Multiple Poses and Thermodynamics of Ligands Targeting Protein Surfaces: The Case of Furosemide Binding to mitoNEET in Aqueous Solution

Human NEET proteins, such as NAF-1 and mitoNEET, are homodimeric, redox iron-sulfur proteins characterized by triple cysteine and one histidine-coordinated [2Fe-2S] cluster. They exist in an oxidized and reduced state. Abnormal release of the cluster is implicated in a variety of diseases, including cancer and neurodegeneration. The computer-aided and structure-based design of ligands affecting cluster release is of paramount importance from a pharmaceutical perspective. Unfortunately, experimental structural information so far is limited to only one ligand/protein complex. This is the X-ray structure of furosemide bound to oxidized mitoNEET. Here we employ an enhanced sampling approach, Localized Volume-based Metadynamics, developed by some of us, to identify binding poses of furosemide to human mitoNEET protein in solution. The binding modes show a high variability within the same shallow binding pocket on the protein surface identified in the X-ray structure. Among the different binding conformations, one of them is in agreement with the crystal structure’s one. This conformation might have been overstabilized in the latter because of the presence of crystal packing interactions, absent in solution. The calculated binding affinity is compatible with experimental data. Our protocol can be used in a straightforward manner in drug design campaigns targeting this pharmaceutically important family of proteins.

PathDetect-SOM: A Neural Network Approach for the Identification of Pathways in Ligand Binding Simulations

Understanding the process of ligand–protein recognition is important to unveil biological mechanisms and to guide drug discovery and design. Enhanced-sampling molecular dynamics is now routinely used to simulate the ligand binding process, resulting in the need for suitable tools for the analysis of large data sets of binding events. Here, we designed, implemented, and tested PathDetect-SOM, a tool based on self-organizing maps to build concise visual models of the ligand binding pathways sampled along single simulations or replicas. The tool performs a geometric clustering of the trajectories and traces the pathways over an easily interpretable 2D map and, using an approximate transition matrix, it can build a graph model of concurrent pathways. The tool was tested on three study cases representing different types of problems and simulation techniques. A clear reconstruction of the sampled pathways was derived in all cases, and useful information on the energetic features of the processes was recovered. The tool is available at https://github.com/MottaStefano/PathDetect-SOM.

Ligand Binding Path Sampling Based on Parallel Cascade Selection Molecular Dynamics: LB-PaCS-MD

Parallel cascade selection molecular dynamics (PaCS-MD) is a rare-event sampling method that generates transition pathways between a reactant and product. To sample the transition pathways, PaCS-MD repeats short-time MD simulations from important configurations as conformational resampling cycles. In this study, PaCS-MD was extended to sample ligand binding pathways toward a target protein, which is referred to as LB-PaCS-MD. In a ligand-concentrated environment, where multiple ligand copies are randomly arranged around the target protein, LB-PaCS-MD allows for the frequent sampling of ligand binding pathways. To select the important configurations, we specified the center of mass (COM) distance between each ligand and the relevant binding site of the target protein, where snapshots generated by the short-time MD simulations were ranked by their COM distance values. From each cycle, snapshots with smaller COM distance values were selected as the important configurations to be resampled using the short-time MD simulations. By repeating conformational resampling cycles, the COM distance values gradually decreased and converged to constants, meaning that a set of ligand binding pathways toward the target protein was sampled by LB-PaCS-MD. To demonstrate relative efficiency, LB-PaCS-MD was applied to several proteins, and their ligand binding pathways were sampled more frequently than conventional MD simulations.

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.

A General Picture of Cucurbit[8]uril Host-Guest Binding

Describing, understanding, and designing complex interaction networks within macromolecular systems remain challenging in modern chemical research. Host-guest systems, despite their relative simplicity in both the structural feature and interaction patterns, still pose problems in theoretical modeling. The barrel-shaped supramolecular container cucurbit[8]uril (CB8) shows promising functionalities in various areas, e.g., catalysis and molecular recognition. It can stably coordinate a series of structurally diverse guests with high affinities. In this work, we examine the binding of seven commonly abused drugs to the CB8 host, aiming at providing a general picture of CB8-guest binding. Extensive sampling of the configurational space of these host-guest systems is performed, and the binding pathway and interaction patterns of CB8-guest complexes are investigated. A thorough comparison of widely used fixed-charge models for drug-like molecules is presented. Iterative refitting of the atomic charges suggests significant conformation dependence of charge generation. The initial model generated at the original conformation could be inaccurate for new conformations explored during conformational search, and the newly fitted charge set improves the prediction-experiment correlation significantly. Our investigations of the configurational space of CB8-drug complexes suggest that the host-guest interactions are more complex than expected. Despite the structural simplicities of these molecules, the conformational fluctuations of the host and the guest molecules and orientations of functional groups lead to the existence of an ensemble of binding modes. The insights of the binding thermodynamics, performance of fixed-charge models, and binding patterns of the CB8-guest systems are useful for studying and elucidating the binding mechanism of other host-guest complexes.

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.

The prediction of protein-ligand unbinding for modern drug discovery

INTRODUCTION

Drug-target thermodynamic and kinetic information have perennially important roles in drug design. The prediction of protein-ligand unbinding, which can provide important kinetic information, in experiments continues to face great challenges. Uncovering protein-ligand unbinding through molecular dynamics simulations has become efficient and inexpensive with the progress and enhancement of computing power and sampling methods.

AREAS COVERED

In this review, various sampling methods for protein-ligand unbinding and their basic principles are firstly briefly introduced. Then, their applications in predicting aspects of protein-ligand unbinding, including unbinding pathways, dissociation rate constants, residence time and binding affinity, are discussed.

EXPERT OPINION

Although various sampling methods have been successfully applied in numerous systems, they still have shortcomings and deficiencies. Most enhanced sampling methods require researchers to possess a wealth of prior knowledge of collective variables or reaction coordinates. In addition, most systems studied at present are relatively simple, and the study of complex systems in real drug research remains greatly challenging. Through the combination of machine learning and enhanced sampling methods, prediction accuracy can be further improved, and some problems encountered in complex systems also may be solved.

All-Atom Simulations Uncover the Molecular Terms of the NKCC1 Transport Mechanism

The secondary-active Na-K-Cl cotransporter 1 (NKCC1), member of the cation-chloride cotransporter (CCC) family, ensures the electroneutral movement of Cl-, Na+, and K+ ions across cellular membranes. NKCC1 regulates Cl- homeostasis and cell volume, handling a pivotal role in transepithelial water transport and neuronal excitability. Aberrant NKCC1 transport is hence implicated in a variety of human diseases (hypertension, renal disorders, neuropathies, and cancer). Building on the newly resolved NKCC1 cryo-EM structure, all-atom enhanced sampling simulations unprecedentedly unlock the mechanism of NKCC1-mediated ion transport, assessing the order and the molecular basis of its interdependent ion translocation. Our outcomes strikingly advance the understanding of the physiological mechanism of CCCs and disclose a key role of CCC-conserved asparagine residues, whose side-chain promiscuity ensures the transport of both negatively and positively charged ions along the same translocation route. This study sets a conceptual basis to devise NKCC-selective inhibitors to treat diseases linked to Cl- dishomeostasis.

Metadynamics-Based Approaches for Modeling the Hypoxia-Inducible Factor 2α Ligand Binding Process

Several methods based on enhanced-sampling molecular dynamics have been proposed for studying ligand binding processes. Here, we developed a protocol that combines the advantages of steered molecular dynamics (SMD) and metadynamics. While SMD is proposed for investigating possible unbinding pathways of the ligand and identifying the preferred one, metadynamics, with the path collective variable (PCV) formalism, is suggested to explore the binding processes along the pathway defined on the basis of SMD, by using only two CVs. We applied our approach to the study of binding of two known ligands to the hypoxia-inducible factor 2α, where the buried binding cavity makes simulation of the process a challenging task. Our approach allowed identification of the preferred entrance pathway for each ligand, highlighted the features of the bound and intermediate states in the free-energy surface, and provided a binding affinity scale in agreement with experimental data. Therefore, it seems to be a suitable tool for elucidating ligand binding processes of similar complex systems.

Modifying the catalytic preference of alpha-amylase toward n-alkanes for bioremediation purposes using in silico strategies

Since the beginning of oil exploration, whole ecosystems have been affected by accidents and bad practices involving petroleum compounds. In this sense, bioremediation stands out as the cheapest and most eco-friendly alternatives to reverse the damage done in oil-impacted areas. However, more efforts must be made to engineer enzymes that could be used in the bioremediation process. Interestingly, a recent work described that α-amylase, one of the most evolutionary conserved enzymes, was able to promiscuously degrade n-alkanes, a class of molecules abundant in the petroleum admixture. Considering that α-amylase is expressed in almost all known organisms, and employed in numerous biotechnological processes, using it can be a great leap toward more efficient applications of enzyme or microorganism-consortia bioremediation approaches. In this work, we employed a strict computational approach to design new α-amylase mutants with potentially enhanced catalytic efficiency toward n-alkanes. Using in silico techniques, such as molecular docking, molecular dynamics, metadynamics, and residue-residue interaction networks, we generated mutants potentially more efficient for degrading n-alkanes, L183Y, and N314A. Our results indicate that the new mutants have an increased binding rate for tetradecane, the longest n-alkane previously tested, which can reside in the catalytic center for more extended periods. Additionally, molecular dynamics and network analysis showed that the new mutations have no negative impact on protein structure than the WT. Our results aid in solidifying this enzyme as one more tool in the petroleum bioremediation toolbox.

Ligand unbinding mechanisms and kinetics for T4 lysozyme mutants from τRAMD simulations

The protein-ligand residence time, τ, influences molecular function in biological networks and has been recognized as an important determinant of drug efficacy. To predict τ, computational methods must overcome the problem that τ often exceeds the timescales accessible to conventional molecular dynamics (MD) simulation. Here, we apply the τ-Random Acceleration Molecular Dynamics (τRAMD) method to a set of kinetically characterized complexes of T4 lysozyme mutants with small, engineered binding cavities. τRAMD yields relative ligand dissociation rates in good accordance with experiments across this diverse set of complexes that differ with regard to measurement temperature, ligand identity, protein mutation and binding cavity. τRAMD thereby allows a comprehensive characterization of the ligand egress routes and determinants of τ. Although ligand dissociation by multiple egress routes is observed, we find that egress via the predominant route determines the value of τ. We also find that the presence of a greater number of metastable states along egress pathways leads to slower protein-ligand dissociation. These physical insights could be exploited in the rational optimization of the kinetic properties of drug candidates.

•

Relative residence times are computed for T4 lysozyme mutant-ligand complexes.

•

τ-Random Acceleration Molecular Dynamics provide efficient sampling of unbinding.

•

Computed dissociation rates show good agreement with all available measured values.

•

Ligand egress via the predominant route determines the value of the residence time.

•

The presence of metastable states along egress pathways slows down dissociation.

Pepsin-like aspartic proteases (PAPs) as model systems for combining biomolecular simulation with biophysical experiments

Pepsin-like aspartic proteases (PAPs) are a class of aspartic proteases which shares tremendous structural similarity with human pepsin. One of the key structural features of PAPs is the presence of a β-hairpin motif otherwise known as flap. The biological function of the PAPs is highly dependent on the conformational dynamics of the flap region. In apo PAPs, the conformational dynamics of the flap is dominated by the rotational degrees of freedom associated with χ1 and χ2 angles of conserved Tyr (or Phe in some cases). However it is plausible that dihedral order parameters associated with several other residues might play crucial roles in the conformational dynamics of apo PAPs. Due to their size, complexities associated with conformational dynamics and clinical significance (drug targets for malaria, Alzheimer's disease etc.), PAPs provide a challenging testing ground for computational and experimental methods focusing on understanding conformational dynamics and molecular recognition in biomolecules. The opening of the flap region is necessary to accommodate substrate/ligand in the active site of the PAPs. The BIG challenge is to gain atomistic details into how reversible ligand binding/unbinding (molecular recognition) affects the conformational dynamics. Recent reports of kinetics (K i, K d) and thermodynamic parameters (ΔH, TΔS, and ΔG) associated with macro-cyclic ligands bound to BACE1 (belongs to PAP family) provide a perfect challenge (how to deal with big ligands with multiple torsional angles and select optimum order parameters to study reversible ligand binding/unbinding) for computational methods to predict binding free energies and kinetics beyond typical test systems e.g. benzamide-trypsin. In this work, i reviewed several order parameters which were proposed to capture the conformational dynamics and molecular recognition in PAPs. I further highlighted how machine learning methods can be used as order parameters in the context of PAPs. I then proposed some open ideas and challenges in the context of molecular simulation and put forward my case on how biophysical experiments e.g. NMR, time-resolved FRET etc. can be used in conjunction with biomolecular simulation to gain complete atomistic insights into the conformational dynamics of PAPs.

Identification of viable TCDD access pathways to human AhR PAS-B ligand binding domain

Unintentionally released in the environment as by-products of industrial activities, dioxins, exemplified by 2,3,7,8-tetrachlorodibenzo-p-dioxin (TCDD), represent a primary concern for human health. Exposure to these chemicals is known to produce a broad spectrum of adverse effects, including cancer. The main mechanism of action of TCDD in humans involves binding to the Aryl hydrocarbon Receptor (AhR). Although qualitatively established, TCDD capture by the AhR remains poorly characterized at the molecular level. Starting from a recently developed structural model of the human AhR PAS-B domain, in this work we attempt the identification of viable TCDD access pathways to the human AhR ligand binding domain by means of molecular dynamics. Based on the result of metadynamics simulations, we identify two main regions that may potentially serve as access paths for TCDD. For each path, we characterize the residues closely interacting with TCDD, thereby suggesting a possible mechanism for TCDD capture. Our results are reviewed and discussed in the light of the available information about Human AhR structure and functions.

Ligand Affinities within the Open-Boundary Molecular Mechanics/Coarse-Grained Framework (I): Alchemical Transformations within the Hamiltonian Adaptive Resolution Scheme

Our recently developed Open-Boundary Molecular Mechanics/Coarse Grained (OB-MM/CG) framework predicts ligand poses in important pharmaceutical targets, such as G-protein Coupled Receptors, even when experimental structural information is lacking. The approach, which is based on GROMOS and AMBER force fields, allows for grand-canonical simulations of protein-ligand complexes by using the Hamiltonian Adaptive Resolution Scheme (H-AdResS) for the solvent. Here, we present a key step toward the estimation of ligand binding affinities for their targets within this approach. This is the implementation of the H-AdResS in the GROMACS code. The accuracy of our implementation is established by calculating hydration free energies of several molecules in water by means of alchemical transformations. The deviations of the GROMOS- and AMBER-based H-AdResS results from the reference fully atomistic simulations are smaller than the accuracy of the force field and/or they are in the range of the published results. Importantly, our predictions are in good agreement with experimental data. The current implementation paves the way to the use of the OB-MM/CG framework for the study of large biological systems.

Thermodynamics and Kinetics of Drug-Target Binding by Molecular Simulation

Computational studies play an increasingly important role in chemistry and biophysics, mainly thanks to improvements in hardware and algorithms. In drug discovery and development, computational studies can reduce the costs and risks of bringing a new medicine to market. Computational simulations are mainly used to optimize promising new compounds by estimating their binding affinity to proteins. This is challenging due to the complexity of the simulated system. To assess the present and future value of simulation for drug discovery, we review key applications of advanced methods for sampling complex free-energy landscapes at near nonergodicity conditions and for estimating the rate coefficients of very slow processes of pharmacological interest. We outline the statistical mechanics and computational background behind this research, including methods such as steered molecular dynamics and metadynamics. We review recent applications to pharmacology and drug discovery and discuss possible guidelines for the practitioner. Recent trends in machine learning are also briefly discussed. Thanks to the rapid development of methods for characterizing and quantifying rare events, simulation's role in drug discovery is likely to expand, making it a valuable complement to experimental and clinical approaches.

SAMPL7 TrimerTrip host-guest binding affinities from extensive alchemical and end-point free energy calculations

The prediction of host-guest binding affinities with computational modelling is still a challenging task. In the 7th statistical assessment of the modeling of proteins and ligands (SAMPL) challenge, a new host named TrimerTrip was synthesized and the thermodynamic parameters of 16 structurally diverse guests binding to the host were characterized. In the TrimerTrip-guest challenge, only structures of the host and the guests are provided, which indicates that the predictions of both the binding poses and the binding affinities are under assessment. In this work, starting from the binding poses obtained from our previous enhanced sampling simulations in the configurational space, we perform extensive alchemical and end-point free energy calculations to calculate the host-guest binding affinities retrospectively. The alchemical predictions with two widely accepted charge schemes (i.e. AM1-BCC and RESP) are in good agreement with the experimental reference, while the end-point estimates perform poorly in reproducing the experimental binding affinities. Aside from the absolute value of the binding affinity, the rank of binding free energies is also crucial in drug design. Surprisingly, the end-point MM/PBSA method seems very powerful in reproducing the experimental rank of binding affinities. Although the length of our simulations is long and the intermediate spacing is dense, the convergence behavior is not very good, which may arise from the flexibility of the host molecule. Enhanced sampling techniques in the configurational space may be required to obtain fully converged sampling. Further, as the length of sampling in alchemical free energy calculations already achieves several hundred ns, performing direct simulations of the binding/unbinding event in the physical space could be more useful and insightful. More details about the binding pathway and mechanism could be obtained in this way. The nonequilibrium method could also be a nice choice if one insists to use the alchemical method, as the intermediate sampling is avoided to some extent.

Enhanced Jarzynski free energy calculations using weighted ensemble

The free energy of transitions between stable states is the key thermodynamic quantity that governs the relative probabilities of the forward and reverse reactions and the ratio of state probabilities at equilibrium. The binding free energy of a drug and its receptor is of particular interest, as it serves as an optimization function for drug design. Over the years, many computational methods have been developed to calculate binding free energies, and while many of these methods have a long history, issues such as convergence of free energy estimates and the projection of a binding process onto order parameters remain. Over 20 years ago, the Jarzynski equality was derived with the promise to calculate equilibrium free energies by measuring the work applied to short nonequilibrium trajectories. However, these calculations were found to be dominated by trajectories with low applied work that occur with extremely low probability. Here, we examine the combination of weighted ensemble algorithms with the Jarzynski equality. In this combined method, an ensemble of nonequilibrium trajectories are run in parallel, and cloning and merging operations are used to preferentially sample low-work trajectories that dominate the free energy calculations. Two additional methods are also examined: (i) a novel weighted ensemble resampler that samples trajectories directly according to their importance to the work of work and (ii) the diffusion Monte Carlo method using the applied work as the selection potential. We thoroughly examine both the accuracy and efficiency of unbinding free energy calculations for a series of model Lennard-Jones atom pairs with interaction strengths ranging from 2 kcal/mol to 20 kcal/mol. We find that weighted ensemble calculations can more efficiently determine accurate binding free energies, especially for deeper Lennard-Jones well depths.

SAMPL7 TrimerTrip host-guest binding poses and binding affinities from spherical-coordinates-biased simulations

Host-guest binding remains a major challenge in modern computational modelling. The newest 7th statistical assessment of the modeling of proteins and ligands (SAMPL) challenge contains a new series of host-guest systems. The TrimerTrip host binds to 16 structurally diverse guests. Previously, we have successfully employed the spherical coordinates as the collective variables coupled with the enhanced sampling technique metadynamics to enhance the sampling of the binding/unbinding event, search for possible binding poses and calculate the binding affinities in all three host-guest binding cases of the 6th SAMPL challenge. In this work, we report a retrospective study on the TrimerTrip host-guest systems by employing the same protocol to investigate the TrimerTrip host in the SAMPL7 challenge. As no binding pose is provided by the SAMPL7 host, our simulations initiate from randomly selected configurations and are proceeded long enough to obtain converged free energy estimates and search for possible binding poses. The calculated binding affinities are in good agreement with the experimental reference, and the obtained binding poses serve as a nice starting point for end-point or alchemical free energy calculations. Note that as the work is performed after the close of the SAMPL7 challenge, we do not participate in the challenge and the results are not formally submitted to the SAMPL7 challenge.

Recent progress in molecular simulation methods for drug binding kinetics

Due to the contribution of drug-target binding kinetics to drug efficacy, there is a high level of interest in developing methods to predict drug-target binding kinetic parameters. During the review period, a wide range of enhanced sampling molecular dynamics simulation-based methods has been developed for computing drug-target binding kinetics and studying binding and unbinding mechanisms. Here, we assess the performance of these methods considering two benchmark systems in detail: mutant T4 lysozyme-ligand complexes and a large set of N-HSP90-inhibitor complexes. The results indicate that some of the simulation methods can already be usefully applied in drug discovery or lead optimization programs but that further studies on more high-quality experimental benchmark datasets are necessary to improve and validate computational methods.

Protein-ligand binding with the coarse-grained Martini model

The detailed understanding of the binding of small molecules to proteins is the key for the development of novel drugs or to increase the acceptance of substrates by enzymes. Nowadays, computer-aided design of protein–ligand binding is an important tool to accomplish this task. Current approaches typically rely on high-throughput docking essays or computationally expensive atomistic molecular dynamics simulations. Here, we present an approach to use the recently re-parametrized coarse-grained Martini model to perform unbiased millisecond sampling of protein–ligand interactions of small drug-like molecules. Remarkably, we achieve high accuracy without the need of any a priori knowledge of binding pockets or pathways. Our approach is applied to a range of systems from the well-characterized T4 lysozyme over members of the GPCR family and nuclear receptors to a variety of enzymes. The presented results open the way to high-throughput screening of ligand libraries or protein mutations using the coarse-grained Martini model.

Computer-aided design of protein-ligand binding is important for the development of novel drugs. Here authors present an approach to use the recently re-parametrized coarse-grained Martini model to perform unbiased millisecond sampling of protein-ligand binding interactions of small drug-like molecules.

Capturing Protein-Ligand Recognition Pathways in Coarse-Grained Simulation

Protein-ligand recognition is dynamic and complex. A key approach in deciphering the mechanism underlying the recognition process is to capture the kinetic process of the ligand in its act of binding to its designated protein cavity. Toward this end, ultralong all-atom molecular dynamics simulation has recently emerged as a popular method of choice because of its ability to record these events at high spatial and temporal resolution. However, success via this route comes at an exorbitant computational cost. Herein, we demonstrate that coarse-grained models of the protein, when systematically optimized to maintain its tertiary fold, can capture the complete process of spontaneous protein-ligand binding from bulk media to the cavity at crystallographic precision and within wall clock time that is orders of magnitude shorter than that of all-atom simulations. The exhaustive sampling of ligand exploration in protein and solvent, harnessed by coarse-grained simulation, leads to elucidation of new ligand recognition pathways and discovery of non-native binding poses.