Benchmarking Quantum Mechanical Levels of Theory for Valence Parametrization in Force Fields.

A wide range of density functional methods and basis sets are available to derive the electronic structure and properties of molecules. Quantum mechanical calculations are too computationally intensive for routine simulation of molecules in the condensed phase, prompting the development of computationally efficient force fields based on quantum mechanical data. Parametrizing general force fields, which cover a vast chemical space, necessitates the generation of sizable quantum mechanical data sets with optimized geometries and torsion scans. To achieve this efficiently, choosing a quantum mechanical method that balances computational cost and accuracy is crucial. In this study, we seek to assess the accuracy of quantum mechanical theory for specific properties such as conformer energies and torsion energetics. To comprehensively evaluate various methods, we focus on a representative set of 59 diverse small molecules, comparing approximately 25 combinations of functional and basis sets against the reference level coupled cluster calculations at the complete basis set limit.

A Polarizable QM/MM Model That Combines the State-Averaged CASSCF and AMOEBA Force Field for Photoreactions in Proteins.

This study presents the polarizable quantum mechanics/molecular mechanics (QM/MM) embedding of the state-averaged complete active space self-consistent field (SA-CASSCF) in the atomic multipole optimized energetics for biomolecular applications (AMOEBA) force field for the purpose of studying photoreactions in protein environments. We describe two extensions of our previous work that combine SA-CASSCF with AMOEBA water models, allowing it to be generalized to AMOEBA models for proteins and other macromolecules. First, we discuss how our QM/MM model accounts for the discrepancy between the direct and polarization electric fields that arises in the AMOEBA description of intramolecular polarization. A second improvement is the incorporation of link atom schemes to treat instances in which the QM/MM boundary goes through covalent bonds. A single-link atom scheme and double-link atom scheme are considered in this work, and we will discuss how electrostatic interaction, van der Waals interaction, and various kinds of valence terms are treated across the boundary. To test the accuracy of the link atom scheme, we will compare QM/MM with full QM calculations and study how the errors in ground state properties, excited state properties, and excitation energies change when tuning the parameters in the link atom scheme. We will also test the new SA-CASSCF/AMOEBA method on an elementary reaction step in NanoLuc, an artificial bioluminescence luciferase. We will show how the reaction mechanism is different when calculated in the gas phase, in polarizable continuum medium (PCM), versus in protein AMOEBA models.

Revealing elusive conformations of sucrose from hydrogen bond J-coupling in H2O: A combined NMR and quantum mechanics study.

Hydrogen bonding is a crucial feature of biomolecules, but its characterization in glycans dissolved in aqueous solutions is challenging due to rapid hydrogen exchange between hydroxyl groups and H2O. In principle, the scalar (J) coupling constant can reveal the relative orientation of the atoms in the molecule. In contrast to J-coupling through H-bonds reported in proteins and nucleic acids, research on J-coupling through H-bonds in glycans dissolved in water is lacking. Here, we use sucrose as a model system for H-bonding studies; its structure, which consists of glucose (Glc) and fructose (Frc), is well-studied, and it is readily available. We apply the in-phase, antiphase-HSQC-TOCSY and quantify previously unreported through H-bond J-values for Frc-OH1-Glc-OH2 in H2O. While earlier reports of Brown and Levy indicate this H-bond as having only a single direction, our reported findings indicate the potential presence of two involving these same atoms, namely, G2OH âž” F1O and F1OH âž” G2O (where F and G stand for Frc and Glc, respectively). The calculated density functional theory J-values for the G2OH âž” F1O agree with the experimental values. Additionally, we detected four other possible H-bonds in sucrose, which require different phi, psi (ϕ, ψ) torsion angles. The ϕ, ψ values are consistent with previous predictions of du Penhoat et al. and Venable et al. Our results will provide new insights into the molecular structure of sucrose and its interactions with proteins.

Temperature-Dependent Fold-Switching Mechanism of the Circadian Clock Protein KaiB.

The oscillator of the cyanobacterial circadian clock relies on the ability of the KaiB protein to switch reversibly between a stable ground-state fold (gsKaiB) and an unstable fold-switched fold (fsKaiB). Rare fold-switching events by KaiB provide a critical delay in the negative feedback loop of this post-translational oscillator. In this study, we experimentally and computationally investigate the temperature dependence of fold switching and its mechanism. We demonstrate that the stability of gsKaiB increases with temperature compared to fsKaiB and that the Q10 value for the gsKaiB → fsKaiB transition is nearly three times smaller than that for the reverse transition. Simulations and native-state hydrogen-deuterium exchange NMR experiments suggest that fold switching can involve both subglobally and near-globally unfolded intermediates. The simulations predict that the transition state for fold switching coincides with isomerization of conserved prolines in the most rapidly exchanging region, and we confirm experimentally that proline isomerization is a rate-limiting step for fold switching. We explore the implications of our results for temperature compensation, a hallmark of circadian clocks, through a kinetic model.

Tuning Potential Functions to Host-Guest Binding Data.

Software to more rapidly and accurately predict protein-ligand binding affinities is of high interest for early-stage drug discovery, and physics-based methods are among the most widely used technologies for this purpose. The accuracy of these methods depends critically on the accuracy of the potential functions that they use. Potential functions are typically trained against a combination of quantum chemical and experimental data. However, although binding affinities are among the most important quantities to predict, experimental binding affinities have not to date been integrated into the experimental data set used to train potential functions. In recent years, the use of host-guest complexes as simple and tractable models of binding thermodynamics has gained popularity due to their small size and simplicity, relative to protein-ligand systems. Host-guest complexes can also avoid ambiguities that arise in protein-ligand systems such as uncertain protonation states. Thus, experimental host-guest binding data are an appealing additional data type to integrate into the experimental data set used to optimize potential functions. Here, we report the extension of the Open Force Field Evaluator framework to enable the systematic calculation of host-guest binding free energies and their gradients with respect to force field parameters, coupled with the curation of 126 host-guest complexes with available experimental binding free energies. As an initial application of this novel infrastructure, we optimized generalized Born (GB) cavity radii for the OBC2 GB implicit solvent model against experimental data for 36 host-guest systems. This refitting led to a dramatic improvement in accuracy for both the training set and a separate test set with 90 additional host-guest systems. The optimized radii also showed encouraging transferability from host-guest systems to 59 protein-ligand systems. However, the new radii are significantly smaller than the baseline radii and lead to excessively favorable hydration free energies (HFEs). Thus, users of the OBC2 GB model currently may choose between GB cavity radii that yield more accurate binding affinities and GB cavity radii that yield more accurate HFEs. We suspect that achieving good accuracy on both will require more far-reaching adjustments to the GB model. We note that binding free-energy calculations using the OBC2 model in OpenMM gain about a 10× speedup relative to corresponding explicit solvent calculations, suggesting a future role for implicit solvent absolute binding free-energy (ABFE) calculations in virtual compound screening. This study proves the principle of using host-guest systems to train potential functions that are transferrable to protein-ligand systems and provides an infrastructure that enables a range of applications.

Copper(II) Affects the Biochemical Behavior of Proinsulin C-peptide by Forming Ternary Complexes with Serum Albumin.

Peptide hormones are essential signaling molecules with therapeutic importance. Identifying regulatory factors that drive their activity gives important insight into their mode of action and clinical development. In this work, we demonstrate the combined impact of Cu(II) and the serum protein albumin on the activity of C-peptide, a 31-mer peptide derived from the same prohormone as insulin. C-peptide exhibits beneficial effects, particularly in diabetic patients, but its clinical use has been hampered by a lack of mechanistic understanding. We show that Cu(II) mediates the formation of ternary complexes between albumin and C-peptide and that the resulting species depend on the order of addition. These ternary complexes notably alter peptide activity, showing differences from the peptide or Cu(II)/peptide complexes alone in redox protection as well as in cellular internalization of the peptide. In standard clinical immunoassays for measuring C-peptide levels, the complexes inflate the quantitation of the peptide, suggesting that such adducts may affect biomarker quantitation. Altogether, our work points to the potential relevance of Cu(II)-linked C-peptide/albumin complexes in the peptide's mechanism of action and application as a biomarker.

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.

The impact of conformational sampling on first-principles calculations of vicinal COCH J-couplings in carbohydrates.

Dihedral angles in organic molecules and biomolecules are vital structural parameters that can be indirectly probed by nuclear magnetic resonance (NMR) measurements of vicinal J-couplings. The empirical relations that map the measured couplings to dihedral angles are typically determined by fitting using static structural models, but this neglects the effects of thermal fluctuations at the finite temperature conditions under which NMR measurements are often taken. In this study, we calculate ensemble-averaged J-couplings for several structurally rigid carbohydrate derivatives using first-principles molecular dynamics simulations to sample the thermally accessible conformations around the minimum energy structure. Our results show that including thermal fluctuation effects significantly shifts the predicted couplings relative to single-point calculations at the energy minima, leading to improved agreement with experiments. This provides evidence that accounting for conformational sampling in first-principles calculations can improve the accuracy of NMR-based structure determination for structurally complex carbohydrates.

Benchmarking Density Functionals, Basis Sets, and Solvent Models in Predicting Thermodynamic Hydricities of Organic Hydrides.

Many renewable energy technologies, such as hydrogen gas synthesis and carbon dioxide reduction, rely on chemical reactions involving hydride anions (H-). When selecting molecules to be used in such applications, an important quantity to consider is the thermodynamic hydricity, which is the free energy required for a species to donate a hydride anion. Theoretical calculations of thermodynamic hydricity depend on several parameters, mainly the density functional, basis set, and solvent model. In order to assess the effects of the above three parameters, we carry out hydricity calculations with different combinations of density functionals, basis sets, and solvent models for a set of organic molecules with known experimental hydricity values. The data are analyzed by comparing the R2 and root-mean-squared error (RMSE) of linear fits with a fixed slope of 1 and using the Akaike Information Criterion to determine statistical significance of the RMSE rank ordering. Based on these results, we quantified the accuracy of theoretical predictions of hydricity and found that the best compromise between accuracy and computational cost was obtained by using the B3LYP-D3 density functional for the geometry optimization and free-energy corrections, either ωB97X-D3 or M06-2X-D3 for single-point energy corrections, combined with a basis set no larger than def-TZVP and the C-PCM ISWIG solvation model. At this level of theory, the RMSEs of hydricity calculations for organic molecules in acetonitrile and dimethyl sulfoxide were found to be <4 and <10 kcal/mol, respectively, for an experimental data set with a dynamic range of 20-150 kcal/mol.

O-Acetyl Migration within the Sialic Acid Side Chain: A Mechanistic Study Using the Ab Initio Nanoreactor.

Many disease-causing viruses target sialic acids on the surface of host cells. Some viruses bind preferentially to sialic acids with O-acetyl modification at the hydroxyl group of C7, C8, or C9 on the glycerol-like side chain. Studies of proteins binding to sialosides containing O-acetylated sialic acids are crucial in understanding the related diseases but experimentally difficult due to the lability of the ester group. We recently showed that O-acetyl migration among hydroxyl groups of C7, C8, and C9 in sialic acids occurs in all directions in a pH-dependent manner. In the current study, we elucidate a full mechanistic pathway for the migration of O-acetyl among C7, C8, and C9. We used an ab initio nanoreactor to explore potential reaction pathways and density functional theory, pKa calculations, and umbrella sampling to investigate elementary steps of interest. We found that when a base is present, migration is easy in any direction and involves three key steps: deprotonation of the hydroxyl group, cyclization between the two carbons, and the migration of the O-acetyl group. This dynamic equilibrium may play a defensive role against pathogens that evolve to gain entry to the cell by binding selectively to one acetylation state.

SARS-CoV-2 and MERS-CoV Spike Protein Binding Studies Support Stable Mimic of Bound 9-O-Acetylated Sialic Acids.

Many disease-causing viruses target sialic acids (Sias), a class of nine-carbon sugars known to coat the surface of many cells, including those in the lungs. Human beta coronaviridae, known for causing respiratory tract diseases, often bind Sias, and some preferentially bind to those with 9-O-Ac-modification. Currently, co-binding of SARS-CoV-2, a beta coronavirus responsible for the COVID-19 pandemic, to human Sias has been reported and its preference towards α2-3-linked Neu5Ac has been shown. Nevertheless, O-acetylated Sias-protein binding studies are difficult to perform, due to the ester lability. We studied the binding free energy differences between Neu5,9Ac2α2-3GalßpNP and its more stable 9-NAc mimic binding to SARS-CoV-2 spike protein using molecular dynamics and alchemical free energy simulations. We identified multiple Sia-binding pockets, including two novel sites, with similar binding affinities to those of MERS-CoV, a known co-binder of sialic acid. In our binding poses, 9-NAc and 9-OAc Sias bind similarly, suggesting an experimentally reasonable mimic to probe viral mechanisms.

Multireference configuration interaction study of the predissociation of C2 via its F1Πu state.

Photodissociation is one of the main destruction pathways for dicarbon (C2) in astronomical environments, such as diffuse interstellar clouds, yet the accuracy of modern astrochemical models is limited by a lack of accurate photodissociation cross sections in the vacuum ultraviolet range. C2 features a strong predissociative F1Πu-X1Σg + electronic transition near 130 nm originally measured in 1969; however, no experimental studies of this transition have been carried out since, and theoretical studies of the F1Πu state are limited. In this work, potential energy curves of excited electronic states of C2 are calculated with the aim of describing the predissociative nature of the F1Πu state and providing new ab initio photodissociation cross sections for astrochemical applications. Accurate electronic calculations of 56 singlet, triplet, and quintet states are carried out at the DW-SA-CASSCF/MRCI+Q level of theory with a CAS(8,12) active space and the aug-cc-pV5Z basis set augmented with additional diffuse functions. Photodissociation cross sections arising from the vibronic ground state to the F1Πu state are calculated by a coupled-channel model. The total integrated cross section through the F1Πu v = 0 and v = 1 bands is 1.198 × 10-13 cm2 cm-1, giving rise to a photodissociation rate of 5.02 × 10-10 s-1 under the standard interstellar radiation field, much larger than the rate in the Leiden photodissociation database. In addition, we report a new 21Σu + state that should be detectable via a strong 21Σu +-X1Σg + band around 116 nm.

Exploration and validation of force field design protocols through QM-to-MM mapping.

The scale of the parameter optimisation problem in traditional molecular mechanics force field construction means that design of a new force field is a long process, and sub-optimal choices made in the early stages can persist for many generations. We hypothesise that careful use of quantum mechanics to inform molecular mechanics parameter derivation (QM-to-MM mapping) should be used to significantly reduce the number of parameters that require fitting to experiment and increase the pace of force field development. Here, we design and train a collection of 15 new protocols for small, organic molecule force field derivation, and test their accuracy against experimental liquid properties. Our best performing model has only seven fitting parameters, yet achieves mean unsigned errors of just 0.031 g cm-3 and 0.69 kcal mol-1 in liquid densities and heats of vaporisation, compared to experiment. The software required to derive the designed force fields is freely available at https://github.com/qubekit/QUBEKit.

Source of Rate Acceleration for Carbocation Cyclization in Biomimetic Supramolecular Cages.

The results of quantum chemical and molecular dynamics calculations reveal that polyanionic gallium-based cages accelerate cyclization reactions of pentadienyl alcohols as a result of substrate cage interactions, preferential binding of reactive conformations of substrate/H3O+ pairs, and increased substrate basicity. However, the increase in basicity dominates. Experimental structure-activity relationship studies in which the metal vertices and overall charge of the cage are varied confirm the model derived via calculations.

Improving Force Field Accuracy by Training against Condensed-Phase Mixture Properties.

Developing a sufficiently accurate classical force field representation of molecules is key to realizing the full potential of molecular simulations as a route to gaining a fundamental insight into a broad spectrum of chemical and biological phenomena. This is only possible, however, if the many complex interactions between molecules of different species in the system are accurately captured by the model. Historically, the intermolecular van der Waals (vdW) interactions have primarily been trained against densities and enthalpies of vaporization of pure (single-component) systems, with occasional usage of hydration free energies. In this study, we demonstrate how including physical property data of binary mixtures can better inform these parameters, encoding more information about the underlying physics of the system in complex chemical mixtures. To demonstrate this, we retrain a select number of Lennard-Jones parameters describing the vdW interactions of the OpenFF 1.0.0 (Parsley) fixed charge force field against training sets composed of densities and enthalpies of mixing for binary liquid mixtures as well as densities and enthalpies of vaporization of pure liquid systems and assess the performance of each of these combinations. We show that retraining against the mixture data improves the force field's ability to reproduce mixture properties, including solvation free energies, correcting some systematic errors that exist when training vdW interactions against properties of pure systems only.

Open Force Field Evaluator: An Automated, Efficient, and Scalable Framework for the Estimation of Physical Properties from Molecular Simulation.

Developing accurate classical force field representations of molecules is key to realizing the full potential of molecular simulations, both as a powerful route to gaining fundamental insights into a broad spectrum of chemical and biological phenomena and for predicting physicochemical and mechanical properties of substances. The Open Force Field Consortium is an industry-funded open science effort to this end, developing open-source tools for rapidly generating new high-quality small-molecule force fields. An integral aspect of this is the parameterization and assessment of force fields against high-quality, condensed-phase physical property data, curated from open data sources such as the NIST ThermoML Archive, alongside quantum chemical data. The quantity of such experimental data in open data archives alone would require an onerous amount of human and computational resources to both curate and estimate manually, especially when estimations must be obtained for numerous sets of force field parameters. Here, we present an entirely automated, highly scalable framework for evaluating physical properties and their gradients in terms of force field parameters. It is written as a modular and extensible Python framework, which employs an intelligent multiscale estimation approach that allows for the automated estimation of properties from simulation and cached simulation data, and a pluggable API for estimation of new properties. In this study, we demonstrate the utility of the framework by benchmarking the OpenFF 1.0.0 small-molecule force field and GAFF 1.8 and GAFF 2.1 force fields against a test set of binary density and enthalpy of mixing measurements curated using the framework utilities. Further, we demonstrate the framework's utility as part of force field optimization by using it alongside ForceBalance, a framework for systematic force field optimization, to retrain a set of nonbonded van der Waals parameters against a training set of density and enthalpy of vaporization measurements.

Proposed Mechanism for the Biosynthesis of the [FeFe] Hydrogenase H-Cluster: Central Roles for the Radical SAM Enzymes HydG and HydE.

Radical S-adenosylmethionine (radical SAM or rSAM) enzymes use their S-adenosylmethionine cofactor bound to a unique Fe of a [4Fe-4S] cluster to generate the "hot" 5'-deoxyadenosyl radical, which drives highly selective radical reactions via specific interactions with a given rSAM enzyme's substrate. This Perspective focuses on the two rSAM enzymes involved in the biosynthesis of the organometallic H-cluster of [FeFe] hydrogenases. We present here a detailed sequential model initiated by HydG, which lyses a tyrosine substrate via a 5'-deoxyadenosyl H atom abstraction from those amino acid's amino group, initially producing dehydroglycine and an oxidobenzyl radical. In this model, two successive radical cascade reactions lead ultimately to the formation of HydG's product, a mononuclear Fe organometallic complex: [Fe(II)(CN)(CO)2(cysteinate)]-, with the iron originating from a unique "dangler" Fe coordinated by a cysteine ligand providing a sulfur bridge to another [4Fe-4S] auxiliary cluster in the enzyme. In turn, in this model, [Fe(II)(CN)(CO)2(cysteinate)]- is the substrate for HydE, the second rSAM enzyme in the biosynthetic pathway, which activates this mononuclear organometallic unit for dimerization, forming a [Fe2S2(CO)4(CN)2] precursor to the [2Fe] H component of the H-cluster, requiring only the completion of the bridging azadithiolate (SCH2NHCH2S) ligand. This model is built upon a foundation of data that incorporates cell-free synthesis, isotope sensitive spectroscopies, and the selective use of synthetic complexes substituting for intermediates in the enzymatic "assembly line". We discuss controversies pertaining to this model and some remaining open issues to be addressed by future work.

Evaluating the Accuracy of the QCEIMS Approach for Computational Prediction of Electron Ionization Mass Spectra of Purines and Pyrimidines.

Mass spectrometry is the most commonly used method for compound annotation in metabolomics. However, most mass spectra in untargeted assays cannot be annotated with specific compound structures because reference mass spectral libraries are far smaller than the complement of known molecules. Theoretically predicted mass spectra might be used as a substitute for experimental spectra especially for compounds that are not commercially available. For example, the Quantum Chemistry Electron Ionization Mass Spectra (QCEIMS) method can predict 70 eV electron ionization mass spectra from any given input molecular structure. In this work, we investigated the accuracy of QCEIMS predictions of electron ionization (EI) mass spectra for 80 purine and pyrimidine derivatives in comparison to experimental data in the NIST 17 database. Similarity scores between every pair of predicted and experimental spectra revealed that 45% of the compounds were found as the correct top hit when QCEIMS predicted spectra were matched against the NIST17 library of >267,000 EI spectra, and 74% of the compounds were found within the top 10 hits. We then investigated the impact of matching, missing, and additional fragment ions in predicted EI mass spectra versus ion abundances in MS similarity scores. We further include detailed studies of fragmentation pathways such as retro Diels-Alder reactions to predict neutral losses of (iso)cyanic acid, hydrogen cyanide, or cyanamide in the mass spectra of purines and pyrimidines. We describe how trends in prediction accuracy correlate with the chemistry of the input compounds to better understand how mechanisms of QCEIMS predictions could be improved in future developments. We conclude that QCEIMS is useful for generating large-scale predicted mass spectral libraries for identification of compounds that are absent from experimental libraries and that are not commercially available.