Weighted-ensemble Brownian dynamics simulations for protein association reactions

Ensemble Adaptive Sampling Scheme: Identifying an Optimal Sampling Strategy via Policy Ranking

Efficient sampling in biomolecular simulations is critical for accurately capturing the complex dynamic behaviors of biological systems. Adaptive sampling techniques aim to improve efficiency by focusing computational resources on the most relevant regions of the phase space. In this work, we present a framework for identifying the optimal sampling policy through metric-driven ranking. Our approach systematically evaluates the policy ensemble and ranks the policies based on their ability to explore the conformational space effectively. Through a series of biomolecular simulation case studies, we demonstrate that the choice of a different adaptive sampling policy at each round significantly outperforms single policy sampling, leading to faster convergence and improved sampling performance. This approach takes an ensemble of adaptive sampling policies and identifies the optimal policy for the next round based on current data. Beyond presenting this ensemble view of adaptive sampling, we also propose two sampling algorithms that approximate this ranking framework on the fly. The modularity of this framework allows incorporation of any adaptive sampling policy, making it versatile and suitable as a comprehensive adaptive sampling scheme.

Path sampling challenges in large biomolecular systems: RETIS and REPPTIS for ABL-imatinib kinetics

Predicting the kinetics of drug-protein interactions is crucial for understanding drug efficacy, particularly in personalized medicine, where protein mutations can significantly alter drug residence times. This study applies replica exchange transition interface sampling and its partial path variant to investigate the dissociation kinetics of imatinib from Abelson nonreceptor tyrosine kinase (ABL) and mutants relevant to chronic myeloid leukemia therapy. These path sampling methods offer a bias-free alternative to conventional approaches requiring qualitative predefined reaction coordinates. Nevertheless, the complex free energy landscape of ABL-imatinib dissociation presents significant challenges. Multiple metastable states and orthogonal barriers lead to parallel unbinding pathways, complicating convergence in transition interface sampling-based methods. Despite employing computational efficiency strategies such as asynchronous replica exchange, full convergence remained elusive. This work provides a critical assessment of path sampling in high-dimensional biological systems, discussing the need for enhanced initialization strategies, advanced Monte Carlo path generation moves, and machine learning-derived reaction coordinates to improve kinetic predictions of drug dissociation with minimal prior knowledge.

Regulatory Interactions between APOBEC3B N- and C-Terminal Domains

APOBEC3B (A3B) is implicated in DNA mutations that facilitate tumor evolution. Although structures of its individual N- and C-terminal domains (NTD and CTD) have been resolved through X-ray crystallography, the full-length A3B (fl-A3B) structure remains elusive, limiting our understanding of its dynamics and mechanisms. In particular, the APOBEC3B C-terminal domain (A3Bctd) is frequently closed in models and structures. In this study, we built several new models of fl-A3B using integrative structural biology methods and selected a top model for further dynamical investigation. We compared the dynamics of the truncated (A3Bctd) to that of the fl-A3B via conventional and Gaussian accelerated molecular dynamics (MD) simulations. Subsequently, we employed weighted ensemble methods to explore the fl-A3B active site opening mechanism, finding that interactions at the NTD-CTD interface enhance the opening frequency of the fl-A3B active site. Our findings shed light on the structural dynamics and potential druggability of fl-A3B, including observations regarding both the active and allosteric sites, which may offer new avenues for therapeutic intervention in cancer.

Unsupervised Learning of Progress Coordinates during Weighted Ensemble Simulations: Application to NTL9 Protein Folding

A major challenge for many rare-event sampling strategies is the identification of progress coordinates that capture the slowest relevant motions. Machine-learning methods that can identify progress coordinates in an unsupervised manner have therefore been of great interest to the simulation community. Here, we developed a general method for identifying progress coordinates "on-the-fly" during weighted ensemble (WE) rare-event sampling via deep learning (DL) of outliers among sampled conformations. Our method identifies outliers in a latent space model of the system's sampled conformations that is periodically trained using a convolutional variational autoencoder. As a proof of principle, we applied our DL-enhanced WE method to simulate the NTL9 protein folding process. To enable rapid tests, our simulations propagated discrete-state synthetic molecular dynamics trajectories using a generative, fine-grained Markov state model. Results revealed that our on-the-fly DL of outliers enhanced the efficiency of WE by >3-fold in estimating the folding rate constant. Our efforts are a significant step forward in the unsupervised learning of slow coordinates during rare event sampling.

A Reflection on the Use of Molecular Simulation to Respond to SARS-CoV-2 Pandemic Threats

Molecular simulations play important roles in understanding the lifecycle of the SARS-CoV-2 virus and contribute to the design and development of antiviral agents and diagnostic tests for COVID. Here, we discuss the insights that such simulations have provided and the challenges involved, focusing on the SARS-CoV-2 main protease (Mpro) and the spike glycoprotein. Mpro is the leading target for antivirals, while the spike glycoprotein is the target for vaccine design. Finally, we reflect on lessons from this pandemic for the simulation community. Data sharing initiatives and collaborations across the international research community contributed to advancing knowledge and should be built on to help in future pandemics and other global challenges such as antimicrobial resistance.

Markov State Models with Weighted Ensemble Simulation: How to Eliminate the Trajectory Merging Bias

The weighted ensemble (WE) algorithm is gaining popularity as a rare event method for studying long timescale processes with molecular dynamics. WE is particularly useful for determining kinetic properties, such as rates of protein (un)folding and ligand (un)binding, where transition rates can be calculated from the flux of trajectories into a target basin of interest. However, this flux depends exponentially on the number of splitting events that a given trajectory experiences before reaching the target state and can vary by orders of magnitude between WE replicates. Markov state models (MSMs) are helpful tools to aggregate information across multiple WE simulations and have previously been shown to provide more accurate transition rates than WE alone. Discrete-time MSMs are models that coarsely describe the evolution of the system from one discrete state to the next using a discrete lag time, τ. When an MSM is built using conventional MD data, longer values of τ typically provide more accurate results. Combining WE simulations with Markov state modeling presents some additional challenges, especially when using a value of τ that exceeds the lag time between resampling steps in the WE algorithm, τWE. Here, we identify a source of bias that occurs when τ > τWE, which we refer to as "merging bias". We also propose an algorithm to eliminate the merging bias, which results in merging bias-corrected MSMs, or "MBC-MSMs". Using a simple model system, as well as a complex biomolecular example, we show that MBC-MSMs significantly outperform both τ = τWE MSMs and uncorrected MSMs at longer lag times.

Molecular simulations reveal intricate coupling between agonist-bound β-adrenergic receptors and G protein

G protein-coupled receptors (GPCRs) and G proteins transmit signals from hormones and neurotransmitters across cell membranes, initiating downstream signaling and modulating cellular behavior. Using advanced computer modeling and simulation, we identified atomistic-level structural, dynamic, and energetic mechanisms of norepinephrine (NE) and stimulatory G protein (Gs) interactions with β-adrenergic receptors (βARs), crucial GPCRs for heart function regulation and major drug targets. Our analysis revealed distinct binding behaviors of NE within β1AR and β2AR despite identical orthosteric binding pockets. β2AR had an additional binding site, explaining variations in NE binding affinities. Simulations showed significant differences in NE dissociation pathways and receptor interactions with the Gs. β1AR binds Gs more strongly, while β2AR induces greater conformational changes in the α subunit of Gs. Furthermore, GTP and GDP binding to Gs may disrupt coupling between NE and βAR, as well as between βAR and Gs. These findings may aid in designing precise βAR-targeted drugs.

Experiment-Guided Refinement of Milestoning Network

Milestoning is an efficient method for calculating rare event kinetics by constructing a continuous-time kinetic network that connects the reactant and product states. Its accuracy depends on both the quality of the underlying force fields and the trajectory sampling. The sampling error can be effectively controlled through various methods. However, the force fields are often not accurate enough, leading to quantitative discrepancies between simulations and experimental data. To address this challenge, we present a refinement approach for Milestoning network based on the maximum caliber (MaxCal), a general variational principle for dynamical systems, to combine simulations and experimental data. The Kullback-Leibler divergence rate between two Milestoning networks is analytically evaluated and minimized as the loss function. Meanwhile, experimental thermodynamic (equilibrium constants) and kinetic (rate constants) data are incorporated as constraints. The use of MaxCal implies that the refined kinetic network is minimally perturbed from the original one while satisfying the experimental constraints. The refined network is expected to align better with available experimental data. The refinement approach is demonstrated using the binding and unbinding dynamics of a series of six small molecule ligands for the model host system, β-cyclodextrin.

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.

WeTICA: A directed search weighted ensemble based enhanced sampling method to estimate rare event kinetics in a reduced dimensional space

Estimating rare event kinetics from molecular dynamics simulations is a non-trivial task despite the great advances in enhanced sampling methods. Weighted Ensemble (WE) simulation, a special class of enhanced sampling techniques, offers a way to directly calculate kinetic rate constants from biased trajectories without the need to modify the underlying energy landscape using bias potentials. Conventional WE algorithms use different binning schemes to partition the collective variable (CV) space separating the two metastable states of interest. In this work, we have developed a new "binless" WE simulation algorithm to bypass the hurdles of optimizing binning procedures. Our proposed protocol (WeTICA) uses a low-dimensional CV space to drive the WE simulation toward the specified target state. We have applied this new algorithm to recover the unfolding kinetics of three proteins: (A) TC5b Trp-cage mutant, (B) TC10b Trp-cage mutant, and (C) Protein G, with unfolding times spanning the range between 3 and 40 μs using projections along predefined fixed Time-lagged Independent Component Analysis (TICA) eigenvectors as CVs. Calculated unfolding times converge to the reported values with good accuracy with more than one order of magnitude less cumulative WE simulation time than the unfolding time scales with or without a priori knowledge of the CVs that can capture unfolding. Our algorithm can be used with other linear CVs, not limited to TICA. Moreover, the new walker selection criteria for resampling employed in this algorithm can be used on more sophisticated nonlinear CV space for further improvements of binless WE methods.

Weighted-ensemble network simulations of the susceptible-infected-susceptible model of epidemics

The presence of erratic or unstable paths in standard kinetic Monte Carlo simulations significantly undermines the accurate simulation and sampling of transition pathways. While typically reliable methods, such as the Gillespie algorithm, are employed to simulate such paths, they encounter challenges in efficiently identifying rare events due to their sequential nature and reliance on exact Monte Carlo sampling. In contrast, the weighted-ensemble method effectively samples rare events and accelerates the exploration of complex reaction pathways by distributing computational resources among multiple replicas, where each replica is assigned a weight reflecting its importance, and evolves independently from the others. Here, we implement the highly efficient and robust weighted-ensemble method to model susceptible-infected-susceptible dynamics on large heterogeneous population networks, and explore the interplay between stochasticity and contact heterogeneity, which ultimately gives rise to disease clearance. Studying a wide variety of networks characterized by fat-tailed asymmetric degree distributions, we are able to compute the mean time to extinction and quasistationary distribution around it in previously inaccessible parameter regimes.

Regulatory interactions between APOBEC3B N- and C-terminal domains

APOBEC3B (A3B) is implicated in DNA mutations that facilitate tumor evolution. Although structures of its individual N- and C-terminal domains (NTD and CTD) have been resolved through X-ray crystallography, the full-length A3B (fl-A3B) structure remains elusive, limiting understanding of its dynamics and mechanisms. In particular, the APOBEC3B C-terminal domain (A3Bctd) active site is frequently closed in models and structures. In this study, we built several new models of fl-A3B using integrative structural biology methods and selected a top model for further dynamical investigation. We compared dynamics of the truncated (A3Bctd) to the fl-A3B via conventional and Gaussian accelerated molecular dynamics (MD) simulations. Subsequently, we employed weighted ensemble methods to explore the fl-A3B active site opening mechanism, finding that interactions at the NTD-CTD interface enhance the opening frequency of the fl-A3B active site. Our findings shed light on the structural dynamics of fl-A3B, which may offer new avenues for therapeutic intervention in cancer.

Augmenting Human Expertise in Weighted Ensemble Simulations through Deep Learning-Based Information Bottleneck

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 chignolin 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.

Augmenting Human Expertise in Weighted Ensemble Simulations through Deep Learning based Information Bottleneck

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.

Exploration of Cryptic Pockets Using Enhanced Sampling Along Normal Modes: A Case Study of KRAS G12D

Identification of cryptic pockets has the potential to open new therapeutic opportunities by discovering ligand binding sites that remain hidden in static apo structures of a target protein. Moreover, allosteric cryptic pockets can become valuable for designing target-selective ligands when the natural ligand binding sites are conserved in variants of a protein. For example, before an allosteric cryptic pocket was discovered, KRAS was considered undruggable due to its smooth surface and conservation of the GDP/GTP binding pocket across the wild type and oncogenic isoforms. Recent identification of the Switch-II cryptic pocket in the KRASG12C mutant and FDA approval of anticancer drugs targeting this site underscores the importance of cryptic pockets in solving pharmaceutical challenges. Here, we present a newly developed approach for the exploration of cryptic pockets using weighted ensemble molecular dynamics simulations with inherent normal modes as progress coordinates applied to the wild type KRAS and the G12D mutant. We performed extensive all-atomic simulations (>400 μs) with and without several cosolvents (xenon, ethanol, benzene), and analyzed trajectories using three distinct methods to search for potential binding pockets. These methods have been applied as a proof-of-concept to KRAS and have shown they can predict known cryptic binding sites. Furthermore, we performed ligand-binding simulations of a known inhibitor (MRTX1133) to shed light on the nature of cryptic pockets in KRASG12D and the role of conformational selection vs induced-fit mechanism in the formation of these cryptic pockets.

Insights into Translocation of Arginine-Rich Cell-Penetrating Peptides across a Model Membrane

It is well-known that membrane deformation and water pores contribute to the spontaneous translocation of arginine-rich cell-penetrating peptides (CPPs). We confirm this through the observation of the spontaneous translocation of single R9 (nona-arginine) and Tat (48-60) peptides across a model membrane using the weighted ensemble (WE) method within all-atom molecular dynamics (MD) simulations. Furthermore, we demonstrate that membrane deformation and the presence of a water pore reduce the effective charge of the CPP and the bending rigidity of the model membrane during translocation. We find that R9 disturbs the model membrane more than Tat (48-60), leading to more efficient translocation of R9 across the model membrane.

Transition-path sampling for run-and-tumble particles

We elaborate and validate a generalization of the renowned transition-path-sampling algorithm for a paradigmatic model of active particles, namely, the run-and-tumble particles. Notwithstanding the nonequilibrium character of these particles, we show how the consequent lack of the microscopical reversibility property, which is usually required by transition-path sampling, can be circumvented by identifying reasonable backward dynamics with a well-defined path-probability density. Our method is then applied to characterize the structure and kinetics of rare transition pathways undergone by run-and-tumble particles having to cross a potential barrier in order to find a target.

Gating residues govern ligand unbinding kinetics from the buried cavity in HIF-2α PAS-B

While transcription factors have been generally perceived as "undruggable," an exception is the HIF-2 hypoxia-inducible transcription factor, which contains an internal cavity that is sufficiently large to accommodate a range of small-molecules, including the therapeutically used inhibitor belzutifan. Given the relatively long ligand residence times of these small molecules and the lack of any experimentally observed pathway connecting the cavity to solvent, there has been great interest in understanding how these drug ligands exit the buried receptor cavity. Here, we focus on the relevant PAS-B domain of hypoxia-inducible factor 2α (HIF-2α) and examine how one such small molecule (THS-017) exits from the buried cavity within this domain on the seconds-timescale using atomistic simulations and ZZ-exchange NMR. To enable the simulations, we applied the weighted ensemble path sampling strategy, which generates continuous pathways for a rare-event process [e.g., ligand (un)binding] with rigorous kinetics in orders of magnitude less computing time compared to conventional simulations. Results reveal the formation of an encounter complex intermediate and two distinct classes of pathways for ligand exit. Based on these pathways, we identified two pairs of conformational gating residues in the receptor: one for the major class (N288 and S304) and another for the minor class (L272 and M309). ZZ-exchange NMR validated the kinetic importance of N288 for ligand unbinding. Our results provide an ideal simulation dataset for rational manipulation of ligand unbinding kinetics.

A Deep Learning-Driven Sampling Technique to Explore the Phase Space of an RNA Stem-Loop

The folding and unfolding of RNA stem-loops are critical biological processes; however, their computational studies are often hampered by the ruggedness of their folding landscape, necessitating long simulation times at the atomistic scale. Here, we adapted DeepDriveMD (DDMD), an advanced deep learning-driven sampling technique originally developed for protein folding, to address the challenges of RNA stem-loop folding. Although tempering- and order parameter-based techniques are commonly used for similar rare-event problems, the computational costs or the need for a priori knowledge about the system often present a challenge in their effective use. DDMD overcomes these challenges by adaptively learning from an ensemble of running MD simulations using generic contact maps as the raw input. DeepDriveMD enables on-the-fly learning of a low-dimensional latent representation and guides the simulation toward the undersampled regions while optimizing the resources to explore the relevant parts of the phase space. We showed that DDMD estimates the free energy landscape of the RNA stem-loop reasonably well at room temperature. Our simulation framework runs at a constant temperature without external biasing potential, hence preserving the information on transition rates, with a computational cost much lower than that of the simulations performed with external biasing potentials. We also introduced a reweighting strategy for obtaining unbiased free energy surfaces and presented a qualitative analysis of the latent space. This analysis showed that the latent space captures the relevant slow degrees of freedom for the RNA folding problem of interest. Finally, throughout the manuscript, we outlined how different parameters are selected and optimized to adapt DDMD for this system. We believe this compendium of decision-making processes will help new users adapt this technique for the rare-event sampling problems of their interest.

System size effects on the free energy landscapes from molecular dynamics of phase-separating bilayers

The "lipid raft" hypothesis proposes that cell membranes contain distinct domains of varying lipid compositions, where "rafts" of ordered lipids and cholesterol coexist with disordered lipid regions. Experimental and theoretical phase diagrams of model membranes have revealed multiple coexisting phases. Molecular dynamics (MD) simulations can also capture spontaneous phase separation of bilayers. However, these methods merely determine the sign of the free energy change upon phase separation-whether or not it is favorable-but not the amplitude. Recently, we developed a workflow to compute the free energy of phase separation from MD simulations using the weighted ensemble method. However, while theoretical treatments generally focus on infinite systems and experimental measurements on mesoscopic to macroscopic systems, MD simulations are comparatively small. Therefore, if we are to put the results of these calculations into the appropriate context, we need to understand the effects the finite size of the simulation has on the computed free energy landscapes. In this study, we investigate this phenomenon by computing free energy profiles for a model phase-separating system as a function of system size, ranging from 324 to 10 110 lipids. The results suggest that, within the limits of statistical uncertainty, bulk-like behavior emerges once the systems contain roughly 4000 lipids.

Rare-Event Sampling using a Reinforcement Learning-Based Weighted Ensemble Method

Despite the power of path sampling strategies in enabling simulations of rare events, such strategies have not reached their full potential. A common challenge that remains is the identification of a progress coordinate that captures the slow relevant motions of a rare event. Here we have developed a weighted ensemble (WE) path sampling strategy that exploits reinforcement learning to automatically identify an effective progress coordinate among a set of potential coordinates during a simulation. We apply our WE strategy with reinforcement learning to three benchmark systems: (i) an egg carton-shaped toy potential, (ii) an S-shaped toy potential, and (iii) a dimer of the HIV-1 capsid protein (C-terminal domain). To enable rapid testing of the latter system at the atomic level, we employed discrete-state synthetic molecular dynamics trajectories using a generative, fine-grained Markov state model that was based on extensive conventional simulations. Our results demonstrate that using concepts from reinforcement learning with a weighted ensemble of trajectories automatically identifies relevant progress co-ordinates among multiple candidates at a given time during a simulation. Due to the rigorous weighting of trajectories, the simulations maintain rigorous kinetics.

Molecular permeation through large pore channels: computational approaches and insights

Computational methods such as molecular dynamics (MD) have illuminated how single-atom ions permeate membrane channels and how selectivity among them is achieved. Much less is understood about molecular permeation through eukaryotic channels that mediate the flux of small molecules (e.g. connexins, pannexins, LRRC8s, CALHMs). Here we describe computational methods that have been profitably employed to explore the movements of molecules through wide pores, revealing mechanistic insights, guiding experiments, and suggesting testable hypotheses. This review illustrates MD techniques such as voltage-driven flux, potential of mean force, and mean first-passage-time calculations, as applied to molecular permeation through wide pores. These techniques have enabled detailed and quantitative modeling of molecular interactions and movement of permeants at the atomic level. We highlight novel contributors to the transit of molecules through these wide pathways. In particular, the flexibility and anisotropic nature of permeant molecules, coupled with the dynamics of pore-lining residues, lead to bespoke permeation dynamics. As more eukaryotic large-pore channel structures and functional data become available, these insights and approaches will be important for understanding the physical principles underlying molecular permeation and as guides for experimental design.

AI in cellular engineering and reprogramming

During the last decade, artificial intelligence (AI) has increasingly been applied in biophysics and related fields, including cellular engineering and reprogramming, offering novel approaches to understand, manipulate, and control cellular function. The potential of AI lies in its ability to analyze complex datasets and generate predictive models. AI algorithms can process large amounts of data from single-cell genomics and multiomic technologies, allowing researchers to gain mechanistic insights into the control of cell identity and function. By integrating and interpreting these complex datasets, AI can help identify key molecular events and regulatory pathways involved in cellular reprogramming. This knowledge can inform the design of precision engineering strategies, such as the development of new transcription factor and signaling molecule cocktails, to manipulate cell identity and drive authentic cell fate across lineage boundaries. Furthermore, when used in combination with computational methods, AI can accelerate and improve the analysis and understanding of the intricate relationships between genes, proteins, and cellular processes. In this review article, we explore the current state of AI applications in biophysics with a specific focus on cellular engineering and reprogramming. Then, we showcase a couple of recent applications where we combined machine learning with experimental and computational techniques. Finally, we briefly discuss the challenges and prospects of AI in cellular engineering and reprogramming, emphasizing the potential of these technologies to revolutionize our ability to engineer cells for a variety of applications, from disease modeling and drug discovery to regenerative medicine and biomanufacturing.

Identification of small-molecule ligand-binding sites on and in the ARNT PAS-B domain

Transcription factors are challenging to target with small-molecule inhibitors due to their structural plasticity and lack of catalytic sites. Notable exceptions include naturally ligand-regulated transcription factors, including our prior work with the hypoxia-inducible factor (HIF)-2 transcription factor, showing that small-molecule binding within an internal pocket of the HIF-2α Per-Aryl hydrocarbon Receptor Nuclear Translocator (ARNT)-Sim (PAS)-B domain can disrupt its interactions with its dimerization partner, ARNT. Here, we explore the feasibility of targeting small molecules to the analogous ARNT PAS-B domain itself, potentially opening a promising route to modulate several ARNT-mediated signaling pathways. Using solution NMR fragment screening, we previously identified several compounds that bind ARNT PAS-B and, in certain cases, antagonize ARNT association with the transforming acidic coiled-coil containing protein 3 transcriptional coactivator. However, these ligands have only modest binding affinities, complicating characterization of their binding sites. We address this challenge by combining NMR, molecular dynamics simulations, and ensemble docking to identify ligand-binding "hotspots" on and within the ARNT PAS-B domain. Our data indicate that the two ARNT/transforming acidic coiled-coil containing protein 3 inhibitors, KG-548 and KG-655, bind to a β-sheet surface implicated in both HIF-2 dimerization and coactivator recruitment. Furthermore, while KG-548 binds exclusively to the β-sheet surface, KG-655 can additionally bind within a water-accessible internal cavity in ARNT PAS-B. Finally, KG-279, while not a coactivator inhibitor, exemplifies ligands that preferentially bind only to the internal cavity. All three ligands promoted ARNT PAS-B homodimerization, albeit to varying degrees. Taken together, our findings provide a comprehensive overview of ARNT PAS-B ligand-binding sites and may guide the development of more potent coactivator inhibitors for cellular and functional studies.

BAD-NEUS: Rapidly converging trajectory stratification

An issue for molecular dynamics simulations is that events of interest often involve timescales that are much longer than the simulation time step, which is set by the fastest timescales of the model. Because of this timescale separation, direct simulation of many events is prohibitively computationally costly. This issue can be overcome by aggregating information from many relatively short simulations that sample segments of trajectories involving events of interest. This is the strategy of Markov state models (MSMs) and related approaches, but such methods suffer from approximation error because the variables defining the states generally do not capture the dynamics fully. By contrast, once converged, the weighted ensemble (WE) method aggregates information from trajectory segments so as to yield unbiased estimates of both thermodynamic and kinetic statistics. Unfortunately, errors decay no faster than unbiased simulation in WE as originally formulated and commonly deployed. Here, we introduce a theoretical framework for describing WE that shows that the introduction of an approximate stationary distribution on top of the stratification, as in nonequilibrium umbrella sampling (NEUS), accelerates convergence. Building on ideas from MSMs and related methods, we generalize the NEUS approach in such a way that the approximation error can be reduced systematically. We show that the improved algorithm can decrease the simulation time required to achieve the desired precision by orders of magnitude.

Simulation-driven design of stabilized SARS-CoV-2 spike S2 immunogens

The full-length prefusion-stabilized SARS-CoV-2 spike (S) is the principal antigen of COVID-19 vaccines. Vaccine efficacy has been impacted by emerging variants of concern that accumulate most of the sequence modifications in the immunodominant S1 subunit. S2, in contrast, is the most evolutionarily conserved region of the spike and can elicit broadly neutralizing and protective antibodies. Yet, S2's usage as an alternative vaccine strategy is hampered by its general instability. Here, we use a simulation-driven approach to design S2-only immunogens stabilized in a closed prefusion conformation. Molecular simulations provide a mechanistic characterization of the S2 trimer's opening, informing the design of tryptophan substitutions that impart kinetic and thermodynamic stabilization. Structural characterization via cryo-EM shows the molecular basis of S2 stabilization in the closed prefusion conformation. Informed by molecular simulations and corroborated by experiments, we report an engineered S2 immunogen that exhibits increased protein expression, superior thermostability, and preserved immunogenicity against sarbecoviruses.

WEDAP: A Python Package for Streamlined Plotting of Molecular Simulation Data

Given the growing interest in path sampling methods for extending the time scales of molecular dynamics (MD) simulations, there has been great interest in software tools that streamline the generation of plots for monitoring the progress of large-scale simulations. Here, we present the WEDAP Python package for simplifying the analysis of data generated from either conventional MD simulations or the weighted ensemble (WE) path sampling method, as implemented in the widely used WESTPA software package. WEDAP facilitates (i) the parsing of WE simulation data stored in highly compressed, hierarchical HDF5 files and (ii) incorporates trajectory weights from WE simulations into all generated plots. Our Python package consists of multiple user-friendly interfaces: a command-line interface, a graphical user interface, and a Python application programming interface. We demonstrate the plotting features of WEDAP through a series of examples using data from WE and conventional MD simulations that focus on the HIV-1 capsid protein's C-terminal domain dimer as a showcase system. The source code for WEDAP is freely available on GitHub at https://github.com/chonglab-pitt/wedap.

Structures, dynamics, complexes, and functions: From classic computation to artificial intelligence

Computational approaches can provide highly detailed insight into the molecular recognition processes that underlie drug binding, the assembly of protein complexes, and the regulation of biological functional processes. Classical simulation methods can bridge a wide range of length- and time-scales typically involved in such processes. Lately, automated learning and artificial intelligence methods have shown the potential to expand the reach of physics-based approaches, ushering in the possibility to model and even design complex protein architectures. The synergy between atomistic simulations and AI methods is an emerging frontier with a huge potential for advances in structural biology. Herein, we explore various examples and frameworks for these approaches, providing select instances and applications that illustrate their impact on fundamental biomolecular problems.

Kinetic Network in Milestoning: Clustering, Reduction, and Transition Path Analysis

We present a reduction of the Milestoning (ReM) algorithm to analyze the high-dimensional Milestoning kinetic network. The algorithm reduces the Milestoning network to low dimensions but preserves essential kinetic information, such as local residence time, exit time, and mean first passage time between any two states. This is achieved in three steps. First, nodes (milestones) in the high-dimensional Milestoning network are grouped into clusters based on the metastability identified by an auxiliary continuous-time Markov chain. Our clustering method is applicable not only to time-reversible networks but also to nonreversible networks generated from practical simulations with statistical fluctuations. Second, a reduced network is established via network transformation, containing only the core sets of clusters as nodes. Finally, transition pathways are analyzed in the reduced network based on the transition path theory. The algorithm is illustrated using a toy model and a solvated alanine dipeptide in two and four dihedral angles.

Combining Transition Path Sampling with Data-Driven Collective Variables through a Reactivity-Biased Shooting Algorithm

Rare event sampling is a central problem in modern computational chemistry research. Among the existing methods, transition path sampling (TPS) can generate unbiased representations of reaction processes. However, its efficiency depends on the ability to generate reactive trial paths, which in turn depends on the quality of the shooting algorithm used. We propose a new algorithm based on the shooting success rate, i.e., reactivity, measured as a function of a reduced set of collective variables (CVs). These variables are extracted with a machine learning approach directly from TPS simulations, using a multitask objective function. Iteratively, this workflow significantly improves the shooting efficiency without any prior knowledge of the process. In addition, the optimized CVs can be used with biased enhanced sampling methodologies to accurately reconstruct the free energy profiles. We tested the method on three different systems: a two-dimensional toy model, conformational transitions of alanine dipeptide, and hydrolysis of acetyl chloride in bulk water. In the latter, we integrated our workflow with an active learning scheme to learn a reactive machine learning-based potential, which allowed us to study the mechanism and free energy profile with an ab initio-like accuracy.

Identification of Small Molecule Ligand Binding Sites On and In the ARNT PAS-B Domain

Transcription factors are generally challenging to target with small molecule inhibitors due to their structural plasticity and lack of catalytic sites. Notable exceptions include several naturally ligand-regulated transcription factors, including our prior work with the heterodimeric HIF-2 transcription factor which showed that small molecule binding within an internal pocket of the HIF-2α PAS-B domain can disrupt its interactions with its dimerization partner, ARNT. Here, we explore the feasibility of similarly targeting small molecules to the analogous ARNT PAS-B domain itself, potentially opening a promising route to simultaneously modulate several ARNT-mediated signaling pathways. Using solution NMR screening of an in-house fragment library, we previously identified several compounds that bind ARNT PAS-B and, in certain cases, antagonize ARNT association with the TACC3 transcriptional coactivator. However, these ligands have only modest binding affinities, complicating characterization of their binding sites. We address this challenge by combining NMR, MD simulations, and ensemble docking to identify ligand-binding 'hotspots' on and within the ARNT PAS-B domain. Our data indicate that the two ARNT/TACC3 inhibitors, KG-548 and KG-655, bind to a β-sheet surface implicated in both HIF-2 dimerization and coactivator recruitment. Furthermore, while KG-548 binds exclusively to the β-sheet surface, KG-655 can additionally bind within a water-accessible internal cavity in ARNT PAS-B. Finally, KG-279, while not a coactivator inhibitor, exemplifies ligands that preferentially bind only to the internal cavity. All three ligands promoted ARNT PAS-B homodimerization, albeit to varying degrees. Taken together, our findings provide a comprehensive overview of ARNT PAS-B ligand-binding sites and may guide the development of more potent coactivator inhibitors for cellular and functional studies.

Dynamical Reweighting for Biased Rare Event Simulations

Dynamical reweighting techniques aim to recover the correct molecular dynamics from a simulation at a modified potential energy surface. They are important for unbiasing enhanced sampling simulations of molecular rare events. Here, we review the theoretical frameworks of dynamical reweighting for modified potentials. Based on an overview of kinetic models with increasing level of detail, we discuss techniques to reweight two-state dynamics, multistate dynamics, and path integrals. We explore the natural link to transition path sampling and how the effect of nonequilibrium forces can be reweighted. We end by providing an outlook on how dynamical reweighting integrates with techniques for optimizing collective variables and with modern potential energy surfaces.

Mechanism of Daunomycin Intercalation into DNA from Enhanced Sampling Simulations

Daunomycin is a widely used anticancer drug, yet the mechanism underlying how it binds to DNA remains contested. 469 all-atom trajectories of daunomycin binding to the DNA oligonucleotide d(GCG CAC GTG CGC) were collected using weighted ensemble (WE)-enhanced sampling. Mechanistic insights were revealed through analysis of the ensemble of trajectories. Initially, the binding process involves a ubiquitous hydrogen bond between the DNA backbone and the NH3+ group on daunomycin. During the binding process, most trajectories exhibited similar structural changes to DNA, including DNA base pair rise, bending, and minor groove width changes. Variability within the ensemble of binding trajectories illuminates differences in the orientation of daunomycin as it initially intercalates; around 10% of trajectories needed minimal rearrangement from intercalation to reaching the fully bound configuration, whereas most needed an additional 1-5 ns to rearrange. The results here emphasize the utility of generating an ensemble of trajectories to discern biomolecular binding mechanisms.

WEDAP: A Python Package for Streamlined Plotting of Molecular Simulation Data

Given the growing interest in path sampling methods for extending the timescales of molecular dynamics (MD) simulations, there has been great interest in software tools that streamline the generation of plots for monitoring the progress of large-scale simulations. Here, we present the WEDAP Python package for simplifying the analysis of data generated from either conventional MD simulations or the weighted ensemble (WE) path sampling method, as implemented in the widely used WESTPA software package. WEDAP facilitates (i) the parsing of WE simulation data stored in highly compressed, hierarchical HDF5 files, and (ii) incorporates trajectory weights from WE simulations into all generated plots. Our Python package consists of multiple user-friendly interfaces: a command-line interface, a graphical user interface, and a Python application programming interface. We demonstrate the plotting features of WEDAP through a series of examples using data from WE and conventional MD simulations that focus on the HIV-1 capsid protein C-terminal domain dimer as a showcase system. The source code for WEDAP is freely available on GitHub at https://github.com/chonglab-pitt/wedap .

Quantitative Prediction of Dissociation Rates of PYK2 Ligands Using Umbrella Sampling and Milestoning

We used umbrella sampling and the milestoning simulation method to study the dissociation of multiple ligands from protein kinase PYK2. The activation barriers obtained from the potential of mean force of the umbrella sampling simulations correlated well with the experimental dissociation rates. Using the zero-temperature string method, we obtained optimized paths along the free-energy surfaces for milestoning simulations of three ligands with a similar structure. The milestoning simulations gave an absolute dissociation rate within 2 orders of magnitude of the experimental value for two ligands but at least 3 orders of magnitude too high for the third. Despite the similarity in their structures, the ligands took different pathways to exit from the binding site of PYK2, making contact with different sets of residues. In addition, the protein experienced different conformational changes for dissociation of the three ligands.

Revisiting Textbook Azide-Clock Reactions: A "Propeller-Crawling" Mechanism Explains Differences in Rates

An ongoing challenge to chemists is the analysis of pathways and kinetics for chemical reactions in solution, including transient structures between the reactants and products that are difficult to resolve using laboratory experiments. Here, we enabled direct molecular dynamics simulations of a textbook series of chemical reactions on the hundreds of ns to μs time scale using the weighted ensemble (WE) path sampling strategy with hybrid quantum mechanical/molecular mechanical (QM/MM) models. We focused on azide-clock reactions involving addition of an azide anion to each of three long-lived trityl cations in an acetonitrile-water solvent mixture. Results reveal a two-step mechanism: (1) diffusional collision of reactants to form an ion-pair intermediate; (2) "activation" or rearrangement of the intermediate to the product. Our simulations yield not only reaction rates that are within error of experiment but also rates for individual steps, indicating the activation step as rate-limiting for all three cations. Further, the trend in reaction rates is due to dynamical effects, i.e., differing extents of the azide anion "crawling" along the cation's phenyl-ring "propellers" during the activation step. Our study demonstrates the power of analyzing pathways and kinetics to gain insights on reaction mechanisms, underscoring the value of including WE and other related path sampling strategies in the modern toolbox for chemists.

PaCS-Toolkit: Optimized Software Utilities for Parallel Cascade Selection Molecular Dynamics (PaCS-MD) Simulations and Subsequent Analyses

Parallel cascade selection molecular dynamics (PaCS-MD) is an enhanced conformational sampling method conducted as a "repetition of time leaps in parallel worlds", comprising cycles of multiple molecular dynamics (MD) simulations performed in parallel and selection of the initial structures of MDs for the next cycle. We developed PaCS-Toolkit, an optimized software utility enabling the use of different MD software and trajectory analysis tools to facilitate the execution of the PaCS-MD simulation and analyze the obtained trajectories, including the preparation for the subsequent construction of the Markov state model. PaCS-Toolkit is coded with Python, is compatible with various computing environments, and allows for easy customization by editing the configuration file and specifying the MD software and analysis tools to be used. We present the software design of PaCS-Toolkit and demonstrate applications of PaCS-MD variations: original targeted PaCS-MD to peptide folding; rmsdPaCS-MD to protein domain motion; and dissociation PaCS-MD to ligand dissociation from adenosine A2A receptor.

Many Roles of Carbohydrates: A Computational Spotlight on the Coronavirus S Protein Binding

Glycosylation is one of the post-translational modifications with more than 50% of human proteins being glycosylated. The exact nature and chemical composition of glycans are inaccessible to X-ray or cryo-electron microscopy imaging techniques. Therefore, computational modeling studies and molecular dynamics must be used as a "computational microscope". The spike (S) protein of SARS-CoV-2 is heavily glycosylated, and a few glycans play a more functional role "beyond shielding". In this mini-review, we discuss computational investigations of the roles of specific S-protein and ACE2 glycans in the overall ACE2-S protein binding. We highlight different functions of specific glycans demonstrated in myriad computational models and simulations in the context of the SARS-CoV-2 virus binding to the receptor. We also discuss interactions between glycocalyx and the S protein, which may be utilized to design prophylactic polysaccharide-based therapeutics targeting the S protein. In addition, we underline the recent emergence of coronavirus variants and their impact on the S protein and its glycans.

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.

Effects of altered backbone composition on the folding kinetics and mechanism of an ultrafast-folding protein

Sequence-encoded protein folding is a ubiquitous biological process that has been successfully engineered in a range of oligomeric molecules with artificial backbone chemical connectivity. A remarkable aspect of protein folding is the contrast between the rapid rates at which most sequences in nature fold and the vast number of conformational states possible in an unfolded chain with hundreds of rotatable bonds. Research efforts spanning several decades have sought to elucidate the fundamental chemical principles that dictate the speed and mechanism of natural protein folding. In contrast, little is known about how protein mimetic entities transition between an unfolded and folded state. Here, we report effects of altered backbone connectivity on the folding kinetics and mechanism of the B domain of Staphylococcal protein A (BdpA), an ultrafast-folding sequence. A combination of experimental biophysical analysis and atomistic molecular dynamics simulations performed on the prototype protein and several heterogeneous-backbone variants reveal the interplay among backbone flexibility, folding rates, and structural details of the transition state ensemble. Collectively, these findings suggest a significant degree of plasticity in the mechanisms that can give rise to ultrafast folding in the BdpA sequence and provide atomic level insights into how protein mimetic chains adopt an ordered folded state.

Thrombin - A Molecular Dynamics Perspective

Thrombin is a crucial enzyme involved in blood coagulation, essential for maintaining circulatory system integrity and preventing excessive bleeding. However, thrombin is also implicated in pathological conditions such as thrombosis and cancer. Despite the application of various experimental techniques, including X-ray crystallography, NMR spectroscopy, and HDXMS, none of these methods can precisely detect thrombin's dynamics and conformational ensembles at high spatial and temporal resolution. Fortunately, molecular dynamics (MD) simulation, a computational technique that allows the investigation of molecular functions and dynamics in atomic detail, can be used to explore thrombin behavior. This review summarizes recent MD simulation studies on thrombin and its interactions with other biomolecules. Specifically, the 17 studies discussed here provide insights into thrombin's switch between 'slow' and 'fast' forms, active and inactive forms, the role of Na+ binding, the effects of light chain mutation, and thrombin's interactions with other biomolecules. The findings of these studies have significant implications for developing new therapies for thrombosis and cancer. By understanding thrombin's complex behavior, researchers can design more effective drugs and treatments that target thrombin.

On committor functions in milestoning

As an optimal one-dimensional reaction coordinate, the committor function not only describes the probability of a trajectory initiated at a phase space point first reaching the product state before reaching the reactant state but also preserves the kinetics when utilized to run a reduced dynamics model. However, calculating the committor function in high-dimensional systems poses significant challenges. In this paper, within the framework of milestoning, exact expressions for committor functions at two levels of coarse graining are given, including committor functions of phase space point to point (CFPP) and milestone to milestone (CFMM). When combined with transition kernels obtained from trajectory analysis, these expressions can be utilized to accurately and efficiently compute the committor functions. Furthermore, based on the calculated committor functions, an adaptive algorithm is developed to gradually refine the transition state region. Finally, two model examples are employed to assess the accuracy of these different formulations of committor functions.

LPATH: A Semiautomated Python Tool for Clustering Molecular Pathways

The pathways by which a molecular process transitions to a target state are highly sought-after as direct views of a transition mechanism. While great strides have been made in the physics-based simulation of such pathways, the analysis of these pathways can be a major challenge due to their diversity and variable lengths. Here, we present the LPATH Python tool, which implements a semiautomated method for linguistics-assisted clustering of pathways into distinct classes (or routes). This method involves three steps: 1) discretizing the configurational space into key states, 2) extracting a text-string sequence of key visited states for each pathway, and 3) pairwise matching of pathways based on a text-string similarity score. To circumvent the prohibitive memory requirements of the first step, we have implemented a general two-stage method for clustering conformational states that exploits machine learning. LPATH is primarily designed for use with the WESTPA software for weighted ensemble simulations; however, the tool can also be applied to conventional simulations. As demonstrated for the C7eq to C7ax conformational transition of the alanine dipeptide, LPATH provides physically reasonable classes of pathways and corresponding probabilities.

Adaptive Sampling Methods for Molecular Dynamics in the Era of Machine Learning

Molecular dynamics (MD) simulations are fundamental computational tools for the study of proteins and their free energy landscapes. However, sampling protein conformational changes through MD simulations is challenging due to the relatively long time scales of these processes. Many enhanced sampling approaches have emerged to tackle this problem, including biased sampling and path-sampling methods. In this Perspective, we focus on adaptive sampling algorithms. These techniques differ from other approaches because the thermodynamic ensemble is preserved and the sampling is enhanced solely by restarting MD trajectories at particularly chosen seeds rather than introducing biasing forces. We begin our treatment with an overview of theoretically transparent methods, where we discuss principles and guidelines for adaptive sampling. Then, we present a brief summary of select methods that have been applied to realistic systems in the past. Finally, we discuss recent advances in adaptive sampling methodology powered by deep learning techniques, as well as their shortcomings.

From Byte to Bench to Bedside: Molecular Dynamics Simulations and Drug Discovery

Molecular dynamics (MD) simulations and computer-aided drug design (CADD) have advanced substantially over the past two decades, thanks to continuous computer hardware and software improvements. Given these advancements, MD simulations are poised to become even more powerful tools for investigating the dynamic interactions between potential small-molecule drugs and their target proteins, with significant implications for pharmacological research.

How Robust Is the Ligand Binding Transition State?

For many drug targets, it has been shown that the kinetics of drug binding (e.g., on rate and off rate) is more predictive of drug efficacy than thermodynamic quantities alone. This motivates the development of predictive computational models that can be used to optimize compounds on the basis of their kinetics. The structural details underpinning these computational models are found not only in the bound state but also in the short-lived ligand binding transition states. Although transition states cannot be directly observed experimentally due to their extremely short lifetimes, recent successes have demonstrated that modeling the ligand binding transition state is possible with the help of enhanced sampling molecular dynamics methods. Previously, we generated unbinding paths for an inhibitor of soluble epoxide hydrolase (sEH) with a residence time of 11 min. Here, we computationally modeled unbinding events with the weighted ensemble method REVO (resampling of ensembles by variation optimization) for five additional inhibitors of sEH with residence times ranging from 14.25 to 31.75 min, with average prediction accuracy within an order of magnitude. The unbinding ensembles are analyzed in detail, focusing on features of the ligand binding transition state ensembles (TSEs). We find that ligands with similar bound poses can show significant differences in their ligand binding TSEs, in terms of their spatial distribution and protein-ligand interactions. However, we also find similarities across the TSEs when examining more general features such as ligand degrees of freedom. Together these findings show significant challenges for rational, kinetics-based drug design.

Understanding the free-energy landscape of phase separation in lipid bilayers using molecular dynamics

Liquid-liquid phase separation inside the cell often results in biological condensates that can critically affect cell homeostasis. Such phase separation events occur in multiple parts of cells, including the cell membranes, where the "lipid raft" hypothesis posits the formation of ordered domains floating in a sea of disordered lipids. The resulting lipid domains often have functional roles. However, the thermodynamics of lipid phase separation and their resulting mechanistic effects on cell function and dysfunction are poorly understood. Understanding such complex phenomena in cell membranes, with their diverse lipid compositions, is exceptionally difficult. For these reasons, simple model systems that can recapitulate similar behavior are widely used to study this phenomenon. Despite these simplifications, the timescale and length scales of domain formation pose a challenge for molecular dynamics (MD) simulations. Thus, most MD studies focus on spontaneous lipid phase separation-essentially measuring the sign (but not the amplitude) of the free-energy change upon separation-rather than directly interrogating the thermodynamics. Here, we propose a proof-of-concept pipeline that can directly measure this free energy by combining coarse-grained MD with enhanced sampling protocols using a novel collective variable. This approach will be a useful tool to help connect the thermodynamics of phase separation with the mechanistic insights already available from MD simulations.

LPATH: A semi-automated Python tool for clustering molecular pathways

The pathways by which a molecular process transitions to a target state are highly sought-after as direct views of a transition mechanism. While great strides have been made in the physics-based simulation of such pathways, the analysis of these pathways can be a major challenge due to their diversity and variable lengths. Here we present the LPATH Python tool, which implements a semi-automated method for linguistics-assisted clustering of pathways into distinct classes (or routes). This method involves three steps: 1) discretizing the configurational space into key states, 2) extracting a text-string sequence of key visited states for each pathway, and 3) pairwise matching of pathways based on a text-string similarity score. To circumvent the prohibitive memory requirements of the first step, we have implemented a general two-stage method for clustering conformational states that exploits machine learning. LPATH is primarily designed for use with the WESTPA software for weighted ensemble simulations; however, the tool can also be applied to conventional simulations. As demonstrated for the C7eq to C7ax conformational transition of alanine dipeptide, LPATH provides physically reasonable classes of pathways and corresponding probabilities.

Approximating First Hitting Point Distribution in Milestoning for Rare Event Kinetics

Milestoning is an efficient method for rare event kinetics calculation using short trajectory parallelization. Mean first passage time (MFPT) is the key kinetic output of Milestoning, whose accuracy crucially depends on the initial distribution of the short trajectory ensemble. The true initial distribution, i.e., the first hitting point distribution (FHPD), has no analytic expression in the general case. Here, we introduce two algorithms, local passage time weighted Milestoning (LPT-M) and Bayesian inference Milestoning (BI-M), to accurately and efficiently approximate FHPD for systems at equilibrium condition. Starting from sampling the Boltzmann distribution on milestones, we calculate the proper weighting factor for the short trajectory ensemble. The methods are tested on two model examples for illustration purpose. Both methods improve significantly over the widely used classical Milestoning (CM) method in terms of the accuracy of MFPT. In particular, BI-M covers the directional Milestoning method as a special case in deterministic Hamiltonian dynamics. LPT-M is especially advantageous in terms of computational costs and robustness with respect to the increasing number of intermediate milestones. Furthermore, a locally iterative correction algorithm for nonequilibrium stationary FHPD is developed for exact MFPT calculation, which can be combined with LPT-M/BI-M and is much cheaper than the exact Milestoning (ExM) method.