Open Force Field BespokeFit: Automating Bespoke Torsion Parametrization at Scale.

The development of accurate transferable force fields is key to realizing the full potential of atomistic modeling in the study of biological processes such as protein-ligand binding for drug discovery. State-of-the-art transferable force fields, such as those produced by the Open Force Field Initiative, use modern software engineering and automation techniques to yield accuracy improvements. However, force field torsion parameters, which must account for many stereoelectronic and steric effects, are considered to be less transferable than other force field parameters and are therefore often targets for bespoke parametrization. Here, we present the Open Force Field QCSubmit and BespokeFit software packages that, when combined, facilitate the fitting of torsion parameters to quantum mechanical reference data at scale. We demonstrate the use of QCSubmit for simplifying the process of creating and archiving large numbers of quantum chemical calculations, by generating a dataset of 671 torsion scans for druglike fragments. We use BespokeFit to derive individual torsion parameters for each of these molecules, thereby reducing the root-mean-square error in the potential energy surface from 1.1 kcal/mol, using the original transferable force field, to 0.4 kcal/mol using the bespoke version. Furthermore, we employ the bespoke force fields to compute the relative binding free energies of a congeneric series of inhibitors of the TYK2 protein, and demonstrate further improvements in accuracy, compared to the base force field (MUE reduced from 0.560.390.77 to 0.420.280.59 kcal/mol and R2 correlation improved from 0.720.350.87 to 0.930.840.97).

Kinetics and free energy of ligand dissociation using weighted ensemble milestoning.

We consider the recently developed weighted ensemble milestoning (WEM) scheme [D. Ray and I. Andricioaei, J. Chem. Phys. 152, 234114 (2020)] and test its capability of simulating ligand-receptor dissociation dynamics. We performed WEM simulations on the following host-guest systems: Na+/Cl- ion pair and 4-hydroxy-2-butanone ligand with FK506 binding protein. As a proof of principle, we show that the WEM formalism reproduces the Na+/Cl- ion pair dissociation timescale and the free energy profile obtained from long conventional MD simulation. To increase the accuracy of WEM calculations applied to kinetics and thermodynamics in protein-ligand binding, we introduced a modified WEM scheme called weighted ensemble milestoning with restraint release (WEM-RR), which can increase the number of starting points per milestone without adding additional computational cost. WEM-RR calculations obtained a ligand residence time and binding free energy in agreement with experimental and previous computational results. Moreover, using the milestoning framework, the binding time and rate constants, dissociation constants, and committor probabilities could also be calculated at a low computational cost. We also present an analytical approach for estimating the association rate constant (kon) when binding is primarily diffusion driven. We show that the WEM method can efficiently calculate multiple experimental observables describing ligand-receptor binding/unbinding and is a promising candidate for computer-aided inhibitor design.

Structure of the Bacillus anthracis dTDP-l-rhamnose biosynthetic pathway enzyme: dTDP-α-d-glucose 4,6-dehydratase, RfbB.

Many bacteria require l-rhamnose as a key cell wall component. This sugar is transferred to the cell wall using an activated donor dTDP-l-rhamnose, which is produced by the dTDP-l-rhamnose biosynthetic pathway. We determined the crystal structure of the second enzyme of this pathway dTDP-α-d-glucose 4,6-dehydratase (RfbB) from Bacillus anthracis. Interestingly, RfbB only crystallized in the presence of the third enzyme of the pathway RfbC; however, RfbC was not present in the crystal. Our work represents the first complete structural characterization of the four proteins of this pathway in a single Gram-positive bacterium.

Naphthablins B and C, Meroterpenoids Identified from the Marine Sediment-Derived Streptomyces sp. CP26-58 Using HeLa Cell-Based Cytological Profiling.

HeLa cell-based cytological profiling (CP) was applied to an extract library of marine sediment-derived actinomycetes to discover new cytotoxic secondary metabolites. Among the hit strains, Streptomyces sp. CP26-58 was selected for further investigation to identify its cytotoxic metabolites. CP revealed that the known ionophore tetronasin (1) was responsible for the cytotoxic effect found in the extract. Furthermore, three naphthoquinone meroterpenoids, naphthablin A (2) and two new derivatives designated as naphthablins B (3) and C (4), were isolated from other cytotoxic fractions. The structures of the new compounds were elucidated based on analysis of their HRESIMS and comprehensive NMR data. The absolute configurations of the new compounds were deduced by simulating ECD spectra and calculating potential energies for the model compounds using density function theory (DFT) calculations. Compound 1 showed a significant cytotoxic effect against HeLa cells with an IC50 value of 0.23 µM, and CP successfully clustered 1 with calcium ionophores.

Evaluation of benzoic acid derivatives as sirtuin inhibitors.

Employing a genetically modified yeast strain as a screening tool, 4-dimethylaminobenzoic acid (5) was isolated from the marine sediment-derived Streptomyces sp. CP27-53 as a weak yeast sirtuin (Sir2p) inhibitor. Using this compound as a scaffold, a series of disubstituted benzene derivatives were evaluated to elucidate the structure activity relationships for Sir2p inhibition. The results suggested that 4-alkyl or 4-alkylaminobenzoic acid is the key structure motif for Sir2p inhibitory activity. The most potent Sir2p inhibitor, 4-tert-butylbenzoic acid (20), among the tested compounds in this study turned out to be a weak but selective SIRT1 inhibitor. The calculated binding free energies between the selected compounds and the catalytic domain of SIRT1 were well correlated to their measured SIRT1 inhibitory activities.

Development and Benchmarking of Open Force Field 2.0.0: The Sage Small Molecule Force Field.

We introduce the Open Force Field (OpenFF) 2.0.0 small molecule force field for drug-like molecules, code-named Sage, which builds upon our previous iteration, Parsley. OpenFF force fields are based on direct chemical perception, which generalizes easily to highly diverse sets of chemistries based on substructure queries. Like the previous OpenFF iterations, the Sage generation of OpenFF force fields was validated in protein-ligand simulations to be compatible with AMBER biopolymer force fields. In this work, we detail the methodology used to develop this force field, as well as the innovations and improvements introduced since the release of Parsley 1.0.0. One particularly significant feature of Sage is a set of improved Lennard-Jones (LJ) parameters retrained against condensed phase mixture data, the first refit of LJ parameters in the OpenFF small molecule force field line. Sage also includes valence parameters refit to a larger database of quantum chemical calculations than previous versions, as well as improvements in how this fitting is performed. Force field benchmarks show improvements in general metrics of performance against quantum chemistry reference data such as root-mean-square deviations (RMSD) of optimized conformer geometries, torsion fingerprint deviations (TFD), and improved relative conformer energetics (ΔΔE). We present a variety of benchmarks for these metrics against our previous force fields as well as in some cases other small molecule force fields. Sage also demonstrates improved performance in estimating physical properties, including comparison against experimental data from various thermodynamic databases for small molecule properties such as ΔHmix, ρ(x), ΔGsolv, and ΔGtrans. Additionally, we benchmarked against protein-ligand binding free energies (ΔGbind), where Sage yields results statistically similar to previous force fields. All the data is made publicly available along with complete details on how to reproduce the training results at https://github.com/openforcefield/openff-sage.

Creation of an HDAC-based yeast screening method for evaluation of marine-derived actinomycetes: discovery of streptosetin A.

A histone deacetylase (HDAC)-based yeast assay employing a URA3 reporter gene was applied as a primary screen to evaluate a marine-derived actinomycete extract library and identify human class III HDAC (SIRT) inhibitors. On the basis of the bioassay-guided purification, a new compound designated as streptosetin A (1) was obtained from one of the active strains identified through the yeast assay. The gross structure of the new compound was elucidated from the 1D and 2D NMR data. The absolute stereostructure of 1 was determined based on X-ray crystal structure analysis and simulation of ECD spectra using time-dependent density functional theory calculations. This compound showed weak inhibitory activity against yeast Sir2p and human SIRT1 and SIRT2.

Development and Benchmarking of Open Force Field v1.0.0-the Parsley Small-Molecule Force Field.

We present a methodology for defining and optimizing a general force field for classical molecular simulations, and we describe its use to derive the Open Force Field 1.0.0 small-molecule force field, codenamed Parsley. Rather than using traditional atom typing, our approach is built on the SMIRKS-native Open Force Field (SMIRNOFF) parameter assignment formalism, which handles increases in the diversity and specificity of the force field definition without needlessly increasing the complexity of the specification. Parameters are optimized with the ForceBalance tool, based on reference quantum chemical data that include torsion potential energy profiles, optimized gas-phase structures, and vibrational frequencies. These quantum reference data are computed and are maintained with QCArchive, an open-source and freely available distributed computing and database software ecosystem. In this initial application of the method, we present essentially a full optimization of all valence parameters and report tests of the resulting force field against compounds and data types outside the training set. These tests show improvements in optimized geometries and conformational energetics and demonstrate that Parsley's accuracy for liquid properties is similar to that of other general force fields, as is accuracy on binding free energies. We find that this initial Parsley force field affords accuracy similar to that of other general force fields when used to calculate relative binding free energies spanning 199 protein-ligand systems. Additionally, the resulting infrastructure allows us to rapidly optimize an entirely new force field with minimal human intervention.

Dimerization misalignment in human glutamate-oxaloacetate transaminase variants is the primary factor for PLP release.

The active form of vitamin B6, pyridoxal 5'-phosphate (PLP), plays an essential role in the catalytic mechanism of various proteins, including human glutamate-oxaloacetate transaminase (hGOT1), an important enzyme in amino acid metabolism. A recent molecular and genetic study showed that the E266K, R267H, and P300L substitutions in aspartate aminotransferase, the Arabidopsis analog of hGOT1, genetically suppress a developmentally arrested Arabidopsis RUS mutant. Furthermore, CD analyses suggested that the variants exist as apo proteins and implicated a possible role of PLP in the regulation of PLP homeostasis and metabolic pathways. In this work, we assessed the stability of PLP bound to hGOT1 for the three variant and wildtype (WT) proteins using a combined 6 µs of molecular dynamics (MD) simulation. For the variants and WT in the holo form, the MD simulations reproduced the "closed-open" transition needed for substrate binding. This conformational transition was associated with the rearrangement of the P15-R32 small domain loop providing substrate access to the R387/R293 binding motif. We also showed that formation of the dimer interface is essential for PLP affinity to the active site. The position of PLP in the WT binding site was stabilized by a unique hydrogen bond network of the phosphate binding cup, which placed the cofactor for formation of the covalent Schiff base linkage with K259 for catalysis. The amino acid substitutions at positions 266, 267, and 300 reduced the structural correlation between PLP and the protein active site and/or integrity of the dimer interface. Principal component analysis and energy decomposition clearly suggested dimer misalignment and dissociation for the three variants tested in our work. The low affinity of PLP in the hGOT1 variants observed in our computational work provided structural rationale for the possible role of vitamin B6 in regulating metabolic pathways.

Cadmium(II) inhibition of human uracil-DNA glycosylase by catalytic water supplantation.

Toxic metals are known to inhibit DNA repair but the underlying mechanisms of inhibition are still not fully understood. DNA repair enzymes such as human uracil-DNA glycosylase (hUNG) perform the initial step in the base excision repair (BER) pathway. In this work, we showed that cadmium [Cd(II)], a known human carcinogen, inhibited all activity of hUNG at 100 µM. Computational analyses based on 2 µs equilibrium, 1.6 µs steered molecular dynamics (SMD), and QM/MM MD determined that Cd(II) ions entered the enzyme active site and formed close contacts with both D145 and H148, effectively replacing the catalytic water normally found in this position. Geometry refinement by density functional theory (DFT) calculations showed that Cd(II) formed a tetrahedral structure with D145, P146, H148, and one water molecule. This work for the first time reports Cd(II) inhibition of hUNG which was due to replacement of the catalytic water by binding the active site D145 and H148 residues. Comparison of the proposed metal binding site to existing structural data showed that D145:H148 followed a general metal binding motif favored by Cd(II). The identified motif offered structural insights into metal inhibition of other DNA repair enzymes and glycosylases.

Engineering trypsin for inhibitor resistance.

The development of effective protease therapeutics requires that the proteases be more resistant to naturally occurring inhibitors while maintaining catalytic activity. A key step in developing inhibitor resistance is the identification of key residues in protease-inhibitor interaction. Given that majority of the protease therapeutics currently in use are trypsin-fold, trypsin itself serves as an ideal model for studying protease-inhibitor interaction. To test the importance of several trypsin-inhibitor interactions on the prime-side binding interface, we created four trypsin single variants Y39A, Y39F, K60A, and K60V and report biochemical sensitivity against bovine pancreatic trypsin inhibitor (BPTI) and M84R ecotin. All variants retained catalytic activity against small, commercially available peptide substrates [kcat /KM = (1.2 ± 0.3) × 10(7) M(-1 ) s(-1) . Compared with wild-type, the K60A and K60V variants showed increased sensitivity to BPTI but less sensitivity to ecotin. The Y39A variant was less sensitive to BPTI and ecotin while the Y39F variant was more sensitive to both. The relative binding free energies between BPTI complexes with WT, Y39F, and Y39A were calculated based on 3.5 µs combined explicit solvent molecular dynamics simulations. The BPTI:Y39F complex resulted in the lowest binding energy, while BPTI:Y39A resulted in the highest. Simulations of Y39F revealed increased conformational rearrangement of F39, which allowed formation of a new hydrogen bond between BPTI R17 and H40 of the variant. All together, these data suggest that positions 39 and 60 are key for inhibitor binding to trypsin, and likely more trypsin-fold proteases.

Conformational dynamics of threonine 195 and the S1 subsite in functional trypsin variants.

Replacing the catalytic serine in trypsin with threonine (S195T variant) leads to a nearly complete loss of catalytic activity, which can be partially restored by eliminating the C42-C58 disulfide bond. The 0.69 µs of combined explicit solvent molecular dynamics (MD) simulations revealed continuous rearrangement of T195 with different conformational preferences between five trypsin variants tested. Among three conformational families observed for the T195 residue, one showed the T195 hydroxyl in a conformation analogous to that of the serine residue in wild-type trypsin, positioning the hydroxyl oxygen atom for attack on the carbonyl carbon of the peptide substrate. MD simulations demonstrated that this conformation was more populated for the C42A/C58V/S195T and C42A/C58A/S195T triple variants than for the catalytically inactive S195T variant and correlated with restored enzymatic activities for triple variants. In addition, observation of the increased motion of the S214-G219 segment in the S195T substituted variants suggested an existence of open and closed conformations for the substrate binding pocket. The closed conformation precludes access to the S1 binding site and could further reduce enzymatic activities for triple variants. Double variants with intact serine residues (C42A/C58A/S195 and C42A/C58V/S195) also showed interchange between closed and open conformations for the S214-G219 segment, but to a lesser extent than the triple variants. The increased conformational flexibility of the S1 subsite, which was not observed for the wild-type, correlated with reduced enzymatic activities and suggested a possible mode of substrate regulation for the trypsin variants tested.