分子动力学模拟是一种以分子的经典力学(分子力场)模型为基础,通过数值求解分子体系的运动方程,研究分子体系的结构与性质的计算机模拟方法。
分子动力学模拟就像用计算机给原子和分子“拍动画片”——根据物理规则(比如原子间的吸引或排斥),计算每个原子每分每秒的位置和速度,把它们的运动过程连起来,就能看到微观世界里的“慢动作镜头”。例如,它能展示药物分子如何找到病毒弱点,或者蛋白与蛋白如何结合到一起,帮助研究者不做实验来预测物质作用机制。
相对于分子对接,分子动力学模拟可通过纳秒级或更长时间尺度的动态模拟来解析受体-配体复合物的构象动态变化、溶剂化效应及氢键网络重构等动态相互作用轨迹。并基于 MM/PBSA 或自由能微扰理论计算结合自由能,可用于验证对接构象的稳定性和准确性,优化受体-配体结合模式,提升受体-配体对接模型和分子识别机制的可靠性。
常用的分子动力学模拟软件包括 AMBER、CHARMM、GROMACS 等,其中 GROMACS 是一款开源软件,可以用于上百万个粒子体系的分子动力学模拟研究,适合大型系统和长时间尺度。
分子动力学模拟的通用流程可分为以下六个核心步骤,每一步都对最终结果的可靠性至关重要:
获取蛋白 3D 结构
获取蛋白的 3D 结构(PDB 格式)。如没有蛋白的 3D 结构,需要对蛋白质进行结构预测。
对化合物进行结构优化,对受体蛋白进行活性位点预测,为后续分子对接提供可靠输入。
从分子对接结果出发,获取受体-配体复合物的初始结合构象,作为动力学模拟的起点。
将蛋白进行预处理并生成拓扑结构;如果是受体-配体体系,也需将配体进行预处理并生成拓扑结构。
依次进行模拟体系溶剂化(如水溶剂体系)与能量最小化,找到局部势能最小值,消除不合理接触。
进行温度平衡(NVT 控温)与压强平衡(NPT 控压),使体系达到目标温度与压强下的稳定状态。
启动分子动力学模拟,采集纳秒级或更长时间尺度的运动轨迹,用于后续结构与能量分析。
根据客户需求进行数据分析,包括 RMSD、RMSF、Rg、SASA、氢键、MM/PBSA 等常规与个性化分析。
📌 如何选择分析类型?
基础分析适用于绝大多数常规项目,可快速评估体系稳定性与结合模式。进阶分析 1 适用于需要高精度结合自由能或相互作用本质解析的项目。进阶分析 2 适用于需要深入理解结合动态过程与结构变化的项目。
均方根偏差(Root mean square deviation, RMSD)表示某一时刻的构象与目标构象所有原子偏差的加和,是衡量体系是否到达稳定的重要依据,可用于量化蛋白质三级结构的稳定性。具体而言,RMSD值较小且稳定表明分子体系在模拟过程中保持了与参考构象相似的结构,显示其三级结构的稳定性。相反,RMSD值的显著变化则表明体系整体与初始状态发生较大的构象改变。
通常情况下,RMSD 值的差值(ΔRMSD)≤ 0.3 nm,是可接受的结果。ΔRMSD 小于 0.2 nm 被广泛认为是验证正确 model 的最有效阈值。
如下研究,45 ns 后 wt 蛋白和翻译后修饰的 PTM 蛋白构象收敛,表明 PTM 和 wt 蛋白在 100 ns 的分子动力学模拟中趋于稳定,适用于进一步分析。
图 2 · 蛋白 RMSD 随时间变化示意图(45 ns 后趋于稳定)
均方根波动(Root Mean Square Fluctuation, RMSF)反映整个动力学模拟中每个氨基酸的平均摆动幅度,通过比较分析各残基的RMSF变化,用于判断受影响区域,从而推断互作区域。通常,RMSF值高的区域代表较高的柔性(loop区或末端),RMSF值较低的区域则更刚性(α-helix,β-sheet)。RMSF分析通过比较残基的波动性,揭示分子的动态变化与功能的关系,用于指导突变,实验优化等。
如下图的RMSF分析结果比较中,PTM和wt的protein在 260–370 左右的氨基酸发生较大的波动,这一段区域是蛋白与配体的结合面,这表明PTM修饰能极大地影响蛋白互作面氨基酸的摆动幅度,从而影响蛋白与配体小分子或配体蛋白的互作。这与实验当中蛋白的翻译后修饰通常影响蛋白与配体结合的结果相一致,能很好地从分子层面解释翻译后修饰影响互作的机制。
图 3 · 蛋白每个残基的 RMSF 值(示意:260–370 区域波动显著)
回转半径(Radius of Gyration, Rg)是衡量蛋白质结构紧密度的重要参数,表示蛋白质中所有原子的质量中心与其几何中心之间的平均平方距离。Rg值越小,说明蛋白质结构越紧密;相反,Rg值越大,说明蛋白质结构越松散。
通过计算和分析Rg值,可以了解蛋白质在动力学模拟过程中整体构象的变化。较低的Rg值通常对应于较为紧密和稳定的结构;较高的Rg值则可能表明蛋白质结构较为松散或展开。Rg的变化趋势能够反映出蛋白质在不同条件下的构象稳定性。
如下图的结果,PTM修饰的蛋白较wt的蛋白有更高的Rg值,表明PTM修饰蛋白整体的结构更松散,与前面RMSF中PTM蛋白和wt蛋白局部改变相印证。PTM修饰蛋白的局部改变对整体球蛋白的结构紧密程度造成影响。
图 4 · 蛋白 Rg 对比示意图(PTM 修饰后结构更松散)
自由能地貌图(Free Energy Landscape,简称FEL)是用于描述分子体系在不同构象或状态之间自由能变化的图形化工具。通过将 Rg 和 RMSD 值这两个参数作为坐标,并通过统计方法计算对应的自由能分布,可以揭示分子在不同构象和状态下的稳定性和动力学信息。
根据图中的能量图可以看到,蓝色区域为体系中具有稳定构象的能量区域。这个区域表示体系在这些构象下自由能最低,显示了PTM蛋白结构在这些状态下的高稳定性。通过提取蓝色区域蛋白的构象,可识别出整个100ns的分子动力学模拟中,稳定的蛋白构象与配体分子的具体互作模式,对理解蛋白质与配体互作的具体模式和蛋白与配体的100ns的分子动力学模拟中的动力学行为具有重要意义。
图 5 · 自由能地貌图(2D和3D)示意:蓝色区域为低自由能稳定构象
对模拟轨迹进行周期消除后,结合前述RMSD等分析,提取蛋白复合物于 20 ns–100 ns 时间段的运动轨迹进行MM/PBSA结合自由能计算,共选取 1000帧 结构进行计算。
在这项分析中,可以计算得到以下几项数据:
ΔvdW (Van der Waals Energy):范德华相互作用能量,反映了分子之间的非极性相互作用;
ΔEEL (Electrostatic Energy):电荷之间的静电相互作用能量;
ΔEPB (Polar Solvation Energy):极性溶剂化能,利用 Poisson-Boltzmann 方程来计算分子在极性溶剂中的溶剂化能;
ΔESA (Non-Polar Solvation Energy):非极性溶剂化能,通常通过溶剂可接触表面积乘以一个经验参数得到;
ΔEtot (ΔGbind):结合自由能总值,用于评估分子之间的结合自由能。
为更深入地分析小分子与蛋白的相互作用及强度,本节计算了每个残基对结合自由能的贡献。
图 6 · 结合自由能残基分解示意(负值越大贡献越显著)
在药物设计过程中,往往需要对化合物的结构进行修饰改造,可以通过这个结合情况图谱来进一步分析化合物各部分的结合情况,从而进一步对化合物结构进行改造优化。与残基分解相似,指纹图谱统计分析了100 ns分子动力学模拟中每一帧的模拟结果中,蛋白与配体互作的氨基酸图谱。可用于辅助判断配体结合蛋白时是否发生“变形”,蛋白与配体互作细节。
图 7 · 小分子化合物与蛋白结合指纹图谱示意(颜色深浅代表互作频率)
以对称匹配微扰理论(Symmetry-Adapted Perturbation Theory, SAPT)为基础,将分子间相互作用能分解为不同物理项,从能量角度更深入地理解分子间相互作用的本质。相比MM/PBSA,SAPT是一种基于量子力学的分子相互作用分解方法,能高精度地分析蛋白与配体之间的相互作用。可用于解释蛋白与药物分子的高精度相互作用分析以及精确解释氢键形成的机制。
图 8 · 蛋白二级结构变化示意(颜色代表不同二级结构类型)
SAPT 基于量子力学,精度更高,可精确分解静电、交换、诱导、色散等能量项;MM/PBSA 基于分子力学,计算速度更快,适合大规模筛选。两者结合可获得更全面的相互作用图像。
在分子动力学模拟中,二级结构分析用于研究和监测蛋白质或核酸等生物大分子中的二级结构变化情况。通过不同的颜色表示每个残基在动力学过程中形成的二级结构变化。
图 8 · 蛋白二级结构变化示意(颜色代表不同二级结构类型)
研究药物分子与靶标蛋白的结合模式和亲和力,基于结合自由能计算和蛋白稳定性计算,预测药物分子结合蛋白的稳定性和蛋白突变耐药性的机制,为药物开发提供分子解释。
研究蛋白的折叠和去折叠的动态过程,研究其发挥功能过程中的构象改变(如酶活性中心的开闭、蛋白信号传导过程中的构象改变、阿尔茨海默病类疾病类的蛋白质异常聚集、不同温度/不同pH/不同盐浓度对蛋白活性的影响等),蛋白突变体对蛋白稳定性的影响等。
研究DNA/RNA的构象变化过程,可用于解释药物对核酸结构的干扰(抑制转录,增殖等药物机制的研究),核酸与蛋白质互作的研究(例如:转录因子与DNA互作研究,核酸损伤后修复酶作用机制),基因编辑分子机制研究(例如:Cas9靶向核酸的机制,表观修饰),核酸药物可行性研究等。
基于分子动力学模拟验证实验设计(蛋白突变体的设计,虚拟筛选分子的可行性),为后续实验提供理论依据与优先级排序。