提升Go?Martini中折叠蛋白质构象系综的扰动优化方法

《Journal of Chemical Theory and Computation》：Improving Conformational Ensembles of Folded Proteins in Go?Martini

本文综述了通过扰动优化方法改进Go?Martini模型，以准确复现原子模拟中折叠蛋白质的构象动力学。针对标准均匀Go?网络在捕获蛋白质本质子空间自由能景观方面的不足，作者提出了一种自动化的PoGo?优化框架。该框架基于微扰理论和粒子群优化算法，通过迭代调整Go?网络中的非均匀相互作用强度，无需在每一步优化都进行模拟，从而显著提升了粗粒化系综与原子参考系综在本质动力学和局部柔性上的一致性。该方法同样适用于优化弹性网络模型，为精确研究蛋白质构象变化提供了高效工具。

　　

引言

蛋白质是生命系统的核心功能分子，其功能常依赖于构象变化以实现不同生物状态间的转换。研究表明，大多数蛋白质动力学可以由少量的集体变量（通过主成分分析(PCA)等技术识别）来有效描述。将高维动力学投影到这些“本质”主成分上并量化低维自由能景观，是理解蛋白质构象动力学和折叠的关键。粗粒化模型通过简化系统表示、牺牲快速自由度的细节，可以显著降低计算成本，从而实现对原子分子动力学在时空尺度上的有力近似。其中，Martini粗粒化力场是模拟生物分子系统最广泛使用的基于物理的力场之一。然而，由于方向性相互作用（如氢键）被平均化，Martini自身无法维持折叠的蛋白质结构。为此，通常将Martini与基于结构的模型，如弹性网络模型(ENM)或Go?模型结合使用。尽管ENM和Go?模型能够保持蛋白质的整体折叠并捕获其柔性，但它们重现构象动力学的能力尚不明确，尤其是在与原子模拟相比时，其在本质子空间中采样到的构象空间有限。

方法

微扰理论

研究采用并扩展了Koyama等人提出的微扰理论框架。对于一个由加和势函数V(q) = Σ_i=1^NL_i(q)描述的分子系统，其正则构型分布为ρ(q) = (1/Z) e^-V(q)。从构成V(q)的N个势能中，选择一个由M个势能L(q) = [L₁(q), ..., L_M(q)]^T组成的子集，并引入微扰系数λ = [λ₁, ..., λ_M]^T，构建一个新的微扰势能：V^λ(q) = V(q) - λ^TL(q) = V(q) - Σ_k=1^Mλ_kL_k(q)。对应的微扰构型分布为ρ^λ(q) = (Z/Z^λ) ρ(q) e^λTL= ρ(q) e^-ψ(λ)+λTL，其中 e^ψ(λ)≡ ?e^λTL?。对于小扰动（∥λ∥ ? 1），可将ψ(λ)在λ=0处进行二阶泰勒展开，得到近似表达式。通过将协方差矩阵C对角化并进行基变换，最终得到一个高效的计算表达式，可以量化一组线性无关的微扰函数所引起的模拟蛋白质构型分布变化。该微扰理论框架为后续的力场优化提供了理论基础，使得无需在每次尝试参数更改时都进行分子动力学模拟，即可分析性地确定网络改进将如何移动构象系综。

系综相似性度量

为了量化测试的原子模型和粗粒化模型在本质子空间上的一致性，研究引入了三个度量指标：均方根内积(RMSIP)、协方差重叠(CO)和切片Wasserstein距离(SWD)。RMSIP衡量由前n个特征向量张成的两个子空间之间的几何相似性。CO与RMSIP类似，但强调沿大方差方向的子空间重叠，同时削弱携带小方差的重叠。与RMSIP或CO等其他结构相似性度量不同，SWD提供了一个真实的距离度量，可以更好地比较采样分布本身。

原子与粗粒化模拟

所有分子动力学模拟均使用GROMACS 2023进行。原子模拟使用Amber99SB*-ILDN力场和TIP3P水模型；粗粒化模拟使用Martini 3力场。Martini拓扑使用martinize2生成。初始弹性网络（κ = 500 kJ/mol/nm²）和Go?网络（ε = 9.4 kJ/mol）使用参考晶体结构构建。对所有系统进行了初始最小化。所有生产模拟均使用蛙跳积分器，原子系统和粗粒化系统的时间步长分别为2 fs和20 fs。使用速度重标度热浴在300 K下维持温度，使用C-rescale气压计在1 bar下维持压力。

结果与讨论

均匀Go?网络在重现本质动力学方面的局限性

研究选取了三个测试系统：T4溶菌酶(T4L)、大肠杆菌核糖结合蛋白(RBP)和大肠杆菌麦芽糖结合蛋白(MBP)。研究产生了5 × 600 ns的原子MD轨迹，以及使用ENM或Go?网络的Martini采样。对前向映射的原子轨迹进行笛卡尔空间主成分分析，揭示了多样的特征运动。将原子和粗粒化轨迹投影到本质子空间显示，均匀的9.4 kJ/mol Go?模型系综明显比原子系综更扩张，且在T4L和MBP中采样了多个最小值。ENM系综虽然部分重叠，但也未能完全复现原子MD模拟在本质子空间内的采样。通过RMSIP、CO和SWD度量进行系统比较发现，虽然两种模型都能合理捕获测试蛋白的本质模式和子空间整体形状，但在采样的构象空间分布上重叠度不佳，特别是Go?模型的SWD值更大。增强Go?网络相互作用强度（ε = 14.1 kJ/mol）或根据原子轨迹调整其平衡距离，仅能带来适度而非最优的改进。这些结果表明，需要优化非均匀的Go?网络以匹配参考系综在本质子空间的行为。

基于微扰的非均匀Go?网络优化

为匹配本质子空间中的参考系综，研究提出并实施了一种基于微扰理论的非均匀Go?网络优化方法。首先，将原子系综映射到Martini粗粒化表示，进行主成分分析并将轨迹投影到所得的本质子空间以获得目标分布ρ_AA(q)。生成一个以均匀Go?网络初始化的粗粒化系综，并投影到同一子空间得到ρ_CG(q)，作为优化的基础。随后构建一个目标函数χ(λ)，该函数衡量应用微扰λ后，预测的粗粒化分布变化与目标分布变化之间的差异。利用微扰理论导出的近似解析表达式，可以在无需运行微扰系综模拟的情况下计算χ(λ)。这使得能够使用快速的粒子群优化(PSO)方法来寻找最小化目标函数的最优微扰λ。然后使用更新后的Go?网络模拟新的粗粒化系综，并迭代重复此过程直至收敛。该方法被命名为PoGo?。

测试表明，该优化算法对所有测试蛋白质都能在约30步内快速收敛，且与初始Go?网络强度的选择无关。优化产生了每种蛋白质独特的Go?网络，部分键增强，部分键减弱。尽管从不同初始强度开始的优化会收敛到不同的参数集（表现为不同的平均相互作用强度），但这些解在目标可观测量（即本质动力学）上是等效的。所有优化方案产生的粗粒化系综与参考原子MD系综在本质动力学上都具有出色的一致性，SWD值接近0.1 nm，比均匀Go?网络报告的值提高了一个数量级以上。

优化的Go?网络重现原子本质动力学并改善局部蛋白质柔性

与未优化的均匀Go?网络相比，优化后的Go?网络产生的系综在视觉上与参考原子系综更相似。对本质子空间的详细分析显示，所有测试蛋白质系统的SWD值显著降低，RMSIP和CO值也有所改善。即使将计算扩展到前五个或十个主成分，相对于未优化的Go?网络仍有改进。作为独立的交叉验证，对骨架珠局部均方根波动(RMSF)的分析表明，优化后的系综与参考原子系综的波动谱一致性显著提高，每个残基的平均无符号偏差大约减少了一半。相比之下，单纯通过均匀增加相互作用强度来改善局部RMSF谱，只能带来与原子本质子空间重叠度的微小提升。这表明，重现局部RMSF并不一定能确保准确的全局动力学，而PoGo?优化方法以更小的总相互作用能扰动，在波动谱和本质子空间上都实现了更好的一致性。

优化的ENM同样可以重现本质动力学

鉴于对ENM中弹簧常数的扰动也是线性的，研究也将优化框架应用于ENM。对于所有三个测试蛋白质，优化仅对原本合理的RMSIP和CO值带来了微小改变，但显著改善了系综重叠和波动谱。这证明了该方法可以同时应用于Martini力场的Go?和ENM组件。

结论

总而言之，研究开发并评估了PoGo?，一种用于Go?网络（此处专门应用于Go?Martini）的基于微扰的优化方法。通过解析地确定网络改进将如何移动构象系综，可以提出单一、信息量最大的更新，避免了在每次尝试参数变化时进行MD模拟的需要。研究在三种不同的蛋白质系统上测试了该方法，在每种情况下都在数十个优化步骤内获得了收敛的Go?网络。由此产生的粗粒化系综在本质子空间中的自由能景观与相应原子参考模拟的结果非常吻合。研究进一步将该框架扩展到优化ENM中的弹簧常数，并取得了类似的成功。此外，研究发现，虽然未明确优化原子波动，但改善沿本质主成分的一致性同时也改善了原子和粗粒化波动谱之间的一致性。简而言之，研究提供了一种完全自动化的方法，用于优化基于Go?Martini的蛋白质模型的本质动力学和波动。