Proof-of-concept studies with a computationally designed Mpro inhibitor as a synergistic combination regimen alternative to Paxlovid

As the SARS-CoV-2 virus continues to spread and mutate, it remains important to focus not only on preventing spread through vaccination but also on treating infection with direct-acting antivirals (DAA). The approval of Paxlovid, a SARS-CoV-2 main protease (Mpro) DAA, has been significant for treatment of patients. A limitation of this DAA, however, is that the antiviral component, nirmatrelvir, is rapidly metabolized and requires inclusion of a CYP450 3A4 metabolic inhibitor, ritonavir, to boost levels of the active drug. Serious drug-drug interactions can occur with Paxlovid for patients who are also taking other medications metabolized by CYP4503A4, particularly transplant or otherwise immunocompromised patients who are most at risk for SARS-CoV-2 infection and the development of severe symptoms. Developing an alternative antiviral with improved pharmacological properties is critical for treatment of these patients. By using a computational and structure-guided approach, we were able to optimize a 100 to 250 μM screening hit to a potent nanomolar inhibitor and lead compound, Mpro61. In this study, we further evaluate Mpro61 as a lead compound, starting with examination of its mode of binding to SARS-CoV-2 Mpro. In vitro pharmacological profiling established a lack of off-target effects, particularly CYP450 3A4 inhibition, as well as potential for synergy with the currently approved alternate antiviral, molnupiravir. Development and subsequent testing of a capsule formulation for oral dosing of Mpro61 in B6-K18-hACE2 mice demonstrated favorable pharmacological properties, efficacy, and synergy with molnupiravir, and complete recovery from subsequent challenge by SARS-CoV-2, establishing Mpro61 as a promising potential preclinical candidate.

Enthalpies and entropies of hydration from Monte Carlo simulations

The changes in free energy, enthalpy, and entropy for transfer of a solute from the gas phase into solution are the fundamental thermodynamic quantities that characterize the solvation process. Owing to the development of methods based on free-energy perturbation theory, computation of free energies of solvation has become routine in conjunction with Monte Carlo (MC) statistical mechanics and molecular dynamics (MD) simulations. Computation of the enthalpy change and by inference the entropy change is more challenging. Two methods are considered in this work corresponding to direct averaging for the solvent and solution and to computing the temperature derivative of the free energy in the van't Hoff approach. The application is for neutral organic solutes in TIP4P water using long MC simulations to improve precision. Definitive results are also provided for pure TIP4P water. While the uncertainty in computed free energies of hydration is ca. 0.05 kcal mol-1, it is ca. 0.4 kcal mol-1 for the enthalpy changes from either van't Hoff plots or the direct method with sampling for 5 billion MC configurations. Partial molar volumes of hydration are also computed by the direct method; they agree well with experimental data with an average deviation of 3 cm3 mol-1. In addition, the results permit breakdown of the errors in the free energy changes from the OPLS-AA force field into their enthalpic and entropic components. The excess hydrophobicity of organic solutes is enthalpic in origin.

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.

Covalent and noncovalent strategies for targeting Lys102 in HIV-1 reverse transcriptase

Reverse transcriptase (RT) is one of three key proteins responsible for the replication cycle of HIV-1 in the host. Several classes of inhibitors have been developed to target the enzyme, with non-nucleoside reverse transcriptase inhibitors forming first-line treatment. Previously, covalent RT inhibitors have been identified and found to bind irreversibly to commonly mutated residues such as Y181C. In this work we aim to circumvent the issue of NNRTI resistance through targeting K102, which has not yet been identified to confer drug resistance. As reported here, 34 compounds were synthesized and characterized biochemically and structurally with wild-type (WT) HIV-1 RT. Two of these inhibitors demonstrate covalent inhibition as evidenced by protein crystallography, enzyme kinetics, mass spectrometry, and antiviral potency in HIV-1 infected human T-cell assays.

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.

Exploring novel HIV-1 reverse transcriptase inhibitors with drug-resistant mutants: A double mutant surprise

HIV-1 reverse transcriptase (RT) remains a key target for HIV drug development. As successful management of the disease requires lifelong treatment, the emergence of resistance mutations is inevitable, making development of new RT inhibitors, which remain effective against resistant variants crucial. To this end, previous computationally guided drug design efforts have resulted in catechol diether compounds, which inhibit wildtype RT with picomolar affinities and appear to be promising preclinical candidates. To confirm that these compounds remain potent against Y181C, a widespread mutation conferring resistance to first generation inhibitors, they were screened against the HIV-1 N119 clinical isolate, reported as a Y181C single mutant. In comparison to a molecular clone with the same mutation, N119 appears less susceptible to inhibition by our preclinical candidate compounds. A more detailed sequencing effort determined that N119 was misidentified and carries V106A in combination with Y181C. While both indolizine and naphthalene substituted catechol diethers are potent against the classical Y181C single mutant, the addition of V106A confers more resistance against the indolizine derivatives than the naphthalene derivatives. Crystal structures presented in this study highlight key features of the naphthyl group, which allow these compounds to remain potent in the double mutant, including stronger interactions with F227 and less reliance on V106 for stabilization of the ethoxy-uracil ring, which makes critical hydrogen bonds with other residues in the binding pocket.

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

A High-Throughput, High-Containment Human Primary Epithelial Airway Organ-on-Chip Platform for SARS-CoV-2 Therapeutic Screening

COVID-19 emerged as a worldwide pandemic in early 2020, and while the rapid development of safe and efficacious vaccines stands as an extraordinary achievement, the identification of effective therapeutics has been less successful. This process has been limited in part by a lack of human-relevant preclinical models compatible with therapeutic screening on the native virus, which requires a high-containment environment. Here, we report SARS-CoV-2 infection and robust viral replication in PREDICT96-ALI, a high-throughput, human primary cell-based organ-on-chip platform. We evaluate unique infection kinetic profiles across lung tissue from three human donors by immunofluorescence, RT-qPCR, and plaque assays over a 6-day infection period. Enabled by the 96 devices/plate throughput of PREDICT96-ALI, we also investigate the efficacy of Remdesivir and MPro61 in a proof-of-concept antiviral study. Both compounds exhibit an antiviral effect against SARS-CoV-2 in the platform. This demonstration of SARS-CoV-2 infection and antiviral dosing in a high-throughput organ-on-chip platform presents a critical capability for disease modeling and therapeutic screening applications in a human physiology-relevant in vitro system.

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.

Progress on the Pharmacological Targeting of Janus Pseudokinases

The Janus kinases (JAKs) are key components of the JAK-STAT signaling pathway and are involved in myriad physiological processes. Though they are the molecular targets of many FDA-approved drugs, these drugs manifest adverse effects due in part to their inhibition of the requisite JAK kinase activity. However, the JAKs uniquely possess an integrated pseudokinase domain (JH2) that regulates the adjacent kinase domain (JH1). The therapeutic targeting of JH2 domains has been less thoroughly explored and may present an avenue to modulate the JAKs without the adverse effects associated with targeting the adjacent JH1 domain. The potential of this strategy was recently demonstrated with the FDA approval of the TYK2 JH2 ligand deucravacitinib for treating plaque psoriasis. In this light, the structure and targetability of the JAK pseudokinases are discussed, in conjunction with the state of development of ligands that bind to these domains.

Design, synthesis, and biological testing of biphenylmethyloxazole inhibitors targeting HIV-1 reverse transcriptase

We report non-nucleoside inhibitors of HIV-1 reverse transcriptase (NNRTIs) using a biphenylmethyloxazole pharmacophore. A crystal structure of benzyloxazole 1 was obtained and suggested the potential viability of biphenyl analogues. In particular, 6a, 6b, and 7 turned out to be potent NNRTIs with low-nanomolar activity in enzyme inhibition and infected T-cell assays, and with low cytotoxicity. Though modeling further suggested that analogues with fluorosulfate and epoxide warheads might provide covalent modification of Tyr188, synthesis and testing did not find evidence for this outcome.

Covalent Modification of the JH2 Domain of Janus Kinase 2

Probe molecules that covalently modify the JAK2 pseudokinase domain (JH2) are reported. Selective targeting of JH2 domains over the kinase (JH1) domains is a necessary feature for ligands intended to evaluate JH2 domains as therapeutic targets. The JH2 domains of three Janus kinases (JAK1, JAK2, and TYK2) possess a cysteine residue in the catalytic loop that does not occur in their JH1 domains. Starting from a non-selective kinase binding molecule, computer-aided design directed attachment of substituents terminating in acrylamide warheads to modify Cys675 of JAK2 JH2. Successful covalent attachment was demonstrated first through observation of enhanced binding with increasing incubation time in fluorescence polarization experiments. Covalent binding also increased selectivity to as much as ca. 30-fold for binding the JAK2 JH2 domain over the JH1 domain after a 20-h incubation. Covalency was confirmed through HPLC electrospray quadrupole time-of-flight HRMS experiments, which revealed the expected mass shifts.

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.

Insights on JAK2 Modulation by Potent, Selective, and Cell-Permeable Pseudokinase-Domain Ligands

JAK2 is a non-receptor tyrosine kinase that regulates hematopoiesis through the JAK-STAT pathway. The pseudokinase domain (JH2) is an important regulator of the activity of the kinase domain (JH1). V617F mutation in JH2 has been associated with the pathogenesis of various myeloproliferative neoplasms, but JAK2 JH2 has been poorly explored as a pharmacological target. In light of this, we aimed to develop JAK2 JH2 binders that could selectively target JH2 over JH1 and test their capacity to modulate JAK2 activity in cells. Toward this goal, we optimized a diaminotriazole lead compound into potent, selective, and cell-permeable JH2 binders leveraging computational design, synthesis, binding affinity measurements for the JH1, JH2 WT, and JH2 V617F domains, permeability measurements, crystallography, and cell assays. Optimized diaminotriazoles are capable of inhibiting STAT5 phosphorylation in both WT and V617F JAK2 in cells.

Conversion of a False Virtual Screen Hit into Selective JAK2 JH2 Domain Binders Using Convergent Design Strategies

The Janus kinase 2 (JAK2) pseudokinase domain (JH2) is an ATP-binding domain that regulates the activity of the catalytic tyrosine kinase domain (JH1). Dysregulation of JAK2 JH1 signaling caused by the V617F mutation in JH2 is implicated in various myeloproliferative neoplasms. To explore if JAK2 activity can be modulated by a small molecule binding to the ATP site in JH2, we have developed several ligand series aimed at selectively targeting the JAK2 JH2 domain. We report here the evolution of a false virtual screen hit into a new JAK2 JH2 series. Optimization guided by computational modeling has yielded analogues with nanomolar affinity for the JAK2 JH2 domain and >100-fold selectivity for the JH2 domain over the JH1 domain. A crystal structure for one of the potent compounds bound to JAK2 JH2 clarifies the origins of the strong binding and selectivity. The compounds expand the platform for seeking molecules to regulate JAK2 signaling, including V617F JAK2 hyperactivation.

Structural Studies and Structure Activity Relationships for Novel Computationally Designed Non-nucleoside Inhibitors and Their Interactions With HIV-1 Reverse Transcriptase

Reverse transcriptase (RT) from the human immunodeficiency virus continues to be an attractive drug target for antiretroviral therapy. June 2022 will commemorate the 30th anniversary of the first Human Immunodeficiency Virus (HIV) RT crystal structure complex that was solved with non-nucleoside reverse transcriptase inhibitor nevirapine. The release of this structure opened opportunities for designing many families of non-nucleoside reverse transcriptase inhibitors (NNRTIs). In paying tribute to the first RT-nevirapine structure, we have developed several compound classes targeting the non-nucleoside inhibitor binding pocket of HIV RT. Extensive analysis of crystal structures of RT in complex with the compounds informed iterations of structure-based drug design. Structures of seven additional complexes were determined and analyzed to summarize key interactions with residues in the non-nucleoside inhibitor binding pocket (NNIBP) of RT. Additional insights comparing structures with antiviral data and results from molecular dynamics simulations elucidate key interactions and dynamics between the nucleotide and non-nucleoside binding sites.

Long-acting and extended-release implant and nanoformulations with a synergistic antiretroviral two-drug combination controls HIV-1 infection in a humanized mouse model

The HIV pandemic has affected over 38 million people worldwide with close to 26 million currently accessing antiretroviral therapy (ART). A major challenge in the long-term treatment of HIV-1 infection is nonadherence to ART. Long-acting antiretroviral (LA-ARV) formulations, that reduce dosing frequency to less than once a day, are an urgent need that could tackle the adherence issue. Here, we have developed two LA-ART interventions, one an injectable nanoformulation, and the other, a removable implant, for the delivery of a synergistic two-drug ARV combination comprising a pre-clinical nonnucleoside reverse transcriptase inhibitor (NNRTI), Compound I, and the nucleoside reverse transcriptase inhibitor (NRTI), 4'-ethynyl-2-fluoro-2'-deoxyadenosine. The nanoformulation is poly(lactide-co-glycolide)-based and the implant is a copolymer of ω-pentadecalactone and p-dioxanone, poly(PDL-co-DO), a novel class of biocompatible, biodegradable materials. Both the interventions, packaged independently with each ARV, released sustained levels of the drugs, maintaining plasma therapeutic indices for over a month, and suppressed viremia in HIV-1-infected humanized mice for up to 42 days with maintenance of CD4+ T cells. These data suggest promise in the use of these new drugs as LA-ART formulations in subdermal implant and injectable mode.

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.

Structure-guided design of a perampanel-derived pharmacophore targeting the SARS-CoV-2 main protease

There is a clinical need for direct-acting antivirals targeting SARS-CoV-2, the coronavirus responsible for the COVID-19 pandemic, to complement current therapeutic strategies. The main protease (Mpro) is an attractive target for antiviral therapy. However, the vast majority of protease inhibitors described thus far are peptidomimetic and bind to the active-site cysteine via a covalent adduct, which is generally pharmacokinetically unfavorable. We have reported the optimization of an existing FDA-approved chemical scaffold, perampanel, to bind to and inhibit Mpro noncovalently with IC50s in the low-nanomolar range and EC50s in the low-micromolar range. Here, we present nine crystal structures of Mpro bound to a series of perampanel analogs, providing detailed structural insights into their mechanism of action and structure-activity relationship. These insights further reveal strategies for pursuing rational inhibitor design efforts in the context of considerable active-site flexibility and potential resistance mechanisms.

Indoloxytriazines as binding molecules for the JAK2 JH2 pseudokinase domain and its V617F variant

Small molecules that selectively bind to the pseudokinase JH2 domain over the JH1 kinase domain of JAK2 kinase are sought. Virtual screening led to the purchase of 17 compounds among which 9 were found to bind to V617F JAK2 JH2 with affinities of 40 - 300 μM in a fluorogenic assay. Ten analogues were then purchased yielding 9 additional active compounds. Aminoanilinyltriazine 22 was particularly notable as it shows no detectable binding to JAK2 JH1, and it has a 65-μM dissociation constant K _d with V617F JAK2 JH2. A crystal structure for 22 in complex with wild-type JAK2 JH2 was obtained to elucidate the binding mode. Additional de novo design led to the synthesis of 19 analogues of 22 with the most potent being 33n with K _d values of 2-3 μM for WT and V617F JAK2 JH2, and with 16-fold selectivity relative to binding with WT JAK2 JH1.

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.

Covalent Inhibition of Wild-Type HIV-1 Reverse Transcriptase Using a Fluorosulfate Warhead

Covalent inhibitors of wild-type HIV-1 reverse transcriptase (CRTIs) are reported. Three compounds derived from catechol diether non-nucleoside inhibitors (NNRTIs) with addition of a fluorosulfate warhead are demonstrated to covalently modify Tyr181 of HIV-RT. X-ray crystal structures for complexes of the CRTIs with the enzyme are provided, which fully demonstrate the covalent attachment, and confirmation is provided by appropriate mass shifts in ESI-TOF mass spectra. The three CRTIs and six noncovalent analogues are found to be potent inhibitors with both IC₅₀ values for in vitro inhibition of WT RT and EC₅₀ values for cytopathic protection of HIV-1-infected human T-cells in the 5-320 nM range.

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.

Targeting the TS dimer interface in bifunctional Cryptosporidium hominis TS-DHFR from parasitic protozoa: Virtual screening identifies novel TS allosteric inhibitors

Effective therapies are lacking to treat gastrointestinal infections caused by the genus Cryptosporidium, which can be fatal in the immunocompromised. One target of interest is Cryptosporidium hominis (C. hominis) thymidylate synthase-dihydrofolate reductase (ChTS-DHFR), a bifunctional enzyme necessary for DNA biosynthesis. Targeting the TS-TS dimer interface is a novel strategy previously used to identify inhibitors against the related bifunctional enzyme in Toxoplasma gondii. In the present study, we target the ChTS dimer interface through homology modeling and high-throughput virtual screening to identifying allosteric, ChTS-specific inhibitors. Our work led to the discovery of methylenedioxyphenyl-aminophenoxypropanol analogues which inhibit ChTS activity in a manner that is both dose-dependent and influenced by the conformation of the enzyme. Preliminary results presented here include an analysis of structure activity relationships and a ChTS-apo crystal structure of ChTS-DHFR supporting the continued development of inhibitors that stabilize a novel pocket formed in the open conformation of ChTS-TS.

Structural investigation of 2-naphthyl phenyl ether inhibitors bound to WT and Y181C reverse transcriptase highlights key features of the NNRTI binding site

Human immunodeficiency virus (HIV)-1 remains as a global health issue that is primarily treated with highly active antiretroviral therapy, a combination of drugs that target the viral life cycle. One class of these drugs are non-nucleoside reverse transcriptase inhibitors (NNRTIs) that target the viral reverse transcriptase (RT). First generation NNRTIs were troubled with poor pharmacological properties and drug resistance, incentivizing the development of improved compounds. One class of developed compounds are the 2-naphthyl phenyl ethers, showing promising efficacy against the Y181C RT mutation. Further biochemical and structural work demonstrated differences in potency against the Y181C mutation and binding mode of the compounds. This work aims to understand the relationship between the binding mode and ability to overcome drug resistance using macromolecular x-ray crystallography. Comparison of 2-naphthyl phenyl ethers bound to Y181C RT reveal that compounds that interact with the invariant W229 are more capable of retaining efficacy against the resistance mutation. Additional modifications to these compounds at the 4-position, computationally designed to compensate for the Y181C mutation, do not demonstrate improved potency. Ultimately, we highlight important considerations for the development of future HIV-1 drugs that are able to combat drug resistance.

Structure-Guided Identification of DNMT3B Inhibitors

Methyltransferase 3 beta (DNMT3B) inhibitors that interfere with cancer growth are emerging possibilities for treatment of melanoma. Herein we identify small molecule inhibitors of DNMT3B starting from a homology model based on a DNMT3A crystal structure. Virtual screening by docking led to purchase of 15 compounds, among which 5 were found to inhibit the activity of DNMT3B with IC₅₀ values of 13-72 μM in a fluorogenic assay. Eight analogues of 7, 10, and 12 were purchased to provide 2 more active compounds. Compound 11 is particularly notable as it shows good selectivity with no inhibition of DNMT1 and 22 μM potency toward DNMT3B. Following additional de novo design, exploratory synthesis of 17 analogues of 11 delivered 5 additional inhibitors of DNMT3B with the most potent being 33h with an IC₅₀ of 8.0 μM. This result was well confirmed in an ultrahigh-performance liquid chromatography (UHPLC)-based analytical assay, which yielded an IC₅₀ of 4.8 μM. Structure-activity data are rationalized based on computed structures for DNMT3B complexes.

Selective Janus Kinase 2 (JAK2) Pseudokinase Ligands with a Diaminotriazole Core

Janus kinases (JAKs) are non-receptor tyrosine kinases that are essential components of the JAK-STAT signaling pathway. Associated aberrant signaling is responsible for many forms of cancer and disorders of the immune system. The present focus is on the discovery of molecules that may regulate the activity of JAK2 by selective binding to the JAK2 pseudokinase domain, JH2. Specifically, the Val617Phe mutation in JH2 stimulates the activity of the adjacent kinase domain (JH1) resulting in myeloproliferative disorders. Starting from a non-selective screening hit, we have achieved the goal of discovering molecules that preferentially bind to the ATP binding site in JH2 instead of JH1. We report the design and synthesis of the compounds and binding results for the JH1, JH2, and JH2 V617F domains, as well as five crystal structures for JH2 complexes. Testing with a selective and non-selective JH2 binder on the autophosphorylation of wild-type and V617F JAK2 is also contrasted.

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.

Molecular and cellular studies evaluating a potent 2-cyanoindolizine catechol diether NNRTI targeting wildtype and Y181C mutant HIV-1 reverse transcriptase

The development of efficacious NNRTIs for HIV/AIDS therapy is commonly met with the emergence of drug resistant strains, including the Y181C variant. Using a computationally-guided approach, we synthesized the catechol diether series of NNRTIs, which display sub-nanomolar potency in cellular assays. Among the most potent were a series of 2-cyanoindolizine substituted catechol diethers, including Compound 1. We present here a thorough evaluation of this compound, including biochemical, cellular, and structural studies. The compound demonstrates low nanomolar potency against both WT and Y181C HIV-1 RT in in vitro and cellular assays. Our crystal structures of both the wildtype and mutant forms of RT in complex with Compound 1 allow the interrogation of this compound's features that allow it to maintain strong efficacy against the drug resistant mutant. Among these are compensatory shifts in the NNRTI binding pocket, persistence of multiple hydrogen bonds, and van der Waals contacts throughout the binding site. Further, the fluorine at the C6 position of the indolizine moiety makes multiple favorable interactions with both RT forms. The present study highlights the indolizine-substituted catechol diether class of NNRTIs as promising therapeutic candidates possessing optimal pharmacological properties and significant potency against multiple RT variants.

Structural and pharmacological evaluation of a novel non-nucleoside reverse transcriptase inhibitor as a promising long acting nanoformulation for treating HIV

Combination antiretroviral therapy (cART) has been proven effective in inhibiting human immunodeficiency virus type 1 (HIV-1) infection and has significantly improved the health outcomes in acquired immune deficiency syndrome (AIDS) patients. The therapeutic benefits of cART have been challenged because of the toxicity and emergence of drug-resistant HIV-1 strains along with lifelong patient compliance resulting in non-adherence. These issues also hinder the clinical benefits of non-nucleoside reverse transcriptase inhibitors (NNRTIs), which are one of the vital components of cART for the treatment of HIV-1 infection. In this study, using a computational and structural based drug design approach, we have discovered an effective HIV -1 NNRTI, compound I (Cmpd I) that is very potent in biochemical assays and which targets key residues in the allosteric binding pocket of wild-type (WT)-RT as revealed by structural studies. Furthermore, Cmpd I exhibited very potent antiviral activity in HIV-1 infected T cells, lacked cytotoxicity (therapeutic index >100,000), and no significant off-target effects were noted in pharmacological assays. To address the issue of non-adherence, we developed a long-acting nanoformulation of Cmpd I (Cmpd I-NP) using poly (lactide-coglycolide) (PLGA) particles. The pharmacokinetic studies of free and nanoformulated Cmpd I were carried out in BALB/c mice. Intraperitoneal administration of Cmpd I and Cmpd I-NP in BALB/c mice revealed prolonged serum residence time of 48 h and 30 days, respectively. The observed serum concentrations of Cmpd I in both cases were sufficient to provide >97% inhibition in HIV-1 infected T-cells. The significant antiviral activity along with favorable pharmacological and pharmacokinetic profile of Cmpd I, provide compelling and critical support for its further development as an anti-HIV therapeutic agent.

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.

Understanding the structural basis of species selective, stereospecific inhibition for Cryptosporidium and human thymidylate synthase

Thymidylate synthase (TS), found in all organisms, is an essential enzyme responsible for the de novo synthesis of deoxythymidine monophosphate. The TS active sites of the protozoal parasite Cryptosporidium hominis and human are relatively conserved. Evaluation of antifolate compound 1 and its R-enantiomer 2 against both enzymes reveals divergent inhibitor selectivity and enzyme stereospecificity. To establish how C. hominis and human TS (ChTS and hTS) selectively discriminate 1 and 2, respectively, we determined crystal structures of ChTS complexed with 2 and hTS complexed with 1 or 2. Coupled with the previously determined structure of ChTS complexed with 1, we discuss a possible mechanism for enzyme stereospecificity and inhibitor selectivity.

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.

Novel allosteric covalent inhibitors of bifunctional Cryptosporidium hominis TS-DHFR from parasitic protozoa identified by virtual screening

Protozoans of the genus Cryptosporidium are the causative agent of the gastrointestinal disease, cryptosporidiosis, which can be fatal in immunocompromised individuals. Cryptosporidium hominis (C. hominis) bifunctional thymidylate synthase-dihydrofolate reductase (TS-DHFR) is an essential enzyme in the folate biosynthesis pathway and a molecular target for inhibitor design. Previous studies have demonstrated the importance of the ChTS-DHFR linker region "crossover helix" to the enzymatic activity and stability of the ChDHFR domain. We conducted a virtual screen of a novel non-active site pocket located at the interface of the ChDHFR domain and crossover helix. From this screen we have identified and characterized a noncompetitive inhibitor, compound 15, a substituted diphenyl thiourea. Through subsequent structure activity relationship studies, we have identified a time-dependent inhibitor lead, compound 15D17, a thiol-substituted 2-hydroxy-N-phenylbenzamide, which is selective for ChTS-DHFR, and whose effects appear to be mediated by covalent bond formation with a non-catalytic cysteine residue adjacent to the non-active site pocket.

Computation of protein-ligand binding free energies using quantum mechanical bespoke force fields

A quantum mechanical bespoke molecular mechanics force field is derived for the L99A mutant of T4 lysozyme and used to compute absolute binding free energies of six benzene analogs to the protein. Promising agreement between theory and experiment highlights the potential for future use of system-specific force fields in computer-aided drug design.

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.

Optimization of Pyrazoles as Phenol Surrogates to Yield Potent Inhibitors of Macrophage Migration Inhibitory Factor

Macrophage migration inhibitory factor (MIF) is a proinflammatory cytokine that is implicated in the regulation of inflammation, cell proliferation, and neurological disorders. MIF is also an enzyme that functions as a keto-enol tautomerase. Most potent MIF tautomerase inhibitors incorporate a phenol, which hydrogen bonds to Asn97 in the active site. Starting from a 113-μm docking hit, we report results of structure-based and computer-aided design that have provided substituted pyrazoles as phenol alternatives with potencies of 60-70 nm. Crystal structures of complexes of MIF with the pyrazoles highlight the contributions of hydrogen bonding with Lys32 and Asn97, and aryl-aryl interactions with Tyr36, Tyr95, and Phe113 to the binding.

Advances and Insights for Small Molecule Inhibition of Macrophage Migration Inhibitory Factor

Macrophage migration inhibitory factor (MIF) is an upstream regulator of the immune response whose dysregulation is tied to a broad spectrum of inflammatory and proliferative disorders. As its complex signaling pathways and pleiotropic nature have been elucidated, it has become an attractive target for drug discovery. Remarkably, MIF is both a cytokine and an enzyme that functions as a keto-enol tautomerase. Strategies including in silico modeling, virtual screening, high-throughput screening, and screening of anti-inflammatory natural products have led to a large and diverse catalogue of MIF inhibitors as well as some understanding of the structure-activity relationships for compounds binding MIF's tautomerase active site. With possible clinical trials of some MIF inhibitors on the horizon, it is an opportune time to review the literature to seek trends, address inconsistencies, and identify promising new avenues of research.

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.

Computationally-guided optimization of small-molecule inhibitors of the Aurora A kinase-TPX2 protein-protein interaction

Free energy perturbation theory, in combination with enhanced sampling of protein-ligand binding modes, is evaluated in the context of fragment-based drug design, and used to design two new small-molecule inhibitors of the Aurora A kinase-TPX2 protein-protein interaction.

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.