2.3.1 Constructing External Metabolites and Reactions

Based on the Yeast-GEM9 model, we added exogenous reactions for ginsenoside synthesis and AHG (3,6-anhydro-L-galactose) conversion, including 13 metabolites and 16 reactions, thereby constructing a Saccharomyces cerevisiae metabolic model capable of utilizing glucose or a mixture of galactose and AHG to synthesize Rh1.

Some required substances and metabolic pathways are not naturally present in yeast. We introduced them into the model. The following table shows the new metabolites used.

Below are the newly used reactions.

The following table shows the upper and lower limits of the relevant reaction fluxes used.

2.3.2 Finding the Target Function Weight Ratio

Cellular resources are limited, requiring a trade-off between growth and product synthesis. To avoid having smaller-weight reactions completely dominated by larger-weight ones, we need to find a balanced weight ratio $\alpha$ for the objective function. In FBA, the biomass synthesis and Rh1 synthesis reactions can be combined into a weighted objective function as follows: $$Maximize(Z)=\alpha * Flux(Biomass)+(1-\alpha)*Flux(Rh1) \tag{1}$$

Let $Z$ be the objective function, $\alpha$ be the weight for strain growth (ranging from 0 to 1), and $ Flux(x)$ represent the flux of reaction $ x $. The carbon source concentration is calculated based on $\alpha$ ranging from 0 to 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)
   

The plot of maximizing objective function $ Z $ is shown below:

Figure 3 Identify the equilibrium parameters of the objective function.

The balanced point $ \alpha $ =0.55 is determined as the point where the objective function value is minimized under the given weight settings $ \alpha $. Choosing this weight ensures the objective function remains in a relatively balanced state, avoiding dominance by a single objective (either biomass or product synthesis).

It can be observed that on the left side of the growth rate curve and the right side of the Rh1 curve, due to imbalance in the growth function, changes in $\alpha$ have no effect, as they are more dominated by the other reaction.

2.3.3 Testing Rh1 Synthesis under Different Carbon Source Strategies

With $ \alpha $ = 0.55, under conditions of maximizing the objective function $ Z $, we calculated the Rh1 synthesis reaction rate at different carbon source concentrations. The results are shown in the figure below:

Figure 4 Ginsenoside Rh1 yield across different carbon source strategies.

The graph shows that, at the same carbon source concentration, the Rh1 synthesis flux using the galactose + AHG mixed carbon source is higher than that using glucose.

2.3.4 Analyzing Metabolic Flux Distribution Differences Under Different Carbon Source Strategies

With $ \alpha $= 0.55 and a carbon source concentration of 1 mmol/gDW/h, under conditions of maximizing objective function $ Z $, we calculated the metabolic flux distribution map (central metabolism, exogenous reactions not shown).

Figure 5 Differences in metabolic flux distribution under different carbon source strategies.

2.3.5 Software Design

Additionally, we found that command-line-based computation is not user-friendly enough. After research, we discovered that existing FBA software lacks a visual interface. Therefore, we developed a software tool called FBA ANALYSIS PLATFORM, featuring a graphical user interface and an intelligent robot assistant for guided operation. It allows users to obtain FBA results without writing any code.

[Click here for more details]

Figure 6 FBA analysis platform: graphically interactive agent-based software.

3.3.1 Tools Used

• ColabFold v1.5.5: A Jupyter Notebook [5] running on Google Colaboratory (Colab) that combines MMseqs2's fast homology search with AlphaFold2 for protein structure prediction.
• AutoDock Vina 1.2.0: An open-source program for molecular docking [6], originally designed by Dr. Oleg Trott from the Scripps Research Institute Molecular Graphics Laboratory (now CCSB).
• ChimeraX 1.10: Provides excellent molecular visualization and analysis functions [7].
• Avogadro: A molecular editor used to prepare ligands for docking studies [8]. It allows building and editing molecular structures, optimizing ligand geometry, and exporting files in formats compatible with docking software. It was used to optimize ligand geometry and assign Gasteiger charges.
• Meeko: Prepares input for AutoDock and processes its output. Co-developed with AutoDock-GPU and AutoDock-Vina, it can parameterize small organic molecules (ligands) as well as proteins and nucleic acids (receptors).
• Consurf Server: A bioinformatics tool that evaluates the evolutionary conservation of amino acid/nucleotide positions in proteins/DNA/RNA based on phylogenetic relationships among homologous sequences [9].

3.3.2 Docking Process

1. After literature review, we obtained the amino acid sequences of the five enzymes.
2. The sequences were input into ColabFold for structure prediction, and results were visualized using USCF ChimeraX. Confidence was color-coded, and predicted aligned error (PAE) plots indicated uncertainty in the predicted distance between any two residues—bluer areas indicating lower prediction error and better relative positioning.

[Click here to download image]

(A) Aga3463 's structure and residue confidence scores.

(B) AqAga 's structure and residue confidence scores.

(C) PdAgaC 's structure and residue confidence scores.

(D) agaNash 's structure and residue confidence scores.

(E) AqNABH 's structure and residue confidence scores.

Figure 8 Protein structures and residue confidence scores generated by alphaFold2 colab.

Ligands Preparation

1. Downloaded the SDF file for a single agar molecule (PubChem CID: 71571511).
2. Used Avogadro to edit and optimize the structure into tetramer and dimer forms, minimizing energy. The figure shows the minimization process for the two ligands.

   [Click here to download images]

   

   Figure 9 (Left Panel) Constructing ligand molecules for agarose tetramer in Avogadro software.

   (Right Panel) Constructing ligand molecules for agarose dimer in Avogadro software.

3. The tetramer was converted to PDBQT format using WGLTools; the dimer used Meeko.

Calculating Active Site

1. Searched the RCSB Protein Data Bank for experimentally determined 3D structures of similar proteins to identify reference active sites [29].
   

   Figure 10 Structural alignment of Aga3463 with pdb_0005u1 in ChimeraX.

2. Consurf identifies functionally important amino acid residues by mapping evolutionary conservation onto 3D structures. Highly conserved regions often reflect importance in maintaining protein structure and/or function, such as catalytic activity or ligand binding [10]. We used the Consurf server to identify highly conserved regions and calculated centroid coordinates using ChimeraX, obtaining at least three usable docking box center coordinates.
   

   Figure 11 Consurf prediction output for Aga3463.

3. Using the above docking parameters, we generated config.txt (using Aga3463 as an example).[Click here to download]

   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
   
     

Performing Molecular Docking

1. Each enzyme underwent at least three different docking simulations, yielding multiple results. For each enzyme, we selected the best binding affinity result and then compared them crosswise.

   Table 4 Docking results.

   Enzyme Type	Name	Binding Energy (kcal/mol)
   		Site 1	Site 2	Site 3	Site 4	Site 5
   Agarase	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
   Neoagarobiose Hydrolase	agaNash	-6.686	-5.065	-7.266	-7.669	-
   	NH852	3.725	-4.691	-7.297	-6.736	-

   Notes: A dash ("-") denotes no valid data.

2. We selected high-scoring docking results—the optimal enzyme combination—and visualized the ligand-receptor complexes in ChimeraX: [Click here to access prediction results / receptor-ligand / parameter files]
   

   Figure 12 (Left Panel) The optimal binding mode of AqAga protein and agarose tetramer ligand.

   (Right Panel) The optimal binding mode of agaNash protein and agarose dimer ligand.

3. We analyze the above results:

   1 AqAga--Tetramer

   1.1 Hydrogen Bond Analysis

   Hydrogen bonds are a key intermolecular interaction, significantly contributing to binding free energy (ΔG) and playing a decisive role in the stability and function of protein-ligand complexes. To deeply understand the specific recognition mechanism between AqAga and the agar tetramer, we systematically analyzed potential hydrogen bonds.

   Hydrogen bond analysis was performed using ChimeraX's built-in “H-Bonds” tool. In the menu, navigate to “Tools” → “Structure Analysis” → “H-Bonds” to open the panel. We set geometric criteria: distance tolerance = 0.4 Å, angle tolerance = 20°.

   

   Figure 13 Hydrogen bonding between AqAga and tetramer.

   Results show that the ligand AqAga forms 9 hydrogen bonds with key residues (TRP 133, ASP 137, THR 150, TRP 155, LYS 177, CYS 307, THR 622, VAL 631, TYR 679), forming a dense interaction network. Typical hydrogen bond distances range from 2.5-3.5 Å. These bonds firmly "anchor" the ligand within the protein's active site.

   Among them, 8 bonds are shorter than 3.5 Å, indicating strong interactions, further confirming this region as a likely key active pocket. These short hydrogen bonds form a coordinated and stable binding interface, significantly enhancing ligand-receptor affinity.

   Notably, the hydrogen bond involving CYS 307 has a length of 3.728 Å, slightly beyond the typical strong hydrogen bond range, and its spatial position deviates from the core of the active pocket, only binding at the ligand's terminus. Thus, its functional contribution may be limited or even negligible.

   1.2 Interaction Analysis

   Besides hydrogen bonds, other interactions such as van der Waals (vdW) forces exist. Since most interactions are distance-dependent, we used the “Contacts” tool to analyze all atomic contacts and indirectly infer possible interactions. Access via “Tools → Structure Analysis → Contacts”.

   

   Figure 14 Contact Analysis bonding between AqAga and tetramer.

   AqAga and the tetramer form 129 contacts, mostly concentrated in the active pocket region identified in hydrogen bond analysis. Notably, the ligand's carbon atoms severely overlap with the two terminal methyl groups (CG2) of VAL 631's side chain (1.026 Å and 0.522 Å overlaps), and closely contact the central atom (CZ) of PHE663's benzene ring (0.477 Å). These strong hydrophobic clashes indicate excessive steric fit, causing significant van der Waals repulsion.

   Meanwhile, complex interactions exist with GLU 677: its OE1 oxygen forms a strong hydrogen bond (2.671 Å) with the ligand, while OE2 oxygen has a 0.574 Å conflict. This "attract-and-repel" situation strongly suggests that GLU 677's side-chain conformation is critical for binding—its precise orientation determines whether a favorable hydrogen bond or unfavorable collision occurs. Additionally, a 0.402 Å conflict between a ligand hydrogen and GLU 468's carbon (CD) represents an unusual and unfavorable interaction, possibly a key site reducing ligand activity.

   In summary, this binding site is a complex environment of hydrophobic and polar interactions. The ligand is anchored by strong hydrophobic effects in the VAL 631/PHE 663 pocket but pays a price of local steric clash. Its ability to form stable hydrogen bonds with GLU 677 depends heavily on the side-chain conformation. Future optimization should focus on moderately reducing the ligand's hydrophobic core volume to relieve clashes and fine-tuning polar interactions with GLU 677 to enhance binding affinity and selectivity.

   1.3 Pocket Analysis

   Substrate channel analysis is a computational and structural biology method to identify possible pathways inside or between large biomolecules (e.g., enzymes) that allow small molecules (substrates, products, ions) to travel from the protein surface to internal active sites. In short, it finds "tunnels" or "channels" in proteins.

   We used CAVER Analyst 2.0, a desktop software for calculating, analyzing, visualizing, and comparing tunnels in proteins, offering a user-friendly interface.

   

   Figure 15 (Left Panel) Surface accessibility of AqAga.

   (Right Panel) Caver substrate channel calculation.

   We evaluated AqAga's surface accessibility, showing the agar tetramer binds within an internal cavity. Further channel analysis indicates the tetramer can traverse the entire active channel (red region), suggesting significant length and depth in the binding pocket, despite some prediction deviations.

   

   Figure 16 Hydrophobicity analysis of active pocket.

   Molecular Lipophilicity Potential (MLP) calculated by ChimeraX shows an alternating pattern of hydrophobic (yellow) and hydrophilic (blue) regions at the binding interface, but overall dominated by hydrophilic amino acids. This aligns with typical polysaccharide-degrading enzymes: the catalytic center is hydrophilic (to bind hydrophilic substrates and accommodate acidic residues like GLU and ASP as proton donors for glycosidic bond cleavage), surrounded by key hydrophobic regions that maintain full catalytic function. Based on MLP and spatial features, we speculate the large hydrophilic region on the right may be the substrate channel entrance, highly hydrophilic and solvent-exposed, favoring initial substrate binding and entry.

   2 agaNash-Dimer

   2.1 Hydrogen Bond Analysis

   

   Figure 17 Hydrogen bonding between agaNash and dimer.

   A total of 8 intermolecular hydrogen bonds formed between agaNash and the dimer: 4 conventional and 4 non-canonical (distance-only). Conventional bonds include the ligand's carbonyl oxygen with HIS205 (ND1, D..A 3.170 Å) and ASN 149 (OD1, D..A 2.987 Å), both within typical range. Non-canonical interactions involve ARG 234 (NH2, D..A 2.820 Å) and HIS263 (NE2, D..A 3.028 Å).

   Overall, multiple hydrogen bonds indicate a broad polar interaction network. Bonds with ASN 149 and ARG 234, with shorter distances, likely contribute more to stability. Longer bonds (e.g., LYS 203 at 3.371 Å) suggest weaker contributions. This multi-site hydrogen bonding, combined with other non-bonded interactions, maintains the complex's binding conformation, suggesting moderate-to-high binding affinity and stability.

   2.2 Interaction Analysis

   

   Figure 18 The interaction between agaNash and dimer.

   A total of 82 contact sites formed between agaNash and the dimer, with several significant spatial conflicts. The most severe clashes occur with PHE125 (CE2, overlap = 0.814 Å) and THR 126 (CG2, overlap = 0.720 Å), indicating strong steric repulsion that may destabilize the complex.

   Moderate overlaps also exist with ARG 234 (NH1, 0.551 Å) and TRP 88 (CH2, 0.413 Å), indicating local geometric mismatch. Notably, some hydrogen-bonding sites (e.g., HIS 205 ND1, ASN 149 OD1) show near-zero or slightly positive overlaps, meaning they form polar interactions but are in tight contact or slight conflict.

   Conversely, many contact sites show ideal negative overlaps, indicating effective hydrophobic matching and vdW interactions. For example, TRP 88 (CE2, -0.249 Å), PHE76 (CE2, -0.259 Å), LYS 203 (NZ, -0.359 Å), and other aromatic/hydrophobic residues form a favorable hydrophobic interface, with overlaps mostly between -0.1 Å and -0.4 Å. These favorable interactions may offset local clashes and provide the main driving force for complex stability.

   In summary, while significant spatial conflicts exist, the extensive network of negative-overlap contacts suggests good steric complementarity. Overall stability may depend on structural plasticity in conflict regions and synergistic hydrophobic effects.

   2.3 Pocket Analysis

   

   Figure 19 (Left Panel) Surface accessibility of AqAga.

   (Right Panel) Analysis of AqAga's pocket hydrophilicity and hydrophobicity.

   The hydrophobic core of the binding pocket is formed by aromatic and aliphatic residues such as TRP 88, PHE 76, PHE 125, and HIS 263, creating a continuous hydrophobic binding interface. Strong hydrophobic interactions (e.g., TRP 88 CE2, PHE 76 CE2, LYS 203 NZ) provide the primary driving force for ligand binding, forming a stable energy foundation.

   Polar interactions are mediated by residues like HIS 205, ASN 149, and ARG 234, forming an 8-hydrogen-bond network with the ligand. HIS 205 (ND1) and ASN 149 (OD1) form conventional hydrogen bonds, while ARG 234 and HIS 263 form non-canonical ones. These polar anchors not only directly participate in binding but also spatially interweave with hydrophobic regions, precisely positioning the ligand and enhancing binding specificity.

3.3.3 Result Analysis

Subsequently, we conducted systematic enzymatic hydrolysis experiments to validate the computational model. Enzyme activity assays showed that AqAga exhibited relatively higher activity among the agarases; AqNABHshowed the highest activity among the neoagarobiose hydrolases (see Figure 21 below).

By comparing prediction results with experimental data, we found that both agarase and neoagarobiose hydrolase activities were highest in the Sq-Ag5 strain. This finding not only confirms the accuracy of the computational model but also provides the optimal enzyme combination for efficient agar degradation. Ultimately, we selected this enzyme combination strain for preliminary fermentation experiments.

Figure 21 Enzyme activity assay of a combination of two enzymes.

4.3.1 Notation Conventions

Other symbols will be explained when referenced later.

4.3.2 Assumptions

1. Glucose, the readily available carbon source, is preferentially consumed and approaches depletion by 48 hours;
2. 3,6-Dehydrated-L-galactose does not accumulate intracellularly and does not affect other carbon source metabolism. It serves solely as an "ineffective product" of agar degradation and requires no additional monitoring.

4.3.3 Parameter Preparation

To characterise promoter expression levels, we quantified the relative expression intensity of each promoter by compiling distinct features from P_TEF, P_GPD, and P_TDH promoters used in wet experiments.

First, measurements of the following indicators of promoter expression strength were collected: mRNA [11] and GFP levels [12].The relative intensity ratios of these three were calculated and normalised to serve as weighted features. Using P_TEF as the baseline, the GFP equivalent intensity value for P_GPD was calculated: $$ P_{\text{GPD,GFP,equ}} = \text{TEF1}_{\text{GFP}} \times \frac{P_{\text{GPD,mRNA}}}{\text{TEF1}_{\text{mRNA}}} \tag{2}$$

The total intensity is calculated using the following formula: $$ S_{\text{total}} = S_{P_{\text{TEF}}} + S_{P_{\text{GPD}}} + S_{P_{\text{TDH}}} = P_{\text{TEF,GFP}} + P_{\text{GPD,GFP,equ}} + P_{\text{TDH3,GFP}} \tag{3}$$

The normalised weight is calculated as: $$ W_j = \frac{S_j}{S_{total}}, \,\, j\in \left\{P_{TEF},P_{GPD},P_{TDH}\right\} \tag{4}$$

Substituting specific values yields the final weight distribution as shown in Table 6.

Promoter combinations may be categorised into homopromoter and heteropromoter combinations, for which the synergistic enhancement coefficient $\alpha$ and average coordination coefficient $\beta$ are defined respectively. The relative expression intensity formula for homopromoter combinations is as follows: $$ \lambda = \alpha \times \lambda_i \tag{5}$$

The relative expression intensity formula for heterologous promoter combinations is as follows: $$ \lambda = \beta \times \frac{\lambda_i + \lambda_j}{2},i \ne j \tag{6}$$

4.3.4 Promoter-regulated Flux Model for Squalene Synthesis

Based on bioreactor kinetics, the carbon source metabolic pathway and key steps of squalene synthesis were decomposed and categorised to construct a system of ordinary differential equations. This model utilises a mixed medium of 25 g/L agarose and 5 g/L glucose as the carbon source, with its metabolic network decomposed into three components: carbon source degradation, intermediate product branching, and squalene synthesis.

1. Carbon Source Degradation Kinetic Model

Based on the classical Monod equation for microbial carbon source metabolism [19], the glucose consumption rate is positively correlated with its concentration ($S_1$). Furthermore, when cell density is too low, synthesis may be limited by insufficient total decomposing enzymes despite a high glucose concentration ($S_1$). Therefore, a cell half-saturation constant ($K_x$) is introduced. When cell density far exceeds the saturation threshold ($K_x$), the cell density regulation coefficient approaches 1. Based on the hypothesis of glucose prioritisation for synthesis and demand analysis, the objective is to maximise squalene yield while satisfying growth requirements. Therefore, a priority weighting parameter is introduced to quantify the stronger pull of synthesis on carbon source consumption, as detailed in the following formula. $$ {\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{7}$$ $$ S_1(t) = 0 \,\,(t>48h) \tag{8}$$

During the early phase, agarase degradation proceeds at a slower rate, primarily supplying nutrients for growth and development. In the later phase, when glucose is depleted, it is mainly utilised for squalene synthesis. Similar to the glucose consumption equation, the synthesis priority coefficient is adjusted due to differing allocations, as detailed in the following formula. $$ \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{9}$$

Subsequently, a dynamic carbon source allocation mechanism is established, centred on maximising squalene yield. A baseline reduction parameter ensures adequate and secure nutrients for early-stage cell growth and development. A cell density coefficient is introduced to provide real-time feedback, preventing excessive nutrient supply when cells are already sufficient. The specific carbon source allocation ratio equation is as follows. $$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{10} $$

2. Cell Growth Kinetic Equation

Observation of experimental data indicates rapid squalene yield growth around 72 hours. Therefore, a growth inhibition coefficient is introduced post-72 hours, increasing over time to ensure carbon sources are prioritised for squalene synthesis pathways once growth requirements are met. The specific formula is as follows. $$ {\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{11}$$

3. Single-Objective Model for Squalene Flavour Synthesis

1. Establishment of the Synthesis Equation

   Under promoter regulation, the squalene synthesis equation is constructed by integrating carbon source allocation and growth development. From an enzymatic and promoter regulatory perspective, higher relative promoter expression strength correlates positively with effective regulatory capacity. This, combined with cell density X(t), provides the foundation for the synthetic pathway by correlating positively with total enzyme expression. Secondly, integrating differential equations from Sections 3.1 and 3.2, the competitive equilibrium between growth and synthesis during the pathway is considered. A higher proportion of synthetic carbon sources enhances conversion efficiency towards squalene, concentrates metabolic flux, and reduces side reactions. The specific equation is as follows. $$\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{12}$$

   Where $k_q$ denotes the catalytic efficiency of enzymes in the MVA pathway, $\lambda_i$ represents the relative expression strength of the promoter, $eta_{q,i}(t)$ indicates the carbon source conversion efficiency.

2. Objective Function and Constraints

   Based on experimental data, yield growth becomes extremely sluggish after 96 hours. Therefore, 96 hours can be set as the final time point for synthetic yield maximisation, denoted as $t_{end}$。

   To prevent premature depletion of agarose during the late fermentation stage, which would halt squalene synthesis due to lack of substrate and hinder sustained yield increase beyond 96 hours, setting constraints on residual carbon sources ensures stable carbon supply throughout the fermentation cycle. Secondly, establishing a requirement that the squalene synthesis rate must not fall below 0.001 mmol/(L·h) after 36 hours prevents prolonged interruptions in synthesis, thereby safeguarding the continuity of squalene production. The single-objective constrained model is formulated as follows. $$ \max_{r(t)} Q_i(t_{\text{end}}) \tag{13}$$ $$ \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{14}$$

4. Solution via Runge-Kutta Algorithm + Constrained Nonlinear Optimisation

Based on the parameter analysis and initial conditions from Section 3.3, the following system of ordinary differential equations is established: $$\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{15}$$

Where, $t \in \left[t_{start},t_{end}\right]=[0.96]$, and the initial conditions are $y_0=y(0)=\left[10.25,0.1,0 \right]^T$

The solutions were obtained using MATLAB's ode45 function, employing the fourth-order Runge-Kutta method[16] to provide candidate solutions, with a fifth-order variable-step numerical solution method controlling the error. $$\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{16}$$

Secondly, parameter boundary constraints are established to formulate the constrained nonlinear error optimisation model as follows: $$\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{17}$$ $$\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{18}$$

Where $p=[p_1,p_2,...,p_{\varGamma}]^T$ denotes the parameter vector to be optimised, and $w_j$represents the time point weights corresponding to 24h, 48h, 72h, and 96h. $lb_k,ub_k$specify the lower and upper bounds for parameters respectively, $g_j(p)$ constitutes the inequality constraint function, and $h_k(p)$ defines the equality constraint function.

Based on experimental data from the six groups—"P_TDH+P_TDH", "P_GPD+P_TDH", "P_TEF+P_GPD", "P_TEF+P_TEF", "P_GPD+P_GPD", and "P_TDH+P_GPD"—as shown in the table, training is conducted to optimise the model's initial parameter settings.

Given that the aforementioned model constitutes a medium-scale constrained nonlinear optimisation problem, the Sequential Quadratic Programming (SQP) algorithm is employed, with analyses conducted on its convergence and stability. Subsequent optimisation and refinement of the algorithm may be undertaken as the problem scale increases.

First, the Lagrange function is defined as follows: $$ \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{19}$$

The KKT conditions must be satisfied, encompassing the stationarity condition, primal feasibility, and dual feasibility conditions, specifically:

Where $p^*$ denotes the optimal solution, $\lambda ^*=[\lambda_{1}^{*},\lambda_{2}^{*},...,\lambda_{m}^{*}]^T$ represents the Lagrange multipliers for inequality constraints, and $\lambda ^*=[\mu_{1}^{*},\mu_{2}^{*},...,\mu_{n}^{*}]^T$ denotes the Lagrange multipliers for equality constraints.

5. Promoter Prediction Based on ODE and Error Optimisation

Predictions were made for the remaining three possible grouping scenarios of the P_TEF, P_GPD, and P_TDH promoters, with experimental testing conducted for each scenario. Given our research objective--to identify optimal promoter combinations for wet experiments while reducing blind exhaustive testing rather than entirely replacing wet experiments. Therefore, focusing solely on the relative ranking of final squalene yields suffices to guide experimental planning. Comparing predicted and experimental squalene yield rankings reveals the following descending order: "P_TEF+P_TDH", "P_GPD+P_TEF", "P_TDH+P_TEF"--consistent with experimental results. The comparison between predicted and experimental outcomes is illustrated below.

Figure 23 Comparison of predicted and experimental results.

6. Model Evaluation

The coefficient of determination R² was calculated, alongside a Kolmogorov-Smirnov (K-S) test, to compare discrepancies between model outputs and experimental data, thereby assessing model accuracy.

The overall coefficient of determination R² for the three sets of tested promoter combinations was 0.7987. Notably, the "P_TDH+P_TEF" combination achieved an R² value as high as 0.895, indicating the model possesses strong interpretability and can effectively characterise the regulatory patterns of promoter combinations on carbon source-dependent olefin production. Further parameter optimisation can be pursued by incorporating additional test data, as detailed in the table below.

The correlation scatter plot is depicted below. It can be observed that the data points cluster near the 1:1 line and fall within the 95% confidence ellipse. This ellipse represents the region with a 95% probability of encompassing all "experimental value - predicted value" data points. The overall centre of the ellipse is slightly lower. Considering the comprehensive nature of the model, this systematic underestimation may stem from numerous constraints and falls within the normal range.

Figure 24 Correlation scatter plot.

The K-S test is employed to verify whether the model-predicted data and experimental data originate from the same distribution, serving as a key indicator for assessing the statistical consistency of the model. We further validated the distributional consistency between predicted and experimental data via the K-S test. Results yielded a significance level p-value of 0.1862 (p > 0.05), accepting the null hypothesis. This indicates that at the 5% significance level, the model predictions and experimental data originate from the same distribution. This avoids scenarios where the model fits well only at isolated points while deviating overall, confirming the model possesses sound statistical consistency and reliability.

4.3.5 Inducible Promoter Optimisation Model

Building upon the aforementioned combinatorial promoter-based model, we further incorporated inducible promoters to more realistically simulate actual biological systems, thereby expanding the model's applicability. Unlike combinatorial promoters, inducible promoters (e.g., P_GAL1, P_GAL7, P_GAL10) exhibit expression levels regulated by external inducer concentrations. This characteristic may influence the balance between cellular growth and product synthesis. The following three promoters were selected for analysis: P_GAL1 [13], P_GAL7 [14], P_GAL10 [15].

1. AHP-Based Evaluation Model for Inducible Promoters

First, we employed the Analytic Hierarchy Process (AHP) to construct an evaluation model for inducible promoters. Key performance parameters were obtained through literature review and experimental analysis, establishing the following target, criterion, and scheme layers.

Using pairwise comparisons via a 1-9 rating scale, construct criterion layer judgment matrix A and calculate weights to determine each criterion's relative importance weighting for the target, along with each scheme's relative merit score under each criterion. Compute the eigenvector of this judgment matrix to obtain criterion weights . Calculate consistency ratio CR < 0.1; validation confirms the judgment matrix's validity.

Next, for each criterion $C_j$, we conducted pairwise comparisons of the promoters, constructed a judgement matrix, and calculated relative scores $s_{ij}$. For each promoter i, its total score $S_i$ was calculated as follows, with results presented in Table 12. The calculated scores for each criterion and the comprehensive score provide preparation for optimising the squalene pathway synthesis model. $$ S_i = \sum_{j=1}^{4} (w_j \times s_{ij}) \tag{21}$$

2. Modification of the Squalene Synthesis Flux Model

Agarose polysaccharides in red algae are broken down into neoagarose disaccharides by agarase, which are then hydrolysed into D-galactose and 3,6-anhydro-L-galactose by neoagarose hydrolase. Inducible promoters differ from combinatorial promoters in that their expression levels are influenced by inducers. Building upon the original model, the previously fixed relative expression intensity has been updated to a dynamic expression intensity. An additional induction efficiency factor has been introduced. Initial values are based on Section 1, with subsequent optimisation.

Secondly, a new galactose consumption kinetic equation was introduced to describe its induction and carbon source functions. $$ \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{22}$$

Here, $k_3$ denotes the galactose consumption rate constant, with an initial value of 0.08 $h^{-1[3]}$。

As part of the galactose is utilised for induction rather than direct metabolism, a galactose carbon source contribution coefficient $c$ is introduced. Due to transcription and translation delays, 3h[13] is required after galactose addition before an increase in promoter activity becomes observable; a delay function $\tau(t)$ is introduced to describe this phenomenon. The revised synthetic equation is as follows. $$ \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{23}$$ $$ \tau(t) = \begin{cases} 0, & t - t_{\text{inducer}} < 3 \\ 1, & t - t_{\text{inducer}} \geq 3 \end{cases} \tag{24}$$

Here, $t_{inducer}$ represents the time for galactose to enter the reaction

Furthermore, based on Section 4.3.4.3, due to the characteristics of inducible promoters, new constraints were introduced: a residual galactose constraint and a promoter activity constraint. These describe the phenomenon where promoter strength abruptly decreases when becomes excessively low.

3. Results and Analysis

Based on experimental data from the six groups--"P_GAL1+P_GAL1, P_GAL1+P_GAL10, P_GAL7+P_GAL7, P_GAL7+P_GAL1, P_GAL7+P_GAL10, P_GAL10+P_GAL7"--the model was trained to optimise its initial parameter settings.

Predictions were made for three grouping scenarios: "P_GAL1+P_GPD", "P_GAL10+P_GAL1", and "P_GAL10+P_GAL10", followed by experimental testing of these configurations. Comparing the predicted results with the experimental outcomes as shown below, it can be observed that the yield decreases in the following order: "P_GAL1+P_GPD", "P_GAL10+P_GAL10", and "P_GAL10+P_GAL1". This ranking aligns with the experimental findings.

Figure 25 Comparison of predicted and experimental results.

4. Model Evaluation

The coefficient of determination R² was calculated, alongside a Kolmogorov-Smirnov (K-S) test, to compare model outcomes with experimental data and assess model accuracy.

The overall R² for the three sets of test promoter combinations was 0.8384, with "P_GAL10+P_GAL1" reaching 0.9149, indicating the model possesses strong interpretability. Details are presented in the table below.

Further validation of the consistency between the distribution of predicted data and experimental data was conducted via the KS test. The results yielded a significance level p-value of 0.7864 (p > 0.05), accepting the null hypothesis. This indicates that, at the 5% significance level, the model's predicted data and experimental data originate from the same distribution, confirming the model's sound statistical consistency and reliability.

4.3.6 Conclusions

Centring on the key regulatory factor in Rh1 ginsenoside biosynthesis from red algae—the promoter combination—we constructed and refined an optimal promoter prediction model based on ODE kinetics. Through model establishment, solution, prediction, and evaluation, alongside the extension of an inducible promoter model, we provide a quantitative tool for studying squalene synthesis mechanisms. The model demonstrates reliable predictive capability, with performance further enhanced upon extension. The overall R² value increased to 0.8384, and the K-S test p-value reached 0.7864. Ranking predictions for combinations such as "P_GAL1+P_GPD" align with experimental results. By predicting high-yield combinations, the model guides the prioritisation of wet-lab experiments, significantly reducing blind trial counts and lowering time and reagent costs. Subsequent research may further enhance model accuracy by incorporating additional experimental data on promoter combinations and optimising parameters such as the synergy coefficient.

5.3.1 Cost Accounting Model

1 Traditional Preparation Methods

1.1 Cost Accounting Boundaries and Classification Framework

The cost structure for preparing ginsenoside Rh1 is categorized into three major components: fixed costs (facility and raw material base investments), variable costs (dynamic expenditures during production), and indirect costs (wastewater/waste residue treatment). Since land costs for industrial production are calculated similarly across multiple methods, they are excluded from this calculation for simplified comparison. The accounting boundary covers the entire process from raw material acquisition to product separation, as detailed in Table 16.

1.2 Production Process Flow

The traditional extraction process for ginsenoside Rh1 primarily involves obtaining it from ginseng and related raw materials or using ginseng and other plant materials as a basis. Common methods include ethanol-acid hydrolysis, enzymatic hydrolysis, and ultrasonic-assisted extraction. Below, we focus on analyzing the “two-step method combining enzymatic conversion with metal ion catalysis.” This approach offers superior product efficiency compared to other methods, with all process parameters supported by clear experimental data, providing scientific backing for cost quantification.

Fresh ginseng procured from origin sites undergoes preliminary screening and drying. After processing, the stems and leaves of fresh ginseng account for 40%-48% of the total ginseng weight [20], with subsequent calculations using an average value of 44%. Through extraction, purification, and concentration/drying processes applied to ginseng stems and leaves, ginsenoside Rh1 can be obtained with a yield of 35.15% [22]. Using Re as the substrate, Rg1 is prepared via a two-step enzymatic conversion process with a yield of 70.5% [23]. Subsequently, Fe³⁺-catalyzed reaction yields the Rh1 isomer with an 80.2% yield [23]. The flowchart is shown in Figure 26.

Figure 26 Flowchart of traditional ginseng extraction methods.

1.3 Itemized Cost Calculation

Assuming the final goal is to obtain 1g of Rh1, first calculate the raw material requirements. Through derivation, the formula for the required fresh ginseng is as follows: $$ n_0 = \frac{1g}{56.64%} \div 35.15\% \div 44\% \approx 11.43g \tag{25}$$

Based on the reaction system and product yield [23], the required amounts of catalyst and reagents for preparing 1g of Rh1 can be calculated. The cost breakdown for the traditional ginseng extraction method is shown in Table 17.

Based on the reaction system and product yield [20], the variable cost formula for raw materials, catalysts, and reagents required to prepare 1 g of Rh1 is as follows: $$ 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 \,\, CNY/g \tag{26}$$

Indirect costs are: $IC_{tradition} = c_{waste} \times n_6 = 0.006$ CNY/g

Assuming a single production run yields 100g of Rh1, with an annual operating period of 250 days, the formula for allocating fixed costs per unit of ginsenoside Rh1 is as follows: $$FC_{tradition} = (10\times S_1 + S_2 \div (5.8 \times 250)+ S_3 \div (10 \times 250))\div 100g = 18.930 \,\, CNY/g \tag{27}$$

Combining Equations (25) to (27), the calculated total cost per unit of Rh1 is as follows: $$C_{tradition} = VC_{traditon} + IC_{traditon} +FC_{tradition} = 21.517+0.006+18.930 = 40.453 \,\, CNY/g \tag{28}$$

2 Cost Model for Red Algae Biosynthesis

In the red algae biosynthesis process we propose, the red algae powder feedstock used is rich in agar components, with an agar content ranging from 61% to 67.3% [21]. For the following calculations, we adopt a value of 64.15%. During the fermentation stage, agarose serves as the key substrate for medium preparation, with a concentration of 25 g/L. Testing revealed a ginsenoside Rh1 yield of 141.78 mg/L in the fermentation broth.

The yield calculation formula for ginsenoside Rh1 from red algae powder is as follows: $$ Y_{red-alga} = (141.78 \div \frac{Agar Quality in Culture Medium}{Agar Content in Red Algae Powder} \times 100) \times 100\% \approx 0.364\% \tag{29}$$

2.1 Laboratory Cost Calculation

First, the cost of preparing 1g of ginsenoside Rh1 at laboratory scale was calculated. The detailed cost breakdown is shown in Table 18.

Based on experimental data, a single batch yields 141.78 mg/L of ginsenoside Rh1. The variable cost calculation formula is as follows: $$ VC_{lab-alga} = (P_{red-alga} + P_{C_{6}H_{12}O_{6}} \times 10 + P_{media} +C_{sep}) \div 0.14178 = 89.724 \,\, CNY/g \tag{30}$$

The indirect cost calculation formula is as follows: $$IC_{lab-alga} = c_{waste} \times 2 \div 0.14178 = 16.720 \,\, CNY/g \tag{31}$$

Based on experimental data, the average production time per batch is 144 hours. Standard equipment operating hours are calculated at 800 hours per year [25]. The formula for allocating fixed costs per unit of ginsenoside Rh1 is as follows: $$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\,\, CNY/g \tag{32}$$

Combining equations (29) to (31), the calculated cost per unit of Rh1 at laboratory scale is: $$C_{lab-alga} = VC_{lab-alga} + IC_{lab-alga} + FC_{lab-alga} = 104.785 \,\, CNY/g \tag{33}$$

2.2 Industrial Cost Calculation

Assuming industrial production using a 100m³ fermentation tank with a single fermentation cycle of 6 days, experimental data indicates a yield of 141.78 mg/L Rh1 per cycle. Rh1 production = 100m³ × 141.78 mg/L = 14.178 kg. The cost breakdown at industrial scale is shown in Table 19:

Calculate variable costs as follows: $$ 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 \,\, CNY/g \tag{34} $$

The formula for indirect costs is as follows: $$IC_{ind-alga} = {C_{waste}}^{\prime} \times 2 \times 100 = 1.673 \,\,CNY/g \tag{35} $$

Assuming an annual working time of 250 days, the fixed cost per unit of Rh1 is calculated as follows: $$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 \,\, CNY/g \tag{36}$$

By combining Equations (29), (34) to (36), the cost per unit of Rh1 in the laboratory is calculated as: $$C_{lab-alga} = VC_{ind-alga} + IC_{ind-alga} +FC_{ind-alga} = 18.656 \,\,CNY/g \tag{37}$$

3 Cost Comparison and Sensitivity Analysis

3.1 Cost Comparison

As shown in Figure 27, when applied industrially, the red algae biosynthesis method demonstrates significantly lower costs compared to traditional ginseng extraction methods, exhibiting strong cost advantages and industrialization potential.

Figure 27 Cost comparison chart.

3.2 Sensitivity Analysis

The following sensitivity analysis examines raw material price fluctuations:

1. If red algae prices increase by 10%, the total cost of laboratory-scale production using the red algae method rises by 2.894% (to 107.818 yuan/g), while the total cost of industrial-scale production using the same method increases by only 0.053% (to 18.666 yuan/g).
2. A 10% increase in fresh ginseng price raises raw material costs by CNY 0.2156/g. The total cost for the traditional method increases by 0.553% (to CNY 40.678/g), which is 10.434 times that of the industrial production using the red algae method.

This demonstrates that the red algae biosynthesis method is more suitable for industrial production compared to traditional fresh ginseng extraction of Rh1. It exhibits lower sensitivity to raw material price fluctuations, stronger cost stability, and greater risk resilience, providing sustainable support for the large-scale mass production of functional cosmetics containing ginsenosides.

Therefore, the red algae biosynthesis method effectively reduces the cost of ginsenoside Rh1 extraction while minimizing production risks associated with raw material price fluctuations. It demonstrates superior cost control capabilities and market competitiveness in industrial production. Within the supply chain, it ensures long-term production cost stability, prevents supply disruptions, and guarantees consistent efficacy concentration of ginsenoside Rh1. This frees up resources for formulation optimization and R&D, significantly advancing the widespread adoption of ginsenosides in cosmetic ingredients.

5.3.2 Carbon Footprint Quantification and Carbon Utilization Rate Calculation Model

1 Methodology and System Boundary

Using the Life Cycle Assessment (LCA) methodology, we quantitatively modeled the carbon footprint and carbon utilization rate for industrial production of ginsenoside Rh1 via traditional ginseng extraction and red algae biosynthesis. By identifying key carbon emission sources and constructing mathematical calculation models, the potential impact of the product on global warming and the composition of its impact across different stages, processes, and spatial locations (expressed in carbon dioxide equivalents) were calculated. Comparing the carbon emission intensity and carbon utilization efficiency of the two pathways provides scientific basis for low-carbon optimization of the Rh1 preparation process and carbon neutrality goals. Concurrently, tools such as convective dispersion equations were employed to enhance the accuracy of accounting for indirect emissions like wastewater.

2 Modeling Assumptions and Principles

1. Transportation Load Assumption: Based on the “Regulations on Road Goods Transportation and Terminal Management” and common industry load standards, the rated single-trip transport capacity of light-duty trucks is assumed to be 6 tons.
2. Transportation Timing Assumption: Referencing the “Technical Conditions for the Safe Operation of Motor Vehicles” (GB 7258) and highway transportation industry operational norms, the average travel speed of diesel light trucks on trunk highways is assumed to be 65 km/h. Considering driver workload and safety management requirements, a rest period is scheduled after every 6-8 hours of continuous driving to account for evaporative emissions during parking.
3. Exclusion Criteria: During data collection, if individual material or energy flows are found to make no substantial contribution to the carbon footprint of a specific unit process, they may be excluded. The cumulative total of excluded flows shall not exceed 5% of the total greenhouse gas emissions for the calculated unit process [26].
4. Data Quality Requirements: Emission factors shall be selected according to the following priority order:
1. Emission factors obtained from actual measurements or mass balances
2. Emission factors provided by suppliers
3. Regional emission factors
4. National emission factors
5. International emission factors

3 LCA-Based Carbon Footprint Quantification Model

To eliminate interference from product yield variations across different process routes in carbon footprint comparisons, the functional unit is uniformly defined as “production of 1g ginsenoside Rh1.”

The life cycle encompasses the continuous and interconnected stages associated with a product, from raw material acquisition or generation from natural resources to end-of-life treatment. Since both traditional preparation and red algae synthesis of Rh1 share identical subsequent processes—product sales, product use, and waste disposal—this model focuses on comparing the environmental benefits of different Rh1 production methods. The system boundary is defined as “cradle-to-gate,” excluding these post-production stages from quantification.

3.1 Carbon Footprint of Traditional Ginseng Extraction

The traditional ginseng-derived Rh1 lifecycle chain comprises cultivation, transportation, and extraction/purification stages. Cultivation is strictly constrained by climatic and soil conditions, concentrating production in Jilin, Liaoning, and Heilongjiang provinces [27], relying on forest soils and artificial fertilization.

Regarding the regional distribution of enterprises along the industrial chain, upstream cultivation and primary processing companies are predominantly clustered in Jilin Province, forming a distinct industrial cluster effect. In contrast, representative listed enterprises in the ginseng industry exhibit relatively dispersed regional distribution, spanning Northeast, Central, and Southern China. Significant differences exist in the positioning of enterprises across these regions: companies in Northeast China primarily focus on ginseng cultivation, harvesting, and primary sales, serving as the supply source at the beginning of the industrial chain; Enterprises in central and southern regions primarily engage in deep processing and R&D, requiring raw material procurement from northeastern production areas. This creates a lengthy transportation chain between cultivation bases and central/southern processing/R&D facilities. Based on enterprise distribution, our model calculates transportation routes using the example of Tonghua, Jilin (northeastern main production area) to Hangzhou, Zhejiang (extraction processing).

According to Equation(2), producing 1g of Rh1 requires 11.43g of fresh ginseng, meaning a single transport corresponds to approximately 524.934kg of Rh1 output. The pollutant emission factors and correction factors for a diesel light-duty truck (China 4 standard) under baseline conditions (temperature 15°C, humidity 50%, speed 30 km/h, load 50%, diesel sulfur content 50 ppm) are shown in Table 20 [28]. The evaporative emission factor during operation is 11.6 g/h, and during parking is 6.5 g/day [28].

$$Exhaust emissions = EF \times \varphi_1 \times \gamma_1 \times \theta_1 \times \lambda_1 \times Transportation distance \times 10^{-6} \tag{38}$$

According to formula (38), the sum of all tailpipe emissions is calculated to be 0.0046 t. Converted to total carbon emissions using GWP: $$\mathrm { 0.0046t \times }\mathrm {(} \frac {\mathrm{44}} {\mathrm{28}} \mathrm{) + 0.00035t \times 25 + (11.6g/h \times 27h + 6.5g/} \mathrm{day}\mathrm{ \times 0.25 }\mathrm{day}\mathrm{)}\mathrm{\times }\mathrm {10^{-6}} \mathrm{\approx 0.01632 t} \tag{39} $$

Transportation Carbon Emissions per Unit Rh1: $$ E_{transport} = 0.01632t \,\,CO_{2}e \div 524.934 \,\,kg Rh1 \approx 0.0312 \,\, kg \,\, CO_{2}e/kg \,\,Rh1 \tag{40} $$

Carbon emissions during the extraction phase primarily stem from the electrical energy consumed by extraction equipment operations. The calculation formula is as follows: $$ E_{extract} = Total power consumption \times Regional Electricity Carbon Emission Factor \tag{41}$$

Referencing chemical instrumentation websites and industry engineering standards, calculate electricity consumption based on equipment power and operating duration, then substitute into the above formula. Carbon emissions calculations for different stages are shown in Table 21.

According to Table 21, the total carbon emissions are calculated as: $$CF_{extract}=E_{plant} +E_{tranport} + E_{extract} = 2.159 \,\,kg \,\,CO_{2}e/g \,\,Rh1 \tag{41}$$

3.2 Carbon Footprint of Rh1 Biosynthesis

The life cycle chain of red algae biosynthetic Rh1 can be divided into cultivation and harvesting, transportation, and fermentation stages. Based on the Rh1 yield, it is calculated that producing 1g of Rh1 requires 274.725g of red algae raw material.

According to data from the China Fishery Statistical Yearbook, China's primary seaweed production regions are Fujian, Shandong, and Liaoning. Among these, Fujian Province is the nation's largest seaweed production area. From a spatial distribution perspective, China's synthetic biology industry chain exhibits relatively comprehensive layouts in provinces and municipalities such as Zhejiang, Jiangsu, Guangdong, Shanghai, and Beijing. This enables a controllable transportation chain from cultivation bases to processing and R&D facilities, forming a “centralized transportation - centralized processing” model that reduces intermediate transfer points. Carbon emissions during the cultivation and harvesting phase primarily stem from the production of cultivation facilities (e.g., polyethylene aquaculture cages) and electricity consumption for harvesting equipment operation. These emissions are calculated using the following formula, incorporating the polyethylene usage in cultivation facilities, the polyethylene emission factor, and the electricity consumption of harvesting equipment. $${E_{F, H}}^{\prime} = Carbon emissions of aquaculture facilities + Carbon emissions of harvesting \tag{19}$$

Similar to the calculation logic in Section 2.3.1, based on Formula (33), the sum of all tail gas emissions is calculated to be 0.0011 t, corresponding to an Rh1 production of approximately 21.840 kg per transport. Converted to total carbon emissions using GWP: $$0.0011t \mathrm { \times (} \frac{44}{28}) + 0.00035t \mathrm { \times} 25 + (11.6g/h \mathrm {\times} 27h + 6.5g / \mathrm{day \times}0.25 \mathrm {) \times } 10^{-6} \approx 0.3161 kg \tag{42}$$

Carbon emissions from transport per unit of Rh1: $${E_{transport}}^{\prime} = 0.3161 + 21.840 = 0.014 \,\, kg\,\,CO_{2}e/g \,\,Rh1 \tag{43}$$

Carbon emissions during the fermentation stage primarily stem from electricity consumption during fermenter operation, calculated using the following formula: $${E_{ferment}}^{\prime} =Total consumption \times Regional electricity carbon emission factor \tag{44}$$

Similar to the calculation logic in Section 2.3.1, the electricity consumption is calculated based on the device power and operating duration and substituted into the above equation. The carbon emissions calculation results for different stages are shown in Table 22.

The total carbon emissions calculated are: $$CF_{red algae} = E_{plant} + E_{transport} + E_{extract} = 0.023 kg\,\,CO_{2}e/g \,\, Rh1 \tag{45}$$

3.3 Comparative Analysis of Method Carbon Footprints

Based on the above comparative calculations, the carbon emissions across each stage of the system boundary are illustrated in Figure 28. Traditional ginseng extraction requires substantial land and resource inputs during cultivation. Furthermore, the dispersion of origin sites and subsequent processing plants, coupled with long transportation distances, results in persistently high carbon emissions. In contrast, red algae sourced from the ocean exhibit lower carbon emissions during harvesting and other stages. Subsequent processing activities are geographically concentrated, leading to an overall carbon footprint significantly smaller than that of ginseng extraction. This disparity points the way for the industry's low-carbon transition. The red algae production model holds greater potential for carbon reduction, helping the sector decrease reliance on high-carbon production methods while advancing green supply chain development and low-carbon production practices.

Figure 28 Carbon emissions of system boundary components under different pathways.

4 Carbon Utilization Rate Calculation Model

Carbon Utilization (CU) serves as a key metric for evaluating the carbon resource efficiency of a preparation process. It calculates the percentage of carbon mass fixed in the target product ginsenoside Rh1 relative to the total carbon emissions (expressed as CO₂ equivalents) across the entire preparation lifecycle. The formula is as follows: $$CU = \frac{C_{Rh1}}{CF_{total}} \times 100\% \tag{46}$$

The chemical molecular formula of ginsenoside Rh1 is $C_{42}H_{72}O_{14}$,Based on the relative atomic masses of the elements, the mass of carbon contained in 1 g of Rh1 is 0.63 g. Based on the carbon footprint calculations for the two preparation pathways in Section 3.3, $CF_{trandition}=2.159 kg\,\,CO_{2}e/g Rh1, CF_{red algae}=0.023 kg\,\,CO_{2}e/g Rh1,$ Substituting these values into the formula yields: $$CU_{traditon} = \frac{0.63}{2159} \times 100\% \approx 0.0292\% \tag{47}$$ $$CU_{red algae} = \frac{0.63}{23} \times 100\% \approx 2.739\% \tag{48}$$

It can be seen that the carbon utilization rate of the red algae biosynthesis method is approximately 94 times that of the traditional ginseng extraction method, indicating its advantage in carbon resource conversion efficiency. Under the same carbon emissions, the red algae method can fix more carbon into the target product Rh1, reducing carbon resource waste.

The primary reason for the disparity in carbon utilization efficiency between the two methods lies in the significantly lower total carbon emissions of the red algae approach. This stems from characteristics such as low carbon emissions during red algae cultivation, short transportation distances, and reduced fermentation energy consumption. The higher carbon utilization efficiency of the red algae method not only aligns with the low-carbon industrialization requirements under the “dual carbon” goals but also reduces potential costs for enterprises under policies like carbon trading and carbon taxes. This further enhances its comprehensive competitiveness for mass production applications in cosmetics.

5.3.3 Sustainability Analysis of Rh1 Biosynthesis from Red Algae

Using “1g ginsenoside Rh1” as the benchmark calculation unit, a “economic cost - environmental benefit” dual-contrast model was constructed. Focusing on three indicators—cost composition, carbon emission intensity, and carbon utilization rate—the model quantitatively compared two preparation pathways, revealing significant differences as shown in Figure 29. The red algae biosynthesis method demonstrates superiority over traditional ginseng extraction in industrial applications across cost, carbon emissions, and carbon utilization rates. This approach reduces procurement costs for cosmetic ingredients, provides flexibility in end-product pricing, and effectively advances the large-scale production of the rare ginsenoside Rh1. Concurrently, it aligns with both the “dual carbon” goals and the industry's demand for green upgrading.

Figure 29 Comparison of costs and carbon utilization rates across different pathway