Computation of Absolute Binding Free Energies for Noncovalent Inhibitors with SARS-CoV-2 Main Protease.

Accurate, routine calculation of absolute binding free energies (ABFEs) for protein-ligand complexes remains a key goal of computer-aided drug design since it can enable screening and optimization of drug candidates. For development and testing of related methods, it is important to have high-quality datasets. To this end, from our own experimental studies, we have selected a set of 16 inhibitors of the SARS-CoV-2 main protease (Mpro) with structural diversity and well-distributed BFEs covering a 5 kcal/mol range. There is also minimal structural uncertainty since X-ray crystal structures have been deposited for 12 of the compounds. For methods testing, we report ABFE results from 2 µs molecular dynamics (MD) simulations using free energy perturbation (FEP) theory. The correlation of experimental and computed results is encouraging, with a Pearson's r2 of 0.58 and a Kendall τ of 0.24. The results indicate that current FEP-based ABFE calculations can be used for identification of active compounds (hits). While their accuracy for lead optimization is not yet sufficient, this activity remains addressable in separate lead series by relative BFE calculations.

Assessing Metadynamics and Docking for Absolute Binding Free Energy Calculations Using Severe Acute Respiratory Syndrome Coronavirus 2 Main Protease Inhibitors.

Absolute binding free energy (ABFE) calculations can be an important part of the drug discovery process by identifying molecules that have the potential to be strong binders for a biomolecular target. Recent work has used free energy perturbation (FEP) theory for these calculations, focusing on a set of 16 inhibitors of the severe acute respiratory syndrome coronavirus 2 main protease (Mpro). Herein, the same data set is evaluated by metadynamics (MetaD), four different docking programs, and molecular mechanics with generalized Born and surface area solvation. MetaD yields a Kendall τ distance of 0.28 and Pearson r2 of 0.49, which reflect somewhat less accuracy than that from the ABFE FEP results. Notably, it is demonstrated that an ensemble docking protocol by which each ligand is docked into the 13 crystal structures in this data set provides improved performance, particularly when docking is carried out with Glide XP (Kendall τ distance = 0.20, Pearson r2 = 0.71), Glide SP (Kendall τ distance = 0.19, Pearson r2 = 0.66), or AutoDock 4 (Kendall τ distance = 0.21, Pearson r2 = 0.55). The best results are obtained with "superconsensus" docking by averaging the 52 results for each compound using the 4 docking protocols and all 13 crystal structures (Kendall τ distance = 0.18, Pearson r2 = 0.73).

Ensemble Geometric Deep Learning of Aqueous Solubility.

Geometric deep learning is one of the main workhorses for harnessing the power of big data to predict molecular properties such as aqueous solubility, which is key to the pharmacokinetic improvement of drug candidates. Two ensembles of graph neural network architectures were built, one based on spectral convolution and the other on spatial convolution. The pretrained models, denoted respectively as SolNet-GCN and SolNet-GAT, significantly outperformed the existing neural networks benchmarked on a validation set of 207 molecules. The SolNet-GCN model demonstrated the best performance on both the training and validation sets, with RMSE values of 0.53 and 0.72 log molar unit and Pearson r2 values of 0.95 and 0.75, respectively. Further, the ranking power of the SolNet models agreed well with a QM-based thermodynamic cycle approach at the PBE-vdW level of theory on a series of benzophenylurea derivatives and a series of benzodiazepine derivatives. Nevertheless, testing the resultant models on a set of inhibitors of the macrophage migration inhibitory factor (MIF) illustrated that the inclusion of atomic attributes to discriminate atoms with a higher tendency to form intermolecular hydrogen bonds in the crystalline state and to identify planar or nonplanar substructures can be beneficial for the prediction of aqueous solubility.

Metadynamics as a Postprocessing Method for Virtual Screening with Application to the Pseudokinase Domain of JAK2.

With standard scoring methods, top-ranked compounds from virtual screening by docking often turn out to be inactive. For this reason, metadynamics, a method used to sample rare events, was studied to further evaluate docking poses with the aim of reducing false positives. Specifically, virtual screening was performed with Glide SP to seek potential molecules to bind to the ATP site in the pseudokinase domain of JAK2 kinase, and promising compounds were selected from the top-ranked 1000 based on visualization. Rescoring with Glide XP, GOLD, and MM/GBSA was unable to differentiate well between active and inactive compounds. Metadynamics was then used to gauge the relative binding affinity from the required time or the potential of mean force needed to dissociate the ligand from the bound complex. With consideration of previously known binders of varying affinities, metadynamics was able to differentiate between the most active compounds and inactive or weakly active ones, and it could identify correctly most of the selected virtual screening compounds as false positives. Thus, metadynamics has the potential to be a viable postprocessing method for virtual screening, minimizing the expense of buying or synthesizing inactive compounds.

QM/MM Calculations for the Cl- + CH3Cl SN2 Reaction in Water Using CM5 Charges and Density Functional Theory.

The prototypical SN2 reaction of chloride ion with methyl chloride has been reinvestigated in aqueous solution using QM/MM methodology featuring MO6-2X/6-31+G(d) calculations with the TIP4P water model, and partial charges were computed with the CM5 method. Though the DFT method yields excellent gas-phase energetics for the reaction, the QM/MM approach is found to yield overestimation of the activation barrier by ca. 12 kcal/mol. The discrepancy is traced to underestimate of the magnitude of the partial charges on the chlorine atoms in the transition structure. When CM1 or CM3 charges based on semiempirical wave functions are used instead, the agreement with experiment is much improved. The findings emphasize the sensitivity of the results of QM/MM calculations to the choice of QM method, the MM force field, and implementation of the QM/MM interface.

LigParGen web server: an automatic OPLS-AA parameter generator for organic ligands.

The accurate calculation of protein/nucleic acid-ligand interactions or condensed phase properties by force field-based methods require a precise description of the energetics of intermolecular interactions. Despite the progress made in force fields, small molecule parameterization remains an open problem due to the magnitude of the chemical space; the most critical issue is the estimation of a balanced set of atomic charges with the ability to reproduce experimental properties. The LigParGen web server provides an intuitive interface for generating OPLS-AA/1.14*CM1A(-LBCC) force field parameters for organic ligands, in the formats of commonly used molecular dynamics and Monte Carlo simulation packages. This server has high value for researchers interested in studying any phenomena based on intermolecular interactions with ligands via molecular mechanics simulations. It is free and open to all at jorgensenresearch.com/ligpargen, and has no login requirements.

Improved Treatment of Nucleosides and Nucleotides in the OPLS-AA Force Field.

DFT calculations have been used to develop improved descriptions of the torsional energetics for nucleosides and nucleotides in the OPLS-AA force field. Scans of nucleotide dihedral angles (γ, χ, and ß) and methyl phosphates provided the bases for the new torsional parameters. In addition, the angle-bending parameters of phosphodiesters and ribose were updated, and adjustments were made to existing carbohydrate torsions to better capture the sugar puckering landscape of ribose. MD simulations of nucleosides with the new parameters demonstrate a significant improvement in the ribose sugar puckering and χ angle distributions. Additionally, energy-minimization of protein-nucleotide crystal structures with the new parameters produced accurate poses.

Molecular dynamics and Monte Carlo simulations for protein-ligand binding and inhibitor design.

BACKGROUND: Non-nucleoside inhibitors of HIV reverse transcriptase are an important component of treatment against HIV infection. Novel inhibitors are sought that increase potency against variants that contain the Tyr181Cys mutation. METHODS: Molecular dynamics based free energy perturbation simulations have been run to study factors that contribute to protein-ligand binding, and the results are compared with those from previous Monte Carlo based simulations and activity data. RESULTS: Predictions of protein-ligand binding modes are very consistent for the two simulation methods; the accord is attributed to the use of an enhanced sampling protocol. The Tyr181Cys binding pocket supports large, hydrophobic substituents, which is in good agreement with experiment. CONCLUSIONS: Although some discrepancies exist between the results of the two simulation methods and experiment, free energy perturbation simulations can be used to rapidly test small molecules for gains in binding affinity. GENERAL SIGNIFICANCE: Free energy perturbation methods show promise in providing fast, reliable and accurate data that can be used to complement experiment in lead optimization projects. This article is part of a Special Issue entitled "Recent developments of molecular dynamics".

Application of a BOSS-Gaussian interface for QM/MM simulations of Henry and methyl transfer reactions.

Hybrid quantum mechanics and molecular mechanics (QM/MM) computer simulations have become an indispensable tool for studying chemical and biological phenomena for systems too large to treat with QM alone. For several decades, semiempirical QM methods have been used in QM/MM simulations. However, with increased computational resources, the introduction of ab initio and density function methods into on-the-fly QM/MM simulations is being increasingly preferred. This adaptation can be accomplished with a program interface that tethers independent QM and MM software packages. This report introduces such an interface for the BOSS and Gaussian programs, featuring modification of BOSS to request QM energies and partial atomic charges from Gaussian. A customizable C-shell linker script facilitates the interprogram communication. The BOSS-Gaussian interface also provides convenient access to Charge Model 5 (CM5) partial atomic charges for multiple purposes including QM/MM studies of reactions. In this report, the BOSS-Gaussian interface is applied to a nitroaldol (Henry) reaction and two methyl transfer reactions in aqueous solution. Improved agreement with experiment is found by determining free-energy surfaces with MP2/CM5 QM/MM simulations than previously reported investigations using semiempirical methods.

Determination of partial molar volumes from free energy perturbation theory.

Partial molar volume is an important thermodynamic property that gives insights into molecular size and intermolecular interactions in solution. Theoretical frameworks for determining the partial molar volume (V°) of a solvated molecule generally apply Scaled Particle Theory or Kirkwood-Buff theory. With the current abilities to perform long molecular dynamics and Monte Carlo simulations, more direct methods are gaining popularity, such as computing V° directly as the difference in computed volume from two simulations, one with a solute present and another without. Thermodynamically, V° can also be determined as the pressure derivative of the free energy of solvation in the limit of infinite dilution. Both approaches are considered herein with the use of free energy perturbation (FEP) calculations to compute the necessary free energies of solvation at elevated pressures. Absolute and relative partial molar volumes are computed for benzene and benzene derivatives using the OPLS-AA force field. The mean unsigned error for all molecules is 2.8 cm(3) mol(-1). The present methodology should find use in many contexts such as the development and testing of force fields for use in computer simulations of organic and biomolecular systems, as a complement to related experimental studies, and to develop a deeper understanding of solute-solvent interactions.

OPLS/2020 Force Field for Unsaturated Hydrocarbons, Alcohols, and Ethers.

The OPLS all-atom force field was updated and applied to modeling unsaturated hydrocarbons, alcohols, and ethers. Testing has included gas-phase conformational energetics, properties of pure liquids, and free energies of hydration. Monte Carlo statistical mechanics (MC) calculations were used to model 60 liquids. In addition, a robust, automated procedure was devised to compute the free energies of hydration with high precision via free-energy perturbation (FEP) calculations using double annihilation. Testing has included larger molecules than in the past, and parameters are reported for the first time for some less common groups including alkynes, allenes, dienes, and acetals. The average errors in comparison with experimental data for the computed properties of the pure liquids were improved with the modified force field (OPLS/2020). For liquid densities and heats of vaporization, the average unsigned errors are 0.01 g/cm3 and 0.2 kcal/mol. The average error and signed error for free energies of hydration are both 1.2 kcal/mol. As noted before, this reflects a systematic overestimate of the hydrophobicity of organic molecules when the parametrization is done to minimize the errors for properties of pure liquids. Implications for the modeling of biomolecular systems with standard force fields are considered.

Optimization of benzyloxazoles as non-nucleoside inhibitors of HIV-1 reverse transcriptase to enhance Y181C potency.

Design of non-nucleoside inhibitors of HIV-1 reverse transcriptase with improved activity towards Tyr181Cys containing variants was pursued with the assistance of free energy perturbation (FEP) calculations. Optimization of the 4-R substituent in 1 led to ethyl and isopropyl analogs 1e and 1f with 1-7 nM potency towards both the wild-type virus and a Tyr181C variant.

Characterization of biaryl torsional energetics and its treatment in OPLS all-atom force fields.

The frequency of biaryl substructures in a database of approved oral drugs has been analyzed. This led to designation of 20 prototypical biaryls plus 10 arylpyridinones for parametrization in the OPLS all-atom force fields. Bond stretching, angle-bending, and torsional parameters were developed to reproduce the MP2 geometries and torsional energy profiles. The transferability of the new parameters was tested through their application to three additional biaryls. The torsional energetics for the 33 biaryl molecules are analyzed and factors leading to preferences for planar and nonplanar geometries are identified. For liquid biphenyl, the computed density and heat of vaporization at the boiling point (255 °C) are also reported.

Refinement of the Optimized Potentials for Liquid Simulations Force Field for Thermodynamics and Dynamics of Liquid Alkanes.

Torsion and Lennard-Jones parameters of the optimized potentials for liquid simulations (OPLS) all-atom force field have been refined for describing thermodynamics and dynamics of a wide range of liquid alkanes. Monte Carlo statistical mechanics (MC) and molecular dynamics (MD) simulations were carried out. For thermodynamics properties, MC simulations with truncated electrostatic interactions performed very closely to MD simulations with a Verlet neighbor list and the particle mesh Ewald algorithm. The average errors in comparison with experimental data for computed properties were improved with the modified force field (OPLS/2020), especially for long-chain alkanes. For liquid densities, heats of vaporization, and free energies of hydration, the average errors are 0.01 g/cm3, 0.2 kcal/mol, and ca. 0.5 kcal/mol, respectively; significant gains were made for relative heats of vaporization of isomeric series. Results for self-diffusion coefficients also reproduce experimental data well for linear alkane liquids up to hexadecane. The new force field is suitable for use in improved modeling of myriad systems of importance in chemistry, biology, and materials science.

Potent Noncovalent Inhibitors of the Main Protease of SARS-CoV-2 from Molecular Sculpting of the Drug Perampanel Guided by Free Energy Perturbation Calculations.

Starting from our previous finding of 14 known drugs as inhibitors of the main protease (Mpro) of SARS-CoV-2, the virus responsible for COVID-19, we have redesigned the weak hit perampanel to yield multiple noncovalent, nonpeptidic inhibitors with ca. 20 nM IC50 values in a kinetic assay. Free-energy perturbation (FEP) calculations for Mpro-ligand complexes provided valuable guidance on beneficial modifications that rapidly delivered the potent analogues. The design efforts were confirmed and augmented by determination of high-resolution X-ray crystal structures for five analogues bound to Mpro. Results of cell-based antiviral assays further demonstrated the potential of the compounds for treatment of COVID-19. In addition to the possible therapeutic significance, the work clearly demonstrates the power of computational chemistry for drug discovery, especially FEP-guided lead optimization.

Optimization of Triarylpyridinone Inhibitors of the Main Protease of SARS-CoV-2 to Low-Nanomolar Antiviral Potency.

Non-covalent inhibitors of the main protease (Mpro) of SARS-CoV-2 having a pyridinone core were previously reported with IC50 values as low as 0.018 µM for inhibition of enzymatic activity and EC50 values as low as 0.8 µM for inhibition of viral replication in Vero E6 cells. The series has now been further advanced by consideration of placement of substituted five-membered-ring heterocycles in the S4 pocket of Mpro and N-methylation of a uracil ring. Free energy perturbation calculations provided guidance on the choice of the heterocycles, and protein crystallography confirmed the desired S4 placement. Here we report inhibitors with EC50 values as low as 0.080 µM, while remdesivir yields values of 0.5-2 µM in side-by-side testing with infectious SARS-CoV-2. A key factor in the improvement is enhanced cell permeability, as reflected in PAMPA measurements. Compounds 19 and 21 are particularly promising as potential therapies for COVID-19, featuring IC50 values of 0.044-0.061 µM, EC50 values of ca. 0.1 µM, good aqueous solubility, and no cytotoxicity.

Quantum mechanical/molecular mechanical modeling finds Diels-Alder reactions are accelerated less on the surface of water than in water.

Quantum and molecular mechanics calculations for the Diels-Alder reactions of cyclopentadiene with 1,4-naphthoquinone, methyl vinyl ketone, and acrylonitrile have been carried out at the vacuum-water interface and in the gas phase. In conjunction with previous studies of these cycloadditions in dilute solution, a more complete picture of aqueous environmental effects emerges with implications for the origin of observed rate accelerations using heterogeneous aqueous suspensions, "on water" conditions. The pure TIP4P water slab maintains the bulk density and hydrogen-bonding properties in central water layers. The bulk region merges to vacuum over a ca. 5 A band with progressive diminution of the density and hydrogen bonding. The relative free energies of activation and transition structures for the reactions at the interface are found to be intermediate between those calculated in the gas phase and in bulk water; i.e., for the reaction with 1,4-naphthoquinone, the DeltaDeltaG(++) values relative to the gas phase are -3.6 and -7.3 kcal/mol at the interface and in bulk water, respectively. Thus, the results do not support the notion that a water surface is more effective than bulk water for catalysis of such pericyclic reactions. The trend is in qualitative agreement with expectations based on density considerations and estimates of experimental rate constants for the gas phase, a heterogeneous aqueous suspension, and a dilute aqueous solution for the reaction of cyclopentadiene with methyl vinyl ketone. Computed energy pair distributions reveal a uniform loss of 0.5-1.0 hydrogen bond for the reactants and transition states in progressing from bulk water to the vacuum-water interface. Orientational effects are apparent at the surface; e.g., the carbonyl group in the methyl vinyl ketone transition structure is preferentially oriented into the surface. Also, the transition structure for the 1,4-naphthoquinone case is buried more in the surface, and the free energy of activation for this reaction is most similar to the result in bulk water.

Explicit Representation of Cation-π Interactions in Force Fields with 1/r4 Nonbonded Terms.

The binding energies for cation-π complexation are underestimated by traditional fixed-charge force fields owing to their lack of explicit treatment of ion-induced dipole interactions. To address this deficiency, an explicit treatment of cation-π interactions has been introduced into the OPLS-AA force field. Following prior work with atomic cations, it is found that cation-π interactions can be handled efficiently by augmenting the usual 12-6 Lennard-Jones potentials with 1/r4 terms. Results are provided for prototypical complexes as well as protein-ligand systems of relevance for drug design. Alkali cation, ammonium, guanidinium, and tetramethylammonium were chosen for the representative cations, while benzene and six heteroaromatic molecules were used as the π systems. The required nonbonded parameters were fit to reproduce structure and interaction energies for gas-phase complexes from density functional theory (DFT) calculations at the ωB97X-D/6-311++G(d,p) level. The impact of the solvent was then examined by computing potentials of mean force (pmfs) in both aqueous and tetrahydrofuran (THF) solutions using the free-energy perturbation (FEP) theory. Further testing was carried out for two cases of strong and one case of weak cation-π interactions between druglike molecules and their protein hosts, namely, the JH2 domain of JAK2 kinase and macrophage migration inhibitory factor. FEP results reveal greater binding by 1.5-4.4 kcal/mol from the addition of the explicit cation-π contributions. Thus, in the absence of such treatment of cation-π interactions, errors for computed binding or inhibition constants of 101-103 are expected.

Identification of 14 Known Drugs as Inhibitors of the Main Protease of SARS-CoV-2.

A consensus virtual screening protocol has been applied to ca. 2000 approved drugs to seek inhibitors of the main protease (M pro ) of SARS-CoV-2, the virus responsible for COVID-19. 42 drugs emerged as top candidates, and after visual analyses of the predicted structures of their complexes with M pro , 17 were chosen for evaluation in a kinetic assay for M pro inhibition. Remarkably 14 of the compounds at 100-µM concentration were found to reduce the enzymatic activity and 5 provided IC 50 values below 40 µM: manidipine (4.8 µM), boceprevir (5.4 µM), lercanidipine (16.2 µM), bedaquiline (18.7 µM), and efonidipine (38.5 µM). Structural analyses reveal a common cloverleaf pattern for the binding of the active compounds to the P1, P1', and P2 pockets of M pro . Further study of the most active compounds in the context of COVID-19 therapy is warranted, while all of the active compounds may provide a foundation for lead optimization to deliver valuable chemotherapeutics to combat the pandemic.

Identification of 14 Known Drugs as Inhibitors of the Main Protease of SARS-CoV-2.

A consensus virtual screening protocol has been applied to ca. 2000 approved drugs to seek inhibitors of the main protease (Mpro) of SARS-CoV-2, the virus responsible for COVID-19. 42 drugs emerged as top candidates, and after visual analyses of the predicted structures of their complexes with Mpro, 17 were chosen for evaluation in a kinetic assay for Mpro inhibition. Remarkably 14 of the compounds at 100-µM concentration were found to reduce the enzymatic activity and 5 provided IC50 values below 40 µM: manidipine (4.8 µM), boceprevir (5.4 µM), lercanidipine (16.2 µM), bedaquiline (18.7 µM), and efonidipine (38.5 µM). Structural analyses reveal a common cloverleaf pattern for the binding of the active compounds to the P1, P1', and P2 pockets of Mpro. Further study of the most active compounds in the context of COVID-19 therapy is warranted, while all of the active compounds may provide a foundation for lead optimization to deliver valuable chemotherapeutics to combat the pandemic.