In the modeling part, we are mainly working on simulating the function of the fermentation system at different levels and serving suggestions for future implementation.

 

We first introduce the pathway of our project.

Figure 1 The biochemical pathway of the project. Escherichia coli Bl21 and Synechococcus PCC 7942 are two selected engineering microorganisms. There are five engineered E. coli Bl21, responsible for the degradation or detection of formaldehyde, degradation or detection of hydrogen sulfide, and sucrose hydrolysis, Synechococcus PCC 7942 for sucrose production. The figure follows the SBGN standard (Novère et al., 2009). Refer to our engineer part for more details.

 

The work consists of two parts, flux balance analysis (FBA) and diffusion model of alginate beads considering bacterial metabolism, specifically, the photosynthesis of Synechococcus PCC 7942, the formaldehyde degradation function, and sucrose hydrolysis by two types of E. coli.

 

  1. ODE system

The Production of FrmR

  1. The transcription process of FrmR,

\[ \frac{d[mRNA_{FrmR}]}{dt} = k_{TC1}DNA_{FrmR}-k_{DegM}[mRNA_{FrmR}] \]

where \(k_{TC1}\) is the transcription rate of of \(FrmR\) and \(k_{DegM}\) is the rate of mRNA degradation.

  1. The translation process of FrmR,

\[ \frac{d[FrmR]}{dt} = k_{TL1}[mRNA_{FrmR}] - k_{DegP}[FrmR] \]

where \(k_{TL1}\) is the rate of FrmR translation and \(k_{DegP}\) is the rate of protein degradation.

The Reaction of FrmR, PfrmR and HCHO

  1. Binding of FrmR protein to the PfrmR promoter region

\[ FrmR + PfrmR \Leftrightarrow FrmR:PfrmR \]

  1. Reaction of FrmR with formaldehyde leading to conformational change

\[FrmR + HCHO \rightarrow FrmR-HCHO\]

  1. Reverse reaction of FrmR-Formaldehyde

\[ FrmR-HCHO \rightarrow FrmR + HCHO \]

  1. Dissociation of FrmR:PfrmR complex due to conformational change

\[ FrmR:PfrmR \rightarrow FrmR + PfrmR \]

Based on the reactions above, we can build a system that represent the process.

About the concentration of \(FrmR\), we have

\[ \frac{d[FrmR]}{dt} = k_{Asso-FP}[FrmR][PfrmR] -k_{Asso-FH}[FrmR][HCHO]\] \[ + k_{Sepe-FH}[FrmR-HCHO] - k_{Sepe-FP}[FrmR:PfrmR] \]

About the concentration of \(PfrmR\), we have

\[ \frac{d[PfrmR]}{dt} = -k_{Asso-FP}[FrmR][PfrmR] + k_{Sepe-FP}[FrmR:PfrmR] \]

About the concentration of \(FrmR:PfrmR\), we have

\[ \frac{d[FrmR:PfrmR]}{dt} = k_{Asso-FP}[FrmR][PfrmR] - k_{Sepe-FP}[FrmR:PfrmR] \]

About the concentration of \(FrmR-HCHO\) , we have

\[ \frac{d[FrmR-HCHO]}{dt} = k_{Asso-FH}[FrmR][HCHO] - k_{Sepe-FH}[FrmR-HCHO] \]

The Production of Fluorescent Protein

  1. The transcription process of Fluorescent Protein “R”,

\[ \frac{d[mRNA_{R}]}{dt} = k_{TC2}[PfrmR] - k_{DegM}[mRNA_{R}]\]

where \(k_{TC2}\) is the transcription rate of \(PfrmR\).

  1. The translation process of “R”,

\[ \frac{d[R]}{dt} = k_{TL2}[PfrmR] - k_{DegP}[R] \]

where \(k_{TL2}\) is the translation rate of \(PfrmR\).

Formaldehyde degradation Equations

  1. \(HCHO_{air} \rightarrow HCHO_{in vivo}\), \(k_0\).

  2. \(HCHO_{in vivo} + GSH \rightarrow HMGS\), \(k_1\).

  3. \(HMGS + NAD^+ + adh \rightarrow mid_1 + NADH + S-formylglutathione + GSH + adh\), \(k_2\)

  4. \(S-formylglutathione + fghA \rightarrow mid_2 + fghA\), \(k_3\)

  5. \(mid_2 + fdh \rightarrow mid_3 + fdh\), \(k_4\)

  6. \(mid_3 \rightarrow CO_2{in vivo}\), \(k_5\)

  7. \(CO_2{in vivo} \rightarrow CO_2{air}\), \(k_6\)

System

\[\frac{d[HCHO_{in\ vivo}]}{dt} = k_0*[HCHO_{air}] - k_1*[HCHO_{in\ vivo}]*[GSH]\]

\[\frac{d[HMGS]}{dt} = k_1*[HCHO_{in\ vivo}][GSH] - k_2[HMGS][NAD^+][adh]\]

\[\frac{d[NADH]}{dt} = k_2*[HMGS][NAD^+][adh] - k_2*[HMGS][NAD^+][adh]\]

\[\frac{d[S-formylglutathione]}{dt} = k_2*[HMGS][NAD^+][adh] - k_3*[S-formylglutathione]*[fghA]\]

\[\frac{d[mid_2]}{dt} = k_3*[S-formylglutathione][fghA] - k_4[mid_2]*[fdh]\]

\[\frac{d[mid_3]}{dt} = k_4*[mid2][fdh] - k_5[mid_3]\]

\[\frac{d[CO2_{in\ vivo}]}{dt} = k_5*[mid_3] - k_6*[CO2_{in\ vivo}]\]

\[\frac{d[CO2_{air}]}{dt} = k_6*[CO2_{in\ vivo}]\]

\[\frac{d[GSH]}{dt} = -k_1*[HCHO_{in\ vivo}][GSH] + k_2[HMGS][NAD^+][adh]\]

\[\frac{d[fghA]}{dt} = -k_3*[S-formylglutathione][fghA] + k_3[S-formylglutathione]*[fghA]\]

\[\frac{d[fdh]}{dt} = -k_4*[mid_2][fdh] + k_4[mid_2]*[fdh]\]

\[\frac{d[adh]}{dt} = -k_2*[HMGS][NAD^+][adh] + k_2*[HMGS][NAD^+][adh]\]

\[\frac{d[NAD^+]}{dt} = -k_2*[HMGS][NAD^+][adh] + k_2*[HMGS][NAD^+][adh]\]

 

 

  1. Flux Balance Analysis (FBA)

 

The first part is about flux balance analysis, an approach for studying a genome-scale metabolic network, containing almost all known metabolic reactions, genes, and related enzymes for an organism (Orth, 2010). The reason we chose this approach as the simulation part is not only because of the constructed published microorganism models, but also through dynamic FBA and some well-developed toolboxes such as COBRA (Heirendt et al., 2019) and COMETs (Dukovski et al., 2021), we can track the changes in metabolite concentrations in the presence of complex reactions within microbial bodies, and even simulate the co-cultivation of blue algae and e. coli.

 

1.1  Verification

First of all, we need to verify the model. We searched for the existing E. coli Bl21 growth experiment data (Christensen & Eriksen, 2002) to compare with the chosen model iB21_1374 (Monk et al., 2013). We adjusted the intake limitation of glucose intake of the model to restrict it from unlimited nutrient intake, and the simulation result is presented below.

 

Table 1. Specific growth rates and yield coefficients for growth of E. coli BL21 on glucose as carbon substrate.

Parameters

Reference value

(Christensen & Eriksen, 2002)

Simulation value

()

0.380.03

0.3796

0.410.01

0.5109

9.60.8

11.22

(), specific growth rates estimated from gram in dry weight,
, yield of biomass per unit of glucose consumed.

 

The equation selected to fit the growth curve is

Figure 2 E.coli simulation growth curve

 

Figure 3 Fitted curve of E.coli’s dry weight against glucose uptake.

Figure 4 Fitted curve of E.coli’s dry weight against NH4+ uptake. The total volume is set to sufficiently high.

Notice that the fit line covers the blue line for the simulation result.

 

Although the result is not perfectly consistent with the reference value, the bias is acceptable and predictable, for the real case has many more hidden variables, such as temperature and thermodynamic process.

 

Meanwhile, we also selected a GEM of Synechococcus PCC 7942, iJB785, published by Broddrick et al. (2016).

 

1.2  Modification

Then, we modify the model according to the metabolite pathway, adding sucrose transportation reaction on the algae model, extracellular sucrose hydrolysis reaction, and formaldehyde degradation pathway on two separate E. coli models.

 

1.3  Dynamic-FBA

After that, we adjusted the parameters to fit the experiment data, and the simulation result of E. colis HCHO consumption is presented below.

 

Figure 5 Left, experiment data obtained in wet lab, LB is medium, Kan is an antibiotic, Cb, 2b are two mutants of BL21. Cb is the improved mutant. Right, simulation results from dynamic FBA on modified GEM. The biomass is set at 300 mg conversion from the OD600 value of Cb at 0.

 

The experiment result of BL21 and the simulation result were similar in trend, which validated the model.

 

1.4  Ratio of coculture system’s components

We were curious about how to set the ratio between algae and E. coli so that the system would function stably.

 

To do that, we first analyzed the maximum production rate of algae.

 

Figure 6 Dry weight growth rate against sucrose production rate in silico, through robust analysis between biomass function and sucrose secretion reaction on GEM of Synechococcus PCC 7942.

 

The maximum sucrose production rate in unit biomass is 0.15 mM/h without a decrease in growth rate. We take this as the chosen value, assuming the cell always tries to maximize its growth rate.

 

Then we set up an optimization model to find the optimal ratio, based on a few assumptions.

 

  1. They are both at exponential growth stage
  2. All sucrose is hydrolyzed as soon as secreted.
  3. Do not distinguish the E. coli for sucrose hydrolysis and formaldehyde degradation.
  4. Other hidden effects are negligible.

 

 

Parameters

Value

()

0.38

()

0.054

0.41

()

0.15

()

180.16

 

We used Python to solve the model, and the optimized result is

We were excited to find that the result is similar to a coculture’s growth data collected from research by Hays et al. (2017).

Figure 7 Photobioreactor continuous co-culture replicates of E. coli W ΔcscR /cscB+ S. elongatus were set up in constant light. The optical density of the entire culture (black scatter plots) and cell counts for the individual species were tracked (green S. elongatus, orange E. coli W ΔcscR) for more than 2 weeks (Hays et al., 2017).

 

The difference may come from multiple reasons, such as some complex interactions and environmental factors. However, this still provided a starting point for future coculture experiments and hardware setup.

 

  1. FEA

 

This part is aimed at giving an insight into how alginate beads affect the metabolite of the microorganism embedded. In our project, we designed three types of beads, one with blue algae for sucrose production, one with E. coli for sucrose degradation, and one with E. coli for formaldehyde degradation. During the development, we realized that we could increase the accuracy by the embedding ODE system, which describes the formaldehyde degradation pathway.

 

The model can be separated into two parts: diffusion and metabolite. We first developed the diffusion model and started a rough simulation, then combined the metabolite part with the diffusion model.

 

2.1  Metabolite function

Before we start, notice that one bead will only have one function.

 

2.1.1        Formaldehyde detection

The detection is based on the FrmR operon regulation on GFP, and the outcome is that when formaldehyde is present, the beads will give green fluorescence under UV light.

 

Based on the research by Woolston et al. (2018), we chose Hill’s equation to fit the experiment data, and the result is presented below.

Figure 8 Curve fitting with Hill’s equation.

Fitted equation:

Where S represents the fluorescence signal, and [F] represents the formaldehyde concentration. The result met our expectations.

Parameters

Estimate

Std. error

P

32.7494

2.5972

5.57e-05 ***

1.6467

0.1996

4.27e-04 ***

 

2.1.2        Photosynthesis

We described photosynthesis with the Michaelis-Menten equation.

 

Parameters

Value(Ritchie, 2008)

36.7

0.062 μmol/s

0.005 μmol/s

Assume
 

This means all glucose from photosynthesis is transformed into sucrose and secreted totally.

2.1.3        Sucrose hydrolysis

This could be summarized as a Michaelis-Menten equation.

 

Parameters

Value(Wanker et al., 1995)

64

0.91 nM/s

Notice that the values are measured from a similar exoenzyme β-fructosidase from the Bacillus subtilis due to the lack of information from Zymomonas mobilis where our SacC gene originated from.

 

2.2  Diffusion model

2.2.1        Model derivation

The diffusion model is based on Fick’s second law

In three-dimension space, the formula can be rewrote with Laplacian operator

Where

Which means for a point, the change of concentration within unit time equals the point’s divergence times original concentration times the diffusion coefficient of the substance in a specific matter, in this case, the alginate is the matter.

 

The diffusion coefficient is related to the type of substance and the concentration of the medium, the values are given below, all at 0.2% of alginate.

 

Diffusion coefficient

Value()

Reference

6.7

(McDonald, 1956)

5.2

(Cussler, 1984)

1.73

(Pu & Yang, 1988)

2.23

Manually, assume
 

10.0

Manually

 

This formula tells us how a substance diffuses into another quantitively, and more importantly, based on former work, we could set up a simulation with the finite element method (FEM), which is suitable for solving partial differential equations and visualization.

 

2.2.2        Light decay

According to Lambert-Beer Law, light will be absorbed when transmitting through a media. So we served an equation to describe the effect.

Where h stands for the distance light traveled in the media, and k is a coefficient.

 

Notice that light is emitted parallelly from the bottom.

 

2.2.3        Finite Element Method (FEM)

We will briefly introduce how it works, to learn more details, please refer to some professional resources.

 

FEM is a numerical method to approximate the solution of continuous, boundary-initial-value problems described by partial differential equations as discrete models (Tekkaya & Soyarslan, 2014), which is commonly applied in practice in the engineering field.

 

Why did we choose FEM? Sometimes obtaining an analysis solution for the model is unnecessary, and we can tolerate the bias, then, instead of solving the model directly, we split the space into thousands even millions of nodes to discretization the space, and update the state of nodes from initial state step by step to discretization the time to simplify the problem.

 

With this idea, we develop the following simulation algorithm.

Figure 9 The numeric simulation algorithm based on diffusion model and FEM.

 

With the simulation, we are able to get insight on how to build and improve the system and prevent some possible traps in advance.

 

2.3  Result

Our simulation will produce a 3-D vector field showing the flux with arrow and concentration with color, a 2-D contour graph showing the concentration distribution on the symmetry plane of the beads parallel to xOy plane in 3-D space.

 

The units are verified, and values are correctly transformed. We picked four representative types to simulate, beads with only diffusion of sucrose, with sucrose hydrolysis, with sucrose production, or with formaldehyde detection.

 

2.3.1        Simple beads

In this case, we do not care about the metabolite term and set it to zero.

Figure 10 Dynamic FEA result of sucrose’s diffusion from pure bead to outside. The environment totally depends on the bead’s diffusion.

 

When this result was obtained, it immediately came to our mind that, if the beads are static, the system's efficiency will be significantly affected due to the limit of substance’s diffusion near the beads.

 

To prove our thought, we refresh the outside concentration with 0 in each cycle to simulate the pump.

 

Figure 11 Dynamic FEA result of sucrose’s diffusion from pure bead to outside. The environment is refreshed in every cycle. The arrow’s length is zero outside the bead, for the constant refresh.

 

The result showed that at 300ms, the concentration at the center is lower than when not refreshed, and the reduction speed is faster, which supports our previous thought and motivates the installation of a pump in our hardware.

 

2.3.2        Consider sucrose hydrolysis

The metabolite term is based on the β-D-fructosidase Michaelis-Menten equation.

 

To show the impact on consumption, we set a small initial value inside the beads, and refresh the environment with a maximum value of 200 μM.

 

Figure 12 Final FEA result of the bead’s sucrose concentration distribution change against time, with an initial value of 2μM inside, and a refresh value of 200μM outside.

 

The result clearly shows that the diffusion will delay the response inside the bead, appearing as a short drop at the beginning, and the system will reach an equilibrium between diffusion and consumption. This confirmed the hypothesis that embedding will affect the system, for the cell density may be higher at the edge than at the center, and the system may not reach its upper limit due to the uneven distribution. The result also provides a direction for optimizing the system, by finding a balance between the size and the alginate concentration of the bead to control the side effect of diffusion.

 

2.3.3        Consider sucrose production

The metabolite term is based on the photosynthesis Michaelis-Menten equation.

 

Figure 13 Final FEA result of the bead’s sucrose concentration distribution change against time in production case. Notice that all nodes inside the beads will produce sucrose. k=0.5

 

The result showed that due to the effect from the diffusion of sucrose, the secreted sucrose accumulated inside the beads. To make the light decay more significantly, we may manually increase the value of k.

 

Figure 14 Final FEA result of the bead’s sucrose concentration distribution change against time in production case with strong light decay. k=5

 

This may occur when cell density is high, thus, we may need to optimize the size of the bead to reduce the side effect of light decay.

 

2.3.4        Consider formaldehyde detection

Now we focused on formaldehyde detection. We decided to simulate the response of the beads toward the decline of the environment formaldehyde concentration. The degradation rate is set as a constant

obtained from the wet lab data when the formaldehyde degradation rate by E. coli is stable. The final result is shown below.

 

Figure 15 Final FEA result of the bead’s response to environment formaldehyde decline.

 

The result showed that there was a slight delay in the center’s response. It is reasonable to assume the hydrogen sulfide detection performed similarly.

 

  1. Conclusion and Discussion

We first set up the general biochemical pathway diagram with the help of wet lab members, which clarifies how the system works. Then, we start research and try to implement FBA to verify the system’s function.

 

When our team decided to build hardware, we started looking for a method to simulate the system and chose to simulate the beads. After that, we developed a model describing an “alginate bead with microorganism embedded”, and set up a mathematical simulation with the finite element method. The result served some insights into how diffusion affects the system, including the necessity of a pump, and the potential effect on cell distribution and growth.

 

There are still many things to be done in the future. The cell growth is limited due to the volume of the beads, and the distribution inside the beads is uneven.

 

All models are wrong, but some are useful”, said George E. P. Box. After a long journey of exploring, meeting, learning, and debugging, we are glad that our work provided some meaningful results and understanding of the system’s function.

 

 

 

Reference

Broddrick, J. T., Rubin, B. E., Welkie, D. G., Du, N., Mih, N., Diamond, S., Lee, J. J., Golden, S. S., & Palsson, B. O. (2016). Unique attributes of cyanobacterial metabolism revealed by improved genome-scale metabolic modeling and essential gene analysis. Proc Natl Acad Sci U S A, 113(51), E8344-e8353. https://doi.org/10.1073/pnas.1613446113

Christensen, M. L., & Eriksen, N. T. (2002). Growth and proton exchange in recombinant Escherichia coli BL21. Enzyme and Microbial Technology, 31(4), 566-574. https://doi.org/https://doi.org/10.1016/S0141-0229(02)00153-9

Cussler, E. L. (1984). Diffusion: Mass Transfer in Fluid Systems.

Dukovski, I., Bajić, D., Chacón, J. M., Quintin, M., Vila, J. C. C., Sulheim, S., Pacheco, A. R., Bernstein, D. B., Riehl, W. J., Korolev, K. S., Sanchez, A., Harcombe, W. R., & Segrè, D. (2021). A metabolic modeling platform for the computation of microbial ecosystems in time and space (COMETS). Nature Protocols, 16(11), 5030-5082. https://doi.org/10.1038/s41596-021-00593-3

Hays, S. G., Yan, L. L. W., Silver, P. A., & Ducat, D. C. (2017). Synthetic photosynthetic consortia define interactions leading to robustness and photoproduction. Journal of Biological Engineering, 11(1). https://doi.org/10.1186/s13036-017-0048-5

Heirendt, L., Arreckx, S., Pfau, T., Mendoza, S. N., Richelle, A., Heinken, A., Haraldsdóttir, H. S., Wachowiak, J., Keating, S. M., Vlasov, V., Magnusdóttir, S., Ng, C. Y., Preciat, G., Žagare, A., Chan, S. H. J., Aurich, M. K., Clancy, C. M., Modamio, J., Sauls, J. T., . . . Fleming, R. M. T. (2019). Creation and analysis of biochemical constraint-based models using the COBRA Toolbox v.3.0. Nature Protocols, 14(3), 639-702. https://doi.org/10.1038/s41596-018-0098-2

McDonald, H. J. (1956). Electrochemistry in Biology and Medicine T. Shedlovsky (Ed.) New York, John Wiley & Sons, Inc., 1955, 369 pp. Clinical Chemistry, 2(1), 56-56. https://doi.org/10.1093/clinchem/2.1.56

Monk, J. M., Charusanti, P., Aziz, R. K., Lerman, J. A., Premyodhin, N., Orth, J. D., Feist, A. M., & Palsson, B. Ø. (2013). Genome-scale metabolic reconstructions of multiple <i>Escherichia coli</i> strains highlight strain-specific adaptations to nutritional environments. Proceedings of the National Academy of Sciences, 110(50), 20338-20343. https://doi.org/10.1073/pnas.1307797110

Novère, N. L., Hucka, M., Mi, H., Moodie, S., Schreiber, F., Sorokin, A., Demir, E., Wegner, K., Aladjem, M. I., Wimalaratne, S. M., Bergman, F. T., Gauges, R., Ghazal, P., Kawaji, H., Li, L., Matsuoka, Y., Villéger, A., Boyd, S. E., Calzone, L., . . . Kitano, H. (2009). The Systems Biology Graphical Notation. Nature Biotechnology, 27(8), 735-741. https://doi.org/10.1038/nbt.1558

Orth, J., Thiele, I. & Palsson, B. (2010). What is flux balance analysis? Nat Biotechnol, 28, 245–248.

Pu, H. T., & Yang, R. Y. (1988). Diffusion of sucrose and yohimbine in calcium alginate gel beads with or without entrapped plant cells. Biotechnol Bioeng, 32(7), 891-896. https://doi.org/10.1002/bit.260320707

Ritchie, R. J. (2008). Fitting light saturation curves were measured using modulated fluorometry. Photosynthesis Research, 96(3), 201-215. https://doi.org/10.1007/s11120-008-9300-7

Tekkaya, A. E., & Soyarslan, C. (2014). Finite Element Method. In (pp. 508-514). Springer Berlin Heidelberg. https://doi.org/10.1007/978-3-642-20617-7_16699

Wanker, E., Huber, A., & Schwab, H. (1995). Purification and characterization of the Bacillus subtilis levanase produced in Escherichia coli. Appl Environ Microbiol, 61(5), 1953-1958. https://doi.org/10.1128/aem.61.5.1953-1958.1995

Woolston, B. M., Roth, T., Kohale, I., Liu, D. R., & Stephanopoulos, G. (2018). Development of a formaldehyde biosensor with application to synthetic methylotrophy. Biotechnol Bioeng, 115(1), 206-215. https://doi.org/10.1002/bit.26455