Flexible docking in solution using metadynamics

Ligand-Binding Calculations with Metadynamics

All-atom molecular dynamics simulations can capture the dynamic degrees of freedom that characterize molecular recognition, the knowledge of which constitutes the cornerstone of rational approaches to drug design and optimization. In particular, enhanced sampling algorithms, such as metadynamics, are powerful tools to dramatically reduce the computational cost required for a mechanistic description of the binding process. Here, we describe the essential details characterizing these simulation strategies, focusing on the critical step of identifying suitable reaction coordinates, as well as on the different analysis algorithms to estimate binding affinity and residence times. We conclude with a survey of published applications that provides explicit examples of successful simulations for several targets.

In Silico Reoptimization of Binding Affinity and Drug-Resistance Circumvention Ability in Kinase Inhibitors: A Case Study with RL-45 and Src Kinase

A major bottleneck in the development of kinase inhibitors has been the onset of drug resistance around the gatekeeper residues of Src kinase. Although recent times have seen the reports of certain second-generation kinase inhibitors which are capable of bypassing the drug resistance by circumventing kinase mutation, their kinase-binding efficacy has remained considerably weaker than that of the classical adenosine 5'-triphosphate-competitive kinase inhibitors. Using a recently synthesized second-generation kinase inhibitor RL-45 as a template, the current work integrates fragment-based drug discovery and quantitative structure-activity relationship study with enhanced molecular dynamics simulation approaches, namely, metadynamics and replica exchange free-energy perturbation, and demonstrates how one can optimally redesign and assess novel Src kinase inhibitors, by minimal introduction of new functional moieties around template kinase inhibitor. Interestingly, unlike many synthetic kinase inhibitors, these in silico optimized small-molecule derivatives of RL-45 are found to be potentially capable of serving dual purposes, crucial for efficacy of an ideal kinase inhibitor: (a) circumventing gatekeeper residue mutation-related drug resistance in Src kinase, unlike many commercial kinase inhibitors and (b) manifesting superior resilience against unbinding from the kinase active site. The computer simulation, boosted by enhanced sampling techniques, further reveals that these designed inhibitors bring about key interactions in the form of significantly long-standing hydrogen bonds and hydrophobic pocket otherwise weak in the template bioactive kinase inhibitor, which enhance the binding efficacy of these newly designed ligands in the kinase-binding pocket.

Parallel cascade selection molecular dynamics to screen for protein complexes generated by rigid docking

We propose a flexible docking simulation based on parallel cascade selection molecular dynamics (PaCS-MD) as a post-processing treatment after a rigid docking simulation. PaCS-MD has been proposed as an enhanced sampling method for generating structural transition pathways from a given reactant to a product. The PaCS-MD cycle consists of the following two steps: (1) selections of important initial structures and (2) their conformational resampling from the selected initial structures. By repeating the conformational resampling from the important initial structures, structural transitions from the reactant to the product are gradually promoted. In the present flexible docking simulation, decoys (protein complexes) are generated by the rigid docking simulation a priori and employed as products of PaCS-MD. Then PaCS-MD is applied to reproduce association processes to the decoys from a reactant (completely separated proteins). To judge whether PaCS-MD found the association processes or not, the root-mean-square deviation measured from decoy (RMSD_decoy) was defined and monitored during the PaCS-MD cycles. By checking the RMSD_decoy values, a set of decoys is screened as a non-near native protein complex. In more detail, PaCS-MD detects near native protein complexes from the generated decoys by imposing a threshold (cutoff) for RMSD_decoy, i.e. RMSD_decoy < cutoff. As a demonstration, the present flexible docking addressed dimerization processes of K48-linked ubiquitin dimer without a covalent bond between its monomers. Finally, PaCS-MD screened out non-near native protein complexes from decoys generated by a rigid docking simulation.

Does the Lipid Bilayer Orchestrate Access and Binding of Ligands to Transmembrane Orthosteric/Allosteric Sites of G Protein-Coupled Receptors?

The ligand-binding sites of many G protein-coupled receptors (GPCRs) are situated around and deeply embedded within the central pocket formed by their seven transmembrane-spanning α-helical domains. Generally, these binding sites are assumed accessible to endogenous ligands from the aqueous phase. Recent advances in the structural biology of GPCRs, along with biophysical and computational studies, suggest that amphiphilic and lipophilic molecules may gain access to these receptors by first partitioning into the membrane and then reaching the binding site via lateral diffusion through the lipid bilayer. In addition, several crystal structures of class A and class B GPCRs bound to their ligands offer unprecedented details on the existence of lipid-facing allosteric binding sites outside the transmembrane helices that can only be reached via lipid pathways. The highly organized structure of the lipid bilayer may direct lipophilic or amphiphilic drugs to a specific depth within the bilayer, changing local concentration of the drug near the binding site and affecting its binding kinetics. Additionally, the constraints of the lipid bilayer, including its composition and biophysical properties, may play a critical role in "pre-organizing" ligand molecules in an optimal orientation and conformation to facilitate receptor binding. Despite its clear involvement in molecular recognition processes, the critical role of the membrane in binding ligands to lipid-exposed transmembrane binding sites remains poorly understood and warrants comprehensive investigation. Understanding the mechanistic basis of the structure-membrane interaction relationship of drugs will not only provide useful insights about receptor binding kinetics but will also enhance our ability to take advantage of the apparent membrane contributions when designing drugs that target transmembrane proteins with improved efficacy and safety. In this minireview, we summarize recent structural and computational studies on membrane contributions to binding processes, elucidating both lipid pathways of ligand access and binding mechanisms for several orthosteric and allosteric ligands of class A and class B GPCRs.

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.

A Multiscale Simulation Approach to Modeling Drug-Protein Binding Kinetics

Drug-target binding kinetics has recently emerged as a sometimes critical determinant of in vivo efficacy and toxicity. Its rational optimization to improve potency or reduce side effects of drugs is, however, extremely difficult. Molecular simulations can play a crucial role in identifying features and properties of small ligands and their protein targets affecting the binding kinetics, but significant challenges include the long time scales involved in (un)binding events and the limited accuracy of empirical atomistic force fields (lacking, e.g., changes in electronic polarization). In an effort to overcome these hurdles, we propose a method that combines state-of-the-art enhanced sampling simulations and quantum mechanics/molecular mechanics (QM/MM) calculations at the BLYP/VDZ level to compute association free energy profiles and characterize the binding kinetics in terms of structure and dynamics of the transition state ensemble. We test our combined approach on the binding of the anticancer drug Imatinib to Src kinase, a well-characterized target for cancer therapy with a complex binding mechanism involving significant conformational changes. The results indicate significant changes in polarization along the binding pathways, which affect the predicted binding kinetics. This is likely to be of widespread importance in binding of ligands to protein targets.

Exploring the PXR ligand binding mechanism with advanced Molecular Dynamics methods

The Pregnane X Receptor (PXR) is a ligand-activated transcription factor belonging to the nuclear receptor family. PXR can bind diverse drugs and environmental toxicants with different binding modes, making it an intriguing target for drug discovery. Here we investigated the binding mechanism of the SR12813 ligand to elucidate the significant steps, from the ligand entrance pathway into the binding cavity, to the ligand-induced conformational changes, and to the exploration of its alternative binding geometries. We used the advanced Molecular Dynamics-based methods implemented in the BiKi suite and developed specific methodological approaches to overcome the complexity induced by the buried and flexible binding cavity. The adopted methods provided a full dynamic description of the binding event and allowed rationalization of the observed multiple binding modes. These results suggest that the same approach could be exploited for the study of other binding processes with similar characteristics.

Intriguing role of water in protein-ligand binding studied by neutron crystallography on trypsin complexes

Hydrogen bonds are key interactions determining protein-ligand binding affinity and therefore fundamental to any biological process. Unfortunately, explicit structural information about hydrogen positions and thus H-bonds in protein-ligand complexes is extremely rare and similarly the important role of water during binding remains poorly understood. Here, we report on neutron structures of trypsin determined at very high resolutions ≤1.5 Å in uncomplexed and inhibited state complemented by X-ray and thermodynamic data and computer simulations. Our structures show the precise geometry of H-bonds between protein and the inhibitors N-amidinopiperidine and benzamidine along with the dynamics of the residual solvation pattern. Prior to binding, the ligand-free binding pocket is occupied by water molecules characterized by a paucity of H-bonds and high mobility resulting in an imperfect hydration of the critical residue Asp189. This phenomenon likely constitutes a key factor fueling ligand binding via water displacement and helps improving our current view on water influencing protein-ligand recognition.

Bridging Molecular Docking to Molecular Dynamics in Exploring Ligand-Protein Recognition Process: An Overview

Computational techniques have been applied in the drug discovery pipeline since the 1980s. Given the low computational resources of the time, the first molecular modeling strategies relied on a rigid view of the ligand-target binding process. During the years, the evolution of hardware technologies has gradually allowed simulating the dynamic nature of the binding event. In this work, we present an overview of the evolution of structure-based drug discovery techniques in the study of ligand-target recognition phenomenon, going from the static molecular docking toward enhanced molecular dynamics strategies.

Collective Variables from Local Fluctuations

We introduce a method to obtain one-dimensional collective variables for studying rarely occurring transitions between two metastable states separated by a high free energy barrier. No previous information, not even approximated, on the path followed during the transition is needed. The only requirement is to know the fluctuations of the system while in the two metastable states. With this information in hand, we build the collective variable using a modified version of Fisher's linear discriminant analysis. The usefulness of this approach is tested on the metadynamics simulation of two representative systems. The first is the freezing of silver iodide into the superionic α-phase, and the second is the study of a classical Diels-Alder reaction. The collective variable works very well in these two diverse cases.

The Molecular Mechanism Underlying Ligand Binding to the Membrane-Embedded Site of a G-Protein-Coupled Receptor

The crystal structure of P2Y₁ receptor (P2Y₁R), a class A GPCR, revealed a special extra-helical site for its antagonist, BPTU, which locates in-between the membrane and the protein. However, due to the limitation of crystallization experiments, the membrane was mimicked by use of detergents, and the information related to the binding of BPTU to the receptor in the membrane environment is rather limited. In the present work, we conducted a total of ∼7.5 μs all-atom simulations in explicit solvent using conventional molecular dynamics and multiple enhanced sampling methods, with models of BPTU and a POPC bilayer, both in the absence and presence of P2Y₁R. Our simulations revealed that BPTU prefers partitioning into the interface of polar/lipophilic region of the lipid bilayer before associating with the receptor. Then, it interacts with the second extracellular loop of the receptor and reaches the binding site through the lipid-receptor interface. In addition, by use of funnel-metadynamics simulations which efficiently enhance the sampling of bound and unbound states, we provide a statistically accurate description of the underlying binding free energy landscape. The calculated absolute ligand-receptor binding affinity is in excellent agreement with the experimental data (Δ G_{b0_theo} = -11.5 kcal mol^-1, Δ G_{b0_exp}= -11.7 kcal mol^-1). Our study broadens the view of the current experimental/theoretical models and our understanding of the protein-ligand recognition mechanism in the lipid environment. The strategy used in this work is potentially applicable to investigate ligands association/dissociation with other membrane-embedded sites, allowing identification of compounds targeting membrane receptors of pharmacological interest.

Fast Dynamic Docking Guided by Adaptive Electrostatic Bias: The MD-Binding Approach

Engineering chemical entities to modify how pharmaceutical targets function, as it is done in drug design, requires a good understanding of molecular recognition and binding. In this context, the limitations of statically describing bimolecular recognition, as done in docking/scoring, call for insightful and efficient dynamical investigations. On the experimental side, the characterization of dynamical binding processes is still in its infancy. Thus, computer simulations, particularly molecular dynamics (MD), are compelled to play a prominent role, allowing a deeper comprehension of the binding process and its causes and thus a more informed compound selection, making more significant the computational contribution to drug discovery (Carlson, H. A. Curr. Opin. Chem. Biol. 2002, 6, 447-452). Unfortunately, MD-based approaches cannot yet describe complex events without incurring prohibitive time and computational costs. Here, we present a new method for fully and dynamically simulating drug-target-complex formations, tested against a real world and pharmaceutically relevant benchmark set. The method, based on an adaptive, electrostatics-inspired bias, envisions a campaign of trivially parallel short MD simulations and a strategy to identify a near native binding pose from the sampled configurations. At an affordable computational cost, this method provided predictions of good accuracy also when the starting protein conformation was different from that of the crystal complex, a known hurdle for traditional molecular docking (Lexa, K. W.; Carlson, H. A. Q. Rev. Biophys. 2012, 45, 301-343). Moreover, along the observed binding routes, it identified some key features also found by much more computationally expensive plain-MD simulations. Overall, this methodology represents significant progress in the description of binding phenomena.

Dynamic Docking: A Paradigm Shift in Computational Drug Discovery

Molecular docking is the methodology of choice for studying in silico protein-ligand binding and for prioritizing compounds to discover new lead candidates. Traditional docking simulations suffer from major limitations, mostly related to the static or semi-flexible treatment of ligands and targets. They also neglect solvation and entropic effects, which strongly limits their predictive power. During the last decade, methods based on full atomistic molecular dynamics (MD) have emerged as a valid alternative for simulating macromolecular complexes. In principle, compared to traditional docking, MD allows the full exploration of drug-target recognition and binding from both the mechanistic and energetic points of view (dynamic docking). Binding and unbinding kinetic constants can also be determined. While dynamic docking is still too computationally expensive to be routinely used in fast-paced drug discovery programs, the advent of faster computing architectures and advanced simulation methodologies are changing this scenario. It is feasible that dynamic docking will replace static docking approaches in the near future, leading to a major paradigm shift in in silico drug discovery. Against this background, we review the key achievements that have paved the way for this progress.

Impact of protein-ligand solvation and desolvation on transition state thermodynamic properties of adenosine A2A ligand binding kinetics

Ligand-protein binding kinetic rates are growing in importance as parameters to consider in drug discovery and lead optimization. In this study we analysed using surface plasmon resonance (SPR) the transition state (TS) properties of a set of six adenosine A2A receptor inhibitors, belonging to both the xanthine and the triazolo-triazine scaffolds. SPR highlighted interesting differences among the ligands in the enthalpic and entropic components of the TS energy barriers for the binding and unbinding events. To better understand at a molecular level these differences, we developed suMetaD, a novel molecular dynamics (MD)-based approach combining supervised MD and metadynamics. This method allows simulation of the ligand unbinding and binding events. It also provides the system conformation corresponding to the highest energy barrier the ligand is required to overcome to reach the final state. For the six ligands evaluated in this study their TS thermodynamic properties were linked in particular to the role of water molecules in solvating/desolvating the pocket and the small molecules. suMetaD identified kinetic bottleneck conformations near the bound state position or in the vestibule area. In the first case the barrier is mainly enthalpic, requiring the breaking of strong interactions with the protein. In the vestibule TS location the kinetic bottleneck is instead mainly of entropic nature, linked to the solvent behaviour.

Self-Avoiding Conformational Sampling Based on Histories of Past Conformational Searches

Self-avoiding conformational sampling (SACS) is proposed as an enhanced conformational sampling method for proteins. In SACS, the following conformational resampling is repeated for a given protein: (1) identification of newly visited states in a subspace and (2) conformational resampling by restarting short-time molecular dynamics (MD) simulations from the newly visited states. To identify the newly visited states, a set of history-dependent histograms projected onto the subspace is used. One is constructed from the trajectories sampled at the current (ith) cycle, and the other is constructed from all of the trajectories accumulated up through the previous ((i - 1)th) cycle. By reference to the history-dependent histograms, the newly visited states appearing at the current (ith) cycle are defined as a difference set between them. By repeating the cycle of conformational resampling, SACS prevents the system from revisiting states that have already been visited for previous cycles, promoting structural transitions via resampling from the newly visited states. To verify the conformational sampling efficiency of SACS, the present method was applied to reveal underlying mechanisms of biologically important domain motions of maltodextrin binding protein in explicit water and successfully reproduced the open-closed transition with a reasonable (nanosecond-order) computational cost.

Lessons learned in induced fit docking and metadynamics in the Drug Design Data Resource Grand Challenge 2

Two of the major ongoing challenges in computational drug discovery are predicting the binding pose and affinity of a compound to a protein. The Drug Design Data Resource Grand Challenge 2 was developed to address these problems and to drive development of new methods. The challenge provided the 2D structures of compounds for which the organizers help blinded data in the form of 35 X-ray crystal structures and 102 binding affinity measurements and challenged participants to predict the binding pose and affinity of the compounds. We tested a number of pose prediction methods as part of the challenge; we found that docking methods that incorporate protein flexibility (Induced Fit Docking) outperformed methods that treated the protein as rigid. We also found that using binding pose metadynamics, a molecular dynamics based method, to score docked poses provided the best predictions of our methods with an average RMSD of 2.01 Å. We tested both structure-based (e.g. docking) and ligand-based methods (e.g. QSAR) in the affinity prediction portion of the competition. We found that our structure-based methods based on docking with Smina (Spearman ρ = 0.614), performed slightly better than our ligand-based methods (ρ = 0.543), and had equivalent performance with the other top methods in the competition. Despite the overall good performance of our methods in comparison to other participants in the challenge, there exists significant room for improvement especially in cases such as these where protein flexibility plays such a large role.

Adaptive simulations, towards interactive protein-ligand modeling

Modeling the dynamic nature of protein-ligand binding with atomistic simulations is one of the main challenges in computational biophysics, with important implications in the drug design process. Although in the past few years hardware and software advances have significantly revamped the use of molecular simulations, we still lack a fast and accurate ab initio description of the binding mechanism in complex systems, available only for up-to-date techniques and requiring several hours or days of heavy computation. Such delay is one of the main limiting factors for a larger penetration of protein dynamics modeling in the pharmaceutical industry. Here we present a game-changing technology, opening up the way for fast reliable simulations of protein dynamics by combining an adaptive reinforcement learning procedure with Monte Carlo sampling in the frame of modern multi-core computational resources. We show remarkable performance in mapping the protein-ligand energy landscape, being able to reproduce the full binding mechanism in less than half an hour, or the active site induced fit in less than 5 minutes. We exemplify our method by studying diverse complex targets, including nuclear hormone receptors and GPCRs, demonstrating the potential of using the new adaptive technique in screening and lead optimization studies.

Thermodynamics of camphor migration in cytochrome P450cam by atomistic simulations

Understanding the mechanisms of ligand binding to enzymes is of paramount importance for the design of new drugs. Here, we report on the use of a novel biased molecular dynamics (MD) methodology to study the mechanism of camphor binding to cytochrome P450cam. Microsecond-long MD simulations allowed us to observe reaction coordinates characterizing ligand diffusion from the active site of cytochrome P450cam to solvent via three egress routes. These atomistic simulations were used to estimate thermodynamic quantities along the reaction coordinates and indicate diverse binding configurations. The results suggest that the diffusion of camphor along the pathway near the substrate recognition site (SRS) is thermodynamically preferred. In addition, we show that the diffusion near the SRS is triggered by a transition from a heterogeneous collection of closed ligand-bound conformers to the basin comprising the open conformations of cytochrome P450cam. The conformational change accompanying this switch is characterized by the retraction of the F and G helices and the disorder of the B' helix. These results are corroborated by experimental studies and provide detailed insight into ligand binding and conformational behavior of the cytochrome family. The presented methodology is general and can be applied to other ligand-protein systems.

Estimation of kinetic and thermodynamic ligand-binding parameters using computational strategies

Kinetic and thermodynamic ligand–protein binding parameters are gaining growing importance as key information to consider in drug discovery. The determination of the molecular structures, using particularly x-ray and NMR techniques, is crucial for understanding how a ligand recognizes its target in the final binding complex. However, for a better understanding of the recognition processes, experimental studies of ligand–protein interactions are needed. Even though several techniques can be used to investigate both thermodynamic and kinetic profiles for a ligand–protein complex, these procedures are very often laborious, time consuming and expensive. In the last 10 years, computational approaches have enormous potential in providing insights into each of the above effects and in parsing their contributions to the changes in both kinetic and thermodynamic binding parameters. The main purpose of this review is to summarize the state of the art of computational strategies for estimating the kinetic and thermodynamic parameters of a ligand–protein binding.

Multisystem altruistic metadynamics-Well-tempered variant

Metadynamics method has been widely used to enhance sampling in molecular simulations. Its original form suffers two major drawbacks, poor convergence in complex (especially biomolecular) systems and its serial nature. The first drawback has been addressed by introduction of a convergent variant known as well-tempered metadynamics. The second was addressed by introduction of a parallel multisystem metadynamics referred to as altruistic metadynamics. Here, we combine both approaches into well-tempered altruistic metadynamics. We provide mathematical arguments and trial simulations to show that it accurately predicts free energy surfaces.

Accelerated flexible protein-ligand docking using Hamiltonian replica exchange with a repulsive biasing potential

A molecular dynamics replica exchange based method has been developed that allows rapid identification of putative ligand binding sites on the surface of biomolecules. The approach employs a set of ambiguity restraints in replica simulations between receptor and ligand that allow close contacts in the reference replica but promotes transient dissociation in higher replicas. This avoids long-lived trapping of the ligand or partner proteins at nonspecific, sticky, sites on the receptor molecule and results in accelerated exploration of the possible binding regions. In contrast to common docking methods that require knowledge of the binding site, exclude solvent and often keep parts of receptor and ligand rigid the approach allows for full flexibility of binding partners. Application to peptide-protein, protein-protein and a drug-receptor system indicate rapid sampling of near-native binding regions even in case of starting far away from the native binding site outperforming continuous MD simulations. An application on a DNA minor groove binding ligand in complex with DNA demonstrates that it can also be used in explicit solvent simulations.

Metadynamics for Perspective Drug Design: Computationally Driven Synthesis of New Protein-Protein Interaction Inhibitors Targeting the EphA2 Receptor

Metadynamics (META-D) is emerging as a powerful method for the computation of the multidimensional free-energy surface (FES) describing the protein-ligand binding process. Herein, the FES of unbinding of the antagonist N-(3α-hydroxy-5β-cholan-24-oyl)-l-β-homotryptophan (UniPR129) from its EphA2 receptor was reconstructed by META-D simulations. The characterization of the free-energy minima identified on this FES proposes a binding mode fully consistent with previously reported and new structure-activity relationship data. To validate this binding mode, new N-(3α-hydroxy-5β-cholan-24-oyl)-l-β-homotryptophan derivatives were designed, synthesized, and tested for their ability to displace ephrin-A1 from the EphA2 receptor. Among them, two antagonists, namely compounds 21 and 22, displayed high affinity versus the EphA2 receptor and resulted endowed with better physicochemical and pharmacokinetic properties than the parent compound. These findings highlight the importance of free-energy calculations in drug design, confirming that META-D simulations can be used to successfully design novel bioactive compounds.

Exploring Multiple Binding Modes Using Confined Replica Exchange Molecular Dynamics

Molecular docking is extensively applied to determine the position of a ligand on its receptor despite the rather poor correspondence between docking scores and experimental binding affinities found in several studies, especially for systems structurally unrelated with those used in the scoring functions' training sets. Here, we present a method for the prediction of binding modes and binding free energies, which uses replica exchange molecular dynamics in combination with a receptor-shaped piecewise potential, confining the ligand in the proximity of the receptor surface and limiting the accessible conformational space of interest. We assess our methodology with a set of protein receptor-ligand test cases. In every case studied, the method is able to locate the ligand on the experimentally known receptor binding site, and it gives as output the binding free energy. The added value of our approach with respect to other available methods is that it quickly performs a conformational space search, providing a set of bound (or unbound) configurations, which can be used to determine phenomenological structural and energetic properties of an experimental binding state as a result of contributions provided by diversified multiple binding poses.

Resolving the problem of trapped water in binding cavities: prediction of host-guest binding free energies in the SAMPL5 challenge by funnel metadynamics

The funnel metadynamics method enables rigorous calculation of the potential of mean force along an arbitrary binding path and thereby evaluation of the absolute binding free energy. A problem of such physical paths is that the mechanism characterizing the binding process is not always obvious. In particular, it might involve reorganization of the solvent in the binding site, which is not easily captured with a few geometrically defined collective variables that can be used for biasing. In this paper, we propose and test a simple method to resolve this trapped-water problem by dividing the process into an artificial host-desolvation step and an actual binding step. We show that, under certain circumstances, the contribution from the desolvation step can be calculated without introducing further statistical errors. We apply the method to the problem of predicting host–guest binding free energies in the SAMPL5 blind challenge, using two octa-acid hosts and six guest molecules. For one of the hosts, well-converged results are obtained and the prediction of relative binding free energies is the best among all the SAMPL5 submissions. For the other host, which has a narrower binding pocket, the statistical uncertainties are slightly higher; longer simulations would therefore be needed to obtain conclusive results.

Toward High-Throughput Predictive Modeling of Protein Binding/Unbinding Kinetics

One of the unaddressed challenges in drug discovery is that drug potency determined in vitro is not a reliable indicator of drug activity in vivo. Accumulated evidence suggests that in vivo activity is more strongly correlated with the binding/unbinding kinetics than the equilibrium thermodynamics of protein-ligand interactions (PLIs). However, existing experimental and computational techniques are insufficient in studying the molecular details of kinetics processes of PLIs on a large scale. Consequently, we not only have limited mechanistic understanding of the kinetic processes but also lack a practical platform for high-throughput screening and optimization of drug leads on the basis of their kinetic properties. For the first time, we address this unmet need by integrating coarse-grained normal mode analysis with multitarget machine learning (MTML). To test our method, HIV-1 protease is used as a model system. We find that computational models based on the residue normal mode directionality displacement of PLIs can not only recapitulate the results from all-atom molecular dynamics simulations but also predict protein-ligand binding/unbinding kinetics accurately. When this is combined with energetic features, the accuracy of combined kon and koff prediction reaches 74.35%. Furthermore, our integrated model provides us with new insights into the molecular determinants of the kinetics of PLIs. We propose that the coherent coupling of conformational dynamics and thermodynamic interactions between the receptor and the ligand may play a critical role in determining the kinetic rate constants of PLIs. In conclusion, we demonstrate that residue normal mode directionality displacement can serve as a kinetic fingerprint to capture long-time-scale conformational dynamics of the binding/unbinding kinetics. When this is coupled with MTML, it is possible to screen and optimize compounds on the basis of their binding/unbinding kinetics in a high-throughput fashion. The further development of such computational tools will bridge one of the critical missing links between in vitro compound screening and in vivo drug activity.

The Dynamic Process of Drug-GPCR Binding at Either Orthosteric or Allosteric Sites Evaluated by Metadynamics

Major advances in G Protein-Coupled Receptor (GPCR) structural biology over the past few years have yielded a significant number of high-resolution crystal structures for several different receptor subtypes. This dramatic increase in GPCR structural information has underscored the use of automated docking algorithms for the discovery of novel ligands that can eventually be developed into improved therapeutics. However, these algorithms are often unable to discriminate between different, yet energetically similar, poses because of their relatively simple scoring functions. Here, we describe a metadynamics-based approach to study the dynamic process of ligand binding to/unbinding from GPCRs with a higher level of accuracy and yet satisfying efficiency.

Advanced molecular dynamics simulation methods for kinase drug discovery

Interest in the application of molecular dynamics (MD) simulations has increased in the field of protein kinase (PK) drug discovery. PKs belong to an important drug target class because they are directly involved in a number of diseases, including cancer. MD methods simulate dynamic biological and chemical events at an atomic level. This information can be combined with other in silico and experimental methods to efficiently target selected receptors. In this review, we present common and advanced methods of MD simulations and we focus on the recent applications of MD-based methodologies that provided significant insights into the elucidation of biological mechanisms involving PKs and into the discovery of novel kinase inhibitors.

Exploring the unbinding of Leishmania (L.) amazonensis CPB derived-epitopes from H2 MHC class I proteins

New strategies to control Leishmania disease demand an extensive knowledge about several aspects of infection including the understanding of its molecular events. In murine models, cysteine proteinase B from Leishmania amazonensis promotes regulation of immune response, and fragments from its C-terminus extension (cyspep) can play a decisive role in the host-parasite interaction. The interaction between cyspep-derived peptides and major histocompatibility complex (MHC) proteins is a crucial factor in Leishmania infections. Seven cyspep-derived peptides, previously identified as capable of interacting with H-2 (murine) MHC class I proteins, were studied in this work. We established a protocol to simulate the unbinding of these peptides from the cleft of H-2 receptors. From the simulations, we estimated the corresponding free energy of dissociation (ΔGd ) and described the molecular events that occur during the exit of peptides from the cleft. To test the reliability of this method, we first applied it to a calibration set of four crystallographic MHC/peptide complexes. Next, we explored the unbinding of the seven complexes mentioned above. Results were consistent with ΔGd values obtained from surface plasmon resonance (SPR) experiments. We also identified some of the primary interactions between peptides and H-2 receptors, and we detected three regions of influence for the interaction. This pattern was systematically observed for the peptides and helped determine a minimum distance for the real interaction between peptides and H-2 proteins occurring at ∼ 25 Å.

Biased Molecular Simulations for Free-Energy Mapping: A Comparison on the KcsA Channel as a Test Case

The calculation of free-energy landscapes in proteins is a challenge for modern numerical simulations. As to the case of potassium ion channels is concerned, it is particularly interesting because of the nanometric dimensions of the selectivity filter, where the complex electrostatics is highly relevant. The present study aims at comparing three different techniques used to bias molecular dynamics simulations, namely Umbrella Sampling, Steered Molecular Dynamics, and Metadynamics, never applied all together in the past to the same channel protein. Our test case is represented by potassium ions permeating the selectivity filter of the KcsA channel.

Free Energy Barriers for the N-Terminal Asparagine to Succinimide Conversion: Quantum Molecular Dynamics Simulations for the Fully Solvated Model

Deamidation of asparagine residues represents one of the main routes for the post-translational modification of protein sequences. We computed the estimates of the free energy barriers for three stages of the deamidation process, deprotonation, cyclization, and deamination, of the conversion of asparagine to the succinimide intermediate within the fully solvated model with explicit water molecules. The Born-Oppenheimer molecular dynamics in the Gaussian and Plane Wave (GPW) approximation as implemented in the CP2K quantum chemistry package was utilized to sample the configurational space. By applying the metadynamics technique, the estimates of the free energy barriers were obtained for three separated stages of the reaction. In agreement with the experimental kinetic measurements, the estimated activation barriers do not exceed 21 kcal/mol. We demonstrate that the use of fully solvated models is the critical issue in theoretical studies of these reactions. We also conclude that more extensive sampling is necessary to obtain full free energy profiles and accurate barriers for the reaction stages.

Molecular view of ligands specificity for CAG repeats in anti-Huntington therapy

Huntington's disease is a fatal and devastating neurodegenerative genetic disorder for which there is currently no cure. It is characterized by Huntingtin protein's mRNA transcripts with 36 or more CAG repeats. Inhibiting the formation of pathological complexes between these expanded transcripts and target proteins may be a valuable strategy against the disease. Yet, the rational design of molecules specifically targeting the expanded CAG repeats is limited by the lack of structural information. Here, we use well-tempered metadynamics-based free energy calculations to investigate pose and affinity of two ligands targeting CAG repeats for which affinities have been previously measured. The first consists of two 4-guanidinophenyl rings linked by an ester group. It is the most potent ligand identified so far, with Kd = 60(30) nM. The second consists of a 4-phenyl dihydroimidazole and 4-1H-indole dihydroimidazole connected by a C-C bond (Kd = 700(80) nM). Our calculations reproduce the experimental affinities and uncover the recognition pattern between ligands' and their RNA target. They also provide a molecular basis for the markedly different affinity of the two ligands for CAG repeats as observed experimentally. These findings may pave the way for a structure-based hit-to-lead optimization to further improve ligand selectivity toward CAG repeat-containing mRNAs.

Carbohydrate-protein interactions: molecular modeling insights

The article reviews the significant contributions to, and the present status of, applications of computational methods for the characterization and prediction of protein-carbohydrate interactions. After a presentation of the specific features of carbohydrate modeling, along with a brief description of the experimental data and general features of carbohydrate-protein interactions, the survey provides a thorough coverage of the available computational methods and tools. At the quantum-mechanical level, the use of both molecular orbitals and density-functional theory is critically assessed. These are followed by a presentation and critical evaluation of the applications of semiempirical and empirical methods: QM/MM, molecular dynamics, free-energy calculations, metadynamics, molecular robotics, and others. The usefulness of molecular docking in structural glycobiology is evaluated by considering recent docking- validation studies on a range of protein targets. The range of applications of these theoretical methods provides insights into the structural, energetic, and mechanistic facets that occur in the course of the recognition processes. Selected examples are provided to exemplify the usefulness and the present limitations of these computational methods in their ability to assist in elucidation of the structural basis underlying the diverse function and biological roles of carbohydrates in their dialogue with proteins. These test cases cover the field of both carbohydrate biosynthesis and glycosyltransferases, as well as glycoside hydrolases. The phenomenon of (macro)molecular recognition is illustrated for the interactions of carbohydrates with such proteins as lectins, monoclonal antibodies, GAG-binding proteins, porins, and viruses.

Studying the binding interactions of allosteric agonists and antagonists of the CXCR4 receptor

Several examples of allosteric modulators of GPCRs have been reported recently in the literature, but understanding their molecular mechanism presents a new challenge for medicinal chemistry. For the specific case of the cellular receptor CXCR4, it is known that pepducins (lipidated fragments of intracellular GPCR loops) such as ATI-2341 modulate CXCR4 activity agonistically via an allosteric mechanism. Moreover, there are also examples of small organic molecules such as AMD11070 and GSK812397 which may also act as allosteric antagonists. However, incomplete knowledge of the ligand-binding sites has hampered a detailed molecular understanding of how these inhibitors work. Here, we attempt to answer this question by analysing the binding interactions between the CXCR4 receptor and the above-mentioned allosteric modulators. We propose two different allosteric binding sites, one located in the intracellular loops 1, 2 and 3 (ICL1, ICL2 and ICL3) which binds the pepducin agonist ATI-2341, and the other at a subsite of the main extracellular orthosteric binding pocket between extracellular loops 1 and 2 and the N-terminus, which binds the antagonists AMD11070 and GSK812397. Allosteric interactions between the CXCR4 and ATI-2341 were predicted by combining different modeling approaches. First, a rotational blind docking search was applied and the best poses were subsequently refined using flexible docking methods and molecular dynamic simulations. For the AMD11070 and GSK812397 antagonists, the entire CXCR4 protein surface was explored by blind docking in order to define the binding region. A second docking analysis by subsites was then performed to refine the allosteric interactions. Finally, we identified the binding residues that appear to be essential for CXCR4 allosteric modulators.