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.

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.

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.

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.

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.

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.

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.

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.

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.

Absolute Free Energy of Binding Calculations for Macrophage Migration Inhibitory Factor in Complex with a Druglike Inhibitor.

Calculation of the absolute free energy of binding (ΔGbind) for a complex in solution is challenging owing to the need for adequate configurational sampling and an accurate energetic description, typically with a force field (FF). In this study, Monte Carlo (MC) simulations with improved side-chain and backbone sampling are used to assess ΔGbind for the complex of a druglike inhibitor (MIF180) with the protein macrophage migration inhibitory factor (MIF) using free energy perturbation (FEP) calculations. For comparison, molecular dynamics (MD) simulations were employed as an alternative sampling method for the same system. With the OPLS-AA/M FF and CM5 atomic charges for the inhibitor, the ΔGbind results from the MC/FEP and MD/FEP simulations, -8.80 ± 0.74 and -8.46 ± 0.85 kcal/mol, agree well with each other and with the experimental value of -8.98 ± 0.28 kcal/mol. The convergence of the results and analysis of the trajectories indicate that sufficient sampling was achieved for both approaches. Repeating the MD/FEP calculations using current versions of the CHARMM and AMBER FFs led to a 6 kcal/mol range of computed ΔGbind. These results show that calculation of accurate ΔGbind for large ligands is both feasible and numerically equivalent, within error limits, using either methodology.

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.

Robust Free Energy Perturbation Protocols for Creating Molecules in Solution.

Accurate methods to estimate free energies play an important role for studying diverse condensed-phase problems in chemistry and biochemistry. The most common methods used in conjunction with molecular dynamics (MD) and Monte Carlo statistical mechanics (MC) simulations are free energy perturbation (FEP) and thermodynamic integration (TI). For common applications featuring the conversion of one molecule to another, simulations are run in stages or multiple "λ-windows" to promote convergence of the results. For computation of absolute free energies of solvation or binding, calculations are needed in which the solute is typically annihilated in the solvent and in the complex. The present work addresses identification of optimal protocols for such calculations, specifically, the creation/annihilation of organic molecules in aqueous solution. As is common practice, decoupling of the perturbations for electrostatic and Lennard-Jones interactions was performed. Consistent with earlier reports, FEP calculations for molecular creations are much more efficient, while annihilations require many more windows and may converge to incorrect values. Strikingly, we find that as few as four windows may be adequate for creation calculations for solutes ranging from argon to ethylbenzene. For a larger druglike molecule, MIF180, which contains 22 non-hydrogen atoms and three rotatable bonds, 10 creation windows are found to be adequate to yield the correct free energy of hydration. Convergence is impeded with procedures that use any sampling in the annihilation direction, and there is no need for postprocessing methods such as the Bennett acceptance ratio (BAR).

Development and Testing of the OPLS-AA/M Force Field for RNA.

Significant improvements have been made to the OPLS-AA force field for modeling RNA. New torsional potentials were optimized based on density functional theory (DFT) scans at the ωB97X-D/6-311++G(d,p) level for potential energy surfaces of the backbone α and γ dihedral angles. In combination with previously reported improvements for the sugar puckering and glycosidic torsion terms, the new force field was validated through diverse molecular dynamics simulations for RNAs in aqueous solution. Results for dinucleotides and tetranucleotides revealed both accurate reproduction of 3 J couplings from NMR and the avoidance of several unphysical states observed with other force fields. Simulations of larger systems with noncanonical motifs showed significant structural improvements over the previous OPLS-AA parameters. The new force field, OPLS-AA/M, is expected to perform competitively with other recent RNA force fields and to be compatible with OPLS-AA models for proteins and small molecules.

Unbinding Dynamics of Non-Nucleoside Inhibitors from HIV-1 Reverse Transcriptase.

Non-nucleoside inhibitors of HIV-1 reverse transcriptase (NNRTIs), which bind to an allosteric site 10-15 Å from the polymerase active site, play a central role in anti-HIV chemotherapy. Though NNRTIs have been known for 30 years, the pathways by which they bind and unbind from HIV-RT have not been characterized. In crystal structures for complexes, three channels are found to extend from the NNRTI binding site to the exterior of the protein, while added mystery comes from the fact that the binding site is collapsed in the unliganded protein. To address this issue, metadynamics simulations have been performed to elucidate the unbinding of four NNRTIs from HIV-RT. A general and transferable collective variable defined by the distance between the center-of-mass (COM) of the binding pocket and COM of the ligand is used to follow the dynamics while minimizing the bias. The metadynamics also allows computation of the barriers to unbinding, which are compared with the observed potencies of the compounds in an antiviral assay.

Enhanced Monte Carlo Methods for Modeling Proteins Including Computation of Absolute Free Energies of Binding.

The generation of a complete ensemble of geometrical configurations is required to obtain reliable estimations of absolute binding free energies by alchemical free energy methods. Molecular dynamics (MD) is the most popular sampling method, but the representation of large biomolecular systems may be incomplete owing to energetic barriers that impede efficient sampling of the configurational space. Monte Carlo (MC) methods can possibly overcome this issue by adapting the attempted movement sizes to facilitate transitions between alternative local-energy minima. In this study, we present an MC statistical mechanics algorithm to explore the protein-ligand conformational space with emphasis on the motions of the protein backbone and side chains. The parameters for each MC move type were optimized to better reproduce conformational distributions of 18 dipeptides and the well-studied T4-lysozyme L99A protein. Next, the performance of the improved MC algorithms was evaluated by computing absolute free energies of binding for L99A lysozyme with benzene and seven analogs. Results for benzene with L99A lysozyme from MD and the optimized MC protocol were found to agree within 0.6 kcal/mol, while a mean unsigned error of 1.2 kcal/mol between MC results and experiment was obtained for the seven benzene analogs. Significant advantages in computation speed are also reported with MC over MD for similar extents of configurational sampling.

Molecular Dynamics Simulations of a Conformationally Mobile Peptide-Based Catalyst for Atroposelective Bromination.

It is widely accepted that structural rigidity is required to achieve high levels of asymmetric induction in catalytic, enantioselective reactions. This fundamental design principle often does not apply to highly selective catalytic peptides that often exhibit conformational heterogeneity. As a result, these complex systems are particularly challenging to study both experimentally and computationally. Herein, we utilize molecular dynamics simulations to investigate the role of conformational mobility on the reactivity and selectivity exhibited by a catalytic, ß-turn-biased peptide in an atroposelective bromination reaction. By means of cluster analysis, multiple distinct conformers of the peptide and a catalyst-substrate complex were identified in the simulations, all of which were corroborated by experimental NMR measurements. The simulations also revealed that a shift in the conformational equilibrium of the peptidic catalyst occurs upon addition of substrate, and the degree of change varies among different substrates. On the basis of these data, we propose a correlation between the composition of the peptide conformational ensemble and its catalytic properties. Moreover, these findings highlight the importance of conformational dynamics in catalytic, asymmetric reactions mediated by oligopeptides, unveiled through high-level, state-of-the-art computational modeling.

Adding a Hydrogen Bond May Not Help: Naphthyridinone vs Quinoline Inhibitors of Macrophage Migration Inhibitory Factor.

Coordination of the ammonium group of Lys32 in the active site of human macrophage migration inhibitory factor (MIF) using a 1,7-naphthyridin-8-one instead of a quinoline is investigated. Both gas- and aqueous-phase DFT calculations for model systems indicate potential benefits for the added hydrogen bond with the lactam carbonyl group, while FEP results are neutral. Three crystal structures are reported for complexes of MIF with 3a, 4a, and 4b, which show that the desired hydrogen bond is formed with O-N distances of 2.8-3.0 Å. Compound 4b is the most potent new MIF inhibitor with Ki and Kd values of 90 and 94 nM; it also has excellent aqueous solubility, 288 µg/mL. Consistent with the FEP results, the naphthyridinones are found to have similar potency as related quinolines in spite of the additional protein-ligand hydrogen bond.

Improved Description of Sulfur Charge Anisotropy in OPLS Force Fields: Model Development and Parameterization.

The atomic point-charge model used in most molecular mechanics force fields does not represent well the electronic anisotropy that is featured in many directional noncovalent interactions. Sulfur participates in several types of such interactions with its lone pairs and σ-holes. The current study develops a new model, via the addition of off-atom charged sites, for a variety of sulfur compounds in the OPLS-AA and OPLS/CM5 force fields to address the lack of charge anisotropy. Parameter optimization is carried out to reproduce liquid-state properties, torsional and noncovalent energetics from reliable quantum mechanical calculations, and electrostatic potentials. Significant improvements are obtained for computed free energies of hydration, reducing the mean unsigned errors from ca. 1.4 to 0.4-0.7 kcal/mol. Enhanced accuracy in directionality and energetics is also obtained for molecular complexes with sulfur-containing hydrogen and halogen bonds. Moreover, the new model reproduces the unusual conformational preferences of sulfur-containing compounds with 1,4-intramolecular chalcogen bonds. Transferability of the new force field parameters to cysteine and methionine is verified via molecular dynamic simulations of blocked dipeptides. The study demonstrates the effectiveness of using off-atom charge sites to address electronic anisotropy, and provides a parametrization methodology that can be applied to other systems.