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.

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.

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.

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.

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.

Deep learning the slow modes for rare events sampling.

The development of enhanced sampling methods has greatly extended the scope of atomistic simulations, allowing long-time phenomena to be studied with accessible computational resources. Many such methods rely on the identification of an appropriate set of collective variables. These are meant to describe the system's modes that most slowly approach equilibrium under the action of the sampling algorithm. Once identified, the equilibration of these modes is accelerated by the enhanced sampling method of choice. An attractive way of determining the collective variables is to relate them to the eigenfunctions and eigenvalues of the transfer operator. Unfortunately, this requires knowing the long-term dynamics of the system beforehand, which is generally not available. However, we have recently shown that it is indeed possible to determine efficient collective variables starting from biased simulations. In this paper, we bring the power of machine learning and the efficiency of the recently developed on the fly probability-enhanced sampling method to bear on this approach. The result is a powerful and robust algorithm that, given an initial enhanced sampling simulation performed with trial collective variables or generalized ensembles, extracts transfer operator eigenfunctions using a neural network ansatz and then accelerates them to promote sampling of rare events. To illustrate the generality of this approach, we apply it to several systems, ranging from the conformational transition of a small molecule to the folding of a miniprotein and the study of materials crystallization.

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.

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.

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.

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.

Integrating NMR and simulations reveals motions in the UUCG tetraloop.

We provide an atomic-level description of the structure and dynamics of the UUCG RNA stem-loop by combining molecular dynamics simulations with experimental data. The integration of simulations with exact nuclear Overhauser enhancements data allowed us to characterize two distinct states of this molecule. The most stable conformation corresponds to the consensus three-dimensional structure. The second state is characterized by the absence of the peculiar non-Watson-Crick interactions in the loop region. By using machine learning techniques we identify a set of experimental measurements that are most sensitive to the presence of non-native states. We find that although our MD ensemble, as well as the consensus UUCG tetraloop structures, are in good agreement with experiments, there are remaining discrepancies. Together, our results show that (i) the MD simulation overstabilize a non-native loop conformation, (ii) eNOE data support its presence with a population of ≈10% and (iii) the structural interpretation of experimental data for dynamic RNAs is highly complex, even for a simple model system such as the UUCG tetraloop.

Path integral molecular dynamics for bosons.

Trapped bosons exhibit fundamental physical phenomena and are at the core of emerging quantum technologies. We present a method for simulating bosons using path integral molecular dynamics. The main difficulty in performing such simulations is enumerating all ring-polymer configurations, which arise due to permutations of identical particles. We show that the potential and forces at each time step can be evaluated by using a recurrence relation which avoids enumerating all permutations, while providing the correct thermal expectation values. The resulting algorithm scales cubically with system size. The method is tested and applied to bosons in a 2-dimensional (2D) trap and agrees with analytical results and numerical diagonalization of the many-body Hamiltonian. An analysis of the role of exchange effects at different temperatures, through the relative probability of different ring-polymer configurations, is also presented.

Neural networks-based variationally enhanced sampling.

Sampling complex free-energy surfaces is one of the main challenges of modern atomistic simulation methods. The presence of kinetic bottlenecks in such surfaces often renders a direct approach useless. A popular strategy is to identify a small number of key collective variables and to introduce a bias potential that is able to favor their fluctuations in order to accelerate sampling. Here, we propose to use machine-learning techniques in conjunction with the recent variationally enhanced sampling method [O. Valsson, M. Parrinello, Phys. Rev. Lett. 113, 090601 (2014)] in order to determine such potential. This is achieved by expressing the bias as a neural network. The parameters are determined in a variational learning scheme aimed at minimizing an appropriate functional. This required the development of a more efficient minimization technique. The expressivity of neural networks allows representing rapidly varying free-energy surfaces, removes boundary effects artifacts, and allows several collective variables to be handled.

Microscopic description of acid-base equilibrium.

Acid-base reactions are ubiquitous in nature. Understanding their mechanisms is crucial in many fields, from biochemistry to industrial catalysis. Unfortunately, experiments give only limited information without much insight into the molecular behavior. Atomistic simulations could complement experiments and shed precious light on microscopic mechanisms. The large free-energy barriers connected to proton dissociation, however, make the use of enhanced sampling methods mandatory. Here we perform an ab initio molecular dynamics (MD) simulation and enhance sampling with the help of metadynamics. This has been made possible by the introduction of descriptors or collective variables (CVs) that are based on a conceptually different outlook on acid-base equilibria. We test successfully our approach on three different aqueous solutions of acetic acid, ammonia, and bicarbonate. These are representative of acid, basic, and amphoteric behavior.

Water-Triggered, Irreversible Conformational Change of SARS-CoV-2 Main Protease on Passing from the Solid State to Aqueous Solution.

The main protease from SARS-CoV-2 is a homodimer. Yet, a recent 0.1-ms-long molecular dynamics simulation performed by D. E. Shaw's research group shows that it readily undergoes a symmetry-breaking event on passing from the solid state to aqueous solution. As a result, the subunits present distinct conformations of the binding pocket. By analyzing this long simulation, we uncover a previously unrecognized role of water molecules in triggering the transition. Interestingly, each subunit presents a different collection of long-lived water molecules. Enhanced sampling simulations performed here, along with machine learning approaches, further establish that the transition to the asymmetric state is essentially irreversible.

Liquid-Liquid Critical Point in Phosphorus.

The study of liquid-liquid phase transitions has attracted considerable attention. One interesting example of this phenomenon is phosphorus, for which the existence of a first-order phase transition between a low density insulating molecular phase and a conducting polymeric phase has been experimentally established. In this Letter, we model this transition by an ab initio quality molecular dynamics simulation and explore a large portion of the liquid section of the phase diagram. We draw the liquid-liquid coexistence curve and determine that it terminates into a second-order critical point. Close to the critical point, large coupled structure and electronic structure fluctuations are observed.

Predicting polymorphism in molecular crystals using orientational entropy.

We introduce a computational method to discover polymorphs in molecular crystals at finite temperature. The method is based on reproducing the crystallization process starting from the liquid and letting the system discover the relevant polymorphs. This idea, however, conflicts with the fact that crystallization has a timescale much longer than that of molecular simulations. To bring the process within affordable simulation time, we enhance the fluctuations of a collective variable by constructing a bias potential with well-tempered metadynamics. We use as a collective variable an entropy surrogate based on an extended pair correlation function that includes the correlation between the orientations of pairs of molecules. We also propose a similarity metric between configurations based on the extended pair correlation function and a generalized Kullback-Leibler divergence. In this way, we automatically classify the configurations as belonging to a given polymorph, using our metric and a hierarchical clustering algorithm. We apply our method to urea and naphthalene. We find different polymorphs for both substances, and one of them is stabilized at finite temperature by entropic effects.

Molecular dynamics simulations of liquid silica crystallization.

Silica is one of the most abundant minerals on Earth and is widely used in many fields. Investigating the crystallization of liquid silica by atomic simulations is of great importance to understand the crystallization mechanism; however, the high crystallization barrier and the tendency of silica to form glasses make such simulations very challenging. Here we have studied liquid silica crystallization to ß-cristobalite with metadynamics, using X-ray diffraction (XRD) peak intensities as collective variables. The frequent transitions between solid and liquid of the biased runs demonstrate the highly successful use of the XRD peak intensities as collective variables, which leads to the convergence of the free-energy surface. By calculating the difference in free energy, we have estimated the melting temperature of ß-cristobalite, which is in good agreement with the literature. The nucleation mechanism during the crystallization of liquid silica can be described by classical nucleation theory.

A metadynamics perspective on the reduction mechanism of the Pt(IV) asplatin prodrug.

Enhanced sampling molecular dynamics has been used to model the reduction mechanism of the antitumoral Asplatin Pt(IV) complex, c,c,t-[PtCl2(NH3)2(OH)(aspirin)] in the presence of l-ascorbic acid as reducing agent. In order to overcome the timescale problem, characteristic of many chemical reactions, we enhanced the sampling of the free energy landscape using Metadynamics. To achieve such a goal, the selection of adequate collective variables is crucial for the application of the method. Recently, a new method called Multi-Class Harmonic Linear Discriminant Analysis (MC-HLDA) has been proposed as a tool for constructing collective variables (CVs) for complex chemical processes. The method reduces the dimensionality of the variable space by generating appropriate linear combinations of several relevant chemical descriptors. The aim of this work is to assess the ability and performance of this method in describing the fundamental features of complex chemical reactions such as the Asplatin reduction mechanism in a compact, simple, and physically transparent manner. © 2019 Wiley Periodicals, Inc.

Metadynamics of Paths.

We present a method to sample reactive pathways via biased molecular dynamics simulations in trajectory space. We show that the use of enhanced sampling techniques enables unconstrained exploration of multiple reaction routes. Time correlation functions are conveniently computed via reweighted averages along a single trajectory and kinetic rates are accessed at no additional cost. These abilities are illustrated analyzing a model potential and the umbrella inversion of NH_{3} in water. The algorithm allows a parallel implementation and promises to be a powerful tool for the study of rare events.