The free energy landscape of small molecule unbinding

Assessment of Kaistella jeonii esterase conformational dynamics in response to poly(ethylene terephthalate) binding

The pervasive presence of plastic in the environment has reached a concerning scale, being identified in many ecosystems. Bioremediation is the cheapest and most eco-friendly alternative to remove this polymer from affected areas. Recent work described that a novel cold-active esterase enzyme extracted from the bacteria Kaistella jeonii could promiscuously degrade PET. Compared to the well-known PETase from Ideonella sakaiensis, this novel esterase presents a low sequence identity yet has a remarkably similar folding. However, enzymatic assays demonstrated a lower catalytic efficiency. In this work, we employed a strict computational approach to investigate the binding mechanism between the esterase and PET. Understanding the underlying mechanism of binding can shed light on the evolutive mechanism of how enzymes have been evolving to degrade these artificial molecules and help develop rational engineering approaches to improve PETase-like enzymes. Our results indicate that this esterase misses a disulfide bridge, keeping the catalytic residues closer and possibly influencing its catalytic efficiency. Moreover, we describe the structural response to the interaction between enzyme and PET, indicating local and global effects. Our results aid in deepening the knowledge behind the mechanism of biological catalysis of PET degradation and as a base for the engineering of novel PETases.

Identification and mechanistic exploration of structural and conformational dynamics of NF-kB inhibitors: rationale insights from in silico and in vitro studies

Increased expression of target genes that code for proinflammatory chemical mediators results from a series of intracellular cascades triggered by activation of dysregulated NF-κB signaling pathway. Dysfunctional NF-kB signaling amplifies and perpetuates autoimmune responses in inflammatory diseases, including psoriasis. This study aimed to identify therapeutically relevant NF-kB inhibitors and elucidate the mechanistic aspects behind NF-kB inhibition. After virtual screening and molecular docking, five hit NF-kB inhibitors opted, and their therapeutic efficacy was examined using cell-based assays in TNF-α stimulated human keratinocyte cells. To investigate the conformational changes of target protein and inhibitor-protein interaction mechanisms, molecular dynamics (MD) simulations, binding free energy calculations together with principal component (PC) analysis, dynamics cross-correlation matrix analysis (DCCM), free energy landscape (FEL) analysis and quantum mechanical calculations were carried out. Among identified NF-kB inhibitors, myricetin and hesperidin significantly scavenged intracellular ROS and inhibited NF-kB activation. Analysis of the MD simulation trajectories of ligand-protein complexes revealed that myricetin and hesperidin formed energetically stabilized complexes with the target protein and were able to lock NF-kB in a closed conformation. Myricetin and hesperidin binding to the target protein significantly impacted conformational changes and internal dynamics of amino acid residues in protein domains. Tyr57, Glu60, Lys144 and Asp239 residues majorly contributed to locking the NF-kB in a closed conformation. The combinatorial approach employing in silico tools integrated with cell-based approaches substantiated the binding mechanism and NF-kB active site inhibition by the lead molecule myricetin, which can be explored as a viable antipsoriatic drug candidate associated with dysregulated NF-kB.Communicated by Ramaswamy H. Sarma.

Reducing the assemblies of amyloid-beta multimers by sodium dodecyl sulfate surfactant at concentrations lower than critical micelle concentration: molecular dynamics simulation exploration

Amyloid-β peptide, the predominant proteinaceous component of senile plaques, is responsible for the incidence of Alzheimer's disease (AD), an age-associated neurodegenerative disorder. Specifically, the amyloid-β(1-42) (Aβ1-42) isoform, known for its high toxicity, is the predominant biomarker for the preliminary diagnosis of AD. The aggregation of the Aβ1-42 peptides can be affected by the components of the cellular medium through changing their structures and molecular interactions. In this study, we investigated the effect of sodium dodecyl sulfate (SDS) at much lower concentrations than the critical micelle concentration (CMC) on Aβ1-42 aggregation. For this purpose, we studied mono-, di-, tri- and tetramers of Aβ1-42 peptide in two different concentrations of SDS molecules (10 and 40 molecules) using a 300 ns molecular dynamics simulation for each system. The distance between the center of mass (COM) of Aβ1-42 peptides confirms that an increase in the number of SDS molecules decreases their aggregation probability due to greater interaction with SDS molecules. Besides, the less compactness parameter reveals the reduced aggregation probability of Aβ1-42 peptides. Based on the energetic FEL landscapes, SDS molecules with the concentration closer to the CMC are an effective inhibitory agent to prevent the formation of Aβ1-42 fibrils. Also, the aggregation direction of the peptide pairs can be predicted by determining the direction of the accumulation-deterrent forces.Communicated by Ramaswamy H. Sarma.

15 Years of molecular simulation of drug-binding kinetics

INTRODUCTION

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

AREAS COVERED

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

EXPERT OPINION

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

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

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

Application of Computational Biology and Artificial Intelligence in Drug Design

Traditional drug design requires a great amount of research time and developmental expense. Booming computational approaches, including computational biology, computer-aided drug design, and artificial intelligence, have the potential to expedite the efficiency of drug discovery by minimizing the time and financial cost. In recent years, computational approaches are being widely used to improve the efficacy and effectiveness of drug discovery and pipeline, leading to the approval of plenty of new drugs for marketing. The present review emphasizes on the applications of these indispensable computational approaches in aiding target identification, lead discovery, and lead optimization. Some challenges of using these approaches for drug design are also discussed. Moreover, we propose a methodology for integrating various computational techniques into new drug discovery and design.

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.

SEEKR2: Versatile Multiscale Milestoning Utilizing the OpenMM Molecular Dynamics Engine

We present SEEKR2 (simulation-enabled estimation of kinetic rates version 2)—the latest iteration in the family of SEEKR programs for using multiscale simulation methods to computationally estimate the kinetics and thermodynamics of molecular processes, in particular, ligand-receptor binding. SEEKR2 generates equivalent, or improved, results compared to the earlier versions of SEEKR but with significant increases in speed and capabilities. SEEKR2 has also been built with greater ease of usability and with extensible features to enable future expansions of the method. Now, in addition to supporting simulations using NAMD, calculations may be run with the fast and extensible OpenMM simulation engine. The Brownian dynamics portion of the calculation has also been upgraded to Browndye 2. Furthermore, this version of SEEKR supports hydrogen mass repartitioning, which significantly reduces computational cost, while showing little, if any, loss of accuracy in the predicted kinetics.

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

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

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

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

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.

Antipyretic activity of Caesalpinia digyna (Rottl.) leaves extract along with phytoconstituent's binding affinity to COX-1, COX-2, and mPGES-1 receptors: In vivo and in silico approaches

Caesalpinia digyna (Rottl.) (Family: Fabaceae) is well known for its numerous medicinal values against several human disorders including fever, senile pruritis, diarrhea, tuberculosis, tonic disorder, diabetes, etc. The current study is intended to investigate the in vivo antipyretic activity of the methanol extract of C. digyna leaves (MECD) and its carbon-tetrachloride (CTCD) and butanol fraction (BTCD). Besides, in silico molecular docking and ADME/T profiling of the selective identified bioactive compounds of C. digyna has been also studied to validate the experimental outcomes and establish a better insight into the possible receptor-ligand interaction affinity. In vivo antipyretic activity of MECD, CTCD and BTCD were evaluated by employing yeast induced pyrexia technique in mice model and in silico analysis of the identified compounds of C. digyna has been implemented using PyRx autodock vina, Discovery Studio 2020, UCSF Chimera software and ADME/T online tools. MECD and BTCD unveiled significant antipyretic activity in dose dependent manner whereas, CTCD failed to exhibit significant antipyretic activity. Comparing to other test sample, MECD (400 mg/kg; b.w) (p < 0.001) displayed maximum inhibition of pyrexia. In molecular docking approach, docking score between −6.60 to −10.20 kcal/mol have been revealed. Besides, in ADME/T analysis, no compound violated the lipiniski’s 5 rules and displayed any toxicity. Biological and computational approaches ascertain the ethno-botanical use of C. digyna as a good agent against pyrexia and the compounds of C. digyna are primarily proved as safe. Hereafter, further analysis is suggested to validate this research.

An Updated Review of Computer-Aided Drug Design and Its Application to COVID-19

The recent outbreak of the deadly coronavirus disease 19 (COVID-19) pandemic poses serious health concerns around the world. The lack of approved drugs or vaccines continues to be a challenge and further necessitates the discovery of new therapeutic molecules. Computer-aided drug design has helped to expedite the drug discovery and development process by minimizing the cost and time. In this review article, we highlight two important categories of computer-aided drug design (CADD), viz., the ligand-based as well as structured-based drug discovery. Various molecular modeling techniques involved in structure-based drug design are molecular docking and molecular dynamic simulation, whereas ligand-based drug design includes pharmacophore modeling, quantitative structure-activity relationship (QSARs), and artificial intelligence (AI). We have briefly discussed the significance of computer-aided drug design in the context of COVID-19 and how the researchers continue to rely on these computational techniques in the rapid identification of promising drug candidate molecules against various drug targets implicated in the pathogenesis of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2). The structural elucidation of pharmacological drug targets and the discovery of preclinical drug candidate molecules have accelerated both structure-based as well as ligand-based drug design. This review article will help the clinicians and researchers to exploit the immense potential of computer-aided drug design in designing and identification of drug molecules and thereby helping in the management of fatal disease.

Deciphering the Agonist Binding Mechanism to the Adenosine A1 Receptor

Despite being among the most characterized G protein-coupled receptors (GPCRs), adenosine receptors (ARs) have always been a difficult target in drug design. To date, no agonist other than the natural effector and the diagnostic regadenoson has been approved for human use. Recently, the structure of the adenosine A1 receptor (A1R) was determined in the active, Gi protein complexed state; this has important repercussions for structure-based drug design. Here, we employed supervised molecular dynamics simulations and mutagenesis experiments to extend the structural knowledge of the binding of selective agonists to A1R. Our results identify new residues involved in the association and dissociation pathway, they suggest the binding mode of N6-cyclopentyladenosine (CPA) related ligands, and they highlight the dramatic effect that chemical modifications can have on the overall binding mechanism, paving the way for the rational development of a structure-kinetics relationship of A1R agonists.

Antidepressant, anxiolytic, antipyretic, and thrombolytic profiling of methanol extract of the aerial part of Piper nigrum: In vivo, in vitro, and in silico approaches

Piper nigrum L. also called black pepper is popular for its numerous uses. The present research is designed to investigate the pharmacological potential of methanol extract of Piper nigrum (MEPN). The antidepressant investigation was performed by using both in vivo forced swimming test (FST) and tail suspension test (TST) methods while the anxiolytic research by hole-board test (HBT) method. Again, the antipyretic analysis was conducted through yeast-induced pyrexia method, whereas clot lysis activity was employed by the thrombolytic method. Furthermore, in silico studies followed by molecular docking analysis of several secondary metabolites, pass prediction, and ADME/T were evaluated with AutoDock Vina, Discovery Studio 2020, UCSF Chimera software PASS online, and ADME/T online tools. The plant extract demonstrated dose-dependent potentiality in antidepressant, anxiolytic, antipyretic, and thrombolytic activities. Induction of MEPN produced a significant (p < .5, p < .001) increase of mobility in FST and TST, and increased the head dipping and decreased the latency of time (p < .01, p < .001) in HBT. MEPN 400 (mg/kg; b.w.; p.o.) lowered the rectal temperature of yeast-induced pyrexia substantially (p < .001). Besides, MEPN produced promising (p < .001) clot lysis activity. In the computational approach, among all the proteins, a docking score was found ranging from -1.0 to -7.90 kcal/mol. Besides, all the compounds were found safe in ADME/T study. The results of our scientific research validate the suitability of this plant as an alternative source of novel therapeutics.

Protein storytelling through physics

Every protein has a story-how it folds, what it binds, its biological actions, and how it misbehaves in aging or disease. Stories are often inferred from a protein's shape (i.e., its structure). But increasingly, stories are told using computational molecular physics (CMP). CMP is rooted in the principled physics of driving forces and reveals granular detail of conformational populations in space and time. Recent advances are accessing longer time scales, larger actions, and blind testing, enabling more of biology's stories to be told in the language of atomistic physics.

Kinetics and free energy of ligand dissociation using weighted ensemble milestoning

We consider the recently developed weighted ensemble milestoning (WEM) scheme [D. Ray and I. Andricioaei, J. Chem. Phys. 152, 234114 (2020)] and test its capability of simulating ligand-receptor dissociation dynamics. We performed WEM simulations on the following host-guest systems: Na⁺/Cl^- ion pair and 4-hydroxy-2-butanone ligand with FK506 binding protein. As a proof of principle, we show that the WEM formalism reproduces the Na⁺/Cl^- ion pair dissociation timescale and the free energy profile obtained from long conventional MD simulation. To increase the accuracy of WEM calculations applied to kinetics and thermodynamics in protein-ligand binding, we introduced a modified WEM scheme called weighted ensemble milestoning with restraint release (WEM-RR), which can increase the number of starting points per milestone without adding additional computational cost. WEM-RR calculations obtained a ligand residence time and binding free energy in agreement with experimental and previous computational results. Moreover, using the milestoning framework, the binding time and rate constants, dissociation constants, and committor probabilities could also be calculated at a low computational cost. We also present an analytical approach for estimating the association rate constant (k_on) when binding is primarily diffusion driven. We show that the WEM method can efficiently calculate multiple experimental observables describing ligand-receptor binding/unbinding and is a promising candidate for computer-aided inhibitor design.

Rapid prediction of crucial hotspot interactions for icosahedral viral capsid self-assembly by energy landscape atlasing validated by mutagenesis

Icosahedral viruses are under a micrometer in diameter, their infectious genome encapsulated by a shell assembled by a multiscale process, starting from an integer multiple of 60 viral capsid or coat protein (VP) monomers. We predict and validate inter-atomic hotspot interactions between VP monomers that are important for the assembly of 3 types of icosahedral viral capsids: Adeno Associated Virus serotype 2 (AAV2) and Minute Virus of Mice (MVM), both T = 1 single stranded DNA viruses, and Bromo Mosaic Virus (BMV), a T = 3 single stranded RNA virus. Experimental validation is by in-vitro, site-directed mutagenesis data found in literature. We combine ab-initio predictions at two scales: at the interface-scale, we predict the importance (cruciality) of an interaction for successful subassembly across each interface between symmetry-related VP monomers; and at the capsid-scale, we predict the cruciality of an interface for successful capsid assembly. At the interface-scale, we measure cruciality by changes in the capsid free-energy landscape partition function when an interaction is removed. The partition function computation uses atlases of interface subassembly landscapes, rapidly generated by a novel geometric method and curated opensource software EASAL (efficient atlasing and search of assembly landscapes). At the capsid-scale, cruciality of an interface for successful assembly of the capsid is based on combinatorial entropy. Our study goes all the way from resource-light, multiscale computational predictions of crucial hotspot inter-atomic interactions to validation using data on site-directed mutagenesis' effect on capsid assembly. By reliably and rapidly narrowing down target interactions, (no more than 1.5 hours per interface on a laptop with Intel Core i5-2500K @ 3.2 Ghz CPU and 8GB of RAM) our predictions can inform and reduce time-consuming in-vitro and in-vivo experiments, or more computationally intensive in-silico analyses.

In silico identification of novel PrfA inhibitors to fight listeriosis: A virtual screening and molecular dynamics studies

Listeria monocytogenes is considered to be one of the most dangerous foodborne pathogens as it can cause listeriosis, a life-threatening human disease. While the incidence of listeriosis is very low its fatality rate is exceptionally high. Because many multi-resistance Listeria monocytogenes strains that do not respond to conventional antibiotic therapy have been recently described, development of new antimicrobials to fight listeriosis is necessary. The positive regulatory factor A (PrfA) is a key homodimeric transcription factor that modulates the transcription of multiple virulence factors which are ultimately responsible of Listeria monocytogenes' pathogenicity. In the present manuscript we describe several new potential PrfA inhibitors that were identified after performing ligand-based virtual screening followed by molecular docking calculations against the wild-type PrfA structure. The three top-scored drug-likeness inhibitors bound to the wild-type PrfA structure were further assessed by Molecular Dynamics (MD) simulations. Besides, the three top-scored inhibitors were docked into a constitutive active apoPrfA mutant structure and the corresponding complexes were also simulated by MD. According to the obtained data, PUBChem 87534955 (P875) and PUBChem 58473762 (P584) may not only bind and inhibit wild-type PrfA but the aforementioned apoPrfA mutant as well. Therefore, P875 and P584 might represent good starting points for the development of a completely new set of antimicrobial agents to treat listeriosis.

The investigation of the G-quadruplex aptamer selectivity to Pb2+ ion: a joint molecular dynamics simulation and density functional theory study

The aptamers with the ability to form a G-quadruplex structure can be stable in the presence of some ions. Hence, study of the interactions between such aptamers and ions can be beneficial to determine the highest selective aptamer toward an ion. In this article, molecular dynamics (MD) simulations and quantum mechanics (QM) calculations have been applied to investigate the selectivity of the T30695 aptamer toward Pb²⁺ in comparison with some ions. The Free Energy Landscape (FEL) analysis indicates that Pb²⁺ has remained inside the aptamer during the MD simulation, while the other ions have left it. The Molecular Mechanics Poisson-Boltzmann Surface Area (MM-PBSA) binding energies prove that the conformational stability of the aptamer is the highest in the presence of Pb²⁺. According to the compaction parameters, the greatest compressed ion-aptamer complex, and hence, the highest ion-aptamer interaction have been induced in the presence of Pb²⁺. The contact maps clarify the closer contacts between the nucleotides of the aptamer in the presence of Pb²⁺. The density functional theory (DFT) results show that Pb²⁺ forms the most stable complex with the aptamer, which is consistent with the MD results. The QM calculations reveal that the N-H bonds and the O…H distances are the longest and the shortest, respectively, in the presence of Pb²⁺. The obtained results verify that the strongest hydrogen bonds (HBs), and hence, the most compressed aptamer structure are induced by Pb²⁺. Besides, atoms in molecules (AIM) and natural bond orbital (NBO) analyses confirm the results.Communicated by Ramaswamy H. Sarma.

Exploring the Ligand Binding/Unbinding Pathway by Selectively Enhanced Sampling of Ligand in a Protein-Ligand Complex

Understanding the protein-ligand binding is of fundamental biological interest and is essential for structure-based drug design. The difficulty in capturing the dynamic process, however, poses a great challenge for current experimental and theoretical simulation techniques. A selective integrated-tempering-sampling molecular dynamics (SITSMD) method offering an option for selectively enhanced sampling of the ligand in a protein-ligand complex was utilized to quantitatively illuminate the binding of benzamidine to the wild-type trypsin protease and its two mutants (S214E and S214K). The SITSMD simulations could produce consistent results as the extensive conventional MD simulation and gave additional insights into the binding pathway for the test protein-ligand complex system using significantly saved computational resource and time, indicating the potential of such a method in investigating protein-ligand binding. Additionally, the simulations identified the different roles of trypsin-benzamidine van der Waals (vdW) and electrostatic interactions in the binding: the former interaction works as the driving force for dragging the benzamidine close to the native binding pocket, and the latter interaction mainly contributes to stabilizing the benzamidine inside the pocket. The S214E mutation introduces more favorable electrostatic interactions, and as a result, both vdW and electrostatic interactions drive the benzamidine binding, lowering the binding and unbinding free energy barrier. In contrast, the S214K mutation prohibits the binding of the benzamidine to the native ligand binding pocket by introducing disliked charge-charge interactions. In summary, these findings suggest that the change in specific residues could modify the protein druggability, including the binding kinetics and thermodynamics.

Dynamic graphical models of molecular kinetics

Most current molecular dynamics simulation and analysis methods rely on the idea that the molecular system can be represented by a single global state (e.g., a Markov state in a Markov state model [MSM]). In this approach, molecules can be extensively sampled and analyzed when they only possess a few metastable states, such as small- to medium-sized proteins. However, this approach breaks down in frustrated systems and in large protein assemblies, where the number of global metastable states may grow exponentially with the system size. To address this problem, we here introduce dynamic graphical models (DGMs) that describe molecules as assemblies of coupled subsystems, akin to how spins interact in the Ising model. The change of each subsystem state is only governed by the states of itself and its neighbors. DGMs require fewer parameters than MSMs or other global state models; in particular, we do not need to observe all global system configurations to characterize them. Therefore, DGMs can predict previously unobserved molecular configurations. As a proof of concept, we demonstrate that DGMs can faithfully describe molecular thermodynamics and kinetics and predict previously unobserved metastable states for Ising models and protein simulations.

REVO: Resampling of ensembles by variation optimization

Conventional molecular dynamics simulations are incapable of sampling many important interactions in biomolecular systems due to their high dimensionality and rough energy landscapes. To observe rare events and calculate transition rates in these systems, enhanced sampling is a necessity. In particular, the study of ligand-protein interactions necessitates a diverse ensemble of protein conformations and transition states, and for many systems, this occurs on prohibitively long time scales. Previous strategies such as WExplore that can be used to determine these types of ensembles are hindered by problems related to the regioning of conformational space. Here, we propose a novel, regionless, enhanced sampling method that is based on the weighted ensemble framework. In this method, a value referred to as "trajectory variation" is optimized after each cycle through cloning and merging operations. This method allows for a more consistent measurement of observables and broader sampling resulting in the efficient exploration of previously unexplored conformations. We demonstrate the performance of this algorithm with the N-dimensional random walk and the unbinding of the trypsin-benzamidine system. The system is analyzed using conformation space networks, the residence time of benzamidine is confirmed, and a new unbinding pathway for the trypsin-benzamidine system is found. We expect that resampling of ensembles by variation optimization will be a useful general tool to broadly explore free energy landscapes.

Kinetics of Drug Binding and Residence Time

The kinetics of drug binding and unbinding is assuming an increasingly crucial role in the long, costly process of bringing a new medicine to patients. For example, the time a drug spends in contact with its biological target is known as residence time (the inverse of the kinetic constant of the drug-target unbinding, 1/ koff). Recent reports suggest that residence time could predict drug efficacy in vivo, perhaps even more effectively than conventional thermodynamic parameters (free energy, enthalpy, entropy). There are many experimental and computational methods for predicting drug-target residence time at an early stage of drug discovery programs. Here, we review and discuss the methodological approaches to estimating drug binding kinetics and residence time. We first introduce the theoretical background of drug binding kinetics from a physicochemical standpoint. We then analyze the recent literature in the field, starting from the experimental methodologies and applications thereof and moving to theoretical and computational approaches to the kinetics of drug binding and unbinding. We acknowledge the central role of molecular dynamics and related methods, which comprise a great number of the computational methods and applications reviewed here. However, we also consider kinetic Monte Carlo. We conclude with the outlook that drug (un)binding kinetics may soon become a go/no go step in the discovery and development of new medicines.

Escape of a Small Molecule from Inside T4 Lysozyme by Multiple Pathways

The T4 lysozyme L99A mutant is often used as a model system to study small-molecule binding to proteins, but pathways for ligand entry and exit from the buried binding site and the associated protein conformational changes have not been fully resolved. Here, molecular dynamics simulations were employed to model benzene exit from its binding cavity using the weighted ensemble (WE) approach to enhance sampling of low-probability unbinding trajectories. Independent WE simulations revealed four pathways for benzene exit, which correspond to transient tunnels spontaneously formed in previous simulations of apo T4 lysozyme. Thus, benzene unbinding occurs through multiple pathways partially created by intrinsic protein structural fluctuations. Motions of several α-helices and side chains were involved in ligand escape from metastable microstates. WE simulations also provided preliminary estimates of rate constants for each exit pathway. These results complement previous works and provide a semiquantitative characterization of pathway heterogeneity for binding of small molecules to proteins.

Kinetics of Ligand-Protein Dissociation from All-Atom Simulations: Are We There Yet?

Large parallel gains in the development of both computational resources and sampling methods have now made it possible to simulate dissociation events in ligand-protein complexes with all-atom resolution. Such encouraging progress, together with the inherent spatiotemporal resolution associated with molecular simulations, has left their use for investigating dissociation processes brimming with potential, both in rational drug design, where it can be an invaluable tool for determining the mechanistic driving forces behind dissociation rate constants, and in force-field development, where it can provide a catalog of transient molecular structures with which to refine force fields. Although much progress has been made in making force fields more accurate, reducing their error for transient structures along a transition path could yet prove to be a critical development helping to make kinetic predictions much more accurate. In what follows, we will provide a state-of-the-art compilation of the enhanced sampling methods based on molecular dynamics (MD) simulations used to investigate the kinetics and mechanisms of ligand-protein dissociation processes. Due to the time scales of such processes being slower than what is accessible using straightforward MD simulations, several ingenious schemes are being devised at a rapid rate to overcome this obstacle. Here we provide an up-to-date compendium of such methods and their achievements and shortcomings in extracting mechanistic insight into ligand-protein dissociation. We conclude with a critical and provocative appraisal attempting to answer the title of this Perspective.

Solvents to Fragments to Drugs: MD Applications in Drug Design

Simulations of molecular dynamics (MD) are playing an increasingly important role in structure-based drug discovery (SBDD). Here we review the use of MD for proteins in aqueous solvation, organic/aqueous mixed solvents (MDmix) and with small ligands, to the classic SBDD problems: Binding mode and binding free energy predictions. The simulation of proteins in their condensed state reveals solvent structures and preferential interaction sites (hot spots) on the protein surface. The information provided by water and its cosolvents can be used very effectively to understand protein ligand recognition and to improve the predictive capability of well-established methods such as molecular docking. The application of MD simulations to the study of the association of proteins with drug-like compounds is currently only possible for specific cases, as it remains computationally very expensive and labor intensive. MDmix simulations on the other hand, can be used systematically to address some of the common tasks in SBDD. With the advent of new tools and faster computers we expect to see an increase in the application of mixed solvent MD simulations to a plethora of protein targets to identify new drug candidates.

Water access and ligand dissociation at the binding site of proteins

Although water is undoubtedly an essential mediator of protein-ligand interactions, whether or not such water molecules are critical for the progress of ligand dissociation remains unclear. To gain a more complete understanding, molecular dynamics simulations are performed with two molecular systems, rigid model binding sites and trypsin-benzamidine. Free-energy landscapes are calculated with a suitably chosen solvent coordinate, which well describes water access to the ligand binding site. The results of free energy provided clear description of water-ligand exchange process, where two different mechanisms appear depending on whether the binding site is buried or not. As the site is more buried, water access is more difficult. When water does not access the site, ligand dissociation produces a large energy barrier, i.e., slow dissociation kinetics. This indicates that control of ligand dissociation kinetics becomes possible with burying the binding site. However, the results also showed that appropriate burying is important because burying reduces not only water access but also ligand binding. The role of the protein structural change is also discussed; it likely plays a similar role to water access because during ligand dissociation, it can make new coordination with the ligand binding site like water. These results contribute to the future pharmaceutical drug design and will be useful for fundamental exploration of various molecular events.

Mapping the Ligand Binding Landscape

The interaction between a ligand and a protein involves a multitude of conformational states. To achieve a particular deeply bound pose, the ligand must search across a rough free-energy landscape with many metastable minima. Creating maps of the ligand binding landscape is a great challenge, as binding and release events typically occur on timescales that are beyond the reach of molecular simulation. The WExplore enhanced sampling method is well suited to build these maps because it is designed to broadly explore free-energy landscapes and is capable of simulating ligand release pathways that occur on timescales as long as minutes. WExplore also uses only unbiased trajectory segments, allowing for the construction of Markov state models (MSMs) and conformation space networks that combine the results of multiple simulations. Here, we use WExplore to study two bromodomain-inhibitor systems using multiple docked starting poses (Brd4-MS436 and Baz2B-ICR7) and synthesize our results using a series of MSMs using time-lagged independent component analysis. Ranking the starting poses by exit rate agrees with the crystal structure pose in both cases. We also predict the most stable pose using the equilibrium populations from the MSM but find that the prediction is not robust as a function of MSM parameters. The simulated trajectories are synthesized into network models that visualize the entire binding landscape for each system, and we examine transition paths between deeply bound stable states. We find that, on average, transitions between deeply bound states convert through the unbound state 81% of the time, implying a trial-and-error approach to ligand binding. We conclude with a discussion of the implications of this result for both kinetics-based drug discovery and virtual screening pipelines that incorporate molecular dynamics.

Estimation of Drug-Target Residence Times by τ-Random Acceleration Molecular Dynamics Simulations

Drug-target residence time (τ), one of the main determinants of drug efficacy, remains highly challenging to predict computationally and, therefore, is usually not considered in the early stages of drug design. Here, we present an efficient computational method, τ-random acceleration molecular dynamics (τRAMD), for the ranking of drug candidates by their residence time and obtaining insights into ligand-target dissociation mechanisms. We assessed τRAMD on a data set of 70 diverse drug-like ligands of the N-terminal domain of HSP90α, a pharmaceutically important target with a highly flexible binding site, obtaining computed relative residence times with an accuracy of about 2.3τ for 78% of the compounds and less than 2.0τ within congeneric series. Analysis of dissociation trajectories reveals features that affect ligand unbinding rates, including transient polar interactions and steric hindrance. These results suggest that τRAMD will be widely applicable as a computationally efficient aid to improving drug residence times during lead optimization.

Protein structure-based drug design: from docking to molecular dynamics

Recent years have witnessed rapid developments of computer-aided drug design methods, which have reached accuracy that allows their routine practical applications in drug discovery campaigns. Protein structure-based methods are useful for the prediction of binding modes of small molecules and their relative affinity. The high-throughput docking of up to 10⁶ small molecules followed by scoring based on implicit-solvent force field can robustly identify micromolar binders using a rigid protein target. Molecular dynamics with explicit solvent is a low-throughput technique for the characterization of flexible binding sites and accurate evaluation of binding pathways, kinetics, and thermodynamics. In this review we highlight recent advancements in applications of ligand docking tools and molecular dynamics simulations to ligand identification and optimization.

Direction-dependent force-induced dissociation dynamics of an entropic-driven lock-and-key assembly

The unbinding dynamics of a nanosized sphere-and-cavity assembly under the pulling of constant force and constant loading rate is explored by dissipative particle dynamics simulations. The formation of this matched lock-and-key pair in a polymer solution is driven by the depletion attraction. The two-dimensional free energy landscape U(x,z) associated with this assembly is constructed. Our results indicate that the unbinding pathway along the orientation of the assembly is unfavorable due to the relatively high energy barrier compared to that along the tortuous minimum path whose energy barrier is not high. It is also found that the dissociation rate depends on the direction of the external force (θ) with respect to the assembly orientation. The presence of the force component perpendicular to the assembly orientation can reduce the bond lifetime significantly by driving the key particle to approach the minimum path. Moreover, the dissociation dynamics can be facilitated even by a pushing force compared to the spontaneous dissociation (without forces). To elucidate the effective pathway under pulling, the escaping position is analyzed and its mean direction with respect to the assembly orientation rises generally with increasing θ, revealing that the presence of the force component along the minimum pathway is helpful. The importance of the direction of the external pulling has been demonstrated in our simple system. Therefore, this effect should be considered in more complicated unbinding experiments.

Multiple Ligand Unbinding Pathways and Ligand-Induced Destabilization Revealed by WExplore

We report simulations of full ligand exit pathways for the trypsin-benzamidine system, generated using the sampling technique WExplore. WExplore is able to observe millisecond-scale unbinding events using many nanosecond-scale trajectories that are run without introducing biasing forces. The algorithm generates rare events by dividing the coordinate space into regions, on-the-fly, and balancing computational effort between regions through cloning and merging steps, as in the weighted ensemble method. The averaged exit flux yields a ligand exit rate of 180 μs, which is within an order of magnitude of the experimental value. We obtain broad sampling of ligand exit pathways, and visualize our findings using conformation space networks. The analysis shows three distinct exit channels, two of which are formed through large, rare motions of the loop regions in trypsin. This broad set of ligand-bound poses is then used to investigate general properties of ligand binding: we observe both a direct stabilizing effect of ligand-protein interactions and an indirect destabilizing effect on intraprotein interactions that is induced by the ligand. Significantly, the crystallographic binding poses are distinguished not only because their ligands induce large stabilizing effects, but also because they induce relatively low indirect destabilizations.

Quantitative Characterization of the Binding and Unbinding of Millimolar Drug Fragments with Molecular Dynamics Simulations

A quantitative characterization of the binding properties of drug fragments to a target protein is an important component of a fragment-based drug discovery program. Fragments typically have a weak binding affinity, however, making it challenging to experimentally characterize key binding properties, including binding sites, poses, and affinities. Direct simulation of the binding equilibrium by molecular dynamics (MD) simulations can provide a computational route to characterize fragment binding, but this approach is so computationally intensive that it has thus far remained relatively unexplored. Here, we perform MD simulations of sufficient length to observe several different fragments spontaneously and repeatedly bind to and unbind from the protein FKBP, allowing the binding affinities, on- and off-rates, and relative occupancies of alternative binding sites and alternative poses within each binding site to be estimated, thereby illustrating the potential of long time scale MD as a quantitative tool for fragment-based drug discovery. The data from the long time scale fragment binding simulations reported here also provide a useful benchmark for testing alternative computational methods aimed at characterizing fragment binding properties. As an example, we calculated binding affinities for the same fragments using a standard free energy perturbation approach and found that the values agreed with those obtained from the fragment binding simulations within statistical error.

How and when does an anticancer drug leave its binding site?

Obtaining atomistic resolution of drug unbinding from a protein is a much sought-after experimental and computational challenge. We report the unbinding dynamics of the anticancer drug dasatinib from c-Src kinase in full atomistic resolution using enhanced sampling molecular dynamics simulations. We obtain multiple unbinding trajectories and determine a residence time in agreement with experiments. We observe coupled protein-water movement through multiple metastable intermediates. The water molecules form a hydrogen bond bridge, elongating a specific, evolutionarily preserved salt bridge and enabling conformation changes essential to ligand unbinding. This water insertion in the salt bridge acts as a molecular switch that controls unbinding. Our findings provide a mechanistic rationale for why it might be difficult to engineer drugs targeting certain specific c-Src kinase conformations to have longer residence times.

Kinetics for Drug Discovery: an industry-driven effort to target drug residence time

A considerable number of approved drugs show non-equilibrium binding characteristics, emphasizing the potential role of drug residence times for in vivo efficacy. Therefore, a detailed understanding of the kinetics of association and dissociation of a target-ligand complex might provide crucial insight into the molecular mechanism-of-action of a compound. This deeper understanding will help to improve decision making in drug discovery, thus leading to a better selection of interesting compounds to be profiled further. In this review, we highlight the contributions of the Kinetics for Drug Discovery (K4DD) Consortium, which targets major open questions related to binding kinetics in an industry-driven public-private partnership.

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

Understanding the structural and energetic requisites of ligand binding toward its molecular target is of paramount relevance in drug design. In recent years, atomistic free energy calculations have proven to be a valid tool to complement experiments in characterizing the thermodynamic and kinetic properties of protein/ligand interaction. Here, we investigate, through a recently developed metadynamics-based protocol, the unbinding mechanism of an inhibitor of the pharmacologically relevant target p38 MAP kinase. We provide a thorough description of the ligand unbinding pathway identifying the most stable binding mode and other thermodynamically relevant poses. From our simulations, we estimated the unbinding rate as k_off = 0.020 ± 0.011 s^-1. This is in good agreement with the experimental value (k_off = 0.14 s^-1). Next, we developed a Markov state model that allowed identifying the rate-limiting step of the ligand unbinding process. Our calculations further show that the solvation of the ligand and that of the active site play crucial roles in the unbinding process. This study paves the way to investigations on the unbinding dynamics of more complex p38 inhibitors and other pharmacologically relevant inhibitors in general, demonstrating that metadynamics can be a powerful tool in designing new drugs with engineered binding/unbinding kinetics.

Protein-ligand (un)binding kinetics as a new paradigm for drug discovery at the crossroad between experiments and modelling

In the last three decades, protein and nucleic acid structure determination and comprehension of the mechanisms, leading to their physiological and pathological functions, have become a cornerstone of biomedical sciences. A deep understanding of the principles governing the fates of cells and tissue at the molecular level has been gained over the years, offering a solid basis for the rational design of drugs aimed at the pharmacological treatment of numerous diseases. Historically, affinity indicators (i.e. Kd and IC50/EC50) have been assumed to be valid indicators of the in vivo efficacy of a drug. However, recent studies pointed out that the kinetics of the drug-receptor binding process could be as important or even more important than affinity in determining the drug efficacy. This eventually led to a growing interest in the characterisation and prediction of the rate constants of protein-ligand association and dissociation. For instance, a drug with a longer residence time can kinetically select a given receptor over another, even if the affinity for both receptors is comparable, thus increasing its therapeutic index. Therefore, understanding the molecular features underlying binding and unbinding processes is of central interest towards the rational control of drug binding kinetics. In this review, we report the theoretical framework behind protein-ligand association and highlight the latest advances in the experimental and computational approaches exploited to investigate the binding kinetics.

Efficient Atomistic Simulation of Pathways and Calculation of Rate Constants for a Protein-Peptide Binding Process: Application to the MDM2 Protein and an Intrinsically Disordered p53 Peptide

The characterization of protein binding processes - with all of the key conformational changes - has been a grand challenge in the field of biophysics. Here, we have used the weighted ensemble path sampling strategy to orchestrate molecular dynamics simulations, yielding atomistic views of protein-peptide binding pathways involving the MDM2 oncoprotein and an intrinsically disordered p53 peptide. A total of 182 independent, continuous binding pathways were generated, yielding a kon that is in good agreement with experiment. These pathways were generated in 15 days using 3500 cores of a supercomputer, substantially faster than would be possible with "brute force" simulations. Many of these pathways involve the anchoring of p53 residue F19 into the MDM2 binding cleft when forming the metastable encounter complex, indicating that F19 may be a kinetically important residue. Our study demonstrates that it is now practical to generate pathways and calculate rate constants for protein binding processes using atomistic simulation on typical computing resources.

Exploration of gated ligand binding recognizes an allosteric site for blocking FABP4-protein interaction

Fatty acid binding protein 4 (FABP4), reversibly binding to fatty acids and other lipids with high affinities, is a potential target for treatment of cancers. The binding site of FABP4 is buried in an interior cavity and thereby ligand binding/unbinding is coupled with opening/closing of FABP4. It is a difficult task both experimentally and computationally to illuminate the entry or exit pathway, especially with the conformational gating. In this report we combine extensive computer simulations, clustering analysis, and the Markov state model to investigate the binding mechanism of FABP4 and troglitazone. Our simulations capture spontaneous binding and unbinding events as well as the conformational transition of FABP4 between the open and closed states. An allosteric binding site on the protein surface is recognized for the development of novel FABP4 inhibitors. The binding affinity is calculated and compared with the experimental value. The kinetic analysis suggests that ligand residence on the protein surface may delay the binding process. Overall, our results provide a comprehensive picture of ligand diffusion on the protein surface, ligand migration into the buried cavity, and the conformational change of FABP4 at an atomic level.