Systematic improvement of empirical energy functions in the era of machine learning.

The impact of targeted replacement of individual terms in empirical force fields is quantitatively assessed for pure water, dichloromethane (CH 2 $$ {}_2 $$ Cl 2 $$ {}_2 $$ ), and solvated K + $$ {}^{+} $$ and Cl - $$ {}^{-} $$ ions. For the electrostatic interactions, point charges (PCs) and machine learning (ML)-based minimally distributed charges (MDCM) fitted to the molecular electrostatic potential are evaluated together with electrostatics based on the Coulomb integral. The impact of explicitly including second-order terms is investigated by adding a fragment molecular orbital (FMO)-derived polarization energy to an existing force field, in this case CHARMM. It is demonstrated that anisotropic electrostatics reduce the RMSE for water (by 1.4 kcal/mol), CH 2 $$ {}_2 $$ Cl 2 $$ {}_2 $$ (by 0.8 kcal/mol) and for solvated Cl - $$ {}^{-} $$ clusters (by 0.4 kcal/mol). An additional polarization term can be neglected for CH 2 $$ {}_2 $$ Cl 2 $$ {}_2 $$ but further improves the models for pure water (by ∼ $$ \sim $$ 1.0 kcal/mol) and hydrated Cl - $$ {}^{-} $$ (by 0.4 kcal/mol), and is key for solvated K + $$ {}^{+} $$ , reducing the RMSE by 2.3 kcal/mol. A 12-6 Lennard-Jones functional form performs satisfactorily with PC and MDCM electrostatics, but is not appropriate for descriptions that account for the electrostatic penetration energy. The importance of many-body contributions is assessed by comparing a strictly 2-body approach with self-consistent reference data. Two-body interactions suffice for CH 2 $$ {}_2 $$ Cl 2 $$ {}_2 $$ whereas water and solvated K + $$ {}^{+} $$ and Cl - $$ {}^{-} $$ ions require explicit many-body corrections. Finally, a many-body-corrected dimer potential energy surface exceeds the accuracy attained using a conventional empirical force field, potentially reaching that of an FMO calculation. The present work systematically quantifies which terms improve the performance of an existing force field and what reference data to use for parametrizing these terms in a tractable fashion for ML fitting of pure and heterogeneous systems.

The first HyDRA challenge for computational vibrational spectroscopy.

Vibrational spectroscopy in supersonic jet expansions is a powerful tool to assess molecular aggregates in close to ideal conditions for the benchmarking of quantum chemical approaches. The low temperatures achieved as well as the absence of environment effects allow for a direct comparison between computed and experimental spectra. This provides potential benchmarking data which can be revisited to hone different computational techniques, and it allows for the critical analysis of procedures under the setting of a blind challenge. In the latter case, the final result is unknown to modellers, providing an unbiased testing opportunity for quantum chemical models. In this work, we present the spectroscopic and computational results for the first HyDRA blind challenge. The latter deals with the prediction of water donor stretching vibrations in monohydrates of organic molecules. This edition features a test set of 10 systems. Experimental water donor OH vibrational wavenumbers for the vacuum-isolated monohydrates of formaldehyde, tetrahydrofuran, pyridine, tetrahydrothiophene, trifluoroethanol, methyl lactate, dimethylimidazolidinone, cyclooctanone, trifluoroacetophenone and 1-phenylcyclohexane-cis-1,2-diol are provided. The results of the challenge show promising predictive properties in both purely quantum mechanical approaches as well as regression and other machine learning strategies.

Interaction at a distance: Xenon migration in Mb.

The transport of ligands, such as NO or O2, through internal cavities is essential for the function of globular proteins, including hemoglobin, myoglobin (Mb), neuroglobin, truncated hemoglobins, or cytoglobin. For Mb, several internal cavities (Xe1 through Xe4) were observed experimentally and they were linked to ligand storage. The present work determines barriers for xenon diffusion and relative stabilization energies for the ligand in the initial and final pocket, linking a transition depending on the occupancy state of the remaining pockets from both biased and unbiased molecular dynamics simulations. It is found that the energetics of a particular ligand migration pathway may depend on the direction in which the transition is followed and the occupancy state of the other cavities. Furthermore, the barrier height for a particular transition can depend in a non-additive fashion on the occupancy of either cavity A or B or simultaneous population of both cavities, A and B. Multiple repeats for the Xe1 → Xe2 transition reveal that the activation barrier is a distribution of barrier heights rather than one single value, which is confirmed by a distribution of transition times for the same transition from unbiased simulations. Dynamic cross correlation maps demonstrate that correlated motions occur between adjacent residues or through space, residue Phe138 is found to be a gate for the Xe1 → Xe2 transition, and the volumes of the internal cavities vary along the diffusion pathway, indicating that there is dynamic communication between the ligand and the protein. These findings suggest that Mb is an allosteric protein.

Structural investigation of CDCA3-Cdh1 protein-protein interactions using in vitro studies and molecular dynamics simulation.

The anaphase-promoting complex/cyclosome (APC/C) ubiquitin ligase and its cofactor, Cdh1, regulate the expression of several cell-cycle proteins and their functions during mitosis. Levels of the protein cell division cycle-associated protein 3 (CDCA3), which is functionally required for mitotic entry, are regulated by APC/CCdh1 . CDCA3 is an intrinsically disordered protein and contains both C-terminal KEN box and D-box recognition motifs, enabling binding to Cdh1. Our previous findings demonstrate that CDCA3 has a phosphorylation-dependent non-canonical ABBA-like motif within the linker region bridging these two recognition motifs and is required for efficient binding to Cdh1. Here, we sought to identify and further characterize additional residues that participate within this ABBA-like motif using detailed in vitro experiments and in silico modeling studies. We identified the role of H-bonds, hydrophobic and ionic interactions across the CDCA3 ABBA-like motif in the linker region between KEN and D-box motifs. This linker region adopts a well-defined structure when bound to Cdh1 in the presence of phosphorylation. Upon alanine mutation, the structure of this region is lost, leading to higher flexibility, and alteration in affinities due to binding to alternate sites on Cdh1. Our findings identify roles for the anchoring residues in the non-canonical ABBA-like motif to promote binding to the APC/CCdh1 and regulation of CDCA3 protein levels.

Uncertainty quantification for predictions of atomistic neural networks.

The value of uncertainty quantification on predictions for trained neural networks (NNs) on quantum chemical reference data is quantitatively explored. For this, the architecture of the PhysNet NN was suitably modified and the resulting model (PhysNet-DER) was evaluated with different metrics to quantify its calibration, the quality of its predictions, and whether prediction error and the predicted uncertainty can be correlated. Training on the QM9 database and evaluating data in the test set within and outside the distribution indicate that error and uncertainty are not linearly related. However, the observed variance provides insight into the quality of the data used for training. Additionally, the influence of the chemical space covered by the training data set was studied by using a biased database. The results clarify that noise and redundancy complicate property prediction for molecules even in cases for which changes - such as double bond migration in two otherwise identical molecules - are small. The model was also applied to a real database of tautomerization reactions. Analysis of the distance between members in feature space in combination with other parameters shows that redundant information in the training dataset can lead to large variances and small errors whereas the presence of similar but unspecific information returns large errors but small variances. This was, e.g., observed for nitro-containing aliphatic chains for which predictions were difficult although the training set contained several examples for nitro groups bound to aromatic molecules. The finding underlines the importance of the composition of the training data and provides chemical insight into how this affects the prediction capabilities of a ML model. Finally, the presented method can be used for information-based improvement of chemical databases for target applications through active learning optimization.

Molecular Dynamics with Conformationally Dependent, Distributed Charges.

Accounting for geometry-induced changes in the electronic distribution in molecular simulation is important for capturing effects such as charge flow, charge anisotropy, and polarization. Multipolar force fields have demonstrated their ability to correctly represent chemically significant features such as anisotropy and sigma holes. It has also been shown that off-center point charges offer a compact alternative with similar accuracy. Here, it is demonstrated that allowing relocation of charges within a minimally distributed charge model (MDCM) with respect to their reference atoms is a viable route to capture changes in the molecular charge distribution depending on geometry, i.e., intramolecular polarization. The approach, referred to as "flexible MDCM" (fMDCM), is validated on a number of small molecules and provides accuracies in the electrostatic potential (ESP) of 0.5 kcal/mol on average compared with reference data from electronic structure calculations, whereas MDCM and point charges have root mean squared errors of a factor of 2 to 5 higher. In addition, MD simulations in the NVE ensemble using fMDCM for a box of flexible water molecules with periodic boundary conditions show a width of 0.1 kcal/mol for the fluctuation around the mean at 300 K on the 10 ns time scale. For water, the equilibrium valence angle in the gas phase is found to increase by 2° for simulations in the condensed phase which is consistent with experiment. The accuracy in capturing the geometry dependence of the ESP together with the long-time stability in energy conserving simulations makes fMDCM a promising tool to introduce advanced electrostatics into atomistic simulations.

Impact of the Characteristics of Quantum Chemical Databases on Machine Learning Prediction of Tautomerization Energies.

An essential aspect for adequate predictions of chemical properties by machine learning models is the database used for training them. However, studies that analyze how the content and structure of the databases used for training impact the prediction quality are scarce. In this work, we analyze and quantify the relationships learned by a machine learning model (Neural Network) trained on five different reference databases (QM9, PC9, ANI-1E, ANI-1, and ANI-1x) to predict tautomerization energies from molecules in Tautobase. For this, characteristics such as the number of heavy atoms in a molecule, number of atoms of a given element, bond composition, or initial geometry on the quality of the predictions are considered. The results indicate that training on a chemically diverse database is crucial for obtaining good results and also that conformational sampling can partly compensate for limited coverage of chemical diversity. The overall best-performing reference database (ANI-1x) performs on average by 1 kcal/mol better than PC9, which, however, contains about 2 orders of magnitude fewer reference structures. On the other hand, PC9 is chemically more diverse by a factor of ∼5 as quantified by the number of atom-in-molecule-based fragments (amons) it contains compared with the ANI family of databases. A quantitative measure for deficiencies is the Kullback-Leibler divergence between reference and target distributions. It is explicitly demonstrated that when certain types of bonds need to be covered in the target database (Tautobase) but are undersampled in the reference databases, the resulting predictions are poor. Examples of this include the poor performance of all databases analyzed to predict C(sp2)-C(sp2) double bonds close to heteroatoms and azoles containing N-N and N-O bonds. Analysis of the results with a Tree MAP algorithm provides deeper understanding of specific deficiencies in predicting tautomerization energies by the reference datasets due to inadequate coverage of chemical space. Capitalizing on this information can be used to either improve existing databases or generate new databases of sufficient diversity for a range of machine learning (ML) applications in chemistry.

Transfer Learning to CCSD(T): Accurate Anharmonic Frequencies from Machine Learning Models.

The calculation of the anharmonic modes of small- to medium-sized molecules for assigning experimentally measured frequencies to the corresponding type of molecular motions is computationally challenging at sufficiently high levels of quantum chemical theory. Here, a practical and affordable way to calculate coupled-cluster quality anharmonic frequencies using second-order vibrational perturbation theory (VPT2) from machine-learned models is presented. The approach, referenced as "NN + VPT2", uses a high-dimensional neural network (PhysNet) to learn potential energy surfaces (PESs) at different levels of theory from which harmonic and VPT2 frequencies can be efficiently determined. The NN + VPT2 approach is applied to eight small- to medium-sized molecules (H2CO, trans-HONO, HCOOH, CH3OH, CH3CHO, CH3NO2, CH3COOH, and CH3CONH2) and frequencies are reported from NN-learned models at the MP2/aug-cc-pVTZ, CCSD(T)/aug-cc-pVTZ, and CCSD(T)-F12/aug-cc-pVTZ-F12 levels of theory. For the largest molecules and at the highest levels of theory, transfer learning (TL) is used to determine the necessary full-dimensional, near-equilibrium PESs. Overall, NN + VPT2 yields anharmonic frequencies to within 20 cm-1 of experimentally determined frequencies for close to 90% of the modes for the highest quality PES available and to within 10 cm-1 for more than 60% of the modes. For the MP2 PESs only ∼60% of the NN + VPT2 frequencies were within 20 cm-1 of the experiment, with outliers up to ∼150 cm-1, compared to the experiment. It is also demonstrated that the approach allows to provide correct assignments for strongly interacting modes such as the OH bending and the OH torsional modes in formic acid monomer and the CO-stretch and OH-bend mode in acetic acid.

Elevating CDCA3 levels in non-small cell lung cancer enhances sensitivity to platinum-based chemotherapy.

Platinum-based chemotherapy remains the cornerstone of treatment for most non-small cell lung cancer (NSCLC) cases either as maintenance therapy or in combination with immunotherapy. However, resistance remains a primary issue. Our findings point to the possibility of exploiting levels of cell division cycle associated protein-3 (CDCA3) to improve response of NSCLC tumours to therapy. We demonstrate that in patients and in vitro analyses, CDCA3 levels correlate with measures of genome instability and platinum sensitivity, whereby CDCA3high tumours are sensitive to cisplatin and carboplatin. In NSCLC, CDCA3 protein levels are regulated by the ubiquitin ligase APC/C and cofactor Cdh1. Here, we identified that the degradation of CDCA3 is modulated by activity of casein kinase 2 (CK2) which promotes an interaction between CDCA3 and Cdh1. Supporting this, pharmacological inhibition of CK2 with CX-4945 disrupts CDCA3 degradation, elevating CDCA3 levels and increasing sensitivity to platinum agents. We propose that combining CK2 inhibitors with platinum-based chemotherapy could enhance platinum efficacy in CDCA3low NSCLC tumours and benefit patients.

Identification of Proteins Deregulated by Platinum-Based Chemotherapy as Novel Biomarkers and Therapeutic Targets in Non-Small Cell Lung Cancer.

Platinum-based chemotherapy remains the cornerstone of treatment for most people with non-small cell lung cancer (NSCLC), either as adjuvant therapy in combination with a second cytotoxic agent or in combination with immunotherapy. Resistance to therapy, either in the form of primary refractory disease or evolutionary resistance, remains a significant issue in the treatment of NSCLC. Hence, predictive biomarkers and novel combinational strategies are required to improve the effectiveness and durability of treatment response 6for people with NSCLC. The aim of this study was to identify novel biomarkers and/or druggable proteins from deregulated protein networks within non-oncogene driven disease that are involved in the cellular response to cisplatin. Following exposure of NSCLC cells to cisplatin, in vitro quantitative mass spectrometry was applied to identify altered protein response networks. A total of 65 proteins were significantly deregulated following cisplatin exposure. These proteins were assessed to determine if they are druggable targets using novel machine learning approaches and to identify whether these proteins might serve as prognosticators of platinum therapy. Our data demonstrate novel candidates and drug-like molecules warranting further investigation to improve response to platinum agents in NSCLC.

GlycoTorch Vina: Docking Designed and Tested for Glycosaminoglycans.

Glycosaminoglycans (GAGs) are a family of anionic carbohydrates that play an essential role in the physiology and pathology of all eukaryotic life forms. Experimental determination of GAG-protein complexes is challenging due to their difficult isolation from biological sources, natural heterogeneity, and conformational flexibility-including possible ring puckering of sulfated iduronic acid from 1C4 to 2SO conformation. To overcome these challenges, we present GlycoTorch Vina (GTV), a molecular docking tool based on the carbohydrate docking program VinaCarb (VC). Our program is unique in that it contains parameters to model 2SO sugars while also supporting glycosidic linkages specific to GAGs. We discuss how crystallographic models of carbohydrates can be biased by the choice of refinement software and structural dictionaries. To overcome these variations, we carefully curated 12 of the best available GAG and GAG-like crystal structures (ranging from tetra- to octasaccharides or longer) obtained from the PDB-REDO server and refined using the same protocol. Both GTV and VC produced pose predictions with a mean root-mean-square deviation (RMSD) of 3.1 Å from the native crystal structure-a statistically significant improvement when compared to AutoDock Vina (4.5 Å) and the commercial software Glide (5.9 Å). Examples of how real-space correlation coefficients can be used to better assess the accuracy of docking pose predictions are given. Comparisons between statistical distributions of empirical "salt bridge" interactions, relevant to GAGs, were compared to density functional theory (DFT) studies of model salt bridges, and water-mediated salt bridges; however, there was generally a poor agreement between these data. Water bridges appear to play an important, yet poorly understood, role in the structures of GAG-protein complexes. To aid in the rapid prototyping of future pose scoring functions, we include a module that allows users to include their own torsional and nonbonded parameters.

Assessing Molecular Docking Tools to Guide Targeted Drug Discovery of CD38 Inhibitors.

A promising protein target for computational drug development, the human cluster of differentiation 38 (CD38), plays a crucial role in many physiological and pathological processes, primarily through the upstream regulation of factors that control cytoplasmic Ca2+ concentrations. Recently, a small-molecule inhibitor of CD38 was shown to slow down pathways relating to aging and DNA damage. We examined the performance of seven docking programs for their ability to model protein-ligand interactions with CD38. A test set of twelve CD38 crystal structures, containing crystallized biologically relevant substrates, were used to assess pose prediction. The rankings for each program based on the median RMSD between the native and predicted were Vina, AD4 > PLANTS, Gold, Glide, Molegro > rDock. Forty-two compounds with known affinities were docked to assess the accuracy of the programs at affinity/ranking predictions. The rankings based on scoring power were: Vina, PLANTS > Glide, Gold > Molegro >> AutoDock 4 >> rDock. Out of the top four performing programs, Glide had the only scoring function that did not appear to show bias towards overpredicting the affinity of the ligand-based on its size. Factors that affect the reliability of pose prediction and scoring are discussed. General limitations and known biases of scoring functions are examined, aided in part by using molecular fingerprints and Random Forest classifiers. This machine learning approach may be used to systematically diagnose molecular features that are correlated with poor scoring accuracy.

Cross-Species Analysis of Glycosaminoglycan Binding Proteins Reveals Some Animal Models Are "More Equal" than Others.

Glycosaminoglycan (GAG) mimetics are synthetic or semi-synthetic analogues of heparin or heparan sulfate, which are designed to interact with GAG binding sites on proteins. The preclinical stages of drug development rely on efficacy and toxicity assessment in animals and aim to apply these findings to clinical studies. However, such data may not always reflect the human situation possibly because the GAG binding site on the protein ligand in animals and humans could differ. Possible inter-species differences in the GAG-binding sites on antithrombin III, heparanase, and chemokines of the CCL and CXCL families were examined by sequence alignments, molecular modelling and assessment of surface electrostatic potentials to determine if one species of laboratory animal is likely to result in more clinically relevant data than another. For each protein, current understanding of GAG binding is reviewed from a protein structure and function perspective. This combinatorial analysis shows chemokine dimers and oligomers can present different GAG binding surfaces for the same target protein, whereas a cleft-like GAG binding site will differently influence the types of GAG structures that bind and the species preferable for preclinical work. Such analyses will allow an informed choice of animal(s) for preclinical studies of GAG mimetic drugs.

Pathway Bifurcation in the (4 + 3)/(5 + 2)-Cycloaddition of Butadiene and Oxidopyrylium Ylides: The Significance of Molecular Orbital Isosymmetry.

By drawing analogies from the dimerization of cyclopentadiene, a novel reaction pathway bifurcation is uncovered in the cycloaddition of oxidopyrylium ylides and butadiene. Analysis of the potential energy surface (at the M06-2X/6-311+G(d,p) level of theory) in combination with Born-Oppenheimer molecular dynamics simulations (M06-2X/6-31+G(d)) demonstrate that both the (4 + 3)- and (5 + 2)-cycloaddition products are accessed from the same transition state. Key indicators of a pathway bifurcation (asynchronous bond formation, and a second transition state for the interconversion of the products) are also observed. The absence of a post-transition state bifurcation in the related oxidopyridinium systems of Krenske and Harmata is rationalized. Finally, the isosymmetry of the oxidopyrylium and cyclopentadiene molecular orbitals as well as the presence of "secondary orbital interactions" are emphasized as the common source of nonstatistical behavior. Application of these principles will allow for the rapid identification of new reaction pathway bifurcations.