位置依赖扩散系数确定的停留时间法：基于偏置分子模拟的提取方法

《Journal of Chemical Theory and Computation》：A Residence-Time Approach for Determining Position-Dependent Diffusivities from Biased Molecular Simulations

　　

位置依赖扩散系数（position-dependent diffusivity，D(z)）是异质环境中分子输运约化随机描述的核心参数，但其从分子动力学（molecular dynamics，MD）模拟中可靠估计仍具挑战。研究人员提出一种停留时间法（residence-time approach，RTA），从偏置模拟中测得的首达（first-exit）统计经平均自由能梯度补偿后提取局域扩散系数。该方法被应用于氧在十六烷/水平板中的扩散、水跨POPC脂质双分子层的渗透，以及水与挥发性有机物跨模型皮肤屏障膜（stratum corneum，SC）的输运。在平板体系中，RTA扩散系数与独立确定的体相参考值一致；在膜体系中，基于RTA所得D(z)的传播子（propagator）预测在无偏MD传播子的较大滞后时间（lag-time）范围内重现，同时揭示某些情况下不存在单一与滞后时间无关的扩散系数剖面可捕获全时间尺度动力学。这些结果支持停留时间统计作为从偏置MD模拟确定有效位置依赖扩散系数的实用途径。

研究背景与动机

在异质体系（膜、纳米孔、多孔材料）中，沿输运坐标z的位形依赖扩散系数D(z)与平均力势（potential of mean force，PMF，F(z)）共同构成一维Smoluchowski描述的核心理念，并经由非均匀溶解-扩散（inhomogeneous solubility–diffusion，ISD）模型桥接至可实验测量的渗透率P。然而从MD轨迹可靠提取D(z)长期受困：基于谐振约束涨落的相关函数法（velocity autocorrelation function，VACF；position autocorrelation function，PACF）需数值积分噪声相关函数且受约束强度影响；基于短滞时间矩、马尔可夫态模型或贝叶斯推断的方法面临离散化、正则化及非马尔可夫记忆效应偏差。更根本地，高维运动投影到单坐标会引入残余记忆，使D(z)估计依赖滞后时间。现有方法常在“自由扩散”与“约束动力学”间给出不一致的有效摩擦。因此亟需一种不依赖谐波约束、不直接处理噪声相关函数积分、且能与增强采样自然衔接的D(z)估计框架。

研究人员以“漂移自由极限”为支点，提出停留时间法（RTA）：在自适应偏力（adaptive biasing force，ABF）模拟中，偏力收敛补偿PMF梯度，使有效自由能F_eff(z)=F(z)+V_b(z)近似平直，约化动力学退化为无漂移项Smoluchowski方程?_tp=D(z)?_z²p。在一维区间Ω=[a,b)内对常扩散系数D作首达时间分析，平均首达时间对初位置z₀为τ(z₀)=(z₀?a)(b?z₀)/(2D)，在平直有效PMF下初位置均匀分布ρ_Ω=1/L，区间平均停留时间τ_r=L²/(12D)，反解即得局域估计D_i=L_i²/(12τ_r,i)。该式精确成立于区间常D、无漂移、近马尔可夫三维投影条件。论文发表于《Journal of Chemical Theory and Computation》。

主要关键技术方法

研究人员采用CHARMM36脂力场与TIP3P水，对十六烷/水平板、POPC/水双分子层、SC多组分脂基质（胆固醇、木蜡酸、神经酰胺NS等比例）三体系在NAMD中跑ABF生产轨（CV为溶质质心z坐标，0.1 ?分箱，边界谐波势），取PMF平稳后的偏置段做停留时间统计；区间宽L=7.5 ?，用Flyvbjerg–Petersen块化法估不确定度。对照用谐波约束轨的VACF/PACF扩散估计。验证分两层：体相区与Einstein-MSD参考比；膜区将RTA的D(z)与ABF的F(z)代入Smoluchowski有限差分解算传播子p(z,t|z₀)，与无偏MD（撤伞形约束释放）传播子在多滞后时间比。

研究结果

4.1 氧在十六烷/水平板中的扩散

ABF的PMF于800 ns后平稳对称，取其后段做RTA。D(z)在水相与十六烷体相出现明确平台，RTA体相值与独立MSD参考值在0.1 ?² ps^?1内重合且误差棒交叠；界面处D(z)平滑过渡。说明RTA在已知体相基准下可复原正确D(z)。

4.2 水跨POPC双分子层渗透

将RTA-D(z)与ABF-F(z)送Smoluchowski解传播子，在从亚百皮秒到数纳秒滞后时间内与无偏MD传播子形态、峰位、展宽一致；与水溶剂中谐波约束VACF/PACF结果趋势吻合但在双分子层中心略低，反映约束动力学改变有效摩擦。证明RTA-D(z)支撑ISD渗透率计算。

4.3 水、丙酮、6-MHO跨SC膜输运

SC体系分七重叠ABF窗口采样。RTA-D(z)显示脂基质内D比水相低约1–2数量级，丙酮与6-MHO剖面随极性变化。传播子验证表明：单一D(z)剖面可在中滞后时间区间重现无偏传播子，但在极短与较长滞后时间偏离，提示投影动力学含滞后时间依赖的有效D(z)——并非RTA失效，而是约化Smoluchowski描述本身在不同尺度由不同有效输运系数描述。

讨论与结论翻译

研究人员讨论指出，RTA避开了谐波约束与相关函数外推，直接由首达统计提取D(z)，实现门槛低且与ABF采样无缝衔接；其代价是D_i应理解为“与停留区间相关的首达有效扩散系数”，而非唯一局域输运系数。在SC等强异质膜中，不存在单一与滞后时间无关的D(z)能全尺度捕捉动力学，这源于高维自由度的非马尔可夫投影。结论部分写明：RTA从偏置MD模拟提取的位置依赖扩散系数在体相参照下准确，在膜体系传播子层面自洽，停留时间统计是从偏置模拟获取有效D(z)的实用路线；未来工作可系统考察区间宽L对有效D(z)滞后时间依赖的定量映射。