Effective data-driven collective variables for free energy calculations from metadynamics of paths.

A variety of enhanced sampling (ES) methods predict multidimensional free energy landscapes associated with biological and other molecular processes as a function of a few selected collective variables (CVs). The accuracy of these methods is crucially dependent on the ability of the chosen CVs to capture the relevant slow degrees of freedom of the system. For complex processes, finding such CVs is the real challenge. Machine learning (ML) CVs offer, in principle, a solution to handle this problem. However, these methods rely on the availability of high-quality datasets-ideally incorporating information about physical pathways and transition states-which are difficult to access, therefore greatly limiting their domain of application. Here, we demonstrate how these datasets can be generated by means of ES simulations in trajectory space via the metadynamics of paths algorithm. The approach is expected to provide a general and efficient way to generate efficient ML-based CVs for the fast prediction of free energy landscapes in ES simulations. We demonstrate our approach with two numerical examples, a 2D model potential and the isomerization of alanine dipeptide, using deep targeted discriminant analysis as our ML-based CV of choice.

Intramolecular and Water Mediated Tautomerism of Solvated Glycine.

Understanding tautomerism and characterizing solvent effects on the dynamic processes pose significant challenges. Using enhanced-sampling molecular dynamics based on state-of-the-art deep learning potentials, we investigated the tautomeric equilibria of glycine in water. We observed that the tautomerism between neutral and zwitterionic glycine can occur through both intramolecular and intermolecular proton transfers. The latter proceeds involving a contact anionic-glycine-hydronium ion pair or separate cationic-glycine-hydroxide ion pair. These pathways with comparable barriers contribute almost equally to the reaction flux.

Structure and polymerization of liquid sulfur across the λ-transition.

The anomalous λ-transition of liquid sulfur, which is supposed to be related to the transformation of eight-membered sulfur rings into long polymeric chains, has attracted considerable attention. However, a detailed description of the underlying dynamical polymerization process is still missing. Here, we study the structures and the mechanism of the polymerization processes of liquid sulfur across the λ-transition as well as its reverse process of formation of the rings. We do so by performing ab initio-quality molecular dynamics simulations thanks to a combination of machine learning potentials and state-of-the-art enhanced sampling techniques. With our approach, we obtain structural results that are in good agreement with the experiments and we report precious dynamical insights into the mechanisms involved in the process.

Data-driven classification of ligand unbinding pathways.

Studying the pathways of ligand-receptor binding is essential to understand the mechanism of target recognition by small molecules. The binding free energy and kinetics of protein-ligand complexes can be computed using molecular dynamics (MD) simulations, often in quantitative agreement with experiments. However, only a qualitative picture of the ligand binding/unbinding paths can be obtained through a conventional analysis of the MD trajectories. Besides, the higher degree of manual effort involved in analyzing pathways limits its applicability in large-scale drug discovery. Here, we address this limitation by introducing an automated approach for analyzing molecular transition paths with a particular focus on protein-ligand dissociation. Our method is based on the dynamic time-warping algorithm, originally designed for speech recognition. We accurately classified molecular trajectories using a very generic descriptor set of contacts or distances. Our approach outperforms manual classification by distinguishing between parallel dissociation channels, within the pathways identified by visual inspection. Most notably, we could compute exit-path-specific ligand-dissociation kinetics. The unbinding timescale along the fastest path agrees with the experimental residence time, providing a physical interpretation to our entirely data-driven protocol. In combination with appropriate enhanced sampling algorithms, this technique can be used for the initial exploration of ligand-dissociation pathways as well as for calculating path-specific thermodynamic and kinetic properties.

Cation Chloride Cotransporter NKCC1 Operates through a Rocking-Bundle Mechanism.

The sodium, potassium, and chloride cotransporter 1 (NKCC1) plays a key role in tightly regulating ion shuttling across cell membranes. Lately, its aberrant expression and function have been linked to numerous neurological disorders and cancers, making it a novel and highly promising pharmacological target for therapeutic interventions. A better understanding of how NKCC1 dynamically operates would therefore have broad implications for ongoing efforts toward its exploitation as a therapeutic target through its modulation. Based on recent structural data on NKCC1, we reveal conformational motions that are key to its function. Using extensive deep-learning-guided atomistic simulations of NKCC1 models embedded into the membrane, we captured complex dynamical transitions between alternate open conformations of the inner and outer vestibules of the cotransporter and demonstrated that NKCC1 has water-permeable states. We found that these previously undefined conformational transitions occur via a rocking-bundle mechanism characterized by the cooperative angular motion of transmembrane helices (TM) 4 and 9, with the contribution of the extracellular tip of TM 10. We found these motions to be critical in modulating ion transportation and in regulating NKCC1's water transporting capabilities. Specifically, we identified interhelical dynamical contacts between TM 10 and TM 6, which we functionally validated through mutagenesis experiments of 4 new targeted NKCC1 mutants. We conclude showing that those 4 residues are highly conserved in most Na+-dependent cation chloride cotransporters (CCCs), which highlights their critical mechanistic implications, opening the way to new strategies for NKCC1's function modulation and thus to potential drug action on selected CCCs.

The role of dynamics in heterogeneous catalysis: Surface diffusivity and N2 decomposition on Fe(111).

Dynamics has long been recognized to play an important role in heterogeneous catalytic processes. However, until recently, it has been impossible to study their dynamical behavior at industry-relevant temperatures. Using a combination of machine learning potentials and advanced simulation techniques, we investigate the cleavage of the N[Formula: see text] triple bond on the Fe(111) surface. We find that at low temperatures our results agree with the well-established picture. However, if we increase the temperature to reach operando conditions, the surface undergoes a global dynamical change and the step structure of the Fe(111) surface is destabilized. The catalytic sites, traditionally associated with this surface, appear and disappear continuously. Our simulations illuminate the danger of extrapolating low-temperature results to operando conditions and indicate that the catalytic activity can only be inferred from calculations that take dynamics fully into account. More than that, they show that it is the transition to this highly fluctuating interfacial environment that drives the catalytic process.

Free energies at QM accuracy from force fields via multimap targeted estimation.

Accurate predictions of ligand binding affinities would greatly accelerate the first stages of drug discovery campaigns. However, using highly accurate interatomic potentials based on quantum mechanics (QM) in free energy methods has been so far largely unfeasible due to their prohibitive computational cost. Here, we present an efficient method to compute QM free energies from simulations using cheap reference potentials, such as force fields (FFs). This task has traditionally been out of reach due to the slow convergence of computing the correction from the FF to the QM potential. To overcome this bottleneck, we generalize targeted free energy methods to employ multiple maps-implemented with normalizing flow neural networks (NNs)-that maximize the overlap between the distributions. Critically, the method requires neither a separate expensive training phase for the NNs nor samples from the QM potential. We further propose a one-epoch learning policy to efficiently avoid overfitting, and we combine our approach with enhanced sampling strategies to overcome the pervasive problem of poor convergence due to slow degrees of freedom. On the drug-like molecules in the HiPen dataset, the method accelerates the calculation of the free energy difference of switching from an FF to a DFTB3 potential by three orders of magnitude compared to standard free energy perturbation and by a factor of eight compared to previously published nonequilibrium calculations. Our results suggest that our method, in combination with efficient QM/MM calculations, may be used in lead optimization campaigns in drug discovery and to study protein-ligand molecular recognition processes.

Kinetics from Metadynamics: Principles, Applications, and Outlook.

Metadynamics is a popular enhanced sampling algorithm for computing the free energy landscape of rare events by using molecular dynamics simulation. Ten years ago, Tiwary and Parrinello introduced the infrequent metadynamics approach for calculating the kinetics of transitions across free energy barriers. Since then, metadynamics-based methods for obtaining rate constants have attracted significant attention in computational molecular science. Such methods have been applied to study a wide range of problems, including protein-ligand binding, protein folding, conformational transitions, chemical reactions, catalysis, and nucleation. Here, we review the principles of elucidating kinetics from metadynamics-like approaches, subsequent methodological developments in this area, and successful applications on chemical, biological, and material systems. We also highlight the challenges of reconstructing accurate kinetics from enhanced sampling simulations and the scope of future developments.

A unified framework for machine learning collective variables for enhanced sampling simulations: mlcolvar.

Identifying a reduced set of collective variables is critical for understanding atomistic simulations and accelerating them through enhanced sampling techniques. Recently, several methods have been proposed to learn these variables directly from atomistic data. Depending on the type of data available, the learning process can be framed as dimensionality reduction, classification of metastable states, or identification of slow modes. Here, we present mlcolvar, a Python library that simplifies the construction of these variables and their use in the context of enhanced sampling through a contributed interface to the PLUMED software. The library is organized modularly to facilitate the extension and cross-contamination of these methodologies. In this spirit, we developed a general multi-task learning framework in which multiple objective functions and data from different simulations can be combined to improve the collective variables. The library's versatility is demonstrated through simple examples that are prototypical of realistic scenarios.

Deep learning collective variables from transition path ensemble.

The study of the rare transitions that take place between long lived metastable states is a major challenge in molecular dynamics simulations. Many of the methods suggested to address this problem rely on the identification of the slow modes of the system, which are referred to as collective variables. Recently, machine learning methods have been used to learn the collective variables as functions of a large number of physical descriptors. Among many such methods, Deep Targeted Discriminant Analysis has proven to be useful. This collective variable is built from data harvested from short unbiased simulations in the metastable basins. Here, we enrich the set of data on which the Deep Targeted Discriminant Analysis collective variable is built by adding data from the transition path ensemble. These are collected from a number of reactive trajectories obtained using the On-the-fly Probability Enhanced Sampling flooding method. The collective variables thus trained lead to more accurate sampling and faster convergence. The performance of these new collective variables is tested on a number of representative examples.

Kinetics of Drug Release from Clay Using Enhanced Sampling Methods.

A key step in the development of a new drug, is the design of drug-excipient complexes that lead to optimal drug release kinetics. Computational chemistry and specifically enhanced sampling molecular dynamics methods can play a key role in this context, by minimizing the need for expensive experiments, and reducing cost and time. Here we show that recent advances in enhanced sampling methodologies can be brought to fruition in this area. We demonstrate the potential of these methodologies by simulating the drug release kinetics of the complex praziquantel-montmorillonite in water. Praziquantel finds promising applications in the treatment of schistosomiasis, but its biopharmaceutical profile needs to be improved, and a cheap material such as the montmorillonite clay would be a very convenient excipient. We simulate the drug release both from surface and interlayer space, and find that the diffusion of the praziquantel inside the interlayer space is the process that limits the rate of drug release.

Enhanced Sampling Aided Design of Molecular Photoswitches.

Advances in the evolving field of atomistic simulations promise important insights for the design and fundamental understanding of novel molecular photoswitches. Here, we use state-of-the-art enhanced simulation techniques to unravel the complex, multistep chemistry of donor-acceptor Stenhouse adducts (DASAs). Our reaction discovery workflow consists of enhanced sampling for efficient chemical space exploration, refinement of newly observed pathways with more accurate ab initio electronic structure calculations, and structural modifications to introduce design principles within future generations of DASAs. We showcase our discovery workflow by not only recovering the full photoswitching mechanism of DASA but also predicting a plethora of new plausible thermal pathways and suggesting a way for their experimental validation. Furthermore, we illustrate the tunability of these newly discovered reactions, leading to a potential avenue for controlling DASA dynamics through multiple external stimuli. Overall, these insights could offer alternative routes to increase the efficiency and control of DASA's photoswitching mechanism, providing new elements to design more complex light-responsive materials.

Rare Event Kinetics from Adaptive Bias Enhanced Sampling.

We introduce a novel enhanced sampling approach named on-the-fly probability enhanced sampling (OPES) flooding for calculating the kinetics of rare events from atomistic molecular dynamics simulation. This method is derived from the OPES approach [Invernizzi and Parrinello, J. Phys. Chem. Lett. 2020, 11, 7, 2731-2736], which has been recently developed for calculating converged free energy surfaces for complex systems. In this paper, we describe the theoretical details of the OPES flooding technique and demonstrate the application on three systems of increasing complexity: barrier crossing in a two-dimensional double-well potential, conformational transition in the alanine dipeptide in the gas phase, and the folding and unfolding of the chignolin polypeptide in an aqueous environment. From extensive tests, we show that the calculation of accurate kinetics not only requires the transition state to be bias-free, but the amount of bias deposited should also not exceed the effective barrier height measured along the chosen collective variables. In this vein, the possibility of computing rates from biasing suboptimal order parameters has also been explored. Furthermore, we describe the choice of optimum parameter combinations for obtaining accurate results from limited computational effort.

Water regulates the residence time of Benzamidine in Trypsin.

The process of ligand-protein unbinding is crucial in biophysics. Water is an essential part of any biological system and yet, many aspects of its role remain elusive. Here, we simulate with state-of-the-art enhanced sampling techniques the binding of Benzamidine to Trypsin which is a much studied and paradigmatic ligand-protein system. We use machine learning methods to determine efficient collective coordinates for the complex non-local network of water. These coordinates are used to perform On-the-fly Probability Enhanced Sampling simulations, which we adapt to calculate also the ligand residence time. Our results, both static and dynamic, are in good agreement with experiments. We find that the presence of a water molecule located at the bottom of the binding pocket allows via a network of hydrogen bonds the ligand to be released into the solution. On a finer scale, even when unbinding is allowed, another water molecule further modulates the exit time.

Characterizing Metastable States with the Help of Machine Learning.

Present-day atomistic simulations generate long trajectories of ever more complex systems. Analyzing these data, discovering metastable states, and uncovering their nature are becoming increasingly challenging. In this paper, we first use the variational approach to conformation dynamics to discover the slowest dynamical modes of the simulations. This allows the different metastable states of the system to be located and organized hierarchically. The physical descriptors that characterize metastable states are discovered by means of a machine learning method. We show in the cases of two proteins, chignolin and bovine pancreatic trypsin inhibitor, how such analysis can be effortlessly performed in a matter of seconds. Another strength of our approach is that it can be applied to the analysis of both unbiased and biased simulations.

Solubility of Organic Salts in Solvent-Antisolvent Mixtures: A Combined Experimental and Molecular Dynamics Simulations Approach.

We combine molecular dynamics simulations with experiments to estimate solubilities of an organic salt in complex growth environments. We predict the solubility by simulations of the growth and dissolution of ions at the crystal surface kink sites at different solution concentrations. Thereby, the solubility is identified as the solution's salt concentration, where the energy of the ion pair dissolved in solution equals the energy of the ion pair crystallized at the kink sites. The simulation methodology is demonstrated for the case of anhydrous sodium acetate crystallized from various solvent-antisolvent mixtures. To validate the predicted solubilities, we have measured the solubilities of sodium acetate in-house, using an experimental setup and measurement protocol that guarantees moisture-free conditions, which is key for a hygroscopic compound like sodium acetate. We observe excellent agreement between the experimental and the computationally evaluated solubilities for sodium acetate in different solvent-antisolvent mixtures. Given the agreement and the rich data the simulations produce, we can use them to complement experimental tasks, which in turn will reduce time and capital in the design of complicated industrial crystallization processes of organic salts.

Exploration vs Convergence Speed in Adaptive-Bias Enhanced Sampling.

In adaptive-bias enhanced sampling methods, a bias potential is added to the system to drive transitions between metastable states. The bias potential is a function of a few collective variables and is gradually modified according to the underlying free energy surface. We show that when the collective variables are suboptimal, there is an exploration-convergence tradeoff, and one must choose between a quickly converging bias that will lead to fewer transitions or a slower to converge bias that can explore the phase space more efficiently but might require a much longer time to produce an accurate free energy estimate. The recently proposed on-the-fly probability enhanced sampling (OPES) method focuses on fast convergence, but there are cases where fast exploration is preferred instead. For this reason, we introduce a new variant of the OPES method that focuses on quickly escaping metastable states at the expense of convergence speed. We illustrate the benefits of this approach in prototypical systems and show that it outperforms the popular metadynamics method.

Computing Rates and Understanding Unbinding Mechanisms in Host-Guest Systems.

The long time scale associated with ligand residence times renders their computation challenging. Therefore, the influence of factors like solvation and steric hindrance on residence times is not fully understood. Here, we demonstrate in a set of model host-guest systems that the recently developed Gaussian mixture based enhanced sampling allows residence times to be computed and enables an understanding of their unbinding mechanism. We observe that guest unbinding often proceeds via a series of intermediate states that can be labeled by the number of water molecules present in the binding cavity. In several cases the residence time is correlated to the water trapping times in the cavity.

Discover, Sample, and Refine: Exploring Chemistry with Enhanced Sampling Techniques.

Over the last few decades, enhanced sampling methods have been continuously improved. Here, we exploit this progress and propose a modular workflow for blind reaction discovery and determination of reaction paths. In a three-step strategy, at first we use a collective variable derived from spectral graph theory in conjunction with the explore variant of the on-the-fly probability enhanced sampling method to drive reaction discovery runs. Once different chemical products are determined, we construct an ad-hoc neural network-based collective variable to improve sampling, and finally we refine the results using the free energy perturbation theory and a more accurate Hamiltonian. We apply this strategy to both intramolecular and intermolecular reactions. Our workflow requires minimal user input and extends the power of ab initio molecular dynamics to explore and characterize the reaction space.

Prediction of a Supersolid Phase in High-Pressure Deuterium.

Supersolid is a mysterious and puzzling state of matter whose possible existence has stirred a vigorous debate among physicists for over 60 years. Its elusive nature stems from the coexistence of two seemingly contradicting properties, long-range order and superfluidity. We report computational evidence of a supersolid phase of deuterium under high pressure (p>800 GPa) and low temperature (T<1.0 K). In our simulations, that are based on bosonic path integral molecular dynamics, we observe a highly concerted exchange of atoms while the system preserves its crystalline order. The exchange processes are favored by the soft core interactions between deuterium atoms that form a densely packed metallic solid. At the zero temperature limit, Bose-Einstein condensation is observed as the permutation probability of N deuterium atoms approaches 1/N with a finite superfluid fraction. Our study provides concrete evidence for the existence of a supersolid phase in high-pressure deuterium and could provide insights on the future investigation of supersolid phases in real materials.