Organophosphorus (OP) compounds are among the most toxic of chemical substances and widely used as insecticides, pesticides, and chemical warfare agents. The most important enzyme inhibited by OP compounds is acetylcholinesterase (AChe). Inactivation of AChe function results in the accumulation of neurotransmitter, leading to death due to serious respiratory disorders. Organophosphorus hydrolase (OPH), also called phosphotriesterase, is a homo-dimeric metalloenzyme that can hydrolyze various OP agents in the circulatory system, resulting in products that are generally of reduced toxicity. The best OPH substrate found to date is the insecticide diethyl p-nitrophenyl phosphate (paraoxon). Most structural and kinetic studies assume that the binding orientation of paraoxon is identical to that of diethyl 4-methylbenzylphosphonate, which is the only substrate analog co-crystallized with OPH. In the current work, we used a combined docking and molecular dynamics (MD) approach to predict the likely binding mode of paraoxon in the OPH active site. We identified a potential binding mode of paraoxon that does not match the binding mode of diethyl 4-methylbenzylphosphonate. Then, we used the predicted binding mode to run MD simulations on the wild type (WT) OPH complexed with paraoxon, and OPH mutants complexed with paraoxon. Additionally, we identified 3 hot-spot residues (D253, H254, and I255) involved in the stability of the OPH active site. To further assess these predictions, we then experimentally assayed single and double mutants involving these residues (D253E, H254S, I255S, D253E-H254R and D253E-I255G) for hydrolytic activity against paraoxon. Computational structural analysis of protein-substrate dynamics shows different hydrogen bonding profiles for mutants involving D253 (D253E, D253E-H254R, and D253E-I255G) compared to WT OPH. Additionally, the binding free energy calculations and the experimental kinetics (particularly, kcat and KM) of the reactions between each OPH mutant and paraoxon show that mutated forms D253E, D253E-H254R, and D253E-I255G exhibit enhanced activity over WT OPH. Interestingly, our experimental results show that the activity of the double mutant D253E-H254R increased by 19-fold compared to WT OPH.

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.

Organophosphorus hydrolase (OPH) is a metalloenzyme that can hydrolyze organophosphorus agents resulting in products that are generally of reduced toxicity. The best OPH substrate found to date is diethyl p-nitrophenyl phosphate (paraoxon). Most structural and kinetic studies assume that the binding orientation of paraoxon is identical to that of diethyl 4-methylbenzylphosphonate, which is the only substrate analog co-crystallized with OPH. In the current work, we used a combined docking and molecular dynamics (MD) approach to predict the likely binding mode of paraoxon. Then, we used the predicted binding mode to run MD simulations on the wild type (WT) OPH complexed with paraoxon, and OPH mutants complexed with paraoxon. Additionally, we identified three hot-spot residues (D253, H254, and I255) involved in the stability of the OPH active site. We then experimentally assayed single and double mutants involving these residues for paraoxon binding affinity. The binding free energy calculations and the experimental kinetics of the reactions between each OPH mutant and paraoxon show that mutated forms D253E, D253E-H254R, and D253E-I255G exhibit enhanced substrate binding affinity over WT OPH. Interestingly, our experimental results show that the substrate binding affinity of the double mutant D253E-H254R increased by 19-fold compared to WT OPH.

The SAMPL challenges focus on testing and driving progress of computational methods to help guide pharmaceutical drug discovery. However, assessment of methods for predicting binding affinities is often hampered by computational challenges such as conformational sampling, protonation state uncertainties, variation in test sets selected, and even lack of high quality experimental data. SAMPL blind challenges have thus frequently included a component focusing on host-guest binding, which removes some of these challenges while still focusing on molecular recognition. Here, we report on the results of the SAMPL7 blind prediction challenge for host-guest affinity prediction. In this study, we focused on three different host-guest categories-a familiar deep cavity cavitand series which has been featured in several prior challenges (where we examine binding of a series of guests to two hosts), a new series of cyclodextrin derivatives which are monofunctionalized around the rim to add amino acid-like functionality (where we examine binding of two guests to a series of hosts), and binding of a series of guests to a new acyclic TrimerTrip host which is related to previous cucurbituril hosts. Many predictions used methods based on molecular simulations, and overall success was mixed, though several methods stood out. As in SAMPL6, we find that one strategy for achieving reasonable accuracy here was to make empirical corrections to binding predictions based on previous data for host categories which have been studied well before, though this can be of limited value when new systems are included. Additionally, we found that alchemical free energy methods using the AMOEBA polarizable force field had considerable success for the two host categories in which they participated. The new TrimerTrip system was also found to introduce some sampling problems, because multiple conformations may be relevant to binding and interconvert only slowly. Overall, results in this challenge tentatively suggest that further investigation of polarizable force fields for these challenges may be warranted.

X-ray crystallography is the gold standard to resolve conformational ensembles that are significant for protein function, ligand discovery, and computational methods development. However, relevant conformational states may be missed at common cryogenic (cryo) data-collection temperatures but can be populated at room temperature. To assess the impact of temperature on making structural and computational discoveries, we systematically investigated protein conformational changes in response to temperature and ligand binding in a structural and computational workhorse, the T4 lysozyme L99A cavity. Despite decades of work on this protein, shifting to RT reveals new global and local structural changes. These include uncovering an apo helix conformation that is hidden at cryo but relevant for ligand binding, and altered side chain and ligand conformations. To evaluate the impact of temperature-induced protein and ligand changes on the utility of structural information in computation, we evaluated how temperature can mislead computational methods that employ cryo structures for validation. We find that when comparing simulated structures just to experimental cryo structures, hidden successes and failures often go unnoticed. When using structural information in ligand binding predictions, both coarse docking and rigorous binding free energy calculations are influenced by temperature effects. The trend that cryo artifacts limit the utility of structures for computation holds across five distinct protein classes. Our results suggest caution when consulting cryogenic structural data alone, as temperature artifacts can conceal errors and prevent successful computational predictions, which can mislead the development and application of computational methods in discovering bioactive molecules.

The Human Immunodeficiency Virus Type 1 nucleocapsid 7 (NCp7) is a multi-functional protein formed by N-terminal and C-terminal domains surrounding two Zn-fingers, linked by a stretch of basic residues, which play a key role in the viral replication. We report the first NCp7 polarizable molecular dynamics (MD) study using the AMOEBA force field complemented by non-polarizable CHARMM simulations. Specif-ically, we compared the relative free-energy stability of two extreme conformations: a compact one having two aromatic residues from each finger, partially stacked, denoted A; and an unfolded one, with the two residues apart, denoted B. Each of these conformations had been previously experimentally advocated to prevail in solution. We compared their theoretical relative free-energy stability using accelerated MD sampling techniques (Steered MD and Umbrella Sampling) and showed that there was a low free energy difference between them. As A and B do not differ in stability by more than 1-1.5 kcal/mol, they should thus coexist in water solution reconciling earlier NMR experimental findings.

The Drug Design Data Resource (D3R) Grand Challenges present an opportunity to assess, in the context of a blind predictive challenge, the accuracy and the limits of tools and methodologies designed to help guide pharmaceutical drug discovery projects. Here, we report the results of our participation in the D3R Grand Challenge 4, which focused on predicting the binding poses and affinity ranking for compounds targeting the beta-amyloid precursor protein (BACE-1). Our ligand similarity-based protocol using HYBRID (OpenEye Scientific Software) successfully identified poses close to the native binding mode for most of the ligands with less than 2 A RMSD accuracy. Furthermore, we compared the performance of our HYBRID-based approach to that of AutoDock Vina and Dock 6 and found that HYBRID performed better here for pose prediction. We also conducted end-point free energy estimates on protein-ligand complexes using molecular mechanics combined with generalized Born surface area method (MM-GBSA). We found that the binding affinity ranking based on MM-GBSA scores have poor correlation with the experimental values. Finally, the main lessons from our participation in D3R Grand Challenge 4 suggest that: i) the generation of the macrocycles conformers is a key step for successful pose prediction, ii) the protonation states of the BACE-1 binding site should be treated carefully, iii) the MM-GBSA method could not discriminate well between different predicted binding poses, and iv) the MM-GBSA method does not perform well at predicting protein-ligand binding affinities here.

The HIV-1 integrase (IN) is a major target for the design of novel anti-HIV inhibitors. Among these, three inhibitors which embody a halobenzene ring derivative (HR) in their structures are presently used in clinics. High-resolution X-ray crystallography of the complexes of the IN-viral DNA transient complex bound to each of the three inhibitors showed in all cases the HR ring to interact within a confined zone of the viral DNA, limited to the highly conserved 5'CpA 3'/5'TpG 3' step. The extension of its extracyclic CX bond is electron-depleted, owing to the existence of the "sigma-hole." It interacts favorably with the electron-rich rings of base G4. We have sought to increase the affinity of HR derivatives for the G4/C16 base pair. We thus designed thirteen novel derivatives and computed their Quantum Chemistry (QC) intermolecular interaction energies (ΔE) with this base-pair. Most compounds had ΔE values significantly more favorable than those of the HR of the most potent halobenzene drug presently used in clinics, Dolutegravir. This should enable the improvement in a modular piece-wise fashion, the affinities of halogenated inhibitors for viral DNA (vDNA). In view of large scale polarizable molecular dynamics simulations on the entirety of the IN-vDNA-inhibitor complexes, validations of the SIBFA polarizable method are also reported, in which the evolution of each ΔE(SIBFA) contribution is compared to its QC counterpart along this series of derivatives.

Using polarizable (AMOEBA) and non-polarizable (CHARMM) force fields, we compare the relative free-energy stability of two extreme conformations of the HIV-1 NCp7 nucleocapsid that had been previously experimentally advocated to prevail in solution. Using accelerated sampling techniques, we show that they differ in stability by no more than 0.75-1.9 kcal/mol depending on the reference protein sequence. While the extended form appears to be the most probable structure, both forms should thus coexist in water explaining the differing NMR findings.

Molecular docking has been successfully used in computer-aided molecular design projects for the identification of ligand poses within protein binding sites. However, relying on docking scores to rank different ligands with respect to their experimental affinities might not be sufficient. It is believed that the binding scores calculated using molecular mechanics combined with the Poisson–Boltzman surface area (MM-PBSA) or generalized Born surface area (MM-GBSA) can predict binding affinities more accurately. In this perspective, we decided to take part in Stage 2 of the Drug Design Data Resource (D3R) Grand Challenge 4 (GC4) to compare the performance of a quick scoring function, AutoDock4, to that of MM-GBSA in predicting the binding affinities of a set of $$\beta$$-Amyloid Cleaving Enzyme 1 (BACE-1) ligands. Our results show that re-scoring docking poses using MM-GBSA did not improve the correlation with experimental affinities. We further did a retrospective analysis of the results and found that our MM-GBSA protocol is sensitive to details in the protein-ligand system: (i) neutral ligands are more adapted to MM-GBSA calculations than charged ligands, (ii) predicted binding affinities depend on the initial conformation of the BACE-1 receptor, (iii) protonating the aspartyl dyad of BACE-1 correctly results in more accurate binding affinity predictions.

Three Integrase (IN) strand transfer inhibitors are in intensive clinical use, raltegravir, elvitegravir anddolutegravir. However, the onset of IN resistance mutations limits their therapeutic efficiency. As put forth earlier, the drug affinity for the intasome could be improved by targeting preferentially the retroviralnucleobases, which are little, if at all, mutation-prone. We report experimental results of anisotropy fluorescence titrations of viral DNA by these three drugs . These show that the ranking of their inhibitory activities of the intasome corresponds to that of their free energies of binding, D Gs,to retroviral DNA, and that such a ranking is only governed by the binding enthalpies, D H, the entropy undergoing marginal variations.This ranking can therefore be directly correlated to that of model Quantum Chemistry (QC) calculations of intermolecular interaction energies of the sole halobenzene ring with the highly conserved retroviral nucleobases G4 and C14, using Density Functional Theory. This DE(QC) ranking is in turn reproduced by the corresponding DE tot values computed with a polarizable molecular mechanics/dynamics procedure, SIBFA (Sum of Interactions Between Fragments Ab initio computed). Such validations should enable polarizable molecular dynamics simulations on more potent inhibitors in their complexes with the complete intasome. Such derivatives should principally encompass modified halobenzene rings.

Molecular docking has been successfully used in computer-aided molecular design projects for the identification of ligand poses within protein binding sites. However, relying on docking scores to rank different ligands with respect to their experimental affinities might not be sufficient. It is believed that the binding scores calculated using molecular mechanics combined with the Poisson-Boltzman surface area (MM-PBSA) or generalized Born surface area (MM-GBSA) can more accurately predict binding affinities. In this perspective, we decided to take part in Stage 2 in the Drug Design Data Resource (D3R) Grand Challenge 4 (GC4) to compare the performance of a quick scoring function, Autodock4, to that of MM-GBSA in predicting the binding affinities of a set of Beta-Amyloid Cleaving Enzyme 1 (BACE-1) ligands. Our results show that re-scoring docking poses using MM-GBSA did not improve the correlation with experimental affinities. We further did a retrospective analysis of the results and found that our MM-GBSA protocol is sensitive to details in the protein-ligand system: i) neutral ligands are more adapted to MM-GBSA calculations than charged ligands, ii) predicted binding affinities depend on the initial conformation of the BACE-1 receptor, iii) protonating the aspartyl dyad of BACE-1 correctly results in more accurate binding pose and affinity predictions.

The Human Immunodeficiency Virus-1 integrase is responsible for the covalent insertion of a newly synthesized double-stranded viral DNA into the host cells, and is an emerging target for antivirus drug design. Raltegravir (RAL) and elvitegravir (EVG) are the first two integrase strand transfer inhibitors used in therapy. However, treated patients eventually develop detrimental resistance mutations. By contrast, a recently approved drug, dolutegravir (DTG), presents a high barrier to resistance. This study aims to understand the increased efficiency of DTG upon focusing on its interaction properties with viral DNA. The results showed DTG to be involved in more extended interactions with viral DNA than EVG. Such interactions involve the halobenzene and scaffold of DTG and EVG and bases 5'G-43', 3'A35'and 3'C45'.

A correct representation of the short-range contributions such as exchange-repulsion (Erep ) and charge-transfer (Ect ) is essential for the soundness of separable, anisotropic polarizable molecular mechanics potentials. Within the context of the SIBFA procedure, this is aimed at by explicit representations of lone pairs in their expressions. It is necessary to account for their anisotropic behaviors upon performing not only in-plane, but also out-of-plane, variations of a probe molecule or cation interacting with a target molecule or molecular fragment. Thus, Erep and Ect have to reproduce satisfactorily the corresponding anisotropies of their quantum chemical (QC) counterparts. A significant improvement of the out-of-plane dependencies was enabled when the sp2 and sp localized lone-pairs are, even though to a limited extent, delocalized on both sides of the plane, above and below the atom bearer but at the closely similar angles as the in-plane lone pair. We report calibration and validation tests on a series of monoligated complexes of a probe Zn(II) cation with several biochemically relevant ligands. Validations are then performed on several polyligated Zn(II) complexes found in the recognition sites of Zn-metalloproteins. Such calibrations and validations are extended to representative monoligated and polyligated complexes of Mg(II) and Ca(II). It is emphasized that the calibration of all three cations was for each ΔE contribution done on a small training set bearing on a limited number of representative N, O, and S monoligated complexes. Owing to the separable nature of ΔE, a secure transferability is enabled to a diversity of polyligated complexes. For these the relative errors with respect to the target ΔE(QC) values are generally < 3%. Overall, the article proposes a full set of benchmarks that could be useful for force field developers. © 2017 Wiley Periodicals, Inc.

In the context of the SIBFA polarizable molecular mechanics/dynamics (PMM/PMD) procedure, we report the calibration and a series of validation tests for the 1,2,4-triazole-3-thione (TZT) heterocycle. TZT acts as the chelating group of inhibitors of dizinc metallo-β-lactamases (MBL), an emerging class of Zn-dependent bacterial enzymes, which by cleaving the β-lactam bond of most β-lactam antibiotics are responsible for the acquired resistance of bacteria to these drugs. Such a study is indispensable prior to performing PMD simulations of complexes of TZT-based inhibitors with MBL's, on account of the anchoring role of TZT in the dizinc MBL recognition site. Calibration was done by comparisons to energy decomposition analyses (EDA) of high-level ab initio QC computations of the TZT complexes with two probes: Zn(II), representative of "soft" dications, and water, representative of dipolar molecules. We performed distance variations of the approach of each probe to each of the two TZT atoms involved in Zn ligation, the S atom and the N atom ortho to it, so that each SIBFA contribution matches its QC counterpart. Validations were obtained by performing in- and out-of-plane angular variations of Zn(II) binding in monoligated Zn(II)-TZT complexes. The most demanding part of this study was then addressed. How well does ΔE(SIBFA) and its individual contributions compare to their QC counterparts in the dizinc binding site of one MBL, L1, whose structure is known from high-resolution X-ray crystallography? Six distinct complexes were considered, namely each separate monozinc site, and the dizinc site, whether ligated or unligated by TZT. Despite the large magnitude of the interaction energies, in all six complexes ΔE(SIBFA) can match ΔE(QC) with relative errors <2% and the proper balance of individual energy contributions. The computations were extended to the dizinc site of another MBL, VIM-2, and its complexes with two other TZT analogues. ΔE(SIBFA) faithfully reproduced ΔE(QC) in terms of magnitude, ranking of the three ligands, and trends of the separate energy contributions. A preliminary extension to correlated calculations is finally presented. All these validations should enable a secure design of a diversity of TZT-containing MBL inhibitors: a structurally and energetically correct anchoring of TZT should enable all other inhibitor groups to in turn optimize their interactions with the other target MBL residues.

We present a short overview of the recent developments and applications of the SIBFA (Sum of Interactions Between Fragments Ab initio computed) polarizable force field.