什么是工业酶研究中的计算化学方法？

本文介绍了工业酶研究，特别是指导酶工程研究中使用的主要计算化学方法，包括分子力学力场和分子动力学模拟、量子力学和量子力学/分子力学联合建模、连续静电建模和分子对接等。文中分别从两个方面对这些方法进行了总结：一是方法本身的基本概念、原始计算结果、适用条件和优缺点；二是通过计算获得的宝贵信息，用于指导突变体和突变体库的设计。

酶的工业应用已有百年历史，酶催化因其高效率、高专一性和高选择性以及环境友好性，被广泛应用于食品、农业、医药和化工等不同行业。 ^[1-2]. 由于工业应用环境与酶在自然界的生存环境相差甚远，天然酶的性质和催化功能与其应用环境往往不匹配或不理想。在这种情况下，有必要借助酶工程来改变酶的天然氨基酸序列，以改善其性能 ^[3]. 最常用的酶工程策略是构建突变库进行筛选，即实验室定向进化 ^[4]. 有效定向进化的必要先决条件之一是，接受筛选的突变体库的库容量（即库中包含的突变体数量）相对于库中潜在有益突变体的比例要足够大。突变体库的大小往往受到筛选方法、可用资源和其他客观条件的限制。关键问题是如何提高突变体库中有效突变体的比例。深入了解酶序列、结构和重要特性之间的关系，有助于识别突变热点，限制突变范围，实现高质量的突变体库设计。计算化学方法是获得这种认识的重要手段。研究表明，与随机突变库相比，基于计算设计的蛋白质突变库可将有效突变体的百分比提高几个数量级 ^[5]. 对于一些高难度的酶工程或蛋白质工程课题，计算所带来的实质性改进可能足以决定课题的最终成败，而这已不再局限于效率的提高。事实上，计算化学和计算生物学方法已经成功地从零开始设计出具有天然酶所不具备的催化功能的人工酶。由于本专辑中的其他综述专门讨论氨基酸序列的自动优化设计方法，本文将重点讨论模拟和分析具有给定氨基酸序列的酶的计算方法。当然，研究人员可以使用这些方法分别研究野生型和突变体，然后比较结果。
长期以来，蛋白质（尤其是酶）研究一直是计算化学研究的一个重要前沿领域 ^[6-8]. 主要方法包括基于经典分子力学力场的分子动力学模拟（经典 MD ^[9], 量子力学 ^[10]和 量子力学/分子力学（QM/MM）组合方法 ^[8,11-12], 分子间复合物预测，即分子对接（Docking） ^{[13]，以及} 可极化连续体模型 (PCM)，以量化静电和溶剂效应[ ^14]. ），如泊松-波尔兹曼模型（PB） ^[14], 以及一些基于几何特性的模型。在本文中，我们将从两个方面对这些方法逐一进行概述：首先是方法本身，包括基本原理、原始计算结果、适用条件和（潜在）优缺点等；其次是如何利用这些方法获取工程相关信息，如深入了解催化相关机理、理论预测或解释催化剂的变化。第二，如何利用这些方法获得与工程学相关的信息，如加深对催化相关机制的理解，从理论上预测或解释不同突变体相对于野生型在性质或功能上的变化，从而指导设计用于定向进化的高质量突变库，或设计高质量突变库，或设计突变库，或根据原始计算结果分析提出特定突变位点和突变类型等。

1

基于经典分子力学力场的分子动力学模拟 (MD)

1.1

方法介绍

我们暂不考虑酶催化过程中的化学变化，只考虑酶的构象变化、酶与反应物（或产物）之间由于分子热运动而形成和解离非共价复合物的过程。在这些过程中，分子的电子状态不会改变（例如，共价键不会断裂或产生），分子力学力场模型适用。所谓分子力学力场是一种经验数学函数，它表达了分子体系的势能与几何构型（即组成分子体系的所有原子的空间坐标）的关系（图 1A）。换句话说，如果我们使用 X 表示所有原子的空间坐标，而 ^VMM (X )来表示分子力场势能，势能随着分子从一种构象到另一种构象的变化而变化。 _X1 到另一个 _X2:
^∆VMM = ^VMM ( _X2 ) – ^VMM ( _X1 ).
根据热力学理论，分子中的原子始终处于热运动状态，即.....、 X 是随时间不断变化的；此外，当我们进行实验观测时，样本总是由大量分子组成（单分子实验除外），不同的分子处于不同的构象状态。因此，从动力学角度看，我们需要考虑构象随时间的变化；从热力学角度看，我们需要考虑不同构象分子的概率分布。分子动力学模拟（MD）是研究这两方面性质最直接的模型（图 1B）。在 MD 模拟中，我们从初始构象开始，根据当前构象和势能函数（力是势能函数相对于原子坐标的负导数）计算每个原子在每个时间点上的作用力，对牛顿运动方程进行数值积分，得到下一个时间点的构象，然后重复这一过程，得到构象随时间演变的轨迹。
在这两者之间，可以使用特殊算法来模拟环境因素（如温度、压力等）对分子运动的影响。根据热力学原理，当时间间隔足够长时，同一分子在不同时间点的构象概率分布和不同分子在热力学平衡态的构象概率分布是相同的（即时间平均相当于系统平均）。因此，如果进行 MD 模拟的时间足够长，模拟得到的构象集就可以作为特定热力学平衡态下分子构象分布的样本。基于这一原理，我们可以根据 MD 得到的时间轨迹分析系统在热力学平衡态下的任意可观测特性。
MD 提供了一种强大的计算工具，可以在原子分辨率下全面分析构象变化的动力学变化过程和重要的微观量子热力学分布，这对于阐明酶等复杂生物大分子机器的设计原理和工作机理尤为重要。由于目前的大分子结构分析实验方法只能提供时空平均的静态结构，因此 MD 模拟在相关研究中具有不可替代的作用。在此前提下，MD 工具本身仍处于不断改进和完善的过程中。在方法论上，MD 的主要局限性来自两个方面：一是分子力场模型的准确性；二是有限的模拟时间难以实现构象空间的全面采样。针对第一个问题，近年来分子力场得到了很大的改进，对生物大分子尤其是蛋白质体系构象平衡的热力学描述精度不断提高，成功模拟了多种蛋白质分子从头部折叠过程到天然结构的构象平衡。 ^[15-16]. 在模拟时间方面，由于计算机硬件和软件的发展，现在可以使用常规计算硬件（如研究小组使用的多核服务器）对一般大小的系统（如水溶液中数百个残基的酶分子）进行微秒级的模拟。在这个时间尺度上，可以观察到结构域或环区的打开和关闭等过程。如果有更多的计算资源，还可以直接模拟底物结合/解离等过程。要研究模拟无法达到的时间尺度上的过程（如异构蛋白的大规模功能变化等），可采用增强型采样方法 ^[17], 前提是用户对 MD 理论有更深入的了解。
目前，大多数 MD 模拟应用的时间尺度从纳秒到微秒不等，构象空间的取样大多局限于初始结构附近（对于单链蛋白质，通常是 3-4 Å 范围内均方根位移（RMS）的结构涨落）。因此，有必要使用合理的初始结构作为 MD 的输入，这样模拟结果才有意义。在大多数情况下，实验确定的晶体结构或基于同源蛋白比较建模的结构都被用作 MD 的初始结构。在模拟酶-底物复合物时，通常需要根据空酶的结构或酶与其他分子的复合物结构来模拟复合物的初始结构，方法是使用分子对接或直接用底物替换晶体结构中的其他小分子（如抑制剂）。MD 模拟还需要构建一个分子力场，以描述系统中的所有化学单元。当要模拟的系统包括作为底物的小分子时，通常情况下 MD 软件包中提供的标准分子力场并不涵盖作为底物的小分子。在这种情况下，可使用可自动生成小分子力场的工具软件 ^[18-19]. 在使用自动生成的力场进行长时间 MD 模拟之前，应先对力场文件进行手动检查并用于短时间的模拟试验。

图 1 分子力学力场 (A) 和 分子动力学模拟 (B)
图 1 分子力学力场（A）和分子动力学模拟（B）

1.2

在酶研究中的应用

从 MD 模拟中获得的信息可通过不同方式用于指导酶工程改造 ^[20]. 例如，通过比较室温和高温 MD 模拟，可以预测酶分子中哪些区域可能对环境温度具有最敏感的结构稳定性。在这些区域引入脯氨酸点突变、二硫键等可能会增强酶的耐热性。 ^[21-24]. 提高稳定性的另一种策略是设计出能形成更多表面氢键和盐键的突变体 ^[25-26]. 在实验验证这些突变体之前，可以同时模拟野生型和突变体，从理论上评估突变是否能达到预期效果 ^[27-28]. 除温度外，MD 还可用于分析环境 pH 值、溶剂等变化对蛋白质构象及其稳定性的影响 ^[29-30].
除了稳定性之外，MD 还被用于预测热点残基，这些残基可能会显著影响与底物结合/产物释放相关的构象动力学，从而为设计突变或突变体库提供基础，这些突变或突变体库可以改变底物选择性、反应选择性、产物释放速率等。 ^[31-32]. 利用 MD 研究底物/反应选择性的方法之一是比较具有不同（初始结构）的酶-底物复合物的模拟结果，并预测具有较高亲和力（或较高反应性）的底物或结构状态。计算亲和性（或反应性）的严格定量方法是自由能计算 ^[33-34]. 由于自由能计算的计算密集性，目前的应用大多采用定性方法进行预测：相对亲和性的定性判别可基于小分子-大分子复合物结构的稳定性、分子间相互作用的平均能量等，而反应性的定性判别则基于催化和反应官能团的相对几何构型分布等。 ^[35]. 这种定性判别的结果可作为设计定向进化序列库的基础。此外，MD 模拟还可用于分析底物结合/产物解离孔周围的热点残基 ^[36-37]. 这类应用涉及模拟小分子与蛋白质的解离途径，如果模拟时间尺度不够，可以使用增强采样技术来克服困难。 ^[38-39].

2

量子力学（QM）和组合量子力学/分子力学（QM/MM）模型

2.1

方法介绍

要模拟酶催化过程中的化学步骤，如共价键的生成和断裂、电子转移以及不同电子状态之间的跃迁，就需要量子力学（QM）模型。目前，计算化学中常用的量子力学模型分为几种类型，如 ab initio ( 自始 )、密度泛函理论 (DFT) 和半经验方法 ^[40]. 其中，半经验方法的计算成本最低。然而，它们都是非第一原理方法，计算结果的可靠性在很大程度上取决于特定的系统和问题。ab initio 和 DFT 方法都是第一原理方法，具有通用性。与 ab initio 计算相比，实用的 DFT 模型可能包含更多的经验理论近似，但 DFT 能够以极高的计算效率处理电子相关能。此外，对于许多化学反应问题，最好的 DFT 模型对反应过程中能量变化等关键参数的计算误差已经小到约 1 kcal/mol 左右，其结果足以作为确定特定催化机理或反应途径的化学合理性的依据。
给定一个分子的几何构型，就可以用 QM 计算出其能量。这被称为单点计算（即只处理几何构型空间中的一点）。QM 模型更常用于分子几何构型的优化，即从初始构型连续变化后找到局部稳定的结构（能量低于相邻结构），或找到连接反应物和生成物的最低能量路径，以及路径上的过渡态。这些计算的计算量很大，因为需要考虑和比较不同的几何构型，通常需要进行数万至数千次单点计算。节省计算量的常见策略是，首先使用精度有限的高效质量管理模型优化大范围的反应路径搜索，然后在搜索到的最低能量构型/路径附近使用精度更高的模型完成构型优化，或者进行单点计算。
目前，将第一原理 QM 方法应用于整个酶分子的计算量很大，基本上仅限于单点计算，仍然缺乏实用性。QM/MM 模型（图 2）通常用于大分子 ^[11]. 在这一模型中，分子系统至少分为两部分：直接参与化学反应的部分用 QM 模型处理，其余部分用分子力学（MM）模型处理。处理 QM-MM 边界和相互作用有几种不同的策略 ^[41]. 在第一原理的 QM/MM 模型中，QM 计算的成本远高于 MM。因此，构象优化方法主要用于 QM 区域，以预测或模拟其几何形状，而分子动力学模拟可用于 MM 部分的采样。 ^[42]. 这意味着计算结果可能对系统 QM 区域的初始结构更为敏感。在这种情况下，需要使用不同的初始结构模型进行计算，以获得可靠的结果。如果半经验方法 ^{[43] 或} 经验价键理论 ^[44-45] 用于质量管理部分 ^, 通过更长时间的 QM/MM MD 采样，有可能更充分地探索构象空间并减少初始结构的影响。

图 2 量子力学 (QM)/Molecular 机械 (MM) 模型
图 2 量子力学（QM）/分子力学（MM）模型。

2.2

方法的应用

^两者 QM 模型 ^[10]和 QM/MM 模型 ^[41] 已被广泛用于酶催化反应化学机制的理论预测和检验。其结果可以帮助我们辨别哪些关键残基参与了化学反应过程，找到反应的限速步骤，建立反应中间体和过渡态的结构模型，分析它们如何与酶环境相互作用等。与 QM 簇模型相比，QM/MM 模型能更真实地模拟化学反应中心所处的酶环境。QM/MM 已被广泛用于理论预测/检验酶催化的化学机制，以及分析和预测环境氨基酸残基对催化过程可能产生的影响。 ^[46]. 原则上，这些结果可用于指导定向进化突变库的设计，目的是提高催化活性和改变特异性或选择性。更具挑战性的研究是根据QM或QM/MM预测的过渡态结构模型，从零开始设计新的活性中心，从而获得全新的人工酶。 ^[47].

3

静电连续介质建模

3.1

该方法的原理

酶催化作用几乎总是在特定的溶液环境中完成的。溶剂效应对酶的特性有着至关重要的影响。用于计算化学处理的溶剂效应的模型可分为两类：显式溶剂模型，例如分子力学力场或量子力学模型，其中每个溶剂分子及其中的每个原子都明确包含在模型中；隐式溶剂或连续介质模型 ^{[48]，其中} 该模型不包括溶剂分子和原子，而是用所谓的 "溶剂平均场 "来处理溶剂效应。显式溶剂模型的优点是能够以完全一致的方式处理溶质和溶剂，真实地模拟特定的相互作用，如溶质-溶剂氢键、盐键等。缺点是溶剂分子数量大，计算量大。此外，溶剂随机上升对系统总能量的影响很大，必须进行长时间的模拟样本平均才能消除上升的影响。隐藏溶剂模型描绘了溶剂的平均效应，并对溶剂的热力学涨落进行了平均。
为简化处理，我们通常在隐藏溶剂模型中将非极性溶剂效应（疏水效应）与极性溶剂效应分开。经验表明，非极性溶质的溶解自由能与其溶剂可及表面积（SASA）成正比。因此，SASA 溶解模型通常用于此部分。该模型的参数包括计算 SASA 所需的原子半径、溶剂分子半径（水分子为 1.4 Å）以及溶解自由能与 SASA 成比例的比例常数。这些参数一般通过拟合小分子溶解自由能的实验值来确定。
考虑极性溶剂效应的最常用模型将溶剂占据的区域视为具有特定介电常数（水的介电常数为 78.4）的连续介质，而将溶质区域视为由具有低介电常数（通常为 2-8）或真空（介电常数为 1）的介质占据（图 3A）。连续介质被溶质区电荷分布产生的静电场极化，由此产生的极化电荷分布反过来又在溶质区产生作用于溶质电荷的静电场。极化电荷产生的电场称为反应场。因此，静电连续介质模型也被称为反应场模型。在溶剂区没有自由离子的连续介质模型中，空间静电势与空间电荷分布之间的关系满足泊松方程。对于含有自由离子的溶液环境，离子的空间分布会受到空间静电势的影响。考虑到这一因素，空间静电势与空间电荷分布之间的关系满足泊松-波尔兹曼方程（Poisson-Boltzmann equation，简称 PB 方程）。PB 方程是关于三维空间中静电势分布与电荷和电介质分布之间关系的偏微分方程，可以用数值方法求解。求解酶等大分子系统的 PB 方程最常用的数值方法是有限差分法（FD），统称为 FDPB 模型（图 3B）。 ^[14]. 利用 FDPB 可以根据溶质的空间电荷分布计算三维静电势，进而计算静电自由能等其他性质。在小分子体系的质量管理计算中，反应场通常等同于分子表面电荷分布产生的电场，相应的模型称为可极化连续介质（PCM）模型。
最常用的极性溶剂效应模型是将溶剂占据的区域视为具有特定介电常数（水为 78.4）的连续介质，而将溶质区域视为具有低介电常数（通常为 2-8）或真空（介电常数为 1）的介质（图 3A）。溶质区域的电荷分布所产生的静电场会使连续介质极化，由此产生的极化电荷分布反过来又会在溶质区域产生作用于溶质电荷的静电场。极化电荷产生的电场称为反应场。因此，静电连续介质模型也被称为反应场模型。在溶剂区没有自由离子的连续介质模型中，空间静电势与空间电荷分布之间的关系满足泊松方程。对于含有自由离子的溶液环境，离子的空间分布会受到空间静电势的影响。考虑到这一因素，空间静电势与空间电荷分布之间的关系满足泊松-波尔兹曼方程（Poisson-Boltzmann equation，简称 PB 方程）。PB 方程是关于三维空间中静电势分布与电荷和电介质分布之间关系的偏微分方程，可以用数值方法求解。求解酶等大分子系统的 PB 方程最常用的数值方法是有限差分法（FD），统称为 FDPB 模型（图 3B）。 ^[14]. 利用 FDPB 可以根据溶质的空间电荷分布计算三维静电势，进而计算静电自由能等其他性质。在小分子系统的质量管理计算中，反应场通常等同于分子表面电荷分布产生的电场，相应的模型称为可极化连续介质（PCM）模型。

图 3 静电连续模型 (A) 和有限差分泊松-波尔兹曼 （发展政策局） 方法 (B)
图 3 静电连续体模型 (A) 和有限差分泊松-波尔兹曼 (FDPB) 方法 (B)。

3.2

方法的应用

连续模型的重要应用之一是研究酶分子中带电氨基酸侧链基团的质子化状态。软件 PROPKA 通过求解 PB 方程来计算不同质子化状态下的静电自由能，从而预测每个可解离基团的 pKa ^[50]. 酶分子的表面静电位分布是影响酶底物选择性的一个重要因素。鉴于酶分子的空间结构和质子化状态，可以通过 FDPB 方法计算出酶分子的表面静电位分布，该方法还能预测氨基酸残基突变或环境 pH 值变化、离子浓度变化等对表面静电位的影响。 ^[14].
使用 QM 簇模型研究酶催化的化学步骤时，通常需要使用 PCM 模型模拟环境对反应区的静电影响。如果反应过程涉及电荷分布的重大变化，不使用连续介质的真空 QM 计算结果就不合理，甚至可能导致错误的定性结论。在 QM/MM 模型中，反应中心一般被可溶性分子包围，这些分子按 MM 方式处理，一般无需考虑连续介质反应场。但是，如果系统的净电荷在反应前后发生变化（如氧化还原电位计算），则很可能需要考虑系统边界外溶液环境对反应自由能的贡献，在这种情况下，可以使用连续介质模型对 QM/MM 结果进行修正。
作为一种兼顾效率和准确性的方法，MM/PBSA 可用于分析蛋白质-蛋白质和蛋白质-小分子复合物的亲和力 ^[49]. 为了消除误差，传统的方法是对复合物进行显式溶剂分子动力学模拟，得到一系列构象；针对每个复合物构象，分别计算复合物整体和组成复合物的每个单体的 MM/PBSA 能量；利用整个构象的整体 MM/PBSA 能量与单体 MM/PBSA 能量之差的平均值来近似计算结合自由能。这种方法可用于分析影响底物亲和力的热点残基，也可用于预测突变体的底物选择性变化。

4

其他方法

4.1

分子对接

对接是指根据单体结构预测复合物结构（和亲和力）的计算过程。小分子-蛋白质对接是基于结构的药物虚拟筛选的核心工具，目前已开发出多种算法 ^[13]. 这些算法和模型也可用于底物-酶复合物的对接。虚拟药物筛选需要考虑大量不同的小分子，出于计算效率的考虑，分子对接计算往往不考虑受体（或仅侧链）的结构变化。与虚拟筛选不同，底物-酶对接研究通常只考虑一种或几种不同的底物，原则上可以更充分地考虑酶的结构变化。要做到这一点，最直接的方法是通过构象取样方法（如 MD）获得不同的酶结构，并分别与底物对接。在底物-酶对接中，通常还可以利用底物和催化功能基团的相对空间排列来筛选/评估对接结果。

4.2

根据几何形状预测小分子孔隙

大量实验研究发现，一些远离活性中心的突变会对酶的催化性能产生很大影响。其中一些位点可能通过改变底物结合/产物释放孔道而发挥作用，孔道大小、孔道周围残基的理化性质等都会改变底物/产物通过率，影响底物选择性。孔道预测方法可用于寻找相关的热残基，并为定向进化文库的设计提供依据。有几种基于几何结构的方法可用于预测蛋白质表面凹坑、内部空腔、连接不同区域的孔等 ^[51-53]. 这些方法使用静态空间结构作为输入，大多采用几何和图论方法，以较高的计算效率实现预测。

4.3

主动式中心比较方法

目前，蛋白质三维结构数据库（PDB）已积累了大量不同结构类型和家族的酶的三维结构数据。如果我们比较不同的酶，就会发现其中一些酶的活性中心具有高度相似性（典型的例子是丝氨酸蛋白酶共有的催化三元活性中心），尽管整体结构序列并不相似。活性中心结构比较法 ^{[54-55] 可以} 可用于自动检索与当前酶的活性中心相似的其他酶的活性中心。在三维空间中将多个相似的活性中心堆叠在一起，分析不同活性中心之间的异同，可为突变位点选择提供有价值的信息。

5

摘要

为了表述清晰，我们对上述方法进行了分类介绍。实际上，不同类型的方法并不相互排斥。它们可以通过多种方式结合使用，以更好地回答我们感兴趣的问题。例如，在酶-底物复合物模拟中，分子对接可用于获得模拟的初始构象；通过 MD 模拟获得的构象集可用于孔预测分析、分子对接、QM/MM 模拟等。如我们前面提到的MM/PBSA方法，是MD和连续介质建模的结合等等。
利用计算化学方法研究蛋白质等生物大分子系统已有 40 多年的历史。这些方法在自身不断发展的同时，在工业酶制剂研究中的应用也越来越广泛。中国在计算化学和工业酶工程领域的研究队伍不断扩大，研究能力迅速提高。随着计算化学与酶工程学科交叉结合的日益紧密，计算化学在酶工程领域的应用也将不断扩大和深化。蛋白质工程、定向进化等技术对工业酶研究产生了巨大影响。未来计算方法的发展，尤其是新型酶设计方法的突破，有望为合成生物学时代的工业酶研究带来新的技术突破。

立即联系我们！

如果您需要价格，请在下表中填写您的联系信息，我们通常会在 24 小时内与您联系。您也可以给我发电子邮件 info@longchangchemical.com 请在工作时间（UTC+8 周一至周六，上午 8:30 至下午 6:00）或使用网站即时聊天工具获得及时回复。

复合葡萄糖淀粉酶	9032-08-0
普鲁兰酶	9075-68-7
木聚糖酶	37278-89-0
纤维素酶	9012-54-8
柚皮苷酶	9068-31-9
β-淀粉酶	9000-91-3
葡萄糖氧化酶	9001-37-0
α-淀粉酶	9000-90-2
果胶酶	9032-75-1
过氧化物酶	9003-99-0
脂肪酶	9001-62-1
过氧化氢酶	9001-05-2
TANNASE	9025-71-2
弹性蛋白酶	39445-21-1
尿素酶	9002-13-5
DEXTRANASE	9025-70-1
L 乳酸脱氢酶	9001-60-9
苹果酸脱氢酶	9001-64-3
胆固醇氧化酶	9028-76-6