Search

The problem of studying rare events is central to many areas of computer simulations. In a recent paper [Kang, P., et al., Nat. Comput. Sci. 4, 451-460, 2024], we have shown that a powerful way of solving this problem passes through the computation of the committor function, and we have demonstrated how the committor can be iteratively computed in a variational way and the transition state ensemble efficiently sampled. Here, we greatly ameliorate this procedure by combining it with a metadynamics-like enhanced sampling approach in which a logarithmic function of the committor is used as a collective variable. This integrated procedure leads to an accurate and balanced sampling of the free energy surface in which transition states and metastable basins are studied with the same thoroughness. We also show that our approach can be used in cases in which competing reactive paths are possible and intermediate metastable are encountered. In addition, we demonstrate how physical insights can be obtained from the optimized committor model and the sampled data, thus providing a full characterization of the rare event under study. We ascribe the success of this approach to the use of a probability-based description of rare events.

Graphs and hypergraphs provide powerful abstractions for modeling interactions among a set of entities of interest and have been attracting a growing interest in the literature thanks to many successful applications in several fields. In particular, they are rapidly expanding in domains such as chemistry and biology, especially in the areas of drug discovery and molecule generation. One of the areas witnessing the fasted growth is the chemical reactions field, where chemical reactions can be naturally encoded as directed hyperedges of a hypergraph. In this paper, we address the chemical reaction classification problem by introducing the notation of a Directed Line Graph (DGL) associated with a given directed hypergraph. On top of it, we build the Directed Line Graph Network (DLGNet), the first spectral-based Graph Neural Network (GNN) expressly designed to operate on a hypergraph via its DLG transformation. The foundation of DLGNet is a novel Hermitian matrix, the Directed Line Graph Laplacian, which compactly encodes the directionality of the interactions taking place within the directed hyperedges of the hypergraph thanks to the DLG representation. The Directed Line Graph Laplacian enjoys many desirable properties, including admitting an eigenvalue decomposition and being positive semidefinite, which make it well-suited for its adoption within a spectral-based GNN. Through extensive experiments on chemical reaction datasets, we show that DGLNet significantly outperforms the existing approaches, achieving on a collection of real-world datasets an average relative-percentage-difference improvement of 33.01%, with a maximum improvement of 37.71%.

Enhanced sampling simulations make the computational study of rare events feasible. A large family of such methods crucially depends on the definition of some collective variables (CVs) that could provide a low-dimensional representation of the relevant physics of the process. Recently, many methods have been proposed to semi-automatize the CV design by using machine learning tools to learn the variables directly from the simulation data. However, most methods are based on feed-forward neural networks and require as input some user-defined physical descriptors. Here, we propose to bypass this step using a graph neural network to directly use the atomic coordinates as input for the CV model. This way, we achieve a fully automatic approach to CV determination that provides variables invariant under the relevant symmetries, especially the permutational one. Furthermore, we provide different analysis tools to favor the physical interpretation of the final CV. We prove the robustness of our approach using different methods from the literature for the optimization of the CV, and we prove its efficacy on several systems, including a small peptide, an ion dissociation in explicit solvent, and a simple chemical reaction.

This letter illustrates the opinion of the molecular dynamics (MD) community on the need to adopt a new FAIR paradigm for the use of molecular simulations. It highlights the necessity of a collaborative effort to create, establish, and sustain a database that allows findability, accessibility, interoperability, and reusability of molecular dynamics simulation data. Such a development would democratize the field and significantly improve the impact of MD simulations on life science research. This will transform our working paradigm, pushing the field to a new frontier. We invite you to support our initiative at the MDDB community (https://mddbr.eu/community/)

We investigate learning the eigenfunctions of evolution operators for time-reversal invariant stochastic processes, a prime example being the Langevin equation used in molecular dynamics. Many physical or chemical processes described by this equation involve transitions between metastable states separated by high potential barriers that can hardly be crossed during a simulation. To overcome this bottleneck, data are collected via biased simulations that explore the state space more rapidly. We propose a framework for learning from biased simulations rooted in the infinitesimal generator of the process and the associated resolvent operator. We contrast our approach to more common ones based on the transfer operator, showing that it can provably learn the spectral properties of the unbiased system from biased data. In experiments, we highlight the advantages of our method over transfer operator approaches and recent developments based on generator learning, demonstrating its effectiveness in estimating eigenfunctions and eigenvalues. Importantly, we show that even with datasets containing only a few relevant transitions due to sub-optimal biasing, our approach recovers relevant information about the transition mechanism.

Determining the kinetic bottlenecks that make transitions between metastable states difficult is key to understanding important physical problems like crystallization, chemical reactions, or protein folding. In all these phenomena, the system spends a considerable amount of time in one metastable state before making a rare but important transition to a new state. The rarity of these events makes their direct simulation challenging, if not impossible. We propose a method to explore the distribution of configurations that the system passes as it translocates from one metastable basin to another. We shall refer to this set of configurations as the transition state ensemble. We base our method on the committor function and the variational principle to which it obeys. We find the minimum of the variational principle via a self-consistent procedure that does not require any input besides the knowledge of the initial and final state. Right from the start, our procedure focuses on sampling the transition state ensemble and allows harnessing a large number of such configurations. With the help of the variational principle, we perform a detailed analysis of the transition state ensemble, ranking quantitatively the degrees of freedom mostly involved in the transition and opening the way for a systematic approach for the interpretation of simulation results and the construction of collective variables.

The understanding of prototropic tautomerism in water and the characterization of solvent effects on protomeric equilibrium pose significant challenges. Using molecular dynamics simulations based on state-of-the-art deep learning potential and enhanced sampling methods, we provide a comprehensive description of all configurational transformations in glycine solvated in water and determine accurate free energy profiles of these processes. We observe that the tautomerism between the neutral and zwitterionic forms of solvated glycine can occur by both intramolecular proton transfer in glycine and intermolecular proton transfer in the contact ion pair (anionic glycine and hydronium ion) or the separated ion pair (cationic glycine and hydroxide ion).

A variety of enhanced sampling 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 enhanced sampling simulations in trajectory space via the metadynamics of paths [arXiv:2002.09281] 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 enhanced sampling simulations. We demonstrate our approach with two numerical examples, a two-dimensional model potential and the isomerization of alanine dipeptide, using deep targeted discriminant analysis as our ML-based CV of choice.

The anomalous $\lambda$-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 $\lambda$-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.

The wind mission Aeolus of the European Space Agency was a groundbreaking achievement for Earth observation. Between 2018 and 2023, the space-borne lidar instrument ALADIN onboard the Aeolus satellite measured atmospheric wind profiles with global coverage which contributed to improving the accuracy of numerical weather prediction. The precision of the wind observations, however, declined over the course of the mission due to a progressive loss of the atmospheric backscatter signal. The analysis of the root cause was supported by the Pierre Auger Observatory in Argentina whose fluorescence detector registered the ultraviolet laser pulses emitted from the instrument in space, thereby offering an estimation of the laser energy at the exit of the instrument for several days in 2019, 2020 and 2021. The reconstruction of the laser beam not only allowed for an independent assessment of the Aeolus performance, but also helped to improve the accuracy in the determination of the laser beam's ground track on single pulse level. The results presented in this paper set a precedent for the monitoring of space lasers by ground-based telescopes and open new possibilities for the calibration of cosmic-ray observatories., Comment: 10 pages, 10 figures

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 (DTW) 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.

This work examines the conformational ensemble involved in $\beta$-hairpin folding by means of advanced molecular dynamics simulations and dimensionality reduction. A fully atomistic description of the protein and the surrounding solvent molecules is used and this complex energy landscape is sampled by means of parallel tempering metadynamics simulations. The ensemble of configurations explored is analysed using the recently proposed sketch-map algorithm. Further simulations allow us to probe how mutations affect the structures adopted by this protein. We find that many of the configurations adopted by a mutant are the same as those adopted by the wild type protein. Furthermore, certain mutations destabilize secondary structure containing configurations by preventing the formation of hydrogen bonds or by promoting the formation of new intramolecular contacts. Our analysis demonstrates that machine-learning techniques can be used to study the energy landscapes of complex molecules and that the visualizations that are generated in this way provide a natural basis for examining how the stabilities of particular configurations of the molecule are affected by factors such as temperature or structural mutations., Comment: Pre-print of doi:10.1021/ct500950z, reproduced with permission from the ACS

Interatomic potentials learned using machine learning methods have been successfully applied to atomistic simulations. However, accurate models require large training datasets, while generating reference calculations is computationally demanding. To bypass this difficulty, we propose a transfer learning algorithm that leverages the ability of graph neural networks (GNNs) to represent chemical environments together with kernel mean embeddings. We extract a feature map from GNNs pre-trained on the OC20 dataset and use it to learn the potential energy surface from system-specific datasets of catalytic processes. Our method is further enhanced by incorporating into the kernel the chemical species information, resulting in improved performance and interpretability. We test our approach on a series of realistic datasets of increasing complexity, showing excellent generalization and transferability performance, and improving on methods that rely on GNNs or ridge regression alone, as well as similar fine-tuning approaches., Comment: 20 pages, 4 figures, 7 tables, published in NeurIPS 2023

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 $\texttt{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.

The phase change compound Ge$_2$Sb$_2$Te$_5$ (GST225) is exploited in advanced non-volatile electronic memories and in neuromorphic devices which both rely on a fast and reversible transition between the crystalline and amorphous phases induced by Joule heating. The crystallization kinetics of GST225 is a key functional feature for the operation of these devices. We report here on the development of a machine-learned interatomic potential for GST225 that allowed us to perform large scale molecular dynamics simulations (over 10000 atoms for over 100 ns) to uncover the details of the crystallization kinetics in a wide range of temperatures of interest for the programming of the devices. The potential is obtained by fitting with a deep neural network (NN) scheme a large quantum-mechanical database generated within Density Functional Theory. The availability of a highly efficient and yet highly accurate NN potential opens the possibility to simulate phase change materials at the length and time scales of the real devices.

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 in short unbiased simulation in the two basins. Here we enrich the set of data on which the Deep Targeted Discriminant Analysis collective variable is built by adding data coming 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 a more accurate sampling and faster convergence. The performance of these new collective variables is tested on a number of representative examples.

This work studies a well-known shared-cache coded caching scenario where each cache can serve an arbitrary number of users, analyzing the case where there is some knowledge about such number of users (i.e., the topology) during the content placement phase. Under the assumption of regular placement and a cumulative cache size that can be optimized across the different caches, we derive the fundamental limits of performance by introducing a novel cache-size optimization and placement scheme and a novel information-theoretic converse. The converse employs new index coding techniques to bypass traditional uniformity requirements, thus finely capturing the heterogeneity of the problem, and it provides a new approach to handle asymmetric settings. The new fundamental limits reveal that heterogeneous topologies can in fact outperform their homogeneous counterparts where each cache is associated to an equal number of users. These results are extended to capture the scenario of topological uncertainty where the perceived/estimated topology does not match the true network topology. This scenario is further elevated to the stochastic setting where the user-to-cache association is random and unknown, and it is shown that the proposed scheme is robust to such noisy or inexact knowledge on the topology., Comment: 39 pages

We present a new method to compute free energies at a quantum mechanical (QM) level of theory from molecular simulations using cheap reference potential energy functions, such as force fields. To overcome the poor overlap between the reference and target distributions, we generalize targeted free energy perturbation (TFEP) to employ multiple configuration maps. While TFEP maps have been obtained before from an expensive training of a normalizing flow neural network (NN), our multimap estimator allows us to use the same set of QM calculations to both optimize the maps and estimate the free energy, thus removing almost completely the overhead due to training. A multimap extension of the multistate Bennett acceptance ratio estimator is also derived for cases where samples from two or more states are available. Furthermore, we propose a one-epoch learning policy that can be used to efficiently avoid overfitting when computing the loss function is expensive compared to generating data. Finally, we show how our multimap approach can be combined with enhanced sampling strategies to overcome the pervasive problem of poor convergence due to slow degrees of freedom. We test our method on the HiPen dataset of drug-like molecules and fragments, and we show that it can accelerate the calculation of the free energy difference of switching from a force field to a DFTB3 potential by about 3 orders of magnitude compared to standard FEP and by a factor of about 8 compared to previously published nonequilibrium calculations., Comment: Main Text: 12 pages, 5 figures, 7 equations. Supplemental Material: 17 pages, 5 figures, 22 equations

Mathematical oncology provides unique and invaluable insights into tumour growth on both the microscopic and macroscopic levels. This review presents state-of-the-art modelling techniques and focuses on their role in understanding glioblastoma, a malignant form of brain cancer. For each approach, we summarise the scope, drawbacks, and assets. We highlight the potential clinical applications of each modelling technique and discuss the connections between the mathematical models and the molecular and imaging data used to inform them. By doing so, we aim to prime cancer researchers with current and emerging computational tools for understanding tumour progression. Finally, by providing an in-depth picture of the different modelling techniques, we also aim to assist researchers who seek to build and develop their own models and the associated inference frameworks., Comment: 30 pages, 3 figures, 3 tables

A key step in the development of a new drug is the design of drug-excipient complexes that lead to an optimal drug release kinetics. Computational chemistry, 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 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 like the montmorillonite clay would be a very convenient excipient. We simulate the drug release both from surface and interlayer space and find that the time it takes to the drug to diffuse from the interlayer space to the solution controls the release time.

We introduce a novel enhanced sampling approach named OPES flooding for calculating the kinetics of rare events from atomistic molecular dynamics simulation. This method is derived from the On-the-fly-Probability-Enhanced-Sampling (OPES) approach [Invernizzi and Parrinello, JPC Lett. 2020], 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 gas phase, and the folding and unfolding of the chignolin polypeptide in 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.

Present-day atomistic simulations generate long trajectories of ever more complex systems. Analyzing these data, discovering metastable states, and uncovering their nature is 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., Comment: Main text: 10 pages, 4 figures. Supplementary Info: 4 pages, 5, figures

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.

We combine molecular dynamics simulations with experiments to estimate solubilities of organic salts 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.

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 on prototypical systems and show that it outperforms the popular metadynamics method.

Multi-antenna coded caching combines a global caching gain, proportional to the total cache size in the network, with an additional spatial multiplexing gain that stems from multiple transmitting antennas. However, classic centralized coded caching schemes are not suitable for dynamic networks as they require prior knowledge of the number of users to indicate what data should be cached at each user during the placement phase. On the other hand, fully decentralized schemes provide comparable gains to their centralized counterparts only when the number of users is very large. In this paper, we propose a novel multi-antenna coded caching scheme for dynamic networks, where instead of defining individual cache contents, we associate users with a limited set of predefined caching profiles. Then, during the delivery phase, we aim at achieving a combined caching and spatial multiplexing gain, comparable to a large extent with the ideal case of fully centralized schemes. The resulting scheme imposes small subpacketization and beamforming overheads, is robust under dynamic network conditions, and incurs small finite-SNR performance loss compared with centralized schemes.

The long timescale associated with ligand residence times renders their computation challenging. Therefore, the influence of factors like solvation and steric hindrance on residence times are 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 understanding their unbinding mechanism. We observe that guest unbinding often proceeds via a series of intermediate states that can be labelled by the number of water molecules present in the binding cavity. And in several cases the residence time is correlated to the water trapping times in the cavity., Comment: 11 pages, including Supplementary Information, 5 figures