Janus kinase type 3 (JAK3), an emerging target for treating autoimmune diseases, possesses a front pocket cysteine that is targeted by covalent modifiers, best represented by the marketed drug ritlecitinib (1). Recently, 2,3-dihydro-1H-inden-1-ylcyanamides have been developed as novel JAK3 inhibitors. Among them, the N-(6-(7H-pyrrolo[2,3-d]pyrimidin-4-yl)-2,3-dihydro-1H-inden-1-yl)cyanamide inhibitor (2) and its methylated analogue (3), while being potent inhibitors, displayed different mechanisms of action (covalent vs noncovalent) and binding modes (Casimiro-Garcia et al., J Med Chem 2018). Prompted by this intriguing behavior, we applied a multiscale approach to characterize the reaction mechanism between the JAK3 front-pocket Cys909 and cyanamide-based inhibitors. Quantum mechanics/molecular mechanics simulations showed that 2 can readily form an isothiourea adduct with the Cys909 only when a conserved water molecule assists the reaction as a proton shuttle and that methylation of the 2,3-dihydro-1H-inden-1-ylcyanamide moiety of 2 hampers the isothiourea formation by displacing this water molecule. Metadynamics and thermodynamic integration simulations were applied to investigate the relative abundance of alternative poses accessible to 2,3-dihydro-1H-inden-1-ylcyanamides, explaining the effect of methylation on the relative binding mode preference. This multiscale approach provides new chemical insights into the mechanism of action of cyanamide inhibitors and emerges as an effective protocol to investigate the interaction between drugs and molecular targets.

The residence time (RT), the time for which a drug remains bound to its biological target, is a critical parameter for drug design. The prediction of this key kinetic property has been proven to be challenging and computationally demanding in the framework of atomistic simulations. In the present work, we setup and applied two distinct metadynamics protocols to estimate the RTs of muscarinic M3 receptor antagonists. In the first method, derived from the conformational flooding approach, the kinetics of unbinding is retrieved from a physics-based parameter known as the acceleration factor α (i.e., the running average over time of the potential deposited in the bound state). Such an approach is expected to recover the absolute RT value for a compound of interest. In the second method, known as the tMETA-D approach, a qualitative estimation of the RT is given by the time of simulation required to drive the ligand from the binding site to the solvent bulk. This approach has been developed to reproduce the change of experimental RTs for compounds targeting the same target. Our analysis shows that both computational protocols are able to rank compounds in agreement with their experimental RTs. Quantitative structure–kinetics relationship (SKR) models can be identified and employed to predict the impact of a chemical modification on the experimental RT once a calibration study has been performed.

We report the fast-track computationally-driven discovery of new SARS-CoV2 Main Protease (Mpro) inhibitors whose potency range from mM for initial non-covalent ligands to high nM for the final covalent compound (IC50=830 +/-50 nM). The project extensively relied on high-resolution all-atom molecular dynamics simulations and absolute binding free energy calculations performed using the polarizable AMOEBA force field. The study is complemented by extensive adaptive sampling simulations used to rationalize different ligands binding poses through the explicit reconstruction of the ligand-protein conformational space. Machine learning predictions are also utilized to predict selected compound properties. Computations were performed on GPU-accelerated supercomputers and high-performance cloud infrastructures to exponentially reduce time-to-solution, and were systematically coupled to nuclear magnetic resonance experiments to drive synthesis and in vitro characterization of compounds. The study highlights the power of in silico strategies that rely on structure-based approaches for drug design and address protein conformational heterogeneity. The proposed scaffolds open a path toward further optimization of Mpro inhibitors with nM affinities.

While a plethora of different protein-ligand docking protocols have been developed over the past twenty years, their performances greatly depend on the provided input protein-ligand pair. In this work we have developed a machine-learning model that uses a combination of convolutional and fully-connected neural networks for the task of predicting the performance of several popular docking protocols given a protein structure and a small compound. We also rigorously evaluate the performance of our model using a widely available database of protein-ligand complexes and different types of data splits. We further open-source all code related to this study so that potential users can make informed guesses on which protocol is best suited for their particular protein-ligand pair.

SkeleDock is a scaffold docking algorithm which uses the structure of a protein-ligand complex as a template to model the binding mode of a chemically similar system. This algorithm was evaluated in the D3R Grand Challenge 4 pose prediction challenge, where it achieved competitive performance. Furthermore, we show that, if crystallized fragments of the target ligand are available, SkeleDock can outperform rDock docking software at predicting the binding mode. This article also addresses the capacity of this algorithm to model macrocycles and deal with scaffold hopping. SkeleDock can be accessed at https://playmolecule.org/SkeleDock/.

The number of entries in the Protein Data Bank (PDB) has doubled in the last decade, and it has increased tenfold in the last twenty years. The availability of an ever-growing number of structures is having a huge impact on the Structure-Based Drug Discovery (SBDD), allowing investigation of new targets and giving the possibility to have multiple structures of the same macromolecule in a complex with different ligands. Such a large resource often implies the choice of the most suitable complex for molecular docking calculation, and this task is complicated by the plethora of possible posing and scoring function algorithms available, which may influence the quality of the outcomes. Here, we report a large benchmark performed on the PDBbind database containing more than four thousand entries and seventeen popular docking protocols. We found that, even in protein families wherein docking protocols generally showed acceptable results, certain ligand-protein complexes are poorly reproduced in the self-docking procedure. Such a trend in certain protein families is more pronounced, and this underlines the importance in identification of a suitable protein-ligand conformation coupled to a well-performing docking protocol.

G-protein coupled receptors (GPCRs) play a pivotal role in transmitting signals at the cellular level. Structural insights can be exploited to support GPCR structure-based drug discovery endeavours. Despite advances in GPCR crystallography, active state structures are scarce. Molecular dynamics (MD) simulations have been used to explore the conformational landscape of GPCRs. Efforts have been made to retrieve active state conformations starting from inactive structures, however to date this has not been possible without using an energy bias. Here, we reconstruct the activation pathways of the apo adenosine receptor (A2A), starting from an inactive conformation, by applying adaptive sampling MD combined with a goal-oriented scoring function. The reconstructed pathways reconcile well with experiments and help deepen our understanding of A2A regulatory mechanisms. Exploration of the apo conformational landscape of A2A reveals the existence of ligand-competent states, active intermediates and state-dependent cholesterol hotspots of relevance for drug discovery. To the best of our knowledge this is the first time an activation process has been elucidated for a GPCR starting from an inactive structure only, using a non-biased MD approach, opening avenues for the study of ligand binding to elusive yet pharmacologically relevant GPCR states.

Drug discovery suffers from high attrition because compounds initially deemed as promising can later show ineffectiveness or toxicity resulting from a poor understanding of their activity profile. In this work, we describe a deep self-normalizing neural network model for the prediction of molecular pathway association and evaluate its performance, showing an AUC ranging from 0.69 to 0.91 on a set of compounds extracted from ChEMBL and from 0.81 to 0.83 on an external data set provided by Novartis. We finally discuss the applicability of the proposed model in the domain of lead discovery. A usable application is available via PlayMolecule.org .

The Cover Feature shows the different behavior of structural water molecules and bulk water. We have developed the AquaMMapS tool, which is able to identify stationary hydration sites close to the protein surface starting from molecular dynamics simulation. Root-mean-square fluctuation (i.e., a 1.4 Å cutoff) is employed to discriminate stationary from nonstationary water molecules. Space is organized into a grid, where cells crossed by stationary water molecules are selected and their occupancy computed along the simulation. Moreover, an empirical scoring function, the AquaMMapScore, has been developed to evaluate the penalty of a ligand displacing a stationary water hot spot. More information can be found in the Full Paper by Stefano Moro et al. on page 522 in Issue 6, 2018 (DOI: 10.1002/cmdc.201700564).

The serine-threonine checkpoint kinase 1 (Chk1) plays a critical role in the cell cycle arrest in response to DNA damage. In the last decade, Chk1 inhibitors have emerged as a novel therapeutic strategy to potentiate the anti-tumour efficacy of cytotoxic chemotherapeutic agents. In the search for new Chk1 inhibitors, a congeneric series of 2-aryl-2 H-pyrazolo[4,3-c]quinolin-3-one (PQ) was evaluated by in-vitro and in-silico approaches for the first time. A total of 30 PQ structures were synthesised in good to excellent yields using conventional or microwave heating, highlighting that 14 of them are new chemical entities. Noteworthy, in this preliminary study two compounds 4e2 and 4h2 have shown a modest but significant reduction in the basal activity of the Chk1 kinase. Starting from these preliminary results, we have designed the second generation of analogous in this class and further studies are in progress in our laboratories.

The reliability of physics-based in-silico studies of protein-ligand complexes highly depends on the quality of available structures and force-field parameters. Both these subjects have been largely addressed by both experimental and computational scientists from industry and academia. Yet, tasks like obtaining an initial structure with the correct protonation states and hydrogen-bond network or accurate force-field parameters for a given ligand can still be out of reach for the non-experts in those particular fields. Here we showcase two software tools that aim at bridging this gap: proteinPrepare [1,2] and parameterize. We show how these softwares can be easily used by the community and how we are integrating these tools within a wider computational pipeline for drug discovery. Stefan Doerr, Toni Giorgino, Gerard Martínez-Rosell, João M. Damas, and Gianni De Fabritiis. High-Throughput Automated Preparation and Simulation of Membrane Proteins with HTMD. Journal of Chemical Theory and Computation 2017 13 (9), 4003-4011. DOI: 10.1021/acs.jctc.7b00480 Gerard Martínez-Rosell, Toni Giorgino, and Gianni De Fabritiis. PlayMolecule ProteinPrepare: A Web Application for Protein Preparation for Molecular Dynamics Simulations. Journal of Chemical Information and Modeling 2017 57 (7), 1511-1516. DOI: 10.1021/acs.jcim.7b00190

Unquestionably, water appears to be an active player in the noncovalent protein-ligand binding process, as it can either bridge interactions between protein and ligand or can be replaced by the bound ligand. Accordingly, in the last decade, alternative computational methodologies have been sought with the aim of predicting the position and thermodynamic profile of water molecules (i.e., hydration sites) in the binding site using either the ligand-bound or ligand-free protein conformation. Herein, we present an alternative approach, named AquaMMapS, that provides a three-dimensional sampling of putative hydration sites. Interestingly, AquaMMapS can post-inspect molecular dynamics (MD) trajectories obtained from different MD engines using indifferently crystallographic or docking-driven structures as a starting point. Moreover, AquaMMapS is naturally integrated into supervised molecular dynamics (SuMD) simulations, presenting the possibility to inspect hydration sites during the ligand-protein association process. Finally, a penalty scoring method, named AquaMMapScoring(AMS), was developed to evaluate the number and nature of the water molecules displaced by a ligand approaching its binding site during the binding event, guiding a medicinal chemist to explore the most suitable regions of a ligand that can be decorated either with or without interfering with the interaction networks mediated by water molecules with specific recognition regions of the protein.

Peptides have gained increased interest as therapeutic agents during recent years. The high specificity and relatively low toxicity of peptide drugs derive from their extremely tight binding to their targets. Indeed, understanding the molecular mechanism of protein-peptide recognition has important implications in the fields of biology, medicine, and pharmaceutical sciences. Even if crystallography and nuclear magnetic resonance are offering valuable atomic insights into the assembling of the protein-peptide complexes, the mechanism of their recognition and binding events remains largely unclear. In this work we report, for the first time, the use of a supervised molecular dynamics approach to explore the possible protein-peptide binding pathways within a timescale reduced up to three orders of magnitude compared with classical molecular dynamics. The better and faster understating of the protein-peptide recognition pathways could be very beneficial in enlarging the applicability of peptide-based drug design approaches in several biotechnological and pharmaceutical fields.

Molecular docking is a powerful tool in the field of computer-aided molecular design. In particular, it is the technique of choice for the prediction of a ligand pose within its target binding site. A multitude of docking methods is available nowadays, whose performance may vary depending on the data set. Therefore, some non-trivial choices should be made before starting a docking simulation. In the same framework, the selection of the target structure to use could be challenging, since the number of available experimental structures is increasing. Both issues have been explored within this work. The pose prediction of a pool of 36 compounds provided by D3R Grand Challenge 2 organizers was preceded by a pipeline to choose the best protein/docking-method couple for each blind ligand. An integrated benchmark approach including ligand shape comparison and cross-docking evaluations was implemented inside our DockBench software. The results are encouraging and show that bringing attention to the choice of the docking simulation fundamental components improves the results of the binding mode predictions.

Molecular recognition is a crucial issue when aiming to interpret the mechanism of known active substances as well as to develop novel active candidates. Unfortunately, simulating the binding process is still a challenging task because it requires classical MD experiments in a long microsecond time scale that are affordable only with a high-level computational capacity. In order to overcome this limiting factor, we have recently implemented an alternative MD approach, named supervised molecular dynamics (SuMD), and successfully applied it to G protein-coupled receptors (GPCRs). SuMD enables the investigation of ligand-receptor binding events independently from the starting position, chemical structure of the ligand, and also from its receptor binding affinity. In this article, we present an extension of the SuMD application domain including different types of proteins in comparison with GPCRs. In particular, we have deeply analyzed the ligand-protein recognition pathways of six different case studies that we grouped into two different classes: globular and membrane proteins. Moreover, we introduce the SuMD-Analyzer tool that we have specifically implemented to help the user in the analysis of the SuMD trajectories. Finally, we emphasize the limit of the SuMD applicability domain as well as its strengths in analyzing the complexity of ligand-protein recognition pathways.

Structure-based drug design (SBDD) has matured within the last two decades as a valuable tool for the optimization of low molecular weight lead compounds to highly potent drugs. The key step in SBDD requires knowledge of the three-dimensional structure of the target-ligand complex, which is usually determined by X-ray crystallography. In the absence of structural information for the complex, SBDD relies on the generation of plausible molecular docking models. However, molecular docking protocols suffer from inaccuracies in the description of the interaction energies between the ligand and the target molecule, and often fail in the prediction of the correct binding mode. In this context, the appropriate selection of the most accurate docking protocol is absolutely relevant for the final molecular docking result, even if addressing this point is absolutely not a trivial task. D3R Grand Challenge 2015 has represented a precious opportunity to test the performance of DockBench, an integrate informatics platform to automatically compare RMDS-based molecular docking performances of different docking/scoring methods. The overall performance resulted in the blind prediction are encouraging in particular for the pose prediction task, in which several complex were predicted with a sufficient accuracy for medicinal chemistry purposes.

In this review, we present a survey of the recent advances carried out by our research groups in the field of ligand-GPCRs recognition process simulations recently implemented at the Molecular Modeling Section (MMS) of the University of Padova. We briefly describe a platform of tools we have tuned to aid the identification of novel GPCRs binders and the better understanding of their binding mechanisms, based on two extensively used computational techniques such as molecular docking and MD simulations. The developed methodologies encompass: (i) the selection of suitable protocols for docking studies, (ii) the exploration of the dynamical evolution of ligand-protein interaction networks, (iii) the detailed investigation of the role of water molecules upon ligand binding, and (iv) a glance at the way the ligand might go through prior reaching the binding site.

The search for G protein-coupled receptors (GPCRs) allosteric modulators represents an active research field in medicinal chemistry. Allosteric modulators usually exert their activity only in the presence of the orthosteric ligand by binding to protein sites topographically different from the orthosteric cleft. They therefore offer potentially therapeutic advantages by selectively influencing tissue responses only when the endogenous agonist is present. The prediction of putative allosteric site location, however, is a challenging task. In facts, they are usually located in regions showing more structural variation among the family members. In the present work, we applied the recently developed Supervised Molecular Dynamics (SuMD) methodology to interpret at the molecular level the positive allosteric modulation mediated by LUF6000 toward the human adenosine A3 receptor (hA3 AR). Our data suggest at least two possible mechanisms to explain the experimental data available. This study represent, to the best of our knowledge, the first case reported of an allosteric recognition mechanism depicted by means of molecular dynamics simulations.

ALK inhibitor crizotinib has shown potent antitumor activity in children with refractory Anaplastic Large Cell Lymphoma (ALCL) and the opportunity to include ALK inhibitors in first-line therapies is oncoming. However, recent studies suggest that crizotinib-resistance mutations may emerge in ALCL patients. In the present study, we analyzed ALK kinase domain mutational status of 36 paediatric ALCL patients at diagnosis to identify point mutations and gene aberrations that could impact on NPM-ALK gene expression, activity and sensitivity to small-molecule inhibitors. Amplicon ultra-deep sequencing of ALK kinase domain detected 2 single point mutations, R335Q and R291Q, in 2 cases, 2 common deletions of exon 23 and 25 in all the patients, and 7 splicing-related INDELs in a variable number of them. The functional impact of missense mutations and INDELs was evaluated. Point mutations were shown to affect protein kinase activity, signalling output and drug sensitivity. INDELs, instead, generated kinase-dead variants with dominant negative effect on NPM-ALK kinase, in virtue of their capacity of forming non-functional heterocomplexes. Consistently, when co-expressed, INDELs increased crizotinib inhibitory activity on NPM-ALK signal processing, as demonstrated by the significant reduction of STAT3 phosphorylation. Functional changes in ALK kinase activity induced by both point mutations and structural rearrangements were resolved by molecular modelling and dynamic simulation analysis, providing novel insights into ALK kinase domain folding and regulation. Therefore, these data suggest that NPM-ALK pre-therapeutic mutations may be found at low frequency in ALCL patients. These mutations occur randomly within the ALK kinase domain and affect protein activity, while preserving responsiveness to crizotinib.