X射线晶体学与QM/MM模拟联手：揭示嘧啶从头合成关键酶OPRTase的催化反应机制 - 附录

本文是主文档的技术附录，详细介绍计算化学方法、模拟参数设置和Supporting Information中的补充结果。

计算方法详解

分子动力学模拟设置

体系构建

起始结构采用PDB ID：6GV9（OPRTase与OA和$\ce{SO4^{2-}}$复合物，1.25 Å分辨率），使用pdb4amber工具去除水分子和硫酸根离子。
质子化状态设置为：所有组氨酸质子化，Asp和Glu去质子化，Lys和Arg质子化。
力场选择方面，蛋白质使用AMBER ff14SB力场，小分子（OA和PRib-PP）使用GAFF力场，电荷由RESP方法在HF/6-31G*水平计算得到。
溶剂化采用TIP3P水模型，八面体盒子，边界距离蛋白质至少10 Å，并加入$\ce{Cl^-}$离子中和体系总电荷。

MD模拟流程

能量最小化：分两步进行，第一步为5000步最速下降加5000步共轭梯度，蛋白质重原子施加$50\,\mathrm{kcal\cdot mol^{-1}\cdot \mathrm{Å}^{-2}}$的约束；第二步为2500步最速下降加2500步共轭梯度，无约束优化
平衡阶段：包括三个步骤，首先NVT升温从0 K至300 K，持续50 ps，施加$10\,\mathrm{kcal\cdot mol^{-1}\cdot \mathrm{Å}^{-2}}$的弱约束；然后NPT平衡在300 K和1 atm下进行500 ps，保持弱约束；最后为NPT系综下的500 ps模拟，无约束
生产模拟：采用Langevin恒温器维持300 K（碰撞频率$2\,\mathrm{ps^{-1}}$），Berendsen控压器维持1 atm（耦合时间2 ps），静电采用PME方法（截断距离10 Å），氢键由SHAKE算法约束，时间步长2 fs，总模拟时长100 ns
轨迹分析：使用CPPTRAJ工具分析关键距离演化，监测活性位点几何构型变化，分析催化环（残基94-110）的柔性

模拟体系概况

整个模拟体系包含约45000个原子（含蛋白质、底物、溶剂和离子），蛋白质部分为213个氨基酸残基，溶剂包含约14000个TIP3P水分子，模拟盒子尺寸约为70 Å × 70 Å × 70 Å。

QM/MM计算细节

QM/MM分区

QM区域（44原子）：包括OA底物完整分子（15个原子）、PRib-PP中的磷酸基团$\ce{PO3^{2-}}$部分（5个原子）、1个显式水分子Wat318（3个原子）
关键催化残基侧链（21原子）：Lys73的NZ-HZ3共3个原子、Asp125的CG-OD1-OD2共4个原子、Arg99*的完整侧链共11个原子、Lys103的NZ-HZ3共3个原子
MM区域：包括蛋白质其余部分、PRib-PP的核糖环部分、所有溶剂分子和抗衡离子
边界处理：采用Link Atom方法处理共价键断裂，在Cα-Cβ键处切割侧链，总共使用5个Link Atom

QM/MM方法学

QM计算采用B3LYP/6-31G*密度泛函理论方法，MM计算使用AMBER ff14SB和GAFF力场，耦合方式为电子嵌入（Electrostatic Embedding），软件包为Gaussian 09结合AMBER 16。

自适应弦方法（Adaptive String Method）

什么是自适应弦方法？

想象你要从山的一边走到另一边，有无数条路可以选择。最小自由能路径（MFEP）就像是找到一条“最容易走”的路——不是最短的直线距离，而是综合考虑了爬坡难度、能量消耗等因素后，最省力、最可能被自然选择的路径。

在分子世界中，化学反应是分子体系在复杂的高维自由能面（能量地形）上“滑动”的过程。自适应弦方法就是用一根“弦”（由多个节点组成的离散化路径）来描绘这条最优路径。“自适应”是指这根弦会自动调整形状，逐步逼近真正的最小自由能路径，而不需要预先定义反应坐标。

为什么需要它？

传统的反应路径优化方法（如NEB, Nudged Elastic Band）通常需要预先定义反应坐标，且在真空或简化模型中进行。但酶促反应的特点是：

高维复杂性：涉及多个原子的协同运动（质子转移、成键断键、蛋白质构象变化）
环境效应：蛋白质和溶剂的动力学涨落显著影响反应路径
多通道性：可能存在多条竞争性反应通道，需要探索哪条能垒最低

自适应弦方法通过在显式溶剂和蛋白质环境中进行MD采样，能够：

自动识别最优反应坐标
考虑环境的动力学效应
准确计算包含熵贡献的自由能

算法原理：迭代优化循环

根据原文描述，自适应弦方法的实现流程如下：

初始化阶段

定义集合变量（CVs）：选择能充分描述反应进程的几何参数
- 成键/断键距离：如 d(N1-C1)、d(C1-O1) 等
- C1原子杂化坐标：描述C1从sp³到sp²再到sp³的转变
- 对于OPRTase的不同机制，使用4-7个CVs组合
- 关键区别：CVs是多维空间的坐标轴（如4-7维），而反应坐标是从反应物到产物的特定路径。传统方法需要预先指定用哪个CV或CV组合作为反应坐标，而自适应弦方法允许在CV空间中自动搜索最优路径
构建初始路径：沿着CVs定义反应物到产物的初始猜测路径
- 选择80个等距节点离散化路径
- 势能profile预探索：在正式的自适应弦方法迭代前，先使用PM6/MM水平对沿CVs的反应路径进行粗略扫描，计算势能剖面（见SI Figure S8）。这比简单的几何线性插值更合理，因为已考虑了体系的能量信息，避免初始路径经过高能区域
- 每个节点代表反应路径上的一个中间构象

迭代优化循环

对于每一轮迭代，执行以下步骤：

独立MD采样（针对每个节点）：对80个节点同时启动独立的MD模拟（时间步1 fs，最多250 ps）。每个节点的模拟受CVs约束，保持在路径上的指定位置。采样该节点附近的构象空间，积累统计力学数据
计算自由能梯度：从每个节点的MD轨迹中提取自由能的近似斜率。这个梯度指示了体系倾向于朝哪个方向演化
节点移动 + 重参数化：每个节点向更低自由能的方向移动，但只能垂直于路径的方向移动（通过拉格朗日乘子去除切向分量），移动后立即重新调整所有节点位置使其等间距。这一步是自适应弦方法的核心：既让路径向MFEP演化（垂直方向往低处走），又防止节点聚集（保持等距约束）
副本交换（增强采样）：每50步尝试一次相邻节点之间的构象交换。帮助克服局部能量势垒，加速收敛
收敛判断：通过测量节点前后位置的平均距离监控收敛。当变化足够小时，弦达到MFEP
重复迭代：重复上述步骤，直到弦最终收敛到MFEP

关于“重参数化”的补充说明

什么是“参数”？弧长坐标s如何计算？

重参数化中的“参数”是指弦上每个节点的弧长坐标 s：

弦是一条连接反应物和产物的曲线，每个节点是弦上的一个点
每个节点$i$对应一个弧长参数 $s_i$，$s_0 = 0$（反应物），$s_{80} = 1$（产物）
s的计算方法：
- 在多维CV空间中，节点$i$和节点$i+1$之间的欧氏距离为：
\[\Delta l_i = \sqrt{\sum_{k=1}^{N_{\text{CV}}} (\text{CV}_k^{i+1} - \text{CV}_k^i)^2}\]
- 从反应物到节点$i$的累积弧长：$L_i = \sum_{j=0}^{i-1} \Delta l_j$
- 归一化的弧长坐标：
\[s_i = \frac{L_i}{L_{\text{total}}}, \quad L_{\text{total}} = \sum_{j=0}^{79} \Delta l_j\]
- 这样确保 $s$ 在0到1之间均匀分布，将多维CV空间投影到一维反应进程坐标
重参数化就是重新调整这些节点在弦上的位置，使得相邻节点之间的弧长间距 $\Delta s = 1/80$ 相等

为什么必须“保持等距”？

防止节点聚集
- 若不约束，节点会自发向低能区聚集（如反应物和产物附近）
- 导致过渡态附近缺乏采样点，无法准确描述能量变化最剧烈的区域
确保算法收敛
- 等距约束是弦方法收敛到正确MFEP的数学必要条件
- 通过拉格朗日乘子去除自由能梯度的切向分量，只保留垂直于路径的分量
实现方式
- 通过三次样条插值重新定义弦的参数化方程
- 在新的等距参数点上重新采样节点位置

需要注意的局限：

节点按弧长（而非能量）均匀分布
能量最高的节点不一定恰好对应过渡态的几何构型
需要额外验证过渡态位置（如通过频率分析）

自由能profile计算

在收敛的路径上进行伞形采样：

在每个节点设置谐振势约束（力常数$10\,\mathrm{kcal\cdot mol^{-1}\cdot \mathrm{Å}^{-2}}$）
每个窗口MD模拟500 ps
使用WHAM（加权直方图分析方法）重构完整的势能均值力（PMF）
得到沿反应进程s的自由能曲线$G(s)$

本研究的具体实现细节

参数	数值
节点数量	80个等距节点
QM方法	PM6（初步探索）+ M06-2X/6-311+G(2df,pd)（精细计算）
MM力场	AMBER ff14SB（蛋白质）+ TIP3P（水）
每节点采样时间	最多250 ps
副本交换频率	每50步尝试一次
收敛标准	梯度 < $0.05\,\mathrm{kcal\cdot mol^{-1}\cdot \AA^{-1}}$
伞形采样窗口	每窗口500 ps，力常数$10\,\mathrm{kcal\cdot mol^{-1}\cdot \mathrm{Å}^{-2}}$

关键理解：初始路径不需要完美。自适应弦方法会在迭代中自动“修正”它，让弦沿着真实的最小自由能路径滑动。这就是“自适应”的含义——算法主动寻找最优路径，而不是死守初始猜测。

自由能微扰（FEP）计算

热力学循环

为了计算OA两种互变异构形式（lactam vs lactim）在酶中的相对稳定性，采用了 Scheme 3 中的热力学循环：

\[\begin{aligned} &\text{OA}_{\text{lactam}}^{\text{gas}} \xrightarrow{\Delta G_{\text{gas}}} \text{OA}_{\text{lactim}}^{\text{gas}}\\ &\quad\downarrow \Delta G_{\text{Amide,p}} \qquad\downarrow \Delta G_{\text{Imidic,p}}\\ &\text{OA}_{\text{lactam}}^{\text{protein}} \xrightarrow{\Delta G_{\text{Protein}}} \text{OA}_{\text{lactim}}^{\text{protein}} \end{aligned}\]

因此：

\[\Delta G_{\text{Protein}} = \Delta G_{\text{gas}} + (\Delta G_{\text{Imidic,p}} - \Delta G_{\text{Amide,p}})\]

模拟参数

λ窗口设置：耦合参数λ从0到1划分为21个窗口（间隔0.05），每个λ窗口包含200 ps的平衡阶段和1 ns的生产阶段，温度维持在300 K，总模拟时间为2 × 21 × 1 ns = 42 ns
Soft-core势函数：参数设置为α = 0.5和σ = 3.0 Å
自由能计算方法：自由能变化ΔG采用Bennett Acceptance Ratio（BAR）方法计算，统计不确定度通过Bootstrap方法进行估计（1000次重采样）

计算结果

气相能量差：使用M06-2X/6-311+G(2df,pd)优化与频率计算得到$\Delta G_{\text{gas}} = 27.5\,\mathrm{kcal\cdot mol^{-1}}$，酰胺形式在真空中最稳定
蛋白质环境相互作用：基于AMBER 16/pmemd.cuda的FEP（21个λ窗口，每窗口1 ns生产段）给出$\Delta G_{\text{Imidic,p}} - \Delta G_{\text{Amide,p}} = -7.6 \pm 0.1\,\mathrm{kcal\cdot mol^{-1}}$，说明活性位点更偏好亚氨酸形式
综合差值：$\Delta G_{\text{Protein}} = 27.5 - 7.6 = 19.9\,\mathrm{kcal\cdot mol^{-1}}$，即便蛋白质提供部分稳定，也不足以翻转互变异构体的能量排序，酰胺形式仍是酶中最稳定并充当反应起点的状态

Q&A

Q1：为什么之前的计算研究未能准确描述OPRTase的反应机制？
A1：以往的计算优化是在真空或简化模型中进行的，忽略了多个关键因素：
- 蛋白质环境的静电效应：保守残基（Lys73、Asp125、Arg99*、Lys103）和$\ce{Mg^{2+}}$对过渡态的静电稳定至关重要
- 蛋白质的灵活性：催化环的开-闭动力学对催化周期至关重要
- 底物互变异构形式的相对稳定性：需要FEP计算才能准确评估酶中酰胺和亚氨酸形式的能量差
- 水分子的作用：活性位点中的水分子作为质子中继，无法在真空计算中体现
- 本研究通过结合高分辨率晶体结构、长时间MD模拟和QM/MM自由能计算，首次全面考虑了这些因素
Q2：自适应弦方法相比传统的反应路径优化有什么优势？
A2：自适应弦方法具有五大优势：
- 自动寻找最小自由能路径（MFEP）：虽需预先选择CVs（如键长、键角），但无需预先指定哪个CV或CV组合是反应坐标，算法在多维CV空间中自动搜索最优路径并投影到一维弧长坐标s
- 考虑熵效应：沿路径进行MD采样，自然包含构象熵
- 路径集合变量（s坐标）：将多维反应空间投影到一维，简化PMF计算
- 副本交换：增强采样效率，加速收敛
- 适用于复杂机制：可处理多步骤、多中间体的复杂反应
- 对于OPRTase这种涉及质子转移、亲核攻击和键断裂的复杂机制，传统方法（如NEB）难以有效处理，而自适应弦方法提供了系统性的解决方案
Q3：为什么水分子作为质子中继比直接质子转移能垒低得多？
A3：能垒差异源于四方面原因：
- 几何约束：N1（OA）到O2A（PRPP）的直接距离较远（约4-5 Å），直接质子转移需要大幅构象重排
- 电荷分离：直接转移产生N1⁻和O2AH的电荷分离态，在低介电环境（蛋白质内部）中能量代价高
- 水分子的双重作用：它作为质子受体和供体减小每步质子转移的距离（约3 Å），形成的$\ce{H3O+}$中间体虽不稳定但寿命足够短，迅速将质子传递给O2A
- 蛋白质环境预组织：MD模拟显示该水分子已预先定位在N1和O2A之间，形成稳定的氢键网络
- 水介导机制利用了格罗特斯机制（Grotthuss mechanism）的优势，通过质子接力显著降低能垒
Q4：如何利用本研究的过渡态信息设计OPRTase抑制剂？
A4：基于过渡态结构的抑制剂设计可采用四种策略：
- 过渡态类似物设计（TSA）：模拟TS几何和电荷分布的小分子
  - C1原子引入部分正电荷或氧碳正离子特征（如用$\ce{CH2+}$或缺电子碳替代）
  - N1-C1键使用部分形成的键长度（约2.3 Å，可用柔性连接模拟）
  - 焦磷酸部分保留负电荷中心以利用Arg99*、Lys103、Lys73的静电相互作用
- 保留关键相互作用
  - 保持与Asp125（通过核糖O2羟基）的氢键
  - 保持与$\ce{Mg^{2+}}$的配位相互作用
  - 保持与催化环残基（Arg99*、Lys103）的多重静电相互作用
- 水分子位点填充：设计能占据关键水分子位置的功能基团，阻断质子转移
- 双底物类似物设计：连接OA和PRPP的结构特征，形成双底物TSA，利用两个底物结合位点的协同效应
- 文献中已报道的一些OPRT抑制剂（如硒代芳香化合物、TSA）可根据本研究的TS结构信息进一步优化
Q5：催化环的开-闭动力学如何影响催化效率和反应选择性？
A5：催化环动力学产生六重影响：
- 底物识别：开放构象允许PRPP进入，只有PRPP结合后催化环才倾向闭合，提供诱导契合机制
- 活性位点隔离：闭合后封闭活性位点排除大部分溶剂水分子，降低介电常数，有利于静电相互作用增强（Lys、Arg与底物）和稳定过渡态电荷分布
- 保持关键水分子：尽管排除大部分水，但闭合时保留参与质子转移的关键水分子
- 防止副反应：封闭环境防止PRPP与其他亲核体（如溶剂水或其他残基）发生非生产性反应
- 产物释放控制：反应后催化环重新打开允许产物释放，Lys103与α-磷酸的相互作用可能帮助引导焦磷酸离去
- 交替位点催化：一个亚基的催化环闭合催化反应时，另一个亚基的环打开释放产物，实现高效的交替催化
- 催化环因此不仅是“盖子”，更是动态调控催化周期各阶段的开关

Supporting Information补充结果

关键距离演化分析

Supporting Information的Figures S3-S7展示了MD模拟过程中活性位点关键距离的时间演化。

图S3：Lys73与OA的相互作用

监测参数为d(NZ(Lys73)-O4(OA))，平均距离为2.8 ± 0.2 Å。该距离在整个模拟中保持稳定，支持Lys73作为质子供体的角色。

图S4：Asp125与PRib-PP的相互作用

监测参数为d(OD1(Asp125)-C1’(PRib-PP))，平均距离为3.2 ± 0.3 Å。距离变化较大，反映催化环的柔性。

图S5：Arg99*与焦磷酸基团的相互作用

监测参数为d(NH1(Arg99*)-Oα(PPi))，平均距离为2.7 ± 0.1 Å。形成稳定的氢键网络，稳定离去基团。

图S6：Lys103与磷酸基团的相互作用

监测参数为d(NZ(Lys103)-Oα(PRib-PP))，平均距离为2.9 ± 0.2 Å。持续的静电相互作用活化磷酸基团。

图S7：水分子Wat318的氢键网络

监测参数包括d(O(Wat318)-O4(OA))为2.8 ± 0.2 Å，以及d(O(Wat318)-OD2(Asp125))为2.7 ± 0.1 Å。水分子稳定地桥接OA和Asp125，支持水介导质子转移机制。

过渡态结构详细分析

图S8：三种机制的过渡态几何构型

该图展示了机制I、II、III在各自过渡态（TS1和TS2）的关键几何参数。

机制I（协同机制）：机制I的TS1几何特征为d(C1’-N1) = 2.1 Å（部分成键）、d(Pα-O) = 2.0 Å（部分断键）、d(N1-H) = 1.3 Å（质子转移进行中）、∠(C1’-N1-C2) = 112°（从平面向四面体过渡），能垒为$16.7\,\mathrm{kcal\cdot mol^{-1}}$。过渡态特征为高度协同，所有化学事件几乎同步发生。

机制II（分步机制，先成键）：机制II的TS1几何（成键步骤）为d(C1’-N1) = 1.9 Å（接近完全成键）、d(Pα-O) = 1.7 Å（尚未断键）、d(N1-H) = 1.1 Å（质子转移完成）。中间体几何为五配位磷原子，不稳定，自由能比反应物高$18.3\,\mathrm{kcal\cdot mol^{-1}}$。TS2几何（断键步骤）的d(Pα-O) = 2.2 Å（断键进行中），总能垒为$21.5\,\mathrm{kcal\cdot mol^{-1}}$（TS2相对反应物）。

机制III（分步机制，先断键）：机制III的TS1几何（断键步骤）为d(Pα-O) = 2.3 Å（接近完全断键）、d(C1’-N1) = 3.5 Å（尚未成键）。中间体为碳正离子（oxocarbenium ion），C1’的电正性极高，由Asp125和周围残基稳定，自由能为+$28.7\,\mathrm{kcal\cdot mol^{-1}}$（相对反应物）。TS2几何（成键步骤）的d(C1’-N1) = 2.0 Å（成键进行中），总能垒为$30.2\,\mathrm{kcal\cdot mol^{-1}}$（过高，不可行）。

三种机制的详细比较

Table S1：机制I、II、III的关键参数对比 | 参数 | 机制I | 机制II | 机制III | |——|——-|——–|———| | 反应路径类型 | 协同 | 分步（先成键） | 分步（先断键） | | TS1能垒 ($\mathrm{kcal\cdot mol^{-1}}$) | 16.7 | 18.3 | 28.7 | | TS2能垒 ($\mathrm{kcal\cdot mol^{-1}}$) | - | 21.5 | 30.2 | | 中间体类型 | 无 | 五配位磷 | 碳正离子 | | 中间体自由能 ($\mathrm{kcal\cdot mol^{-1}}$) | - | +18.3 | +28.7 | | 关键质子供体 | Lys73 | Lys73 | Lys73 | | 质子转移时机 | 与成键同步 | 成键前 | 断键后 | | 实验$k_{\text{cat}}$对应能垒 ($\mathrm{kcal\cdot mol^{-1}}$) | 15.5 | 15.5 | 15.5 | | 计算误差 ($\mathrm{kcal\cdot mol^{-1}}$) | +1.2 | +6.0 | +14.7 | | 机制可行性 | ✓ 最优 | ✗ 能垒偏高 | ✗ 能垒过高 |

结论:

机制I（协同机制）与实验数据吻合最好，计算能垒（$16.7\,\mathrm{kcal\cdot mol^{-1}}$）接近实验值（$15.5\,\mathrm{kcal\cdot mol^{-1}}$）
机制II和III的能垒显著偏高，与实验观测到的高效催化不符
协同机制避免了形成高能中间体，降低了反应能垒

计算资源与软件

使用的主要软件包

AMBER 16：MD模拟和FEP计算
Gaussian 09：QM/MM计算
CPPTRAJ：轨迹分析
VMD 1.9.3：结构可视化
PyMOL 2.0：作图和结构分析
WHAM：伞形采样数据分析

计算资源配置

MD模拟使用NVIDIA Tesla V100 GPU加速，QM/MM计算使用48核Intel Xeon处理器，总计算时间约50000 CPU小时。

本附录详细介绍了OPRTase反应机制研究中使用的计算化学方法和补充结果，这些技术细节对于理解主文档的结论、评估研究质量以及为类似研究提供方法学参考具有重要价值。

Mendelevium

Contact