Simulating rare events using a weighted ensemble-based string method

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.

Finite temperature string by K-means clustering sampling with order parameters as collective variables for molecular crystals: application to polymorphic transformation between β-CL-20 and ε-CL-20

Polymorphic transformation of molecular crystals is a fundamental phase transition process, and it is important practically in the chemical, material, biopharmaceutical, and energy storage industries. However, understanding of the transformation mechanism at the molecular level is poor due to the extreme simulating challenges in enhanced sampling and formulating order parameters (OPs) as the collective variables that can distinguish polymorphs with quite similar and complicated structures so as to describe the reaction coordinate. In this work, two kinds of OPs for CL-20 were constructed by the bond distances, bond orientations and relative orientations. A K-means clustering algorithm based on the Euclidean distance and sample weight was used to smooth the initial finite temperature string (FTS), and the minimum free energy path connecting β-CL-20 and ε-CL-20 was sketched by the string method in collective variables, and the free energy profile along the path and the nucleation kinetics were obtained by Markovian milestoning with Voronoi tessellations. In comparison with the average-based sampling, the K-means clustering algorithm provided an improved convergence rate of FTS. The simulation of transformation was independent of OP types but was affected greatly by finite-size effects. A surface-mediated local nucleation mechanism was confirmed and the configuration located at the shoulder of potential of mean force, rather than overall maximum, was confirmed to be the critical nucleus formed by the cooperative effect of the intermolecular interactions. This work provides an effective way to explore the polymorphic transformation of caged molecular crystals at the molecular level.

Controlling Uncertainty of Empirical First-Passage Times in the Small-Sample Regime

We derive general bounds on the probability that the empirical first-passage time τ[over ¯]_{n}≡∑_{i=1}^{n}τ_{i}/n of a reversible ergodic Markov process inferred from a sample of n independent realizations deviates from the true mean first-passage time by more than any given amount in either direction. We construct nonasymptotic confidence intervals that hold in the elusive small-sample regime and thus fill the gap between asymptotic methods and the Bayesian approach that is known to be sensitive to prior belief and tends to underestimate uncertainty in the small-sample setting. We prove sharp bounds on extreme first-passage times that control uncertainty even in cases where the mean alone does not sufficiently characterize the statistics. Our concentration-of-measure-based results allow for model-free error control and reliable error estimation in kinetic inference, and are thus important for the analysis of experimental and simulation data in the presence of limited sampling.

Weighted ensemble: Recent mathematical developments

Weighted ensemble (WE) is an enhanced sampling method based on periodically replicating and pruning trajectories generated in parallel. WE has grown increasingly popular for computational biochemistry problems due, in part, to improved hardware and accessible software implementations. Algorithmic and analytical improvements have played an important role, and progress has accelerated in recent years. Here, we discuss and elaborate on the WE method from a mathematical perspective, highlighting recent results that enhance the computational efficiency. The mathematical theory reveals a new strategy for optimizing trajectory management that approaches the best possible variance while generalizing to systems of arbitrary dimension.

Identifying a Feasible Transition Pathway between Two Conformational States for a Protein

Proteins usually need to transit between different conformational states to fulfill their biological functions. In the mechanistic study of such transition processes by molecular dynamics simulations, identification of the minimum free energy path (MFEP) can substantially reduce the sampling space, thus enabling rigorous thermodynamic evaluation of the process. Conventionally, the MFEP is derived by iterative local optimization from an initial path, which is typically generated by simple brute force techniques like the targeted molecular dynamics (tMD). Therefore, the quality of the initial path determines the successfulness of MFEP estimation. In this work, we propose a method to improve derivation of the initial path. Through iterative relaxation-biasing simulations in a bidirectional manner, this method can construct a feasible transition pathway connecting two known states for a protein. Evaluation on small, fast-folding proteins against long equilibrium trajectories supports the good sampling efficiency of our method. When applied to larger proteins including the catalytic domain of human c-Src kinase as well as the converter domain of myosin VI, the paths generated by our method deviate significantly from those computed with the generic tMD approach. More importantly, free energy profiles and intermediate states obtained from our paths exhibit remarkable improvements over those from tMD paths with respect to both physical rationality and consistency with a priori knowledge.

A Companion Guide to the String Method with Swarms of Trajectories: Characterization, Performance, and Pitfalls

The string method with swarms of trajectories (SMwST) is an algorithm that identifies a physically meaningful transition pathway─a one-dimensional curve, embedded within a high-dimensional space of selected collective variables. The SMwST algorithm leans on a series of short, unbiased molecular dynamics simulations spawned at different locations of the discretized path, from whence an average dynamic drift is determined to evolve the string toward an optimal pathway. However conceptually simple in both its theoretical formulation and practical implementation, the SMwST algorithm is computationally intensive and requires a careful choice of parameters for optimal cost-effectiveness in applications to challenging problems in chemistry and biology. In this contribution, the SMwST algorithm is presented in a self-contained manner, discussing with a critical eye its theoretical underpinnings, applicability, inherent limitations, and use in the context of path-following free-energy calculations and their possible extension to kinetics modeling. Through multiple simulations of a prototypical polypeptide, combining the search of the transition pathway and the computation of the potential of mean force along it, several practical aspects of the methodology are examined with the objective of optimizing the computational effort, yet without sacrificing accuracy. In light of the results reported here, we propose some general guidelines aimed at improving the efficiency and reliability of the computed pathways and free-energy profiles underlying the conformational transitions at hand.

WESTPA 2.0: High-Performance Upgrades for Weighted Ensemble Simulations and Analysis of Longer-Timescale Applications

The weighted ensemble (WE) family of methods is one of several statistical mechanics-based path sampling strategies that can provide estimates of key observables (rate constants and pathways) using a fraction of the time required by direct simulation methods such as molecular dynamics or discrete-state stochastic algorithms. WE methods oversee numerous parallel trajectories using intermittent overhead operations at fixed time intervals, enabling facile interoperability with any dynamics engine. Here, we report on the major upgrades to the WESTPA software package, an open-source, high-performance framework that implements both basic and recently developed WE methods. These upgrades offer substantial improvements over traditional WE methods. The key features of the new WESTPA 2.0 software enhance the efficiency and ease of use: an adaptive binning scheme for more efficient surmounting of large free energy barriers, streamlined handling of large simulation data sets, exponentially improved analysis of kinetics, and developer-friendly tools for creating new WE methods, including a Python API and resampler module for implementing both binned and "binless" WE strategies.

A Minimal, Adaptive Binning Scheme for Weighted Ensemble Simulations

A promising approach for simulating rare events with rigorous kinetics is the weighted ensemble path sampling strategy. One challenge of this strategy is the division of configurational space into bins for sampling. Here we present a minimal adaptive binning (MAB) scheme for the automated, adaptive placement of bins along a progress coordinate within the framework of the weighted ensemble strategy. Results reveal that the MAB binning scheme, despite its simplicity, is more efficient than a manual, fixed binning scheme in generating transitions over large free energy barriers, generating a diversity of pathways, estimating rate constants, and sampling conformations. The scheme is general and extensible to any rare-events sampling strategy that employs progress coordinates.

Wepy: A Flexible Software Framework for Simulating Rare Events with Weighted Ensemble Resampling

Here, we introduce the open-source software framework wepy (https://github.com/ADicksonLab/wepy) which is a toolkit for running and analyzing weighted ensemble (WE) simulations. The wepy toolkit is in pure Python and as such is highly portable and extensible, making it an excellent platform to develop and use new WE resampling algorithms such as WExplore, REVO, and others while leveraging the entire Python ecosystem. In addition, wepy simplifies WE-specific analyses by defining out-of-core tree-like data structures using the cross-platform HDF5 file format. In this paper, we discuss the motivations and challenges for simulating rare events in biomolecular systems. As has previously been shown, high-dimensional WE resampling algorithms such as WExplore and REVO have been successful at these tasks, especially for rare events that are difficult to describe by one or two collective variables. We explain in detail how wepy facilitates implementation of these algorithms, as well as aids in analyzing the unique structure of WE simulation results. To explain how wepy and WE work in general, we describe the mathematical formalism of WE, an overview of the architecture of wepy, and provide code examples of how to construct, run, and analyze simulation results for a protein-ligand system (T4 Lysozyme in an implicit solvent). This paper is written with a variety of readers in mind, including (1) those curious about how to leverage WE rare-event simulations for their domain, (2) current WE users who want to begin using new high-dimensional resamplers such as WExplore and REVO, and (3) expert users who would like to prototype or implement their own algorithms that can be easily adopted by others.

Automatic error control during forward flux sampling of rare events in master equation models

Enhanced sampling methods, such as forward flux sampling (FFS), have great capacity for accelerating stochastic simulations of nonequilibrium biochemical systems involving rare events. However, the description of the tradeoffs between simulation efficiency and error in FFS remains incomplete. We present a novel and mathematically rigorous analysis of the errors in FFS that, for the first time, covers the contribution of every phase of the simulation. We derive a closed form expression for the optimally efficient count of samples to take in each FFS phase in terms of a fixed constraint on sampling error. We introduce a new method, forward flux pilot sampling (FFPilot), that is designed to take full advantage of our optimizing equation without prior information or assumptions about the phase weights and costs along the transition path. In simulations of both single and multidimensional gene regulatory networks, FFPilot is able to completely control sampling error. We then discuss how memory effects can introduce additional error when relaxation along the transition path is slow. This extra error can be traced to correlations between the FFS phases and can be controlled by monitoring the covariance between them. Finally, we show that, in sets of simulations with matched error, FFPilot is on the order of tens-to-hundreds of times faster than direct sampling and noticeably more efficient than previous FFS methods.

A Suite of Tutorials for the WESTPA Rare-Events Sampling Software [Article v1.0]

The weighted ensemble (WE) strategy has been demonstrated to be highly efficient in generating pathways and rate constants for rare events such as protein folding and protein binding using atomistic molecular dynamics simulations. Here we present five tutorials instructing users in the best practices for preparing, carrying out, and analyzing WE simulations for various applications using the WESTPA software. Users are expected to already have significant experience with running standard molecular dynamics simulations using the underlying dynamics engine of interest (e.g. Amber, Gromacs, OpenMM). The tutorials range from a molecular association process in explicit solvent to more complex processes such as host-guest association, peptide conformational sampling, and protein folding.

Protein-protein binding pathways and calculations of rate constants using fully-continuous, explicit-solvent simulations

A grand challenge in the field of biophysics has been the complete characterization of protein-protein binding processes at atomic resolution. This characterization requires the direct simulation of binding pathways starting from the initial, unbound state and proceeding through states that are too transient to be captured by experiment. Here, we applied the weighted ensemble path sampling strategy to orchestrate atomistic simulation of protein-protein binding pathways. Our simulation generated 203 fully-continuous and independent pathways along with rate constants for the binding process involving the barnase and barstar proteins. Results reveal multiple binding pathways along a "funnel-like" free energy landscape in which the formation of the "encounter complex" intermediate is rate-limiting followed by a relatively rapid rearrangement of the encounter complex to the bound state. Among all diffusional collisions, only ∼11% were productive. In the most probable binding pathways, the proteins rotated to a large extent (likely via electrostatic steering) in order to collide productively followed by "rolling" of the proteins along each other's binding interfaces to reach the bound state. Consistent with experiment, R59 was identified as the most kinetically important barnase residue for the binding process. Furthermore, protein desolvation occurs late in the binding process during the rearrangement of the encounter complex to the bound state. Notably, the positions of crystallographic water molecules that bridge hydrogen bonds between barnase and barstar are occupied in the bound-state ensemble. Our simulation was completed in a month using 1600 CPU cores at a time, demonstrating that it is now practical to carry out atomistic simulations of protein-protein binding.

Efficiently sampling conformations and pathways using the concurrent adaptive sampling (CAS) algorithm

Molecular dynamics simulations are useful in obtaining thermodynamic and kinetic properties of bio-molecules, but they are limited by the time scale barrier. That is, we may not obtain properties' efficiently because we need to run microseconds or longer simulations using femtosecond time steps. To overcome this time scale barrier, we can use the weighted ensemble (WE) method, a powerful enhanced sampling method that efficiently samples thermodynamic and kinetic properties. However, the WE method requires an appropriate partitioning of phase space into discrete macrostates, which can be problematic when we have a high-dimensional collective space or when little is known a priori about the molecular system. Hence, we developed a new WE-based method, called the "Concurrent Adaptive Sampling (CAS) algorithm," to tackle these issues. The CAS algorithm is not constrained to use only one or two collective variables, unlike most reaction coordinate-dependent methods. Instead, it can use a large number of collective variables and adaptive macrostates to enhance the sampling in the high-dimensional space. This is especially useful for systems in which we do not know what the right reaction coordinates are, in which case we can use many collective variables to sample conformations and pathways. In addition, a clustering technique based on the committor function is used to accelerate sampling the slowest process in the molecular system. In this paper, we introduce the new method and show results from two-dimensional models and bio-molecules, specifically penta-alanine and a triazine trimer.

Large enhancement of response times of a protein conformational switch by computational design

The design of protein conformational switches—or proteins that change conformations in response to a signal such as ligand binding—has great potential for developing novel biosensors, diagnostic tools, and therapeutic agents. Among the defining properties of such switches, the response time has been the most challenging to optimize. Here we apply a computational design strategy in synergistic combination with biophysical experiments to rationally improve the response time of an engineered protein-based Ca²⁺-sensor in which the switching process occurs via mutually exclusive folding of two alternate frames. Notably, our strategy identifies mutations that increase switching rates by as much as 32-fold, achieving response times on the order of fast physiological Ca²⁺ fluctuations. Our computational design strategy is general and may aid in optimizing the kinetics of other protein conformational switches.

The rational optimization of response times of protein conformational switches is a major challenge for biomolecular switch design. Here the authors present a generally applicable computational design strategy that in combination with biophysical experiments can improve response times using a Ca²⁺-sensor as an example.

Computational membrane biophysics: From ion channel interactions with drugs to cellular function

The rapid development of experimental and computational techniques has changed fundamentally our understanding of cellular-membrane transport. The advent of powerful computers and refined force-fields for proteins, ions, and lipids has expanded the applicability of Molecular Dynamics (MD) simulations. A myriad of cellular responses is modulated through the binding of endogenous and exogenous ligands (e.g. neurotransmitters and drugs, respectively) to ion channels. Deciphering the thermodynamics and kinetics of the ligand binding processes to these membrane proteins is at the heart of modern drug development. The ever-increasing computational power has already provided insightful data on the thermodynamics and kinetics of drug-target interactions, free energies of solvation, and partitioning into lipid bilayers for drugs. This review aims to provide a brief summary about modeling approaches to map out crucial binding pathways with intermediate conformations and free-energy surfaces for drug-ion channel binding mechanisms that are responsible for multiple effects on cellular functions. We will discuss post-processing analysis of simulation-generated data, which are then transformed to kinetic models to better understand the molecular underpinning of the experimental observables under the influence of drugs or mutations in ion channels. This review highlights crucial mathematical frameworks and perspectives on bridging different well-established computational techniques to connect the dynamics and timescales from all-atom MD and free energy simulations of ion channels to the physiology of action potentials in cellular models. This article is part of a Special Issue entitled: Biophysics in Canada, edited by Lewis Kay, John Baenziger, Albert Berghuis and Peter Tieleman.

Ab initio sampling of transition paths by conditioned Langevin dynamics

We propose a novel stochastic method to generate Brownian paths conditioned to start at an initial point and end at a given final point during a fixed time t_f under a given potential U(x). These paths are sampled with a probability given by the overdamped Langevin dynamics. We show that these paths can be exactly generated by a local stochastic partial differential equation. This equation cannot be solved in general but we present several approximations that are valid either in the low temperature regime or in the presence of barrier crossing. We show that this method warrants the generation of statistically independent transition paths. It is computationally very efficient. We illustrate the method first on two simple potentials, the two-dimensional Mueller potential and the Mexican hat potential, and then on the multi-dimensional problem of conformational transitions in proteins using the "Mixed Elastic Network Model" as a benchmark.

Weighted Ensemble Simulation: Review of Methodology, Applications, and Software

The weighted ensemble (WE) methodology orchestrates quasi-independent parallel simulations run with intermittent communication that can enhance sampling of rare events such as protein conformational changes, folding, and binding. The WE strategy can achieve superlinear scaling-the unbiased estimation of key observables such as rate constants and equilibrium state populations to greater precision than would be possible with ordinary parallel simulation. WE software can be used to control any dynamics engine, such as standard molecular dynamics and cell-modeling packages. This article reviews the theoretical basis of WE and goes on to describe successful applications to a number of complex biological processes-protein conformational transitions, (un)binding, and assembly processes, as well as cell-scale processes in systems biology. We furthermore discuss the challenges that need to be overcome in the next phase of WE methodological development. Overall, the combined advances in WE methodology and software have enabled the simulation of long-timescale processes that would otherwise not be practical on typical computing resources using standard simulation.

Markov State Models of gene regulatory networks

Background

Gene regulatory networks with dynamics characterized by multiple stable states underlie cell fate-decisions. Quantitative models that can link molecular-level knowledge of gene regulation to a global understanding of network dynamics have the potential to guide cell-reprogramming strategies. Networks are often modeled by the stochastic Chemical Master Equation, but methods for systematic identification of key properties of the global dynamics are currently lacking.

Results

The method identifies the number, phenotypes, and lifetimes of long-lived states for a set of common gene regulatory network models. Application of transition path theory to the constructed Markov State Model decomposes global dynamics into a set of dominant transition paths and associated relative probabilities for stochastic state-switching.

Conclusions

In this proof-of-concept study, we found that the Markov State Model provides a general framework for analyzing and visualizing stochastic multistability and state-transitions in gene networks. Our results suggest that this framework—adopted from the field of atomistic Molecular Dynamics—can be a useful tool for quantitative Systems Biology at the network scale.

Electronic supplementary material

The online version of this article (doi:10.1186/s12918-017-0394-4) contains supplementary material, which is available to authorized users.

Path-sampling strategies for simulating rare events in biomolecular systems

Despite more than three decades of effort with molecular dynamics simulations, long-timescale (ms and beyond) biologically relevant phenomena remain out of reach in most systems of interest. This is largely because important transitions, such as conformational changes and (un)binding events, tend to be rare for conventional simulations (<10μs). That is, conventional simulations will predominantly dwell in metastable states instead of making large transitions in complex biomolecular energy landscapes. In contrast, path sampling approaches focus computing effort specifically on transitions of interest. Such approaches have been in use for nearly 20 years in biomolecular systems and enabled the generation of pathways and calculation of rate constants for ms processes, including large protein conformational changes, protein folding, and protein (un)binding.

Simulating Current-Voltage Relationships for a Narrow Ion Channel Using the Weighted Ensemble Method

Ion channels are responsible for a myriad of fundamental biological processes via their role in controlling the flow of ions through water-filled membrane-spanning pores in response to environmental cues. Molecular simulation has played an important role in elucidating the mechanism of ion conduction, but connecting atomistically detailed structural models of the protein to electrophysiological measurements remains a broad challenge due to the computational cost of reaching the necessary time scales. Here, we introduce an enhanced sampling method for simulating the conduction properties of narrow ion channels using the Weighted ensemble (WE) sampling approach. We demonstrate the application of this method to calculate the current–voltage relationship as well as the nonequilibrium ion distribution at steady-state of a simple model ion channel. By direct comparisons with long brute force simulations, we show that the WE simulations rigorously reproduce the correct long-time scale kinetics of the system and are capable of determining these quantities using significantly less aggregate simulation time under conditions where permeation events are rare.

Efficient Atomistic Simulation of Pathways and Calculation of Rate Constants for a Protein-Peptide Binding Process: Application to the MDM2 Protein and an Intrinsically Disordered p53 Peptide

The characterization of protein binding processes - with all of the key conformational changes - has been a grand challenge in the field of biophysics. Here, we have used the weighted ensemble path sampling strategy to orchestrate molecular dynamics simulations, yielding atomistic views of protein-peptide binding pathways involving the MDM2 oncoprotein and an intrinsically disordered p53 peptide. A total of 182 independent, continuous binding pathways were generated, yielding a kon that is in good agreement with experiment. These pathways were generated in 15 days using 3500 cores of a supercomputer, substantially faster than would be possible with "brute force" simulations. Many of these pathways involve the anchoring of p53 residue F19 into the MDM2 binding cleft when forming the metastable encounter complex, indicating that F19 may be a kinetically important residue. Our study demonstrates that it is now practical to generate pathways and calculate rate constants for protein binding processes using atomistic simulation on typical computing resources.

Principles and Overview of Sampling Methods for Modeling Macromolecular Structure and Dynamics

Investigation of macromolecular structure and dynamics is fundamental to understanding how macromolecules carry out their functions in the cell. Significant advances have been made toward this end in silico, with a growing number of computational methods proposed yearly to study and simulate various aspects of macromolecular structure and dynamics. This review aims to provide an overview of recent advances, focusing primarily on methods proposed for exploring the structure space of macromolecules in isolation and in assemblies for the purpose of characterizing equilibrium structure and dynamics. In addition to surveying recent applications that showcase current capabilities of computational methods, this review highlights state-of-the-art algorithmic techniques proposed to overcome challenges posed in silico by the disparate spatial and time scales accessed by dynamic macromolecules. This review is not meant to be exhaustive, as such an endeavor is impossible, but rather aims to balance breadth and depth of strategies for modeling macromolecular structure and dynamics for a broad audience of novices and experts.

Unbiased Rare Event Sampling in Spatial Stochastic Systems Biology Models Using a Weighted Ensemble of Trajectories

The long-term goal of connecting scales in biological simulation can be facilitated by scale-agnostic methods. We demonstrate that the weighted ensemble (WE) strategy, initially developed for molecular simulations, applies effectively to spatially resolved cell-scale simulations. The WE approach runs an ensemble of parallel trajectories with assigned weights and uses a statistical resampling strategy of replicating and pruning trajectories to focus computational effort on difficult-to-sample regions. The method can also generate unbiased estimates of non-equilibrium and equilibrium observables, sometimes with significantly less aggregate computing time than would be possible using standard parallelization. Here, we use WE to orchestrate particle-based kinetic Monte Carlo simulations, which include spatial geometry (e.g., of organelles, plasma membrane) and biochemical interactions among mobile molecular species. We study a series of models exhibiting spatial, temporal and biochemical complexity and show that although WE has important limitations, it can achieve performance significantly exceeding standard parallel simulation—by orders of magnitude for some observables.

Stochastic simulations (simulations where randomness plays a role) of even simple biological systems are often so computationally intensive that it is impossible, in practice, to simulate them exhaustively and gather good statistics about the likelihood of different outcomes. The difficulty is compounded for the observation of rare events in these simulations; unfortunately, rare events, such as state transitions and barrier crossings, are often those of particular interest. Using the weighted ensemble (WE) method, we are able to enhance the characterization of rare events in cell biology simulations, but in such a way that the statistics for these events remain unbiased. The histogram of outcomes that WE produces has the same shape as a naive one, but the resolution of events in the tails of the histogram is greatly improved. This improved resolution in rare event statistics can be used to infer unbiased estimates of long timescale dynamics from short simulations, and we show that using a weighted ensemble can result in a reduction in total simulation time needed to sample certain events of interest in spatial, stochastic models of biological systems.

Estimating first-passage time distributions from weighted ensemble simulations and non-Markovian analyses

First-passage times (FPTs) are widely used to characterize stochastic processes such as chemical reactions, protein folding, diffusion processes or triggering a stock option. In previous work (Suarez et al., JCTC 2014;10:2658-2667), we demonstrated a non-Markovian analysis approach that, with a sufficient subset of history information, yields unbiased mean first-passage times from weighted-ensemble (WE) simulations. The estimation of the distribution of the first-passage times is, however, a more ambitious goal since it cannot be obtained by direct observation in WE trajectories. Likewise, a large number of events would be required to make a good estimation of the distribution from a regular "brute force" simulation. Here, we show how the previously developed non-Markovian analysis can generate approximate, but highly accurate, FPT distributions from WE data. The analysis can also be applied to any other unbiased trajectories, such as from standard molecular dynamics simulations. The present study employs a range of systems with independent verification of the distributions to demonstrate the success and limitations of the approach. By comparison to a standard Markov analysis, the non-Markovian approach is less sensitive to the user-defined discretization of configuration space.

Conformational transition paths harbor structures useful for aiding drug discovery and understanding enzymatic mechanisms in protein kinases

This short article examines the usefulness of fast simulations of conformational transition paths in elucidating enzymatic mechanisms and guiding drug discovery for protein kinases. It applies the transition path method in the MOIL software package to simulate the paths of conformational transitions between six pairs of structures from the Protein Data Bank. The structures along the transition paths were found to resemble experimental structures that mimic transient structures believed to form during enzymatic catalysis or conformational transitions, or structures that have drug candidates bound. These findings suggest that such simulations could provide quick initial insights into the enzymatic mechanisms or pathways of conformational transitions of proteins kinases, or could provide structures useful for aiding structure-based drug design.

Highly Efficient Computation of the Basal kon using Direct Simulation of Protein-Protein Association with Flexible Molecular Models

An essential baseline for determining the extent to which electrostatic interactions enhance the kinetics of protein–protein association is the “basal” k_on, which is the rate constant for association in the absence of electrostatic interactions. However, since such association events are beyond the milliseconds time scale, it has not been practical to compute the basal k_on by directly simulating the association with flexible models. Here, we computed the basal k_on for barnase and barstar, two of the most rapidly associating proteins, using highly efficient, flexible molecular simulations. These simulations involved (a) pseudoatomic protein models that reproduce the molecular shapes, electrostatic, and diffusion properties of all-atom models, and (b) application of the weighted ensemble path sampling strategy, which enhanced the efficiency of generating association events by >130-fold. We also examined the extent to which the computed basal k_on is affected by inclusion of intermolecular hydrodynamic interactions in the simulations.

DNA-Binding Kinetics Determines the Mechanism of Noise-Induced Switching in Gene Networks

Gene regulatory networks are multistable dynamical systems in which attractor states represent cell phenotypes. Spontaneous, noise-induced transitions between these states are thought to underlie critical cellular processes, including cell developmental fate decisions, phenotypic plasticity in fluctuating environments, and carcinogenesis. As such, there is increasing interest in the development of theoretical and computational approaches that can shed light on the dynamics of these stochastic state transitions in multistable gene networks. We applied a numerical rare-event sampling algorithm to study transition paths of spontaneous noise-induced switching for a ubiquitous gene regulatory network motif, the bistable toggle switch, in which two mutually repressive genes compete for dominant expression. We find that the method can efficiently uncover detailed switching mechanisms that involve fluctuations both in occupancies of DNA regulatory sites and copy numbers of protein products. In addition, we show that the rate parameters governing binding and unbinding of regulatory proteins to DNA strongly influence the switching mechanism. In a regime of slow DNA-binding/unbinding kinetics, spontaneous switching occurs relatively frequently and is driven primarily by fluctuations in DNA-site occupancies. In contrast, in a regime of fast DNA-binding/unbinding kinetics, switching occurs rarely and is driven by fluctuations in levels of expressed protein. Our results demonstrate how spontaneous cell phenotype transitions involve collective behavior of both regulatory proteins and DNA. Computational approaches capable of simulating dynamics over many system variables are thus well suited to exploring dynamic mechanisms in gene networks.

WESTPA: an interoperable, highly scalable software package for weighted ensemble simulation and analysis

The weighted ensemble (WE) path sampling approach orchestrates an ensemble of parallel calculations with intermittent communication to enhance the sampling of rare events, such as molecular associations or conformational changes in proteins or peptides. Trajectories are replicated and pruned in a way that focuses computational effort on underexplored regions of configuration space while maintaining rigorous kinetics. To enable the simulation of rare events at any scale (e.g., atomistic, cellular), we have developed an open-source, interoperable, and highly scalable software package for the execution and analysis of WE simulations: WESTPA (The Weighted Ensemble Simulation Toolkit with Parallelization and Analysis). WESTPA scales to thousands of CPU cores and includes a suite of analysis tools that have been implemented in a massively parallel fashion. The software has been designed to interface conveniently with any dynamics engine and has already been used with a variety of molecular dynamics (e.g., GROMACS, NAMD, OpenMM, AMBER) and cell-modeling packages (e.g., BioNetGen, MCell). WESTPA has been in production use for over a year, and its utility has been demonstrated for a broad set of problems, ranging from atomically detailed host–guest associations to nonspatial chemical kinetics of cellular signaling networks. The following describes the design and features of WESTPA, including the facilities it provides for running WE simulations and storing and analyzing WE simulation data, as well as examples of input and output.

Investigations of α-helix↔β-sheet transition pathways in a miniprotein using the finite-temperature string method

A parallel implementation of the finite-temperature string method is described, which takes into account the invariance of coordinates with respect to rigid-body motions. The method is applied to the complex α-helix↔β-sheet transition in a β-hairpin miniprotein in implicit solvent, which exhibits much of the complexity of conformational changes in proteins. Two transition paths are considered, one derived from a linear interpolant between the endpoint structures and the other derived from a targeted dynamics simulation. Two methods for computing the conformational free energy (FE) along the string are compared, a restrained method, and a tessellation method introduced by E. Vanden-Eijnden and M. Venturoli [J. Chem. Phys. 130, 194103 (2009)]. It is found that obtaining meaningful free energy profiles using the present atom-based coordinates requires restricting sampling to a vicinity of the converged path, where the hyperplanar approximation to the isocommittor surface is sufficiently accurate. This sampling restriction can be easily achieved using restraints or constraints. The endpoint FE differences computed from the FE profiles are validated by comparison with previous calculations using a path-independent confinement method. The FE profiles are decomposed into the enthalpic and entropic contributions, and it is shown that the entropy difference contribution can be as large as 10 kcal/mol for intermediate regions along the path, compared to 15-20 kcal/mol for the enthalpy contribution. This result demonstrates that enthalpic barriers for transitions are offset by entropic contributions arising from the existence of different paths across a barrier. The possibility of using systematically coarse-grained representations of amino acids, in the spirit of multiple interaction site residue models, is proposed as a means to avoid ad hoc sampling restrictions to narrow transition tubes.

Simultaneous Computation of Dynamical and Equilibrium Information Using a Weighted Ensemble of Trajectories

Equilibrium formally can be represented as an ensemble of uncoupled systems undergoing unbiased dynamics in which detailed balance is maintained. Many nonequilibrium processes can be described by suitable subsets of the equilibrium ensemble. Here, we employ the “weighted ensemble” (WE) simulation protocol [Huber and Kim, Biophys. J.1996, 70, 97–110] to generate equilibrium trajectory ensembles and extract nonequilibrium subsets for computing kinetic quantities. States do not need to be chosen in advance. The procedure formally allows estimation of kinetic rates between arbitrary states chosen after the simulation, along with their equilibrium populations. We also describe a related history-dependent matrix procedure for estimating equilibrium and nonequilibrium observables when phase space has been divided into arbitrary non-Markovian regions, whether in WE or ordinary simulation. In this proof-of-principle study, these methods are successfully applied and validated on two molecular systems: explicitly solvated methane association and the implicitly solvated Ala4 peptide. We comment on challenges remaining in WE calculations.

WExplore: hierarchical exploration of high-dimensional spaces using the weighted ensemble algorithm

As most relevant motions in biomolecular systems are inaccessible to conventional molecular dynamics simulations, algorithms that enhance sampling of rare events are indispensable. Increasing interest in intrinsically disordered systems and the desire to target ensembles of protein conformations (rather than single structures) in drug development motivate the need for enhanced sampling algorithms that are not limited to "two-basin" problems, and can efficiently determine structural ensembles. For systems that are not well-studied, this must often be done with little or no information about the dynamics of interest. Here we present a novel strategy to determine structural ensembles that uses dynamically defined sampling regions that are organized in a hierarchical framework. It is based on the weighted ensemble algorithm, where an ensemble of copies of the system ("replicas") is directed to new regions of configuration space through merging and cloning operations. The sampling hierarchy allows for a large number of regions to be defined, while using only a small number of replicas that can be balanced over multiple length scales. We demonstrate this algorithm on two model systems that are analytically solvable and examine the 10-residue peptide chignolin in explicit solvent. The latter system is analyzed using a configuration space network, and novel hydrogen bonds are found that facilitate folding.

Efficient stochastic simulation of chemical kinetics networks using a weighted ensemble of trajectories

We apply the "weighted ensemble" (WE) simulation strategy, previously employed in the context of molecular dynamics simulations, to a series of systems-biology models that range in complexity from a one-dimensional system to a system with 354 species and 3680 reactions. WE is relatively easy to implement, does not require extensive hand-tuning of parameters, does not depend on the details of the simulation algorithm, and can facilitate the simulation of extremely rare events. For the coupled stochastic reaction systems we study, WE is able to produce accurate and efficient approximations of the joint probability distribution for all chemical species for all time t. WE is also able to efficiently extract mean first passage times for the systems, via the construction of a steady-state condition with feedback. In all cases studied here, WE results agree with independent "brute-force" calculations, but significantly enhance the precision with which rare or slow processes can be characterized. Speedups over "brute-force" in sampling rare events via the Gillespie direct Stochastic Simulation Algorithm range from ~10(12) to ~10(18) for characterizing rare states in a distribution, and ~10(2) to ~10(4) for finding mean first passage times.