High-Dimensional Enhanced Sampling via Regularized Path-Dependent McKean--Vlasov Dynamics using Tensor Density Approximation

This paper proposes a scalable, regularized path-dependent McKean-Vlasov framework for high-dimensional enhanced sampling that improves statistical stability through path-history measures and achieves efficient numerical realization via optimization-free tensor density approximation, enabling effective exploration of complex energy landscapes with collective variable dimensions up to 64.

The Gist

Imagine you are trying to explore a vast, foggy mountain range to find every hidden valley and peak. This mountain range represents the "energy landscape" of a molecule. In a standard simulation, the molecule is like a hiker who gets stuck in one deep valley (a "metastable state") because the mountains around it are too high to climb over. The hiker just walks around in that one valley for a long time, never seeing the rest of the world.

Scientists want to see the whole map, but the hiker is too slow and the mountains are too high. This is the problem of sampling: getting a complete picture of a complex system without waiting for an impossible amount of time.

Here is how this paper solves that problem, using simple analogies:

1. The Old Way: The "Instantaneous" Map

Previous methods tried to help the hiker by drawing a map of where they have been right now and telling them, "Go where you haven't been yet!"

• The Problem: If you only have a few hikers (which is usually the case in computer simulations), the map they draw is very shaky and full of holes. It's like trying to draw a detailed map of a city based on the path of a single person walking for five minutes. The map is too noisy, and the instructions become confusing.
• The Math Issue: To make the map smooth enough to follow, old methods had to do a massive amount of complex math (called "convolution") that becomes impossible to calculate when the mountain range has many dimensions (like 64 different directions to move).

2. The New Solution: The "Memory" Hiker

The authors propose a new way to guide the hiker. Instead of looking at where the hiker is right this second, they look at the entire history of the hiker's journey.

• The Memory Trick: Imagine the hiker carries a backpack that remembers every step they've taken over the last hour. The guide looks at this full history to decide where to push the hiker next.
• Why it helps: Even if you only have a few hikers, their history is long. By averaging over time (the path) rather than just counting how many hikers are in a spot right now, the map becomes much smoother and more reliable. This allows the simulation to work well even with a small number of computer "walkers."

3. The "Smart" Compass (Regularization)

The new method also fixes a "roughness" problem. If the hiker's history shows a tiny, empty spot, the old math might get confused and say, "Go there!" or "Don't go there!" in a jerky, unpredictable way.

• The Fix: The authors added a "smoothing filter" (called regularization). Think of it like a smart compass that refuses to give a direction if the data is too shaky. It gently nudges the hiker away from crowded areas and toward empty ones, but it does so smoothly so the hiker doesn't get jolted around. This makes the math stable and prevents the simulation from crashing.

4. The "Folding" Map (Tensor Density)

The biggest challenge is that the mountain range has 64 dimensions. Imagine trying to draw a map of a city where you need to track 64 different variables at once (temperature, wind, humidity, traffic, etc., all at the same time). A normal grid map would require more paper than exists in the universe to draw this.

• The Solution: The authors use a technique called Functional Hierarchical Tensor (FHT).
• The Analogy: Instead of trying to draw the whole 64-dimensional map on one giant sheet of paper, they break the map down into smaller, connected pieces that can be "folded" together efficiently. It's like packing a complex 3D object into a flat suitcase by folding it in a specific, smart pattern. This allows them to store and calculate the map of the 64-dimensional world without needing a supercomputer to run out of memory.

5. The Results: Exploring the Unexplored

The team tested this method on several "mountain ranges":

• Simple Hills: A 2D test case where they could see the whole map.
• Peptides: Small protein chains with 3 to 9 moving parts.
• Proteins: Real biological molecules.
  • Chignolin: A small protein with 16 moving parts.
  • Villin Headpiece: A slightly larger protein with 64 moving parts.

The Outcome:
In standard simulations, the hiker would get stuck in the "native" folded shape of the protein and never unfold. With this new method, the hiker successfully explored the entire landscape, finding the folded state, the intermediate states (half-folded), and the fully unfolded states. They were able to do this even with 64 dimensions, a scale that was previously considered too difficult for these types of adaptive sampling methods.

Summary

The paper introduces a new way to simulate molecules by:

1. Using memory: Looking at the whole journey history instead of just the current moment to get a smoother, more reliable guide.
2. Smoothing the path: Adding a filter to prevent the guide from giving confusing instructions in empty areas.
3. Folding the map: Using a smart mathematical "folding" technique to handle maps with up to 64 dimensions, which was previously impossible.

This allows scientists to see the full "mountain range" of complex molecules much faster and more accurately than before.

Technical Summary

Technical Summary: High-Dimensional Enhanced Sampling via Regularized Path-Dependent McKean–Vlasov Dynamics using Tensor Density Approximation

Problem Statement
Sampling from high-dimensional Gibbs measures associated with complex energy landscapes is a central challenge in computational statistical mechanics. When the potential energy $U(\mathbf{x})$ contains multiple well-separated local minima, standard Langevin dynamics often fails to explore the configuration space within feasible simulation times due to rare barrier-crossing events. Enhanced-sampling methods address this by introducing adaptive biasing potentials along prescribed collective variables (CVs) to flatten the free-energy landscape. However, existing approaches face two critical scalability limitations in high-dimensional CV spaces ($m \sim 10\text{--}100$):

1. Numerical Intractability of Regularization: Recent Wasserstein-gradient-flow formulations (e.g., Lelièvre, Lin, and Monmarché [19]) require a nested convolution to regularize the free-energy functional. The outer convolution over the CV domain necessitates a grid-based representation that becomes computationally prohibitive as the CV dimension $m$ increases.
2. Finite-Ensemble Instability: Practical molecular dynamics (MD) simulations typically employ a small number of replicas ($M \sim \mathcal{O}(10)$). The instantaneous empirical marginal density derived from such small ensembles is noisy and unstable, leading to unreliable biasing drifts when approximating the mean-field limit.

Methodology
The authors propose a reformulation of the adaptive biasing process as a regularized, path-dependent McKean–Vlasov stochastic differential equation (SDE). This approach decouples the bias construction from the exact Wasserstein gradient flow structure to achieve numerical tractability and stability.

1. Direct Regularization of the Drift: Instead of regularizing the free-energy functional (which induces the nested convolution), the method directly regularizes the CV marginal density entering the biasing drift.

   • The CV marginal density is first smoothed via mollification ($\varphi_\delta * \mu_\xi$).
   • A second regularization using a Softplus function, $K_{\tau, \epsilon}(r) = \epsilon + \tau \text{Softplus}(r/\tau)$, ensures a uniform positive lower bound.
   • This yields a drift term of the form $\nabla \log K_{\tau, \epsilon}((\varphi_\delta * \mu_\xi)(\xi(\mathbf{x})))$, which is globally Lipschitz and avoids the outer convolution over the CV space.

2. Path-Dependent Formulation: To mitigate the noise inherent in small-replica ensembles, the instantaneous law $\mu_t$ in the McKean–Vlasov drift is replaced by a weighted history measure $\mu^q_t$.

   • $\mu^q_t$ accumulates the trajectory history up to time $t$ weighted by a probability measure $q$ on $[0,1]$.
   • This transforms the dynamics into a path-distribution-dependent SDE, effectively trading a noisy ensemble average for a time average along the sample path.

3. Tensor-Based Density Estimation: The history-averaged CV marginal density is approximated using an optimization-free Functional Hierarchical Tensor (FHT) representation.

   • The density is expanded in a tensor-product basis of Gaussian functions.
   • The coefficient tensor is compressed into a hierarchical low-rank format (binary tree structure) using sketching techniques.
   • This allows the density and its score function to be evaluated efficiently with linear scaling in the CV dimension $m$, provided the density admits a low-rank structure.

4. Theoretical Guarantees: The authors establish the well-posedness (existence and uniqueness of strong solutions) for both the regularized McKean–Vlasov dynamics and its path-dependent generalization under suitable regularity assumptions. Furthermore, under strong dissipativity conditions, they demonstrate that the path-dependent dynamics is asymptotically consistent with the invariant measure of the corresponding instantaneous-law dynamics.

Key Results
The method was validated on benchmark systems with CV dimensions ranging from 2 to 64:

• Müller Potential (2D): The method successfully facilitated transitions between metastable basins. The reweighted free-energy surface (FES) accurately reproduced the reference landscape, with errors concentrated only in sparsely sampled high-energy regions.
• Alanine Systems (2D and 4D): For alanine dipeptide (2D) and ACE-(ALA)$_2$-NME (4D), the adaptive bias significantly accelerated convergence compared to unbiased MD. The method recovered principal wells and barriers with reduced run-to-run variability.
• Peptoid Systems (3D and 9D): Applied to s-(1)-phenylethyl peptoid (3D) and its trimer (9D), the method demonstrated robust exploration of the torsional landscape. The FHT estimator produced smooth, multimodal free-energy profiles without spurious oscillations, even in 9 dimensions.
• Protein Benchmarks (16D and 64D):
  • Chignolin (16 CVs): The method drove the system from the native $\beta$-hairpin to unfolded states, resolving four distinct metastable basins (folded, intermediate, and two unfolded states) in the reconstructed FES.
  • Villin Headpiece (64 CVs): This represents a significant scalability test. The method successfully constructed an adaptive bias in a 64-dimensional space, resolving the folded, intermediate, and unfolded basins of the 36-residue protein. The reweighted FES projections clearly separated these states, demonstrating the method's ability to handle realistic, high-dimensional protein folding landscapes.

Significance and Claims
The paper claims to provide a scalable adaptive enhanced-sampling framework capable of handling CV dimensions (up to 64) that are typically beyond the reach of standard grid-based or mean-field enhanced sampling methods.

• Novelty: The work introduces a path-dependent McKean–Vlasov formulation that replaces the instantaneous law with a history measure, specifically addressing the "small-replica" regime common in molecular simulations.
• Efficiency: By avoiding the outer convolution through direct regularization and utilizing FHT approximations, the method achieves computational tractability in high dimensions.
• Distinction: Unlike tensor-compressed metadynamics (TT-metadynamics), which compresses accumulated Gaussian hills, this method uses low-rank tensors as a density estimator to close a regularized score-driven stochastic dynamics.
• Limitations and Future Work: The authors modestly note that the direct regularization sacrifices the exact Wasserstein gradient-flow structure, meaning existing variational convergence theories do not directly apply. They identify the need for future work on convergence theory for realistic molecular potentials and the adaptive selection of tensor ranks and history weights. The FHT approximation is explicitly described as a surrogate for constructing the bias, not as the final unbiased FES estimator, which is recovered via reweighting or MBAR.