How Robust Is the Ligand Binding Transition State?

For many drug targets, it has been shown that the kinetics of drug binding (e.g., on rate and off rate) is more predictive of drug efficacy than thermodynamic quantities alone. This motivates the development of predictive computational models that can be used to optimize compounds on the basis of their kinetics. The structural details underpinning these computational models are found not only in the bound state but also in the short-lived ligand binding transition states. Although transition states cannot be directly observed experimentally due to their extremely short lifetimes, recent successes have demonstrated that modeling the ligand binding transition state is possible with the help of enhanced sampling molecular dynamics methods. Previously, we generated unbinding paths for an inhibitor of soluble epoxide hydrolase (sEH) with a residence time of 11 min. Here, we computationally modeled unbinding events with the weighted ensemble method REVO (resampling of ensembles by variation optimization) for five additional inhibitors of sEH with residence times ranging from 14.25 to 31.75 min, with average prediction accuracy within an order of magnitude. The unbinding ensembles are analyzed in detail, focusing on features of the ligand binding transition state ensembles (TSEs). We find that ligands with similar bound poses can show significant differences in their ligand binding TSEs, in terms of their spatial distribution and protein-ligand interactions. However, we also find similarities across the TSEs when examining more general features such as ligand degrees of freedom. Together these findings show significant challenges for rational, kinetics-based drug design.

Flexible Topology: A Dynamic Model of a Continuous Chemical Space.

Ligand design problems involve searching chemical space for a molecule with a set of desired properties. As chemical space is discrete, this search must be conducted in a pointwise manner, separately investigating one molecule at a time, which can be inefficient. We propose a method called "Flexible Topology", where a ligand is composed of a set of shapeshifting "ghost" atoms, whose atomic identities and connectivity can dynamically change over the course of a simulation. Ghost atoms are guided toward their target positions using a translation-, rotation-, and index-invariant restraint potential. This is the first step toward a continuous model of chemical space, where a dynamic simulation can move from one molecule to another by following gradients of a potential energy function. This builds on a substantial history of alchemy in the field of molecular dynamics simulation, including the Lambda dynamics method developed by Brooks and co-workers [X. Kong and C.L. Brooks III, J. Chem. Phys. 105, 2414 (1996)], but takes it to an extreme by associating a set of four dynamical attributes with each shapeshifting ghost atom that control not only its presence but also its atomic identity. Here, we outline the theoretical details of this method, its implementation using the OpenMM simulation package, and some preliminary studies of ghost particle assembly simulations in vacuum. We examine a set of 10 small molecules, ranging in size from 6 to 50 atoms, and show that Flexible Topology is able to consistently assemble all of these molecules to high accuracy, beginning from randomly initialized positions and attributes.

Adrenomedullin 2/intermedin is a slow off-rate, long-acting endogenous agonist of the adrenomedullin2 G protein-coupled receptor.

Adrenomedullin 2/intermedin (AM2/IMD), adrenomedullin (AM), and calcitonin gene-related peptide (CGRP) have functions in the cardiovascular, lymphatic, and nervous systems by activating three heterodimeric receptors comprising the class B GPCR CLR and a RAMP1, -2, or -3 modulatory subunit. CGRP and AM prefer the RAMP1 and RAMP2/3 complexes, respectively, whereas AM2/IMD is thought to be relatively nonselective. Accordingly, AM2/IMD exhibits overlapping actions with CGRP and AM, so the rationale for this third agonist for the CLR-RAMP complexes is unclear. Here, we report that AM2/IMD is kinetically selective for CLR-RAMP3, known as the AM2R, and we define the structural basis for its distinct kinetics. In live cell biosensor assays, AM2/IMD-AM2R elicited longer-duration cAMP signaling than the other peptide-receptor combinations. AM2/IMD and AM bound the AM2R with similar equilibrium affinities, but AM2/IMD had a slower off-rate and longer receptor residence time, thus explaining its prolonged signaling capacity. Peptide and receptor chimeras and mutagenesis were used to map the regions responsible for the distinct binding and signaling kinetics to the AM2/IMD mid-region and the RAMP3 extracellular domain (ECD). Molecular dynamics simulations revealed how the former forms stable interactions at the CLR ECD-transmembrane domain interface and how the latter augments the CLR ECD binding pocket to anchor the AM2/IMD C terminus. These strong binding components only combine in the AM2R. Our findings uncover AM2/IMD-AM2R as a cognate pair with unique temporal features, reveal how AM2/IMD and RAMP3 collaborate to shape CLR signaling, and have significant implications for AM2/IMD biology.

Adrenomedullin 2/intermedin is a slow off-rate, long-acting endogenous agonist of the adrenomedullin 2 G protein-coupled receptor.

The signaling peptides adrenomedullin 2/intermedin (AM2/IMD), adrenomedullin (AM), and CGRP have overlapping and distinct functions in the cardiovascular, lymphatic, and nervous systems by activating three shared receptors comprised of the class B GPCR CLR in complex with a RAMP1, -2, or -3 modulatory subunit. Here, we report that AM2/IMD, which is thought to be a non-selective agonist, is kinetically selective for CLR-RAMP3, known as the AM 2 R. AM2/IMD-AM 2 R elicited substantially longer duration cAMP signaling than the eight other peptide-receptor combinations due to AM2/IMD slow off-rate binding kinetics. The regions responsible for the slow off-rate were mapped to the AM2/IMD mid-region and the RAMP3 extracellular domain. MD simulations revealed how these bestow enhanced stability to the complex. Our results uncover AM2/IMD-AM 2 R as a cognate pair with unique temporal features, define the mechanism of kinetic selectivity, and explain how AM2/IMD and RAMP3 collaborate to shape the signaling output of a clinically important GPCR.

The Small Molecule Antagonist KCI807 Disrupts Association of the Amino-Terminal Domain of the Androgen Receptor with ELK1 by Modulating the Adjacent DNA Binding Domain.

The androgen receptor (AR) is a crucial coactivator of ELK1 for prostate cancer (PCa) growth, associating with ELK1 through two peptide segments (358-457 and 514-557) within the amino-terminal domain (NTD) of AR. The small-molecule antagonist 5-hydroxy-2-(3-hydroxyphenyl)chromen-4-one (KCI807) binds to AR, blocking ELK1 binding and inhibiting PCa growth. We investigated the mode of interaction of KCI807 with AR using systematic mutagenesis coupled with ELK1 coactivation assays, testing polypeptide binding and Raman spectroscopy. In full-length AR, deletion of neither ELK1 binding segment affected sensitivity of residual ELK1 coactivation to KCI807. Although the NTD is sufficient for association of AR with ELK1, interaction of the isolated NTD with ELK1 was insensitive to KCI807. In contrast, coactivation of ELK1 by the AR-V7 splice variant, comprising the NTD and the DNA binding domain (DBD), was sensitive to KCI807. Deletions and point mutations within DBD segment 558-595, adjacent to the NTD, interfered with coactivation of ELK1, and residual ELK1 coactivation by the mutants was insensitive to KCI807. In a glutathione S-transferase pull-down assay, KCI807 inhibited ELK1 binding to an AR polypeptide that included the two ELK1 binding segments and the DBD but did not affect ELK1 binding to a similar AR segment that lacked the sequence downstream of residue 566. Raman spectroscopy detected KCI807-induced conformational change in the DBD. The data point to a putative KCI807 binding pocket within the crystal structure of the DBD and indicate that either mutations or binding of KCI807 at this site will induce conformational changes that disrupt ELK1 binding to the NTD. SIGNIFICANCE STATEMENT: The small-molecule antagonist KCI807 disrupts association of the androgen receptor (AR) with ELK1, serving as a prototype for the development of small molecules for a novel type of therapeutic intervention in drug-resistant prostate cancer. This study provides basic information needed for rational KCI807-based drug design by identifying a putative binding pocket in the DNA binding domain of AR through which KCI807 modulates the amino-terminal domain to inhibit ELK1 binding.

Quality over quantity: Sampling high probability rare events with the weighted ensemble algorithm.

The prediction of (un)binding rates and free energies is of great significance to the drug design process. Although many enhanced sampling algorithms and approaches have been developed, there is not yet a reliable workflow to predict these quantities. Previously we have shown that free energies and transition rates can be calculated by directly simulating the binding and unbinding processes with our variant of the WE algorithm "Resampling of Ensembles by Variation Optimization", or "REVO". Here, we calculate binding free energies retrospectively for three SAMPL6 host-guest systems and prospectively for a SAMPL9 system to test a modification of REVO that restricts its cloning behavior in quasi-unbound states. Specifically, trajectories cannot clone if they meet a physical requirement that represents a high likelihood of unbinding, which in the case of this work is a center-of-mass to center-of-mass distance. The overall effect of this change was difficult to predict, as it results in fewer unbinding events each of which with a much higher statistical weight. For all four systems tested, this new strategy produced either more accurate unbinding free energies or more consistent results between simulations than the standard REVO algorithm. This approach is highly flexible, and any feature of interest for a system can be used to determine cloning eligibility. These findings thus constitute an important improvement in the calculation of transition rates and binding free energies with the weighted ensemble method.

Predicting the structural basis of targeted protein degradation by integrating molecular dynamics simulations with structural mass spectrometry.

Targeted protein degradation (TPD) is a promising approach in drug discovery for degrading proteins implicated in diseases. A key step in this process is the formation of a ternary complex where a heterobifunctional molecule induces proximity of an E3 ligase to a protein of interest (POI), thus facilitating ubiquitin transfer to the POI. In this work, we characterize 3 steps in the TPD process. (1) We simulate the ternary complex formation of SMARCA2 bromodomain and VHL E3 ligase by combining hydrogen-deuterium exchange mass spectrometry with weighted ensemble molecular dynamics (MD). (2) We characterize the conformational heterogeneity of the ternary complex using Hamiltonian replica exchange simulations and small-angle X-ray scattering. (3) We assess the ubiquitination of the POI in the context of the full Cullin-RING Ligase, confirming experimental ubiquitinomics results. Differences in degradation efficiency can be explained by the proximity of lysine residues on the POI relative to ubiquitin.

A biosensor of protein foldedness identifies increased "holdase" activity of chaperones in the nucleus following increased cytosolic protein aggregation.

Chaperones and other quality control machinery guard proteins from inappropriate aggregation, which is a hallmark of neurodegenerative diseases. However, how the systems that regulate the "foldedness" of the proteome remain buffered under stress conditions and in different cellular compartments remains incompletely understood. In this study, we applied a FRET-based strategy to explore how well quality control machinery protects against the misfolding and aggregation of "bait" biosensor proteins, made from the prokaryotic ribonuclease barnase, in the nucleus and cytosol of human embryonic kidney 293T cells. We found that those barnase biosensors were prone to misfolding, were less engaged by quality control machinery, and more prone to inappropriate aggregation in the nucleus as compared with the cytosol, and that these effects could be regulated by chaperone Hsp70-related machinery. Furthermore, aggregation of mutant huntingtin exon 1 protein (Httex1) in the cytosol appeared to outcompete and thus prevented the engagement of quality control machinery with the biosensor in the cytosol. This effect correlated with reduced levels of DNAJB1 and HSPA1A chaperones in the cell outside those sequestered to the aggregates, particularly in the nucleus. Unexpectedly, we found Httex1 aggregation also increased the apparent engagement of the barnase biosensor with quality control machinery in the nucleus suggesting an independent implementation of "holdase" activity of chaperones other than DNAJB1 and HSPA1A. Collectively, these results suggest that proteostasis stress can trigger a rebalancing of chaperone abundance in different subcellular compartments through a dynamic network involving different chaperone-client interactions.

Local Ion Densities can Influence Transition Paths of Molecular Binding.

Improper reaction coordinates can pose significant problems for path-based binding free energy calculations. Particularly, omission of long timescale motions can lead to over-estimation of the energetic barriers between the bound and unbound states. Many methods exist to construct the optimal reaction coordinate using a pre-defined basis set of features. Although simulations are typically conducted in explicit solvent, the solvent atoms are often excluded by these feature sets-resulting in little being known about their role in reaction coordinates, and ultimately, their role in determining (un)binding rates and free energies. In this work, analysis is done on an extensive set of host-guest unbinding trajectories, working to characterize differences between high and low probability unbinding trajectories with a focus on solvent-based features, including host-ion interactions, guest-ion interactions and location-dependent ion densities. We find that differences in ion densities as well as guest-ion interactions strongly correlate with differences in the probabilities of reactive paths that are used to determine free energies of (un)binding and play a significant role in the unbinding process.

Creating Maps of the Ligand Binding Landscape for Kinetics-Based Drug Discovery.

Simulations of ligand-protein interactions can be very useful for drug design and to gain biological insight. Full pathways of ligand-protein binding can be used to get information about ligand binding transition states, which form the rate-limiting step of the binding and release processes. However, these simulations are typically limited by the presence of large energy barriers that separate stable poses of interest. Here we describe a simulation protocol for exploring and analyzing landscapes of ligand-protein interactions that makes use of molecular docking, enhanced molecular simulation with the weighted ensemble algorithm, and network analysis. It can be accomplished using a modest cluster of graphics processing units and freely accessible software. This protocol focuses on the construction and analysis of a network model of ligand binding poses and provides links to resources that describe the other steps in more detail. The end result of this protocol is a map of the ligand-protein binding landscape that identifies transition states of the ligand binding pathway, as well as alternative bound poses that could be stabilized with modifications to the ligand.

Identification and Characterization of Two Structurally Related Dipeptides that Enhance Catalytic Efficiency of Neurolysin.

Neurolysin (Nln) is a recently recognized endogenous mechanism functioning to preserve the brain from ischemic injury. To further understand the pathophysiological function of this peptidase in stroke and other neurologic disorders, the present study was designed to identify small molecule activators of Nln. Using a computational approach, the structure of Nln was explored, which was followed by docking and in silico screening of ∼140,000 molecules from the National Cancer Institute Developmental Therapeutics Program database. Top ranking compounds were evaluated in an Nln enzymatic assay, and two hit histidine-dipeptides were further studied in detail. The identified dipeptides enhanced the rate of synthetic substrate hydrolysis by recombinant (human and rat) and mouse brain-purified Nln in a concentration-dependent manner (micromolar A50 and Amax ≥ 300%) but had negligible effect on activity of closely related peptidases. Both dipeptides also enhanced hydrolysis of Nln endogenous substrates neurotensin, angiotensin I, and bradykinin and increased efficiency of the synthetic substrate hydrolysis (Vmax/Km ratio) in a concentration-dependent manner. The dipeptides and competitive inhibitor dynorphin A (1-13) did not affect each other's affinity for Nln, suggesting differing nature of their respective binding sites. Lastly, drug affinity responsive target stability (DARTS) and differential scanning fluorimetry (DSF) assays confirmed concentration-dependent interaction of Nln with the activator molecule. This is the first study demonstrating that Nln activity can be enhanced by small molecules, although the peptidic nature and low potency of the activators limit their application. The identified dipeptides provide a chemical scaffold to develop high-potency, drug-like molecules as research tools and potential drug leads. SIGNIFICANCE STATEMENT: This study describes discovery of two molecules that selectively enhance activity of peptidase Nln-a newly recognized cerebroprotective mechanism in the poststroke brain. The identified molecules will serve as a chemical scaffold for development of drug-like molecules to further study Nln and may become lead structures for a new class of drugs. In addition, our conceptual and methodological framework and research findings might be used for other peptidases and enzymes, the activation of which bears therapeutic potential.

Perturbation of ACE2 Structural Ensembles by SARS-CoV-2 Spike Protein Binding.

The human ACE2 enzyme serves as a critical first recognition point of coronaviruses, including SARS-CoV-2. In particular, the extracellular domain of ACE2 interacts directly with the S1 tailspike protein of the SARS-CoV-2 virion through a broad protein-protein interface. Although this interaction has been characterized by X-ray crystallography, these structures do not reveal significant differences in the ACE2 structure upon S1 protein binding. In this work, using several all-atom molecular dynamics simulations, we show persistent differences in the ACE2 structure upon binding. These differences are determined with the linear discriminant analysis (LDA) machine learning method and validated using independent training and testing datasets, including long trajectories generated by D. E. Shaw Research on the Anton 2 supercomputer. In addition, long trajectories for 78 potent ACE2-binding compounds, also generated by D. E. Shaw Research, were projected onto the LDA classification vector in order to determine whether the ligand-bound ACE2 structures were compatible with S1 protein binding. This allows us to predict which compounds are "apo-like" versus "complex-like" and to pinpoint long-range ligand-induced allosteric changes in the ACE2 structure.

Predicting partition coefficients for the SAMPL7 physical property challenge using the ClassicalGSG method.

The prediction of [Formula: see text] values is one part of the statistical assessment of the modeling of proteins and ligands (SAMPL) blind challenges. Here, we use a molecular graph representation method called Geometric Scattering for Graphs (GSG) to transform atomic attributes to molecular features. The atomic attributes used here are parameters from classical molecular force fields including partial charges and Lennard-Jones interaction parameters. The molecular features from GSG are used as inputs to neural networks that are trained using a "master" dataset comprised of over 41,000 unique [Formula: see text] values. The specific molecular targets in the SAMPL7 [Formula: see text] prediction challenge were unique in that they all contained a sulfonyl moeity. This motivated a set of ClassicalGSG submissions where predictors were trained on different subsets of the master dataset that are filtered according to chemical types and/or the presence of the sulfonyl moeity. We find that our ranked prediction obtained 5th place with an RMSE of 0.77 [Formula: see text] units and an MAE of 0.62, while one of our non-ranked predictions achieved first place among all submissions with an RMSE of 0.55 and an MAE of 0.44. After the conclusion of the challenge we also examined the performance of open-source force field parameters that allow for an end-to-end [Formula: see text] predictor model: General AMBER Force Field (GAFF), Universal Force Field (UFF), Merck Molecular Force Field 94 (MMFF94) and Ghemical. We find that ClassicalGSG models trained with atomic attributes from MMFF94 can yield more accurate predictions compared to those trained with CGenFF atomic attributes.

ClassicalGSG: Prediction of log P using classical molecular force fields and geometric scattering for graphs.

This work examines methods for predicting the partition coefficient (log P) for a dataset of small molecules. Here, we use atomic attributes such as radius and partial charge, which are typically used as force field parameters in classical molecular dynamics simulations. These atomic attributes are transformed into index-invariant molecular features using a recently developed method called geometric scattering for graphs (GSG). We call this approach "ClassicalGSG" and examine its performance under a broad range of conditions and hyperparameters. We train ClassicalGSG log P predictors with neural networks using 10,722 molecules from the OpenChem dataset and apply them to predict the log P values from four independent test sets. The ClassicalGSG method's performance is compared to a baseline model that employs graph convolutional networks. Our results show that the best prediction accuracies are obtained using atomic attributes generated with the CHARMM generalized force field and 2D molecular structures.

Membrane-Mediated Ligand Unbinding of the PK-11195 Ligand from TSPO.

The translocator protein (TSPO), previously known as the peripheral benzodiazepine receptor, is of longstanding medical interest as both a biomarker for neuroinjury and a potential drug target for neuroinflammation and other disorders. Recently, it was shown that ligand residence time is a key factor determining steroidogenic efficacy of TSPO-binding compounds. This spurs interest in simulations of (un)binding pathways of TSPO ligands, which could reveal the molecular interactions governing ligand residence time. In this study, we use a weighted ensemble algorithm to determine the unbinding pathway for different poses of PK-11195, a TSPO ligand used in neuroimaging. In contrast with previous studies, our results show that PK-11195 does not dissociate directly into the solvent but instead dissociates via the lipid membrane by going between the transmembrane helices. We analyze this path ensemble in detail, constructing descriptors that can facilitate a general understanding of membrane-mediated ligand binding. We construct a set of Markov state models augmented with additional straightforward simulations to determine pose-specific ligand residence times. Together, we combine over 40 µs of trajectory data to form a coherent picture of the ligand binding landscape. We find that multiple starting poses yield residence times that roughly agree with the experimental quantity. The ligand binding transition states predicted by these Markov state models occur when PK-11195 is already in the membrane and involves only minimal ligand-protein interactions. This has implications for the design of new long-residence-time TSPO ligands.

Wepy: A Flexible Software Framework for Simulating Rare Events with Weighted Ensemble Resampling.

Here, we introduce the open-source software framework wepy (https://github.com/ADicksonLab/wepy) which is a toolkit for running and analyzing weighted ensemble (WE) simulations. The wepy toolkit is in pure Python and as such is highly portable and extensible, making it an excellent platform to develop and use new WE resampling algorithms such as WExplore, REVO, and others while leveraging the entire Python ecosystem. In addition, wepy simplifies WE-specific analyses by defining out-of-core tree-like data structures using the cross-platform HDF5 file format. In this paper, we discuss the motivations and challenges for simulating rare events in biomolecular systems. As has previously been shown, high-dimensional WE resampling algorithms such as WExplore and REVO have been successful at these tasks, especially for rare events that are difficult to describe by one or two collective variables. We explain in detail how wepy facilitates implementation of these algorithms, as well as aids in analyzing the unique structure of WE simulation results. To explain how wepy and WE work in general, we describe the mathematical formalism of WE, an overview of the architecture of wepy, and provide code examples of how to construct, run, and analyze simulation results for a protein-ligand system (T4 Lysozyme in an implicit solvent). This paper is written with a variety of readers in mind, including (1) those curious about how to leverage WE rare-event simulations for their domain, (2) current WE users who want to begin using new high-dimensional resamplers such as WExplore and REVO, and (3) expert users who would like to prototype or implement their own algorithms that can be easily adopted by others.

Enhanced Jarzynski free energy calculations using weighted ensemble.

The free energy of transitions between stable states is the key thermodynamic quantity that governs the relative probabilities of the forward and reverse reactions and the ratio of state probabilities at equilibrium. The binding free energy of a drug and its receptor is of particular interest, as it serves as an optimization function for drug design. Over the years, many computational methods have been developed to calculate binding free energies, and while many of these methods have a long history, issues such as convergence of free energy estimates and the projection of a binding process onto order parameters remain. Over 20 years ago, the Jarzynski equality was derived with the promise to calculate equilibrium free energies by measuring the work applied to short nonequilibrium trajectories. However, these calculations were found to be dominated by trajectories with low applied work that occur with extremely low probability. Here, we examine the combination of weighted ensemble algorithms with the Jarzynski equality. In this combined method, an ensemble of nonequilibrium trajectories are run in parallel, and cloning and merging operations are used to preferentially sample low-work trajectories that dominate the free energy calculations. Two additional methods are also examined: (i) a novel weighted ensemble resampler that samples trajectories directly according to their importance to the work of work and (ii) the diffusion Monte Carlo method using the applied work as the selection potential. We thoroughly examine both the accuracy and efficiency of unbinding free energy calculations for a series of model Lennard-Jones atom pairs with interaction strengths ranging from 2 kcal/mol to 20 kcal/mol. We find that weighted ensemble calculations can more efficiently determine accurate binding free energies, especially for deeper Lennard-Jones well depths.

On Calculating Free Energy Differences Using Ensembles of Transition Paths.

The free energy of a process is the fundamental quantity that determines its spontaneity or propensity at a given temperature. In particular, the binding free energy of a drug candidate to its biomolecular target is used as an objective quantity in drug design. Recently, binding kinetics-rates of association (k on) and dissociation (k off)-have also demonstrated utility for their ability to predict efficacy and in some cases have been shown to be more predictive than the binding free energy alone. Some methods exist to calculate binding kinetics from molecular simulations, although these are typically more difficult to calculate than the binding affinity as they depend on details of the transition path ensemble. Assessing these rate constants can be difficult, due to uncertainty in the definition of the bound and unbound states, large error bars and the lack of experimental data. As an additional consistency check, rate constants from simulation can be used to calculate free energies (using the log of their ratio) which can then be compared to free energies obtained experimentally or using alchemical free energy perturbation. However, in this calculation it is not straightforward to account for common, practical details such as the finite simulation volume or the particular definition of the "bound" and "unbound" states. Here we derive a set of correction terms that can be applied to calculations of binding free energies using full reactive trajectories. We apply these correction terms to revisit the calculation of binding free energies from rate constants for a host-guest system that was part of a blind prediction challenge, where significant deviations were observed between free energies calculated with rate ratios and those calculated from alchemical perturbation. The correction terms combine to significantly decrease the error with respect to computational benchmarks, from 3.4 to 0.76 kcal/mol. Although these terms were derived with weighted ensemble simulations in mind, some of the correction terms are generally applicable to free energies calculated using physical pathways via methods such as Markov state modeling, metadynamics, milestoning, or umbrella sampling.

The SAMPL6 SAMPLing challenge: assessing the reliability and efficiency of binding free energy calculations.

Approaches for computing small molecule binding free energies based on molecular simulations are now regularly being employed by academic and industry practitioners to study receptor-ligand systems and prioritize the synthesis of small molecules for ligand design. Given the variety of methods and implementations available, it is natural to ask how the convergence rates and final predictions of these methods compare. In this study, we describe the concept and results for the SAMPL6 SAMPLing challenge, the first challenge from the SAMPL series focusing on the assessment of convergence properties and reproducibility of binding free energy methodologies. We provided parameter files, partial charges, and multiple initial geometries for two octa-acid (OA) and one cucurbit[8]uril (CB8) host-guest systems. Participants submitted binding free energy predictions as a function of the number of force and energy evaluations for seven different alchemical and physical-pathway (i.e., potential of mean force and weighted ensemble of trajectories) methodologies implemented with the GROMACS, AMBER, NAMD, or OpenMM simulation engines. To rank the methods, we developed an efficiency statistic based on bias and variance of the free energy estimates. For the two small OA binders, the free energy estimates computed with alchemical and potential of mean force approaches show relatively similar variance and bias as a function of the number of energy/force evaluations, with the attach-pull-release (APR), GROMACS expanded ensemble, and NAMD double decoupling submissions obtaining the greatest efficiency. The differences between the methods increase when analyzing the CB8-quinine system, where both the guest size and correlation times for system dynamics are greater. For this system, nonequilibrium switching (GROMACS/NS-DS/SB) obtained the overall highest efficiency. Surprisingly, the results suggest that specifying force field parameters and partial charges is insufficient to generally ensure reproducibility, and we observe differences between seemingly converged predictions ranging approximately from 0.3 to 1.0 kcal/mol, even with almost identical simulations parameters and system setup (e.g., Lennard-Jones cutoff, ionic composition). Further work will be required to completely identify the exact source of these discrepancies. Among the conclusions emerging from the data, we found that Hamiltonian replica exchange-while displaying very small variance-can be affected by a slowly-decaying bias that depends on the initial population of the replicas, that bidirectional estimators are significantly more efficient than unidirectional estimators for nonequilibrium free energy calculations for systems considered, and that the Berendsen barostat introduces non-negligible artifacts in expanded ensemble simulations.

Consistently dated Atlantic sediment cores over the last 40 thousand years.

Rapid changes in ocean circulation and climate have been observed in marine-sediment and ice cores over the last glacial period and deglaciation, highlighting the non-linear character of the climate system and underlining the possibility of rapid climate shifts in response to anthropogenic greenhouse gas forcing. To date, these rapid changes in climate and ocean circulation are still not fully explained. One obstacle hindering progress in our understanding of the interactions between past ocean circulation and climate changes is the difficulty of accurately dating marine cores. Here, we present a set of 92 marine sediment cores from the Atlantic Ocean for which we have established age-depth models that are consistent with the Greenland GICC05 ice core chronology, and computed the associated dating uncertainties, using a new deposition modeling technique. This is the first set of consistently dated marine sediment cores enabling paleoclimate scientists to evaluate leads/lags between circulation and climate changes over vast regions of the Atlantic Ocean. Moreover, this data set is of direct use in paleoclimate modeling studies.