Search

The recent surge in Generative Artificial Intelligence (AI) has introduced exciting possibilities for computational chemistry. Generative AI methods have made significant progress in sampling molecular structures across chemical species, developing force fields, and speeding up simulations. This Perspective offers a structured overview, beginning with the fundamental theoretical concepts in both Generative AI and computational chemistry. It then covers widely used Generative AI methods, including autoencoders, generative adversarial networks, reinforcement learning, flow models and language models, and highlights their selected applications in diverse areas including force field development, and protein/RNA structure prediction. A key focus is on the challenges these methods face before they become truly predictive, particularly in predicting emergent chemical phenomena. We believe that the ultimate goal of a simulation method or theory is to predict phenomena not seen before, and that Generative AI should be subject to these same standards before it is deemed useful for chemistry. We suggest that to overcome these challenges, future AI models need to integrate core chemical principles, especially from statistical mechanics.

The weighted ensemble (WE) method stands out as a widely used segment-based sampling technique renowned for its rigorous treatment of kinetics. The WE framework typically involves initially mapping the configuration space onto a low-dimensional collective variable (CV) space and then partitioning it into bins. The efficacy of WE simulations heavily depends on the selection of CVs and binning schemes. The recently proposed State Predictive Information Bottleneck (SPIB) method has emerged as a promising tool for automatically constructing CVs from data and guiding enhanced sampling through an iterative manner. In this work, we advance this data-driven pipeline by incorporating prior expert knowledge. Our hybrid approach combines SPIB-learned CVs to enhance sampling in explored regions with expert-based CVs to guide exploration in regions of interest, synergizing the strengths of both methods. Through benchmarking on alanine dipeptide and chignoin systems, we demonstrate that our hybrid approach effectively guides WE simulations to sample states of interest, and reduces run-to-run variances. Moreover, our integration of the SPIB model also enhances the analysis and interpretation of WE simulation data by effectively identifying metastable states and pathways, and offering direct visualization of dynamics.

In this work we examine the nucleation from NaCl aqueous solutions within nano-confined environments, employing enhanced sampling molecular dynamics simulations integrated with machine learning-derived reaction coordinates. Through our simulations, we successfully induce phase transitions between solid, liquid, and a hydrated phase, typically observed at lower temperatures in bulk environments. Interestingly, nano-confinement serves to stabilize the solid phase and elevate melting points. Our simulations explain these findings by underscoring the significant role of water, alongside ion aggregation and subtle, anistropic dielectric behavior, in driving nucleation within nano-constrained environments. This letter thus provides a framework for sampling, analyzing and understanding nucleation processes under nano-confinement.

Ligand unbinding is mediated by the free energy change, which has intertwined contributions from both energy and entropy. It is important but not easy to quantify their individual contributions. We model hydrophobic ligand unbinding for two systems, a methane particle and a C60 fullerene, both unbinding from hydrophobic pockets in all-atom water. By using a modified deep learning framework, we learn a thermodynamically optimized reaction coordinate to describe hydrophobic ligand dissociation for both systems. Interpretation of these reaction coordinates reveals the roles of entropic and enthalpic forces as ligand and pocket sizes change. Irrespective of the contrasting roles of energy and entropy, we also find that for both the systems the transition from the bound to unbound states is driven primarily by solvation of the pocket and ligand, independent of ligand size. Our framework thus gives useful thermodynamic insight into hydrophobic ligand dissociation problems that are otherwise difficult to glean., Comment: 27 pages; 5 figures

Even though nucleation is ubiquitous in different science and engineering problems, investigating nucleation is extremely difficult due to the complicated ranges of time and length scales involved. In this work, we simulate NaCl nucleation in both molten and aqueous environments using enhanced sampling all-atom molecular dynamics with deep learning-based estimation of reaction coordinates. By incorporating various structural order parameters and learning the reaction coordinate as a function thereof, we achieve significantly improved sampling relative to traditional ad hoc descriptions of what drives nucleation, particularly in the aqueous medium. Our results reveal a one-step nucleation mechanism in both environments, with reaction coordinate analysis highlighting the importance of local ion density in distinguishing solid and liquid states. However, while fluctuations in the local ion density are necessary to drive nucleation, they are not sufficient. Our analysis shows that near the transition states, descriptors such as enthalpy and local structure become crucial. Our protocol proposed here enables robust nucleation analysis and phase sampling, and could offer insights into nucleation mechanisms for generic small molecules in different environments.

Kinases compose one of the largest fractions of the human proteome, and their misfunction is implicated in many diseases, in particular cancers. The ubiquitousness and structural similarities of kinases makes specific and effective drug design difficult. In particular, conformational variability due to the evolutionarily conserved DFG motif adopting in and out conformations and the relative stabilities thereof are key in structure-based drug design for ATP competitive drugs. These relative conformational stabilities are extremely sensitive to small changes in sequence, and provide an important problem for sampling method development. Since the invention of AlphaFold2, the world of structure-based drug design has noticably changed. In spite of it being limited to crystal-like structure prediction, several methods have also leveraged its underlying architecture to improve dynamics and enhanced sampling of conformational ensembles, including AlphaFold2-RAVE. Here, we extend AlphaFold2-RAVE and apply it to a set of kinases: the wild type DDR1 sequence and three mutants with single point mutations that are known to behave drastically differently. We show that AlphaFold2-RAVE is able to efficiently recover the changes in relative stability using transferable learnt order parameters and potentials, thereby supplementing AlphaFold2 as a tool for exploration of Boltzmann-weighted protein conformations.

Phase transitions are ubiquitous across life, yet hard to quantify and describe accurately. In this work, we develop an approach for characterizing generic attributes of phase transitions from very limited observations made deep within different phases' domains of stability. Our approach is called Thermodynamic Maps, which combines statistical mechanics and molecular simulations with score-based generative models. Thermodynamic Maps enable learning the temperature dependence of arbitrary thermodynamic observables across a wide range of temperatures. We show its usefulness by calculating phase transition attributes such as melting temperature, temperature-dependent heat capacities, and critical exponents. For instance, we demonstrate the ability of thermodynamic maps to infer the ferromagnetic phase transition of the Ising model, including temperature-dependent heat capacity and critical exponents, despite never having seen samples from the transition region. In addition, we efficiently characterize the temperature-dependent conformational ensemble and compute melting curves of the two RNA systems GCAA tetraloop and HIV-TAR, which are notoriously hard to sample due to glassy-like landscapes.

Molecular dynamics (MD) enables the study of physical systems with excellent spatiotemporal resolution but suffers from severe time-scale limitations. To address this, enhanced sampling methods have been developed to improve exploration of configurational space. However, implementing these is challenging and requires domain expertise. In recent years, integration of machine learning (ML) techniques in different domains has shown promise, prompting their adoption in enhanced sampling as well. Although ML is often employed in various fields primarily due to its data-driven nature, its integration with enhanced sampling is more natural with many common underlying synergies. This review explores the merging of ML and enhanced MD by presenting different shared viewpoints. It offers a comprehensive overview of this rapidly evolving field, which can be difficult to stay updated on. We highlight successful strategies like dimensionality reduction, reinforcement learning, and flow-based methods. Finally, we discuss open problems at the exciting ML-enhanced MD interface., Comment: Submitted as invited article to Annual Review of Physical Chemistry vol 75; updated formatting issues

While representation learning has been central to the rise of machine learning and artificial intelligence, a key problem remains in making the learned representations meaningful. For this, the typical approach is to regularize the learned representation through prior probability distributions. However, such priors are usually unavailable or are ad hoc. To deal with this, recent efforts have shifted towards leveraging the insights from physical principles to guide the learning process. In this spirit, we propose a purely dynamics-constrained representation learning framework. Instead of relying on predefined probabilities, we restrict the latent representation to follow overdamped Langevin dynamics with a learnable transition density - a prior driven by statistical mechanics. We show this is a more natural constraint for representation learning in stochastic dynamical systems, with the crucial ability to uniquely identify the ground truth representation. We validate our framework for different systems including a real-world fluorescent DNA movie dataset. We show that our algorithm can uniquely identify orthogonal, isometric and meaningful latent representations.

When examining dynamics occurring at non-zero temperatures, both energy and entropy must be taken into account while describing activated barrier crossing events. Furthermore, good reaction coordinates need to be constructed to describe different metastable states and the transition mechanisms between them. Here we use a physics-based machine learning method called the State Predictive Information Bottleneck (SPIB) to find non-linear reaction coordinates for three systems of varying complexity. The SPIB is able to predict correctly an entropic bottleneck for an analytical flat-energy double-well system and identify the entropy- and energy-dominated pathways for an analytical four-well system. Finally, for a simulation of benzoic acid permeation through a lipid bilayer, SPIB is able to discover the entropic and energetic barriers to the permeation process. Given these results, we thus establish that SPIB is a reasonable and robust method for finding the important entropy and energy/enthalpy barriers in physical systems, which can then be used for enhanced understanding and sampling of different activated mechanisms.

We study NaCl ion-pair dissociation in a dilute aqueous solution using computer simulations both for the full system with long range Coulomb interactions and for a well chosen reference system with short range intermolecular interactions. Analyzing results using concepts from Local Molecular Field (LMF) theory and the recently proposed AI-based analysis tool "State predictive information bottleneck" (SPIB) we show that the system with short range interactions can accurately reproduce the transition rate for the dissociation process, the dynamics for moving between the underlying metastable states, and the transition state ensemble. Contributions from long range interactions can be largely neglected for these processes because long range forces from the direct interionic Coulomb interactions are almost completely canceled ($>90\%$) by those from solvent interactions over the length scale where the transition takes place. Thus for this important monovalent ion-pair system, short range forces alone are able to capture detailed consequences of the collective solvent motion, allowing the use of physically suggestive and computationally efficient short range models for the disassociation event. We believe that the framework here should be applicable to disentangling mechanisms for more complex processes such as multivalent ion disassociation, where previous work has suggested that long range contributions may be more important.

Understanding the driving forces behind the nucleation of different polymorphs is of great importance for material sciences and the pharmaceutical industry. This includes understanding the reaction coordinate that governs the nucleation process as well as correctly calculating the relative free energies of different polymorphs. Here we demonstrate, for the prototypical case of urea nucleation from melt, how one can learn such a 1-dimensional reaction coordinate as a function of pre-specified order parameters, and use it to perform efficient biased all-atom molecular dynamics simulations. The reaction coordinate is learnt as a function of generic thermodynamic and structural order parameters using the "Spectral Gap Optimization of Order Parameters (SGOOP)" approach [P. Tiwary and B. J. Berne, Proc. Natl. Acad. Sci. (2016)], and is biased using well-tempered metadynamics simulations. The reaction coordinate gives insight into the role played by different structural and thermodynamics order parameters, and the biased simulations obtain accurate relative free energies for different polymorphs. This includes accurate prediction of the approximate pressure at which urea undergoes a phase transition and one of the metastable polymorphs becomes the most stable conformation. We believe the ideas demonstrated in thus work will facilitate efficient sampling of nucleation in complex, generic systems.

The ability to make sense of the massive amounts of high-dimensional data generated from molecular dynamics (MD) simulations is heavily dependent on the knowledge of a low dimensional manifold (parameterized by a reaction coordinate or RC) that typically distinguishes between relevant metastable states and which captures the relevant slow dynamics of interest. Methods based on machine learning and artificial intelligence have been proposed over the years to deal with learning such low-dimensional manifolds, but they are often criticized for a disconnect from more traditional and physically interpretable approaches. To deal with such concerns, in this work, we propose a deep learning based State Predictive Information Bottleneck (SPIB) approach to learn the RC from high dimensional molecular simulation trajectories. We demonstrate analytically and numerically how the RC learnt in this approach is deeply connected to the committor in chemical physics, and can be used to accurately identify transition states. A crucial hyperparameter in this approach is the time-delay, or how far into the future the algorithm should make predictions about. Through careful comparisons for benchmark systems, we demonstrate that this hyperparameter choice gives useful control over how coarse-grained we want the metastable state classification of the system to be. We thus believe that this work represents a step forward in systematic application of deep learning based ideas to molecular simulations in a way that bridges the gap between artificial intelligence and traditional chemical physics., Comment: 11 pages, 13 figures

Recurrent neural networks (RNNs) have led to breakthroughs in natural language processing and speech recognition, wherein hundreds of millions of people use such tools on a daily basis through smartphones, email servers and other avenues. In this work, we show such RNNs, specifically Long Short-Term Memory (LSTM) neural networks can also be applied to capturing the temporal evolution of typical trajectories arising in chemical and biological physics. Specifically, we use a character-level language model based on LSTM. This learns a probabilistic model from 1-dimensional stochastic trajectories generated from molecular dynamics simulations of a higher dimensional system. We show that the model can not only capture the Boltzmann statistics of the system but it also reproduce kinetics at a large spectrum of timescales. We demonstrate how the embedding layer, introduced originally for representing the contextual meaning of words or characters, exhibits here a nontrivial connectivity between different metastable states in the underlying physical system. We demonstrate the reliability of our model and interpretations through different benchmark systems and a single molecule force spectroscopy trajectory for multi-state riboswitch. We anticipate that our work represents a stepping stone in the understanding and use of RNNs for modeling and predicting dynamics of complex stochastic molecular systems.

In this work, we revisit our recent iterative machine learning (ML) -- molecular dynamics (MD) technique "Reweighted autoencoded variational Bayes for enhanced sampling (RAVE)" (Ribeiro, Bravo, Wang, Tiwary, J. Chem. Phys. 149 072301 (2018) and Wang, Ribeiro, Tiwary, Nature Commun. 10 3573 (2019)) and analyze as well as formalize some of its approximations. These including: (a) the choice of a predictive time-delay, or how far into the future should the ML try to predict the state of a given system output from MD, and (b) for short time-delays, how much of an error is made in approximating the biased propagator for the dynamics as the unbiased propagator. We demonstrate through a master equation framework as to why the exact choice of time-delay is irrelevant as long as a small non-zero value is adopted. We also derive a correction to reweight the biased propagator, and somewhat to our dissatisfaction but also to our reassurance, find that it barely makes a difference to the intuitive picture we had previously derived and used.

Molecular dynamics (MD) has become a powerful tool for studying biophysical systems, due to increasing computational power and availability of software. Although MD has made many contributions to better understanding these complex biophysical systems, there remain methodological difficulties to be surmounted. First, how to make the deluge of data generated in running even a microsecond long MD simulation human comprehensible. Second, how to efficiently sample the underlying free energy surface and kinetics. In this short perspective, we summarize machine learning based ideas that are solving both of these limitations, with a focus on their key theoretical underpinnings and remaining challenges.

The ability to predict accurate thermodynamic and kinetic properties in biomolecular systems is of both scientific and practical utility. While both remain very difficult, predictions of kinetics are particularly difficult because rates, in contrast to free energies, depend on the route taken and are thus not amenable to all enhanced sampling methods. It has recently been demonstrated that it is possible to recover kinetics through so called `infrequent metadynamics' simulations, where the simulations are biased in a way that minimally corrupts the dynamics of moving between metastable states. This method, however, requires the bias to be added slowly, thus hampering applications to processes with only modest separations of timescales. Here we present a frequency-adaptive strategy which bridges normal and infrequent metadynamics. We show that this strategy can improve the precision and accuracy of rate calculations at fixed computational cost, and should be able to extend rate calculations for much slower kinetic processes., Comment: 15 pages, 2 figures, 2 tables

Here we propose the Reweighted Autoencoded Variational Bayes for Enhanced Sampling (RAVE) method, a new iterative scheme that uses the deep learning framework of variational autoencoders to enhance sampling in molecular simulations. RAVE involves iterations between molecular simulations and deep learning in order to produce an increasingly accurate probability distribution along a low-dimensional latent space that captures the key features of the molecular simulation trajectory. Using the Kullback-Leibler divergence between this latent space distribution and the distribution of various trial reaction coordinates sampled from the molecular simulation, RAVE determines an optimum, yet nonetheless physically interpretable, reaction coordinate and optimum probability distribution. Both then directly serve as the biasing protocol for a new biased simulation, which is once again fed into the deep learning module with appropriate weights accounting for the bias, the procedure continuing until estimates of desirable thermodynamic observables are converged. Unlike recent methods using deep learning for enhanced sampling purposes, RAVE stands out in that (a) it naturally produces a physically interpretable reaction coordinate, (b) is independent of existing enhanced sampling protocols to enhance the fluctuations along the latent space identified via deep learning, and (c) it provides the ability to easily filter out spurious solutions learned by the deep learning procedure. The usefulness and reliability of RAVE is demonstrated by applying it to model potentials of increasing complexity, including computation of the binding free energy profile for a hydrophobic ligand-substrate system in explicit water with dissociation time of more than three minutes, in computer time at least twenty times less than that needed for umbrella sampling or metadynamics.

We consider a range of model potentials with metastable states undergoing molecular dynamics coupled to a thermal bath in the high friction regime, and consider how the optimal reaction coordinate depends on the diffusion anisotropy. For this we use our recently proposed method 'Spectral gap optimization of order parameters (SGOOP)' (Tiwary and Berne, Proc. Natl. Acad. Sci. 113 2839 2016). We show how available information about dynamical observables in addition to static information can be incorporated into SGOOP, which can then be used to accurately determine the 'best' reaction coordinate for arbitrary anisotropies. We compare our results with transmission coefficient calculations and published benchmarks where applicable or available respectively.

We use a recently proposed method called Spectral Gap Optimization of Order Parameters (SGOOP) (Tiwary and Berne, Proc. Natl. Acad. Sci 2016, 113, 2839 (2016)), to determine an optimal 1-dimensional reaction coordinate (RC) for the unbinding of a bucky-ball from a pocket in explicit water. This RC is estimated as a linear combination of the multiple available order parameters that collectively can be used to distinguish the various stable states relevant for unbinding. We pay special attention to determining and quantifying the degree to which water molecules should be included in the RC. Using SGOOP with under-sampled biased simulations, we predict that water plays a distinct role in the reaction coordinate for unbinding in the case when the ligand is sterically constrained to move along an axis of symmetry. This prediction is validated through extensive calculations of the unbinding times through metadynamics, and by comparison through detailed balance with unbiased molecular dynamics estimate of the binding time. However when the steric constraint is removed, we find that the role of water in the reaction coordinate diminishes. Here instead SGOOP identifies a good one-dimensional RC involving various motional degrees of freedom., Comment: 7 pages, 5 figures

We consider the rate of transition for a particle between two metastable states coupled to a thermal environment for various magnitudes of the coupling strength, using the recently proposed infrequent metadynamics approach (Tiwary and Parrinello, Phys. Rev. Lett. 111, 230602 (2013)). We are interested in understanding how this approach for obtaining rate constants performs as the dynamics regime changes from energy diffusion to spatial diffusion. Reassuringly, we find that the approach works remarkably well for various coupling strengths in the strong coupling regime, and to some extent even in the weak coupling regime., Comment: 3 pages, 1 figure, submitted to J. Chem. Phys

Metadynamics is a popular enhanced sampling scheme wherein by periodic application of a repulsive bias, one can surmount high free energy barriers and explore complex landscapes. Recently metadynamics was shown to be mathematically well founded, in the sense that the biasing procedure is guaranteed to converge to the true free energy surface in the long time limit irrespective of the precise choice of biasing parameters. A differential equation governing the post-transient convergence behavior of metadynamics was also derived. In this short communication, we revisit this differential equation, expressing it in a convenient and elegant Riccati-like form. A perturbative solution scheme is then developed for solving this differential equation, which is valid for any generic biasing kernel. The solution clearly demonstrates the robustness of metadynamics to choice of biasing parameters and gives further confidence in the widely used method., Comment: submitted to J. Chem. Phys

In modern day simulations of many-body systems much of the computational complexity is shifted to the identification of slowly changing molecular order parameters called collective variables (CV) or reaction coordinates. A vast array of enhanced sampling methods are based on the identification and biasing of these low-dimensional order parameters, whose fluctuations are important in driving rare events of interest. Here describe a new algorithm for finding optimal low-dimensional collective variables for use in enhanced sampling biasing methods like umbrella sampling, metadynamics and related methods, when limited prior static and dynamic information is known about the system, and a much larger set of candidate CVs is specified. The algorithm involves estimating the best combination of these candidate CVs, as quantified by a maximum path entropy estimate of the spectral gap for dynamics viewed as a function of that CV. Through multiple practical examples, we show how this post-processing procedure can lead to optimization of CV and several orders of magnitude improvement in the convergence of the free energy calculated through metadynamics, essentially giving the ability to extract useful information even from unsuccessful metadynamics runs., Comment: 7 pages, 4 figures; corrected missing figure number and added a reference

A key factor influencing a drug's efficacy is its residence time in the binding pocket of the host protein. Using atomistic computer simulation to predict this residence time and the associated dissociation process is a desirable but extremely difficult task due to the long timescales involved. This gets further complicated by the presence of biophysical factors such as steric and solvation effects. In this work, we perform molecular dynamics (MD) simulations of the unbinding of a popular prototypical hydrophobic cavity-ligand system using a metadynamics based approach that allows direct assessment of kinetic pathways and parameters. When constrained to move in an axial manner, we find the unbinding time to be on the order of 4000 sec. In accordance with previous studies, we find that the ligand must pass through a region of sharp dewetting transition manifested by sudden and high fluctuations in solvent density in the cavity. When we remove the steric constraints on ligand, the unbinding happens predominantly by an alternate pathway, where the unbinding becomes 20 times faster, and the sharp dewetting transition instead becomes continuous. We validate the unbinding timescales from metadynamics through a Poisson analysis, and by comparison through detailed balance to binding timescale estimates from unbiased MD. This work demonstrates that enhanced sampling can be used to perform explicit solvent molecular dynamics studies at timescales previously unattainable, obtaining direct and reliable pictures of the underlying physio-chemical factors including free energies and rate constants., Comment: 7 pages, 4 figures, supplementary PDF file, submitted

Metadynamics is a commonly used and successful enhanced sampling method. By the introduction of a history dependent bias which depends on a restricted number of collective variables(CVs) it can explore complex free energy surfaces characterized by several metastable states separated by large free energy barriers. Here we extend its scope by introducing a simple yet powerful method for calculating the rates of transition between different metastable states. The method does not rely on a previous knowledge of the transition states or reaction co-ordinates, as long as CVs are known that can distinguish between the various stable minima in free energy space. We demonstrate that our method recovers the correct escape rates out of these stable states and also preserves the correct sequence of state-to-state transitions, with minimal extra computational effort needed over ordinary metadynamics. We apply the formalism to three different problems and in each case find excellent agreement with the results of long unbiased molecular dynamics runs., Comment: 4 pages, 2 figures, 1 supplemental file