Weighted-ensemble Brownian dynamics simulations for protein association reactions

Improving double-ended transition state searches for soft-matter systems

Transitions between different stable configurations of biomolecules are important in understanding disease mechanisms, structure-function relations, and novel molecular-scale engineering. The corresponding pathways can be characterized efficiently using geometry optimization schemes based on double-ended transition state searches. An interpolation is first constructed between the known states and then refined, yielding a band that contains transition state candidates. Here, we analyze an example where various interpolation schemes lead to bands with a single step transition, but the correct pathway actually proceeds via an intervening, low-energy minimum. We compare a number of different interpolation schemes for this problem. We systematically alter the number of discrete images in the interpolations and the spring constants used in the optimization and test two schemes for adjusting the spring constants and image distribution, resulting in a total of 2760 different connection attempts. Our results confirm that optimized bands are not necessarily a good description of the transition pathways in themselves, and further refinement to actually converge transition states and establish their connectivity is required. We see an improvement in the optimized bands if we employ the adjustment of spring constants with doubly-nudged elastic band and a smaller improvement from the image redistribution. The example we consider is representative of numerous cases we have encountered in a wide variety of molecular and condensed matter systems.

Using Coarse-Grained Simulations to Characterize the Mechanisms of Protein-Protein Association

The formation of functionally versatile protein complexes underlies almost every biological process. The estimation of how fast these complexes can be formed has broad implications for unravelling the mechanism of biomolecular recognition. This kinetic property is traditionally quantified by association rates, which can be measured through various experimental techniques. To complement these time-consuming and labor-intensive approaches, we developed a coarse-grained simulation approach to study the physical processes of protein–protein association. We systematically calibrated our simulation method against a large-scale benchmark set. By combining a physics-based force field with a statistically-derived potential in the simulation, we found that the association rates of more than 80% of protein complexes can be correctly predicted within one order of magnitude relative to their experimental measurements. We further showed that a mixture of force fields derived from complementary sources was able to describe the process of protein–protein association with mechanistic details. For instance, we show that association of a protein complex contains multiple steps in which proteins continuously search their local binding orientations and form non-native-like intermediates through repeated dissociation and re-association. Moreover, with an ensemble of loosely bound encounter complexes observed around their native conformation, we suggest that the transition states of protein–protein association could be highly diverse on the structural level. Our study also supports the idea in which the association of a protein complex is driven by a “funnel-like” energy landscape. In summary, these results shed light on our understanding of how protein–protein recognition is kinetically modulated, and our coarse-grained simulation approach can serve as a useful addition to the existing experimental approaches that measure protein–protein association rates.

Frontier Expansion Sampling: A Method to Accelerate Conformational Search by Identifying Novel Seed Structures for Restart

Traditional molecular dynamics (MD) simulations have difficulties in tracking the slow molecular motions, at least partially due to the waste of sampling in already sampled regions. Here, we proposed a new enhanced sampling method, frontier expansion sampling (FEXS), to improve the sampling efficiency of molecular simulations by iteratively selecting seed structures diversely distributed at the "frontier" of an already sampled region to initiate new simulations. Different from other enhanced sampling methods, FEXS identifies the "frontier" seeds by integrating the Gaussian mixture model and the convex hull algorithm, which effectively improves the structural variation among the selected seeds and thus the descendant simulations. Validation in three protein systems, including the folding of chignolin, open-to-closed transition of maltodextrin binding protein, and internal conformational change of bovine pancreatic trypsin inhibitor, confirmed the effectiveness of this novel method in enhancing the sampling of conventional MD simulations to observe the large-scale protein conformational changes. When compared with other enhanced sampling methods like the structural dissimilarity sampling (SDS), FEXS reached at least the same level of sampling efficiency but was capable of providing complementary information in the three tested protein systems.

On Calculating Free Energy Differences Using Ensembles of Transition Paths

The free energy of a process is the fundamental quantity that determines its spontaneity or propensity at a given temperature. In particular, the binding free energy of a drug candidate to its biomolecular target is used as an objective quantity in drug design. Recently, binding kinetics—rates of association (k_on) and dissociation (k_off)—have also demonstrated utility for their ability to predict efficacy and in some cases have been shown to be more predictive than the binding free energy alone. Some methods exist to calculate binding kinetics from molecular simulations, although these are typically more difficult to calculate than the binding affinity as they depend on details of the transition path ensemble. Assessing these rate constants can be difficult, due to uncertainty in the definition of the bound and unbound states, large error bars and the lack of experimental data. As an additional consistency check, rate constants from simulation can be used to calculate free energies (using the log of their ratio) which can then be compared to free energies obtained experimentally or using alchemical free energy perturbation. However, in this calculation it is not straightforward to account for common, practical details such as the finite simulation volume or the particular definition of the “bound” and “unbound” states. Here we derive a set of correction terms that can be applied to calculations of binding free energies using full reactive trajectories. We apply these correction terms to revisit the calculation of binding free energies from rate constants for a host-guest system that was part of a blind prediction challenge, where significant deviations were observed between free energies calculated with rate ratios and those calculated from alchemical perturbation. The correction terms combine to significantly decrease the error with respect to computational benchmarks, from 3.4 to 0.76 kcal/mol. Although these terms were derived with weighted ensemble simulations in mind, some of the correction terms are generally applicable to free energies calculated using physical pathways via methods such as Markov state modeling, metadynamics, milestoning, or umbrella sampling.

OPTIMIZING WEIGHTED ENSEMBLE SAMPLING OF STEADY STATES

We propose parameter optimization techniques for weighted ensemble sampling of Markov chains in the steady-state regime. Weighted ensemble consists of replicas of a Markov chain, each carrying a weight, that are periodically resampled according to their weights inside of each of a number of bins that partition state space. We derive, from first principles, strategies for optimizing the choices of weighted ensemble parameters, in particular the choice of bins and the number of replicas to maintain in each bin. In a simple numerical example, we compare our new strategies with more traditional ones and with direct Monte Carlo.

ScMile: A Script to Investigate Kinetics with Short Time Molecular Dynamics Trajectories and the Milestoning Theory

Studies of complex and rare events in condensed phase systems continue to attract considerable attention. Milestoning is a useful theory and algorithm to investigate the long-time dynamics of activated molecular events. It is based on launching a large number of short trajectories and statistical analysis of the outcome. The implementation of the theory in a computer script is described that enables more efficient Milestoning calculation, reducing user time and errors, and automating a significant fraction of the algorithm. The script exploits a molecular dynamics engine, which at present is NAMD, to run the short trajectories. However, since the script is external to the engine, the script can be easily adapted to different molecular dynamics codes. The outcomes of the short trajectories are analyzed to obtain a kinetic and thermodynamic description of the entire process. While many examples of Milestoning were published in the past, we provide two simple examples (a conformational transition of alanine dipeptide in a vacuum and aqueous solution) to illustrate the use of the script.

Transition path sampling for non-equilibrium dynamics without predefined reaction coordinates

We develop two novel transition path sampling (TPS) algorithms for harvesting ensembles of rare event trajectories using non-equilibrium dynamics. These methods have the advantage that no predefined reaction coordinate is needed. Instead, an instantaneous reaction coordinate is based on the current path. Constituting a Monte Carlo random walk in trajectory space, the algorithms can be viewed as bridging between the original TPS methodology and the Rosenbluth based forward flux sampling methodology. We illustrate the new methods on toy models undergoing equilibrium and non-equilibrium dynamics, including an active Brownian particle system. For the latter, we find that transitions between steady states occur via states that are locally ordered but globally disordered.

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.

Extinction risk of a metapopulation under bistable local dynamics

We study the extinction risk of a fragmented population residing on a network of patches coupled by migration, where the local patch dynamics includes deterministic bistability. Mixing between patches is shown to dramatically influence the population's viability. We demonstrate that slow migration always increases the population's global extinction risk compared to the isolated case, while at fast migration synchrony between patches minimizes the population's extinction risk. Moreover, we discover a critical migration rate that maximizes the extinction risk of the population, and identify an early-warning signal when approaching this state. Our theoretical results are confirmed via the highly efficient weighted ensemble method. Notably, our theoretical formalism can also be applied to studying switching in gene regulatory networks with multiple transcriptional states.

Measurement of the real part of the Clausius-Mossotti factor of dielectrophoresis for Brownian particles

A method is proposed for measuring the real part of the Clausius-Mossotti factor ( K r ) of dielectrophoresis for Brownian particles based on a solution of the Smoluchowski equation using a designed polydimethysilloxane microchannel with planar hyperbolic electrodes on its glass substrate. An approximate two-dimensional spring-like dielectrophoretic force is generated in the device, and the data necessarily measured is the time evolution of the in-plane particle displacement undergoing confined Brownian motion. Validity of the measurement was checked against the zeta potentials in the literature based on the classical theory of surface conductance using polystyrene particles of size of one micron. As the dielectrophoretic force depends on K r , which is usually unknown for bio-particles and some engineered particles, and is seldom measured; this study is important from the academic point of view and could be helpful for the manipulation and characterization of sub-micron particles using dielectrophoresis. Extension of the method to the measurement of permanent dipole moment and total polarizability of particle was developed theoretically and discussed by incorporating an optical tweezer into the present device.

Transient probability currents provide upper and lower bounds on non-equilibrium steady-state currents in the Smoluchowski picture

Probability currents are fundamental in characterizing the kinetics of nonequilibrium processes. Notably, the steady-state current Jss for a source-sink system can provide the exact mean-first-passage time (MFPT) for the transition from the source to sink. Because transient nonequilibrium behavior is quantified in some modern path sampling approaches, such as the "weighted ensemble" strategy, there is strong motivation to determine bounds on Jss-and hence on the MFPT-as the system evolves in time. Here, we show that Jss is bounded from above and below by the maximum and minimum, respectively, of the current as a function of the spatial coordinate at any time t for one-dimensional systems undergoing overdamped Langevin (i.e., Smoluchowski) dynamics and for higher-dimensional Smoluchowski systems satisfying certain assumptions when projected onto a single dimension. These bounds become tighter with time, making them of potential practical utility in a scheme for estimating Jss and the long time scale kinetics of complex systems. Conceptually, the bounds result from the fact that extrema of the transient currents relax toward the steady-state current.

Large-Scale Conformational Changes and Protein Function: Breaking the in silico Barrier

Large-scale conformational changes are essential to link protein structures with their function at the cell and organism scale, but have been elusive both experimentally and computationally. Over the past few years developments in cryo-electron microscopy and crystallography techniques have started to reveal multiple snapshots of increasingly large and flexible systems, deemed impossible only short time ago. As structural information accumulates, theoretical methods become central to understand how different conformers interconvert to mediate biological function. Here we briefly survey current in silico methods to tackle large conformational changes, reviewing recent examples of cross-validation of experiments and computational predictions, which show how the integration of different scale simulations with biological information is already starting to break the barriers between the in silico, in vitro, and in vivo worlds, shedding new light onto complex biological problems inaccessible so far.

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.

Galerkin approximation of dynamical quantities using trajectory data

Understanding chemical mechanisms requires estimating dynamical statistics such as expected hitting times, reaction rates, and committors. Here, we present a general framework for calculating these dynamical quantities by approximating boundary value problems using dynamical operators with a Galerkin expansion. A specific choice of basis set in the expansion corresponds to the estimation of dynamical quantities using a Markov state model. More generally, the boundary conditions impose restrictions on the choice of basis sets. We demonstrate how an alternative basis can be constructed using ideas from diffusion maps. In our numerical experiments, this basis gives results of comparable or better accuracy to Markov state models. Additionally, we show that delay embedding can reduce the information lost when projecting the system's dynamics for model construction; this improves estimates of dynamical statistics considerably over the standard practice of increasing the lag time.

REVO: Resampling of ensembles by variation optimization

Conventional molecular dynamics simulations are incapable of sampling many important interactions in biomolecular systems due to their high dimensionality and rough energy landscapes. To observe rare events and calculate transition rates in these systems, enhanced sampling is a necessity. In particular, the study of ligand-protein interactions necessitates a diverse ensemble of protein conformations and transition states, and for many systems, this occurs on prohibitively long time scales. Previous strategies such as WExplore that can be used to determine these types of ensembles are hindered by problems related to the regioning of conformational space. Here, we propose a novel, regionless, enhanced sampling method that is based on the weighted ensemble framework. In this method, a value referred to as "trajectory variation" is optimized after each cycle through cloning and merging operations. This method allows for a more consistent measurement of observables and broader sampling resulting in the efficient exploration of previously unexplored conformations. We demonstrate the performance of this algorithm with the N-dimensional random walk and the unbinding of the trypsin-benzamidine system. The system is analyzed using conformation space networks, the residence time of benzamidine is confirmed, and a new unbinding pathway for the trypsin-benzamidine system is found. We expect that resampling of ensembles by variation optimization will be a useful general tool to broadly explore free energy landscapes.

Statistical Uncertainty Analysis for Small-Sample, High Log-Variance Data: Cautions for Bootstrapping and Bayesian Bootstrapping

Recent advances in molecular simulations allow the evaluation of previously unattainable observables, such as rate constants for protein folding. However, these calculations are usually computationally expensive, and even significant computing resources may result in a small number of independent estimates spread over many orders of magnitude. Such small-sample, high "log-variance" data are not readily amenable to analysis using the standard uncertainty (i.e., "standard error of the mean") because unphysical negative limits of confidence intervals result. Bootstrapping, a natural alternative guaranteed to yield a confidence interval within the minimum and maximum values, also exhibits a striking systematic bias of the lower confidence limit in log space. As we show, bootstrapping artifactually assigns high probability to improbably low mean values. A second alternative, the Bayesian bootstrap strategy, does not suffer from the same deficit and is more logically consistent with the type of confidence interval desired. The Bayesian bootstrap provides uncertainty intervals that are more reliable than those from the standard bootstrap method but must be used with caution nevertheless. Neither standard nor Bayesian bootstrapping can overcome the intrinsic challenge of underestimating the mean from small-size, high log-variance samples. Our conclusions are based on extensive analysis of model distributions and reanalysis of multiple independent atomistic simulations. Although we only analyze rate constants, similar considerations will apply to related calculations, potentially including highly nonlinear averages like the Jarzynski relation.

Practical rare event sampling for extreme mesoscale weather

Extreme mesoscale weather, including tropical cyclones, squall lines, and floods, can be enormously damaging and yet challenging to simulate; hence, there is a pressing need for more efficient simulation strategies. Here, we present a new rare event sampling algorithm called quantile diffusion Monte Carlo (quantile DMC). Quantile DMC is a simple-to-use algorithm that can sample extreme tail behavior for a wide class of processes. We demonstrate the advantages of quantile DMC compared to other sampling methods and discuss practical aspects of implementing quantile DMC. To test the feasibility of quantile DMC for extreme mesoscale weather, we sample extremely intense realizations of two historical tropical cyclones, 2010 Hurricane Earl and 2015 Hurricane Joaquin. Our results demonstrate quantile DMC's potential to provide low-variance extreme weather statistics while highlighting the work that is necessary for quantile DMC to attain greater efficiency in future applications.

Computational Estimation of Microsecond to Second Atomistic Folding Times

Despite the development of massively parallel computing hardware including inexpensive graphics processing units (GPUs), it has remained infeasible to simulate the folding of atomistic proteins at room temperature using conventional molecular dynamics (MD) beyond the microsecond scale. Here, we report the folding of atomistic, implicitly solvated protein systems with folding times τ ranging from ∼10 μs to ∼100 ms using the weighted ensemble (WE) strategy in combination with GPU computing. Starting from an initial structure or set of structures, WE organizes an ensemble of GPU-accelerated MD trajectory segments via intermittent pruning and replication events to generate statistically unbiased estimates of rate constants for rare events such as folding; no biasing forces are used. Although the variance among atomistic WE folding runs is significant, multiple independent runs are used to reduce and quantify statistical uncertainty. Folding times are estimated directly from WE probability flux and from history-augmented Markov analysis of the WE data. Three systems were examined: NTL9 at low solvent viscosity (yielding τf = 0.8-9 μs), NTL9 at water-like viscosity (τf = 0.2-2 ms), and Protein G at low viscosity (τf = 3-200 ms). In all cases, the folding time, uncertainty, and ensemble properties could be estimated from WE simulation; for Protein G, this characterization required significantly less overall computing than would be required to observe a single folding event with conventional MD simulations. Our results suggest that the use and calibration of force fields and solvent models for precise estimation of kinetic quantities is becoming feasible.

On the removal of initial state bias from simulation data

Classical atomistic simulations of biomolecules play an increasingly important role in molecular life science. The structure of current computing architectures favors methods that run multiple trajectories at once without requiring extensive communication between them. Many advanced sampling strategies in the field fit this mold. These approaches often rely on an adaptive logic and create ensembles of comparatively short trajectories whose starting points are not distributed according to the correct Boltzmann weights. This type of bias is notoriously difficult to remove, and Markov state models (MSMs) are one of the few strategies available for recovering the correct kinetics and thermodynamics from these ensembles of trajectories. In this contribution, we analyze the performance of MSMs in the thermodynamic reweighting task for a hierarchical set of systems. We show that MSMs can be rigorous tools to recover the correct equilibrium distribution for systems of sufficiently low dimensionality. This is conditional upon not tampering with local flux imbalances found in the data. For a real-world application, we find that a pure likelihood-based inference of the transition matrix produces the best results. The removal of the bias is incomplete, however, and for this system, all tested MSMs are outperformed by an alternative albeit less general approach rooted in the ideas of statistical resampling. We conclude by formulating some recommendations for how to address the reweighting issue in practice.

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.

On approximating a weak Markovian process as Markovian: Are we justified when discarding longtime correlations

The effect for removing weak longtime correlation is studied using a model system that contains a driven atom at liquid density under strong thermal fluctuations. The force that drives the tagged particle is about 1% of the average random force experienced by the particle. The tagged particle is allowed to assume a range of masses from 1/8 to 80 times that of a surrounding particle to study the effects of inertia. The driving force is indefinitely correlated but much weaker than "random" fluctuations from the environment. From this study, it is shown that the environmental influence is not fully random leading to the force autocorrelation function being a poor metric for detecting the correlated driving force. Although the velocity autocorrelation function shows stronger correlation for systems with higher inertia, the velocity autocorrelation function decays to a very small value of 2.5×10^-3 even for the most massive driven particle. For systems with small inertia, our study reveals that discarding longtime correlation has negligible influence on the first passage time (FPT) estimate, whereas for particles with large inertia, the deviation can indeed be appreciable. It is interesting that the Markov State Model (MSM) still produces reasonable estimates on the FPT even when a very short lag time that clearly violates the Markovianity assumption is used. This is likely a result of favorable error cancellations when the MSM transition probability matrices were constructed using trajectories that are much longer than the lag time.

A Multiscale Computational Model for Simulating the Kinetics of Protein Complex Assembly

Proteins fulfill versatile biological functions by interacting with each other and forming high-order complexes. Although the order in which protein subunits assemble is important for the biological function of their final complex, this kinetic information has received comparatively little attention in recent years. Here we describe a multiscale framework that can be used to simulate the kinetics of protein complex assembly. There are two levels of models in the framework. The structural details of a protein complex are reflected by the residue-based model, while a lower-resolution model uses a rigid-body (RB) representation to simulate the process of complex assembly. These two levels of models are integrated together, so that we are able to provide the kinetic information about complex assembly with both structural details and computational efficiency.

Escape of a Small Molecule from Inside T4 Lysozyme by Multiple Pathways

The T4 lysozyme L99A mutant is often used as a model system to study small-molecule binding to proteins, but pathways for ligand entry and exit from the buried binding site and the associated protein conformational changes have not been fully resolved. Here, molecular dynamics simulations were employed to model benzene exit from its binding cavity using the weighted ensemble (WE) approach to enhance sampling of low-probability unbinding trajectories. Independent WE simulations revealed four pathways for benzene exit, which correspond to transient tunnels spontaneously formed in previous simulations of apo T4 lysozyme. Thus, benzene unbinding occurs through multiple pathways partially created by intrinsic protein structural fluctuations. Motions of several α-helices and side chains were involved in ligand escape from metastable microstates. WE simulations also provided preliminary estimates of rate constants for each exit pathway. These results complement previous works and provide a semiquantitative characterization of pathway heterogeneity for binding of small molecules to proteins.

Dissociation Process of a MDM2/p53 Complex Investigated by Parallel Cascade Selection Molecular Dynamics and the Markov State Model

Recently, we efficiently generated dissociation pathways of a protein-ligand complex without applying force bias with parallel cascade selection molecular dynamics (PaCS-MD) and showed that PaCS-MD in combination with the Markov state model (MSM) yielded a binding free energy comparable to experimental values. In this work, we applied the same procedure to a complex of MDM2 protein and the transactivation domain of p53 protein (TAD-p53), the latter of which is known to be very flexible in the unbound state. Using 30 independent MD simulations in PaCS-MD, we successfully generated 25 dissociation pathways of the complex, which showed complete or partial unfolding of the helical region of TAD-p53 during the dissociation process within an average simulation time of 154.8 ± 46.4 ns. The standard binding free energy obtained in combination with one-dimensional-, three-dimensional (3D)- or C_α-MSM was in good agreement with those determined experimentally. Using 3D-MSM based on the center of mass position of TAD-p53 relative to MDM2, the dissociation rate constant was calculated, which was comparable to those measured experimentally. C_α-MSM based on all C_α coordinates of TAD-p53 reproduced the experimentally measured standard binding free energy, and dissociation and association rate constants. We conclude that the combination of PaCS-MD and MSM offers an efficient computational procedure to calculate binding free energies and kinetic rates.

OpenPathSampling: A Python Framework for Path Sampling Simulations. 1. Basics

Transition path sampling techniques allow molecular dynamics simulations of complex systems to focus on rare dynamical events, providing insight into mechanisms and the ability to calculate rates inaccessible by ordinary dynamics simulations. While path sampling algorithms are conceptually as simple as importance sampling Monte Carlo, the technical complexity of their implementation has kept these techniques out of reach of the broad community. Here, we introduce an easy-to-use Python framework called OpenPathSampling (OPS) that facilitates path sampling for (bio)molecular systems with minimal effort and yet is still extensible. Interfaces to OpenMM and an internal dynamics engine for simple models are provided in the initial release, but new molecular simulation packages can easily be added. Multiple ready-to-use transition path sampling methodologies are implemented, including standard transition path sampling (TPS) between reactant and product states and transition interface sampling (TIS) and its replica exchange variant (RETIS), as well as recent multistate and multiset extensions of transition interface sampling (MSTIS, MISTIS). In addition, tools are provided to facilitate the implementation of new path sampling schemes built on basic path sampling components. In this paper, we give an overview of the design of this framework and illustrate the simplicity of applying the available path sampling algorithms to a variety of benchmark problems.

Kinetics of Ligand-Protein Dissociation from All-Atom Simulations: Are We There Yet?

Large parallel gains in the development of both computational resources and sampling methods have now made it possible to simulate dissociation events in ligand-protein complexes with all-atom resolution. Such encouraging progress, together with the inherent spatiotemporal resolution associated with molecular simulations, has left their use for investigating dissociation processes brimming with potential, both in rational drug design, where it can be an invaluable tool for determining the mechanistic driving forces behind dissociation rate constants, and in force-field development, where it can provide a catalog of transient molecular structures with which to refine force fields. Although much progress has been made in making force fields more accurate, reducing their error for transient structures along a transition path could yet prove to be a critical development helping to make kinetic predictions much more accurate. In what follows, we will provide a state-of-the-art compilation of the enhanced sampling methods based on molecular dynamics (MD) simulations used to investigate the kinetics and mechanisms of ligand-protein dissociation processes. Due to the time scales of such processes being slower than what is accessible using straightforward MD simulations, several ingenious schemes are being devised at a rapid rate to overcome this obstacle. Here we provide an up-to-date compendium of such methods and their achievements and shortcomings in extracting mechanistic insight into ligand-protein dissociation. We conclude with a critical and provocative appraisal attempting to answer the title of this Perspective.

Trajectory stratification of stochastic dynamics

We present a general mathematical framework for trajectory stratification for simulating rare events. Trajectory stratification involves decomposing trajectories of the underlying process into fragments limited to restricted regions of state space (strata), computing averages over the distributions of the trajectory fragments within the strata with minimal communication between them, and combining those averages with appropriate weights to yield averages with respect to the original underlying process. Our framework reveals the full generality and flexibility of trajectory stratification, and it illuminates a common mathematical structure shared by existing algorithms for sampling rare events. We demonstrate the power of the framework by defining strata in terms of both points in time and path-dependent variables for efficiently estimating averages that were not previously tractable.

Best Practices for Quantification of Uncertainty and Sampling Quality in Molecular Simulations [Article v1.0]

The quantitative assessment of uncertainty and sampling quality is essential in molecular simulation. Many systems of interest are highly complex, often at the edge of current computational capabilities. Modelers must therefore analyze and communicate statistical uncertainties so that "consumers" of simulated data understand its significance and limitations. This article covers key analyses appropriate for trajectory data generated by conventional simulation methods such as molecular dynamics and (single Markov chain) Monte Carlo. It also provides guidance for analyzing some 'enhanced' sampling approaches. We do not discuss systematic errors arising, e.g., from inaccuracy in the chosen model or force field.

Conformational change of a biomolecule studied by the weighted ensemble method: Use of the diffusion map method to extract reaction coordinates

We simulate the nonequilibrium ensemble dynamics of a biomolecule using the weighted ensemble method, which was introduced in molecular dynamics simulations by Huber and Kim and further developed by Zuckerman and co-workers. As the order parameters to characterize its conformational change, we here use the coordinates derived from the diffusion map (DM) method, one of the manifold learning techniques. As a concrete example, we study the kinetic properties of a small peptide, chignolin in explicit water, and calculate the conformational change between the folded and misfolded states in a nonequilibrium way. We find that the transition time scales thus obtained are comparable to those using previously employed hydrogen-bond distances as the order parameters. Since the DM method only uses the 3D Cartesian coordinates of a peptide, this shows that the DM method can extract the important distance information of the peptide without relying on chemical intuition. The time scales are compared well with the previous results using different techniques, non-Markovian analysis and core-set milestoning for a single long trajectory. We also find that the most significant DM coordinate turns out to extract a dihedral angle of glycine, and the previously studied relaxation modes are well correlated with the most significant DM coordinates.

Exploring Configuration Space and Path Space of Biomolecules Using Enhanced Sampling Techniques-Searching for Mechanism and Kinetics of Biomolecular Functions

To understand functions of biomolecules such as proteins, not only structures but their conformational change and kinetics need to be characterized, but its atomistic details are hard to obtain both experimentally and computationally. Here, we review our recent computational studies using novel enhanced sampling techniques for conformational sampling of biomolecules and calculations of their kinetics. For efficiently characterizing the free energy landscape of a biomolecule, we introduce the multiscale enhanced sampling method, which uses a combined system of atomistic and coarse-grained models. Based on the idea of Hamiltonian replica exchange, we can recover the statistical properties of the atomistic model without any biases. We next introduce the string method as a path search method to calculate the minimum free energy pathways along a multidimensional curve in high dimensional space. Finally we introduce novel methods to calculate kinetics of biomolecules based on the ideas of path sampling: one is the Onsager⁻Machlup action method, and the other is the weighted ensemble method. Some applications of the above methods to biomolecular systems are also discussed and illustrated.

Mapping the Ligand Binding Landscape

The interaction between a ligand and a protein involves a multitude of conformational states. To achieve a particular deeply bound pose, the ligand must search across a rough free-energy landscape with many metastable minima. Creating maps of the ligand binding landscape is a great challenge, as binding and release events typically occur on timescales that are beyond the reach of molecular simulation. The WExplore enhanced sampling method is well suited to build these maps because it is designed to broadly explore free-energy landscapes and is capable of simulating ligand release pathways that occur on timescales as long as minutes. WExplore also uses only unbiased trajectory segments, allowing for the construction of Markov state models (MSMs) and conformation space networks that combine the results of multiple simulations. Here, we use WExplore to study two bromodomain-inhibitor systems using multiple docked starting poses (Brd4-MS436 and Baz2B-ICR7) and synthesize our results using a series of MSMs using time-lagged independent component analysis. Ranking the starting poses by exit rate agrees with the crystal structure pose in both cases. We also predict the most stable pose using the equilibrium populations from the MSM but find that the prediction is not robust as a function of MSM parameters. The simulated trajectories are synthesized into network models that visualize the entire binding landscape for each system, and we examine transition paths between deeply bound stable states. We find that, on average, transitions between deeply bound states convert through the unbound state 81% of the time, implying a trial-and-error approach to ligand binding. We conclude with a discussion of the implications of this result for both kinetics-based drug discovery and virtual screening pipelines that incorporate molecular dynamics.

Predicting ligand binding affinity using on- and off-rates for the SAMPL6 SAMPLing challenge

Interest in ligand binding kinetics has been growing rapidly, as it is being discovered in more and more systems that ligand residence time is the crucial factor governing drug efficacy. Many enhanced sampling methods have been developed with the goal of predicting ligand binding rates ([Formula: see text]) and/or ligand unbinding rates ([Formula: see text]) through explicit simulation of ligand binding pathways, and these methods work by very different mechanisms. Although there is not yet a blind challenge for ligand binding kinetics, here we take advantage of experimental measurements and rigorously computed benchmarks to compare estimates of [Formula: see text] calculated as the ratio of two rates: [Formula: see text]. These rates were determined using a new enhanced sampling method based on the weighted ensemble framework that we call "REVO": Reweighting of Ensembles by Variance Optimization. This is a further development of the WExplore enhanced sampling method, in which trajectory cloning and merging steps are guided not by the definition of sampling regions, but by maximizing trajectory variance. Here we obtain estimates of [Formula: see text] and [Formula: see text] that are consistent across multiple simulations, with an average log10-scale standard deviation of 0.28 for on-rates and 0.56 for off-rates, which is well within an order of magnitude and far better than previously observed for previous applications of the WExplore algorithm. Our rank ordering of the three host-guest pairs agrees with the reference calculations, however our predicted [Formula: see text] values were systematically lower than the reference by an average of 4.2 kcal/mol. Using tree network visualizations of the trajectories in the REVO algorithm, and conformation space networks for each system, we analyze the results of our sampling, and hypothesize sources of discrepancy between our [Formula: see text] values and the reference. We also motivate the direct inclusion of [Formula: see text] and [Formula: see text] challenges in future iterations of SAMPL, to further develop the field of ligand binding kinetics prediction and modeling.

Multiscale enhanced sampling of glucokinase: Regulation of the enzymatic reaction via a large scale domain motion

Enhanced sampling yields a comprehensive structural ensemble or a free energy landscape, which is beyond the capability of a conventional molecular dynamics simulation. Our recently developed multiscale enhanced sampling (MSES) method employs a coarse-grained model coupled with the target physical system for the efficient acceleration of the dynamics. MSES has demonstrated applicability to large protein systems in solution, such as intrinsically disordered proteins and protein-protein and protein-ligand interactions. Here, we applied the MSES simulation to an important drug discovery target, glucokinase (GCK), to elucidate the structural basis of the positive cooperativity of the enzymatic reaction at an atomistic resolution. MSES enabled us to compare two sets of the free energy landscapes of GCK, for the glucose-bound and glucose-unbound forms, and thus demonstrated the drastic change of the free energy surface depending on the glucose concentration. In the glucose-bound form, we found two distinct basins separated by a high energy barrier originating from the domain motion and the folding/unfolding of the α13 helix. By contrast, in the glucose-unbound form, a single flat basin extended to the open and super-open states. These features illustrated the two distinct phases achieving the cooperativity, the fast reaction cycle staying in the closed state at a high glucose concentration and the slow cycle primarily in the open/super-open state at a low concentration. The weighted ensemble simulations revealed the kinetics of the structural changes in GCK with the synergetic use of the MSES results; the rate constant of the transition between the closed state and the open/super-open states, k_C/O = 1.1 ms^-1, is on the same order as the experimental catalytic rate, k_cat = 0.22 ms^-1. Finally, we discuss the pharmacological activities of GCK activators (small molecular drugs modulating the GCK activity) in terms of the slight changes in the domain motion, depending on their chemical structures as regulators. The present study demonstrated the capability of the enhanced sampling and the associated kinetic calculations for understanding the atomistic structural dynamics of protein systems in physiological environments.

Rare-event sampling of epigenetic landscapes and phenotype transitions

Stochastic simulation has been a powerful tool for studying the dynamics of gene regulatory networks, particularly in terms of understanding how cell-phenotype stability and fate-transitions are impacted by noisy gene expression. However, gene networks often have dynamics characterized by multiple attractors. Stochastic simulation is often inefficient for such systems, because most of the simulation time is spent waiting for rare, barrier-crossing events to occur. We present a rare-event simulation-based method for computing epigenetic landscapes and phenotype-transitions in metastable gene networks. Our computational pipeline was inspired by studies of metastability and barrier-crossing in protein folding, and provides an automated means of computing and visualizing essential stationary and dynamic information that is generally inaccessible to conventional simulation. Applied to a network model of pluripotency in Embryonic Stem Cells, our simulations revealed rare phenotypes and approximately Markovian transitions among phenotype-states, occurring with a broad range of timescales. The relative probabilities of phenotypes and the transition paths linking pluripotency and differentiation are sensitive to global kinetic parameters governing transcription factor-DNA binding kinetics. Our approach significantly expands the capability of stochastic simulation to investigate gene regulatory network dynamics, which may help guide rational cell reprogramming strategies. Our approach is also generalizable to other types of molecular networks and stochastic dynamics frameworks.

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.