Reducing Weighted Ensemble Variance With Optimal Trajectory Management

TL;DR

This paper tackles the challenge of high-variance MFPT estimates in weighted ensemble simulations of biophysical processes. It extends a previous local MFPT–driven binning framework by leveraging history-augmented MSMs (haMSMs) to optimally partition phase space and allocate sampling toward high-variance regions, validated on synthetic Trp-cage dynamics and all-atom NTL9 folding in different frictions. Across synMD and high-friction atomistic models, the optimized binning consistently reduces run-to-run MFPT variance and improves the reliability of kinetic estimates, with the most dramatic gains in the slow-relaxation, high-friction regime. The approach offers a practical pathway to apply WE with principled parameterization to complex biomolecular kinetics, supported by open-source tools and a scalable workflow for bin construction and haMSM analysis.

Abstract

Weighted ensemble (WE) is an enhanced path-sampling method that is conceptually simple, widely applicable, and statistically exact. In a WE simulation, an ensemble of trajectories is periodically pruned or replicated to enhance sampling of rare transitions and improve estimation of mean first passage times (MFPTs). However, poor choices of the parameters governing pruning and replication can lead to high-variance MFPT estimates. Our previous work [J. Chem. Phys. 158, 014108 (2023)] presented an optimal WE parameterization strategy and applied it in low-dimensional example systems. The strategy harnesses estimated local MFPTs from different initial configurations to a single target state. In the present work, we apply the optimal parameterization strategy to more challenging, high-dimensional molecular models, namely, synthetic molecular dynamics (MD) models of Trp-cage folding and unfolding, as well as atomistic MD models of NTL9 folding in high-friction and low-friction continuum solvents. In each system we use WE to estimate the MFPT for folding or unfolding events. We show that the optimal parameterization reduces the variance of MFPT estimates in three of four systems, with dramatic improvement in the most challenging atomistic system. Overall, the parameterization strategy improves the accuracy and reliability of WE estimates for the kinetics of biophysical processes.

Figures (8)

• Figure 1: Rate constant estimation using WE (schematic). (A) The rate constant is estimated by the "flux" or trajectory weight entering the target state per unit time (red arrows). (B) Idealized flux curve over time for rate constant estimation from a WE run, which plateaus to a steady-state value equal to the reciprocal MFPT. (C) Realistic flux curves exhibit significant variance among different WE runs (run-to-run variance).
• Figure 2: Optimizing WE bins. Training data for optimization is obtained from initial WE simulations using arbitrary bins (blue bin boundaries). After classification of phase space into states (purple and blue boundaries), a history-augmented Markov state model (haMSM) is constructed from the training data and used to estimate the local MFPT values for each state. The states are subsequently combined in an optimal way (green boundaries) to minimize variance in estimating the "global" MFPT from initial to target state.
• Figure 3: 1D double-well potential function $V(x)$ and equilibrium density ${\rm e}^{-\beta V(x)} / \int {\rm e}^{-\beta V(y)} dy$, plotted for $\beta = 1$. The axes are scaled for illustration purposes.
• Figure 4: One-dimensional double well and functions used in variance minimization, plotted for $\beta = 1$. (A) Nonequilibrium steady-state distribution. The distribution $\pi$ (blue) is negligible near the sink, in contrast to the equilibrium steady-state distribution (Fig. \ref{['fig:1d_V_and_eq']}). (B) Discrepancy. The discrepancy $h$ (orange) illustrates that the local MFPT is larger near the source and smaller near the sink. (C) Variance. The variance $v^2$ (green) peaks near the energy barrier where trajectories might move to either side. (D) Optimal allocation. The optimal allocation $\pi v$ (red) is maximal on the ascending side of the barrier closer to the source. All functions are scaled to have minimum value 0 and maximum value 1. For clarity, the potential is not plotted to scale.
• Figure 5: Variance reduction applied to WE simulations of synMD Trp-cage unfolding. The discrepancy functions in panels (A) to (C) have been rescaled to lie in the interval [0, 1] for visibility. (A) Unoptimized WE bins indicated as colored bands with synMD states (dots) plotted using the RMSD and discrepancy function $h$ as coordinates. (B) Exactly optimized WE bins. (C) haMSM optimized WE bins with cluster centers indicated as dots. (D) Representative instantaneous flux profiles compared to the exact reference value (dotted line). (E) Representative cumulative flux profiles. (F) Statistical summary of variation in flux values, based on 20 runs for each scheme. In the box-and-whisker plots, the horizontal orange line indicates the median while the top and bottom boundaries indicate Q3 and Q1 values, and the height of the box corresponds to the interquartile range (IQR). The top and bottom whiskers correspond to Q3 + 1.5 IQR and Q1 - 1.5 IQR, and the black circles correspond to outliers outside the whiskers.