多肽的分子动力学模拟：力场选择、溶剂模型与构象采样的实战经验

多肽的分子动力学模拟处于蛋白质模拟和有机小分子模拟之间的一条窄路上。蛋白质模拟有成熟的全原子力场(AMBER ff19SB、CHARMM36m)和粗粒化方案(Martini)作为底层支撑，有机小分子有GAFF、CGenFF可以补位，多肽恰恰是两者的过渡地带——短肽在溶液中高度柔性，构象空间比折叠蛋白大几个数量级，而力场的参数化精度在多肽的非标准氨基酸残基和非天然侧链上打了折扣。

力场选择的物理限制

AMBER和CHARMM是两条主要路线。AMBER系列(ff14SB、ff19SB)的多肽主链二面角参数经过了量子力学(QM)计算的重新拟合，对α-螺旋和β-折叠的倾向性描述明显优于早期版本(ff99SB)。CHARMM36m在2016年更新了多肽骨架的CMAP修正项，解决了早期CHARMM22/CMAP版本中”螺旋偏好”——某些序列明明实验上是无规卷曲(random coil)，模拟里却自发折叠成了α-螺旋。

模拟短肽(5-15个残基)时遇到过一个经典问题：AMBER ff14SB对丙氨酸二肽(alanine dipeptide)的Ramachandran分布与溶液NMR J-coupling数据吻合得很好，但对含D-氨基酸或N-甲基化残基的多肽，力场参数的外推能力骤降。原因是AMBER力场的二面角参数是针对天然L-氨基酸拟合的，手性反转或甲基化修饰改变了二面角的旋转能垒，而力场里的参数没有重新训练。在模拟含有非天然氨基酸的短肽时，最好先跑一轮QM优化(HF/6-31G*或DFT/B3LYP/6-31G*水平)校准关键二面角的旋转能垒，然后用修正后的参数启动MD。这不是标准流程里的必选项，但在精度敏感的构象分析中，跳过这一步后面对构象分布的讨论都站不住。

CHARMM36m对多肽侧链的处理方式略有不同——它的侧链扭转参数经过了更细致的QM势能面扫描，对Tyr、Trp等大侧链的旋转异构体分布描述更准确。代价是CHARMM的参数格式和GROMACS的拓扑文件转换过程可能引入格式适配上的小错误——特别是Urey-Bradley项和CMAP项在不同的GROMACS版本中支持程度不一样，需要手工检查转换后的itp文件。

溶剂模型：显式还是隐式，不只是速度之争

多肽的分子动力学模拟中，溶剂模型的选择直接影响构象的稳定性。两派各有道理：

**显式溶剂模型**(TIP3P、TIP4P-Ew、OPC)的物理图像直接——水分子一个一个是存在的，氢键网络、疏水作用、溶剂介电屏蔽都可以在模拟中自然涌现。代价是计算量：一个20残基多肽 + TIP3P水盒子 + 150mM NaCl(约8000~15000个水分子)，系统总原子数3万~5万，在GPU上一天大约跑100~200ns(GROMACS/AMBER，RTX 3090级别)，一周不到1μs。对于多肽来说，1μs不一定够——短肽的α-螺旋-无规卷曲转变的特征时间在微秒甚至毫秒量级，1μs里的构象变化可能是局部的而不是全局的。

**隐式溶剂模型**(GBSA、PBSA)把溶剂的静电屏蔽效应用连续介质模型代替，氢键网络用解析函数近似，疏水作用通过溶剂可及表面积(SASA)来估算。它的速度优势巨大——没有水分子，系统原子数从3~5万降到几百到几千，一天能跑几微秒甚至几十微秒。但代价也很清楚：GB模型对带电残基(Lys、Arg、Asp、Glu)之间的盐桥稳定性描述远不如显式溶剂——GB模型中两个带相反电荷的侧链会过度稳定，导致构象被锁在非物理的紧密结构中出不来。

一个折中方案在几个项目中被验证有效：先用GB隐式溶剂跑一轮快速构象采样(1~5μs)覆盖构象空间的主要盆地，筛选出5~10个代表性构象(聚类时RMSD截断取2~3Å)，然后对这组构象分别用显式溶剂跑200~500ns的精算模拟——隐式给广度，显式给精度。这个策略把隐式溶剂的速度优势和显式溶剂的物理精度用在了它们各自擅长的地方。

增强采样：当你的模拟时间不够时该怎么办

多肽的分子动力学模拟中最头疼的问题不是参数，是采样——构象空间太大，模拟时间永远不够。增强采样方法可以在不依赖更长时长的情况下加速构象探索：

**REMD(副本交换分子动力学)**：用多个副本在不同温度下并行模拟，每隔一定步数尝试交换相邻副本的温度。高温副本能翻越构象能垒(因为热运动更强)，低温副本能在能量极小值附近精细采样。REMD对多肽的效果优于对折叠蛋白——因为多肽的能垒通常更矮更宽，温度复制之间构象重叠更好，交换接受率更高(一般>20%就有效)。代价是需要多个副本并行——一般16~48个副本覆盖300~450K的温度范围，GPU资源配置翻倍。

**Metadynamics**：通过在自由能面上逐步”填坑”(添加高斯势)来迫使系统离开已经访问过的构象盆地。对于多肽的φ/ψ二面角空间探索，选1~2个关键的二面角做集合变量(CV)，用Metadynamics可以在几百纳秒内画出完整的自由能面。风险在于如果CV选错了——比如用端到端距离做CV来推动构象变化，但实际的速率限制步是某个特定二面角旋转——Metadynamics会在不相关的方向上浪费计算资源，自由能面不会显著改变。

**aMD(加速分子动力学)**：在二面角势能面上施加一个”提升势”——当系统处于低能构象时抬高能量面，减小能垒高度，加速构象转变。aMD的参数只需要两个(E_thresh和α)，比Metadynamics简单，但自由能矫正(reweighting)的误差偏大，特别是当提升势在几个kJ/mol以上时。

三个方法的取舍：REMD给最接近物理的自由能估计(因为没有引入偏置势)，代价是计算量最大;Metadynamics适合采样有限的集体变量空间，但对CV选择敏感;aMD最易上手、不需要选CV，但自由能结果的定量精度偏低。多肽模拟选哪个，取决于最终目标是”获得正确的自由能数值”(选REMD)还是”快速排除不可能的构象区”(选aMD或Metadynamics)。

多肽聚集：从单分子到多分子多肽的步进

单个多肽在溶液中的构象行为，和多个多肽在一起时的聚集行为，是两个不同尺度的物理问题。

多肽的分子动力学模拟从单分子跨到多分子多肽聚集时，需要考虑浓度的物理对应。模拟盒子里放了4个多肽 + 水，折算的浓度大约在50~100mM——这远高于生理浓度(通常μM级别)，这意味着模拟里看到的聚集趋势会被高估——多肽分子之间碰撞的频率被浓度放大了10³~10⁴倍。

一个补救措施是用coarse-grained(CG)模型在低浓度下跑大规模系统——Martini力场能将4个重原子合并成一个作用珠(bead)，系统尺寸可以从几十个残基的多肽扩展到数百个多肽的系统，在微秒级模拟中观察到聚集的成核-生长过程。CG模型的缺点是丢失原子分辨率——氢键的几何定向性、π-π堆叠的精确角度这些驱动聚集的结构细节被粗粒化后无法准确再现。所以CG模拟适合回答”会不会聚集、形成什么形态的聚集体”，不适合回答”聚集界面上哪个残基是锚点”。

全原子和粗粒化结合，是这个领域的常见打法：CG先跑一遍大型系统确认是否有聚集倾向和聚集形态，再把CG代表性构象反向映射(backmapping)到全原子模型进行精细采样，确认关键残基间的相互作用模式。

分析方法：不要只有RMSD和RMSF

多肽的分子动力学模拟的结果分析如果只停留在RMSD(root mean square deviation)和RMSF(root mean square fluctuation)上，就浪费了模拟产生的大量信息。

**聚类分析**：用RMSD或二面角距离将轨迹按构象相似性聚成若干个簇，每个簇对应一个构象子态。只报平均结构没有意义——多肽在溶液中不是一个构象，而是一个构象系综(conformational ensemble)，平均结构可能是一个物理上不存在的虚拟构象。

**二级结构时间线**：用DSSP(Define Secondary Structure of Proteins)算法对每一帧的残基进行二级结构分类。在多肽的构象变化分析中，更关注的是二级结构的”寿命”——某个α-helix在轨迹中存在时间的长短，以及它消失-重现之间的间隔。

**氢键网络分析**：多肽分子内氢键的形成和断裂是构象转变的直接驱动力。统计哪些残基对之间的氢键占比较高(occpancy>30%)、哪些是短寿命的瞬时氢键(<10%)，能从分子层面上解释”为什么这个多肽偏向于α-螺旋”——因为i→i+4主链氢键的占有率远高于i→i+3(3₁₀-helix)或i→i+5(π-helix)。

**溶剂可及表面积(SASA)**：SASA的时间序列能捕捉疏水塌缩——疏水残基(Leu、Ile、Val、Phe)从暴露→包埋的转变在多肽折叠中起关键作用。SASA下降的快慢反映了疏水核形成的时间尺度。

多肽的分子动力学模拟从选择力场开始，到解读分析结果结束，每一个决策节点的选择都在约束最终结论的可靠性。这份可靠性不能靠一篇Nature Methods的benchmark论文来担保——它应该靠对自己系统的针对性验证来建立。