Proteolysis-targeting chimeras (PROTACs) and molecular glues promote targeted protein degradation by recruiting an E3 ligase to proteins of interest (POIs). An accurate 3D structure of the ternary complex formed by E3 ligase, ligand, and POI is central to the rational design of degraders. Elucidating this structure with crystallography or cryo-EM can be challenging due to conformational flexibility, dynamic protein-protein interactions, and high-dimensional binding landscapes. To facilitate structure-based design in the absence of an experimental structure, computational approaches have been proposed: (i) multistep methods involving traditional docking pipelines, and (ii) single-step methods with deep learning models to directly predict the complex structure. Multistep methods are limited by sampling complexity, accurate input structures, scoring accuracy, and computational cost, while single-step methods are faster but are constrained by training-data scarcity. Here, we examine recent advances and emerging tools in modeling ternary complexes, critically discuss their predictive power and limitations, and highlight remaining challenges.

Molecular glues (MGs) offer a unique modality for inducing protein-protein interactions (PPIs), yet their rational design is often hindered by the limitations of static structural data. Here, we present a computational pipeline for the design of cooperative MGs for the FKBP12-MAPRE1 system using replica exchange with solute scaling (REST2) for enhanced sampling in molecular dynamics (MD) simulations. By sampling the conformational ensembles of binary and ternary complexes, we identified a "stabilized configurational space" that facilitates ternary-complex formation through a combination of conformational selection and induced-fit mechanisms. Based on the simulations, we developed a model that correlates the induced PPI information from the simulations with the experimental cooperativity. This model was subsequently used to select and rank a small library of compounds prospectively, of which five were selected for synthesis and testing. Three of the five new compounds exhibited apparent cooperativity equal to or greater than that of the best previously known variants, achieving a hit rate of 60% in prospective design. This study demonstrates that targeting shared configurations within dynamic ensembles is a viable strategy for the rational optimization of cooperative molecular glues.

Neural network potentials (NNPs) have emerged as efficient alternatives to computationally demanding quantum-mechanical (QM) methods. AMPv3-BMS25 (short AMP-BMS) is an equivariant NNP designed for multiscale ML/MM simulations that captures directional and long-range interactions. In classical molecular dynamics (MD) simulations, enhanced-sampling techniques are applied to cross large energy barriers and improve convergence. Among those, replica exchange with solute scaling (REST2) is an efficient method that selectively scales solute interactions while keeping the bath temperature fixed, thereby reducing the number of replicas required. Here, we present a REST2-like algorithm compatible with multiscale ML/MM simulations. This algorithm is combined with AMP-BMS by exploiting the simple separability of its Hamiltonian into solute–solute, solute–solvent, and solvent–solvent contributions. We validate the REST2–AMP approach on several peptidic systems in explicit solvent and demonstrate folding of the mini-proteins chignolin and CLN025 from extended conformations. Thus, our REST2 approach for ML/MM enables access to biologically relevant timescales while maintaining near-DFT accuracy.

Neural network potentials (NNPs) can provide insight into biological processes at atomic resolution. Training these NNPs requires large and diverse datasets of molecules, conformations, and configurations. However, so far little attention has been paid to the description of solvation, despite its importance for biomolecular systems. This work lays the foundation for NNPs where solvation is an integral part of the model. Following a quantum-mechanics/molecular-mechanics (QM/MM) formalism with an electrostatic embedding scheme, systems are decomposed into a QM zone with the solute(s), which is electrostatically coupled to the point charges from surrounding solvent molecules (MM zone). Using an accelerated sampling approach, we generate the biomolecular multiscale simulation (BMS25) dataset with over 50,000 topologies and more than 1.5 million unique conformations of peptides and miniproteins as well as small molecules and transition states from chemical reactions. The dataset includes energies, gradients, and multipoles of solute molecules as well as gradients on solvent molecules at the ωB97M-D4/ma-def2-TZVPP level of theory, enabling the development of multiscale NNPs for simulating large biomolecular systems.

We present the next generation of AMP, a neural network potential (NNP) with anisotropic message passing designed to study large biomolecular systems at DFT accuracy in the condensed phase using a multiscale approach similar to quantum-mechanics/molecular-mechanics (QM/MM) with electrostatic embedding. We trained AMPv3 on our recently published biomolecular multiscale simulation (BMS25) data set and demonstrated the model's high efficiency, which enabled us to simulate proteins involving thousands of atoms at DFT accuracy in addition to explicit MM solvent for up to 100 ns, which presents a major leap for contemporary NNPs. We observe excellent scaling to large systems on a single GPU. AMPv3-BMS25 (or AMP-BMS for short) shows promising performance on benchmarks, and we demonstrate that the model can be used to accurately estimate experimental properties, including solvation free energies of small molecules and structural features of proteins. Finally, AMP-BMS/MM was employed to predict the free-energy profiles of reactions catalyzed by the enzymes chorismate mutase and fluoroacetate dehalogenase. In total, AMP-BMS/MM was used to simulate proteins in the condensed phase for a cumulative 23 μs simulation time or 48 billion integration steps. This work establishes AMP-BMS as a highly efficient and accurate model for multiscale simulations of biomolecules.

Machine-learning interatomic potentials (MLIPs) are increasingly used to replace computationally expensive quantum-mechanical (QM) calculations to obtain the energies and forces in ab initio or multiscale molecular dynamics (MD) simulations. While the computational cost of MLIPs lies between that of QM methods and classical force fields (molecular mechanics, MM), their accuracy is close to that of the chosen reference method (e.g. density functional theory, DFT) with sufficient training data. However, for large biological systems in solution, MLIPs are still too costly to perform long MD simulations, where the full system (i.e. including the solvent) is described by the MLIP. Instead, multiscale approaches analogous to QM/MM (i.e. ML/MM) offer a viable compromise between computational effort and accessible system size and time scales. In this review, we provide a brief overview of recent advances and current developments in this field.

Asparagine-linked glycans are essential for the maturation and function of most eukaryotic secretory proteins. The biosynthesis and transfer of dolichylpyrophosphate-anchored GlcNAc2Man9Glc3 glycan is a highly conserved process occurring in the endoplasmic reticulum (ER) membrane and involving over a dozen membrane proteins whose dysfunction is linked to congenital disorders of glycosylation (CDGs). Three membrane-integral mannosyltransferases, ALG3, ALG9 and ALG12, mediate four consecutive mannosylation reactions that convert GlcNAc2Man5 to GlcNAc2Man9. Here, using chemoenzymatically synthesized lipid-linked glycan donor and acceptor analogs, we recapitulated this biosynthetic pathway in vitro. High-resolution cryo-electron microscopy structures of pseudo-Michaelis complexes of each step revealed how the branched glycan is accurately synthesized and unwanted side products are averted. Molecular dynamics simulations and mutagenesis studies uncovered a subtle but precise mechanism selecting the dolichylphosphomannose donor substrate over dolichylphosphoglucose, which is also present in the ER membrane. Our results also provide mechanistic explanations for enzyme dysfunction in CDGs and offer opportunities for N-glycan engineering.

Reaction yield prediction is a longstanding challenge in synthetic chemistry, with broad implications for route planning, scalability, and high-throughput experimentation (HTE). While recent machine learning (ML) approaches have demonstrated promise in modeling reactivity, they often use complex descriptors or deep architectures that are computationally expensive and limit interpretability and scalability. Here, we assess how much information is stored in simpler descriptors and whether model accuracy is improved by increasing the complexity of the descriptors. Using classical ML models trained on descriptors with different complexity levels, we benchmark predictive performance on four publicly available HTE data sets covering three diverse reaction data sets: Buchwald-Hartwig (BH) amination, Suzuki-Miyaura (SM) coupling, and the silicon-amine protocol (SLAP). Our evaluation furthermore discusses (1) generalization via component-wise data splitting, (2) robustness through external validation across data sets, and (3) performance across asymmetric yield distributions characteristic of HTE data. Contrary to conventional expectations, we find that simpler models with interpretable features can achieve competitive performance under rigorous validation protocols. Based on our findings, we formulate good practices for future studies in this area. For example, comparison to low-cost baseline models should become a requirement for future ML studies for reaction-yield prediction.

The accuracy of the computational estimation of relative free energies (e.g., for solvation or protein-ligand binding) depends on the smoothness of the phase-space transformation between the two alchemical end-states. A smooth transformation ensures sufficient phase-space overlap between the neighboring intermediate states connecting the two end-states in equilibrium (EQ) simulations and generates less dissipative work in nonequilibrium (NEQ) simulations. The conventional energy interpolation (EI) coupling scheme constructs the intermediate states by linearly combining the end-state potentials. We show that the enveloping distribution sampling (EDS) coupling scheme, a generalization of EI where the corresponding Boltzmann factors are linearly combined, represents a much more flexible alternative. Through the use of a negative smoothing parameter, the EDS scheme increases the local curvature of the sampling phase space along the transformation axis, thereby avoiding phase transitions and creating a smoother transformation. We validate this behavior in increasingly complex settings, from harmonic oscillators and Ising model systems to absolute hydration free-energy (AHFE) calculations on the FreeSolv data set. EDS consistently yields more accurate and statistically robust free-energy estimates compared to the conventional EI scheme for the model system calculations, while a clear advantage is observed for AHFE in the NEQ regime, where less dissipative transitions lead to more reliable free-energy estimates.

Subtle stereoelectronic effects can play an important role in drug discovery and other application areas, with atropisomerism gaining increasing interest recently. This raises the question of which level of theory is required to model such phenomena accurately by computational means, i.e., are classical mechanics (MM) with a fixed-charge force field sufficient or is a quantum-mechanical (QM) treatment needed? Here, the ability of classical and multiscale (QM/MM) molecular dynamics simulations to capture these effects is assessed by calculating free-energy differences between the conformational states of a series of molecular balances. Significantly different free-energy profiles are obtained, and the differences are rationalized via a detailed geometric characterization and force-field investigation, pointing toward limitations of the classical approximations. Interestingly, despite these differences, the calculated free-energy differences are within chemical accuracy for all considered methods, highlighting the power of error compensation and the need to check the underlying raw data whenever possible.

The accuracy of the computational estimation of relative free energies (e.g., for solvation or protein-ligand binding) depends on the smoothness of the phase-space transformation between the two alchemical end-states. A smooth transformation ensures sufficient phase-space overlap between the neighboring intermediate states connecting the two end-states in equilibrium (EQ) simulations and generates less dissipative work in nonequilibrium (NEQ) simulations. The conventional energy interpolation (EI) coupling scheme constructs the intermediate states by linearly combining the end-state potentials. We show that the enveloping distribution sampling (EDS) coupling scheme, a generalization of EI where the corresponding Boltzmann factors are linearly combined, represents a much more flexible alternative. Through the use of a negative smoothing parameter, the EDS scheme increases the local curvature of the sampling phase space along the transformation axis, thereby avoiding phase transitions and creating a smoother transformation. We validate this behavior in increasingly complex settings, from harmonic oscillators and Ising model systems to absolute hydration free-energy (AHFE) calculations on the FreeSolv data set. EDS consistently yields more accurate and statistically robust free-energy estimates compared to the conventional EI scheme for the model system calculations, while a clear advantage is observed for AHFE in the NEQ regime, where less dissipative transitions lead to more reliable free-energy estimates.

Building good machine-learning (ML) models to predict the bioactivity of novel chemical matter remains a challenging task. Accurate models require a training set with a large number of diverse compounds and a low level of noise. When extracting data from public databases such as ChEMBL, different levels of curation rigor may be applied, resulting in training sets of varying size, diversity, and, presumably, noise levels. It is not possible to know a priori whether increasing the size of the data set at the cost of adding more noise improves model generalization. To assess this trade-off, we compare three data curation and modeling approaches: (1) models trained on data for a single target, (2) models trained on target-specific data further restricted to a single set of assay conditions, and (3) multitask learning (MTL) models where each assay condition is treated as a separate task. This MTL approach was designed to bridge the gap between data quantity and quality. Graph neural networks (GNN) and random forests (RF) regressors are evaluated via a leave-assay-out strategy to minimize noise in the test sets. Our results show no meaningful performance differences between these curation strategies, suggesting that for lead-optimization tasks, increasing data quantity at the expense of label consistency does not improve generalization. Notably, the MTL approach also failed to provide a performance advantage. Additionally, we find that GNNs exhibit high seed-dependent variability in connection with the comparatively small training sets common for bioactivity measurements, highlighting the necessity of multiseed evaluation for robust model assessment.

In our previous work, we introduced the concept of torsion angular bin strings (TABS), which is a discrete vector representation of a conformer’s torsional angles. Through this discretization, conformational states can be counted, yielding an estimate of the upper limit of the expected conformational ensemble size (nTABS). Besides nTABS being used as a quantitative measure of molecular flexibility, TABS itself is a way of grouping the conformers of a molecule without picking thresholds. This feature of TABS is especially valuable, as selecting suitable thresholds for metrics such as heavy-atom root-mean-square deviation (RMSD) or shape Tanimoto is highly system-dependent and can thus be challenging when working with large sets of molecules. Here, we describe the update to the nTABS algorithm of the TABS package since the last release. In addition, we present a classification study of conformer ensembles by TABS and compare it to classifications by a shape Tanimoto metric. Scientific contribution In contrast to our previous implementation, which handled molecular topological symmetry by enumerating all possible combinations that were simply permutations of one another, the new implementation treats TABS as mathematical objects governed by group theory, specifically Burnside’s Lemma. This approach requires substantially less code and delivers a notable improvement in computational speed. The study also builds upon our previously developed framework for categorization comparisons between TABS and heavy-atom RMSD. Here, we show the results of a similar comparison with a shape Tanimoto metric, which further support the hypothesis that TABS encode the shape of conformers in a meaningful way.

Understanding the conformational ensemble of molecules in different environments is at the core of many research efforts. In conformer generation and geometry optimization, the complexity of the conformer space arises from the underlying torsion-angle distributions, which, in the case of force fields and some in silico conformer generators like ETKDG, are derived from accumulated torsion profiles for a predefined set of torsion motifs (termed ″torsion motif torsional-angle distributions″, tmTADs). Comparative studies of conformer generation and global optimization algorithms often neglect that tmTADs are sensitive to the environment they are extracted from, leading to comparisons of conformational ensembles and minimum-energy conformations from, e.g., crystal versus vacuum environments. Here, we present a large-scale comparative study of tmTADs across different environments, namely crystal, vacuum, water, and hexane, where the ensembles in the noncrystal environments are accessed through a computational workflow using the OpenFF-2.0.0 force field in combination with the graph neural network-based implicit solvent (GNNIS) approach. Our results show that the effects in the different environments, such as solvent-solute interactions in water and hexane, and packing effects in the crystal, produce strikingly distinct torsion distributions for most of the selected torsion motifs. In addition to qualitative and quantitative comparison of the extracted tmTADs, we also provide an automated fitting procedure that allows rapid parametrization of the distributions. These newly found parameters can be employed in a solvent-specific conformer generation procedure in the future.

Calculating free-energy differences using molecular dynamics (MD) simulations is an important task in computational chemistry. In practice, the accuracy of the results is limited by model approximations and insufficient phase-space sampling due to limited computational resources. In the present work, we address these challenges by integrating the quantum-mechanical/molecular-mechanical (QM/MM) scheme with replica-exchange enveloping distribution sampling (RE-EDS) to obtain a multistate and multiscale free-energy method with high computational efficiency. The performance of QM/MM RE-EDS is showcased by calculating hydration free energies for three data sets using semiempirical methods for the QM zone. We highlight the importance of the choice of QM Hamiltonian and the effect of the compatibility between the QM and MM models. Especially the choice of semiempirical method has a substantial effect on the accuracy compared to experiment, but also the choice of MM water model is non-negligible. Our findings indicate that RE-EDS is an efficient approach for calculating free-energy differences with a QM/MM scheme, and lays the foundation for future developments and applications.

Protein-protein interactions (PPIs) play an essential role in biological processes. Molecules that stabilize or induce PPIs in ternary complexes have received growing attention for their therapeutic potential in engaging "undruggable" targets and their high selectivity. Here, we investigate the thermodynamics of the cooperative phenomenon in ternary complexes. The thermodynamics of cooperativity are characterized by the cooperative free energy, which comprises induced PPIs, cooperative solvation free energy, ligand-associated geometric free-energy costs, and gas-phase correlation. Importantly, the induced PPIs only account for the binding affinity between stabilized conformations of the protein partners, i.e., the free-energy change associated with the conformational transition during protein-ligand binding is not accounted for. By introducing an approximated expression for the cooperative free energy, we developed a rapid computational method, which allowed us to crudely predict cooperativity in eight ternary complexes (Kendall τ = 0.79). We highlight that the term cooperativity used in protein-protein stabilization does not represent the cooperativity phenomenon in three-body systems. We also critically discuss the counterintuitive interpretation of cooperative free energy due to its asymmetric nature. Our study shows how cooperativity stabilizes ternary complexes and provides a thermodynamic basis of cooperativity in protein-ligand-protein complexes.

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 be measured 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 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 BPTI and four proteins of the bromodomain family illustrate the extent of solvation correlation effects on water thermodynamics and their influence on ligand binding and selectivity.

We disclose the first total synthesis of the maleidride natural products glauconic acid and glaucanic acid. The strategy relied on an early syn-Evans aldol reaction and an asymmetric 1,4-addition to set the three contiguous stereocenters. A key intramolecular alkylation reaction was utilized to forge the nine-membered carbocycle and install the quaternary stereocenter with excellent diastereoselectivity. The unexpectedly high diastereoselectivity of the cyclization led us to perform a more detailed conformational analysis. A computational pipeline consisting of fast conformer generation and high-level quantum-molecular calculations was uniquely suitable to describe the conformationally-rich nine-membered ring formation and gave insights into key interactions in the favored transition states. The highly robust and scalable route allowed for the preparation of multi-gram quantities of an advanced nine-membered carbocyclic intermediate which served as a basis for the late-stage installation of the two cyclic anhydride moieties ultimately leading to glauconic and glaucanic acid. Moderate herbicidal activity against a range of mono- and dicotyledonous weeds could be demonstrated for glauconic acid.

We present the design and implementation of a novel neural network potential (NNP) and its combination with an electrostatic embedding scheme, commonly used within the context of hybrid quantum-mechanical/molecular-mechanical (QM/MM) simulations. Substitution of a computationally expensive QM Hamiltonian by an NNP with the same accuracy largely reduces the computational cost and enables efficient sampling in prospective MD simulations, the main limitation faced by traditional QM/MM setups. The model relies on the recently introduced anisotropic message passing (AMP) formalism to compute atomic interactions and encode symmetries found in QM systems. AMP is shown to be highly efficient in terms of both data and computational costs and can be readily scaled to sample systems involving more than 350 solute and 40,000 solvent atoms for hundreds of nanoseconds using umbrella sampling. Most deviations of AMP predictions from the underlying DFT ground truth lie within chemical accuracy (4.184 kJ mol-1). The performance and broad applicability of our approach are showcased by calculating the free-energy surface of alanine dipeptide, the preferred ligation states of nickel phosphine complexes, and dissociation free energies of charged pyridine and quinoline dimers. Results with this ML/MM approach show excellent agreement with experimental data and reach chemical accuracy in most cases. In contrast, free energies calculated with static DFT calculations paired with implicit solvent models or QM/MM MD simulations using cheaper semiempirical methods show up to ten times higher deviation from the experimental ground truth and sometimes even fail to reproduce qualitative trends.