ME丨天津科技大学王敏教授团队基于人工智能的高精度酶约束模型的自动构建

在系统生物学中，高质量的基因组规模代谢模型（Gems）对微生物细胞工厂的合理设计至关重要，但传统模型构建方法在非模式工业微生物上存在挑战。这主要是因为非模式菌株存在蛋白质注释不充分以及酶动力学参数缺失的问题。尽管有一些基因功能预测工具，但它们依赖现有数据库和已知功能注释进行预测，对于未被广泛研究的蛋白质可能不适用，且在准确注释未知蛋白质功能方面仍存在困难。随着“下一代细胞工厂”概念的出现，在生物炼制应用中展现出独特优势。G. thermoglucosidasius作为一种具有高温耐受性、生长迅速且污染风险低等优点的微生物，同时还具备高效的基因编辑工具，是极具潜力的下一代细胞工厂。然而，由于缺乏代谢信息，其代谢工程需求难以通过传统方法得到满足。该研究提出基于人工智能技术的BIM策略，通过一系列工具和模块完成G. thermoglucosidasius高精度酶约束模型的构建及优化。该研究为快速构建非模式菌株模型提供了新的策略和方法，展示了人工智能在代谢工程领域的巨大潜力，对推动非模式微生物在工业生产中的应用具有重要意义。

基因组规模代谢网络模型的改进

基于现有的G.thermoglucosidasius基因组规模代谢网络p-thermo模型覆盖酶数量有限的问题。该研究首先通过EnzymeExtractor程序从其基因组中提取相关信息并整理为FASTA文件，再用CLEAN进行预测分析，得到3641个蛋白质的功能注释，经处理和去重后确定1356个独特功能蛋白质。对比这些数据与已知蛋白注释及p-thermo模型中的酶，通过将预测的功能蛋白与p-thermo模型中现有的酶进行比较，发现112种酶的预测与基因组注释一致，但未包括在模型中。此外，从预测中鉴定出置信区间大于0.9的205种酶，从而从预测中鉴定出317种新酶。去除重复后，添加249种酶以增强模型完整性（图1和图2）。为避免在模型构建中纳入非特征性反应，研究系统地汇编了目标属特异性的相关基因、酶、化合物和反应，为G. thermoglucosidasius NCIMB 11955创建了一个全面的生化反应数据库，并与目标属特异性数据库进行交叉引用。同时，开发了MetaPatchM模块来自动化补充基因组规模代谢网络模型。该模块优先纳入注释基因，结合CLEAN预测信息，保留未注释且置信区间大于0.9的酶，并将其与同属反应数据库匹配，最终确定308个反应来更新模型。更新后的模型包含1484个反应和1253种代谢物，并依据该模型构建了代谢通量矩阵。

图1 基因组注释结果分析

图2 p-thermo模型、G. thermoglucosidasius NCIMB 11955基因组和CLEAN预测结果的比较

酶约束模型的构建

基因组规模模型的预处理。研究首先以G. thermoglucosidasius G-ther11955模型为基础，采用GECKO方法，利用酶动力学参数（k_cat）和酶丰度数据对反应施加通量约束。接着，将模型转换为raven格式，以便后续对代谢网络模型进行更便捷的构建、分析和可视化操作。

可逆反应和同工酶反应的拆分。为解决酶与反应之间的复杂关系（如同工酶、多功能酶以及单个酶催化可逆反应等情况导致的约束复杂问题），研究在酶约束模型中对可逆反应采用特殊处理技术。将边界参数为（-1000，1000）的可逆反应R拆分成边界参数为（0，1000）的R和R_rev两个不可逆反应，G-ther11955模型中最初的644个可逆反应拆分后变为2128个反应，代谢物和基因数量不变。同时，对涉及同工酶等的反应引入过渡反应并进行拆分，168个涉及同工酶的反应被分成614个反应，模型反应数最终扩展到2574个。

酶动力学参数的积分。为更准确地预测细胞代谢过程，研究人员首先将酶当作底物，把kcat的倒数作为系数纳入模型。接着，采用e_i等参数描述模型中反应所需酶量与其他生物参数关系。因缺乏非模式生物蛋白质丰度数据，通过敏感性分析确定e_pool值为94 mg/g DW，参考GECKO 3.0设定为0.5，将酶作为伪代谢物引入模型模拟酶循环。

酶约束代谢网络模型的校准。首先，该研究通过整合酶动力学参数和丰度信息建立了酶约束模型，并开发了Enzyme-insert程序来自动化此过程。为解决酶约束代谢网络模型构建中因细胞简化表述及模型与实际细胞代谢差异导致的过约束问题，该研究对Enzyme-insert模块进行优化。优化后的模块在添加酶约束后对模型进行全局分析，通过执行通量平衡分析（FBA）评估每个反应，特别关注FBA计算中通量为零或接近零的反应，这些反应被认为受到过度约束。然后识别并消除这些约束，以提高模型的准确性和功能性。经过全局扫描模型，确定了359个抑制性反应，包括62个同工酶和54个多功能酶反应，这些反应导致了通路阻塞。此外，该研究根据酶的基因组表示情况添加连接反应并纳入总酶池约束，使模型反应总数增加到3345个。最后，研究应用Enzyme-insert工具去除359个过约束反应对应的酶连接，得到包含2986个反应的酶约束模型ec_G-ther11955，涵盖2986个化学反应和1834种代谢物（图3）。新增代谢物涉及700个酶和由168个同工酶引入产生的过渡代谢物，这些酶广泛参与多种代谢途径，得到168个同工酶和443个多功能酶（表1）。

图3 模型构建过程中的代谢信息

表1 ec_G-ther11955模型特征描述

G.thermoglucosidasius NCIMB 11955酶限制模型的特征分析

研究首先将ec_G-ther11955模型与p-thermo模型比较，发现二者拓扑结构相似（表2），但也存在差异。ec_G-ther11955模型全局聚类系数更低、平均最短路径更长、网络结构更分散、复杂且连通性更高（图4）。此外，采用MEMOTE工具测试仅含化学计量数据的模型，去除相关反应并标准化名称后，评估分数为94%，略高于P-thermo模型的90%。研究采用两种模型预测核黄素生产的优化目标，通过OptForce和PKT算法进行预测。ec_G-ther11955模型预测得到purAGtg基因是目标基因，而P-thermo模型不能预测核黄素为目标基因。研究人员通过为P-thermo模型添加核黄素出口反应使其能预测核黄素为目标基因后，发现该模型预测的微生物生长和核黄素生产情况与实验结果不符（图4）。接着采用ec_G-ther11955模型预测预测核黄素产生的敲除靶标，鉴定了4个基因敲除靶标（AOT13_03190、AOT13_05975、AOT13_05975和AOT13_05975），其中三个已被报道对提高核黄素产量有效，准确率至少达到75%。

表2 原始代谢模型和酶约束模型的连通性度量

图4 原始代谢p-thermo模型和酶约束ec_G-ther11955模型的节点度值直方图

L-缬氨酸修饰靶点的预测

核黄素生产案例初步验证了ec_G-ther11955模型在识别G. thermoglucosidasius代谢途径修饰的基因敲除靶点方面的潜力。在此基础上，选择缬氨酸合成作为代谢工程的目标来进一步验证模型。研究采用OptForce、PKT和FSEOF算法预测在维持正常细胞生长的同时实现高产的靶点。初始预测得到9个基因敲除靶点和8个基因过表达靶点（图5和表3）。

图5 利用G. thermoglucosidasius NCIMB 11955改造的枯草芽孢杆菌生产L-缬氨酸的改进策略

表3 高产缬氨酸单靶点预测

在此基础上，进一步模拟该模型以获得组合修饰靶点。从模拟后细胞存活开始并考虑维持正常细胞生长的角度出发，预测了7个组合修饰靶标（表4）。使用已知的调节和代谢信息分析这些靶标的代谢特征，包括减少前体消耗、切断竞争途径、最小化碳损失、平衡辅因子、增加NADPH水平和将碳流重定向到所需产物。最后，该研究进行了有针对性的遗传修饰，包括基因敲除和过表达，以提高菌株的L-缬氨酸产量。

表4 高产L-缬氨酸联合靶点预测

提高L-缬氨酸产量的菌株代谢工程

该研究首先对基因敲除菌株和过表达菌株进行摇瓶发酵实验。敲除基因实验中，除ΔsdhB菌株外，多数敲除菌株L-缬氨酸产量显著提高，8个敲除基因（gnd、fum、akgdh、pflB、pdhA、pyc、mdh和ldhA）验证了模型88.9%的预测准确性并发现新靶点（图6）。

图6 预测靶标敲除对PKT和OptForce生产L-缬氨酸的影响

过表达菌株实验中，过表达菌株OilvBNC实现了L-缬氨酸产量为112 mg/L，相较于携带空载体pucG3.8的对照菌株提高了100.6%。其他过表达基因如ilvD、gpsA等分别使L-缬氨酸产量有不同程度的提高，而glus基因过表达对产量无显著影响，glpk基因过表达使产量降低约35.2%（图7）。

图7 基于强制目标通量的通量扫描预测靶标过表达对L-缬氨酸生产的影响

基于单靶点修饰的结果，选择了四个有效的基因敲除靶点进行组合修饰实验，包括mdh、ldhA、pdhA和sdhA基因。接着，进行摇瓶发酵实验验证组合修饰的效果。实验结果显示，组合敲除的菌株在L-缬氨酸产量上有显著提升，相较于野生型菌株和单靶点敲除菌株，产量进一步提高，验证了组合修饰策略的有效性（图8）。

图8 预测靶点对L-缬氨酸生产的影响