Experimental and Theoretical Studies on the Reactions of Aliphatic Imines with Isocyanates.

In the context of a project aiming at the replacement of the 3-substituted ß-lactam ring in classical ß-lactam antibiotics by an N(3)-acyl-1,3-diazetidinone moiety, we have investigated the reaction of isocyanates with imines derived from allyl glycinate and differently substituted propionaldehydes. Imines of aromatic aldehydes with anilines have been reported to react with acyl isocyanates to give 1,3-diazetidinones or 2,3-dihydro-4H-1,3,5-oxadiazin-4-ones, via [2+2] or [4+2] cycloaddition, respectively. However, neither of these products was formed with imines derived from allyl glycinate and 2-(mono)methyl propionaldehydes. α,α-Dimethylation of the imine enabled the [4+2] cycloaddition pathway, but the desired 1,3-diazetidinone products were not observed. Surprisingly, the imines obtained from thioesters of 2,2-dimethyl 3-oxo propionic acid reacted with aryl isocyanates or with benzyl isocyanate to give 5,5-dimethyl-2,4-dioxo-6-(aryl-/alkylthio)tetrahydropyrimidines, via thiol displacement and re-addition to a putative six-membered iminium intermediate. These experimental results obtained for the reactions could be rationalized by DFT calculations. In addition, we have shown that N(3)-acyl-1,3-diazetidinone and 2,3-dihydro-4H-1,3,5-oxadiazin-4-one products can be distinguished based on experimental IR data in combination with theoretical reference spectra employing the IR spectra alignment (IRSA) algorithm. This discrimination was not possible by means of 1 H, 13 C, or 15 N NMR spectroscopy.

Combining IC50 or Ki Values from Different Sources Is a Source of Significant Noise.

As part of the ongoing quest to find or construct large data sets for use in validating new machine learning (ML) approaches for bioactivity prediction, it has become distressingly common for researchers to combine literature IC50 data generated using different assays into a single data set. It is well-known that there are many situations where this is a scientifically risky thing to do, even when the assays are against exactly the same target, but the risks of assays being incompatible are even higher when pulling data from large collections of literature data like ChEMBL. Here, we estimate the amount of noise present in combined data sets using cases where measurements for the same compound are reported in multiple assays against the same target. This approach shows that IC50 assays selected using minimal curation settings have poor agreement with each other: almost 65% of the points differ by more than 0.3 log units, 27% differ by more than one log unit, and the correlation between the assays, as measured by Kendall's τ, is only 0.51. Requiring that most of the assay metadata in ChEMBL matches ("maximal curation") in order to combine two assays improves the situation (48% of the points differ by more than 0.3 log units, 13% by more than one log unit, and Kendall's τ is 0.71) at the expense of having smaller data sets. Surprisingly, our analysis shows similar amounts of noise when combining data from different literature Ki assays. We suggest that good scientific practice requires careful curation when combining data sets from different assays and hope that our maximal curation strategy will help to improve the quality of the data that are being used to build and validate ML models for bioactivity prediction. To help achieve this, the code and ChEMBL queries that we used for the maximal curation approach are available as open-source software in our GitHub repository, https://github.com/rinikerlab/overlapping_assays.

Validating Small-Molecule Force Fields for Macrocyclic Compounds Using NMR Data in Different Solvents.

Macrocycles are a promising class of compounds as therapeutics for difficult drug targets due to a favorable combination of properties: They often exhibit improved binding affinity compared to their linear counterparts due to their reduced conformational flexibility, while still being able to adapt to environments of different polarity. To assist in the rational design of macrocyclic drugs, there is need for computational methods that can accurately predict conformational ensembles of macrocycles in different environments. Molecular dynamics (MD) simulations remain one of the most accurate methods to predict ensembles quantitatively, although the accuracy is governed by the underlying force field. In this work, we benchmark four different force fields for their application to macrocycles by performing replica exchange with solute tempering (REST2) simulations of 11 macrocyclic compounds and comparing the obtained conformational ensembles to nuclear Overhauser effect (NOE) upper distance bounds from NMR experiments. Especially, the modern force fields OpenFF 2.0 and XFF yield good results, outperforming force fields like GAFF2 and OPLS/AA. We conclude that REST2 in combination with modern force fields can often produce accurate ensembles of macrocyclic compounds. However, we also highlight examples for which all examined force fields fail to produce ensembles that fulfill the experimental constraints.

Understanding and Quantifying Molecular Flexibility: Torsion Angular Bin Strings.

Molecular flexibility is a commonly used, but not easily quantified term. It is at the core of understanding composition and size of a conformational ensemble and contributes to many molecular properties. For many computational workflows, it is necessary to reduce a conformational ensemble to meaningful representatives, however defining them and guaranteeing the ensemble's completeness is difficult. We introduce the concepts of torsion angular bin strings (TABS) as a discrete vector representation of a conformer's dihedral angles and the number of possible TABS (nTABS) as an estimation for the ensemble size of a molecule, respectively. Here, we show that nTABS corresponds to an upper limit for the size of the conformational space of small molecules and compare the classification of conformer ensembles by TABS with classifications by RMSD. Overcoming known drawbacks like the molecular size dependency and threshold picking of the RMSD measure, TABS is shown to meaningfully discretize the conformational space and hence allows e.g. for fast checks of the coverage of the conformational space. The current proof-of-concept implementation is based on the ETKDGv3 conformer generator as implemented in the RDKit and known torsion preferences extracted from small-molecule crystallographic data.

lwreg: A Lightweight System for Chemical Registration and Data Storage.

Here, we present lwreg, a lightweight, yet flexible chemical registration system supporting the capture of both two-dimensional molecular structures (topologies) and three-dimensional conformers. lwreg is open source, with a simple Python API, and is designed to be easily integrated into computational workflows. In addition to lwreg itself, we also introduce a straightforward schema for storing experimental data and metadata in the registration database. This direct connection between compound structural information and data generated using those structures creates a powerful tool for data analysis and experimental reproducibility. The software is available at and installable directly from https://github.com/rinikerlab/lightweight-registration.

Unraveling motion in proteins by combining NMR relaxometry and molecular dynamics simulations: A case study on ubiquitin.

Nuclear magnetic resonance (NMR) relaxation experiments shine light onto the dynamics of molecular systems in the picosecond to millisecond timescales. As these methods cannot provide an atomically resolved view of the motion of atoms, functional groups, or domains giving rise to such signals, relaxation techniques have been combined with molecular dynamics (MD) simulations to obtain mechanistic descriptions and gain insights into the functional role of side chain or domain motion. In this work, we present a comparison of five computational methods that permit the joint analysis of MD simulations and NMR relaxation experiments. We discuss their relative strengths and areas of applicability and demonstrate how they may be utilized to interpret the dynamics in MD simulations with the small protein ubiquitin as a test system. We focus on the aliphatic side chains given the rigidity of the backbone of this protein. We find encouraging agreement between experiment, Markov state models built in the χ1/χ2 rotamer space of isoleucine residues, explicit rotamer jump models, and a decomposition of the motion using ROMANCE. These methods allow us to ascribe the dynamics to specific rotamer jumps. Simulations with eight different combinations of force field and water model highlight how the different metrics may be employed to pinpoint force field deficiencies. Furthermore, the presented comparison offers a perspective on the utility of NMR relaxation to serve as validation data for the prediction of kinetics by state-of-the-art biomolecular force fields.

DASH properties: Estimating atomic and molecular properties from a dynamic attention-based substructure hierarchy.

Recently, we presented a method to assign atomic partial charges based on the DASH (dynamic attention-based substructure hierarchy) tree with high efficiency and quantum mechanical (QM)-like accuracy. In addition, the approach can be considered "rule based"-where the rules are derived from the attention values of a graph neural network-and thus, each assignment is fully explainable by visualizing the underlying molecular substructures. In this work, we demonstrate that these hierarchically sorted substructures capture the key features of the local environment of an atom and allow us to predict different atomic properties with high accuracy without building a new DASH tree for each property. The fast prediction of atomic properties in molecules with the DASH tree can, for example, be used as an efficient way to generate feature vectors for machine learning without the need for expensive QM calculations. The final DASH tree with the different atomic properties as well as the complete dataset with wave functions is made freely available.

Accounting for Solvation Correlation Effects on the Thermodynamics of Water Networks in Protein Cavities.

Macromolecular recognition and ligand binding are at the core of biological function and drug discovery efforts. Water molecules play a significant role in mediating the protein-ligand interaction, acting as more than just the surrounding medium by affecting the thermodynamics and thus the outcome of the binding process. As individual water contributions are impossible to measure experimentally, a range of computational methods have emerged to identify hydration sites in protein pockets and characterize their energetic contributions for drug discovery applications. Even though several methods model solvation effects explicitly, they focus on determining the stability of specific water sites independently and neglect solvation correlation effects upon replacement of clusters of water molecules, which typically happens in hit-to-lead optimization. In this work, we rigorously determine the conjoint effects of replacing all combinations of water molecules in protein binding pockets through the use of the RE-EDS multistate free-energy method, which combines Hamiltonian replica exchange (RE) and enveloping distribution sampling (EDS). Applications on the small bovine pancreatic trypsin inhibitor and four proteins of the bromodomain family illustrate the extent of solvation correlation effects on water thermodynamics, with the favorability of replacement of the water sites by pharmacophore probes highly dependent on the composition of the water network and the pocket environment. Given the ubiquity of water networks in biologically relevant protein targets, we believe our approach can be helpful for computer-aided drug discovery by providing a pocket-specific and a priori systematic consideration of solvation effects on ligand binding and selectivity.

Accelerating Alchemical Free Energy Prediction Using a Multistate Method: Application to Multiple Kinases.

Alchemical free-energy methods based on molecular dynamics (MD) simulations have become important tools to identify modifications of small organic molecules that improve their protein binding affinity during lead optimization. The routine application of pairwise free-energy methods to rank potential binders from best to worst is impacted by the combinatorial increase in calculations to perform when the number of molecules to assess grows. To address this fundamental limitation, our group has developed replica-exchange enveloping distribution sampling (RE-EDS), a pathway-independent multistate method, enabling the calculation of alchemical free-energy differences between multiple ligands (N > 2) from a single MD simulation. In this work, we apply the method to a set of four kinases with diverse binding pockets and their corresponding inhibitors (42 in total), chosen to showcase the general applicability of RE-EDS in prospective drug design campaigns. We show that for the targets studied, RE-EDS is able to model up to 13 ligands simultaneously with high sampling efficiency, leading to a substantial decrease in computational cost when compared to pairwise methods.

DASH: Dynamic Attention-Based Substructure Hierarchy for Partial Charge Assignment.

We present a robust and computationally efficient approach for assigning partial charges of atoms in molecules. The method is based on a hierarchical tree constructed from attention values extracted from a graph neural network (GNN), which was trained to predict atomic partial charges from accurate quantum-mechanical (QM) calculations. The resulting dynamic attention-based substructure hierarchy (DASH) approach provides fast assignment of partial charges with the same accuracy as the GNN itself, is software-independent, and can easily be integrated in existing parametrization pipelines, as shown for the Open force field (OpenFF). The implementation of the DASH workflow, the final DASH tree, and the training set are available as open source/open data from public repositories.

Improving the IR spectra alignment algorithm with spectra deconvolution and combination with Raman or VCD spectroscopy.

The relative stereochemistry of organic molecules can be determined by comparing theoretical and experimental infrared (IR) spectra of all isomers and assessing the best match. For this purpose, we have recently developed the IR spectra alignment (IRSA) algorithm for automated optimal alignment. IRSA provides a set of quantitative metrics to identify the candidate structure that agrees best with the experimental spectrum. While the correct diastereomer could be determined for the tested sets of rigid and flexible molecules, two issues were identified with more complex compounds that triggered further development. First, strongly overlapping peaks in the IR spectrum are not treated adequately in the current IRSA implementation. Second, the alignment of multiple spectra from different sources (e.g. IR and VCD or Raman) can be improved. In this study, we present an in-depth discussion of these points, followed by the description of modifications to the IRSA algorithm to address them. In particular, we introduce the concept of deconvolution of the experimental and theoretical spectra with a set of pseudo-Voigt bands. The pseudo-Voigt bands have a set of parameters, which can be employed in the alignment algorithm, leading to improved scoring functions. We test the modified algorithm on two data sets. The first set contains compounds with IR and Raman spectra measured in this study, and the second set contains compounds with IR and VCD spectra available in the literature. We show that the algorithm is able to determine the correct diastereomer in all cases. The results highlight that vibrational spectroscopy can be a valuable alternative or complementary method to inform about the stereochemistry of compounds, and the performance of the updated IRSA algorithm suggests that it is a powerful tool for quantitative-based spectral assignments in academia and industry.

Influence of the Fluorophore Mobility on Distance Measurements by Gas-Phase FRET.

Gas-phase Förster resonance energy transfer (FRET) combines mass spectrometry and fluorescence spectroscopy for the conformational analysis of mass-selected biomolecular ions. In FRET, fluorophore pairs are typically covalently attached to a biomolecule using short linkers, which affect the mobility of the dye and the relative orientation of the transition dipole moments of the donor and acceptor. Intramolecular interactions may further influence the range of motion. Yet, little is known about this factor, despite the importance of intramolecular interactions in the absence of a solvent. In this study, we applied transition metal ion FRET (tmFRET) to probe the mobility of a single chromophore pair (Rhodamine 110 and Cu2+) as a function of linker lengths to assess the relevance of intramolecular interactions. Increasing FRET efficiencies were observed with increasing linker length, ranging from 5% (2 atoms) to 28% (13 atoms). To rationalize this trend, we profiled the conformational landscape of each model system using molecular dynamics (MD) simulations. We captured intramolecular interactions that promote a population shift toward smaller donor-acceptor separation for longer linker lengths and induce a significant increase in the acceptor's transition dipole moment. The presented methodology is a first step toward the explicit consideration of a fluorophore's range of motion in the interpretation of gas-phase FRET experiments.

Energy-based clustering: Fast and robust clustering of data with known likelihood functions.

Clustering has become an indispensable tool in the presence of increasingly large and complex datasets. Most clustering algorithms depend, either explicitly or implicitly, on the sampled density. However, estimated densities are fragile due to the curse of dimensionality and finite sampling effects, for instance, in molecular dynamics simulations. To avoid the dependence on estimated densities, an energy-based clustering (EBC) algorithm based on the Metropolis acceptance criterion is developed in this work. In the proposed formulation, EBC can be considered a generalization of spectral clustering in the limit of large temperatures. Taking the potential energy of a sample explicitly into account alleviates requirements regarding the distribution of the data. In addition, it permits the subsampling of densely sampled regions, which can result in significant speed-ups and sublinear scaling. The algorithm is validated on a range of test systems including molecular dynamics trajectories of alanine dipeptide and the Trp-cage miniprotein. Our results show that including information about the potential-energy surface can largely decouple clustering from the sampled density.

Implicit solvent approach based on generalized Born and transferable graph neural networks for molecular dynamics simulations.

Molecular dynamics simulations enable the study of the motion of small and large (bio)molecules and the estimation of their conformational ensembles. The description of the environment (solvent) has, therefore, a large impact. Implicit solvent representations are efficient but, in many cases, not accurate enough (especially for polar solvents, such as water). More accurate but also computationally more expensive is the explicit treatment of the solvent molecules. Recently, machine learning has been proposed to bridge the gap and simulate, in an implicit manner, explicit solvation effects. However, the current approaches rely on prior knowledge of the entire conformational space, limiting their application in practice. Here, we introduce a graph neural network based implicit solvent that is capable of describing explicit solvent effects for peptides with different compositions than those contained in the training set.

Simulation of aqueous solutes using the adaptive solvent-scaling (AdSoS) scheme.

The Adaptive Solvent-Scaling (AdSoS) scheme [J. Chem. Phys. 155 (2021) 094107] is an adaptive-resolution approach for performing simulations of a solute embedded in a fine-grained (FG) solvent region surrounded by a coarse-grained (CG) solvent region, with a continuous FG ↔ CG switching of the solvent resolution across a buffer layer. Instead of relying on a distinct CG solvent model, AdSoS is based on CG models defined by a dimensional scaling of the FG solvent by a factor s, accompanied by the s-dependent modulation of its mass and interaction parameters. The latter changes are designed to achieve an isomorphism between the dynamics of the FG and CG models, and to preserve the dispersive and dielectric solvation properties of the solvent with respect to a solute at FG resolution. As a result, the AdSoS scheme minimizes the thermodynamic mismatch between different regions of the adaptive-resolution system. The present article generalizes the scheme initially introduced for a pure atomic liquid in slab geometry to more practically relevant situations involving (i) a molecular dipolar solvent (e.g., water); (ii) a radial geometry (i.e., spherical rather than planar layers); and (iii) the inclusion of a solute (e.g., water molecule, dipeptide, ion, or ion pair).

Predicting Collagen Triple Helix Stability through Additive Effects of Terminal Residues and Caps.

Collagen model peptides (CMPs) consisting of proline-(2S,4R)-hydroxyproline-glycine (POG) repeats have provided a breadth of knowledge of the triple helical structure of collagen, the most abundant protein in mammals. Predictive tools for triple helix stability have, however, lagged behind since the effect of CMPs with different frames ([POG]n , [OGP]n , or [GPO]n ) and capped or uncapped termini have so far been underestimated. Here, we elucidated the impact of the frame, terminal functional group and its charge on the stability of collagen triple helices. Combined experimental and theoretical studies with frame-shifted, capped and uncapped CMPs revealed that electrostatic interactions, strand preorganization, interstrand H-bonding, and steric repulsion at the termini contribute to triple helix stability. We show that these individual contributions are additive and allow for the prediction of the melting temperatures of CMP trimers.

Frame Shifts Affect the Stability of Collagen Triple Helices.

Collagen model peptides (CMPs), composed of proline-(2S,4R)-hydroxyproline-glycine (POG) repeat units, have been extensively used to study the structure and stability of triple-helical collagen─the dominant structural protein in mammals─at the molecular level. Despite the more than 50-year history of CMPs and numerous studies on the relationship between the composition of single-stranded CMPs and the thermal stability of the assembled triple helices, little attention has been paid to the effects arising from their terminal residues. Here, we show that frame-shifted CMPs, which share POG repeat units but terminate with P, O, or G, form triple helices with vastly different thermal stabilities. A melting temperature difference as high as 16 °C was found for triple helices from 20-mers Ac-OG[POG]6-NH2 and Ac-[POG]6PO-NH2, and triple helices of the constitutional isomers Ac-[POG]7-NH2 and Ac-[GPO]7-NH2 melt 10 °C apart. A combination of thermal denaturation, circular dichroism and NMR spectroscopic studies, and molecular dynamics simulations revealed that the stability differences originate from the propensity of the peptide termini to preorganize into a polyproline-II helical structure. Our results advise that care must be taken when designing peptide mimics of structural proteins, as subtle changes in the terminal residues can significantly affect their properties. Our findings also provide a general and straightforward tool for tuning the stability of CMPs for applications as synthetic materials and biological probes.

Replica-Exchange Enveloping Distribution Sampling Using Generalized AMBER Force-Field Topologies: Application to Relative Hydration Free-Energy Calculations for Large Sets of Molecules.

Free-energy differences between pairs of end-states can be estimated based on molecular dynamics (MD) simulations using standard pathway-dependent methods such as thermodynamic integration (TI), free-energy perturbation, or Bennett's acceptance ratio. Replica-exchange enveloping distribution sampling (RE-EDS), on the other hand, allows for the sampling of multiple end-states in a single simulation without the specification of any pathways. In this work, we use the RE-EDS method as implemented in GROMOS together with generalized AMBER force-field (GAFF) topologies, converted to a GROMOS-compatible format with a newly developed GROMOS++ program amber2gromos, to compute relative hydration free energies for a series of benzene derivatives. The results obtained with RE-EDS are compared to the experimental data as well as calculated values from the literature. In addition, the estimated free-energy differences in water and in vacuum are compared to values from TI calculations carried out with GROMACS. The hydration free energies obtained using RE-EDS for multiple molecules are found to be in good agreement with both the experimental data and the results calculated using other free-energy methods. While all considered free-energy methods delivered accurate results, the RE-EDS calculations required the least amount of total simulation time. This work serves as a validation for the use of GAFF topologies with the GROMOS simulation package and the RE-EDS approach. Furthermore, the performance of RE-EDS for a large set of 28 end-states is assessed with promising results.

Incorporating NOE-Derived Distances in Conformer Generation of Cyclic Peptides with Distance Geometry.

Nuclear magnetic resonance (NMR) data from NOESY (nuclear Overhauser enhancement spectroscopy) and ROESY (rotating frame Overhauser enhancement spectroscopy) experiments can easily be combined with distance geometry (DG) based conformer generators by modifying the molecular distance bounds matrix. In this work, we extend the modern DG based conformer generator ETKDG, which has been shown to reproduce experimental crystal structures from small molecules to large macrocycles well, to include NOE-derived interproton distances. In noeETKDG, the experimentally derived interproton distances are incorporated into the distance bounds matrix as loose upper (or lower) bounds to generate large conformer sets. Various subselection techniques can subsequently be applied to yield a conformer bundle that best reproduces the NOE data. The approach is benchmarked using a set of 24 (mostly) cyclic peptides for which NOE-derived distances as well as reference solution structures obtained by other software are available. With respect to other packages currently available, the advantages of noeETKDG are its speed and that no prior force-field parametrization is required, which is especially useful for peptides with unnatural amino acids. The resulting conformer bundles can be further processed with the use of structural refinement techniques to improve the modeling of the intramolecular nonbonded interactions. The noeETKDG code is released as a fully open-source software package available at www.github.com/rinikerlab/customETKDG.

RestraintMaker: a graph-based approach to select distance restraints in free-energy calculations with dual topology.

The calculation of relative binding free energies (RBFE) involves the choice of the end-state/system representation, of a sampling approach, and of a free-energy estimator. System representations are usually termed "single topology" or "dual topology". As the terminology is often used ambiguously in the literature, a systematic categorization of the system representations is proposed here. In the dual-topology approach, the molecules are simulated as separate molecules. Such an approach is relatively easy to automate for high-throughput RBFE calculations compared to the single-topology approach. Distance restraints are commonly applied to prevent the molecules from drifting apart, thereby improving the sampling efficiency. In this study, we introduce the program RestraintMaker, which relies on a greedy algorithm to find (locally) optimal distance restraints between pairs of atoms based on geometric measures. The algorithm is further extended for multi-state methods such as enveloping distribution sampling (EDS) or multi-site [Formula: see text]-dynamics. The performance of RestraintMaker is demonstrated for toy models and for the calculation of relative hydration free energies. The Python program can be used in script form or through an interactive GUI within PyMol. The selected distance restraints can be written out in GROMOS or GROMACS file formats. Additionally, the program provides a human-readable JSON format that can easily be parsed and processed further. The code of RestraintMaker is freely available on GitHub https://github.com/rinikerlab/restraintmaker.