1 引言
工业酶催化过程中，高温操作（60-100 °C）可以大幅降低体系黏度、抑制微生物污染和提升传质效率与反应速率，是生物炼制或纤维素降解过程的核心需求
[1,2]。但多数来源于中温微生物的工业酶在此温度范围内面临构象失稳和活性衰减问题，导致酶制剂补加频繁，生产成本上涨。例如糖苷水解酶家族1（glycoside hydrolase family 1，GH1）的β-葡萄糖苷酶是纤维素水解途径的关键限速酶，负责催化纤维素水解途径的最终步骤，但是该酶的热稳定性不足以达到第二代生物燃料经济性水平
[3,4]。
长期以来，蛋白质热稳定性改造主要依赖随机诱变和定向进化等实验策略
[5,6]。这类方法通过构建突变体文库并结合高通量筛选获得性能改良的酶，在蛋白质工程领域取得了显著成功。但是这种方法往往需要大量的实验研究，耗时、昂贵，且对高通量筛选依赖性强，难以高效、低成本地系统解析热稳定性提升的分子机制
[7]。此外，由于突变是随机造成的，因此定向进化的迭代次数较多，限制了其在复杂酶体系中的应用效率。这些局限性促使研究者寻求更具理性设计能力的计算模拟方法。
计算模拟方法的兴起，使酶的热稳定性能够合理地进行设计。分子动力学模拟（molecular dynamics，MD）可以通过原子分辨率捕获高温诱导构象涨落和氢键网络破裂
[8,9]；FoldX、Rosetta以及自由能微扰（free energy perturbation，FEP）等方法可评估突变对折叠稳定性的贡献
[10,11]；机器学习（machine learning，ML）和蛋白语言模型（protein language model，PLM）则可用于大规模候选突变排序
[12,13]。这些方法联合应用将有望进一步推动蛋白质工程由经验驱动向计算辅助设计转变
[10,14,15]。
GH1家族糖苷酶具有典型的磷酸丙糖异构酶（triosephosphate isomerase，TIM）桶结构，为计算模拟提供了结构规则、功能复杂的研究模型：保守桶状骨架有利于结构建模，而底物通道环区、表面静电网络和寡聚界面的差异又对构象采样、自由能计算和跨亚型模型迁移提出挑战
[16-20]。近年来，GH1酶的计算设计已由单一工具应用发展为“动态识别—能量筛选—智能组合—实验反馈”的多尺度策略
[21-23]。
本文聚焦计算模拟方法学本身，系统梳理了其在GH1家族糖苷酶高温稳定性设计中的技术路径、整合策略与发展前沿，以期计算模拟能大幅缩短研发周期、降低实验成本，推动从经验优化向理性优化的转变，为解决GH1酶热稳定性问题、工业酶快速改造和高温工业生物催化提供思路指导。
2 GH1家族酶的结构特征与热稳定性
2.1 典型的(β/α)₈ TIM桶状结构
GH1家族酶的三维结构核心是一个高度保守的TIM桶，由8个β‑α重复单元构成，内部为疏水β‑折叠核心，外部为亲水α‑螺旋（
图1）
[16,17]。从计算建模角度，TIM桶的高度对称性为结构预测与分子模拟提供了理想框架。而在GH1糖苷酶上，活性部位正好处于该桶的C端开口处
[18,19]。这种精密的排布不仅为催化反应提供了稳定的支架，还通过其独特的几何对称性确保了底物的准确定位。
连接β链和α螺旋之间的环区形成了底物通道、结合口袋以及水分子网络，是稳定性与催化功能发生耦联的关键区域
[20]。由于环区序列的不同，不同的GH1同源酶和功能亚型之间存在较大的变化，因而经常作为定点突变和动力学筛选的目标。例如，在Tamaki等
[19]针对GH1 β-葡萄糖苷酶的研究中，通过改变活性中心周围非保守氨基酸位点可以大幅度提高酶的催化效率、底物结合能力以及热稳定性。在Sengupta等
[21]的研究中，他们进一步对来源于土壤杆菌
Agrobacterium tumefaciens的GH1 β-葡萄糖苷酶H0HC94的非保守口袋残基进行了结构导向改造。Q319A 使纤维二糖比活提高约1.8倍、催化常数（turnover number，k
cat）提高约1.5倍、催化效率（k
cat/K
m，其中K
m为米氏常数）提高约3倍，但其熔解温度（melting temperature，T
m）为53.9 °C，与野生型的53.7 °C基本一致。该结果说明，环区或口袋突变可以显著改善底物可达性和催化效率，却不一定提高整体热稳定性，稳定性评价必须与动力学指标分开验证。
2.2 热稳定性的计算热力学描述
蛋白质热稳定性的计算与评估必须建立在严谨的热力学框架之上。其热力学本质是折叠态（folded state）与展开态（unfolded state）之间的焓变（ΔH）与熵变（ΔS）的平衡
[20]。折叠自由能ΔG = G
folded − G
unfolded 以及其随温度变化的规律，是进行稳定性预测的理论基础
[20]。在具体计算中，突变对稳定性的影响会通过自由能变化量ΔΔG来进行量化，其中ΔΔG定义为突变体与野生型蛋白折叠自由能之差（ΔΔG
fold = ΔG
fold,mut - ΔG
fold,WT），其中负值通常代表突变有利于稳定性的提高
[10]。从热力学角度讲，增加蛋白稳定性可以通过增加折叠态的有利相互作用（如氢键、疏水作用和静电相互作用）或降低展开态的构象自由度来实现
[24,25]。
2.3 GH1家族热稳定性的共性规律
在这些普遍稳定化因素的基础上，GH1酶还表现出与TIM桶桶口环区、催化位点及底物通道动态密切相关的结构特异性。GH1家族酶的稳定性调控具有结构特异性，其热适应机制与高度保守的TIM桶折叠密切相关。近年来大量的结构比较、MD模拟以及自由能分析表明，GH1酶热稳定性的增强一般并非仅仅依赖于单一的结构因素，而是由“TIM桶核心刚性—环区动态柔性—表面静电网络”三类因素共同决定的
[19,21,26]。本文在综合Tamaki等
[19]、Sengupta等
[21]、Niu等
[22]及Bian等
[23]等文献实验证据的基础上，将上述三类因素提炼为“TIM桶核心刚性—环区动态柔性—表面静电网络”三层协同稳定（Three-Tier Cooperative Stabilization，3TCS）概念框架。
首先，TIM桶核心区域疏水压实是保持GH1酶高温稳定性的基础。嗜热的GH1酶通常拥有更加紧密的疏水核心，β折叠之间形成了更加致密的疏水堆积，从而减少在高温下桶状骨架的构象波动；中温GH1酶往往具有更为广泛的空腔或者更大的柔性间隙，在β7-β8附近以及桶口附近都存在着较大的空腔或者柔性间隙，较容易在高温条件下发生局部解折叠
[27]。
其次，环区动态柔性调控是GH1酶与许多传统工业酶不同的重要特征。研究发现环区3、环区6和β8-α8连接区参与底物通道开闭、糖环定位和产物释放；适度降低这些区域的无序涨落可能提高熔解温度（Tm），但过度刚化会限制催化所需的构象转换。因此，环区设计的目标不是使均方根波动最小，而是保留能够连接底物进入、催化和产物释放的功能构象系综。
此外，表面盐桥网络和电荷重排在GH1高温适应上亦有重要意义。嗜热GH1酶往往具有更稳定的长程静电网络，降低了高温诱导下的局部结构涨落
[28]。
2.4 GH1家族亚型差异对稳定性设计与模型迁移的影响
GH1家族具有多种催化功能，鉴于现有计算工程研究和工业应用案例的分布，本文选取β-葡萄糖苷酶、β-半乳糖苷酶和β-甘露糖苷酶作为代表性亚型，讨论其突变设计特点及跨亚型模型迁移问题。
β-葡萄糖苷酶的工程研究最为丰富。H0HC94的Q319A突变通过扩大底物通道显著提高了纤维二糖的催化效率，但其T
m基本没有变化
[21]。Niu等
[22]利用计算重设计构建了TrCel1b-H13的人工三聚体，使其最适温度提高了55 °C，T
m提高了20 °C，在65 °C下的半衰期和水解活性分别提高416倍和3223倍。说明即使同属于β-葡萄糖苷酶，也需根据目标性质选择不同的结构层级进行优化。
β-半乳糖苷酶的突变设计更加注重半乳糖基供体与受体的定位匹配以及转糖基反应竞争关系的调控。Xing等
[29]对来源于嗜热古菌
Vulcanisaeta thermophila的GH1 β-半乳糖苷酶VthBga进行口袋突变，获得了A266G/F441T双突变体；在最优条件下，该突变体使酪醇半乳糖苷产率提高至44.5%，体现了底物通道与活性口袋几何重塑的价值。该案例主要优化的是转糖基活性而非热稳定性，也说明在不同工程目标之间不能将活性提升简单等同于稳定性增强。
β-甘露糖苷酶的工程研究主要集中于底物特异性调控，其关键在于活性位点对底物C2位羟基构型的识别。Kaper等
[30]通过交换来自
Pyrococcus horikoshii的β-甘露糖苷酶BglB与来自
Pyrococcus furiosus的β-葡萄糖苷酶CelB在活性位点附近的差异残基Gln77/Arg77和Asp206/Asn206，显著改变了酶对甘露糖苷和葡萄糖苷底物的相对催化偏好。该案例说明β-甘露糖苷酶的改造更依赖活性口袋微环境与底物识别模式的重塑，而不能简单套用以整体稳定化为目标的设计策略。由于这些位点同时影响局部氢键、水分子定位和环区构象，稳定性模型若只编码序列保守性，容易把“底物特异性位点”误判为可直接移植的稳定化位点。
现有计算模型跨亚型迁移面临了一些困难。实验标签和条件不一致，不同研究人员使用了Tm、半失活温度T50、半衰期、残余活性等不同标签，测定温度和底物也不一致；变长环区和插入/缺失导致序列和结构对齐不稳定，GH1酶为(β/α)8 TIM桶折叠，核心桶保守，但环区（尤其是活性口袋附近）和插入/缺失变异大；pH依赖的质子化状态，活性口袋pKa差异、保守水分子网络、寡聚化会显著影响局部构象动态；非葡萄糖苷酶亚型及多位点组合数据稀缺，GH1家族中β-葡萄糖苷酶数据最多，β-半乳糖苷酶、β-甘露糖苷酶等亚型实验数据少得多，导致模型在这些亚型上的不确定性高，无法有效校准。
未来可从数据、特征和验证策略三个方面改进模型迁移能力：建立按功能亚型、实验标签和测定条件分层的多任务数据集；在蛋白语言模型或图神经网络中引入亚型嵌入、口袋几何、局部氢键网络、保守水分子和寡聚状态等特征；利用少量目标亚型实验数据进行迁移学习微调，并结合MD构象系综生成动态描述符，而非仅依赖单一静态结构。
3 计算模拟方法在热稳定性研究中的应用
3.1 分子动力学模拟
MD模拟是研究GH1酶高温动态行为的最重要工具，可在原子尺度捕捉TIM桶核心波动、环区开闭以及高温诱导的氢键破裂过程
[8,9]，常用分析包括均方根偏差（root-mean-square deviation，RMSD）、均方根波动（root-mean-square fluctuation，RMSF）、主成分分析（principal component analysis，PCA）和动态交叉相关矩阵（dynamic cross-correlation matrix，DCCM），分别用于评价整体稳定性、局部柔性、主要构象运动和远程协同变化。
针对GH1(β/α)8 TIM桶结构，常用的技术参数设置包括AMBER ff14SB或CHARMM36力场、TIP3P/TIP4P水模型和周期性边界条件（periodic boundary conditions，PBC）。生产阶段应采用多个独立重复并根据研究目标选择更长时间尺度；高温模拟可设置为60-100 °C，但应与实验条件和力场适用范围相匹配。PBC盒边距、离子强度、质子化状态和寡聚体完整性均可能影响表面静电与桶口环区的结果。
对于构象跃迁较慢的环区，可采用副本交换分子动力学（replica-exchange molecular dynamics，REMD）、高斯加速分子动力学（Gaussian accelerated molecular dynamics，GaMD）或元动力学（metadynamics）增强采样。GH1环区的高柔性容易造成常规MD采样不充分，因此应通过多重复轨迹、自由能面重建和收敛检验区分真实构象变化与采样噪声。
3.2 稳定性相关能量评估与自由能计算
与传统依赖经验筛选的酶工程策略相比，稳定性相关能量评分和严格自由能计算可从不同层次评估突变效应。其中，FoldX和Rosetta主要用于经验性候选排序，MM-PBSA可辅助比较同一体系内的结合或界面相互作用，而基于热力学循环的炼金自由能方法可对少量候选进行更严格的相对自由能估算，这种方法已成为GH1家族糖苷酶理性热稳定性设计中的核心计算手段之一
[10,11]。蛋白质热稳定性的本质来源于折叠态与展开态之间自由能差的变化，而突变对稳定性的影响通常使用自由能变化值进行描述。常用方法包括基于经验势函数的FoldX、Rosetta以及分子力学/泊松-玻尔兹曼表面积法（molecular mechanics/Poisson-Boltzmann surface area，MM-PBSA）和自由能微扰（FEP）等，这些方法在多种工业酶稳定性改造研究中取得了良好应用效果
[10,31,32]。
3.2.1 FoldX
FoldX是一个以经验能量函数为基础的快速自由能计算方法软件，被广泛用于预测点突变对蛋白质折叠稳定性的影响
[33]。该方法通过对范德华作用、氢键、静电相互作用、溶剂化效应以及构象熵等能量项进行参数化建模，可高效计算突变引起的自由能变化量。在GH1家族酶的研究中，FoldX特别适用于筛选位于TIM桶外围结构、表面残基或二级结构边缘区域的潜在稳定化突变
[10]。
FoldX适合快速排除明显不稳定化的单点突变，但其静态构象近似不宜被解释为实验T
m的替代指标。以H0HC94为例，Q319A显著提高底物可达性和催化效率，而T
m与野生型基本一致
[21]。这一结果提醒：对于活性中心环区或口袋突变，应把FoldX排序与构象采样、动力学指标和实验热稳定性共同评价，不能仅凭单一自由能变化量判断“稳活双升”。
3.2.2 Rosetta
Rosetta是最为常用的蛋白质结构建模及理性设计平台之一，采用包含物理项和统计经验项的加权评分函数，并通过构象采样算法探索突变后的结构重排
[34]。Rosetta的自由能变化量计算方法能够在一定程度上考虑突变位点附近的局部主链松弛和侧链构象适配，并可对多位点组合进行整体打分
[35,36]。
Rosetta Flex_ddG通过主链采样与侧链重排，比刚性静态模型更加充分地评估环区突变及局部构象适配
[35,36]。Peccati等
[37]进一步表明，使用AlphaFold构象系综而非单一结构能够改善酶稳定化突变预测。对于GH1多位点组合，Rosetta仍主要提供候选排序而非可靠的高阶上位效应定量；关键组合应结合MD系综、显式溶剂自由能计算和实验验证。
3.2.3 FEP与MM-PBSA
自由能微扰通过热力学循环和多个λ窗口来计算高精度相对自由能。典型流程包括拓扑构建、数十至数百 ns的平衡MD、多窗口采样，并使用Bennett接受比（Bennett acceptance ratio，BAR）或热力学积分（thermodynamic integration，TI）估算自由能；收敛性应通过正反向滞回、窗口重叠和块平均误差进行评估。
MM-PBSA适用于从充分采样的MD轨迹中提取构象并进行中通量相对排序，其结果由分子力学能、极性/非极性溶剂化和熵项组成。对GH1环区和寡聚界面突变，构象熵、保守水和质子化状态的不确定性往往是主要误差来源，因此更适合作为同一体系内的相对筛选，而非跨亚型绝对比较。
3.2.4 GH1体系中自由能计算方法的综合比较
FoldX、Rosetta、MM-PBSA和FEP在GH1糖苷酶热稳定性设计中各有侧重（
表1）。FoldX计算速度快，适合TIM桶外围表面位点和疏水核心压实的单点突变初筛，但因主要基于静态结构，对GH1高度柔性的环区预测能力有限。Rosetta可通过局部主链采样和侧链重排处理一定程度的构象适配，因而更适合环区及多位点候选组合的结构精修，但其对高阶上位效应的定量预测仍受到采样范围和评分函数误差的限制。MM-PBSA适合中通量突变体排序，而FEP具有较严格的热力学基础，可用于少量关键候选的精细计算，但其可靠性取决于热力学循环构建、采样充分性和收敛性。总体而言，推荐采用“FoldX快速初筛 + Rosetta环区精修 + MM-PBSA批量排序 + FEP关键验证”的分层策略，以平衡计算效率、预测精度与GH1体系的动态特性需求。
3.3 机器学习与人工智能方法
3.3.1 基本原理与模型类型
机器学习和人工智能（Artificial Intelligence，AI）技术在蛋白质工程领域快速发展，是研究酶热稳定性预测及优化的重要手段之一。其核心思想就是通过从大量实验数据、序列信息和结构特征中学习模式，建立预测模型来评估突变对蛋白质稳定性的影响
[12]。ML方法与传统计算方法不同，它能够从大量实验数据中自动学习序列—结构—稳定性之间的复杂非线性关系，无需显式建模分子层面的物理相互作用，在处理高维度序列空间和大规模突变筛选方面具有明显优势
[38]。
在热稳定性研究中，ML模型一般基于监督学习框架，输入包括氨基酸序列、进化信息、结构描述符（如二级结构比例、氢键网络），输出为热稳定性指标如T
m或自由能变化量
[39]。近年来，深度学习和蛋白语言模型的发展，进一步提升了模型对复杂蛋白质动态的捕捉能力
[13]。
3.3.2 GH1家族酶的ML应用案例
Pudžiuvelytė等
[40]开发的TemStaPro可直接根据氨基酸序列预测蛋白质在多个温度阈值下的稳定性类别。该模型主要基于已知最适生长温度的微生物蛋白质组数据进行训练，并在CRISPR-Cas II类效应蛋白中进行了应用验证。对于GH1家族，该工具可用于候选同源酶的初步筛选，但其在不同GH1功能亚型中的预测精度仍需通过独立实验数据验证。
Erkanli等
[41]构建了一个序列与温度联合驱动的三模块机器学习框架，用于预测β-葡萄糖苷酶在不同温度下的k
cat/K
m。该模型反映的是温度依赖的催化性能，而非蛋白质热力学稳定性。因此，该工作可用于辅助评价酶的催化热适应性，但不能替代T
m、T
50或半衰期等热稳定性实验。
Rapp等
[42]开发的自驱动蛋白质景观探索平台（Self-driving Autonomous Machines for Protein Landscape Exploration，SAMPLE）整合了机器学习与自动化实验，实现了蛋白质热稳定性的自主优化。该平台部署了4个智能代理，以工程改造糖苷水解酶的热耐受性为目标。每个代理每轮选择3个序列进行测试，共运行20轮设计—测试—学习—分析循环。尽管各代理的搜索行为存在差异，但均能快速收敛于热稳定酶。该平台展示了AI驱动的闭环优化在酶热稳定性工程中的潜力。
3.3.3 局限性与改进方向
尽管ML和AI在GH1糖苷酶热稳定性预测中取得了显著进展，但目前相关研究仍面临明显局限。
首先，训练数据偏差是限制模型泛化能力的核心问题。当前多数公开稳定性数据库中，去稳定化突变远多于稳定化突变，而GH1家族高质量实验数据本身数量有限，因此模型容易学习到“突变通常降低稳定性”的统计偏差。此外，不同GH1亚型之间在底物特异性、环区长度及活性位点微环境上的差异，也进一步增加了跨家族泛化难度
[43]。
其次，虽然深度学习模型能够实现较高预测精度，但往往无法解释为何某些突变会提高热稳定性。对于GH1酶这种具有复杂催化动态的体系，仅获得预测结果而缺乏机制解释，难以真正指导工业酶工程设计
[38]。
此外，对于罕见突变或者多个位点协同效应的处理，现有深度学习模型往往需要大量计算资源进行训练，且在小样本数据集上易过拟合
[44]。
近年来，越来越多研究开始尝试将MD轨迹、Rosetta能量项以及进化耦联信息作为多模态输入，以提高模型对动态构象及长程上位效应的学习能力。此外，从GH1亚型特异性角度，当前ML模型的训练集主要来自β-葡萄糖苷酶，对β-半乳糖苷酶和β-甘露糖苷酶的预测能力尚未得到系统性验证。鉴于β-甘露糖苷酶底物结合口袋几何差异较大，建议在将当前ML模型应用于非β-葡萄糖苷酶亚型时，优先进行迁移学习微调（fine-tuning）而非直接外推。未来，结合物理约束与动态信息“Physics-informed AI”可能成为GH1热稳定性预测的重要发展方向
[45]。
3.4 多尺度集成计算策略
3.4.1 策略框架与核心原理
多尺度集成计算策略（Multi-scale Computational Strategies）是蛋白质热稳定性理性设计的一种高级策略，通过整合分子动力学模拟、自由能计算、机器学习/蛋白语言模型以及AlphaFold结构预测等方法，构建多尺度、多模态的集成计算框架，从而显著提升稳定性预测与突变设计的效率和可靠性
[23,46,47]。其核心原理在于充分发挥各方法的互补优势：AlphaFold等深度学习工具为无实验结构的目标酶提供可靠的初始模型
[12,37]；MD模拟提供原子级动态构象信息和识别热敏感柔性区域；Rosetta和FoldX可用于突变稳定性变化的近似评估、候选排序及局部结构优化；机器学习实现高通量模式识别、上位效应预测和自主迭代设计。最后以T
m、半衰期、活性及工业耐受性实验验证并更新模型（
图2）。
3.4.2 GH1家族酶的多尺度集成案例
3.4.2.1 多种MD方法解析单点稳定化机制
Rocha等
[48]以
Paenibacillus polymyxa GH1 β-葡萄糖苷酶BglA的E96K和M416I稳定化突变为对象，联合不同温度下平衡MD、PCA、分子对接、元动力学、加速MD、泊松-玻尔兹曼表面分析、网格非均匀溶剂化理论（grid inhomogeneous solvation theory）和构象熵计算，并与实验表型相互印证。结果表明，E96K主要通过表面静电重排及疏水核心溶剂隔离增强稳定性，M416I则主要通过重新分配活性位点环区动力学发挥作用。该研究展示了“实验突变—多轨迹动力学—自由能景观—机制解释”的完整链条，但其候选位点并非由AI组合模型生成。
3.4.2.2 计算重设计寡聚界面实现稳定性与活性协同提升
Niu等
[22]通过计算蛋白质界面设计，在GH1 β-葡萄糖苷酶TrCel1b-H13中构建人工三聚体界面，并结合结构分析和实验验证评估了界面装配效率及寡聚体热稳定性。经过改造后最适温度提高了55 °C，T
m提高20 °C；在65 °C下，半衰期和水解活性分别提高416倍和3223倍。该案例说明，除单点残基筛选外，通过理性设计寡聚化界面建立跨亚基相互作用，可产生幅度更大的稳定化效应。
3.4.2.3 活性口袋改造揭示稳活目标解耦
Sengupta等
[21]基于结构与序列保守性选择H0HC94口袋位点，并以Q319A进行实验验证。该突变使纤维二糖催化效率提高约3倍，而T
m几乎不变，半衰期仅轻微提升（1.4倍）。虽然该研究没有覆盖FEP和AI组合突变，但它提供了多尺度流程中不可忽略的负面校准：能量评分或活性改善不能代替T
m和半衰期测试，稳定性与活性多目标优化应在筛选阶段分别建模。
3.4.2.4 PLM组合突变与动态上位效应分析的迁移范式
Bian等
[23]在肌酸酶体系中把低阶组合突变数据、温度引导PLM迭代微调、全组合空间排序、两轮实验验证、MD和DCCM机制分析连接起来，获得13位点高稳定突变体。该研究不是GH1案例，但其流程最接近"AI组合突变—实验反馈—动态上位效应解释"，可作为GH1在积累足够单点和低阶组合数据后的迁移模板。
3.4.3 现存挑战
虽然多尺度集成计算方法大大提高了热稳定性的预测精度以及突变设计成功率，但仍存在若干典型的方法边界与失效场景，有必要系统归因。
力场偏差与采样不足会在高柔性环区、带电表面和保守水网络中被放大。Peccati等
[37]表明，使用AlphaFold构象系综可改善Rosetta稳定化预测，反向说明单一静态结构难以覆盖突变诱导的构象重排。对GH1环区应采用多初始构象、多重复轨迹和构象分层评分。
多位点组合不能用单点ΔΔG简单相加。溶剂介导的静电耦联、疏水核心重排、寡聚界面装配和熵补偿均可产生符号或协同上位效应
[23,44,49]。因此，组合突变应先用低成本模型缩小空间，再对高不确定性组合开展Rosetta Flex_ddG、MD系综或FEP复核。
针对上述边界，推荐采用“FoldX初筛→Rosetta Flex_ddG精修→FEP验证关键位点”分层策略降低预测偏差。此外，跨尺度数据融合与参数优化仍面临诸多挑战。一方面多方法集成往往需要跨尺度数据融合与复杂参数调整，这对计算资源和算法稳定性都有较高的需求
[50]。另一方面不同尺度方法之间的结果整合和权重分配也没有统一的标准，可能引入系统性偏差。特别是在模拟极端工业条件时，由于实验数据不足或训练样本有限，模型预测可能出现过拟合或泛化能力下降的问题
[51]。此外，对于新发现的GH1同源酶，当前深度学习和机器学习模型在跨物种或跨功能背景下的稳定性预测能力仍有待进一步提升
[50]。
4 挑战、工业转化与未来展望
4.1 当前面临的关键挑战
4.1.1 数据偏差与模型泛化
训练数据的系统性偏差是ML模型泛化能力的重要限制因素。当前的标准数据库（如ProTherm、FireProtDB）中去稳定化突变通常占据多数，而稳定化突变往往不足三分之一
[51]。同时，多位点组合突变可能产生显著的上位效应（epistasis），使其稳定性变化偏离各单点突变效应的简单加和，这使得基于单点突变数据训练的模型在预测组合突变体时可靠性显著下降
[44,52]。对于GH1家族而言，不同亚型（β-葡萄糖苷酶、β-半乳糖苷酶、β-甘露糖苷酶）之间的序列和结构差异进一步加剧了泛化难度。
4.1.2 上位效应的结构来源与计算优化
GH1多位点突变的上位效应来源于TIM桶内部和外围的多层耦联。首先，相邻β-α单元形成贯穿桶核心的疏水堆积网络，一个位点引起的空腔或侧链重排可跨重复单元传播；其次，β3-α3、β6-α6环区及C-末端β8-α8连接区通过协同运动控制底物通道，远距离突变可改变其开闭概率；再次，表面盐桥、保守水分子、催化残基质子化状态及寡聚界面构成环境依赖的静电网络。催化谷氨酸及周围水化结构还会使突变效应随底物结合状态而变化。上述因素使组合突变效应具有非加和性，甚至出现单点有利而组合失活的符号上位效应。
现阶段可采用的五类计算优化包括：构建单点—双点—低阶组合的分层训练集，显式学习二阶和高阶交互项；在PLM或图神经网络中加入残基距离、接触图、寡聚界面和亚型信息；用MD构象系综、DCCM和互信息特征表征动态长程耦联；按精度分层复核，首先利用低成本序列模型、FoldX或Rosetta稳定性评分缩小候选空间，再通过MD构象系综和相应能量分析复核高不确定性组合；对于少量关键突变，可进一步采用严格热力学循环方法进行精细计算，并最终通过实验验证。通过主动学习或贝叶斯优化选择最有信息量的组合进行实验。Bian等
[23]以低阶组合数据微调PLM并用DCCM解释长程上位效应，为GH1提供了可迁移范式，但仍需建立GH1专属的低阶组合突变数据集，以验证该范式的预测精度与机制解释能力。
4.1.3 稳定性—活性权衡
“稳定性—活性权衡”是GH1家族糖苷酶热稳定性设计中最核心的问题之一，也是限制工业化应用的重要瓶颈。
传统蛋白稳定化理论通常认为，提高整体刚性能够增强热稳定性。然而对于GH1酶而言，其催化功能本身高度依赖局部动态构象变化，因此过度刚化往往会破坏催化所需的功能柔性。
GH1酶属于典型的retaining glycosidase
[18]，其催化过程涉及底物精确定位、糖基中间体形成、水分子重排以及产物释放等多个动态步骤。在这一过程中，活性中心附近环区需要发生精细的开闭运动，以实现底物进入与催化残基协同定位。因此，GH1酶催化体系本质上是一种“动态催化系统”。
GH1酶的热稳定性优化需在核心刚性与功能柔性间取得平衡。通过选择性刚化策略—即稳定TIM桶骨架的同时保留活性中心动态自由度—可实现稳定性与催化活性的协同提升
[22,23,45,46,53]。
此外，多点协同突变设计也开始突破传统“稳定性—活性”对立模式。一些AI驱动设计研究表明，通过合理协调核心刚性、环区柔性以及表面静电网络，可以实现稳定性与催化活性的同步提升。这说明GH1酶的热稳定性优化并非绝对受限于活性损失，而是依赖于更加精细的动态构象调控。
因此，未来GH1酶工程设计的关键方向将不再是“全面刚化”，而是通过选择性刚化策略实现局部稳定化与功能动态之间的协同平衡，即在增强整体热稳定性的同时保留催化所需的动态自由度。
4.1.4 计算成本与效率瓶颈
多尺度集成策略虽然预测精度高，但计算成本巨大。一个完整的MD + Rosetta + ML流程可能需要数万CPU·h，限制了其在工业规模酶改造中的应用。例如，对一个典型的GH1酶进行全面的Rosetta饱和突变扫描，即使采用并行计算也需要数周时间。开发更高效的近似算法（如基于图神经网络的快速自由能变化量预测）和云计算平台是未来的重要方向。
4.2 新兴计算工具与发展趋势
4.2.1 物理-AI融合模型
近年来，物理约束与人工智能的融合成为蛋白质设计的新趋势。Gelman等
[54]开发的突变效应迁移学习（mutational effect transfer learning，METL）框架通过引入生物物理先验（如能量函数、构象约束），在模拟数据上预训练并结合实验功能数据微调，在小样本泛化和残基位置外推等任务中显著优于传统基于进化信息的模型。Li和Luo
[55]提出的基于重连策略的稳定性预测（stability prediction using a rewired strategy，SPURS）框架将进化尺度建模（Evolutionary Scale Modeling，ESM）系列的第二代蛋白语言模型ESM-2与反向折叠模型ProteinMPNN相结合，在大规模热稳定性数据上进行监督微调，旨在实现对未见蛋白及未见突变的自由能变化量预测。这些工具对训练数据有限的GH1糖苷酶体系具有重要意义。
4.2.2 生成式AI与从头设计
近年来，生成式人工智能为蛋白质结构与序列的协同设计提供了新工具。主链扩散模型RFdiffusion可通过迭代去噪生成满足特定拓扑、对称性、结合位点或功能基序约束的蛋白质骨架，而ProteinMPNN可根据给定骨架开展逆折叠序列设计，二者共同构成“骨架生成—序列优化”的设计流程
[56,57]。这一策略为探索GH1样TIM桶骨架和高稳定候选序列提供了新路径
[58,59]。但GH1酶的催化功能不仅取决于整体骨架稳定性，还受到催化残基几何排布、过渡态稳定、保守水分子网络及底物通道动态的共同制约。因此，现有生成式工具尚不能直接保证所设计序列具有预期的糖苷水解活性和底物特异性，相关候选仍需结合催化构象建模、分子模拟及实验验证，现阶段更适合作为中长期研究方向。
4.2.3 多目标工业适应性优化
工业GH1糖苷酶在保证热稳定性的同时还要兼顾葡萄糖耐受性、有机溶剂稳定性及宽pH适应性等复杂性能指标
[60,61]。推动计算设计由单一稳定性优化向多维性能协同优化转变，需要开发多目标优化算法和相应的实验验证平台。
4.2.4 动态构象系综设计
传统的蛋白质设计主要基于单一静态结构，而越来越多的证据表明，蛋白质的“动态构象系综”而非单一构象决定了其功能。未来的计算设计将更多考虑蛋白质的动态特性，通过稳定“功能性构象”而非简单增加刚性来实现稳定性与活性的双赢。这需要发展新的计算工具，能够同时预测蛋白质的动态行为和热力学稳定性。
针对GH1酶特有的工业痛点，基于当前技术成熟度，以下方向具有近期突破潜力。
葡萄糖耐受性与热稳定性的耦联设计—工业β-葡萄糖苷酶在高温下同时面临热失稳和产物（葡萄糖）抑制的双重制约，构建能够同时预测Tm变化与葡萄糖抑制常数变化的双目标自由能变化量模型将是近期重要突破口。
GH1多亚型通用计算模型构建—开发基于泛GH1结构保守域的多任务学习框架，覆盖β-葡萄糖苷酶、β-半乳糖苷酶和β-甘露糖苷酶三类亚型，以克服单亚型模型的泛化瓶颈，显著扩展计算设计的工业应用范围。
5 结论与展望
5.1 结论
GH1糖苷酶的高温稳定性设计本质上是保守TIM桶骨架、功能环区和表面/界面相互作用的协同优化。MD与增强采样适于识别热敏感构象和远程耦联，自由能方法适于分层排序，ML与PLM适于扩大搜索和组合突变；但任何单一工具都不能同时可靠处理跨亚型分布偏移、动态上位效应和稳定性—活性权衡。本文提出的3TCS框架将核心刚性、环区动态柔性和表面静电网络统一为可计算的机制约束，并通过代表性案例说明结构层级和目标性质必须与工具选择相匹配。
5.2 展望
针对数据偏差，应建设GH1专属标准化数据库，统一记录序列、结构、亚型、底物、pH、温度、时间、Tm、T50、半衰期和活性，并纳入失败突变及重复实验；模型评价应采用家族外和亚型外测试，而非仅随机划分。
针对上位效应，应系统积累双点和低阶组合数据，发展融合接触图、MD构象系综和DCCM特征的结构感知PLM/图模型，并以不确定性驱动的主动学习选择下一轮实验；高阶组合需用显式溶剂方法和实验进行最终验证。
针对稳定性—活性权衡，应把Tm、半衰期、kcat/Km、底物特异性和工业耐受性纳入Pareto多目标优化，优先稳定TIM桶核心并抑制非功能性表面的过度波动，同时保留底物通道环区的功能构象系综。
针对计算成本，应建立“PLM/FoldX初筛—Rosetta构象系综精修—MD/MM-PBSA排序—FEP严格验证—实验验证”的分层流程，并结合代理模型、GPU计算和自驱动实验平台形成可迭代闭环。通过这四条路线，GH1酶工程有望由分散的计算辅助逐步转向数据标准化、机制可解释和成本可控的理性设计。
PDF
