A Residence-Time Approach for Determining Position-Dependent Diffusivities from Biased Molecular Simulations

Abstract: We introduce a residence-time approach (RTA) for determining position-dependent diffusivities from biased molecular dynamics simulations. The method is formulated for trajectory segments in which the effective drift along the transport coordinate is negligible, as realized here using adaptive biasing force simulations. In this regime, local diffusivities are obtained directly from mean first-exit times out of finite spatial intervals. Unlike conventional fluctuation-based approaches, the RTA does not require dedicated harmonically restrained simulations or numerical integration of noisy time-correlation functions. We assess the method for oxygen diffusion across a hexadecane slab, water permeation across a lipid bilayer, and permeation of water and selected volatile organic compounds through a model skin-barrier membrane. In the slab system, the RTA reproduces independently determined bulk diffusivities within statistical uncertainty. In the membrane systems, the inferred diffusivity profiles are supported by propagator-level validation. These results establish the RTA as a practical approach for extracting position-dependent diffusivities from biased molecular simulations.

• The paper introduces a novel residence-time approach (RTA) for determining local diffusivities from biased ABF simulations.
• It validates RTA across bulk hexadecane/water systems, POPC bilayers, and SC membranes, demonstrating strong agreement with MD-derived propagators.
• The method overcomes limitations of conventional estimators by providing robust statistical estimation and improved permeability predictions.

Residence-Time Approach for Local Diffusivity Estimation in Biased Molecular Simulations

Introduction and Context

Diffusivity profiles are indispensable for quantifying transport phenomena in heterogeneous environments such as lipid membranes and barrier systems. Reliable extraction of position-dependent diffusivities, $D(z)$, from MD simulations provides the thermodynamic and kinetic parameters necessary to bridge molecular-level phenomena with macroscopic permeability. Despite the development of a suite of methods—equilibrium fluctuation-based estimators, propagator inference, kinetic reconstruction, and moment-based techniques—numerical and conceptual challenges persist. Many approaches demand extensive sampling, are sensitive to analysis parameters, or are fundamentally restricted by the presence of noise, memory effects, and model assumptions.

This work introduces a residence-time approach (RTA) for extracting $D(z)$ from biased molecular simulations, exploiting the drift-free dynamics realized in adaptive biasing force (ABF) regimes. Unlike VACF or PACF methods, RTA obviates the need for restrained simulations and direct correlation function integration. Validation is performed with systematic analysis of prototypical systems: oxygen in a hexadecane/water slab, water permeation through a POPC bilayer, and the diffusion of water and volatile organics through a model skin-barrier membrane.

Theoretical Foundation: Residence-Time Estimator

The RTA is framed within the Smoluchowski formalism for projected stochastic dynamics, focusing on the regime where ABF has converged to a flat effective free energy. In this limit, projected dynamics along $z$ become drift-free,

$$\partial_t p(z,t) = \partial_z [ D(z) \partial_z p(z,t) ].$$

The key identity leveraged in RTA connects the mean first-exit time (MFET) from an interval $\Omega = [a, b)$ with width $L$ to the local diffusivity (assuming $D(z)$ is approximately constant over $L$):

$\tau_r = \frac{L^2}{12 D}.$

Given a time series $z(t)$, RTA samples starting frames within each interval and averages the empirically determined exit times, delivering a direct, statistically robust estimator:

$D(z)$0

This construction is valid over intervals in which both $D(z)$1 and the effective drift are approximately constant and negligible, respectively, and where Markovian behavior dominates.

Methodological Implementation

The model systems encompass increasing levels of complexity and heterogeneity, ranging from simple phase-separated slabs to biologically relevant, multicomponent lipid matrices. ABF simulations (NAMD3, Colvars, CHARMM36/CGenFF/TIP3P) are employed to drive convergence to flat PMFs, enabling robust segmentation into drift-free trajectories. Statistical reliability is addressed with blocking procedures for variance estimation. The practical choice of interval width ($D(z)$27.5 $D(z)$3) is justified via propagator-level validation, balancing local homogeneity with spatial resolution.

Validation and System-Specific Results

Bulk System: Oxygen Diffusion in Hexadecane/Water

ABF provides rapid PMF convergence (Figure 1), enabling local residence-time analysis. The RTA diffusivity profile shows well-defined plateaus in both aqueous and hydrocarbon phases.

Figure 1: Convergence of PMF profiles from ABF simulations of oxygen diffusion across the hexadecane/water slab.

Figure 2: Diffusivity profiles for oxygen permeation across the hexadecane/water slab. Symbols denote computed diffusivities; lines are guides to the eye.

Direct comparison with bulk MSD-based diffusivities establishes statistical compatibility, with deviations within 0.1 $D(z)$4—well inside practical tolerances. This anchors the RTA as quantitatively precise when reference values exist.

Membrane System: Water in POPC Bilayer

ABF provides kinetic acceleration and convergent PMFs (Figure 3). RTA-derived diffusivities are compared to VACF and PACF estimates (Figure 4). The RTA yields intermediate values, specifically at the membrane center.

Figure 3: Convergence of PMF profiles from ABF simulations of water permeation across a POPC bilayer.

Figure 4: Diffusivity profiles across the POPC lipid bilayer obtained from the VACF, PACF, and residence-time approaches (RTA).

Propagator-level analysis (solving the Smoluchowski equation with inferred $D(z)$5 and $D(z)$6) provides stringent, lag-time-resolved validation. RTA-based propagators best reproduce central-membrane MD propagators at intermediate lag times ($D(z)$7 ps), while VACF and PACF estimates under- and over-broaden, respectively (Figure 5). At very short or long lags, discrepancies signal residual non-Markovianity and the limitations of one-dimensional descriptions.

Figure 5: Comparison of MD-derived propagators with model predictions obtained using diffusivity profiles from the VACF, PACF, and residence-time approaches (RTA) for water permeation across the POPC bilayer ($D(z)$8); representative lag times shown.

Barrier Membrane: SC Model with Water and VOCs

The SC system’s compositional heterogeneity and slow relaxation pose challenges for diffusivity estimators. RTA analysis is performed only on well-converged ABF windows; the resulting diffusivity profiles demonstrate that RTA tracks PACF values closely for acetone/6-MHO, and consistently lies between PACF and VACF for water (Figure 6).

Figure 6: Position-dependent diffusivity profiles across the SC membrane for (a) water, (b) acetone, and (c) 6-MHO. Symbols denote computed diffusivities; lines are guides to the eye.

Propagation analysis reveals that RTA-based $D(z)$9 produces propagators in close agreement with unbiased MD for water, especially at the membrane center and over all lag times tested, outperforming alternative estimators (Figure 7).

Figure 7: Comparison of MD-derived propagators with model predictions for water permeation across the SC membrane, using VACF, PACF, and RTA-based diffusivity profiles ($z$0).

Transport Coefficients: Permeability Predictions

Permeability coefficients computed from inhomogeneous solubility–diffusion theory reflect integrated PMF--diffusivity pairings (Figure 8). VACF-based permeabilities are always maximal; the interplay between RTA and PACF varies by system and solute. For water in POPC, the hierarchy is $z$1; for SC-acetone, PMF differences exert dominant influence.

Figure 8: Permeability coefficients computed using diffusivity profiles from VACF, PACF, and RTA methods.

Implications and Outlook

RTA provides a robust, practical, and conceptually transparent method for extracting $z$2 from biased MD simulations, exploiting drift-free segments of ABF trajectories. It circumvents the restraint-strength and lag-time dependencies intrinsic to correlation-based estimators, eliminates the need for dedicated sampling, and supports statistically rigorous uncertainty quantification. The presented applications demonstrate that RTA:

• Reproduces bulk reference diffusivities in homogeneous phases within uncertainty bounds,
• Yields intermediate and often more physically accurate descriptions in lipid bilayers relative to established methods,
• Provides best-in-class agreement with MD propagators across biologically relevant time scales in complex barriers.

Importantly, the results highlight that a single lag-time-independent diffusivity profile cannot generally capture projected transport across all dynamical regimes; system-specific and solute-specific non-Markovian effects remain a limiting factor in reduced descriptions.

Conclusion

The residence-time approach stands as a practically attractive, computationally efficient estimator for position-dependent diffusivities in biased MD studies. Across a series of validation systems—ranging from simple slabs to highly heterogeneous, ordered lipid matrices—RTA yields diffusivity and permeability values that are consistent with, or superior to, established estimators when judged at the propagator level.

Immediate practical recommendations include adoption of RTA in studies where precise permeability estimation or physical interpretability of $z$3 is paramount, and its coupling with robust ABF protocols. Future directions should systematically scrutinize interval-width dependence, extend validation to more diverse solute classes, and develop hybridized models that explicitly address residual, lag-dependent non-Markovianity. This approach establishes a foundation for both methodological innovation and the routine, robust extraction of kinetic parameters from trajectory data in complex soft-matter systems.

Reference: "A Residence-Time Approach for Determining Position-Dependent Diffusivities from Biased Molecular Simulations" (2604.01940).