Exploration vs Convergence Speed in Adaptive-Bias Enhanced Sampling

Lambda-ABF-OPES: Faster Convergence with High Accuracy in Alchemical Free Energy Calculations

Predicting the binding affinity between small molecules and target macromolecules while combining both speed and accuracy is a cornerstone of modern computational drug discovery, which is critical for accelerating therapeutic development. Despite recent progress in molecular dynamics (MD) simulations, such as advanced polarizable force fields and enhanced sampling techniques, estimating absolute binding free energies (ABFEs) remains computationally challenging. To overcome these difficulties, we introduce a highly efficient hybrid methodology that couples the Lambda-adaptive biasing force (Lambda-ABF) scheme with on-the-fly probability enhanced sampling (OPES). This approach achieves up to a 9-fold improvement in sampling efficiency and computational speed compared to the original Lambda-ABF when used in conjunction with the AMOEBA polarizable force field, yielding converged results at a fraction of the cost of standard techniques.

Everything everywhere all at once: a probability-based enhanced sampling approach to rare events

The problem of studying rare events is central to many areas of computer simulations. We recently proposed an approach to solving this problem that involves computing the committor function, showing how it can be iteratively computed in a variational way while efficiently sampling the transition state ensemble. Here we greatly improve 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 procedure leads to an accurate sampling of the free energy surface in which transition states and metastable basins are studied with the same thoroughness. We show that our approach can be used in cases with the possibility of competing reactive paths and metastable intermediates. 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.

Data efficient learning of molecular slow modes from nonequilibrium metadynamics

Enhanced sampling simulations help overcome free energy barriers and explore molecular conformational space by applying external bias potential along suitable collective variables (CVs). However, identifying optimal CVs that align with the slow modes of complex molecular systems with many coupled degrees of freedom can be a significant challenge. Deep time-lagged independent component analysis (Deep-TICA) addresses this issue by employing an artificial neural network that generates non-linear combinations of molecular descriptors to learn the slowest degrees of freedom. Training Deep-TICA CVs, however, typically requires long equilibrium simulations that can sample multiple recrossing events across various metastable conformations of the molecule. This requirement can often be prohibitively expensive, thereby limiting its widespread application. In this study, we present an algorithm that enables the training of Deep-TICA CVs using a limited amount of trajectory data obtained from short, non-equilibrium metadynamics simulations that only sample one forward transition from the initial to the final state. We achieve this by utilizing the variational Koopman algorithm, which reweights short off-equilibrium trajectories to reflect the equilibrium probability densities. We demonstrate that enhanced sampling simulations conducted along the Koopman reweighted Deep-TICA CV can accurately and efficiently converge the free energy surface for systems such as the Müller-Brown potential, alanine dipeptide, and the chignolin mini-protein. Our approach, therefore, addresses the key challenge of inferring slow modes from limited trajectory data, making it more feasible to use deep learning CVs to study molecular processes of practical relevance.

Slow dynamical modes from static averages

In recent times, efforts have been made to describe the evolution of a complex system not through long trajectories but via the study of probability distribution evolution. This more collective approach can be made practical using the transfer operator formalism and its associated dynamics generator. Here, we reformulate in a more transparent way the result of Devergne et al. [Adv. Neural Inform. Process. Syst. 37, 75495-75521 (2024)] and show that the lowest eigenfunctions and eigenvalues of the dynamics generator can be efficiently computed using data easily obtainable from biased simulations. We also show explicitly that the long time dynamics can be reconstructed by using the spectral decomposition of the dynamics operator.

PLUMED Tutorials: A collaborative, community-driven learning ecosystem

In computational physics, chemistry, and biology, the implementation of new techniques in shared and open-source software lowers barriers to entry and promotes rapid scientific progress. However, effectively training new software users presents several challenges. Common methods like direct knowledge transfer and in-person workshops are limited in reach and comprehensiveness. Furthermore, while the COVID-19 pandemic highlighted the benefits of online training, traditional online tutorials can quickly become outdated and may not cover all the software's functionalities. To address these issues, here we introduce "PLUMED Tutorials," a collaborative model for developing, sharing, and updating online tutorials. This initiative utilizes repository management and continuous integration to ensure compatibility with software updates. Moreover, the tutorials are interconnected to form a structured learning path and are enriched with automatic annotations to provide broader context. This paper illustrates the development, features, and advantages of PLUMED Tutorials, aiming to foster an open community for creating and sharing educational resources.

Long-Range Electrostatics in Serine Proteases: Machine Learning-Driven Reaction Sampling Yields Insights for Enzyme Design

Computational enzyme design is a promising technique for producing novel enzymes for industrial and clinical needs. A key challenge that this technique faces is to consistently achieve the desired activity. Fundamental studies of natural enzymes revealed critical contributions from second-shell - and even more distant - residues to their remarkable efficiency. In particular, such residues organize the internal electrostatic field to promote the reaction. Engineering such fields computationally proved to be a promising strategy, which, however, has some limitations. Charged residues necessarily form specific patterns of local interactions that may be exploited for structural integrity. As a result, it is impossible to probe the electrostatic field alone by substituting amino acids. We hypothesize that an approach that isolates the influences of residues' charges from other influences could yield deeper insights. We use molecular modeling with AI-enhanced QM/MM reaction sampling to implement such an approach and apply it to a model serine protease subtilisin. We find that the negative charge 8 Å away from the catalytic site is crucial to achieving the enzyme's catalytic efficiency, contributing more than 2 kcal/mol to lowering the barrier. In contrast, a positive charge from the second-closest charged residue opposes the efficiency of the reaction by raising the barrier by 0.8 kcal/mol. This result invites discussion into the role of this residue and trade-offs that might have taken place in the evolution of such enzymes. Our approach is transferable and can help investigate the evolution of electrostatic preorganization in other enzymes. We believe that the study and engineering of electrostatic fields in enzymes is a promising direction to advance both fundamental and applied enzymology and lead to the design of new powerful biocatalysts.

Propensity of Water Self-Ions at Air(Oil)-Water Interfaces Revealed by Deep Potential Molecular Dynamics with Enhanced Sampling

The preference of water self-ions (hydronium and hydroxide) toward air/oil-water interfaces is one of the hottest topics in water research due to its importance for understanding properties, phenomena, and reactions of interfaces. In this work, we performed enhanced-sampling molecular dynamics simulations based on state-of-the-art neural network potentials with approximate M06-2X accuracy to investigate the propensity of hydronium and hydroxide ions at air/oil(decane)-water interfaces, which can simultaneously describe well the water autoionization process forming these ions, the recombination of ions, and the ionic distribution along the normal distance to the interface by employing a set of appropriate Voronoi collective variables. A stable ionic double-layer distribution is observed near the air-water interface, while the distribution is different at oil-water interfaces, where hydronium tends to be repelled from the interface into the bulk water, whereas hydroxide, with an interfacial stabilization free energy of -0.6 kcal/mol, is enriched in the interfacial layer. Through simulations of oil droplets in water, we further reveal that the interfacial propensity of hydroxide ions is caused by the positive charge distribution of the oil-water interface contributed by hydrogens of the dangling OH bonds of the interfacial water layer and the outermost layer decane molecules lying flat on the droplet. The present results may aid in understanding the acid-base nature of water interfaces with wide applications.

Computational advances in discovering cryptic pockets for drug discovery

A number of promising therapeutic target proteins have been considered "undruggable" due to the lack of well-defined ligandable pockets. Substantial research in protein dynamics has elucidated the existence of "cryptic" pockets that only exist transiently and become favorable for binding in the presence of a ligand. These pockets provide an avenue to target challenging proteins, inspiring the development of multiple computational methods. This review highlights established cryptic pocket modeling approaches like mixed solvent molecular dynamics and recent applications of enhanced sampling and AI-based methods in therapeutically relevant proteins.

Accelerating Molecular Dynamics through Informed Resetting

We present a procedure for enhanced sampling of molecular dynamics simulations through informed stochastic resetting. Many phenomena, such as protein folding and crystal nucleation, occur over time scales inaccessible in standard simulations. We recently showed that stochastic resetting can accelerate molecular simulations that exhibit broad transition time distributions. However, standard stochastic resetting does not exploit any information about the reaction progress. For a model system and chignolin in explicit water, we demonstrate that an informed resetting protocol leads to greater accelerations than standard stochastic resetting in molecular dynamics and Metadynamics simulations. This is achieved by resetting only when a certain condition is met, e.g., when the distance from the target along the reaction coordinate is larger than some threshold. We use these accelerated simulations to infer important kinetic observables such as the unbiased mean first-passage time and direct transit time. For the latter, Metadynamics with informed resetting leads to speedups of 2-3 orders of magnitude over unbiased simulations with relative errors of only ∼35-70%. Our work significantly extends the applicability of stochastic resetting for enhanced sampling of molecular simulations.

Molecular dynamics-guided reaction discovery reveals endoperoxide-to-alkoxy radical isomerization as key branching point in α-pinene ozonolysis

Secondary organic aerosols (SOAs) significantly impact Earth's climate and human health. Although the oxidation of volatile organic compounds (VOCs) has been recognized as the major contributor to the atmospheric SOA budget, the mechanisms by which this process produces SOA-forming highly oxygenated organic molecules (HOMs) remain unclear. A major challenge is navigating the complex chemical landscape of these transformations, which traditional hypothesis-driven methods fail to thoroughly investigate. Here, we explore the oxidation of α-pinene, a critical atmospheric biogenic VOC, using a novel reaction discovery approach based on molecular dynamics and state-of-the-art enhanced sampling techniques. Our approach successfully identifies all established reaction pathways of α-pinene ozonolysis, as well as discovers multiple novel species and pathways without relying on a priori chemical knowledge. In particular, we unveil a key branching point that leads to the rapid formation of alkoxy radicals, whose high and diverse reactivity help to explain hitherto unexplained oxidation pathways suggested by mass spectral peaks observed in α-pinene ozonolysis experiments. This branching point is likely prevalent across a variety of atmospheric VOCs and could be crucial in establishing the missing link to SOA-forming HOMs.

Enhanced Sampling with Suboptimal Collective Variables: Reconciling Accuracy and Convergence Speed

We introduce an enhanced sampling algorithm to obtain converged free energy landscapes of molecular rare events, even when the collective variable (CV) used for biasing is not optimal. Our approach samples a time-dependent target distribution by combining the on-the-fly probability enhanced sampling and its exploratory variant, OPES Explore. This promotes more transitions between the relevant metastable states and accelerates the convergence speed of the free energy estimate. We demonstrate the successful application of this combined algorithm on the two-dimensional Wolfe-Quapp potential, millisecond time-scale ligand-receptor binding in the trypsin-benzamidine complex, and folding-unfolding transitions in chignolin mini-protein. Our proposed algorithm can compute accurate free energies at an affordable computational cost and is robust in terms of the choice of CVs, making it particularly promising for the simulation of complex biomolecular systems.

Free Energy-Based Refinement of DNA Force Field via Modification of Multiple Nonbonding Energy Terms

The amber-OL21 force field (ff) was developed to better describe noncanonical DNA, including Z-DNA. Despite its improvements for DNA simulations, this study found that OL21's scope of application was limited by embedded ff artifacts. In a benchmark set of seven DNA molecules, including two double-stranded DNAs transitioning between B- and Z-DNA and five single-stranded DNAs folding into mini-dumbbell or G-quadruplex structures, the free energy landscapes obtained using OL21 revealed several issues: Z-DNA was overly stabilized; misfolded states in mini-dumbbell DNAs were most stable; DNA GQ folding was consistently biased toward an antiparallel topology. To address these issues, a simple van der Waals (vdW) correction scheme, referred to as vdW5, was proposed for OL21. This involved revising multiple nonbonding energy terms to improve the overall quality of the free energy landscapes. The vdW5 correction effectively eliminated the artifacts in OL21, providing significantly improved free energy representations for DNAs tested. The vdW5-level revision substantially enhanced the predictive power of DNA simulations, thereby extending the scope of application for amber-OL21.

Probing the Conformational Ensemble of the Amyloid Beta 16-22 Fragment with Parallel-Bias Metadynamics

Aβ(16-22) is a segment of the Alzheimer's-related β-amyloid peptide that plays a crucial role in its aggregation. This study applies well-tempered parallel-bias metadynamics to investigate the impact of several denaturants and osmolytes on the conformational ensembles of both termini-capped and uncapped Aβ(16-22) monomers. Comparison of the different sets of collective variables in the metadynamics bias shows that using the set of backbone torsional angles results in better and faster convergence of simulations than employing more general structural characteristics of the short peptide. The equilibrium conformational ensembles of the peptides are characterized in pure water and in the presence of TMAO, urea, guanidinium chloride, and trifluoroethanol. In particular, trifluoroethanol and TMAO are found to increase the population of compact peptide conformations, whereas urea and guanidinium chloride favor extended structures. The analysis of the free energy surfaces in the presence of various substances with a comparison of the behavior of the capped and uncapped peptide forms reveals the role of different types of intrapeptide interactions such as salt bridges, hydrophobic contacts, and hydrogen bonds in stabilization of the compact or extended structures. As compounds reducing the abundance of the compact states of Aβ(16-22) and other disordered peptides are likely to suppress their amyloid fibril formation, simulations in the systems with this short peptide may be useful for the virtual screening of such compounds.

Binding Modes and Water-Mediation of Polyelectrolyte Adsorption to a Neutral CaCO3 Surface

Aqueous polyelectrolytes are effective mineralization inhibitors due to their ability to template onto crystal surfaces and chelate ions in solution. These additives have been shown to alter the morphology of calcium carbonate crystals, making them promising candidates for biological and industrial applications. However, while key to designing more effective mineralization inhibitors, the molecular mechanisms governing the interactions between polyelectrolytes and crystal surfaces remain poorly understood. In this study, we investigate the adsorption of poly(acrylic acid) (PAA) on the dominant calcite ( 101̅ 4 ) cleavage plane using all-atom molecular dynamics simulations. Although the calcite slab is electrostatically neutral, its charge distribution induces a strong electrostatic potential in an aqueous solution, which leads to significant water structuring at the interface. We observe a very favorable adsorption affinity of the polyelectrolyte chain to the surface, yet the structure of the interfacial water is not significantly affected. Direct interactions between the monomers on the polyelectrolyte and the calcite surface are infrequent, despite variations in chain length, charge density of the polyelectrolyte, and solution conditions. Intriguingly, the polyelectrolyte interaction with the calcite surface is dominantly mediated through bridging hydrogen bond interactions. As the polyelectrolyte adsorbs to the surface, the chain conformation adapts to the interfacial water structure by increasing polyelectrolyte-water contacts and integrates into pre-existing hydrogen bond networks. We found that water-mediated interactions are more dominant than direct interactions between the polyelectrolyte and the surface. This suggests an alternative pathway to the widely accepted notion that entropic effects due to water reorganization are the primary driving force. These results suggest that the polyelectrolyte binding affinity can be tuned by altering the polymer chain interactions with the interfacial water structure in addition to the surface itself.

Host-Guest Binding Free Energies à la Carte: An Automated OneOPES Protocol

Estimating absolute binding free energies from molecular simulations is a key step in computer-aided drug design pipelines, but the agreement between computational results and experiments is still very inconsistent. Both the accuracy of the computational model and the quality of the statistical sampling contribute to this discrepancy, yet disentangling the two remains a challenge. In this study, we present an automated protocol based on OneOPES, an enhanced sampling method that exploits replica exchange and can accelerate several collective variables to address the sampling problem. We apply this protocol to 37 host-guest systems. The simplicity of setting up the simulations and producing well-converged binding free energy estimates without the need to optimize simulation parameters provides a reliable solution to the sampling problem. This, in turn, allows for a systematic force field comparison and ranking according to the correlation between simulations and experiments, which can inform the selection of an appropriate model. The protocol can be readily adapted to test more force field combinations and study more complex protein-ligand systems, where the choice of an appropriate physical model is often based on heuristic considerations rather than systematic optimization.

Investigating Ligand-Mediated Conformational Dynamics of Pre-miR21: A Machine-Learning-Aided Enhanced Sampling Study

MicroRNAs (miRs) are short, noncoding RNA strands that regulate the activity of mRNAs by affecting the repression of protein translation, and their dysregulation has been implicated in several pathologies. miR21 in particular has been implicated in tumorigenesis and anticancer drug resistance, making it a critical target for drug design. miR21 biogenesis involves precise biochemical pathways, including the cleavage of its precursor, pre-miR21, by the enzyme Dicer. The present work investigates the conformational dynamics of pre-miR21, focusing on the role of adenine29 in switching between Dicer-binding-prone and inactive states. We also investigated the effect of L50, a cyclic peptide binder of pre-miR21 and a weak inhibitor of its processing. Using time series data and our novel collective variable-based enhanced sampling technique, OneOPES, we simulated these conformational changes and assessed the effect of L50 on the conformational plasticity of pre-miR21. Our results provide insight into peptide-induced conformational changes and pave the way for the development of a computational platform for the screening of inhibitors of pre-miR21 processing that considers RNA flexibility, a stepping stone for effective structure-based drug design, with potentially broad applications in drug discovery.

Overcoming Sampling Issues and Improving Computational Efficiency in Collective-Variable-Based Enhanced-Sampling Simulations: A Tutorial

This tutorial is designed to help users overcome sampling challenges and improve computational efficiency in collective-variable (CV)-based enhanced-sampling, or importance-sampling, simulations. Toward this end, we introduce well-tempered metadynamics-extended adaptive biasing force (WTM-eABF) and its integration with Gaussian accelerated molecular dynamics (GaMD). Additionally, use will be made of a method for identifying the least-free-energy pathway (LFEP) and multiple concurrent pathways on high-dimensional free-energy surfaces. We illustrate these sampling techniques with the conformational equilibria of trialanine and chignolin in aqueous solution as test cases. This tutorial assumes that the user has prior experience with molecular dynamics (MD) simulations, in general, with the popular program NAMD, and to some extent with Colvars, the module for CV-based calculations. This tutorial can, however, in large measure be used in conjunction with alternate MD engines that support the Colvars module such as GROMACS, LAMMPS, and Tinker-HP.

Enhanced Sampling of Biomolecular Slow Conformational Transitions Using Adaptive Sampling and Machine Learning

Biomolecular simulations often suffer from the "time scale problem", hindering the study of rare events occurring over extended time scales. Enhanced sampling techniques aim to alleviate this issue by accelerating conformational transitions, yet they typically necessitate well-defined collective variables (CVs), posing a significant challenge. Machine learning offers promising solutions but typically requires rich training data encompassing the entire free energy surface (FES). In this work, we introduce an automated iterative pipeline designed to mitigate these limitations. Our protocol first utilizes a CV-free count-based adaptive sampling method to generate a data set rich in rare events. From this data set, slow modes are identified using Koopman-reweighted time-lagged independent component analysis (KTICA), which are subsequently leveraged by on-the-fly probability enhanced sampling (OPES) to efficiently explore the FES. The effectiveness of our pipeline is demonstrated and further compared with the common Markov State Model (MSM) approach on two model systems with increasing complexity: alanine dipeptide (Ala2) and deca-alanine (Ala10), underscoring its applicability across diverse biomolecular simulations.

A Mode Evolution Metric to Extract Reaction Coordinates for Biomolecular Conformational Transitions

The complex, multidimensional energy landscape of biomolecules makes the extraction of suitable, nonintuitive collective variables (CVs) that describe their conformational transitions challenging. At present, dimensionality reduction approaches and machine learning (ML) schemes are employed to obtain CVs from molecular dynamics (MD)/Monte Carlo (MC) trajectories or structural databanks for biomolecules. However, minimum sampling conditions to generate reliable CVs that accurately describe the underlying energy landscape remain unclear. Here, we address this issue by developing a Mode evolution Metric (MeM) to extract CVs that can pinpoint new states and describe local transitions in the vicinity of a reference minimum from nonequilibrated MD/MC trajectories. We present a general mathematical formulation of MeM for both statistical dimensionality reduction and machine learning approaches. Application of MeM to MC trajectories of model potential energy landscapes and MD trajectories of solvated alanine dipeptide reveals that the principal components which locate new states in the vicinity of a reference minimum emerge well before the trajectories locally equilibrate between the associated states. Finally, we demonstrate a possible application of MeM in designing efficient biased sampling schemes to construct accurate energy landscape slices that link transitions between states. MeM can help speed up the search for new minima around a biomolecular conformational state and enable the accurate estimation of thermodynamics for states lying on the energy landscape and the description of associated transitions.

Absolute Binding Free Energies with OneOPES

The calculation of absolute binding free energies (ABFEs) for protein-ligand systems has long been a challenge. Recently, refined force fields and algorithms have improved the quality of the ABFE calculations. However, achieving the level of accuracy required to inform drug discovery efforts remains difficult. Here, we present a transferable enhanced sampling strategy to accurately calculate absolute binding free energies using OneOPES with simple geometric collective variables. We tested the strategy on two protein targets, BRD4 and Hsp90, complexed with a total of 17 chemically diverse ligands, including both molecular fragments and drug-like molecules. Our results show that OneOPES accurately predicts protein-ligand binding affinities with a mean unsigned error within 1 kcal mol-1 of experimentally determined free energies, without the need to tailor the collective variables to each system. Furthermore, our strategy effectively samples different ligand binding modes and consistently matches the experimentally determined structures regardless of the initial protein-ligand configuration. Our results suggest that the proposed OneOPES strategy can be used to inform lead optimization campaigns in drug discovery and to study protein-ligand binding and unbinding mechanisms.

A molecular dynamics study on the ion-mediated self-assembly of monolayer-protected nanoclusters

We studied the effects of metal and molecular cations on the aggregation of atomically precise monolayer-protected nanoclusters (MPCs) in an explicit solvent using atomistic molecular dynamics simulations. While divalent cations such as Zn2+ and Cd2+ promote aggregation by forming ligand-cation-ligand bridges between the MPCs, molecular cations such as tetraethylammonium and cholinium inhibit their aggregation by getting adsorbed into the MPC's ligand shell and reducing the ligand's motion. Here, we studied the aggregation of Au25(SR)18 nanoclusters with two types of ligands, para-mercaptobenzoic acid and D-penicillamine, as prototypical examples.

Analogies and Differences in the Photoactivation Mechanism of Bathy and Canonical Bacteriophytochromes Revealed by Multiscale Modeling

Bacteriophytochromes are light-sensing biological machines that switch between two photoreversible states, Pr and Pfr. Their relative stability is opposite in canonical and bathy bacteriophytochromes, but in both cases the switch between them is triggered by the photoisomerization of an embedded bilin chromophore. We applied an integrated multiscale strategy of excited-state QM/MM nonadiabatic dynamics and (QM/)MM molecular dynamics simulations with enhanced sampling techniques to the Agrobacterium fabrum bathy phytochrome and compared the results with those obtained for the canonical phytochrome Deinococcus radiodurans. Contrary to what recently suggested, we found that photoactivation in both phytochromes is triggered by the same hula-twist motion of the bilin chromophore. However, only in the bathy phytochrome, the bilin reaches the final rotated structure already in the first intermediate. This allows a reorientation of the binding pocket in a microsecond time scale, which can propagate through the entire protein causing the spine to tilt.

Peptidomimetics design and characterization: Bridging experimental and computer-based approaches

Peptidomimetics, designed to mimic peptide biological activity with more drug-like properties, are increasingly pivotal in medicinal chemistry. They offer enhanced systemic delivery, cell penetration, target specificity, and protection against peptidases when compared to their native peptide counterparts. Already utilized in treating diverse diseases like neurodegenerative disorders, cancer and infectious diseases, their future in medicine seems bright, with many peptidomimetics in clinical trials or development stages. Peptidomimetics are well-suited for addressing disturbed protein-protein interactions (PPIs), which often underlie various pathologies. Structural biology and computational methods like molecular dynamics simulations facilitate rational design, whereas machine learning algorithms accelerate protein structure prediction, enabling efficient drug development. Experimental validation via various spectroscopic, biophysical, and biochemical assays confirms computational predictions and guides further optimization. Peptidomimetics, with their tailored constrained structures, represent a frontier in drug design focused on targeting PPIs. In this overview, we present a comprehensive landscape of peptidomimetics, encompassing perspectives on involvement in pathologies, chemical strategies, and methodologies for their characterization, spanning in silico, in vitro and in cell approaches. With increasing interest from pharmaceutical sectors, peptidomimetics hold promise for revolutionizing therapeutic approaches, marking a new era of precision drug discovery.

Replica Exchange of Expanded Ensembles: A Generalized Ensemble Approach with Enhanced Flexibility and Parallelizability

Generalized ensemble methods such as Hamiltonian replica exchange (HREX) and expanded ensemble (EE) have been shown effective in free energy calculations for various contexts, given their ability to circumvent free energy barriers via nonphysical pathways defined by states with different modified Hamiltonians. However, both HREX and EE methods come with drawbacks, such as limited flexibility in parameter specification or the lack of parallelizability for more complicated applications. To address this challenge, we present the method of replica exchange of expanded ensembles (REXEE), which integrates the principles of HREX and EE methods by periodically exchanging coordinates of EE replicas sampling different yet overlapping sets of alchemical states. With the solvation free energy calculation of anthracene and binding free energy calculation of the CB7-10 binding complex, we show that the REXEE method achieves the same level of accuracy in free energy calculations as the HREX and EE methods, while offering enhanced flexibility and parallelizability. Additionally, we examined REXEE simulations with various setups to understand how different exchange frequencies and replica configurations influence the sampling efficiency in the fixed-weight phase and the weight convergence in the weight-updating phase. The REXEE approach can be further extended to support asynchronous parallelization schemes, allowing looser communications between larger numbers of loosely coupled processors such as cloud computing and therefore promising much more scalable and adaptive executions of alchemical free energy calculations. All algorithms for the REXEE method are available in the Python package ensemble_md, which offers an interface for REXEE simulation management without modifying the source code in GROMACS.

Insights into the mechanism of peptide fibril growth on gold surface

Understanding the formation of β-fibrils over the gold surface is of paramount interest in nano-bio-medicinal Chemistry. The intricate mechanism of self-assembly of neurofibrillogenic peptides and their growth over the gold surface remains elusive, as experiments are limited in unveiling the microscopic dynamic details, in particular, at the early stage of the peptide aggregation. In this work, we carried out equilibrium molecular dynamics and enhanced sampling simulations to elucidate the underlying mechanism of the growth of an amyloid-forming sequence of tau fragments over the gold surface. Our results disclose that the collective intermolecular interactions between the peptide chains and peptides with the gold surface facilitate the peptide adsorption, followed by integration, finally leading to the fibril formation.

Isomerization Pathways of a Mismatched Base Pair of A:8OG in Free Duplex DNA

The A:8OG base pair (bp) is the outcome of DNA replication of the mismatched C:8OG bp. A high A:8OG bp population increases the C/G to A/T transversion mutation, which is responsible for various diseases. MutY is an important enzyme in the error-proof cycle and reverts A:8OG to C:8OG bp by cleaving adenine from the A:8OG bp. Several X-ray crystallography studies have determined the structure of MutY during the lesion scanning and lesion recognition stages. Interestingly, glycosidic bond (χ) angles of A:8OG bp in those two lesion recognition structures were found to differ, which implies that χ-torsion isomerization should occur during the lesion recognition process. In this study, as a first step to understanding this isomerization process, we characterized the intrinsic dynamic features of A:8OG in free DNAs by a free energy landscape simulation at the all-atom level. In this study, four isomerization states were assigned in the order of abundance: Aanti:8OGsyn > Aanti:8OGanti > Asyn:8OGanti ≈ Asyn:8OGsyn. Of these bp states, only 8OG in Asyn:8OGanti was located in the extrahelical space, whereas the purine bases (A and 8OG) in the other bp states remained inside the DNA helix. Also, free energy landscapes showed that the isomerization processes connecting these four bp states proceeded mostly in the intrahelical space via successive single glycosidic bond rotations of either A or 8OG.

Short-Time Infrequent Metadynamics for Improved Kinetics Inference

Infrequent Metadynamics is a popular method to obtain the rates of long time-scale processes from accelerated simulations. The inference procedure is based on rescaling the first-passage times of the Metadynamics trajectories using a bias-dependent acceleration factor. While useful in many cases, it is limited to Poisson kinetics, and a reliable estimation of the unbiased rate requires slow bias deposition and prior knowledge of efficient collective variables. Here, we propose an improved inference scheme, which is based on two key observations: (1) the time-independent rate of Poisson processes can be estimated using short trajectories only. (2) Short trajectories experience minimal bias, and their rescaled first-passage times follow the unbiased distribution even for relatively high deposition rates and suboptimal collective variables. Therefore, by basing the inference procedure on short time scales, we obtain an improved trade-off between speedup and accuracy at no additional computational cost, especially when employing suboptimal collective variables. We demonstrate the improved inference scheme for a model system and two molecular systems.

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.

SWISH-X, an Expanded Approach to Detect Cryptic Pockets in Proteins and at Protein-Protein Interfaces

Protein-protein interactions mediate most molecular processes in the cell, offering a significant opportunity to expand the set of known druggable targets. Unfortunately, targeting these interactions can be challenging due to their typically flat and featureless interaction surfaces, which often change as the complex forms. Such surface changes may reveal hidden (cryptic) druggable pockets. Here, we analyze a set of well-characterized protein-protein interactions harboring cryptic pockets and investigate the predictive power of current computational methods. Based on our observations, we developed a new computational strategy, SWISH-X (SWISH Expanded), which combines the established cryptic pocket identification capabilities of SWISH with the rapid temperature range exploration of OPES MultiThermal. SWISH-X is able to reliably identify cryptic pockets at protein-protein interfaces while retaining its predictive power for revealing cryptic pockets in isolated proteins, such as TEM-1 β-lactamase.

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.

Three-State Diffusion Model of DNA Glycosylase Translocation along Stretched DNA as Revealed by Free Energy Landscapes at the All-Atom Level

DNA glycosylases play key roles in the maintenance of genomic integrity. These enzymes effectively find rare damaged sites in DNA and participate in subsequent base excision repair. Single-molecule and ensemble experiments have revealed key aspects of this damage-site searching mechanism and the involvement of facilitated diffusion. In this study, we describe free energy landscapes of enzyme translocation along nonspecific DNA obtained using a fully atomistic molecular dynamics (MD) simulation of a well-known DNA glycosylase, human 8-oxoguanine DNA glycosylase 1 (hOGG1). Based on an analysis of simulated free energy profiles, we propose a three-state model for the damage-site searching mechanism. In the three states, named the L1, L2, and L3 states, the L1 state is a helical sliding mode in close contact with DNA, whereas the L2 state is a major- or minor-groove tracking mode in loose contact with DNA and the L3 state is a two-dimensional freely diffusing mode during which hOGG1 is somewhat removed from the DNA surface (∼24 Å away from the surface). This three-state model well describes key experimental findings obtained from single-molecule and ensemble experiments and provides a unified molecular picture of the DNA lesion-searching mechanism of hOGG1.

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.

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.

Collective Variable-Based Enhanced Sampling: From Human Learning to Machine Learning

Enhanced-sampling algorithms relying on collective variables (CVs) are extensively employed to study complex (bio)chemical processes that are not amenable to brute-force molecular simulations. The selection of appropriate CVs characterizing the slow movement modes is of paramount importance for reliable and efficient enhanced-sampling simulations. In this Perspective, we first review the application and limitations of CVs obtained from chemical and geometrical intuition. We also introduce path-sampling algorithms, which can identify path-like CVs in a high-dimensional free-energy space. Machine-learning algorithms offer a viable approach to finding suitable CVs by analyzing trajectories from preliminary simulations. We discuss both the performance of machine-learning-derived CVs in enhanced-sampling simulations of experimental models and the challenges involved in applying these CVs to realistic, complex molecular assemblies. Moreover, we provide a prospective view of the potential advancements of machine-learning algorithms for the development of CVs in the field of enhanced-sampling simulations.

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.

Combining stochastic resetting with Metadynamics to speed-up molecular dynamics simulations

Metadynamics is a powerful method to accelerate molecular dynamics simulations, but its efficiency critically depends on the identification of collective variables that capture the slow modes of the process. Unfortunately, collective variables are usually not known a priori and finding them can be very challenging. We recently presented a collective variables-free approach to enhanced sampling using stochastic resetting. Here, we combine the two methods, showing that it can lead to greater acceleration than either of them separately. We also demonstrate that resetting Metadynamics simulations performed with suboptimal collective variables can lead to speedups comparable with those obtained with optimal collective variables. Therefore, applying stochastic resetting can be an alternative to the challenging task of improving suboptimal collective variables, at almost no additional computational cost. Finally, we propose a method to extract unbiased mean first-passage times from Metadynamics simulations with resetting, resulting in an improved tradeoff between speedup and accuracy. This work enables combining stochastic resetting with other enhanced sampling methods to accelerate a broad range of molecular simulations.

Standard Binding Free-Energy Calculations: How Far Are We from Automation?

Recent success stories suggest that in silico protein-ligand binding free-energy calculations are approaching chemical accuracy. However, their widespread application remains limited by the extensive human intervention required, posing challenges for the neophyte. As such, it is critical to develop automated workflows for estimating protein-ligand binding affinities with minimum personal involvement. Key human efforts include setting up and tuning enhanced-sampling or alchemical-transformation algorithms as a preamble to computational binding free-energy estimations. Additionally, preparing input files, bookkeeping, and postprocessing represent nontrivial tasks. In this Perspective, we discuss recent progress in automating standard binding free-energy calculations, featuring the development of adaptive or parameter-free algorithms, standardization of binding free-energy calculation workflows, and the implementation of user-friendly software. We also assess the current state of automated standard binding free-energy calculations and evaluate the limitations of existing methods. Last, we outline the requirements for future algorithms and workflows to facilitate automated free-energy calculations for diverse protein-ligand complexes.

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.

Computational Probing of the Folding Mechanism of Human Telomeric G-Quadruplex DNA

The human telomeric (htel) sequences in the terminal regions of human telomeres form diverse G-quadruplex (GQ) structures. Despite much experimental efforts to elucidate the folding pathways of htel GQ, no comprehensive model of htel GQ folding has been presented. Here, we describe folding pathways of the htel GQ determined by state-of-the-art enhanced sampling molecular dynamics simulation at the all-atom level. Briefly, GQ folding is initiated by the formation of a single-hairpin and then followed by the formation of double-hairpins, which then branch via distinct folding pathways to produce different GQ topologies (antiparallel chair, antiparallel basket, hybrids 1 and 2, and parallel propeller). In addition to these double-hairpin states, three-triad and two-tetrad structures in antiparallel backbone alignment serve as key intermediates that connect the GQ folding and transition between two different GQs.

Vendi sampling for molecular simulations: Diversity as a force for faster convergence and better exploration

Molecular dynamics (MD) is the method of choice for understanding the structure, function, and interactions of molecules. However, MD simulations are limited by the strong metastability of many molecules, which traps them in a single conformation basin for an extended amount of time. Enhanced sampling techniques, such as metadynamics and replica exchange, have been developed to overcome this limitation and accelerate the exploration of complex free energy landscapes. In this paper, we propose Vendi Sampling, a replica-based algorithm for increasing the efficiency and efficacy of the exploration of molecular conformation spaces. In Vendi sampling, replicas are simulated in parallel and coupled via a global statistical measure, the Vendi Score, to enhance diversity. Vendi sampling allows for the recovery of unbiased sampling statistics and dramatically improves sampling efficiency. We demonstrate the effectiveness of Vendi sampling in improving molecular dynamics simulations by showing significant improvements in coverage and mixing between metastable states and convergence of free energy estimates for four common benchmarks, including Alanine Dipeptide and Chignolin.

Perspectives on Ligand/Protein Binding Kinetics Simulations: Force Fields, Machine Learning, Sampling, and User-Friendliness

Computational techniques applied to drug discovery have gained considerable popularity for their ability to filter potentially active drugs from inactive ones, reducing the time scale and costs of preclinical investigations. The main focus of these studies has historically been the search for compounds endowed with high affinity for a specific molecular target to ensure the formation of stable and long-lasting complexes. Recent evidence has also correlated the in vivo drug efficacy with its binding kinetics, thus opening new fascinating scenarios for ligand/protein binding kinetic simulations in drug discovery. The present article examines the state of the art in the field, providing a brief summary of the most popular and advanced ligand/protein binding kinetics techniques and evaluating their current limitations and the potential solutions to reach more accurate kinetic models. Particular emphasis is put on the need for a paradigm change in the present methodologies toward ligand and protein parametrization, the force field problem, characterization of the transition states, the sampling issue, and algorithms' performance, user-friendliness, and data openness.

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.

OneOPES, a Combined Enhanced Sampling Method to Rule Them All

Enhanced sampling techniques have revolutionized molecular dynamics (MD) simulations, enabling the study of rare events and the calculation of free energy differences in complex systems. One of the main families of enhanced sampling techniques uses physical degrees of freedom called collective variables (CVs) to accelerate a system's dynamics and recover the original system's statistics. However, encoding all the relevant degrees of freedom in a limited number of CVs is challenging, particularly in large biophysical systems. Another category of techniques, such as parallel tempering, simulates multiple replicas of the system in parallel, without requiring CVs. However, these methods may explore less relevant high-energy portions of the phase space and become computationally expensive for large systems. To overcome the limitations of both approaches, we propose a replica exchange method called OneOPES that combines the power of multireplica simulations and CV-based enhanced sampling. This method efficiently accelerates the phase space sampling without the need for ideal CVs, extensive parameters fine tuning nor the use of a large number of replicas, as demonstrated by its successful applications to protein-ligand binding and protein folding benchmark systems. Our approach shows promise as a new direction in the development of enhanced sampling techniques for molecular dynamics simulations, providing an efficient and robust framework for the study of complex and unexplored problems.

Understanding Molecular Aggregation of Ligand-Protected Atomically-Precise Metal Nanoclusters

Monolayer-protected atomically precise nanoclusters (MPCs) are an important class of molecules due to their unique structural features and diverse applications, including bioimaging, sensors, and drug carriers. Understanding the atomistic and dynamical details of their self-assembly process is crucial for designing system-specific applications. Here, we applied molecular dynamics and on-the-fly probability-based enhanced sampling simulations to study the aggregation of Au₂₅(pMBA)₁₈ MPCs in aqueous and methanol solutions. The MPCs interact via both hydrogen bonds and π-stacks between the aromatic ligands to form stable dimers, oligomers, and crystals. The dimerization free energy profiles reveal a pivotal role of the ligand charged state and solvent mediating the molecular aggregation. Furthermore, MPCs' ligands exhibit suppressed conformational flexibility in the solid phase due to facile intercluster hydrogen bonds and π-stacks. Our work provides unprecedented molecular-level dynamical details of the aggregation process and conformational dynamics of MPCs ligands in solution and crystalline phases.

Improving Speed and Affordability without Compromising Accuracy: Standard Binding Free-Energy Calculations Using an Enhanced Sampling Algorithm, Multiple-Time Stepping, and Hydrogen Mass Repartitioning

Accurate evaluation of protein-ligand binding free energies in silico is of paramount importance for understanding the mechanisms of biological regulation and providing a theoretical basis for drug design and discovery. Based on a series of atomistic molecular dynamics simulations in an explicit solvent, using well-tempered metadynamics extended adaptive biasing force (WTM-eABF) as an enhanced sampling algorithm, the so-called "geometrical route" offers a rigorous theoretical framework for binding affinity calculations that match experimental values. However, although robust, this strategy remains expensive, requiring substantial computational time to achieve convergence of the simulations. Improving the efficiency of the geometrical route, while preserving its reliability through improved ergodic sampling, is, therefore, highly desirable. In this contribution, having identified the computational bottleneck of the geometrical route, to accelerate the calculations we combine (i) a longer time step for the integration of the equations of motion with hydrogen-mass repartitioning (HMR), and (ii) multiple time-stepping (MTS) for collective-variable and biasing-force evaluation. Altogether, we performed 50 independent WTM-eABF simulations in triplicate for the "physical" separation of the Abl kinase-SH3 domain:p41 complex, following different HMR and MTS schemes, while tuning, in distinct protocols, the parameters of the enhanced-sampling algorithm. To demonstrate the consistency and reliability of the results obtained with the best-performing setups, we carried out quintuple simulations. Furthermore, we demonstrated the transferability of our method to other complexes by triplicating a 200 ns separation simulation of nine chosen protocols for the MDM2-p53:NVP-CGM097 complex. [Holzer et al. J. Med. Chem. 2015, 58, 6348-6358.] Our results, based on an aggregate simulation time of 14.4 μs, allowed an optimal set of parameters to be identified, able to accelerate convergence by a factor of three without any noticeable loss of accuracy.

Alchemical Metadynamics: Adding Alchemical Variables to Metadynamics to Enhance Sampling in Free Energy Calculations

Performing alchemical transformations, in which one molecular system is nonphysically changed to another system, is a popular approach adopted in performing free energy calculations associated with various biophysical processes, such as protein-ligand binding or the transfer of a molecule between environments. While the sampling of alchemical intermediate states in either parallel (e.g., Hamiltonian replica exchange) or serial manner (e.g., expanded ensemble) can bridge the high-probability regions in the configurational space between two end states of interest, alchemical methods can fail in scenarios where the most important slow degrees of freedom in the configurational space are, in large part, orthogonal to the alchemical variable, or if the system gets trapped in a deep basin extending in both the configurational and alchemical space. To alleviate these issues, we propose to use alchemical variables as an additional dimension in metadynamics, making it possible to both sample collective variables and to enhance sampling in free energy calculations simultaneously. In this study, we validate our implementation of "alchemical metadynamics" in PLUMED with test systems and alchemical processes with varying complexities and dimensionalities of collective variable space, including the interconversion between the torsional metastable states of a toy system and the methylation of a nucleoside both in the isolated form and in a duplex. We show that multidimensional alchemical metadynamics can address the challenges mentioned above and further accelerate sampling by introducing configurational collective variables. The method can trivially be combined with other metadynamics-based algorithms implemented in PLUMED. The necessary PLUMED code changes have already been released for general use in PLUMED 2.8.

Skipping the Replica Exchange Ladder with Normalizing Flows

We combine replica exchange (parallel tempering) with normalizing flows, a class of deep generative models. These two sampling strategies complement each other, resulting in an efficient method for sampling molecular systems characterized by rare events, which we call learned replica exchange (LREX). In LREX, a normalizing flow is trained to map the configurations of the fastest-mixing replica into configurations belonging to the target distribution, allowing direct exchanges between the two without the need to simulate intermediate replicas. This can significantly reduce the computational cost compared to standard replica exchange. The proposed method also offers several advantages with respect to Boltzmann generators that directly use normalizing flows to sample the target distribution. We apply LREX to some prototypical molecular dynamics systems, highlighting the improvements over previous methods.

Enhancing sampling with free-energy calculations

In recent years, considerable progress has been made to enhance sampling and help address biological questions, including, but not limited to conformational transitions in biomolecules and protein-ligand reversible binding, hitherto intractable by brute-force computer simulations. Many of these advances result from the development of a palette of methods aimed at exploring rare events through reliable free-energy calculations. The advent of new, often conceptually related methods has also rendered difficult the choice of the best suited option for a given problem. Here, we focus on geometrical transformations and algorithms designed to enhance sampling along adequately chosen progress variables, tracing their theoretical foundations, and showing how they are connected and can be blended together for improved performance.

Extending and validating bubble nucleation rate predictions in a Lennard-Jones fluid with enhanced sampling methods and transition state theory

We calculate bubble nucleation rates in a Lennard-Jones fluid through explicit molecular dynamics simulations. Our approach—based on a recent free energy method (dubbed reweighted Jarzynski sampling), transition state theory, and a simple recrossing correction—allows us to probe a fairly wide range of rates in several superheated and cavitation regimes in a consistent manner. Rate predictions from this approach bridge disparate independent literature studies on the same model system. As such, we find that rate predictions based on classical nucleation theory, direct brute force molecular dynamics simulations, and seeding are consistent with our approach and one another. Published rates derived from forward flux sampling simulations are, however, found to be outliers. This study serves two purposes: First, we validate the reliability of common modeling techniques and extrapolation approaches on a paradigmatic problem in materials science and chemical physics. Second, we further test our highly generic recipe for rate calculations, and establish its applicability to nucleation processes.

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.