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.

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.

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).

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/cm³, 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 (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₅₀ 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.

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/r⁴ 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 10¹-10³ are expected.

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.

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 ₅₀ 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.

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 (ΔG_bind) 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 ΔG_bind 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 ΔG_bind 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 ΔG_bind. These results show that calculation of accurate ΔG_bind 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 S_N2 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).

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.

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.

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.

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.

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 K_i and K_d 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.

1.14*CM1A-LBCC: Localized Bond-Charge Corrected CM1A Charges for Condensed-Phase Simulations

The quality of the 1.14*CM1A and 1.20*CM5 charge models was evaluated for calculations of free energies of hydration. For a set of 426 neutral molecules, 1.14*CM1A and 1.20*CM5 yield MADs of 1.26 and 1.21 kcal/mol, respectively. The 1.14*CM1A charges, which can be readily obtained for large systems, exhibit large deviations only for a subset of functional groups. The results for these cases were systematically improved using localized bond-charge corrections (LBCC) by which offsetting adjustments are made to the partial charges for atoms in specified bond types. Only 19 LBCCs were needed to yield 1.14*CM1A-LBCC charges that reduce the errors for the 426 ΔG_hyd values to only 0.61 kcal/mol. The modified charge method was also tested in computation of heats of vaporization and densities for pure organic liquids, yielding average errors of 1.40 kcal/mol and 0.024 g/cm³, similar to those for 1.14*CM1A.

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.

Performance of Protein-Ligand Force Fields for the Flavodoxin-Flavin Mononucleotide System

The ability to accurately perform molecular dynamics and free energy perturbation calculations for protein-ligand systems is of broad interest to the biophysical and pharmaceutical sciences. In this work, several popular force fields are evaluated for reproducing experimental properties of the flavodoxin/flavin mononucleotide system. Calculated (3)J couplings from molecular dynamics simulations probing φ and χ1 dihedral angles are compared to over 1000 experimental measurements. Free energy perturbation calculations were also executed between different protein mutants for comparison with experimental data for relative free energies of binding. Newer versions of popular protein force fields reproduced (3)J backbone and side chain couplings with good accuracy, with RMSD values near or below one hertz in most cases. OPLS-AA/M paired with CM5 charges for the ligand performed particularly well, both for the (3)J couplings and FEP results, with a mean unsigned error for relative free energies of binding of 0.36 kcal/mol.

Biomolecular Force Field Parameterization via Atoms-in-Molecule Electron Density Partitioning

Molecular mechanics force fields, which are commonly used in biomolecular modeling and computer-aided drug design, typically treat nonbonded interactions using a limited library of empirical parameters that are developed for small molecules. This approach does not account for polarization in larger molecules or proteins, and the parametrization process is labor-intensive. Using linear-scaling density functional theory and atoms-in-molecule electron density partitioning, environment-specific charges and Lennard-Jones parameters are derived directly from quantum mechanical calculations for use in biomolecular modeling of organic and biomolecular systems. The proposed methods significantly reduce the number of empirical parameters needed to construct molecular mechanics force fields, naturally include polarization effects in charge and Lennard-Jones parameters, and scale well to systems comprised of thousands of atoms, including proteins. The feasibility and benefits of this approach are demonstrated by computing free energies of hydration, properties of pure liquids, and the relative binding free energies of indole and benzofuran to the L99A mutant of T4 lysozyme.

Hydration Properties and Solvent Effects for All-Atom Solutes in Polarizable Coarse-Grained Water

Due to the importance of water in chemical and biological systems, a coarse-grained representation of the solvent can greatly simplify the description of the system while retaining key thermodynamic properties of the medium. A multiscale solvation model that couples all-atom solutes and polarizable Martini coarse-grained water (AAX/CGS) is developed to reproduce free energies of hydration of organic solutes. Using Monte Carlo/free energy perturbation (MC/FEP) calculations, results from multiscale and all-atom simulations are compared. Improved accuracy is obtained with the AAX/CGS approach for hydrophobic and sulfur- or halogen-containing solutes, but larger deviations are found for polar solute molecules where hydrogen bonding is featured. Furthermore, solvent effects on conformational and tautomeric equilibria of AA solutes were investigated using AA, CG, and GB/SA solvent models. It is found that the CG solvent model can reproduce well the medium effects from experiment and AA simulations; however, the GB/SA solvent model fails in some cases. A 7-30-fold reduction in computational cost is found for the present AAX/CGS multiscale simulations compared to the AA alternative.

Performance of B3LYP Density Functional Methods for a Large Set of Organic Molecules

Testing of the commonly used hybrid density functional B3LYP with the 6-31G(d), 6-31G(d,p), and 6-31+G(d,p) basis sets has been carried out for 622 neutral, closed-shell organic compounds containing the elements C, H, N, and O. The focus is comparison of computed and experimental heats of formation and isomerization energies. In addition, the effect of an empirical dispersion correction term has been evaluated and found to improve agreement with the experimental data. For the 622 compounds, the mean absolute errors (MAE) in the heats of formation are 3.1, 2.6, 2.7, and 2.4 kcal/mol for B3LYP/6-31G(d), B3LYP/6-31G(d,p), B3LYP/6-31+G(d,p), and B3LYP/6-31+G(d,p) with the dispersion correction. A diverse set of 34 isomerizations highlights specific issues of general interest, such as performance on differences in steric effects, conjugation, and bonding. The corresponding MAEs for the isomerizations are 2.7, 2.4, 2.2, and 1.9 kcal/mol. Improvement is obtained for isomerizations of amines and alcohols when both polarization and diffuse functions are used, but the overstabilization of linear alkanes compared to branched isomers can be relieved only with the dispersion correction. Besides the insights on DFT methods, the study also aimed to quantify the gains in accuracy that can be achieved by replacing energetics from NDO-based semiempirical methods with DFT results. Since the MAEs obtained with the PDDG/PM3 method for the 622 heats of formation and 34 isomerizations are 2.8 and 2.3 kcal/mol, negligible advantage in accuracy for the B3LYP-based methods emerged in the absence of the dispersion corrections.

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.

Improved Peptide and Protein Torsional Energetics with the OPLSAA Force Field

The development and validation of new peptide dihedral parameters are reported for the OPLS-AA force field. High accuracy quantum chemical methods were used to scan φ, ψ, χ1, and χ2 potential energy surfaces for blocked dipeptides. New Fourier coefficients for the dihedral angle terms of the OPLS-AA force field were fit to these surfaces, utilizing a Boltzmann-weighted error function and systematically examining the effects of weighting temperature. To prevent overfitting to the available data, a minimal number of new residue-specific and peptide-specific torsion terms were developed. Extensive experimental solution-phase and quantum chemical gas-phase benchmarks were used to assess the quality of the new parameters, named OPLS-AA/M, demonstrating significant improvement over previous OPLS-AA force fields. A Boltzmann weighting temperature of 2000 K was determined to be optimal for fitting the new Fourier coefficients for dihedral angle parameters. Conclusions are drawn from the results for best practices for developing new torsion parameters for protein force fields.

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.

Evaluation of CM5 Charges for Condensed-Phase Modeling

The recently developed Charge Model 5 (CM5) is tested for its utility in condensed-phase simulations. The CM5 approach, which derives partial atomic charges from Hirshfeld population analyses, provides excellent results for gas-phase dipole moments and is applicable to all elements of the periodic table. Herein, the adequacy of scaled CM5 charges for use in modeling aqueous solutions has been evaluated by computing free energies of hydration (ΔGhyd) for 42 neutral organic molecules via Monte Carlo statistical mechanics. An optimal scaling factor for the CM5 charges was determined to be 1.27, resulting in a mean unsigned error (MUE) of 1.1 kcal/mol for the free energies of hydration. Testing for an additional 20 molecules gave an MUE of 1.3 kcal/mol. The high precision of the results is confirmed by free energy calculations using both sequential perturbations and complete molecular annihilation. Performance for specific functional groups is discussed; sulfur-containing molecules yield the largest errors. In addition, the scaling factor of 1.27 is shown to be appropriate for CM5 charges derived from a variety of density functional methods and basis sets. Though the average errors from the 1.27*CM5 results are only slightly lower than those using 1.14*CM1A charges, the broader applicability and easier access to CM5 charges via the Gaussian program are additional attractive features. The 1.27*CM5 charge model can be used for an enormous variety of applications in conjunction with many fixed-charge force fields and molecular modeling programs.

Enhanced Monte Carlo Sampling through Replica Exchange with Solute Tempering

With a view to improving the consistency of free energy perturbation calculations in Monte Carlo simulations of protein–ligand complexes, we have implemented the replica exchange with solute tempering (REST) method in the MCPRO software. By augmenting the standard REST approach with regular attempted jumps in selected dihedral angles, our combined method facilitates sampling of ligand binding modes that are separated by high free energy barriers and ensures that computed free energy changes are considerably less dependent on the starting conditions and the chosen mutation pathway than those calculated with standard Monte Carlo sampling. We have applied the enhanced sampling method to the calculation of the activities of seven non-nucleoside inhibitors of HIV-1 reverse transcriptase, and its Tyr181Cys variant, and have shown that a range of binding orientations is possible depending on the nature of the ligand and the presence of mutations at the binding site.

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.

Bifunctional inhibition of human immunodeficiency virus type 1 reverse transcriptase: mechanism and proof-of-concept as a novel therapeutic design strategy

Human immunodeficiency virus type 1 reverse transcriptase (HIV-1 RT) is a major target for currently approved anti-HIV drugs. These drugs are divided into two classes: nucleoside and non-nucleoside reverse transcriptase inhibitors (NRTIs and NNRTIs). This study illustrates the synthesis and biochemical evaluation of a novel bifunctional RT inhibitor utilizing d4T (NRTI) and a TMC-derivative (a diarylpyrimidine NNRTI) linked via a poly(ethylene glycol) (PEG) linker. HIV-1 RT successfully incorporates the triphosphate of d4T-4PEG-TMC bifunctional inhibitor in a base-specific manner. Moreover, this inhibitor demonstrates low nanomolar potency that has 4.3-fold and 4300-fold enhancement of polymerization inhibition in vitro relative to the parent TMC-derivative and d4T, respectively. This study serves as a proof-of-concept for the development and optimization of bifunctional RT inhibitors as potent inhibitors of HIV-1 viral replication.

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.

Virtual screening and optimization yield low-nanomolar inhibitors of the tautomerase activity of Plasmodium falciparum macrophage migration inhibitory factor

The Plasmodium falciparum orthologue of the human cytokine, macrophage migratory inhibitory factor (PfMIF), is produced by the parasite during malaria infection and modulates the host's immune response. As for other MIF orthologues, PfMIF has tautomerase activity, whose inhibition may influence the cytokine activity. To identify small-molecule inhibitors of the tautomerase activity of PfMIF, virtual screening has been performed by docking 2.1 million compounds into the enzymatic site. Assaying of 17 compounds identified four as active. Substructure search for the most potent of these compounds, a 4-phenoxypyridine analogue, identified four additional compounds that were purchased and also shown to be active. Thirty-one additional analogues were then designed, synthesized, and assayed. Three were found to be potent PfMIF tautomerase inhibitors with K(i) of ∼40 nM; they are also highly selective with K(i) > 100 μM for human MIF.

Methyl effects on protein-ligand binding

The effects of addition of a methyl group to a lead compound on biological activity are examined. A literature analysis of >2000 cases reveals that an activity boost of a factor of 10 or more is found with an 8% frequency, and a 100-fold boost is a 1 in 200 event. Four cases in the latter category are analyzed in depth to elucidate any unusual aspects of the protein-ligand binding, distribution of water molecules, and changes in conformational energetics. The analyses include Monte Carlo/free-energy perturbation (MC/FEP) calculations for methyl replacements in inhibitor series for p38α MAP kinase, ACK1, PTP1B, and thrombin. Methyl substitutions ortho to an aryl ring can be particularly effective at improving activity by inducing a propitious conformational change. The greatest improvements in activity arise from coupling the conformational gain with the burial of the methyl group in a hydrophobic region of the protein.

Improving MM-GB/SA Scoring through the Application of the Variable Dielectric Model

A variable dielectric model based on residue types for better description of protein-ligand electrostatics in MM-GBSA scoring is reported. The variable dielectric approach provides better correlation with binding data and reduces the score dynamic range, typically observed in the standard MM-GB/SA method. The latter supports the view that exaggerated enthalpic separation between weak and potent compounds due to the lack of shielding effects in the model is greatly responsible for the wide scoring spread.

Effects of Water Placement on Predictions of Binding Affinities for p38α MAP Kinase Inhibitors

Monte Carlo free energy perturbation (MC/FEP) calculations have been applied to compute the relative binding affinities of 17 congeneric pyridazo-pyrimidinone inhibitors of the protein p38α MAP kinase. Overall correlation with experiment was found to be modest when the complexes were hydrated using a traditional procedure with a stored solvent box. Significant improvements in accuracy were obtained when the MC/FEP calculations were repeated using initial solvent distributions optimized by the water placement algorithm JAWS. The results underscore the importance of accurate placement of water molecules in a ligand binding site for the reliable prediction of relative free energies of binding.

Synthesis and evaluation of selected key methyl ether derivatives of vancomycin aglycon

A select series of methyl ether derivatives of vancomcyin aglycon were prepared and examined for antimicrobial activity against vancomycin-sensitive Staphylococcus aureus and vancomycin-resistant Enterococci faecalis as well as their binding affinity for D-Ala-D-Ala and D-Ala-D-Lac. The intent of the study was to elucidate the role selected key methyl groups may play in the improvement of the in vitro antimicrobial profile of the tetra methyl ether derivative of vancomycin aglycon against vancomycin-resistant Enterococci faecalis previously reported. In these studies, methodology for selective derivatization of the A-, B-, and D-ring was developed that defines the relative reactivity of the four phenols of vancomycin aglycon, providing a foundation for future efforts for site-directed modification of the vancomycin aglycon core.

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.

Energetics of displacing water molecules from protein binding sites: consequences for ligand optimization

A strategy in drug design is to consider enhancing the affinity of lead molecules with structural modifications that displace water molecules from a protein binding site. Because success of the approach is uncertain, clarification of the associated energetics was sought in cases where similar structural modifications yield qualitatively different outcomes. Specifically, free-energy perturbation calculations were carried out in the context of Monte Carlo statistical mechanics simulations to investigate ligand series that feature displacement of ordered water molecules in the binding sites of scytalone dehydratase, p38-alphaMAP kinase, and EGFR kinase. The change in affinity for a ligand modification is found to correlate with the ease of displacement of the ordered water molecule. However, as in the EGFR example, the binding affinity may diminish if the free-energy increase due to the removal of the bound water molecule is not more than compensated by the additional interactions of the water-displacing moiety. For accurate computation of the effects of ligand modifications, a complete thermodynamic analysis is shown to be needed. It requires identification of the location of water molecules in the protein-ligand interface and evaluation of the free-energy changes associated with their removal and with the introduction of the ligand modification. Direct modification of the ligand in free-energy calculations is likely to trap the ordered molecule and provide misleading guidance for lead optimization.

Prediction of the water content in protein binding sites

An efficient molecular simulation methodology has been developed to determine the positioning of water molecules in the binding site of a protein or protein-ligand complex. Occupancies and absolute binding free energies of water molecules are computed using a statistical thermodynamics approach. The methodology, referred to as Just Add Water Molecules (JAWS), features "theta-water" molecules that can appear and disappear on a binding-site grid. Key approximations render the technique far more efficient than conventional free energy simulations. Testing of JAWS on five diverse examples (neuraminidase, scytalone dehydratase, major urinary protein 1, beta-lactoglobulin, and COX-2) demonstrates its accuracy in locating hydration sites in comparison to results from high-resolution crystal structures. Possible applications include aid in refinement of protein crystal structures, drug lead optimization, setup of docking calculations, and simulations of protein-ligand complexes.

Vancomycin analogs: Seeking improved binding of d-Ala-d-Ala and d-Ala-d-Lac peptides by side-chain and backbone modifications

In order to seek vancomycin analogs with improved performance against VanA and VanB resistant bacterial strains, extensive computational investigations have been performed to examine the effects of side-chain and backbone modifications. Changes in binding affinities for tripeptide cell-wall precursor mimics, Ac(2)-l-Lys-d-Ala-d-Ala (3) and Ac(2)-l-Lys-d-Ala-d-Lac (4), with vancomycin analogs were computed with Monte Carlo/free energy perturbation (MC/FEP) calculations. Replacements of the 3-hydroxyl group in residue 7 with small alkyl or alkoxy groups, which improve contacts with the methyl side chain of the ligands'd-Ala residue, are predicted to be the most promising to enhance binding for both ligands. The previously reported amine backbone modification as in 5 is shown to complement the hydrophobic modifications for binding monoacetylated tripeptides. In addition, replacement of the hydroxyl groups in residues 5 and 7 by fluorine is computed to have negligible impact on binding the tripeptides, though it may be pharmacologically advantageous.

In Silico Improvement of beta3-peptide inhibitors of p53 x hDM2 and p53 x hDMX

There is great interest in molecules capable of inhibiting the interactions between p53 and its negative regulators hDM2 and hDMX, as these molecules have validated potential against cancers in which one or both oncoproteins are overexpressed. We reported previously that appropriately substituted beta(3)-peptides inhibit these interactions and, more recently, that minimally cationic beta(3)-peptides are sufficiently cell permeable to upregulate p53-dependent genes in live cells. These observations, coupled with the known stability of beta-peptides in a cellular environment, and the recently reported structures of hDM2 and hDMX, motivated us to exploit computational modeling to identify beta-peptides with improved potency and/or selectivity. This exercise successfully identified a new beta(3)-peptide, beta53-16, that possesses the highly desirable attribute of high affinity for both hDM2 and hDMX and identifies the 3,4-dichlorophenyl moiety as a novel determinant of hDMX affinity.

Discovery of wild-type and Y181C mutant non-nucleoside HIV-1 reverse transcriptase inhibitors using virtual screening with multiple protein structures

To discover non-nucleoside inhibitors of HIV-1 reverse transcriptase (NNRTIs) that are effective against both wild-type (WT) virus and variants that encode the clinically troublesome Tyr181Cys (Y181C) RT mutation, virtual screening by docking was carried out using three RT structures and more than 2 million commercially available compounds. Two of the structures are for WT-virus with different conformations of Tyr181, while the third structure incorporates the Y181C modification. Eventually nine compounds were purchased and assayed. Three of the compounds show low-micromolar antiviral activity toward either or both the wild-type and Y181C HIV-1 strains. The study illustrates a viable protocol to seek anti-HIV agents with enhanced resistance profiles.

Exploiting structural analysis, in silico screening, and serendipity to identify novel inhibitors of drug-resistant falciparum malaria

Plasmodium falciparum thymidylate synthase-dihydrofolate reductase (TS-DHFR) is an essential enzyme in folate biosynthesis and a major malarial drug target. This bifunctional enzyme thus presents different design approaches for developing novel inhibitors against drug-resistant mutants. We performed a high-throughput in silico screen of a database of diverse, drug-like molecules against a non-active-site pocket of TS-DHFR. The top compounds from this virtual screen were evaluated by in vitro enzymatic and cellular culture studies. Three compounds active to 20 microM IC(50)'s in both wildtype and antifolate-resistant P. falciparum parasites were identified; moreover, no inhibition of human DHFR enzyme was observed, indicating that the inhibitory effects appeared to be parasite-specific. Notably, all three compounds had a biguanide scaffold. However, relative free energy of binding calculations suggested that the compounds might preferentially interact with the active site over the screened non-active-site region. To resolve the two possible modes of binding, co-crystallization studies of the compounds complexed with TS-DHFR enzyme were performed. Surprisingly, the structural analysis revealed that these novel, biguanide compounds do indeed bind at the active site of DHFR and additionally revealed the molecular basis by which they overcome drug resistance. To our knowledge, these are the first co-crystal structures of novel, biguanide, non-WR99210 compounds that are active against folate-resistant malaria parasites in cell culture.

Vancomycin resistance: modeling backbone variants with D-Ala-D-Ala and D-Ala-D-Lac peptides

To seek vancomycin analogs with broader antibacterial activity, effects of backbone modifications for the agylcon 2 on binding with D-Ala-D-Ala- and D-Ala-D-Lac-containing peptides were investigated by Monte Carlo/free energy perturbation (MC/FEP) calculations. The experimental trend in binding affinities for 2 with three tripeptides was well reproduced. Possible modifications of the peptide bond between residues 4 and 5 were then considered, specifically for conversion of the OCNH linkage to CH(2)NH(2)(+) (6), FCCH (7), HCCH (8), and HNCO (9). The MC/FEP results did not yield binding improvements for 7, 8, and 9, though the fluorovinyl replacement is relatively benign. The previously reported analog 6 remains as the only variant that exhibits improved affinity for the D-Ala-D-Lac sequence and acceptable affinity for the D-Ala-D-Ala sequence.