Weighted Ensemble Milestoning

• WEM is a computational strategy that combines weighted ensemble sampling with milestoning to accurately compute rare-event kinetics and thermodynamics in complex molecular systems.
• It partitions configurational space into milestones and employs resampling to reduce simulation time while estimating MFPTs, rate constants, and free energy profiles.
• Extensions like WEM-RR and analytical corrections for diffusion enhance sampling diversity and adjust rates, making the approach valuable for drug discovery and biomolecular studies.

Weighted Ensemble Milestoning (WEM) is a computational strategy that synthesizes the milestoning approach with the weighted ensemble (WE) path sampling method to compute both thermodynamic and kinetic observables for rare event processes in complex molecular systems. By integrating the master-equation rigor of milestoning with the statistical efficiency of the WE algorithm, WEM achieves significant reductions in computational cost for simulating transitions that typically occur on timescales inaccessible to standard molecular dynamics (MD). This approach provides estimates for quantities such as mean first-passage times (MFPTs), committor probabilities, free energy profiles, rate constants, and time-correlation functions, with demonstrated applicability to both model systems and realistic biomolecular complexes (Ray et al., 2020, Ray et al., 2019).

1. Theoretical Framework

WEM partitions configurational space along a user-chosen reaction coordinate (RC) $\xi(x)$ into a set of non-overlapping hypersurfaces referred to as milestones $\{\xi_0, \xi_1, ..., \xi_M\}$. Trajectories are initiated at each milestone and propagated until they reach a neighboring milestone. Two core statistical entities are accumulated:

• Milestone-to-milestone transition kernel $K_{ij}$: Probability that a trajectory starting on milestone $i$ first exits on neighboring milestone $j$ ($j = i \pm 1$),

$$K_{ij} = \sum_{k \in \Gamma(i \to j)} w_k, \quad K_{ij}=0 \;\text{for}\; j \neq i\pm1$$

where $w_k$ are statistical weights from WE, and $\Gamma(i \to j)$ is the set of trajectories starting on $i$ and first hitting $j$.

• Mean lifetime at milestone $i$, $\overline{T}_i$:

$$\overline{T}_i = \sum_{k\in\Gamma(i\to\{i-1,i+1\})} w_k\,t_k$$

with $t_k$ the dwell time before escape.

The stationary distribution of fluxes $q_i$ is determined via:

$q^T K = q^T,$

subject to appropriate boundary conditions. The equilibrium occupancy at milestone $i$ is:

$P_{eq,i} = q_i \overline{T}_i$

and the free energy profile referenced to milestone $0$ is:

$\Delta G_i = -k_B T \ln\frac{P_{eq,i}}{P_{eq,0}}$

The MFPT from reactant ($0$) to product ($M$) milestones follows:

$$\langle \tau \rangle = \frac{\sum_i q_i \overline{T}_i}{q_f}$$

where $q_f$ is the steady-state flux through the product milestone (Ray et al., 2020, Ray et al., 2019).

2. Algorithmic Workflow

WEM proceeds in the following main steps:

1. Milestone anchoring: Save a representative configuration $x_i^0$ from initialization MD at each milestone $\xi_i$.
2. WE binning: Subdivide the segment between adjacent milestones into WE bins.
3. WE iteration: Launch $N$ trajectories from $x_i^0$, each with weight $1/N$. Propagate for a WE cycle duration $\delta t$ under unbiased dynamics.
4. Resampling: After each $\delta t$, retain a fixed number of trajectories in each bin by splitting and merging, conserving total bin weight.
5. Milestone crossing accounting: If a trajectory exits to a neighboring milestone, record its exit statistics ($t_k, w_k$) and remove it.
6. Assembly: Upon convergence, aggregate $K_{ij}$ and $\overline{T}_i$ for all adjacent milestone pairs, solve the master equation, and extract kinetic/thermodynamic observables.

This workflow is highly parallelizable, as WE sampling within each milestone segment is independent (Ray et al., 2020, Ray et al., 2019).

3. WEM-RR: Restraint-Release Extension

The WEM-Restraint-Release (WEM-RR) protocol addresses insufficient initial condition sampling from a single anchor per milestone in standard WEM. It imposes a harmonic restraint on the RC at each milestone for $L$ WE iterations to generate a set of $L \times N$ distinct configurations $\{x_i^\ell\}$. After discarding the restraint, full WE runs are started from each anchor, dramatically increasing sampling diversity without additional atomic-scale cost. This reduces variance in both the transition kernel $K_{ij}$ and mean lifetime estimates $\overline{T}_i$.

Pseudocode for WEM-RR anchoring:

for each milestone i: restrain ξ(x) ≈ ξ_i; run L WE iterations → anchors {x_i^ℓ} for each ℓ in 1…L: initialize N trajectories at x_i^ℓ with w=1/(N L) run WE until K_{i,i±1}, T_i converge assemble global K, T → solve for kinetics & ΔG

(Ray et al., 2020).

4. Analytical Corrections for Diffusion-Limited Binding

For association events dominated by diffusive arrival at a boundary milestone $\xi_r$ with a nearly flat free energy profile, WEM supports an analytical correction to the computed $k_\text{on}$. This accounts for the effective fraction $\alpha$ of the outer milestone's surface that leads inward (computed via Monte-Carlo sampling), the ligand diffusion coefficient $D$, and the transition kernel element $K_{r, r-1}$:

$$k_\text{on}^\text{diff} = 4\pi D r \alpha K_{r, r-1} N_A$$

where $N_A$ is Avogadro's number. This correction enables accurate estimation of association rates in the diffusion-limited regime (Ray et al., 2020).

5. Benchmark Applications and Validation

WEM and WEM-RR have been validated on both model and biomolecular systems:

System	Observables	Result accuracy (vs reference)
Na$^+$/Cl$^-$ ion	Free energy barrier,	MFPT$_{unbind}$ ≈ 1 ns (1.25 ns brute-force MD);
pair in water	MFPT, $k_\text{off}$, $k_\text{on}$	$k_\text{off} \approx 1\times10^9$ s$^{-1}$, $k_\text{on} \approx 1\times10^9$ M$^{-1}$s$^{-1}$, $\Delta G \approx 0$ $k_BT$ (within <1 $k_BT$ of reference)
FKBP–BUT ligand	$\Delta G$, MFPT,	Standard WEM: large variance; WEM-RR ($L=10$): $\langle\tau_{off}\rangle \approx 9$ ns (cf. 7.2$\pm$3.0 ns MD), $\Delta G^{\ddagger}\approx8.3\pm0.3\,k_BT$, $k_\text{on}^{\text{diff}} \approx 0.7\times10^9$ M$^{-1}$s$^{-1}$, $K_d \approx 0.17$ mM, $\Delta G\approx1.8\,k_BT$ (closely matching long Anton MD and metadynamics)

In all cases, WEM achieves order-of-magnitude improvements in computational effort compared to conventional or WE-only simulations, with final results in agreement with reference MD and, where available, experiment (Ray et al., 2020).

6. Implementation, Efficiency, and Applicability

WEM's computational steps can be summarized as follows:

• Milestone selection: Milestones should be placed orthogonally to the slow reaction coordinate; finer spacing in high free energy barrier regions is beneficial.
• WE bins: Within each milestone cell, 5–20 bins are typical; binning enables granular resampling needed for efficient WE.
• Convergence: $|K_{i,i-1} + K_{i,i+1} - 1| < \epsilon_1$ is routinely used; stability in $K_{ij}$ and $\overline{T}_i$ needed before assembly.

Computational efficiency is high due to the embarrassingly parallel nature of individual WE runs between milestones. For example, FKBP–BUT dissociation calculations (both kinetics and thermodynamics) required less than 100 ns of aggregate simulation time across all milestones, in contrast to tens of μs for brute-force MD (Ray et al., 2020). No biasing forces are introduced, so physical mechanisms are preserved. Quantities accessible include committor probabilities, MFPTs, all rate constants (on/off), $K_d$, free energy landscapes, and full time-correlation functions.

WEM and WEM-RR are especially suited to high-throughput sampling in drug discovery, screening, and any application requiring rare-event kinetics well beyond the microsecond scale, as long as an appropriate low-dimensional RC can be identified (Ray et al., 2020, Ray et al., 2019).

7. Relation to Other Enhanced Sampling Strategies

WEM is distinguished from classical milestoning by its use of WE to efficiently estimate transition statistics within the milestone formalism. The synergy enables accurate thermodynamics and kinetics in fewer total MD steps versus either method alone. The master-equation-based assembly of observables ensures quantitative rigor, while WE acceleration mitigates the sampling bottlenecks endemic to direct trajectory-based methods. In contrast to single-trajectory brute-force MD and equilibrium WE methods, WEM leverages both statistical resampling and partitioned configuration space to realize faster convergence and accessibility to both equilibrium and dynamical quantities (Ray et al., 2019).

---

References:

• [Kinetics and Free Energy of Ligand Dissociation Using Weighted Ensemble Milestoning, (Ray et al., 2020)]
• [Weighted Ensemble Milestoning (WEM): A Combined Approach for Rare Event Simulations, (Ray et al., 2019)]