Daniel M. Zuckerman

Equilibrium Sampling in Biomolecular Simulations

Equilibrium sampling of biomolecules remains an unmet challenge after more than 30 years of atomistic simulation. Efforts to enhance sampling capability, which are reviewed here, range from the development of new algorithms to parallelization to novel uses of hardware. Special focus is placed on classifying algorithms--most of which are underpinned by a few key ideas--in order to understand their fundamental strengths and limitations. Although algorithms have proliferated, progress resulting from novel hardware use appears to be more clear-cut than from algorithms alone, due partly to the lack of widely used sampling measures.

Efficient and verified simulation of a path ensemble for conformational change in a united-residue model of calmodulin

The computational sampling of rare, large-scale, conformational transitions in proteins is a well appreciated challenge—for which a number of potentially efficient path-sampling methodologies have been proposed. Here, we study a large-scale transition in a united-residue model of calmodulin using the “weighted ensemble” (WE) approach of Huber and Kim. Because of the models relative simplicity, we are able to compare our results with brute-force simulations. The comparison indicates that the WE approach quantitatively reproduces the brute-force results, as assessed by considering (i) the reaction rate, (ii) the distribution of event durations, and (iii) structural distributions describing the heterogeneity of the paths. Importantly, the WE method is readily applied to more chemically accurate models, and by studying a series of lower temperatures, our results suggest that the WE method can increase efficiency by orders of magnitude in more challenging systems.

Dynamic reaction paths and rates through importance-sampled stochastic dynamics

We extend a previously developed method, based on Wagner’s stochastic formulation of importance sampling, to the calculation of reaction rates and to a simple quantitative description of finite-temperature, average dynamic paths. Only the initial and final states are required as input—no information on transition state(s) is necessary. We demonstrate the method for a single particle moving on the two-dimensional Muller–Brown potential surface. Beyond computing the forward and reverse rates for this surface, we determine the average path, which exhibits “saddle point avoidance.” The method may be generalized to arbitrary numbers of degrees of freedom and to arbitrary types of stochastic dynamics.

Efficient use of nonequilibrium measurement to estimate free energy differences for molecular systems

A promising method for calculating free energy differences ΔF is to generate nonequilibrium data via “fast‐growth” simulations or by experiments—and then use Jarzynskis equality. However, a difficulty with using Jarzynskis equality is that ΔF estimates converge very slowly and unreliably due to the nonlinear nature of the calculation—thus requiring large, costly data sets. The purpose of the work presented here is to determine the best estimate for ΔF given a (finite) set of work values previously generated by simulation or experiment. Exploiting statistical properties of Jarzynskis equality, we present two fully automated analyses of nonequilibrium data from a toy model, and various simulated molecular systems. Both schemes remove at least several kBT of bias from ΔF estimates, compared to direct application of Jarzynskis equality, for modest sized data sets (100 work values), in all tested systems. Results from one of the new methods suggest that good estimates of ΔF can be obtained using 5–40‐fold less data than was previously possible. Extending previous work, the new results exploit the systematic behavior of bias due to finite sample size. A key innovation is better use of the more statistically reliable information available from the raw data.

Steady-state simulations using weighted ensemble path sampling

We extend the weighted ensemble (WE) path sampling method to perform rigorous statistical sampling for systems at steady state. A straightforward steady-state implementation of WE is directly practical for simple landscapes, but not when significant metastable intermediates states are present. We therefore develop an enhanced WE scheme, building on existing ideas, which accelerates attainment of steady state in complex systems. We apply both WE approaches to several model systems, confirming their correctness and efficiency by comparison with brute-force results. The enhanced version is significantly faster than the brute force and straightforward WE for systems with WE bins that accurately reflect the reaction coordinate(s). The new WE methods can also be applied to equilibrium sampling, since equilibrium is a steady state.

Simulations of the Alternating Access Mechanism of the Sodium Symporter Mhp1

Sodium coupled cotransporters of the five-helix inverted repeat (5HIR) superfamily use an alternating access mechanism to transport a myriad of small molecules across the cell membrane. One of the primary steps in this mechanism is the conformational transition from a state poised to bind extracellular substrates to a state that is competent to deliver substrate to the cytoplasm. Here, we construct a coarse-grained model of the 5HIR benzylhydantoin transporter Mhp1 that incorporates experimental structures of the outward- and inward-open states to investigate the mechanism of this conformational change. Using the weighted ensemble path-sampling method, we rigorously sample the outward- to inward-facing transition path ensemble. The transition path ensemble reveals a heterogeneous set of pathways connecting the two states and identifies two modes of transport: one consistent with a strict alternating access mechanism and another where decoupling of the inner and outer gates causes the transient formation of a continuous permeation pathway through the transporter. We also show that the conformational switch between the outward- and inward-open states results from rigid body motions of the hash motif relative to the substrate bundle, supporting the rocking bundle hypothesis. Finally, our methodology provides the groundwork for more chemically detailed investigations of the alternating mechanism.

Transition-event durations in one-dimensional activated processes

Despite their importance in activated processes, transition-event durations--which are much shorter than first passage times--have not received a complete theoretical treatment. The authors therefore study the distribution rhob(t) of durations of transition events over a barrier in a one-dimensional system undergoing overdamped Langevin dynamics. The authors show that rhob(t) is determined by a Fokker-Planck equation with absorbing boundary conditions and obtain a number of results, including (i) the analytic form of the asymptotic short-time transient behavior, which is universal and independent of the potential function; (ii) the first nonuniversal correction to the short-time behavior leading to an estimate of a key physical time scale; (iii) following previous work, a recursive formulation for calculating, exactly, all moments of rhob based solely on the potential function-along with approximations for the distribution based on a small number of moments; and (iv) a high-barrier approximation to the long-time (t-->infinity) behavior of rhob(t). The authors also find that the mean event duration does not depend simply on the barrier-top frequency (curvature) but is sensitive to details of the potential. All of the analytic results are confirmed by transition-path-sampling simulations implemented in a novel way. Finally, the authors discuss which aspects of the duration distribution are expected to be general for more complex systems.

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.

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.

Transition events in butane simulations: Similarities across models

From a variety of long simulations of all-atom butane using both stochastic and fully solved molecular dynamics, we have uncovered striking generic behavior that also occurs in one-dimensional systems. We find an apparently universal distribution of transition event durations, as well as a characteristic speed profile along the reaction coordinate. An approximate analytic distribution of event durations, derived from a one-dimensional model, correctly predicts the asymptotic behavior of the universal distribution for both short and long durations.