2.3.1 构建外源代谢物与外源反应

基于Yeat-GEM9模型，添加人参皂苷合成和AHG转化的外援反应，共13种代谢物，16种反应，从而构建出一个可利用葡萄糖或半乳糖和AHG，合成Rh1的酿酒酵母代谢模型。

部分所需物质和代谢途径并非酵母中天然具有，我们额外引入到模型中，下表为使用的新代谢物。

ID	Name	Cellular Localization	Kegg_id
s_5000	Dammarenediol_c	Cytoplasm	C19829
s_5001	Dammarenediol_er	Endoplasmic reticulum	C19829
s_5002	Protopanaxadiol_c	Cytoplasm	C20715
s_5003	Protopanaxadiol_er	Endoplasmic reticulum	C20715
s_5004	Protopanaxatriol_c	Cytoplasm	C20715
s_5005	Protopanaxatriol_er	Endoplasmic reticulum	C20715
s_5006	Ginsenoside Rh1_c	Cytoplasm	C22129
s_5007	Ginsenoside Rh1_er	Endoplasmic reticulum	C22129
s_5008	Ginsenoside Rh1	Extracellular	C22129
s_5009	3,6-Anhydro-D-Galactose	Cytoplasm	C06474
s_5010	3,6-Anhydro-D-Galactonate	Cytoplasm	-
s_5011	2-Dehydro-3-Deoxy-D-Galactonate	Cytoplasm	-
s_5012	2-Keto-3-Deoxy-6-Phosphohexonate	Cytoplasm	-

下表为使用的新反应。

ID	Name	Reaction	Kegg_id
r_5000	DM Synthesis(c)	(S)-2,3-Epoxysqualene_c + H2O_c <-> Dammarenediol_c	R09701
r_5001	DM Synthesis(er)	(S)-2,3-Epoxysqualene_er + H2O_er <-> Dammarenediol_er	R09701
r_5002	PPD Synthesis(c)	Dammarenediol_c + NADPH_c + Oxygen_c <-> Protopanaxadiol_c + NADP_c + H2O_c	R10540
r_5003	PPD Synthesis(er)	Dammarenediol_er + NADPH_er + Oxygen_er <-> Protopanaxadiol_er + NADP_er + H2O_er	R10540
r_5004	PPT Synthesis(c)	Protopanaxadiol_c + NADPH_c + Oxygen_c <-> Protopanaxatriol_c + NADP_c + H2O_c	R10591
r_5005	PPT Synthesis(er)	Protopanaxadiol_er + NADPH_er + Oxygen_er <-> Protopanaxatriol_er + NADP_er + H2O_er	R10591
r_5006	Rh1 Synthesis(c)	Protopanaxatriol_c + UDP-D-glucose_c <-> Ginsenoside Rh1_c + UDP_cytoplasm	R12388
r_5007	Rh1 Synthesis(er)	Protopanaxatriol_er + UDP-D-glucose_er <-> Ginsenoside Rh1_er + UDP_er	R12388
r_5008	Rh1 Transport	Ginsenoside Rh1_c <-> Ginsenoside Rh1_er	-
r_5009	Rh1 Synthesis	Ginsenoside Rh1_c ->Ginsenoside Rh1	-
r_5010	AHG degradation	3,6-Anhydro-D-galactose + NADP(+) -> 3,6-Anhydro- D-galactonate + NADPH	-
r_5011	AHGA degradation	3,6-Anhydro-D-galactonate -> 2-Dehydro-3-Deoxy- D-Galactonate	-
r_5012	KDGal degradation	2-Dehydro-3-Deoxy-D-Galactonate -> 2-Keto-3- Deoxy-6-Phosphohexonate	-
r_5013	KDPG degradation	2-Keto-3-Deoxy-6-Phosphohexonate -> Pyruvate + Glyceraldehyde 3-Phosphate	-
r_5014	Rh1 Exchange	Ginsenoside Rh1->	-
r_5015	AHG Exchange	AHG ->	-

下表为采用的相关反应通量的上下限设置。

Attribute	Reaction(Rxn)	Description	Lower Bound	Upper Bound	Units
Key Variables	DM Synthesis(c)	Provided Oxygen	a	0	mmol/gDW/h
	DM Synthesis(er)	Provided Glucose	b	0
	PPD Synthesis(c)	Provided Galactose	c	0
Irrelevant Variable	EX_nh4_e	Unlimited Supply	-	-
	EX_fe2_e		-	-
	EX_pi_e		-	-
	EX_k_e		-	-
	EX_na1_e		-	-
	EX_so4_e		-	-
	EX_cl_e		-	-
	EX_cu2_e		-	-
	EX_mn2_e		-	-
	EX_zn2_e		-	-
	EX_mg2_e		-	-
	EX_ca2_e		-	-

2.3.2 寻找目标函数比例

细胞资源有限，其必须在生长与产物合成之间进行调配。为了避免权重较小的反应被权重大的一方完全夺去模型的注意力，我们需要找到反应体系下的平衡权重 $ \alpha $ . 在FBA中，可将生物量合成和Rh1合成这两个反应构建带权重的目标函数。公式如下： $$Maximize(Z)=\alpha * Flux(Biomass)+(1-\alpha)*Flux(Rh1) \tag{1}$$

$ Z $ 是目标函数，$\alpha$ 是菌株生长的权重（范围为0-1），$ Flux(x)$ 代表 $ x $ 反应的通量。同时设定碳源浓度根据 $\alpha$ 从0-1，进行计算。

                
list_growth=[]
for x in range(0,1010,10): weight = (x/1000, 1-x/1000) 
  # Oxygen input
  oxygen_exchange_rxn_id='r_1992'
  oxygen_exchange_rxn = model.reactions.get_by_id(oxygen_exchange_rxn_id) 
  oxygen_exchange_rxn.bounds=(-1000,0)

  # Glucose input 
  glucose_exchange_rxn_id='r_1714'
  glucose_exchange_rxn = model.reactions.get_by_id(glucose_exchange_rxn_id) 
  glucose_exchange_rxn.bounds=(0,0)

  # Galactose input 
  galactose_exchange_rxn_id = 'r_1710' 
  galactose_exchange_rxn = model.reactions.get_by_id(galactose_exchange_rxn_id) 
  galactose_exchange_rxn.bounds = (-0.5216, 0)

  # AHG input
  AHG_rxn_id = 'r_5015'  # AHG exchange
  AHG_rxn = model.reactions.get_by_id(AHG_rxn_id) 
  AHG_rxn.bounds = (-0.5216, 0)

  # Ginsenoside biosynthesi 
  rh1_synthesis_rxn_id='r_5014'
  rh1_synthesis_rxn = model.reactions.get_by_id(rh1_synthesis_rxn_id) 
  rh1_synthesis_rxn.bounds=(0,1000)

  # Objective
  biomass_rxn_id = 'r_2111'# biomass r_2111/r_4041/r_4048 
  biomass_rxn = model.reactions.get_by_id(biomass_rxn_id) 
  biomass_rxn.bounds=(0,1000)

  etoh_rxn_id='r_1761'
  etoh_rxn = model.reactions.get_by_id(etoh_rxn_id) 
  etoh_rxn.bounds=(0,1000)

  # Weighted solution r_5015 
  model.objective = weight[0]*biomass_rxn.flux_expression+weight[1]*rh1_synthesis_rxn.flux_expression
  model.solver.objective_direction = 'max'

  # FBA optimizing 
  solution = model.optimize()

  list_growth.append((
      f'{x / 1000}:{round(1 - x / 1000, 3)}', 
      solution.fluxes["r_2111"], # biomass flux
      solution.fluxes["r_5014"],  # Rh1 flux
      weight[0] * solution.fluxes["r_2111"] + weight[1] * solution.fluxes["r_5014"] # objective value
    ))
  print(list_growth)
                
              

作出最大化目标函数 $ Z $ 的图像如下：

同时求解出平衡 $ \alpha $ =0.55。平衡 $ \alpha $ 在权重设置下，目标函数值最小的点。选择该点对应的权重值作为目标函数的权重，可以让目标函数处于一个相对平衡的状态，从而避免目标函数的不平衡状态，避免单一目标（生物量或产物合成）过于主导。

观察到生长速率曲线左侧和Rh1曲线右侧，由于生长函数的不平衡状态，$ \alpha $ 的变化对其未有影响，更多受另一反应的主导影响。

2.3.3 测试不同碳源策略下，Rh1的合成量

$ \alpha $ = 0.55 ，在不同浓度的碳源下，最大化目标函数 $ Z $ 条件下，计算Rh1合成反应的速率，结果如下图所示

从图中可以看出，在同一碳源浓度下，半乳糖+AHG混合碳源比葡萄糖碳源的Rh1合成反应通量更大。

2.3.4 分析不同碳源策略下，代谢通量分布差异

在$ \alpha $= 0.55 ，碳源浓度为 1 mmol/gDW /h 下，最大化目标函数 $ Z $ 条件，计算代谢的通量分布图（中心代谢，未显示外源反应）

2.3.5 软件设计

另外，在过程中我们认为，使用命令行代码的计算方式并不够用户友好。调研后我们发现，市面上已有的的FBA软件都未提供可视化操作界面，因此我们开发了一款具有可视化界面，并且提供智能机器人辅助操作计算的软件FBA ANALYSIS PLATFORM，无需撰写代码，即可使用FBA算法获得预期结果。

点击此处详细了解

3.3.1 使用的工具

• ColabFold v1.5.5：在 Google Colaboratory（Colab）平台上运行的 Jupyter Notebook，将 MMseqs2的快速同源性搜索与AlphaFold2相结合，进行蛋白质预测。
• AutoDock Vina 1.2.0 [6]：一个用于进行分子对接的开源程序。它最初由斯克里普斯研究所分子图形实验室（现为CCSB）的Oleg Trott 博士设计和实现。
• ChimeraX 1.10 [7]：提供了优秀的分子可视化和分析功能。
• Avogadro [8]：一个分子编辑器，也是一款用于制备对接研究配体的可视化工具。它允许构建和编辑分子结构、优化配体几何构型，并导出与对接软件兼容格式的文件。它被用来优化配体的几何构型并对其施加Gastieger电荷。
• Meeko：负责为 AutoDock 准备输入并处理其输出。它与 AutoDock-GPU 和 AutoDock-Vina 共同开发。可以参数化小有机分子（配体）以及蛋白质和核酸（受体）。
• Consurf服务器 [9]：一个生物信息学工具，用于基于同源序列之间的系统发育关系来评估蛋白质/DNA/RNA 分子中氨基酸/核酸位置的进化保守性。

3.3.2 对接过程

受体结构预测

1. 经过文献查找，我们获取五种酶的氨基酸序列
2. 将获取的序列输入到ColabFold中进行结构预测，使用USCF ChimeraX进行结果的可视化。根据置信度进行颜色编码，同时还有predicted aligned error PAE图，代表预测结构中任意两个残基之间距离的不确定性，同样也是越蓝预测误差越小，两个残基之间相对位置预测得越好。

点击此处获得图片下载链接

配体制作准备

1. 从PubChem下载CID为71571511的琼脂单分子sdf文件
2. 使用Avogadro编辑优化为四聚体和二聚体的结构文件，同时使其能量最小化，下图为两个配体最小化过程。

   点击此处获得图片下载链接

   

   Figure 1(left) FBA原理图

   Figure 3(right) Yeast-GEM9示意图

3. 其中四聚体使用WGLTooLs进行pdbqt格式转化，二聚体使用Meeko进行

计算活性位点

1. 在RCSB蛋白质数据库中搜索近似的实验确定的蛋白质3D结构，从中寻找可以作为参考的活性位点。[29]
   

   图为Aga3463与pdb_0005u1在ChimeraX中进行结构匹配的结果

2. Consurf 的主要用途是识别功能重要的氨基酸残基，映射到3D结构上以识别潜在的功能位点。进化上高度保守的区域，往往能反映该位点在维持蛋白质结构和/或功能中的重要性，即 催化活性、配体结合等。[10] 我们使用Consurf服务器对氨基酸残基进行识别，通过选择高度保守的区域，使用ChimeraX进行取点质心计算，一次获得至少三个可用的对接盒子中心坐标.
   

   图为Aga3463在Consurf里预测的输出结果

3. 计算以上对接参数，获得config.txt（这里以Aga3463为例），完整参数文件请点击此处下载

   config contents

   

                             
   %%=====config_01=====%%
   receptor = /mnt/d/Season_04/iGEM/proteinBuilding/thirdBuild/Aga3463/Aga3463.pdbqt
   ligand = /mnt/d/Season_04/iGEM/proteinBuilding/thirdBuild/Aga3463/agar.pdbqt
   center_x = -46.04
   center_y = -5.99
   center_z = 45.07
   size_x = 50
   size_y = 30
   size_z = 30
   exhaustiveness = 8
   num_modes = 9
   energy_range = 3
   out = /mnt/d/Season_04/iGEM/proteinBuilding/thirdBuild/Aga3463/output_01.pdbqt
   
   %%=====config_02=====%%
   receptor = /mnt/d/Season_04/iGEM/proteinBuilding/thirdBuild/Aga3463/Aga3463.pdbqt
   ligand = /mnt/d/Season_04/iGEM/proteinBuilding/thirdBuild/Aga3463/agar.pdbqt
   center_x = -43.65
   center_y = -4.61
   center_z = 38.18
   size_x = 50
   size_y = 30
   size_z = 30
   exhaustiveness = 8
   num_modes = 9
   energy_range = 3
   out = /mnt/d/Season_04/iGEM/proteinBuilding/thirdBuild/Aga3463/output_02.pdbqt
   
   %%=====config_03=====%%
   receptor = /mnt/d/Season_04/iGEM/proteinBuilding/thirdBuild/Aga3463/Aga3463.pdbqt
   ligand = /mnt/d/Season_04/iGEM/proteinBuilding/thirdBuild/Aga3463/agar.pdbqt
   center_x = 33.08
   center_y = 2.16
   center_z = -28.34
   size_x = 50
   size_y = 30
   size_z = 30
   exhaustiveness = 8
   num_modes = 9
   energy_range = 3
   out = /mnt/d/Season_04/iGEM/proteinBuilding/thirdBuild/Aga3463/output_03.pdbqt
   
   %%=====config_04=====%%
   receptor = /mnt/d/Season_04/iGEM/proteinBuilding/thirdBuild/Aga3463/Aga3463.pdbqt
   ligand = /mnt/d/Season_04/iGEM/proteinBuilding/thirdBuild/Aga3463/agar.pdbqt
   center_x = -45.38
   center_y = -6.85
   center_z = 46.32
   size_x = 50
   size_y = 30
   size_z = 30
   exhaustiveness = 8
   num_modes = 9
   energy_range = 3
   out = /mnt/d/Season_04/iGEM/proteinBuilding/thirdBuild/Aga3463/output_04.pdbqt
   
   %%=====config_05=====%%
   receptor = /mnt/d/Season_04/iGEM/proteinBuilding/thirdBuild/Aga3463/Aga3463.pdbqt
   ligand = /mnt/d/Season_04/iGEM/proteinBuilding/thirdBuild/Aga3463/agar.pdbqt
   center_x = 35.03
   center_y = 3.57
   center_z = -26.91
   size_x = 50
   size_y = 30
   size_z = 30
   exhaustiveness = 8
   num_modes = 9
   energy_range = 3
   out = /mnt/d/Season_04/iGEM/proteinBuilding/thirdBuild/Aga3463/output_05.pdbqt
   
   
   
     

进行分子对接

1. 每一种酶都进行了至少三种不同对接参数的模拟，计算结果同样会是多个。我们对于每一个酶都单独选出其最佳的结合亲和力，接着再与其他酶的结果进行横向对比。

   酶类型	名称	结合能 kcal/mol
   		位点1	位点2	位点3	位点4	位点5
   琼胶酶	Aga3464	-6.23	-7.205	-7.766	-6.077	-7.413
   	AqAga	-8.275	-10.054	-7.361	-6.079	-7.636
   	PdAgaC	-7.396	-7.440	-8.574	-6.722	-8.647
   新琼二糖水解酶	agaNash	-6.686	-5.065	-7.266	-7.669	-
   	NH852	3.725	-4.691	-7.297	-6.736	-

2. 我们选择出在分子对接中的高分结果，即最佳的酶组合，将配体和受体一同在ChimeraX软件中展示：点击此处获取预测结果/受体-配体/参数文件
   

   Figure 1(left) AqAga

   Figure 3(right) agaNash

3. 我们对以上结果进行分析：

   1 AqAga-四聚体

   1.1 氢键计算

   氢键作为分子间相互作用的关键类型，是结合自由能（ΔG）的重要贡献者，对蛋白质-配体复合物的稳定性与功能行使具有决定性作用。为深入理解AqAga与琼胶四聚体之间的特异性识别机制，我们首先对该复合物体系中可能形成的氢键进行了系统分析。

   氢键计算通过ChimeraX软件内置的“H-Bonds”工具完成。具体操作如下：在菜单栏中依次选择“Tools” → “Structure Analysis” → “H-Bonds”，以打开氢键分析面板。为准确捕捉潜在的氢键相互作用，我们设置了以下几何判据：键长距离容差（Distance tolerance）为 0.4 Å，键角角度容差（Angle tolerance）为 20°。

   

   AqAga和四聚体形成的氢键

   分析结果显示，配体分子AqAga与蛋白质四聚体在多个关键残基位点（包括 TRP133、ASP137、THR150、TRP155、LYS177、CYS307、THR622、VAL631 和 TYR679）之间形成了9个氢键，构成了一个密集的相互作用网络。氢键的典型作用距离范围在2.5–3.5 Å之间，而该配体通过与上述多个重要残基形成氢键，被稳固地“锚定”在蛋白质的活性位点区域。

   在这些氢键中，有8个键的键长均小于3.5 Å，显示出较强的作用力，进一步印证该区域极有可能为关键的活性口袋。这些短距离氢键共同构建了一个协调而稳定的结合界面，显著增强了配体与受体之间的亲和力。

   值得注意的是，307Cys所参与的氢键键长为3.728 Å，略高于典型强氢键的范围，并且其空间位置偏离活性口袋的核心区域，仅与配体分子末端结合。因此，该氢键在实际结合中的贡献可能相对有限，其功能性作用较弱，甚至可能不参与关键结合过程。

   1.2 相互作用分析

   除了氢键外，蛋白质配体复合物中还存在其他相互作用，例如范德华相互作用（van der Waals interactions, vdW）等。由于大多数相互作用主要和距离相关，此时可以使用“Contacts”工具分析所有的原子接触，来间接推测可能存在的相互作用。使用方法也基本与氢键分析类似，在菜单栏点击“Tools... Structure Analysis... Contacts”打开“Contacts”面板；

   

   AqAga和四聚体的相互作用

   AqAga和四聚体形成总共129 个contacts，其中大多数分布在氢键分析中推测的活性口袋区域。分析显示，配体的碳原子与VAL 631侧链的两个末端甲基（CG2）分别产生了1.026 Å和0.522 Å的严重空间重叠，同时与PHE 663苯环的中心原子（CZ）也存在0.477 Å的紧密接触。这一系列强烈的疏水区域冲突表明，配体与该口袋的契合度过高，导致了显著的范德华排斥。

   与此同时，配体与GLU 677的侧链存复杂的相互作用：其氧原子（OE1）与配体形成了一个距离仅为2.671 Å的强氢键，而其另一个氧原子（OE2）却与配体碳原子存在0.574 Å的冲突。这种“既吸引又排斥”的局面强烈暗示GLU 677的侧链构象对于配体结合至关重要，其羧基的精确取向决定了它与配体发生的是有利的氢键还是不利的原子碰撞。此外，配体一个氢原子与GLU 468的羧基碳（CD）存在0.402 Å的冲突，这是一种不寻常且不利的相互作用，可能表明该区域是降低配体活性的一个关键位点，需要通过调整配体相应基团的朝向或化学结构来避免。

   综上所述，该结合位点是一个疏水作用与极性作用交织的复杂环境。配体通过强疏水效应被锚定在VAL631/PHE663口袋中，但付出了局部位阻过大的代价；其与GLU677能否形成稳定氢键则高度依赖于对方侧链的构象。未来的优化策略应着重于适度减小配体疏水核心的体积以缓解冲突，并精细调整配体与GLU677间的极性相互作用，从而在减少不利冲突的同时巩固有利的氢键，实现结合亲和力与选择性的提升。

   1.3 口袋分析

   底物通道分析是一种计算和结构生物学方法，用于研究大型生物分子（如酶或蛋白质复合物）内部或之间的可能路径，这些路径允许小分子底物、产物或离子从蛋白质表面运输到其内部的活性位点，或在不同活性位点之间穿梭。简单来说，它就是寻找蛋白质内部的“隧道”或“通道”。

   我们使用CAVER Analyst 2.0来进行底物通道分析。这是一款专门用于计算、分析、可视化和比较蛋白质（或其他生物大分子）中隧道和通道的桌面软件。它是广受欢迎，提供了一个用户友好的界面，使得原本需要命令行操作的分析变得直观易懂。

   

   Figure 1(left) AqAga的表面可及性

   Figure 3(right) Caver 底物通道计算

   我们评估了AqAga的表面可及性，结果显示琼胶四聚体结合于酶的内部空腔中。进一步的底物通道分析表明，该四聚体底物能够贯穿AqAga的整个活性通道（红色区域），提示其结合口袋具有显著的长度与深度，尽管部分通道区域的预测仍存在一定偏差。

   

   活性口袋疏水性分析

   利用ChimeraX计算的分子亲脂性势（Molecular Lipophilicity Potential, MLP）显示，四聚体与蛋白质的结合区域呈现疏水区（黄色）与亲水区（蓝色）交替分布的格局，但整体以亲水氨基酸为主导。这一特征与典型多糖降解酶的活性口袋特征一致：尽管催化中心以亲水环境为主（便于结合亲水性多糖底物及容纳酸性残基如GLU和ASP作为质子供体参与糖苷键的断裂），其周围仍分布着若干关键疏水区域，共同维持完整的催化功能。基于MLP结果和空间结构特征，我们推测右侧大范围的亲水区域可能为底物通道的入口，该区域不仅具有高度亲水性，而且直接暴露于溶剂环境，有利于底物的初始结合与进入。

   2 AqNABH-二聚体

   2.1 氢键计算

   

   AqNABH和二聚体形成的氢键

   对配体分子AqNABH与蛋白质二聚体之间的相互作用进行了分析。结果表明，二者仅在蛋白质亚基的ILE 424残基主链羰基氧处形成一个分子间氢键，其供体-受体距离（D..A dist）为2.863 Å，氢-受体距离（D-H..A dist）为2.132 Å。该键长略大于标准氢键范围，表明此氢键强度较弱，但仍处于可识别范围内。从整体结合模式来看，仅依靠单个弱氢键维持的相互作用体系可能处于一种结合亲和力较低、稳定性欠佳的状态，这提示该复合物的结合可能还依赖于其他非键相互作用（如疏水作用、范德华力等）作为主要稳定因素。

   2.2 相互作用分析

   

   AqNABH和二聚体的相互作用

   AqNABH与二聚体之间共形成45个接触位点，涉及关键残基包括ILE 227、VAL 229、PHE 228、MET 308、TYR 310、LYS 311、ASN 367、ASN 369、ASP 370、ARG 420和ILE 424。其中，ARG 420的CD原子与配体之间的相互作用最为显著（overlap = 0.399 Å，距离 = 3.361 Å），同时配体氢原子与ASP 370主链CA原子之间也存在明显空间重叠（overlap = 0.307 Å，距离 = 2.573 Å）。这类冲突表明，在当前结合构象下，配体与上述残基之间发生了显著的空间排斥，可能不利于复合物的结构稳定性。

   除上述严重冲突外，部分接触位点（如配体与LYS 311的CA原子、ARG 420的CG原子之间）表现出接近零或略负的overlap值，属于紧密但处于边缘状态的接触。这类接触在几何构象发生微小扰动时可能转变为冲突。

   另一方面，多个疏水性和芳香族残基——如PHE 228、TYR 310、VAL 229和ILE 424——与配体之间形成了较为理想的负overlap接触（范围约 -0.1 Å 至 -0.35 Å），表明配体在结合口袋中实现了有效的疏水匹配，这些相互作用可能对复合物的结合稳定性产生积极贡献。

   2.3 口袋分析

   

   AqAga的表面可及性;AqAga的通道推测;AqAga的口袋亲水性分析

   琼胶二聚体稳定镶嵌于AqAga蛋白表面的活性口袋中。根据结构特征与催化残基的空间排布推测，底物分子可能首先通过由疏水残基ILE 424等构成的区域进入活性口袋。随后，位于催化中心的酸性氨基酸残基GLU 422发挥催化功能，特异性地切割琼胶二聚体分子中的糖苷键，使其发生水解断裂。催化反应生成的小分子产物，包括半乳糖和3,6-脱水-L-半乳糖，随后可能经由蛋白表面通道（其出口可能与ILE 153等残基相邻）释放至溶剂环境中，完成整个催化循环。

   进一步对AqAga活性口袋的理化环境进行亲水性分析表明，该口袋主要由亲水性氨基酸残基所构成。这些残基在空间上形成了高度极性的微环境，有助于稳定底物及反应中间体中带电或极性基团的取向。此外，亲水残基的集中分布可能促进了水分子在催化过程中的参与，不仅为水解反应提供必要条件，也有利于产物分子的溶剂化与后续释放。

3.3.3 结果分析

之后，我们开展了系统的组合酶解实验，验证计算模型。结果显示，琼胶酶活力测定实验中，AqAga的活性相对其他酶较强；在新琼二糖水解酶活力测定实验中，AqNABH的活性相对其他酶较强（见下6图）。

通过对比预测结果与实验数据，我们发现无论是琼胶酶活性还是新琼二糖水解酶活性，Sq- Ag5菌株的总酶活均达到最高。这一发现不仅证实了计算模型的准确性，也为琼胶高效生物降解提供了最优的酶组合方案。最终我们选择该酶组合菌株用于初步发酵实验.

4.3.1 符号约定

符号	单位	说明
$S_{1}(t)$	mmol/L	胞外葡萄糖浓度
$S_{2}(t)$	mmol/L	胞外琼胶浓度
$\lambda_i$	\	第i种组合型启动子组合的相对表达强度
$C_{i,j}$	\	第j种诱导性启动子的第i准则
$Q_{i}(t)$	mmol/L	第 i 种启动子组合的角鲨烯产量

其他符号，在后文提及时说明。

4.3.2 条件假设

1. 优先消耗易利用碳源葡萄糖，48h时接近耗尽；
2. 3,6 - 脱水 - L - 半乳糖完全不积累于胞内，且不影响其他碳源代谢，仅作为琼胶分解的 “无效产物”，无需额外监测。

4.3.3 参数准备

为表征各启动子表达水平，我们通过收集湿实验使用的PTEF、GPD、PTDH启动子的不同特征，从而量化每个启动子的相对表达强度。

首先，收集衡量启动子表达强度指标mRNA[11]、GFP水平[12]的测量值来计算三者的相对强度比例，并进行归一化，作为权重特征。以PTEF为基准，计算GPD的GFP等效强度值： $$ GPD_{GFP,equ} = TEF1_{GFP} \times \frac{GPD_{mRNA}}{TEF1_{mRNA}} \tag{1}$$

计算总强度，具体公式如下： $$ S_{total} = S_{PTEF} + S_{GPD} + S_{PTDH} = PTEF_{GFP} + GPD_{GFP,equ} + TDH3_{GFP} \tag{2}$$

计算归一化权重为： $$ W_j = \frac{S_j}{S_{total}}, \,\, j\in \left\{PTEF,GPD,PTDH \right\} \tag{3}$$

代入具体数值，计算得到最终权重分配如表1

启动子	相对表达强度取值
PTEF	0.265
GPD	0.418
PTDH	0.317

对于启动子组合，可以分为同启动子组合与异启动子组合两类，分别为其定义协同增强系数$\alpha$与平均协调系数$\beta$，同启动子组合相对表达强度公式如下： $$ \lambda = \alpha \times \lambda_i \tag{4}$$

异启动子组合相对表达强度公式如下： $$ \lambda = \beta \times \frac{\lambda_i + \lambda_j}{2},i \ne j \tag{5}$$

4.3.4 基于启动子调控的角鲨烯合成通量模型

基于生物反应动力学，对碳源代谢路径与角鲨烯合成的关键环节进行拆解划分，构建常微分方程组。本模型以25g/L 琼胶与 10g/L 葡萄糖混合体系为碳源，其代谢网络可拆解为碳源分解、中间产物分流及角鲨烯合成三部分。

1. 碳源分解动力学模型

基于微生物碳源代谢的经典动力学Monod方程[19]，由于葡萄糖消耗速率与自身浓度正相关，同时当细胞量过低，尽管较大也可能因分解酶总量不足受到限制，引入细胞半饱和常数,当细胞量远大于时，细胞密度调节系数趋近于1。基于葡萄糖优先供合成的假设及需求分析，在满足生长发育的前提下以最大化角鲨烯产量为目标，因此引入优先级权重参数，量化合成对于碳源消耗更强的拉力作用，具体公式如下。 $$ {\frac{d S_{1}}{d t}}=-\ k_{1}S_{1}\cdot{\frac{X(t)}{K_{x}+X(t)}}\left[r(t),{\frac{\mu(t)}{\mu_{m a x}}}+1.4.(1-r(t))^{\cdot}{\frac{Q_{t}(t)}{Q_{max}}}\right]\quad(0\leq t\leq48h) \tag{6}$$ $$ S_1(t) = 0 \,\,(t>48h) \tag{7}$$

前期，琼胶酶分解速度较缓，主要用于生长发育营养的提供，后期葡萄糖耗尽，主要用于角鲨烯的合成。与葡萄糖消耗方程类似，由于分配不同，调整合成优先级系数，具体公式如下。 $$ \frac{dS_2}{dt} = \begin{cases} -k_2 \cdot S_2 \cdot \frac{X(t)}{K_x + X(t)} \cdot \left[ r(t) \cdot \frac{\mu(t)}{\mu_{\max}} + 0.8 \cdot (1 - r(t)) \cdot \frac{Q_i(t)}{Q_{\max}} \right], & 0 \leq t \leq 48h \\ -k_2 \cdot S_2 \cdot \frac{X(t)}{K_x + X(t)} \cdot 1.2 \cdot (1 - r(t)) \cdot \frac{Q_i(t)}{Q_{\max}}, & t > 48h \end{cases} \tag{8}$$

接下来，建立动态碳源分配机制，以角鲨烯产量最大化为核心，通过基础下降参数，保障前期细胞生长发育营养充足安全，引入细胞密度系数，实时反馈避免细胞足够但仍提供大量营养，具体碳源分配比例方程如下。 $$r(t)=\operatorname*{max}\left\{0.55,{\mathrm{min}}\left\{0.75,{\mathrm{ro-}}\,k_{r,i}\cdot{\mathrm{t}}-0.1\cdot{\mathrm{tanh}}\left(2\cdot\left({\frac{X(t)}{X_{m a x}}}-0.6\right)\right)\right\}\right\} \tag{9} $$

2. 细胞生长动力学方程

通过观察实验数据，角鲨烯产量在72h前后增长速度快，因此，我们在72h后引入生长抑制系数，随时间而增强，确保碳源在生长发育充足时优先流向合成角鲨烯通道。具体公式如下。 $$ {\frac{d X}{d t}}=Y\cdot r(t)\cdot\left(-{\frac{d S_{1}}{d t}}-{\frac{d S_{2}}{d t}}\right)\cdot\left(1-{\frac{X(t)}{X_{m a x}}}\right)\cdot\exp\left(-0.2\cdot\left({\frac{t}{72}}-1\right)^{2}\right) \tag{10}$$

3. 角鲨烯通量合成单目标模型

1. 3.1 合成方程的建立

   在启动子调控下，联合碳源供给分配、生长发育，构建角鲨烯合成方程。从酶促与启动子调控上，启动子相对表达强度越高，有效调控能力越强，结合细胞密度X(t)，与酶总表达量成正相关关系，为合成通道提供基础。其次，结合3.1、3.2微分方程，考虑合成过程中生长发育与合成的竞争平衡。当合成碳源比例越高时，碳源向角鲨烯的转换效率越强，代谢流更为集中，副反应减少。具体方程如下。 $$\frac{dQ_i}{dt} = k_q \cdot \lambda_i \cdot X(t) \cdot (1 - r(t)) \cdot \left(0.7 \cdot \left(-\frac{dS_1}{dt}\right) + \left(-\frac{dS_2}{dt}\right)\right) \cdot \eta_{q,i}(t) \tag{11}$$

   其中，$k_q$为MVA途径酶催化效率，$\lambda_i$为启动子相对表达强度，$n_{q,i}(t)$为碳源转换效率。

2. 3.2目标函数及约束条件

   根据实验数据，96h后产量增速极缓，因此可设96h为合成产量的最终时间，以最大化该时刻为目标，记为$t_{end}$。

   为避免发酵后期琼胶过早耗尽，角鲨烯合成将因缺乏物质基础而停滞，导致角鲨烯产量无法持续增长到 96 h，对于碳源剩余设置约束可保障整个发酵周期内碳源供给的稳定性。其次，设置36 h后角鲨烯合成速率需不低于0.001mmol/(L·h)，防止合成长时间中断情况，保障角鲨烯合成的连续性。建立单目标约束模型如下。 $$ \max_{r(t)} Q_i(t_{\text{end}}) \tag{12}$$ $$ \left\{ \begin{array}{l} 0.65 \leq r(t) \leq 0.75 \quad (\forall t \in [0,96]) \\ X(t) \geq 0.4, X_{\max} = 2.0 \quad (\forall t \in [48,96]) \\ S_2(t_{\text{end}}) \geq 0.4 \\ \frac{dQ_i}{dt} \geq 0.001 \quad (\forall t \in [36,96]) \end{array} \right. \tag{13}$$

4. 基于龙格-库塔算法 +有约束非线性优化的求解

根据3的参数分析与初始条件，建立常微分方程组如下： $$\frac{dy}{dt} = \begin{bmatrix} \frac{dS_1}{dt} \\ \frac{dS_{12}}{dt} \\ \frac{dX}{dt} \\ \frac{dQ_i}{dt} \end{bmatrix} = f(t, y) \tag{14}$$

其中，$t \in \left[t_{start},t_{end}\right]=[0.96]$，初始条件为$y_0=y(0)=\left[10.25,0.1,0 \right]^T$

使用MATLAB的ode45求解，四阶龙格-库塔方法[16]提供候选解，五阶方法控制误差的变步长数值解法。 $$\left\{ \begin{aligned} k_1 &= h \cdot f(t_n, y_n) \\ k_2 &= h \cdot f\left(t_n + \frac{h}{2}, y_n + \frac{k_1}{2}\right) \\ k_3 &= h \cdot f\left(t_n + \frac{h}{2}, y_n + \frac{k_2}{2}\right) \\ k_4 &= h \cdot f(t_n + h, y_n + k_3) \\ y_{n+1} &= y_n + \frac{1}{6}(k_1 + 2k_2 + 2k_3 + k_4) \\ t_{n+1} &= t_n + h \end{aligned} \right. \tag{15}$$

其次，设置参数边界约束，建立有约束非线性误差优化模型如下： $$\min_{p} f(p) = \sum_{i=1}^{6} \sum_{j=1}^{4} w_j \left[ Q_{\text{pred}}^{(i)}(t_j, p) - Q_{\text{exp}}^{(i)}(t_j) \right]^2 \tag{16}$$ $$\left\{ \begin{array}{l} lb_k \leq p_k \leq ub_k, \quad k = 1, 2, \ldots, 7 \\ g_j(\mathbf{p}) \leq 0, \quad j = 1, 2, \ldots, m \\ h_k(\mathbf{p}) = 0, \quad k = 1, 2, \ldots, n \end{array} \right. \tag{17}$$

其中，$p=[p_1,p_2,...,p_{\varGamma}]^T$为待优化参数向量，$w_j$为时间点权重，分别对应 24h, 48h, 72h, 96h 时刻。$lb_k,ub_k$分别为参数的下界和上界，$g_j(p)$为不等式约束函数，$h_k(p)$为等式约束函数。

根据“PTDH+PTDH”,“PGPD+PTDH”,“PTEF+PGPD”,“PTEF+PTEF”,“PGPD+PGPD”,“PTDH+PGPD”这六组的实验数据，如表所示，进行训练，优化模型初步设定的参数值。

启动子组合	24h-平均	48h-平均	72h-平均	96h-平均
PTEF+PTEF	92.8400	187.3833	286.7200	365.1000
PTEF+PGPD	92.7467	173.3533	301.0733	363.7467
PGPD+PGPD	57.4267	117.3567	229.1567	326.0167
PGPD+PTDH	72.8400	153.6000	296.9133	325.0200
PTDH+PTDH	79.3933	172.8233	347.7667	424.4700
PTDH+PGPD	97.6783	211.3600	337.9067	439.2200

考虑到上述模型为中等规模的有约束非线性优化问题，采用序列二次规划(SQP)算法，并进行收敛性及稳定性分析。后续可以根据问题规模的扩大优化改进算法。

首先，设拉格朗日函数如下： $$ \mathcal{L}(p, \lambda, \mu) = f(p) + \sum_{j=1}^{m} \lambda_j g_j(p) + \sum_{k=1}^{n} \mu_k h_k(p) \tag{18}$$

必须满足KKT条件，包含平稳性条件、原始可行性及对偶可行性条件，具体如下：

其中，$p^*$为最优解，$\lambda ^*=[\lambda_{1}^{*},\lambda_{2}^{*},...,\lambda_{m}^{*}]^T$为不等式约束的拉格朗日乘子，$\lambda ^*=[\mu_{1}^{*},\mu_{2}^{*},...,\mu_{n}^{*}]^T$为等式约束的拉格朗日乘子。

5. 基于ODE及误差优化的启动子预测

对PTEF、GPD、PTDH这三个启动子的剩余三种可能分组情况进行预测，并对这三种情况进行实验测试。基于我们的研究目标，即为湿实验探寻最佳启动子组合，减少盲目穷举实验，并非完全取代湿实验。因此，只需重点关注角鲨烯最终产量的相对排名，就能够起到为实验探路的作用，对比预测与实验的角鲨烯产量排名，可以看出产量由大到小依次为：“PTEF+PTDH”、“PGPD+PTEF”、“PTDH+PTEF”，大小排比与实验保持一致。将预测结果与实验结果进行对比如下图。

启动子组合	24h-平均	48h-平均	72h-平均	96h-平均
PTEF+PTDH	93.23007	187.3833	299.0133	365.1000
PGPD+PTEF	62.5100	173.3533	302.3167	363.7467
PTDH+PTEF	52.5500	117.3567	237.4233	326.0167

6. 模型评估

计算决定系数R²，同时进行 Kolmogorov - Smirnov（K - S）检验，比较模型结果与实验数据的差异，评估模型准确性。

三组测试启动子组合的整体决定系数R²为0.7987，其中，“PTDH+PTEF”组合高达0.895，说明模型具有较强的可解释能力，能够有效刻画启动子组合对碳源烯产量的调控规律，后续可以进一步通过增加测试数据优化参数，具体如下表。

启动子组合	$R^2$
PTEF+PTDH	0.6348
PGPD+PTEF	0.8373
PTDH+PTEF	0.8945
整体	0.7987

绘制相关性散点图如下，可以看出数据点分布在1:1线附近，落在95%的置信椭圆内，即有95%概率能包含所有“实验值-预测值”数据点的区域，椭圆的中心整体略低，出于模型全面考虑可能存在较多的约束造成系统性的低估属于正常范围。

K-S 检验用于验证模型预测数据与实验数据是否来自同一分布，是评估模型统计一致性的关键指标。我们进一步通过KS检验对预测数据与实验数据的分布一致性进行验证，结果显示显著性水平p值为0.1862（p > 0.05），接受原假设。这表明在5%显著性水平下，模型预测数据与实验数据来自同一分布，避免了模型仅在个别点拟合良好而整体偏离的情况，验证了模型具有良好的统计一致性和可靠性。

4.3.5 诱导性启动子优化模型

在上述基于组合型启动子的模型基础上，我们进一步引入诱导型启动子，以更真实地模拟实际生物系统，进一步拓展了模型的应用范围。诱导型启动子（如GAL1、GAL7、GAL10等）不同于组合型启动子，其表达强度受外界诱导物的浓度调控，这一特性有可能影响协调细胞生长与产物合成之间的平衡，下面选取 GAL1[13]、GAL7[14]、GAL10[15]这三个启动子进行分析。

1. 基于AHP的诱导性启动子评价模型

首先，我们使用层次分析法（AHP）构建诱导性启动子评价模型，通过查找文献及实验分析获取启动子的关键性能参数，建立目标层、准则层和方案层如下。

评价体系	说明
目标层	选择最佳启动子
准则层$C_j$	$C_1$表达强度
	$C_2$诱导倍数
	$C_3$调控特点
	$C_4$系统稳定性
方案层	待评估启动子：GAL1、 GAL7、GAL10

通过1-9标度法两两比较的方式，构建准则层判断矩阵A与权重计算，确定每个准则对目标的相对重要性权重，及每个方案在每个准则下的相对优劣得分。通过计算该判断矩阵的特征向量，得到各准则的权重。计算一致性比例CR < 0.1，通过检验，证明判断矩阵合理。

\	$C_1$	$C_2$	$C_3$	$C_4$
$C_1$	1	3	5	7
$C_2$	$\frac{1}{3}$	1	3	5
$C_3$	$\frac{1}{5}$	$\frac{1}{3}$	1	3
$C_4$	$\frac{1}{7}$	$\frac{1}{5}$	$\frac{1}{3}$	1

接下来，分别针对每个准则$C_j$，我们对启动子进行两两比较，构造判断矩阵及计算相对得分$s_{ij}$。对于每一个启动子i，其总得分$S_i$计算公式如下，计算结果如表5。计算得到各准则得分及综合得分为角鲨烯通道合成模型的优化提供准备。 $$ S_i = \sum_{j=1}^{4} (w_j \times s_{ij}) \tag{20}$$

启动子	$C_1$	$C_2$	$C_3$	$C_4$	综合得分
GAL1	0.373	0.090	0.039	0.006	0.502
GAL7	0.135	0.090	0.039	0.014	0.218
GAL10	0.049	0.090	0.039	0.036	0.280

2. 角鲨烯合成通量模型修正

红藻中的琼胶多糖在琼胶酶的作用下分解为新琼二糖，接着在新琼二糖水解酶的作用下，分解为 D - 半乳糖和 3,6 - 脱水 - L - 半乳糖。诱导型启动子不同于组合型启动子，其表达强度受到诱导剂的影响。基于原有模型，将原先固定的相对表达强度更新为动态表达强度，并新增诱导效率因子，根据1初始取值，后续优化

其次，新增半乳糖消耗动力学方程，描述其诱导及碳源功能。 $$ \frac{dS_3}{dt} = -k_3 \cdot S_3 \cdot \frac{X(t)}{K_x + X(t)} \cdot \left[ \alpha \cdot r(t) \cdot \frac{\mu(t)}{\mu_{\max}} + \beta \cdot (1 - r(t)) \cdot \frac{Q_i(t)}{Q_{\max}} \right] \tag{21}$$

其中，$k_3$为半乳糖消耗速率常数, 初始取0.08 $h^{-1[3]}$。

由于部分半乳糖用于诱导而非直接代谢，引入半乳糖碳源贡献系数c。由于转录及翻译延迟，半乳糖添加后需才能观察到启动子活性上升, 引入延迟函数描述这一现象。修正后的合成方程如下。 $$ \frac{dQ_i}{dt} = k_q \cdot \lambda_{\text{inducer}} \cdot \varepsilon \cdot \tau(t) \cdot X(t) \cdot (1 - r(t)) \cdot \left[ 0.7 \left( -\frac{dS_1}{dt} \right) + \left( -\frac{dS_2}{dt} \right) + c \left( -\frac{dS_3}{dt} \right) \right] \cdot \eta_{q,i}(t) \tag{22}$$ $$ \tau(t) = \begin{cases} 0, & t - t_{\text{inducer}} < 3 \\ 1, & t - t_{\text{inducer}} \geq 3 \end{cases} \tag{24}$$

其中，$t_{inducer}$为半乳糖进入反应时间。

此外，基于4.3.4.3，由于诱导型启动子的特征，新增半乳糖残留约束及启动子活性约束，描述当过低时，启动子强度骤降的现象。

3. 求解结果与分析

根据“PGall+PGall、PGall+GAL10、PGal7+PGal7、PGal7+PGall、PGal7+GAL10、GAL10+GAL7”这六组的实验数据，进行训练，优化模型初步设定的参数值。

启动子组合	24h-平均	48h-平均	72h-平均	96h-平均
PGall+PGall	94.5067	191.6800	347.4600	427.7900
PGall+GAL10	126.9300	250.5367	391.5667	505.2600
PGal7+PGal7	131.7667	287.3533	373.1600	477.8167
PGal7+PGall	121.7433	274.8033	354.6333	475.7767
PGal7+GAL10	117.0783	246.8233	427.1033	500.3700
GAL10+GAL7	152.1233	318.4700	469.2700	594.6800

对“PGal1+PGPD” 、“GAL10+PGal1”、“GAL10+GAL10”三种分组情况进行预测，并对这三种情况进行实验测试。将预测结果与实验结果进行对比如下图，可以看出产量由大到小依次为：“PGal1+PGPD”、“GAL10+GAL10”、“GAL10+PGal1”，大小排比与实验保持一致。

启动子组合	24h-平均	48h-平均	72h-平均	96h-平均
PGall+PGPD	113.1333	187.3833	435.4800	549.6800
GAL10+GAL10	118.7550	252.1567	396.1333	481.0667
GAL10+PGall	125.8067	267.4500	411.8667	505.5833

4. 模型评估

计算决定系数R²，同时进行 Kolmogorov - Smirnov（K - S）检验，比较模型结果与实验数据的差异，评估模型准确性。

三组测试启动子组合的整体R²为0.8384，其中，“GAL10+PGal1”高达0.9149，说明模型具有较强的可解释能力，具体如下表。

启动子组合	$R^2$
PGal1+PGPD	0.8598
GAL10+GAL10	0.7406
GAL10+PGal1	0.9149
整体	0.8384

进一步通过KS检验对预测数据与实验数据的分布一致性进行验证，结果显示显著性水平p值为0.7864（p > 0.05），接受原假设。这表明在5%显著性水平下，模型预测数据与实验数据来自同一分布，验证了模型具有良好的统计一致性和可靠性。

4.3.6 结论

围绕人参皂苷 Rh1 红藻生物合成中角鲨烯产量的关键调控因素 —— 启动子组合，我们构建并完善了基于 ODE 动力学的最佳启动子预测模型，通过模型建立、求解、预测与评估，以及诱导性启动子模型的扩展，为角鲨烯合成机制研究提供了量化工具。模型预测能力可靠，在扩展后性能进一步提升，整体 R² 升至 0.8384，K-S 检验 p 值达 0.7864，对 “PGal1+PGPD”等组合的排名预测与实验一致。通过模型预测高产量组合，可指导湿实验开展的优先级排序，显著减少盲目实验次数，降低时间与试剂成本。后续研究可以进一步通过增加更多启动子组合的实验数据、优化协同系数等参数，进一步提升模型的精度。

1 成本核算模型

1.1 传统制备方法

1.1.1 成本核算边界与分类框架

将人参皂苷 Rh1 制备成本划分为固定成本（设施与原料基础投入）、可变成本（生产过程动态支出）、间接成本（废水 / 废渣处理）三大类，由于多种方法工业生产用地成本计算类似，为简化对比，不纳入计算。核算边界覆盖从原料获取到产物分离的全流程，具体如下表。

成本类别	参数符号	核心构成	核算逻辑	数据来源
固定成本	FC	提取设备折旧、配套设施成本等	按使用年限均摊，参考工业用地与设备投资标准	行业工程定额、企业报价、市场监管文件
可变成本	IC	原料采购、辅料消耗、能耗、人工、税收减免	按单位产量动态计算，结合市场价格	市场报价
间接成本	VC	废水 / 废渣处理、原料种植 / 养殖环境影响	参考环境治理市场价	环保报告
总成本	C	FC+IC+VC	\	\

1.1.2 工艺流程

人参皂苷 Rh1 的传统提取工艺主要是从人参及其相关原料中获取或者基于人参等植物原料，常见的有乙醇酸水解法、生物酶法、超声波辅助提取。下面选取“酶转化与金属离子催化联用的二步法”进行具体分析，该方法相比其他方法在产物效率上具有优势，其工艺参数均有清晰实验数据为成本量化核算提科学支撑。

从产地采购的需先经过预处理筛选，鲜参烘干加工后其中人参茎叶可达到人参的40%～48% [20]，后续计算取平均值44%。经提取、精制纯化、浓缩干燥等工艺处理人参茎叶，可获得收率达35.15％[21]的人参皂苷Re。将Re作为底物，使用二步法转化工艺，通过酶催化反应制备Rg1，得率为70.5%[22]，再以 Fe3+为催化剂进行反应,最终得到人参皂苷Rh1组异构体，得率为80.2%[22]。 流程图如下图所示

1.1.3 工艺流程

假设最终需获得1g Rh1，首先计算原料用量，经推导，所需鲜参计算公式如下： $$ n_0 = \frac{1g}{56.64%} \div 35.15\% \div 44\% \approx 11.43g \tag{1}$$

根据反应体系及产物用量[20]，可计算出制备 1g Rh1 需催化剂与试剂的用量，传统人参提取法成本明细如表2所示。

成本大类	分类	参数符号	参数名称	取值	单位	用量	数据来源
可变成本 $VC_{tradition}$	原料成本	$P_{ginseng}$	鲜参成本	0.196	元 / g	11.43	产地价
	加工成本	$C_{process}$	鲜参清洗、烘干、切片	0.007	元 / g	\	康美中药城
	催化剂与试剂成本 $C_{catalyst}$	$c_{1}$	麦芽汁培养基	0.26	元 / g	$n_1$ = 54.64g/g Rh1	网购平台
		$c_{2}$	硫酸铵	0.001	元 / g	$n_2$ = 29.71g/g Rh1	网购平台
		$c_{3}$	AB - 8 大孔吸附树脂	0.316	元 / g	$n_3$ = 10g/g Rh1	网购平台
		$c_{4}$	乙醇	0.005	元 / g	$n_4$ = 21.49g/g Rh1	网购平台
		$c_{5}$	六水合三氯化铁	0.122	元 / g	$n_5$ = 13.88g/g Rh1	网购平台
间接成本 $IC_{traditon}$	环境污染成本	$ c_{waste} $	废水处理	1.186	元/立方米	$n_6$ = 5L/g Rh1	国发改委[24]
固定成本 $FC_{traditon}$	仪器成本	$ S_{1} $	透析袋	180	元/m	1m/ L 酶液	网购平台
		$ S_{2} $	离心机	100000	元/5.8 年	\	网购平台、 化工仪器网
		$ S_{3} $	旋转蒸发仪	60000	元/10 年	\	网购平台、 售后维修公司

根据反应体系及产物用量[20]，制备 1g Rh1 需原料、催化剂与试剂成本计算公式如下： $$ VC_{traditon} = P_{ginseng} \times n_0 + c_1 \times n_1 + c_2 \times n_2 + c_3 \times n_3 + c_4 \times n_4 + c_5 \times n_5 = +21.517 \,\, 元/g \tag{2}$$

间接成本为：$IC_{tradition} = c_{waste} \times n_6 = 0.006$ 元/g

假设单次生产100g Rh1，每年工作时间为250天，单位人参皂苷 Rh1 的固定成本分摊计算公式如下： $$FC_{tradition} = (10\times S_1 + S_2 \div (5.8 \times 250)+ S_3 \div (10 \times 250))\div 100g = 18.930 \,\, 元/g \tag{3}$$

联立（1）-（3）式，计算得到的单位 Rh1 总成本为： $$C_{tradition} = VC_{traditon} + IC_{traditon} +FC_{tradition} = 21.517+0.006+18.930 = 40.453 \,\, 元/g \tag{4}$$

1.2 红藻生物合成法

在我们提出的红藻生物合成制备工艺中，所采用的红藻粉原料富含琼胶成分，其琼胶含量区间为61%-67.3%[21]，下面计算取64.15%。发酵阶段以琼胶作为关键底物配置培养基，设定培养基中琼胶的浓度为25g/L，经检测，发酵液中人参皂苷Rh1的产量达到141.78mg/L。

红藻粉制备人参皂苷 Rh1 的得率计算公式如下： $$ Y_{red-alga} = (141.78 \div \frac{培养基中琼胶质量}{红藻粉琼胶含量} \times 100) \times 100\% \approx 0.364\% \tag{5}$$

1.2.1 实验室成本计算

首先对实验室规模下制备 1g 人参皂苷 Rh1 的成本进行核算，具体成本明细如表3所示。

成本大类	分类	参数符号	参数名称	取值	单位	数据来源
可变成本 $VC_{lab-alga}$	原料成本	$P_{red-alga}$	红藻成本	4.301	元 / L	网购平台
		$P_{C_{6}H_{12}O_{6}}$	葡萄糖成本	0.003	元 / g	葡萄糖 10g/L 网购平台
		$P_{media}$	其他培养基	0.39	元 / L	网购平台
	分离纯化成本	$c_{sep}$	萃取试剂 + 硅胶柱耗材	8	元 / L	网购平台
间接成本$IC_{lab-alga}$	环境污染成本	$c_{waste}$	发酵废水处理	1.186	元/立方米	发改委[24] (总量平均为2L)
固定成本$FC_{lab-alga}$	仪器成本	${S_1}^{\prime}$	实验室发酵罐	80000	元/5年	网购平台
		${S_2}^{\prime}$	实验室离心机	25632	元/5 年	网购平台、化工仪器网
		${S_3}^{\prime}$	实验室旋转蒸发仪	6200	元/10 年	网购平台、售后维修公司

根据实验数据，单次实验可生产 141.78mg/L 人参皂苷 Rh1，可变成本计算公式如下： $$ VC_{lab-alga} = (P_{red-alga} + P_{C_{6}H_{12}O_{6}} \times 10 + P_{media} +C_{sep}) \div 0.14178 = 89.724 \,\, 元/g \tag{6}$$

间接成本计算公式如下： $$IV_{lab-alga} = c_{waste} \times 2 \div 0.14178 = 16.720 \,\, 元/g \tag{7}$$

根据实验数据，单次生产平均耗时为144小时，一般设备额定机时按每年 800 小时计算[25],单位人参皂苷 Rh1 的固定成本分摊计算公式如下： $$FC_{lab-alga} = ({S_3}^{\prime} \div (5\times 800) + {S_2}^{\prime} \div (5\times 800) + {S_2}^{\prime} \div (5\times 800)) \div 144 \div 0.14178 = 1.331\,\, 元/g \tag{8}$$

联立(7)-(9)式，计算得到实验室规模下单位 Rh1 成本为： $$C_{lab-alga} = VC_{lab-alga} + IV_{lab-alga} + FC_{lab-alga} = 104.785 \,\, 元/g \tag{9}$$

1.2.2 工业成本计算

假设采用 100m³ 发酵罐进行工业生产，单次发酵周期为6 天，根据实验数据，单次生产141.78mg/L Rh1，Rh1 产量 = 100m³×141.78mg/L=14.178kg，工业规模下的成本明细如表4所示

成本大类	分类	参数符号	参数名称	取值	单位	数据来源
可变成本 ${VC_{ind-alga}}^{\prime}$	原料成本	${P_{red-alga}}^{\prime}$	红藻成本	1373.626	元 /100m³	网购平台
		${P_{C_{6}H_{12}O_{6}}}^{\prime}$	葡萄糖成本	0.003	元 / g	网购平台 （葡萄糖 10g/L）
		${P_{media}}^{\prime}$	其他培养基	39000.000	元 /100m³	网购平台
	分离纯化成本	${c_{sep}}^{\prime}$	萃取试剂 + 硅胶柱耗材	200000.000	元 /100m³	网购平台
间接成本${IC_{lab-alga}}^{\prime}$	环境污染成本	${c_{waste}}^{\prime}$	发酵废水处理	118.600	元/立方米	发改委[24] (总量平均为2L)
固定成本${FC_{lab-alga}}^{\prime}$	仪器成本	${S_1}^{\prime \prime}$	发酵罐	16800.000	元/8年	网购平台
		${S_2}^{\prime \prime}$	离心机	100000.000	元/5.8 年	网购平台、化工仪器网
		${S_3}^{\prime \prime}$	旋转蒸发仪	60000.000	元/8年	网购平台、售后维修公司

计算可变成本如下: $$ VC_{ind-alga} = ({P_{red-alga}}^{\prime} + {P_{C_{6}H_{12}O_{6}}}^{\prime} \times 10 + {C_{sep}}^{\prime}) \div 14.178 \div 1000 = 16.975 \,\, 元/g \tag{10} $$

间接成本计算公式如下： $$IC_{ind-alga} = {C_{waste}}^{\prime} \times 2 \times 100 = 1.673 \,\,元/g \tag{11} $$

假设每年工作时间为250天，计算分摊得到单位 Rh1 固定成本为： $$FC_{ind-alga}=({S_3}^{\prime\prime} \div (8 \times 250) + {S_2}^{\prime \prime} \div (5.8 \times 250)+ {S_1}^{\prime \prime} \div (8 \times 250)) \div 14.78 \div 1000= 0.008 \,\, 元/g \tag{12}$$

联立(5)式、(10)-(12)式，计算得到实验室中单位 Rh1 成本为： $$C_{lab-alga} = VC_{ind-alga} + IC_{ind-alga} +FC_{ind-alga} = 18.656 \,\,元/g \tag{13}$$

1.3 成本对比与敏感性分析

1.3.1 成本对比

由下图，红藻生物合成法制备应用到工业时，成本显著低于传统人参提取的工业制备成本，展现出良好的成本优势与工业化潜力。

1.3.2 敏感性分析

下面针对原料价格，进行敏感性分析：

1. 若红藻价格上涨 10%，红藻法实验室制备总成本上升 2.894%（至 107.818元 /g），而红藻法工业制备总成本仅上升0.053 %（18.666元 /g）。
2. 鲜参价格上涨 10%，原料成本上涨为0.2156元 /g，传统法总成本上升0.553%（至40.678元 /g），为红藻法工业制备的10.434倍

这表明，红藻生物合成法相较于传统鲜参提取 Rh1更适合工业化生产，因对原料价格波动敏感度低、成本稳定性与抗风险能力更强，能为后续应用于人参皂苷类功效型化妆品规模化量产提供持久性支撑。

因此，红藻生物合成法能有效降低人参皂苷 Rh1 提取的成本且受原料定价波动的生产风险小，在工业化生产中展现出更强的成本控制能力与市场竞争力。在供应链端，它能保障生产成本长期稳定，避免生产断供情况发生，进而确保人参皂苷 Rh1 功效浓度稳定达标，为企业配方优化与研发腾出空间，有力推动人参皂苷在化妆品原料领域的广泛应用。

2 碳足迹量化与碳利用率计算模型

2.1 方法介绍及系统边界范围

基于生命周期法（LCA），对传统人参提取法与红藻生物合成法工业制备人参皂苷 Rh1 的碳足迹和碳利用率进行量化建模。通过识别关键碳排放源、构建数学计算模型，计算产品对全球变暖的潜在影响以及在不同阶段、不同过程、不同空间位置的影响构成 (以二氧化碳当量表示)。对比两种路径的碳排放强度与碳利用效率，为 Rh1制备工艺的低碳化优化及碳中和目标提供科学依据，同时结合对流弥散方程等工具完善废水等间接排放的核算精度。

2.2 模型假设及原则

1. 运输载重假设：依据《道路货物运输及站场管理规定》及行业常见载重标准，假设轻型载货汽车单次额定运输量为 6t；
2. 运输时效假设：参考《机动车运行安全技术条件》（GB 7258）及公路运输行业运营规范，假设柴油轻型货车在干线公路运输中的平均行驶速度为 65km/h；结合驾驶员劳动强度及安全管理要求，设定每连续行驶 6-8h 进行一次休息，用于核算运输过程中的驻车蒸发排放。
3. 取舍准则：在数据收集中若发现个别物质流或能量流对某一单元过程的碳足迹无实质性贡献，则可以将其排除在外，累计排除总量不得超过本次计算单元过程总温室气体排放量的５％[26]；
4. 数据质量要求： 排放因子的选择按以下优先顺序:
1. 实际测量或质量平衡获得的排放因子
2. 供应商提供的排放因子
3. 区域排放因子
4. 国内排放因子
5. 国际排放因子

2.3 基于LCA的碳足迹量化模型

消除不同工艺路线下产品产量差异对碳足迹对比的干扰，统一设定功能单位为“生产 1g 人参皂苷 Rh1”。

生命周期是产品相关的连续且相互连接的阶段包括原材料获取或从自然资源中生成原材料至生命末期处理。由于无论传统制备还是红藻合成Rh1，后续流程一致，均为产品销售、 产品使用和废弃处置。而本模型目的在于对比不同方法下制备Rh1的环境效益，主要考虑系统边界为“从摇篮到大门”的形式，暂不将其纳入量化范围。

2.3.1 传统人参提取法

传统人参合成Rh1生命周期链条可划分为种植阶段、运输阶段与提取纯化阶段。从种植环节来看，受气候、土壤等生长环境条件的严格限制，人参主产区集中于吉林、辽宁、黑龙江三省 [27]，依赖森林土壤与人工施肥。

从产业链企业的区域分布特征来看，上游种植及初加工企业主要聚集于吉林省，形成了明显的产业集群效应；而人参产业中的代表性上市企业，区域分布则相对分散，覆盖东北、中部及南部等多个地区，且不同区域的企业定位存在显著差异：东北地区的企业多聚焦于人参的种植、采收与初级销售，承担产业链源头供给职能；中南部地区的企业则大多以人参深加工及研发为核心业务，需从东北产区采购原料，这就使得人参从种植基地到中南部加工研发环节之间，形成了较长距离的运输链条。依照企业分布，我们模型计算以运输路线：吉林通化（东北主产区）至浙江杭州为例（提取加工）。

由(2)式，生产 1g Rh1 需 11.43g 鲜参 ，单次运输对应 Rh1 产量约 524.934kg。柴油轻型载货汽车（国 4 标准）在基准条件（温度 15℃、湿度 50%、速度 30km/h、负载 50%、柴油硫含量 50ppm）下的污染物排放因子及修正因子[28]如表5。行驶过程蒸发排放因子为11.6 g/h，驻车过程蒸发排放因子为6.5 g / 天[28]。

气体	尾气基准排放因子 EF(g/km)	气象修正因子（φ₁）	速度修正因子（γ₁）	负载修正因子（θ₁）	劣化修正因子（λ₁）
CO	1.48	1.33	0.7	1.23	1.3
HC	0.186	1.33	0.64	1.00	1.3
$NO_x$	2.636	0.88	0.60	1.29	1.3
PM	0.058	1	0.65	1.18	1.3

$$尾气排放量 = EF \times \varphi_1 \times \gamma_1 \times \theta_1 \times \lambda_1 \times 运输距离 \times 10^{-6} \tag{14}$$

根据公式(14)，计算得到各尾气排放量之和为0.0046t。按照GWP换算成总碳排放量为： $$\mathrm { 0.0046t \times }\mathrm {(} \frac {\mathrm{44}} {\mathrm{28}} \mathrm{) + 0.00035t \times 25 + (11.6g/h \times 27h + 6.5g/} \mathrm{天}\mathrm{ \times 0.25 }\mathrm{天}\mathrm{)}\mathrm{\times }\mathrm {10^{-6}} \mathrm{\approx 0.01632 t} \tag{15} $$

故单位 Rh1 运输碳排放：​ $$ E_{运输} = 0.01632t \,\,CO_{2}e \div 524.934 \,\,kg Rh1 \approx 0.0312 \,\, kg \,\, CO_{2}e/kg \,\,Rh1 \tag{16} $$

提取阶段碳排放主要来源于提取设备运行消耗的电能，计算公式为： $$ E_{提取} = 提取阶段总耗电量 \times 区域电力碳排放因子 \tag{17}$$

参考化工仪器网、行业工程定额，根据设备功率与运行时长计算出耗电量代入上式中。不同阶段碳排放计算结果如表6所示。

碳排放阶段	参数符号	结果($kg\,\,CO_{2}e/g$ Rh1)	说明	数据来源
种植	$E_{种植}$	1.850	每 100g 鲜参对应化肥用量0.8g	白山市统计局
运输	$E_{运输}$	0.031	公路运输距离约1755km	参考实际物流干线里程
提取	$E_{提取}$	0.278	东北电力平均排放因子0.5564 kg $CO_{2}$/kWh	生态环境部前瞻网行业报告

根据表6，计算总碳排放为： $$CF_{提取}=E_{种植} +E_{运输} + E_{提取} = 2.159 \,\,kg \,\,CO_{2}e/g \,\,Rh1 \tag{18}$$

2.3.2 红藻生物合成法

红藻生物合成Rh1生命周期链条可划分为养殖采收、运输、发酵阶段。根据Rh1得率可以计算出生成1g Rh1 需红藻原料274.725g。

根据中国渔业统计年鉴数据显示，中国海藻主要产区为福建、山东和辽宁。其中，福建省是我国海藻最大的产区。从产业链空间分布来看，我国合成生物学产业链在浙江、江苏、广东、上海、北京等省市具有较为完善的布局。这就使得红藻从养殖基地到加工研发环节之间的运输链条可控，形成 “集中运输 - 集中加工” 的模式，减少中间转运节点。养殖采收阶段碳排放主要来源于养殖设施（如聚乙烯养殖网箱）生产与采收设备运行耗电，将养殖设施聚乙烯用量、聚乙烯排放因子、采收设备耗电量代入以下式子计算。 $${E_{养殖采收}}^{\prime} = 养殖设施排放量 + 采收设备耗电碳排放量 \tag{19}$$

与2.3.1计算逻辑类似，根据公式(19)，计算得到各尾气排放量之和为0.0011 t, 单次运输对应 Rh1 产量约 21.840 kg。按照GWP换算成总碳排放量为： $$0.0011t \mathrm { \times (} \frac{44}{28}) + 0.00035t \mathrm { \times} 25 + (11.6g/h \mathrm {\times} 27h + 6.5g / \mathrm{天 \times}0.25 \mathrm {) \times } 10^{-6} \approx 0.3161 kg \tag{20}$$

故单位 Rh1 运输碳排放: $${E_{运输}}^{\prime} = 0.3161 + 21.840 = 0.014 \,\, kg\,\,CO_{2}e/g \,\,Rh1 \tag{21}$$

发酵阶段碳排放主要来源于发酵罐运行耗电，计算公式为 $${E_{发酵}}^{\prime} =发酵阶段总耗电量 \times 区域电子碳排放因子 \tag{22}$$

与2.3.1计算逻辑类似，根据设备功率与运行时长计算出耗电量代入上式中。不同阶段碳排放计算结果如表7所示。

碳排放阶段	参数符号	结果($kg \,\,CO_{2}e/g$ Rh1)	说明	数据来源
养殖采收	${E_{养殖采收}}^{\prime}$	0.003	聚乙烯排放因子0.06 kg $CO_{2}/kWh$	环保局
运输	${E_{运输}}^{\prime}$	0.014	假设由福建红藻产区至广东发酵厂，运输距离约400km	参考实际里程
发酵	${E_{发酵}}^{\prime}$	0.006	电力平均排放因子0.5366 kg $CO_{2}/kWh$	生态环境部

计算得到总碳排放为： $$CF_{提取} = E_{种植} + E_{运输} + E_{提取} = 0.023 kg\,\,CO_{2}e/g \,\, Rh1 \tag{23}$$

2.3.3 方法碳足迹对比分析

根据上述计算对比，系统边界各环节碳排放量如图3所示，传统人参提取因种植环节需大量土地、资源投入，并且原产地与后续工艺加工厂分布较为分散，运输路程远，导致碳排放居高不下；而红藻取自海洋，采收等环节碳排放较低，后续加工等环节地理位置集中，整体碳足迹远小于人参提取。这种差异为产业低碳转型指明方向，红藻生产模式在减碳上更具潜力，能助力行业降低对高碳生产方式的依赖，推动绿色供应链构建与低碳生产实践。

2.4 碳足迹利用率计算模型

碳利用率（Carbon Utilization, CU）是衡量制备工艺碳资源利用效率的关键指标，即计算目标产物人参皂苷 Rh1中固定的碳元素质量占制备全生命周期总碳排放量（以 CO₂当量计）的百分比，公式如下： $$CU = \frac{C_{Rh1}}{CF_{总}} \times 100\% \tag{24}$$

人参皂苷 Rh1 的化学分子式为$C_{42}H_{72}O_{14}$,据元素相对原子质量计算1g Rh1 所含碳元素质量为0.63 g 。根据3.3两种制备路径的碳足迹计算结果，$CF_{传统}=2.159 kg\,\,CO_{2}e/g Rh1, CF_{红藻}=0.023 kg\,\,CO_{2}e/g Rh1,$ 代入公式得到： $$CU_{传统} = \frac{0.63}{2159} \times 100\% \approx 0.0292\% \tag{25}$$ $$CU_{红藻} = \frac{0.63}{23} \times 100\% \approx 2.739\% \tag{26}$$

由此可知，红藻生物合成法的碳利用率约为传统人参提取法的 94 倍，表明其在碳资源转化效率上具有优势，相同碳排放量下，红藻法能将更多碳固定到目标产物 Rh1 中，减少碳资源浪费。

两种方法碳利用率差异的主要原因是红藻法的总碳排放量远低于传统法，这与红藻养殖碳排放低、运输距离短、发酵能耗低等特点相关。红藻法更高的碳利用率不仅符合 “双碳” 目标下的产业低碳化要求，还能降低企业在碳交易、碳税等政策下的潜在成本，进一步强化其化妆品量产应用的综合竞争力。

3 红藻生物合成Rh1可持续性分析

以 “1g 人参皂苷 Rh1” 为基准核算单位，构建 “经济成本 - 环境效益” 二元对比模型，聚焦成本构成、碳排放强度、碳利用率三项指标，建模量化对比两种制备路径，差异显著，具体如图4所示。红藻生物合成法在工业应用上无论是成本或是碳排放、碳利用率均优于传统的人参提取法，可降低化妆品原料采购成本，为终端产品定价提供灵活空间，有效推动人参中的稀有人参皂苷Rh1规模量产，同时兼顾符合 “双碳” 目标与产业绿色升级需求。