Modification of GEM for EcN

In this project, we try to introduce a pathway for utilizing erythritol into EcN in order to improve the colonization of EcN in the ileal, also to secrete beneficial metabolites such as PQQ at the same time.

In this section, we will describe how we provide theoretical support for the feasibility of our subject based on GEM, through the adjustment of the metabolic network, the combination with enzyme kinetics, and the optimization of details.

The framework of our model, the construction and modification of each reaction and the coupling between the components will be described below. The units of the models are volume (L), time (hr), concentration (mM or mmol/L), and biomass (g/L).

Before reading this section, we would like you to gain some knowledge about GEM.

Modification of EcN model in anoxic environment

The goal of the project is to colonize the ileal, so we also need to modify EcN under anoxic conditions. (see codes below)

It should be noted that the modification we mean the adjustment of the bound of certain reactions, for example, the limit of the oxygen absorption reaction will be adjusted from 10 mmol/(gDw*h) to 0.6 mmol/(gDw*h). These numerical changes are based on data from Kim, D., Kim, Y., Yoon, S.H., 2021..

EX_o2_e = EcN_p.reactions.get_by_id('EX_o2_e')
EX_o2_e.lower_bound = -0.6

DTARTD = EcN_p.reactions.get_by_id('DTARTD')
DTARTD.lower_bound = -1000
DTARTD.upper_bound = 1000

LCARS = EcN_p.reactions.get_by_id('LCARS')
LCARS.lower_bound = -1000
LCARS.upper_bound = 1000

SUCTARTtpp = EcN_p.reactions.get_by_id('SUCTARTtpp')
SUCTARTtpp.lower_bound = -1000
SUCTARTtpp.upper_bound = 1000

TARTD = EcN_p.reactions.get_by_id('TARTD')
TARTD.lower_bound = -1000
TARTD.upper_bound = 1000

TARTRt7pp = EcN_p.reactions.get_by_id('TARTRt7pp')
TARTRt7pp.lower_bound = -1000
TARTRt7pp.upper_bound = 1000

ARGAGMt7pp = EcN_p.reactions.get_by_id('ARGAGMt7pp')
ARGAGMt7pp.lower_bound = -1000
ARGAGMt7pp.upper_bound = 1000

FHL = EcN_p.reactions.get_by_id('FHL')
FHL.lower_bound = -1000
FHL.upper_bound = 1000

EX_gal_e = EcN_p.reactions.get_by_id('EX_gal_e')
EX_gal_e.lower_bound = -20

EX_acnam_e = EcN_p.reactions.get_by_id('EX_acnam_e')
EX_acnam_e.lower_bound = -20

EX_acgam_e = EcN_p.reactions.get_by_id('EX_acgam_e')
EX_acgam_e.lower_bound = -20

EX_glc__D_e = EcN_p.reactions.get_by_id('EX_glc__D_e')
EX_glc__D_e.lower_bound = -20

EX_acgal_e = EcN_p.reactions.get_by_id('EX_acgal_e')
EX_acgal_e.lower_bound = -20

EX_gal_e = EcN_p.reactions.get_by_id('EX_gal_e')
EX_gal_e.lower_bound = -20

EX_glcn_e = EcN_p.reactions.get_by_id('EX_glcn_e')
EX_glcn_e.lower_bound = -20

EX_srb__L_e = EcN_p.reactions.get_by_id('EX_srb__L_e')
EX_srb__L_e.lower_bound = -20

EX_drib_e = EcN_p.reactions.get_by_id('EX_drib_e')
EX_drib_e.lower_bound = -20

EX_fuc__L_e = EcN_p.reactions.get_by_id('EX_fuc__L_e')
EX_fuc__L_e.lower_bound = -20

EX_glyc_e = EcN_p.reactions.get_by_id('EX_glyc_e')
EX_glyc_e.lower_bound = -20

EX_akg_e = EcN_p.reactions.get_by_id('EX_akg_e')
EX_akg_e.lower_bound = -20

EX_lac__L_e = EcN_p.reactions.get_by_id('EX_lac__L_e')
EX_lac__L_e.lower_bound = -20
                        

Since the anoxic transformation is the most fundamental part, in the following description, "EcN before transformation" refers to the EcN model that has only been modified by an anoxic environment.

Introduction of metabolic pathways

According to the requirements of this year's project, we will transform the EcN model in the following three aspects:

1. Secretion of PQQ:

There is no PQQ in the original GEM of EcN, so we needed to introduce PQQ into the GEM of EcN first.

PQQ secretion consists of two steps. The first step is the intracellular production of PQQ reaction (Stites, T.E. et al., 2000):

\begin{equation} Glu+Tyr+3O2\longrightarrow 2H_2O_2+H_2O+PQQ_e \label{(1)} \end{equation}

Glu = EcN.metabolites.get_by_id('glu__L_c')
Tyr = EcN.metabolites.get_by_id('tyr__L_c')
O2 = EcN.metabolites.get_by_id('o2_c')
H2O2 = EcN.metabolites.get_by_id('h2o2_c')
H2O = EcN.metabolites.get_by_id('h2o_c')
PQQ = Metabolite(
    'PQQ_c',
    formula='C14H6N2O8',
    name='pyrroloquinoline quinone',
    compartment='c'
)

PQQ_syn = Reaction('PQQS')
PQQ_syn.add_metabolites({
    Glu: -1,
    Tyr: -1,
    O2: -3,
    H2O2: 2,
    H2O: 1,
    PQQ: 1
})

EX_PQQ = Reaction('EX_PQQ')
EX_PQQ.add_metabolites({
    PQQ: -1
})
EcN.add_reactions({PQQ_syn, EX_PQQ})
                        

After introducing the first step of PQQ secretion into GEM, we found that we also needed to describe the process of secretion of PQQ, which is the second step. Considering the bit rough of GEM model, we simulated the second step by enzyme kinetics (see in "Construction of ODEs system")

2.Oxidation of erythritol into erythrulose

This part is divided into the following three processes (Zou, X. et al., 2017): (1) Erythritol is absorbed into EcN under the action of ABC transporter (2) erythritol is converted to L-Erythrulose, (3) L-Erythrulose is converted to L-Erythrose:

\begin{equation} Erythritol_{e}\longrightarrow Erythritol_{c} \label{(2)} \end{equation} \begin{equation} Erythritol_{c} \longrightarrow L-Erythrulose \label{(3)} \end{equation} \begin{equation} L-Erythrulose \longrightarrow L-Erythrose \label{(4)} \end{equation}

Since the above three substances (Erythritol, L-Erythrulose and L-Erythrose) are not present in the GEM of EcN, in order to minimize errors and to simulate realistic conditions, we simulated this whole part using enzyme kinetics (see "Construction of ODEs system"). ).

3. Erythritol enters the pentose phosphate pathway (PPP):

The first step of erythritol entry into the PPP pathway is still the aspiration of erythritol into EcN under the action of ABC transporter, a step shared with "erythritol production from erythritol".

After inhalation of erythritol, the utilization of erythritol follows the following steps respectively (Barbier, T. et al., 2014):

\begin{equation} Erythritol_{c} \longrightarrow L-Erythritol-4-P \label{(5)} \end{equation} \begin{equation} L-Erythritol-4-P \longrightarrow L-3-tetrulose-4-P \label{(6)} \end{equation} \begin{equation} L-3-tetrulose-4-P \longrightarrow D-3-tetrulose-4-P \label{(7)} \end{equation} \begin{equation} D-3-tetrulose-4-P \longrightarrow D-erythrulose-4-P \label{(8)} \end{equation} \begin{equation} D-Erythrulose-4-P \longrightarrow D-Erythrose-4-P \label{(9)} \end{equation} \begin{equation} D-Erythrose-4-P \longrightarrow Fructose-6-P \label{(10)} \end{equation} \begin{equation} D-Erythrose \longrightarrow Glyceraldehyde-3-P \label{(11)} \end{equation}

The reactions represented by Eqs. 10-11 are already present in the GEM of EcN, so all we need to do is to construct the enzyme kinetic equations corresponding to Eqs. 5-9 and to dynamically couple the enzyme kinetic equations corresponding to Eqs. 7-9 to the GEM (see "Coupling of GEM and Michaelis equations").

L_tet=Metabolite(
    'L_tet4p',
    formula='C4H8O7P',
    name='L-tetrulose-4-p',
    compartment='c'
)
D_tet=Metabolite(
    'D_tet4p',
    formula='C4H8O7P',
    name='D-tetrulose-4-p',
    compartment='c'
)
D_erythrulose4p=Metabolite(
    'D_eryth4p',
    formula='C4H8O7P',
    name='D-erythrulose-4-p',
    compartment='c'
)


EryC = Reaction('EryC')
EryC.name = 'iso of tetrulose '
EryC.subsystem = 'Cell Envelope Biosynthesis'
EryC.lower_bound = -1000. 
EryC.upper_bound = 1000.  


EryC.add_metabolites({
    D_tet: 1.0
})
EryC


TpiA2 = Reaction('TpiA2')
TpiA2.name = 'tetrulose to erythrulose '
TpiA2.subsystem = 'Cell Envelope Biosynthesis'
TpiA2.lower_bound = -1000.  # This is the default
TpiA2.upper_bound = 1000.  # This is the default

TpiA2.add_metabolites({
    D_tet: -1.0,
    D_erythrulose4p: 1.0
})
TpiA2


D_erythrose4p=EcN.metabolites.get_by_id('e4p_c')
RpiB = Reaction('RpiB')
RpiB.name = 'erythrulose to erythrose'
RpiB.subsystem = 'Cell Envelope Biosynthesis'
RpiB.lower_bound = -1000.  # This is the default
RpiB.upper_bound = 1000.  # This is the default

RpiB.add_metabolites({
    D_erythrulose4p: -1.0,
    D_erythrose4p: 1.0
})
RpiB
EcN.add_reactions({RpiB, TpiA2, EryC})
                        

After this series of reactions, the obtained Fructose-6-P and Glyceraldehyde-3-P could enter the PPP pathway and participate in the construction of EcN biomass, thus enhancing the biomass of EcN.

After introducing equation (1) into GEM, we found that GEM was not friendly to the intracellular and extracellular exchange reactions of the newly introduced substances. Due to the features of GEM, if the secretion reaction of PQQ is introduced directly into the GEM of EcN, then all the PQQ will be secreted to the extracellular as soon it is produced, which is not practical and not what we want.

In addition, we also needed to consider the concentration conversion of intracellular and extracellular substances due to the difference in volume, assuming that the cells are uniformly distributed in the environment, after converting the units, the result is shown below.

\begin{equation} {dCout}=\frac{Biomass*Vout*V{single}}{W_{single}*Vout}dCin \label{(12)} \end{equation}

Thus, we constructed the enzyme kinetic equation for PQQ secretion:

\begin{equation} \frac{dPQQ_{e}}{dt}=-\frac{Biomass*Vout*V_{single}}{Wsingle*Vout}*\frac{Kcat_{PQQ_e}*Enzyme*PQQ_{e}}{Km_{PQQ_{e}}+PQQ_{e}} \label{(13)} \end{equation}

Besides, we also constructed the enzyme kinetic equation related to Eqs. 2-6:

\begin{equation} \frac{dErythritol_e}{dt}= -\frac{Kcat_{Erythritol_e}*Enzyme*Erythritol_e}{Erythritol_e+Km_{Erythritol_e}}*\frac{Biomass*Vout*V_{single}}{Wsingle*Vout} \label{(14)} \end{equation} \begin{equation} \frac{dErythritol_c}{dt}=\frac{dErythritol_e}{dt}-\frac{Kcat_{sld}*Enzyme*Erythritol_c}{Km_{sld}+Erythritol_c}-\frac{Km_{EryA}*Enzyme*Erythritol}{Km_{EryA}+Erythritol_c} \label{(15)} \end{equation} \begin{equation} \frac{dL-Erythrulose}{dt}=\frac{Km_{sld}*Enzyme*Erythritol_c}{Km_{sld}+Erythritol_c}-\frac{Kcat_{L-RI}*Enzyme*Erythrulose}{Km_{L-RI}+Erythrulose} \label{(16)} \end{equation} \begin{equation} \frac{L-Erythrose}{dt}=\frac{Kcat_{L-RI}*Enryme*Erythrulose}{Km_{L-RI}+Erythrulose} \label{(17)} \end{equation} \begin{equation} \frac{L-Erythritol-4-P}{dt}= \frac{Kcat_{EryA}*Enzyme*Erythritol}{Km_{EryA}+Erythlose}-\frac{Kcat_{EryB}*Enzyme*L-Erythritol-4-P}{Km_{EryB}+L-Erythritol-4-P} \label{(18)} \end{equation}

In addition, considering the substance and production of Eqs. 7-9 are the same in molecular formula. For parametric considerations, we merge Eqs. 7-9 into a one-step enzyme kinetic reaction.

Coupling of GEM and Michaelis equations

After constructing the GEM and Michaelis equation models, we now need to combine these two parts, which is a crucial step for the part of "Erythritol into PPP pathway".

The GEM and the Michaelis equation model are two separate models, and we would use the dFBA model for our dynamic simulations, so we planned to incorporate the results of the Michaelis equation into the GEM to achieve a coupling effect:

We created a new pathway in the GEM and we wanted to couple the two models by changing the "flux" of this pathway:

\begin{equation} None \longrightarrow D-tetrulose \label{(19)} \end{equation}

What we want to do is quantitative, but the unit of "flux" in GEM (mmol/(g*h)) is not the same as the unit of the enzyme kinetic equation (mmol/(L*h)), so we have the following equation for the conversion:

\begin{equation} V_{intracellular\; metabolites}=flux*Biomass*\frac{V_{out}}{V_{cells}} \label{(20)} \end{equation} \begin{equation} V_{cells}=\frac{Biomass}{W_{single cell}}*V_{out}*V_{single cell} \label{(21)} \end{equation}

We queried some parameters of EcN: Volume of cell ≈ 1*10-15L, Cell dry weight ≈ 3*10-13g (Milo, R., Phillips, R., 2016). After bringing in the parameters, we obtained the "flux to rate conversion equation" that fits EcN:

\begin{equation} V_{intracellular\; metabolites}=300flux \label{(22)} \end{equation}

\begin{equation} flux_{L-3-tetrulose}=\frac{Kcat*Enzyme*L-3-tetrulose-4-P}{300(Km_{L-3-tetrulose-4-P}+L-3-tetrulose-4-P)} \label{(23)} \end{equation} \begin{equation} \frac{d{L-3-tetrulose}}{dt}=\frac{Kcat_{EryB}*Enzyme*L-Erythritol-4-P}{Km_{EryB}+L-Erythritol-4-P}-300flux_{L-3-tetrulose} \label{(24)} \end{equation}

At the end of the model building, after converting the units, here is the rate of change of Biomass (g/L) in this model:

\begin{equation} \frac{dBiomass}{dt}=flux_{Biomass}*Biomass-D*Biomass \label{(25)} \end{equation}

After finishing the coupling, the initial construction of the model is completed, and various computational simulations can be performed by this EcN model.

Verification of the effectiveness of the modification by dFBA

In the previous section, we have built the modified EcN model by GEM and enzyme kinetics model, which is a work worth celebrating, but we still needed to do some validation of the modified model by dFBA model simulation:

The change rate of key substances in the model:

In order to verify the accuracy of the model, we first need to observe whether the key substances involved in the model are consistent with the cognitive and experimental rules, so we simulated the changes of each key substance by dFBA when Erythritol (10mm) was given in addition to the nutrients required for basic survival:

From the above figure, it is clear that the pattern of change of each substance is in line with our conjecture and close to the reality, which proves that our modification is realistic.

Secretion of PQQ:

Next, we modelled for the secretion of PQQ in the modified EcN model with sufficient carbon sources such as glucose and Erythritol to verify the validity of our model. The results of the simulations are as follows:

According to the graph, most of the PQQ is exported, which is basically consistent with the literature that 98% of PQQ will be secreted to the extracellular. (Velterop, J.S., et al., 1995). Also, the concentration of intracellular PQQ will maintain a certain level, which fit the feature of a coenzyme.

Comparison between before and after modification:

Another key part of the model is the uptake, conversion and utilization of Erythritol. The most important pathway for the uptake of Erythritol into the intracellular compartment is the progressive conversion of Erythritol into various compounds and then into the PPP pathway, thus enhancing the Biomass of EcN and achieving our goal of enhancing EcN colonization.

So we simulated the co-culture of EcN model before and after modification under the set initial conditions (10mm Erythritol and 10mm glucose):

From the figure above, it is observed that the modified EcN has stronger growth and material utilization ability than the unmodified EcN in the early stage of bacterial growth, and this advantage is more obvious in the late stage of material deprivation.

We can conclude that the modified EcN is able to increase its growth using erythritol, and the modified EcN is superior to the unmodified EcN in terms of viability! The modified EcN has a competitive advantage in the genus.

Sensitivity analysis

Sensitivity analysis of enzyme concentration

Enzyme concentration is a crucial parameter in the whole model, and is the most important part of the sensitivity analysis. Therefore, we will set a gradient concentration of Enzyme and observe the change of Biomass in different cases. The results of the simulations are as follows:

The results are very satisfactory, and we can see that the change of enzyme concentration has a relatively small effect on the trend and size of Biomass, which further proves the correctness and robustness of the model.

Sensitivity analysis of enzyme kinetics parameters

Since in the previous model construction, we simplified Eqs. 7-9 as a one-step response and made an estimate for this parameter, in order to verify the correctness and robustness of the model, we also needed to do a sensitivity analysis for Kcat_L3t4p. The effect of Kcat_L3t4p on the model was negligible.

In addition, we also did a sensitivity analysis on Kcat_EryB, and we found that the effect of EryB on the model is also extremely small, what does this mean? Let's start our next section from this point.

Combination of model and experiment and implications

Which Ery gene is more critical in the erythritol metabolic pathway? erythritol import, EryA, EryB or EryC?

Observing the effect of erythritol on the biomass of the modified EcN, we can divide the whole process into five steps, named erythritol import, EryA, EryB, EryC(EryC includes the process of Eqs. 7-9) respectively. We hope to guide the experiment with our model, for example, how these five steps relate to each other, how these five steps together affect the final biomass, and which of these five steps is more important.

The meaning of our results would be to change the activity of the target enzyme by regulating the expression of the corresponding gene. We investigate the relationship between these four steps, namely erythritol import, EryA ,EryB and EryC, as well as the effect on biomass by regulating enzyme concentrations (to simulate the enhancement of this gene)

We take biomass flux as our dependent variable and the enzyme concentrations of the four reactions as the independent variables, and the results are shown below.

By observing the results above, we can conclude that the degree of influence of the four steps on the whole process is ranked in erythritol import > EryC > EryA> EryB.

As Fig 9 shown, you can see the different roles of 4 genes play in the erythritol utilization process. If the ability of importing erythritol stay constant, no matter how the rest three genes express(the correlated width of pipe increase), the biomass will not change. Again, this conclusion applies to the relationship between EryC to EryB and EryC to EryA.

The first part of this result is very convincing in terms of its relevance to reality, and the second part of this result can also be used as a guide for the experiment:

The second half of this result can also provide some guidance to the experiment: When increase the Biomass by using Erythritol in the PPP pathway, the adequacy and importance of the EryA part is higher than that of the EryB part, which is consistent with the results of the sensitivity analysis of EryB in our previous paper. This is in line with the results of our sensitivity analysis of EryB in the previous section, thus increasing the credibility of the model.

Then, we summarize as follows: If we want to enhance the absorption and the use of Erythritol by EcN, it is necessary to enhance the intensity of gene expression in the part of erythritol import. Secondly, considering that the EryC part is composed of multiple steps and the complexity of the experiment, enhancing the gene expression intensity of the EryA step should be the second most necessary. The EryB part is not as influential as the previous parts. This conclusion reveals how the wet experiment part regulate the expression of each gene in order to work better.

Reference

Barbier, T., Collard, F., Zúñiga-Ripa, A., Moriyón, I., Godard, T., Becker, J., Wittmann, C., Van Schaftingen, E., Letesson, J.-J., 2014. Erythritol feeds the pentose phosphate pathway via three new isomerases leading to D-erythrose-4-phosphate in Brucella. Proc. Natl. Acad. Sci. U.S.A. 111, 17815–17820. https://doi.org/10.1073/pnas.1414622111

Burger C, Kessler C, Gruber S, Ehrenreich A, Liebl W, Weuster-Botz D. L-Erythrulose production with a multideletion strain of Gluconobacter oxydans. Appl Microbiol Biotechnol. 2019 Jun;103(11):4393-4404. doi: 10.1007/s00253-019-09824-w. Epub 2019 Apr 18. PMID: 31001743.

Hosie AH, Allaway D, Galloway CS, Dunsby HA, Poole PS. Rhizobium leguminosarum has a second general amino acid permease with unusually broad substrate specificity and high similarity to branched-chain amino acid transporters (Bra/LIV) of the ABC family. J Bacteriol. 2002 Aug;184(15):4071-80. doi: 10.1128/JB.184.15.4071-4080.2002. PMID: 12107123; PMCID: PMC135202.

Kim, D., Kim, Y., Yoon, S.H., 2021. Development of a Genome-Scale Metabolic Model and Phenome Analysis of the Probiotic Escherichia coli Strain Nissle 1917. IJMS 22, 2122. https://doi.org/10.3390/ijms22042122

Lillo AM, Tetzlaff CN, Sangari FJ, Cane DE. Functional expression and characterization of EryA, the erythritol kinase of Brucella abortus, and enzymatic synthesis of L-erythritol-4-phosphate. Bioorg Med Chem Lett. 2003 Feb 24;13(4):737-9. doi: 10.1016/s0960-894x(02)01032-6. PMID: 12639570.

Lodeiro A, Melgarejo A. Robustness in Escherichia coli glutamate and glutamine synthesis studied by a kinetic model. J Biol Phys. 2008 Apr;34(1-2):91-106. doi: 10.1007/s10867-008-9109-9. Epub 2008 Aug 30. PMID: 19669495; PMCID: PMC2577740.

Mahmood S, Iqbal MW, Zhang W, Mu W. A review on l-ribose isomerases for the biocatalytic production of l-ribose and l-ribulose. Food Res Int. 2021 Jul;145:110409. doi: 10.1016/j.foodres.2021.110409. Epub 2021 May 11. PMID: 34112412.

Milo, R., Phillips, R., 2016. Cell biology by the numbers. Garland Science, Taylor & Francis Group, New York, NY.

Schiller C, Fröhlich CP, Giessmann T, et al. Intestinal fluid volumes and transit of dosage forms as assessed by magnetic resonance imaging. Aliment Pharmacol Ther. 2005;22(10):971-979. doi:10.1111/j.1365-2036.2005.02683.x

Stites, T.E., Sih, T.R., Rucker, R.B., 2000. Synthesis of [14C]pyrroloquinoline quinone (PQQ) in E. coli using genes for PQQ synthesis from K. pneumoniae. Biochimica et Biophysica Acta (BBA) - General Subjects 1524, 247–252. https://doi.org/10.1016/S0304-4165(00)00166-5

Velterop, J.S., Sellink, E., Meulenberg, J.J., David, S., Bulder, I., Postma, P.W., 1995. Synthesis of pyrroloquinoline quinone in vivo and in vitro and detection of an intermediate in the biosynthetic pathway. J Bacteriol 177, 5088–5098. https://doi.org/10.1128/jb.177.17.5088-5098.1995

Wang L, Johnson ZL, Wasserman MR, Levring J, Chen J, Liu S. Characterization of the kinetic cycle of an ABC transporter by single-molecule and cryo-EM analyses. Elife. 2020 May 27;9:e56451. doi: 10.7554/eLife.56451.

Webb AJ, Hosie AH. A member of the second carbohydrate uptake subfamily of ATP-binding cassette transporters is responsible for ribonucleoside uptake in Streptococcus mutans. J Bacteriol. 2006 Dec;188(23):8005-12. doi: 10.1128/JB.01101-06. Epub 2006 Sep 22. PMID: 16997965; PMCID: PMC1698198.

Zou, X., Lin, J., Mao, X., Zhao, S., Ren, Y., 2017. Biosynthesis of L-Erythrose by Assembly of Two Key Enzymes in Gluconobacter oxydans. J. Agric. Food Chem. 65, 7721–7725. https://doi.org/10.1021/acs.jafc.7b02201