核酸结构的分子动力学模拟：从双螺旋到配体结合的动态路径

做核酸分子动力学模拟的人，第一次建体系的时候大概率会被离子中和这一步卡住。DNA在溶液里浑身带负电，磷酸骨架每两个核苷酸就贡献一个负电荷，一条12碱基的双链DNA净电荷轻易超过-20e。直接扔进水盒子里跑MD，系统撑不过100 ps——静电排斥让双链直接散架。

这不是你的模拟设置错了，是力场在提醒你：真实溶液里这些负电荷是被抗衡离子屏蔽的，你的模拟盒子也必须把它们加进去。

力场选择：AMBER还是CHARMM

核酸力场的水比蛋白质力场还深。AMBER的DNA力场从ff94一路迭代到最新的OL15，每个版本都在修正糖苷键二面角的偏好。CHARMM的DNA力场（CHARMM36）对A-form和B-form的平衡有更好的描述。

实际选型有个简单的判据：如果你的模拟里DNA主要是标准B-form双链，AMBER OL15是目前最稳的选择；如果涉及非标准构象（Z-DNA、三链、G-四链体），CHARMM36的覆盖度更好。

RNA的情况更复杂。RNA的构象空间比DNA大得多，碱基堆积和氢键模式也更丰富。AMBER的RNA力场（OL3）是目前最广泛验证的，但仍有已知的偏差——比如U-turn构象的能量评估偏高，需要手动修正 glycosidic torsion 的参数。

建系统时的关键步骤

第一步：结构准备。从PDB下载的结构，先检查有没有缺失的碱基或修饰。很多晶体结构里只给了电子密度可见的碱基，末端几个碱基密度弱，模型里直接砍掉了。MD模拟里缺碱基会导致局部结构不稳，需要先用建模工具（LEaP for AMBER，psfgen for CHARMM）补全。

第二步：溶剂化和离子中和。SPC/E和TIP3P水模型都能用，TIP3P和AMBER力场配套更好。离子中和有个细节：不能只加够电中性的Na+或Cl-，还要加够生理盐浓度（~150 mM NaCl）。只做电中性中和，DNA的沟渠里离子分布不真实，会影响后续配体结合的自由能计算精度。

第三步：能量最小化。先做500-1000步的最陡下降，再做共轭梯度，直到最大力小于1000 kJ/mol·nm。如果体系里有晶体水，最小化的步数要加倍，晶体水的位置往往和MD平衡结构有偏差，需要慢慢弛豫。

第四步：平衡。NVT平衡先把温度拉到指定值（通常300K），用Langevin恒温器或者v-rescale。NPT平衡再拉压强到1 atm。核酸体系建议各做500 ps，步长2 fs（用constraint的话可以4 fs）。

轨迹分析能告诉你什么

RMSD和RMSF：RMSD看整个结构是否稳定，核酸的RMSD通常在2-4 Å之间波动，超过6 Å就要怀疑结构是否散了。RMSF看哪个碱基柔性大，通常末端碱基的RMSF高，结合口袋里的碱基如果RMSF异常高，说明这个位点构象不稳定。

氢键分析：核酸的氢键分两种——Watson-Crick配对氢键和沟渠里的水分子介导氢键。GROMACS的gmx hbond可以直接算，但更细的分析需要自己写脚本：每个碱基对的氢键数目随时间变化，能告诉你配体结合是否破坏了局部的氢键网络。

Minor groove宽度：这个量对配体识别非常关键。AT-rich区域的minor groove天然比较窄，小分子如果锚定在这里，通常会和沟底的一串水分子形成有序结构。计算minor groove宽度可以用mdanalysis或cpptraj，沿着螺旋轴扫描每个碱基步的P-P距离。

自由能计算：MM-PBSA和热力学积分

配体结合的ΔG，最常用MM-PBSA近似。GROMACS里有gmx_MMPBSA工具，AMBER里有全套的mm_pbsa.pl和MMPBSA.py。

MM-PBSA的好处是快，从MD轨迹里取样几百帧就能给出结合自由能的近似值。缺点是误差大（通常3-5 kcal/mol），而且对熵的贡献估算粗糙。如果要做决策级的结合能比较，需要做热力学积分（Thermodynamic Integration）或者FEP，算起来贵10-100倍，但精度可信得多。

常见坑点

坑一：模拟时间不够。核酸的构象变化时间尺度在微秒到毫秒级，100 ns的模拟往往只采样了局部构象。如果关心的是构象转变（比如DNA弯曲、配体诱导的构象变化），至少需要微秒级采样，或者用加速采样方法（Gaussian accelerated MD，GaMD）。

坑二：力场漂移。有些AMBER老版本力场（ff99）对DNA的A-form有系统偏差，会在模拟里把B-form DNA慢慢扭成A-form。确认你用的力场版本，OL15和bsc1是这个领域目前最稳的选择。

坑三：忽略离子特异性。Na+和K+在核酸体系里的表现不一样，K+更容易稳定G-四链体。如果用错了抗衡离子，G-四链体的拓扑结构在模拟里可能维持不住，会直接解折叠。

什么时候MD能回答你的问题

MD模拟最适合回答动态问题：配体结合和解离的路径是什么？哪个碱基的突变会破坏结构稳定性？药物分子在沟渠里是怎么调整的？

不适合回答的问题：绝对的亲和力度数（用实验测，MD只能给相对比较）；突变后的折叠自由能变化（用专用的ΔΔG预测工具，比如Rosetta）；大于毫秒的时间尺度事件（用Markov State Model降维，但建模成本很高）。

核酸分子动力学模拟的价值，不在于算出某个精确数值，而在于让你看到晶体结构里看不到的那个维度——时间。结构给了你一张静态的地图，MD告诉你这张地图上哪些路是通的、哪些是死的。