Overview

Our model serves two main purposes:

Quantitative Description of Project Design: Our dual-microbe system involves complex processes spanning lignin degradation, metabolic conversion, and microbial population dynamics. Through mathematical modeling, we can quantitatively describe these intricate interactions, providing crucial insights into system behavior that would be challenging to obtain through wet lab experiments alone.

Computational Methods for Project Engineering: The models help determine key parameters for system optimization, reduce experimental trial-and-error, connect different wet lab modules, and mathematically encapsulate biological components for predictable engineering.

Our modeling framework consists of three interconnected parts:

1. ODE models - Simulate the kinetic processes of lignin decomposition and metabolic conversion, providing dynamic insights into reaction rates.
2. Flux Balance Analysis (FBA) - Predicts optimal metabolic fluxes and theoretical succinate yields at genome scale.
3. Dual-Species System Dynamics - Explore population interactions and system stability.

These models provide a comprehensive understanding of our project and yield valuable computational results for both fundamental research and industrial application.

Model for Lignin Decomposition and Succinate Synthesis

I. Model Design

We examine the reaction equations and enzyme parameters of the four main parts in the project design: “Lignin decomposition by laccase”, “Metabolic entry of decomposition products into TCA”, “TCA cycle”, and “Succinate extraction”, and characterize the system’s metabolism through ordinary differential equations. We hope to deepen our understanding of the project system and provide guidance for further improvements in wet lab experiments.

II. Pre-assumptions

• The system is well-mixed in compartments, ignoring diffusion limitations inside the culture medium and cytoplasm, but considering the separation between culture medium and cytoplasm.
• Oxygen concentration is suitable and not a limiting factor.
• pH and temperature remain constant during the reaction.
• The total concentration of various enzymes remains constant.
• Products of non-equilibrium reactions do not inhibit the reactions.
• Complex free radical chain reactions are simplified as depolymerization.
• The baseline production of various substances in the engineered bacteria is stable.
• Parameters for some enzymes that were not queried are given by referring to parameters of similar enzymes.
• The total cellular protein concentration is estimated as 250g/L.

III. Lignin Decomposition by Laccase

1.Model Variables

2.Metabolic Reactions

Reaction 1: Direct oxidation of phenolic lignin (L_p) by laccase

$$E + L_p \xrightarrow{k_1} E:L_p \xrightarrow{k_{cat,p}} E_{red} + L_p^{\bullet}\quad (phenolic\text{ }lignin\text{ }radical)$$ $$E_{red} + O_2 \xrightarrow{k_{ox}} E + H_2O\quad (enzyme\text{ }regeneration)$$ $$L_p^{\bullet} \xrightarrow{k_{depoly,p}} a\cdot Syr+b\cdot Van+c\cdot pHB\quad(radical-initiated\text{ }depolymerization)$$ Lignin is composed of three types of monomers: S, G, and H. When lignin is completely decomposed by laccase, these three monomers form syringate, vanillic, and p-Hydroxybenzoic, respectively.

Assuming enzyme regeneration and radical depolymerization are relatively fast, the entire process can be viewed as the consumption of L_p. Since L_p and M_red compete for the enzyme E, the consumption rate v₁ of L_p is described using the competitive inhibition Michaelis-Menten equation: $$v_1 = \frac{V_{max,p} [L_p]}{K_{m,p}(1 + \frac{[M_{red}]}{K_{m,m}}) + [L_p]}$$

Reaction 2: Oxidation of mediator (M_red) by laccase

$$E + M_{red} \xrightarrow{k_2} E:M_{red} \xrightarrow{k_{cat,m}} E_{red} + M_{ox}$$ $$E_{red} + O_2 \xrightarrow{k_{ox}} E + H_2O\quad(enzyme\text{ }regeneration)$$

Similarly, considering the competition between L_p and M_red, the generation rate of M_ox (i.e., the consumption rate of M_red) v₂ is described as: $$v_2 = \frac{V_{max,m} [M_{red}]}{K_{m,m}(1 + \frac{[L_p]}{K_{m,p}}) + [M_{red}]}$$

Reaction 3: Oxidation of non-phenolic lignin (L_np) by oxidized mediator (M_ox)

$$M_{ox} + L_{np} \xrightarrow{k_3} M_{red} + L_{np}^{\bullet}\quad(non-phenolic\text{ }lignin\text{ }radical)$$ $$L_{np}^{\bullet} \xrightarrow{k_{depoly,np}} d\cdot Syr+e\cdot Van+f\cdot pHB\quad(radical-initiated\text{ }depolymerization)$$

Simplification: Assuming radical depolymerization is fast, the reaction rate v₃ is mainly determined by the concentrations of M_ox and L_np: $$\begin{align}v_3 = k_3 [M_{ox}] [L_{np}]\end{align}$$

3.Parameter List:

4.Ordinary Differential Equations:

$$\frac{d[L_p]}{dt} = v_{L_{p}}-v_1 $$ $$\frac{d[L_{np}]}{dt} = v_{L_{np}}-v_3 $$ $$\frac{d[M_{ox}]}{dt} = v_2 - v_3 $$ $$ \begin{align} [M_{red}] = [M_{total}] - [M_{ox}] \end{align} $$ $$\frac{d[\text{Syr}]}{dt}=a\cdot v_1+d\cdot v_3$$ $$\frac{d[\text{Van}]}{dt}=b\cdot v_1+e\cdot v_3$$ $$\frac{d[\text{pHB}]}{dt}=c\cdot v_1+f\cdot v_3$$

5.Result Presentation:

IV.Metabolic Entry of Decomposition Products into TCA

1.Model Variables

2.Metabolic Reactions

Reaction 4: Entry of Syringate into the TCA cycle

Syringate will be gradually metabolized into pyruvate and oxaloacetate via the gallate metabolic pathway in Pseudomonas putida. $$\begin{align}Syr \longrightarrow Pyr + OAA\end{align}$$ To maintain model interpretability, we equate the metabolic rate of Syringate inside Pseudomonas putida to the rate catalyzed by a certain step enzyme, thus using the Michaelis-Menten equation for description: $$v_4=\frac{V_{max}^{Syr}\cdot [Syr]}{K_m^{Syr}+[Syr]}$$ In the fermentation system, since the laccase decomposition process occurs in the total liquid, subsequent metabolic steps occur in the bacterial cells. However, the bacterial volume only accounts for a small part of the total liquid volume, so the consumption rate of Syringate is corrected by multiplying by a coefficient rt. But the generation rates of pyruvate and oxaloacetate are not multiplied by the coefficient rt, as they describe concentrations inside the engineered bacteria, not in the fermentation system.

Reaction 5: Entry of Vanillic into the TCA cycle

Vanillic will be gradually metabolized into succinate and acetyl-CoA via the Protocatechuate metabolic pathway in Pseudomonas putida. $$Van \xrightarrow[+CoA]{} Succ + Ac_{CoA}$$ To maintain model interpretability, we use the simplification method of Reaction 4 for description, and similarly multiply the consumption rate of vanillic by the coefficient rt, but not the generation rates of succinate and acetyl-CoA: $$v_5=\frac{V_{max}^{Van}\cdot [Van]}{K_m^{Van}+[Van]}$$

Reaction 6: Entry of p-Hydroxybenzoic into the TCA cycle

p-Hydroxybenzoic will also be gradually metabolized into succinate and acetyl-CoA via the Protocatechuate metabolic pathway in Pseudomonas putida. $$pHB \xrightarrow[+CoA]{} Succ + Ac_{CoA}$$ To maintain model interpretability, we use the simplification method of Reaction 4 for description, and similarly multiply the consumption rate of p-Hydroxybenzoic by the coefficient rt, but not the generation rates of succinate and acetyl-CoA: $$v_6=\frac{V_{max}^{pHB}\cdot [pHB]}{K_m^{pHB}+[pHB]}$$

3.Parameter List:

4.Ordinary Differential Equations

$$\frac{d[\text{Syr}]}{dt}=a\cdot v_1+d\cdot v_3-v_4\cdot rt$$ $$\frac{d[\text{Van}]}{dt}=b\cdot v_1+e\cdot v_3-v_5\cdot rt$$ $$\frac{d[\text{pHB}]}{dt}=c\cdot v_1+f\cdot v_3-v_6\cdot rt$$

$$\frac{d[\text{OAA}]}{dt}=v_{OAA}+v_4$$ $$\frac{d[\text{Pyr}]}{dt}=v_{Pyr}+v_4$$ $$\frac{d[\text{Ac}_\text{{CoA}}]}{dt}=v_{Ac_{CoA}}+v_5+v_6$$ $$\frac{d[\text{Succ}]}{dt}=v_{Succ}+v_5+v_6$$

5.Result Presentation:

V.TCA Cycle:

1.Model Variables

2.Metabolic Reactions

Reaction 7: Conversion of Pyruvate to Acetyl-CoA

$$Pyr \xrightarrow[+CoA-2.5ATP]{PDHc} Ac_{CoA}$$ The reaction mechanism is a ping-pong mechanism, so its reaction rate is described using the corresponding Michaelis-Menten equation variant: $$v_7=\frac{V_{max}^{PDHc}\cdot [Pyr]\cdot [CoA]}{K_m^{PDHc,Pyr}\cdot [CoA]+K_m^{PDHc,CoA}\cdot [Pyr]+[Pyr]\cdot [CoA]}$$

Reaction 8: TCA cycle (part1)

$$Ac_{CoA} + OAA \xrightarrow[-CoA]{CS} Cit \xleftrightarrow{ACO} Icit \xrightarrow[-2.5ATP]{IDH} \alpha KG \xrightarrow[+CoA-2.5ATP]{\alpha KGDHC} Succinyl_{CoA} \xrightarrow[-CoA-ATP]{SCS} Succ $$ For most enzyme-catalyzed reactions in the TCA cycle, we use the Michaelis-Menten equation to describe their reaction rates. Specifically, for the step $Ac_{CoA} + OAA \xrightarrow[-CoA]{CS} Cit$, the reaction mechanism is an ordered sequential mechanism; for the step $\alpha KG \xrightarrow[+CoA-2.5ATP]{\alpha KGDHC} Succinyl_{CoA}$, the reaction mechanism is a ping-pong mechanism; their reaction rates are described using the corresponding Michaelis-Menten equation variants. For the step $Cit \xleftrightarrow{ACO} Icit$, the forward and reverse reaction rates are both described using the standard Michaelis-Menten equation. $$v_{8.1}= \frac{V_{max}^{CS}\cdot [OAA]\cdot [Ac_{CoA}]}{K_m^{CS,Ac_{CoA}}\cdot [OAA] + K_m^{CS,OAA}\cdot [Ac_{CoA}] + [OAA]\cdot [Ac_{CoA}] + K_s^{CS,OAA}\cdot K_m^{CS,Ac_{CoA}}}$$ $$v_{8.2+}=\frac{V_{max}^{+ACO}\cdot [Cit]}{K_m^{+ACO}+[Cit]}$$ $$v_{8.2-}=\frac{V_{max}^{-ACO}\cdot [Icit]}{K_m^{-ACO}+[Icit]}$$ $$v_{8.3}=\frac{V_{max}^{IDH}\cdot [Icit]}{K_m^{IDH}+[Icit]}$$ $$v_{8.4}=\frac{V_{max}^{\alpha KGDHC}\cdot [\alpha KG]\cdot [CoA]}{K_m^{\alpha KGDHC,CoA}\cdot [\alpha KG] + K_m^{\alpha KGDHC,\alpha KG}\cdot [CoA] + [\alpha KG]\cdot [CoA]}$$ $$v_{8.5}=\frac{V_{max}^{SCS}\cdot [Succinyl_{CoA}]}{K_m^{SCS}+[Succinyl_{CoA}]}$$

Reaction 9: TCA cycle (part2)

$$OAA \xleftrightarrow[+2.5ATP]{MDH} Mal \xleftrightarrow{FH} Fum \xrightarrow[+1.5ATP]{FRD} Succ$$ For the step $Fum \xrightarrow[+1.5ATP]{FRD} Succ$, the Michaelis-Menten equation is used to describe its reaction rate. For the steps $OAA \xleftrightarrow[+2.5ATP]{MDH} Mal$ and $Mal \xleftrightarrow{FH} Fum$, the forward and reverse reaction rates are both described using the standard Michaelis-Menten equation. $$v_{9.1+}=\frac{V_{max}^{+MDH}\cdot [OAA]}{K_m^{+MDH}+[OAA]}$$ $$v_{9.1-}=\frac{V_{max}^{-MDH}\cdot [Mal]}{K_m^{-MDH}+[Mal]}$$ $$v_{9.2+}=\frac{V_{max}^{+FH}\cdot [Mal]}{K_m^{+FH}+[Mal]}$$ $$v_{9.2-}=\frac{V_{max}^{-FH}\cdot [Fum]}{K_m^{-FH}+[Fum]}$$ $$v_{9.3}=\frac{V_{max}^{FRD}\cdot [Fum]}{K_m^{FRD}+[Fum]}$$

Reaction 10: Regeneration of OAA

$$Pyr\xrightarrow[+ATP]{PC} OAA$$ Due to the break in the TCA cycle, OAA cannot be regenerated through the TCA cycle, so the pyruvate carboxylation catalyzed by PC becomes the main generation pathway for OAA. The standard Michaelis-Menten equation is used to describe its reaction rate. $$v_{10}=\frac{V_{max}^{PC}\cdot [Pyr]}{K_m^{PC}+[Pyr]}$$

3.Parameter List:

4.Ordinary Differential Equations

$$\frac{d[\text{ATP}]}{dt}=v_{ATP}+v_{8.5}+2.5\cdot (v_7+ v_{8.3}+ v_{8.4}+v_{9.1-}-v_{9.1+})-1.5\cdot v_{9.3}$$ $$\frac{d[\text{OAA}]}{dt}=v_{OAA}+v_4-v_{8.1}-v_{9.1+}+v_{9.1-}$$ $$\frac{d[\text{Pyr}]}{dt}=v_{Pyr}+v_4-v_7$$ $$\frac{d[\text{Ac}_\text{{CoA}}]}{dt}=v_{Ac_{CoA}}+v_5+v_6+v_7-v_{8.1}$$ $$\frac{d[\text{Cit}]}{dt}=v_{Cit}+v_{8.1}-v_{8.2+}+v_{8.2-}$$ $$\frac{d[\text{Icit}]}{dt}=v_{Icit}+v_{8.2+}-v_{8.2-}-v_{8.3}$$ $$\frac{d[\alpha \text{KG}]}{dt}=v_{\alpha KG}+v_{8.3}-v_{8.4}$$ $$\frac{d[\text{Succinyl}_\text{CoA}]}{dt}=v_{Succinyl_{CoA}}+v_{8.4}-v_{8.5}$$ $$\frac{d[\text{Succ}]}{dt}=v_{Succ}+v_5+v_6+v_{8.5}+v_{9.3}$$ $$\frac{d[\text{Fum}]}{dt}=v_{Fum}+v_{9.2+}-v_{9.2-}-v_{9.3}$$ $$\frac{d[\text{Mal}]}{dt}=v_{Mal}+v_{9.1+}-v_{9.1-}-v_{9.2+}+v_{9.2-}$$ $$\frac{d[\text{CoA}]}{dt}=v_{CoA}-v_5-v_6-v_7+v_{8.1}-v_{8.4}+v_{8.5}$$

5.Result Presentation:

VI.Succinate Extraction

1.Model Variables

2.Metabolic Reactions

Reaction 11: Extraction of Succinate

$$\begin{align}Succ \rightarrow Produc\end{align}$$ Since we have introduced a plasmid containing the succinate efflux carrier protein gene into the engineered bacteria, we assume that the produced Succinate is stably and quickly excreted from the engineered bacteria and extracted in the system. The extraction rate is analogous to a first-order reaction rate formula. Meanwhile, since the bacterial volume only accounts for a small part of the total liquid volume, the Succinate concentration is multiplied by the coefficient rt to represent the Succinate concentration in the fermentation system. $$\begin{align}v_{11}=K_{Succ}\cdot [Succ]\cdot rt\end{align}$$

3.Parameter List:

4.Ordinary Differential Equations:

$$\frac{d[\text{Succ}]}{dt}=v_{Succ}+v_5+v_6+v_{8.5}+v_{9.3}-v_{11}$$ $$\frac{d[\text{Produc}]}{dt}=v_{11}$$

5.Result Presentation:

VII.Model Conclusions and Guidance for Wet Lab Experiments:

1. The rate of laccase decomposition of phenolic lignin is faster than that of non-phenolic lignin.
2. Syringate has the highest proportion among lignin decomposition products, but its further metabolism is more difficult.
3. After further metabolism, most decomposition products enter the TCA cycle via oxaloacetate and pyruvate, while a small part enters via Acetyl-CoA and directly generates the product Succinate.
4. For the broken TCA cycle, Succinate is mainly produced through the equilibrium reaction in branch two (oxaloacetate, malate, fumarate) consuming energy.
5. The production rate of Succinate is approximately 7mmol/L * min, where /L corresponds to per liter of fermentation system.
6. ATP is consumed at a rate of about 20mmol/L * min, where /L corresponds to per liter of engineered bacterial cytoplasm, requiring supply of appropriate energy substances to maintain fermentation.
7. Based on the rate bottlenecks revealed by the model, methods such as increasing the content of Syringate degradation enzymes, TCA cycle branch two enzymes, or constructing pathways to recycle Acetyl-CoA can be considered to further improve the system’s Succinate production.

VIII.Complete List of Variables, Parameters, and Ordinary Differential Equations:

IX.Complete MATLAB code：

clear
clc

%%

first=[0.700;0.0008;0; 0;0;0; 0;0;0;
       0;0;0; 0;0;0; 0;0;-100; 0];

% 1 :Lp
% 2 :Lnp
% 3 :Mox
% 4 :Syr
% 5 :Van
% 6 :pHB
% 7 :ATP
% 8 :OAA
% 9 :Pyr
% 10:AcCoA
% 11:Cit
% 12:Icit
% 13:akG
% 14:SucCoA
% 15:Succ
% 16:Fum
% 17:Mal
% 18:CoA
% 19:Produc

function d=fuc(t,y)

Etotal  =2;
Mtotal  =1;
Kp      =0.8;
Km      =0.32;
Vp      =0.025;
Vm      =0.25;
k3      =60;
a       =0.4;
b       =0.5;
c       =0.1;
d       =0.85;
e       =0.15;
f       =0;
rt      =0.01;

KSyr    =15;
VSyr    =10;
KVan    =15;
VVan    =100;
KpHB    =15;
VpHB    =50;

KPDHcPyr    =0.037;
KPDHcCoA    =0;
VPDHc       =5.62;
KCSAcCoA    =0.0013;
KCSOAAm     =25.5;
KCSOAAs     =140;
VCS         =0.108;%%%%%%%
KACO1       =11;
VACO1       =2.088;
KACO0       =0.015;
VACO0       =3.045;
KIDH        =0.0114;
VIDH        =34.56;
KaKGCoA     =0;
KaKGaKG     =0.22;
VaKG        =8.738;
KSUCL       =0.06;
VSUCL       =0.0138;
KMDH1       =0.049;
VMDH1       =137.16;
KMDH0       =2.6;
VMDH0       =32;
KFH1        =1.1;
VFH1        =23.777;
KFH0        =0.15;
VFH0        =2.6598;
KFRD        =0.028;
VFRD        =426.9;
KPC         =1.36;
VPC         =89.25;
vATP        =0;
vOAA        =0;
vPyr        =0;
vAcCoA      =-0.9*y(10);
vCit        =0;
vIcit       =0;
vaKG        =0;
vSucCoA     =0;
vSucc       =0;
vFum        =0;
vMal        =0;
vCoA        =0;
KSucc       =0.5;

v1=Vp*y(1)/(Kp*(1+(Mtotal-y(3))/Km)+y(1));
v2=Vm*(Mtotal-y(3))/(Km*(1+y(1)/Kp)+Mtotal-y(3));
v3=k3*y(3)*y(2);

v4=VSyr*y(4)/(KSyr+y(4));
v5=VVan*y(5)/(KVan+y(5));
v6=VpHB*y(6)/(KpHB+y(6));
v7=VPDHc*y(9)*y(18)/( KPDHcPyr*y(18)+ KPDHcCoA*y(9)+ y(9)*y(18) );

v81=VCS*y(8)*y(10)/( KCSAcCoA*y(8)+ KCSOAAm*y(10)+ y(8)*y(10)+ KCSOAAs*KCSAcCoA );
v821=VACO1*y(11)/(KACO1+y(11));
v820=VACO0*y(12)/(KACO0+y(12));
v83=VIDH*y(12)/(KIDH+y(12));
v84=VaKG*y(13)*y(18)/( KaKGCoA*y(13)+ KaKGaKG*y(18)+ y(13)*y(18) );
v85=VSUCL*y(14)/(KSUCL+y(14));

v911=VMDH1*y(8)/(KMDH1+y(8));
v910=VMDH0*y(17)/(KMDH0+y(17));
v921=VFH1*y(17)/(KFH1+y(17));
v920=VFH0*y(16)/(KFH0+y(16));
v93=VFRD*y(16)/(KFRD+y(16));

v10=VPC*y(9)/(KPC+y(9));

v11=KSucc*y(15)*rt;

Lp=-v1+0.01;
Lnp=-v3+0.05;
Mox=v2-v3;
Syr=a*v1+d*v3-v4*rt;
Van=b*v1+e*v3-v5*rt;
pHB=c*v1+f*v3-v6*rt;
ATP=vATP+v85+2.5*(v7+v83+v84+v910-v911)-1.5*v93-v10;
OAA=vOAA+v4-v81-v911+v910+v10;
Pyr=vPyr+v4-v7-v10;
AcCoA=vAcCoA+v5+v6+v7-v81;
Cit=vCit+v81-v821+v820;
Icit=vIcit+v821-v820-v83;
aKG=vaKG+v83-v84;
SucCoA=vSucCoA+v84-v85;
Succ=vSucc+v5+v6+v85+v93-v11;
Fum=vFum+v921-v920-v93;
Mal=vMal+v911-v910-v921+v920;
CoA=vCoA-v5-v6-v7+v81-v84+v85;
Produc=v11;

d=[Lp;Lnp;Mox;Syr;Van;pHB;ATP;OAA;
    Pyr;AcCoA;Cit;Icit;aKG;SucCoA;
    Succ;Fum;Mal;CoA;Produc];
end

options = odeset('NonNegative', [1,2,3,4,5,6,8,9,10,11,12,13,14,15,16,17,19]);

[t0,y0] = ode45(@fuc,[0 900],first,options);

t=t0(1:10000:end, :);
y=y0(1:10000:end, :);

%% Obtain the process of variable changes
Etotal  =2;
Mtotal  =1;
Kp      =0.8;
Km      =0.05;
Vp      =0.025;
Vm      =0.25;
k3      =60;
a       =0.4;
b       =0.5;
c       =0.1;
d       =0.85;
e       =0.15;
f       =0;
rt      =0.01;

KSyr    =15;
VSyr    =10;
KVan    =15;
VVan    =100;
KpHB    =15;
VpHB    =50;

KPDHcPyr    =0.037;
KPDHcCoA    =0;
VPDHc       =5.62;
KCSAcCoA    =0.0013;
KCSOAAm     =25.5;
KCSOAAs     =140;
VCS         =0.108;
KACO1       =11;
VACO1       =2.088;
KACO0       =0.015;
VACO0       =3.045;
KIDH        =0.0114;
VIDH        =34.56;
KaKGCoA     =0;
KaKGaKG     =0.22;
VaKG        =8.738;
KSUCL       =0.06;
VSUCL       =0.0138;
KMDH1       =0.049;
VMDH1       =137.16;
KMDH0       =2.6;
VMDH0       =32;
KFH1        =1.1;
VFH1        =23.777;
KFH0        =0.15;
VFH0        =2.6598;
KFRD        =0.028;
VFRD        =426.9;
KPC         =1.36;
VPC         =89.25;
vATP        =0;
vOAA        =0;
vPyr        =0;
vAcCoA      =-0.9*y(10);
vCit        =0;
vIcit       =0;
vaKG        =0;
vSucCoA     =0;
vSucc       =0;
vFum        =0;
vMal        =0;
vCoA        =100;
KSucc       =0.5;

v1=Vp.*y(:,1)./(Kp*(1+(Mtotal-y(:,3))/Km)+y(:,1));
v2=Vm*(Mtotal-y(:,3))./(Km*(1+y(:,1)/Kp)+Mtotal-y(:,3));
v3=k3.*y(:,3).*y(:,2);

v4=VSyr.*y(:,4)./(KSyr+y(:,4));
v5=VVan.*y(:,5)./(KVan+y(:,5));
v6=VpHB.*y(:,6)./(KpHB+y(:,6));
v7=VPDHc.*y(:,9).*y(:,18)./( KPDHcPyr.*y(:,18)+ KPDHcCoA.*y(:,9)+ y(:,9).*y(:,18) );

v81=VCS.*y(:,8).*y(:,10)./( KCSAcCoA.*y(:,8)+ KCSOAAm.*y(:,10)+ y(:,8).*y(:,10)+ KCSOAAs*KCSAcCoA );
v821=VACO1.*y(:,11)./(KACO1+y(:,11));
v820=VACO0.*y(:,12)./(KACO0+y(:,12));
v83=VIDH.*y(:,12)./(KIDH+y(:,12));
v84=VaKG.*y(:,13).*y(:,18)./( KaKGCoA.*y(:,13)+ KaKGaKG.*y(:,18)+ y(:,13).*y(:,18) );
v85=VSUCL.*y(:,14)./(KSUCL+y(:,14));

v911=VMDH1.*y(:,8)./(KMDH1+y(:,8));
v910=VMDH0.*y(:,17)./(KMDH0+y(:,17));
v921=VFH1.*y(:,17)./(KFH1+y(:,17));
v920=VFH0.*y(:,16)./(KFH0+y(:,16));
v93=VFRD.*y(:,16)./(KFRD+y(:,16));

v10=VPC.*y(:,9)./(KPC+y(:,9));

v11=KSucc.*y(:,15)*rt;

Lp=-v1+0.01;
Lnp=-v3+0.05;
Mox=v2-v3;
Syr=a*v1+d*v3-v4*rt;
Van=b*v1+e*v3-v5*rt;
pHB=c*v1+f*v3-v6*rt;
ATP=vATP+v85+2.5*(v7+v83+v84+v910-v911)-1.5*v93-v10;
OAA=vOAA+v4-v81-v911+v910+v10;
Pyr=vPyr+v4-v7-v10;
AcCoA=vAcCoA+v5+v6+v7-v81;
Cit=vCit+v81-v821+v820;
Icit=vIcit+v821-v820-v83;
aKG=vaKG+v83-v84;
SucCoA=vSucCoA+v84-v85;
Succ=vSucc+v5+v6+v85+v93-v11;
Fum0=vFum+v921-v920-v93;
Fum = movmean(Fum0, 3);
Mal=vMal+v911-v910-v921+v920;
CoA=vCoA-v5-v6-v7+v81-v84+v85;
Produc=v11;


%% Figure sum:
figure;
vP=v1+v3;
vin=vOAA+v4+vSucc+v5+v6;
vout=v11;
vA=vATP+v85+2.5*(v7+v83+v84+v910-v911)-1.5*v93-v10;
vP=normalize(vP,'range');
vin=normalize(vin,'range');
vout=normalize(vout,'range');
vA=normalize(-vA,'range');
hold on;
plot(t,vP,'LineWidth',2);
plot(t,vin,'LineWidth',2);
plot(t,vout,'LineWidth',2);
plot(t,vA,'LineWidth',2);
legend('Lignin','Metabolic Entry into TCA','Succinate production','ATP consumption');
title('Rate after normalize');
xlabel('Time /min');
ylabel('Velocity');
hold off;

%% Other plottings have been omitted

FBA Model: Model for Succinate Yield Prediction

The model aims at predicting the succinate production with the dual-microbes system in industry. According to the model, each 1 g biomass can produce 5.343 g succinate each hour.

Principle & Design

The model is based on the COnstraint-Based Reconstruction and Analysis (COBRA) methods which has been wildly used to model the metabolic networks in both prokaryotes and eukaryotes on the genome-scale. This model does not require difficult-to-measure kinetic parameters. Instead, it lists the stoichiometric coefficients of each metabolic reaction in the form of a numerical matrix as constraints. Each row represents a metabolite and each column represents a reaction. Under steady state, the flux of each reaction is given by Sv = 0, which defines a linear equation system. Then, linear programming is used to determine the flux distribution, which is optimized or minimized by maximizing or minimizing the objective function within the allowable flux space defined by the mass balance equation and the reaction boundaries^[24].

With the model, the hourly succinate production and biomass by 1 g biomass can be roughly calculated. Then, with integration, we might get the total production required in the industrial design. With proper construction of the matebolism network, it might be a method to offer opimized prediction in absence of experimental data on a large scale.

Realization

The constrcution of the model started from a exist model downloaded from BIGG for wild type P.Putida. Although exist model for T.reesei was unvailable, the core metabolism was similar. Therefore, this model served as the base for the dual-microbes system and extra reactions of both engineered microbses were added to it later^[25].

# import model of wild P.Putida from IBGG
import cobra
model = cobra.io.read_sbml_model('iJN746.xml')

# create a new reaction for succinate output
from cobra import Model, Reaction, Metabolite

reaction = Reaction('SUCCt')
reaction.name = 'succinate transport'
reaction.subsystem = 'succinate transport'
reaction.lower_bound = 0.0
reaction.upper_bound = 999999.9

succ_c = Metabolite(
    'succ_c',
    formula = 'C4H4O4',
    name='Succinate',
    compartment='c')
succ_e = Metabolite(
    'succ_e',
    formula = 'C4H4O4',
    name='Succinate',
    compartment='c')

reaction.add_metabolites({
    succ_c: -1.0,
    succ_e: 1.0,
})  

reaction.reaction

model.add_reactions([reaction])

Table.1. the Reactions added in the final model

Then, the bounds of the reactions controled by the genes to be blocked were set as “(0,0)”. In this model, the reactions where succinate is transferd to fumarate was blocked. After that, the medium was set according to the composition of the broth after the first fermentation.

# block transformation of succinate to fumarate
model.reactions.get_by_id('SUCD1').bounds = (-9999.9,0)
model.reactions.get_by_id('SUCDi').bounds = (0,0)

Table.2. the medium composition

# change your medium
medium['EX_glc__D_e'] = 1.01
model.medium = medium

Result & Discussion

Solve the model with the optimizing projective as the reaction where succinate is produced in the Kreb’s cycle, and then the uptake and secrection can be print out as listed below. These data were used in the prediction of succinate production which can be found on the page for hardware.

# solve the model and output
model.objective = "SUCOAS"
model.optimize()
model.summary()

The model also has something to be improved or noticed:

• The basic model of wild types isn’t the same as what they should be in live microbes. There’s surely some reactions missing.
• The time and the conversion rate of each reaction isn’t considered.
• A steady state is assumed so that it’s not the same as the real fermentation process.
• The model can only provide the ideal case, and for industry, a large-scale experiment is always neccesary.

Dual-Species Microbial System: Mathematical Modeling and Simulation

Model Description

This document describes a mathematical model of a dual-species microbial system consisting of Trichoderma reesei and Pseudomonas putida. The model captures the synergistic interaction where T. reesei degrades lignin to produce aromatic compounds that serve as growth substrates for P. putida, which in turn produces succinate as a valuable metabolic product^[26].

The system represents a typical syntrophic relationship in microbial communities, where the metabolic activities of one species create favorable conditions for another. This model provides a quantitative framework for understanding and optimizing such cooperative systems for industrial applications like lignocellulosic biomass conversion^[27].

Basic Assumptions

1. Spatial Homogeneity: The system is well-mixed with uniform distribution of all components.

2. Monod Kinetics: Microbial growth follows Monod kinetics, where growth rate depends on limiting substrate concentration.

3. Enzyme-Mediated Reactions: Lignin degradation is catalyzed by laccase enzymes produced by T. reesei, following Michaelis-Menten kinetics.

4. Product Inhibition: High concentrations of metabolic products (particularly succinate) can inhibit microbial growth.

5. Mass Conservation: Substrate consumption, biomass production, and product formation are connected through yield coefficients.

6. Dynamic Enzyme Production: Enzyme synthesis is regulated by substrate induction and follows first-order kinetics.

7. Non-growth Associated Maintenance: Microbial populations experience natural decay rates.

Model Equations

State Variables