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

Increased throughput in methods for simulating protein ligand binding and unbinding

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

Comprehensive Survey of Consensus Docking for High-Throughput Virtual Screening

The rapid advances of 3D techniques for the structural determination of proteins and the development of numerous computational methods and strategies have led to identifying highly active compounds in computer drug design. Molecular docking is a method widely used in high-throughput virtual screening campaigns to filter potential ligands targeted to proteins. A great variety of docking programs are currently available, which differ in the algorithms and approaches used to predict the binding mode and the affinity of the ligand. All programs heavily rely on scoring functions to accurately predict ligand binding affinity, and despite differences in performance, none of these docking programs is preferable to the others. To overcome this problem, consensus scoring methods improve the outcome of virtual screening by averaging the rank or score of individual molecules obtained from different docking programs. The successful application of consensus docking in high-throughput virtual screening highlights the need to optimize the predictive power of molecular docking methods.

Using Selectively Scaled Molecular Dynamics Simulations to Assess Ligand Poses in RNA Aptamers

Predicting the structure (or pose) of RNA-ligand complexes is an important problem in RNA structural biology. Although one could use computational docking to rapidly sample putative poses of RNA-ligand complexes, accurately discriminating the native-like poses from non-native, decoy poses remains a formidable challenge. Here, we started from the assumption that native-like RNA-ligand poses are less likely to dissociate during molecular dynamics simulations, and then we used enhanced simulations to promote ligand unbinding for diverse poses of a handful of RNA aptamer-ligand complexes. By fitting unbinding profiles obtained from the simulations to a single exponential, we identified the characteristic decay time (τ) as particularly effective at resolving native poses from decoys. We also found that a simple regression model trained to predict the simulation-derived parameters directly from structure could also discriminate ligand poses for similar RNA aptamers. Characterizing the unbinding properties of individual poses may thus be an effective strategy for enhancing pose prediction for ligands interacting with RNA aptamers. A similar strategy might be applicable to other ligandable RNAs; however, further analysis will be required to confirm this hypothesis.

Spatial organization of hydrophobic and charged residues affects protein thermal stability and binding affinity

What are the molecular determinants of protein–protein binding affinity and whether they are similar to those regulating fold stability are two major questions of molecular biology, whose answers bring important implications both from a theoretical and applicative point of view. Here, we analyze chemical and physical features on a large dataset of protein–protein complexes with reliable experimental binding affinity data and compare them with a set of monomeric proteins for which melting temperature data was available. In particular, we probed the spatial organization of protein (1) intramolecular and intermolecular interaction energies among residues, (2) amino acidic composition, and (3) their hydropathy features. Analyzing the interaction energies, we found that strong Coulombic interactions are preferentially associated with a high protein thermal stability, while strong intermolecular van der Waals energies correlate with stronger protein–protein binding affinity. Statistical analysis of amino acids abundances, exposed to the molecular surface and/or in interaction with the molecular partner, confirmed that hydrophobic residues present on the protein surfaces are preferentially located in the binding regions, while charged residues behave oppositely. Leveraging on the important role of van der Waals interface interactions in binding affinity, we focused on the molecular surfaces in the binding regions and evaluated their shape complementarity, decomposing the molecular patches in the 2D Zernike basis. For the first time, we quantified the correlation between local shape complementarity and binding affinity via the Zernike formalism. In addition, considering the solvent interactions via the residue hydropathy, we found that the hydrophobicity of the binding regions dictates their shape complementary as much as the correlation between van der Waals energy and binding affinity. In turn, these relationships pave the way to the fast and accurate prediction and design of optimal binding regions as the 2D Zernike formalism allows a rapid and superposition-free comparison between possible binding surfaces.

Practical Protocols for Efficient Sampling of Kinase-Inhibitor Binding Pathways Using Two-Dimensional Replica-Exchange Molecular Dynamics

Molecular dynamics (MD) simulations are increasingly used to study various biological processes such as protein folding, conformational changes, and ligand binding. These processes generally involve slow dynamics that occur on the millisecond or longer timescale, which are difficult to simulate by conventional atomistic MD. Recently, we applied a two-dimensional (2D) replica-exchange MD (REMD) method, which combines the generalized replica exchange with solute tempering (gREST) with the replica-exchange umbrella sampling (REUS) in kinase-inhibitor binding simulations, and successfully observed multiple ligand binding/unbinding events. To efficiently apply the gREST/REUS method to other kinase-inhibitor systems, we establish modified, practical protocols with non-trivial simulation parameter tuning. The current gREST/REUS simulation protocols are tested for three kinase-inhibitor systems: c-Src kinase with PP1, c-Src kinase with Dasatinib, and c-Abl kinase with Imatinib. We optimized the definition of kinase-ligand distance as a collective variable (CV), the solute temperatures in gREST, and replica distributions and umbrella forces in the REUS simulations. Also, the initial structures of each replica in the 2D replica space were prepared carefully by pulling each ligand from and toward the protein binding sites for keeping stable kinase conformations. These optimizations were carried out individually in multiple short MD simulations. The current gREST/REUS simulation protocol ensures good random walks in 2D replica spaces, which are required for enhanced sampling of inhibitor dynamics around a target kinase.

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

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

Enhanced Molecular Dynamics Method to Efficiently Increase the Discrimination Capability of Computational Protein-Protein Docking

Protein–protein docking typically consists of the generation of putative binding conformations, which are subsequently ranked by fast heuristic scoring functions. The simplicity of these functions allows for computational efficiency but has severe repercussions on their discrimination capabilities. In this work, we show the effectiveness of suitable descriptors calculated along short scaled molecular dynamics runs in recognizing the nearest-native bound conformation among a set of putative structures generated by the HADDOCK tool for eight protein–protein systems.

Immunoproteasome and Non-Covalent Inhibition: Exploration by Advanced Molecular Dynamics and Docking Methods

The selective inhibition of immunoproteasome is a valuable strategy to treat autoimmune, inflammatory diseases, and hematologic malignancies. Recently, a new series of amide derivatives as non-covalent inhibitors of the β1i subunit with K_i values in the low/submicromolar ranges have been identified. Here, we investigated the binding mechanism of the most potent and selective inhibitor, N-benzyl-2-(2-oxopyridin-1(2H)-yl)propanamide (1), to elucidate the steps from the ligand entrance into the binding pocket to the ligand-induced conformational changes. We carried out a total of 400 ns of MD-binding analyses, followed by 200 ns of plain MD. The trajectories clustering allowed identifying three representative poses evidencing new key interactions with Phe31 and Lys33 together in a flipped orientation of a representative pose. Further, Binding Pose MetaDynamics (BPMD) studies were performed to evaluate the binding stability, comparing 1 with four other inhibitors of the β1i subunit: N-benzyl-2-(2-oxopyridin-1(2H)-yl)acetamide (2), N-cyclohexyl-3-(2-oxopyridin-1(2H)-yl)propenamide (3), N-butyl-3-(2-oxopyridin-1(2H)-yl)propanamide (4), and (S)-2-(2-oxopyridin-1(2H)-yl)-N,4-diphenylbutanamide (5). The obtained results in terms of free binding energy were consistent with the experimental values of inhibition, confirming 1 as a lead compound of this series. The adopted methods provided a full dynamic description of the binding events, and the information obtained could be exploited for the rational design of new and more active inhibitors.

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

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

Observing spontaneous, accelerated substrate binding in molecular dynamics simulations of glutamate transporters

Glutamate transporters are essential for removing the neurotransmitter glutamate from the synaptic cleft. Glutamate transport across the membrane is associated with elevator-like structural changes of the transport domain. These structural changes require initial binding of the organic substrate to the transporter. Studying the binding pathway of ligands to their protein binding sites using molecular dynamics (MD) simulations requires micro-second level simulation times. Here, we used three methods to accelerate aspartate binding to the glutamate transporter homologue Gltph and to investigate the binding pathway. 1) Two methods using user-defined forces to prevent the substrate from diffusing too far from the binding site. 2) Conventional MD simulations using very high substrate concentrations in the 0.1 M range. The final, substrate bound states from these methods are comparable to the binding pose observed in crystallographic studies, although they show more flexibility in the side chain carboxylate function. We also captured an intermediate on the binding pathway, where conserved residues D390 and D394 stabilize the aspartate molecule. Finally, we investigated glutamate binding to the mammalian glutamate transporter, excitatory amino acid transporter 1 (EAAT1), for which a crystal structure is known, but not in the glutamate-bound state. Overall, the results obtained in this study reveal new insights into the pathway of substrate binding to glutamate transporters, highlighting intermediates on the binding pathway and flexible conformational states of the side chain, which most likely become locked in once the hairpin loop 2 closes to occlude the substrate.

Tuning Local Hydration Enables a Deeper Understanding of Protein-Ligand Binding: The PP1-Src Kinase Case

Water plays a key role in biomolecular recognition and binding. Despite the development of several computational and experimental approaches, it is still challenging to comprehensively characterize water-mediated effects on the binding process. Here, we investigate how water affects the binding of Src kinase to one of its inhibitors, PP1. Src kinase is a target for treating several diseases, including cancer. We use biased molecular dynamics simulations, where the hydration of predetermined regions is tuned at will. This computational technique efficiently accelerates the SRC-PP1 binding simulation and allows us to identify several key and yet unexplored aspects of the solvent's role. This study provides a further perspective on the binding phenomenon, which may advance the current drug design approaches for the development of new kinase inhibitors.

Thermodynamics and Kinetics of Drug-Target Binding by Molecular Simulation

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

Combined strategies in structure-based virtual screening

The identification and optimization of lead compounds are inalienable components in drug design and discovery pipelines. As a powerful computational approach for the identification of hits with novel structural scaffolds, structure-based virtual screening (SBVS) has exhibited a remarkably increasing influence in the early stages of drug discovery. During the past decade, a variety of techniques and algorithms have been proposed and tested with different purposes in the scope of SBVS. Although SBVS has been a common and proven technology, it still shows some challenges and problems that are needed to be addressed, where the negative influence regardless of protein flexibility and the inaccurate prediction of binding affinity are the two major challenges. Here, focusing on these difficulties, we summarize a series of combined strategies or workflows developed by our group and others. Furthermore, several representative successful applications from recent publications are also discussed to demonstrate the effectiveness of the combined SBVS strategies in drug discovery campaigns.

DFMD: Fast and Effective DelPhiForce Steered Molecular Dynamics Approach to Model Ligand Approach Toward a Receptor: Application to Spermine Synthase Enzyme

Here we report a novel approach, the DelPhiForce Molecular Dynamics (DFMD) method, for steered molecular dynamics simulations to model receptor-ligand association involving charged species. The main purpose of developing DFMD is to simulate ligand's trajectory toward the receptor and thus to predict the "entrance" of the binding pocket and conformational changes associated with the binding. We demonstrate that the DFMD is superior compared with molecular dynamics simulations applying standard cut-offs, provides correct binding forces, allows for modeling the ligand approach at long distances and thus guides the ligand toward the correct binding spot, and it is very fast (frequently the binding is completed in <1 ns). The DFMD is applied to model the binding of two ligands to a receptor (spermine synthase) and it is demonstrated that it guides the ligands toward the corresponding pockets despite of the initial ligand's position with respect to the receptor. Predicted conformational changes and the order of ligand binding are experimentally verified.

Encounter complexes and hidden poses of kinase-inhibitor binding on the free-energy landscape

Modern drug discovery increasingly focuses on the drug-target binding kinetics which depend on drug (un)binding pathways. The conventional molecular dynamics simulation can observe only a few binding events even using the fastest supercomputer. Here, we develop 2D gREST/REUS simulation with enhanced flexibility of the ligand and the protein binding site. Simulation (43 μs in total) applied to an inhibitor binding to c-Src kinase covers 100 binding and unbinding events. On the statistically converged free-energy landscapes, we succeed in predicting the X-ray binding structure, including water positions. Furthermore, we characterize hidden semibound poses and transient encounter complexes on the free-energy landscapes. Regulatory residues distant from the catalytic core are responsible for the initial inhibitor uptake and regulation of subsequent bindings, which was unresolved by experiments. Stabilizing/blocking of either the semibound poses or the encounter complexes can be an effective strategy to optimize drug-target residence time.

Natural Alkaloid Berberine Activity against Pseudomonas aeruginosa MexXY-Mediated Aminoglycoside Resistance: In Silico and in Vitro Studies

The multidrug efflux system MexXY-OprM, inside the resistance-nodulation-division family, is a major determinant of aminoglycoside resistance in Pseudomonas aeruginosa. In the fight aimed to identify potential efflux pump inhibitors among natural compounds, the alkaloid berberine emerged as a putative inhibitor of MexXY-OprM. In this work, we elucidated its interaction with the extrusor protein MexY and assessed its synergistic activity with aminoglycosides. In particular, we built an in silico model for the MexY protein in its trimeric association using both AcrB (E. coli) and MexB (P. aeruginosa) as 3D templates. This model has been stabilized in the bacterial cytoplasmic membrane using a molecular dynamics approach and used for ensemble docking to obtain the binding site mapping. Then, through dynamic docking, we assessed its binding affinity and its synergism with aminoglycosides focusing on tobramycin, which is widely used in the treatment of pulmonary infections. In vitro assays validated the data obtained: the results showed a 2-fold increase of the inhibitory activity and 2-4 log increase of the killing activity of the association berberine-tobramycin compared to those of tobramycin alone against 13/28 tested P. aeruginosa clinical isolates. From hemolytic assays, we preliminarily assessed berberine's low toxicity.

Ligand Pathways in Nuclear Receptors

Nuclear receptors (NRs) are ligand-inducible transcription factors that play an essential role in a multitude of physiological processes as well as diseases, rendering them attractive drug targets. Crystal structures revealed the binding site of NRs to be buried in the core of the protein, with no obvious route for ligands to access this cavity. The process of ligand binding is known to be an often-neglected contributor to the efficacy of drug candidates and is thought to influence the selectivity and specificity of NRs. While experimental methods generally fail to highlight the dynamic processes of ligand access or egress on the atomistic scale, computational methods have provided fundamental insight into the pathways connecting the buried binding pocket to the surrounding environment. Methods based on molecular dynamics (MD) and Monte Carlo simulations have been applied to identify pathways and quantify their capability to transport ligands. Here, we systematically review findings of more than 20 years of research in the field, including the applied methodology and controversies. Further, we establish a unified nomenclature to describe the pathways with respect to their location relative to protein secondary structure elements and summarize findings relevant to drug design. Lastly, we discuss the effect of NR interaction partners such as coactivators and corepressors, as well as mutations on the pathways.

Discovering New Casein Kinase 1d Inhibitors with an Innovative Molecular Dynamics Enabled Virtual Screening Workflow

The value of including protein flexibility in structure-based drug design (SBDD) is widely documented, and currently, molecular dynamics (MD) simulations represent a powerful tool to investigate protein dynamics. Yet, the inclusion of MD-derived information in pre-existing SBDD workflows is still far from trivial. We recently published an integrated MD-FLAP (Fingerprints for Ligands and Proteins) approach combining MD, clustering and Linear Discriminant Analysis (LDA) for enhancing accuracy, efficacy, and for protein conformational selection in virtual screening (VS) campaigns. Here we prospectively applied the MD-FLAP workflow to discover novel chemotypes inhibiting the Casein Kinase 1 delta (CSNK1D) enzyme. We first obtained a VS model able to separate active from inactive compounds, with a global AUC of 0.9 and a partial ROC enrichment at 0.5% of 0.18, and use it to mine the internal Pfizer screening database. Seven active molecules sharing a phenyl-indazole scaffold, not yet reported among CSNK1D inhibitors, were found. The most potent inhibitor showed an IC₅₀ of 134 nM.

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.

Binding kinetics of cariprazine and aripiprazole at the dopamine D3 receptor

The dissociation behaviours of aripiprazole and cariprazine at the human D₂ and D₃ receptor are evaluated. A potential correlation between kinetics and in vivo profiles, especially cariprazine’s action on negative symptoms in schizophrenia, is investigated. The binding kinetics of four ligands were indirectly evaluated. After the receptor preparations were pre-incubated with the unlabelled ligands, the dissociation was initiated with an excess of [³H]spiperone. Slow dissociation kinetics characterizes aripiprazole and cariprazine at the D₂ receptor. At the D₃ receptor, aripiprazole exhibits a slow monophasic dissociation, while cariprazine displays a rapid biphasic behaviour. Functional ß-arrestin assays and molecular dynamics simulations at the D₃ receptor confirm a biphasic binding behaviour of cariprazine. This may influence its in vivo action, as the partial agonist could react rapidly to variations in the dopamine levels of schizophrenic patients and the ligand will not quantitatively dissociate from the receptor in one single step. With these findings novel agents may be developed that display rapid, biphasic dissociation from the D₃R to further investigate this effect on in vivo profiles.

Predicting how drug molecules bind to their protein targets

There have been substantial advances in the application of molecular modelling and simulation to drug discovery in recent years, as massive increases in computer power are coupled with continued development in the underlying methods and understanding of how to apply them. Here, we survey recent advances in one particular area-predicting how a known ligand binds to a particular protein. We focus on the four contributing classes of calculation: predicting where a binding site is on a protein; characterizing where chemical functional groups will bind to that site; molecular docking to generate a binding mode for a ligand and dynamics simulations to refine that pose and allow for protein conformation change. Examples of successful application are provided for each class.

BiKi Life Sciences: A New Suite for Molecular Dynamics and Related Methods in Drug Discovery

In this paper, we introduce the BiKi Life Sciences suite. This software makes it easy for computational medicinal chemists to run ad hoc molecular dynamics protocols in a novel and task-oriented environment; as a notebook, BiKi (acronym of Binding Kinetics) keeps memory of any activity together with dependencies among them. It offers unique accelerated protein-ligand binding/unbinding methods and other useful tools to gain actionable knowledge from molecular dynamics simulations and to simplify the drug discovery process.