Implementation

Our digital twin model was implemented using the MATLAB COBRA toolbox for genome-scale metabolic modeling. The toolbox was accessed using the COBRA package for python (Ebrahim et al., 2013).

Reaction List Implementation

The initial metabolic reaction list ‘ICW1057’ for P. fluorescens SBW25 was obtained from the only genome-scale metabolic model for P. fluorescens (Haung and Lin, 2019). ICW1057 reaction list was used as the foundation for helper and main strain, however, this model showed poor results as key reactions like glucose uptake were missing. By stimulating growth on minimal media with known P. fluorescens, further missing reactions and metabolites were identified (https://bacdive.dsmz.de/strain/12851). 57 reactions and 30 metabolites were added to the ICW1057 model to obtain the growth behavior for the ReMixHD strain P. fluorescens DSM50090 (see Table 2.1 and 2.2). The improved reaction list for the P. fluorescens DSM50090 model is available here.

Metabolites added to the P. fluorescens DSM 50090 model
ID Name Compartments Strain
0 cpd10001_c0 Terephtalic_Acid_c0 c0 Main
1 cpd10002_c0 Poly_Ethylene_c0 c0 Helper
2 cpd10003_c0 Ethylene_Glycol_c0 c0 Helper
3 cpd10004_c0 DCD_c0 c0 Main
4 cpd10005_c0 Octanol_c0 c0 Main, Helper
5 cpd10006_c0 Hexanol_c0 c0 Main, Helper
6 cpd10007_c0 Decanol_c0 c0 Main, Helper
7 cpd10008_c0 Heptanol_c0 c0 Main, Helper
8 cpd10009_c0 Pentanol_c0 c0 Main, Helper
9 cpd10010_c0 Nonanol_c0 c0 Main, Helper
10 cpd10011_c0 Undecanol_c0 c0 Main, Helper
11 cpd10012_c0 Dodecanol_c0 c0 Main, Helper
12 cpd10013_c0 Dodecanoyl_CoA_c0 c0 Main, Helper
13 cpd10014_c0 Undecanoyl_CoA_c0 c0 Main, Helper
14 cpd10015_c0 Decanoyl_CoA_c0 c0 Main, Helper
15 cpd10016_c0 Nonanoyl_CoA_c0 c0 Main, Helper
16 cpd10017_c0 Octanoyl_CoA_c0 c0 Main, Helper
17 cpd10018_c0 Heptanoyl_CoA_c0 c0 Main, Helper
18 cpd10019_c0 Hexanoyl_CoA_c0 c0 Main, Helper
19 cpd10020_c0 Pentanoyl_CoA_c0 c0 Main, Helper
20 cpd10021_c0 Butanoyl_CoA_c0 c0 Main, Helper
21 cpd00029_e0 Acetate_e0 e0 Main, Helper
22 cpd10003_e0 Ethylene_Glycol_e0 e0 Main, Helper
23 cpd10002_e0 Poly_Ethylene_e0 e0 Helper
24 cpd00224_e0 L_Arabinose_e0 e0 Main, Helper
25 cpd00036_e0 Succinate_e0 e0 Main, Helper
26 cpd00161_e0 Threonine_e0 e0 Main, Helper
27 cpd00363_e0 Ethanol_e0 e0 Main, Helper
28 cpd10001_e0 Terephtalic_Acid_e0 e0 Main
29 cpd10022_e0 PET_e0 e0 Main
Reactions added to the P. fluorescens DSM 50090 model
ID Name Compartments Strain
0 EX_cpd00013_1_e0 EX_NH3_1_e0 -> e0 Main, Helper
1 IM_cpd00027 IM_Glucose c0 -> e0 Main, Helper
2 IM_cpd00033 IM_Glycine c0 -> e0 Main, Helper
3 IM_cpd00078 IM_Tryptophane c0 -> e0 Main, Helper
4 IM_cpd00100 IM_Glycerol c0 -> e0 Main, Helper
5 EX_cpd00029_e0 EX_Acetate_e0 -> e0 Main, Helper
6 IM_cpd00029 IM_Acetate c0 -> e0 Main, Helper
7 EX_cpd10003_e0 EX_Ethylene_Glycol_e0 -> e0 Main, Helper
8 IM_cpd10003 IM_Ethylene_Glycol c0 -> e0 Main, Helper
9 EX_cpd10002_e0 EX_Poly_Ethylene_e0 -> e0 Helper
10 IM_cpd10002 IM_Poly_Ethylene c0 -> e0 Helper
11 EX_cpd00224_e0 EX_L_Arabinose_e0 -> e0 Main, Helper
12 IM_cpd00224 IM_L_Arabinose c0 -> e0 Main, Helper
13 EX_cpd00036_e0 EX_Succinate_e0 -> e0 Main, Helper
14 IM_cpd00036 IM_Succinate c0 -> e0 Main, Helper
15 EX_cpd00161_e0 EX_Threonine_e0 -> e0 Main, Helper
16 IM_cpd00161 IM_Threonine c0 -> e0 Main, Helper
17 EX_cpd00363_e0 EX_Ethanol_e0 -> e0 Main, Helper
18 IM_cpd00363 IM_Ethanol c0 -> e0 Main, Helper
19 EX_cpd10001_e0 EX_Terephtalic_Acid_e0 -> e0 Main
20 IM_cpd10001 IM_Terephtalic_Acid c0 -> e0 Main
21 EX_PET_e0 -> PET_e0 -> e0 Main
22 rxn20047_c0 PET_e0 -> Terephtalic acid + Ethylene_glycol -> e0 Main
23 rxn20001_c0 Terephtalic_Acid_c0 -> DCD_c0 -> c0 Main
24 rxn20002_c0 DCD_c0 -> Protocatechuate_c0 -> c0 Main, Helper
25 rxn20003_c0 4_Carboxymuconolactone_c0 -> 3_oxoadipate_enol_lactone_c0 -> c0 Main, Helper
26 rxn20004_c0 Ethylene_Glycol_c0 -> Glycolaldehyde_c0 -> c0 Main, Helper
27 rxn20005_c0 Glycolaldehyde_c0 -> Glycolate_c0 -> c0 Main, Helper
28 rxn20006_c0 Poly_Ethylene -> relative splitting -> c0 Main, Helper
29 rxn20007_c0 Hexanol_c0 -> Hexanoate_c0 -> c0 Main, Helper
30 rxn20008_c0 Octanol_c0 -> octanoate_c0 -> c0 Main, Helper
31 rxn20009_c0 Decanol_c0 -> Decanoate_c0 -> c0 Main, Helper
32 rxn20010_c0 Hexanoate_c0 -> Hexanoyl_CoA_c0 -> c0 Main, Helper
33 rxn20011_c0 octanoate_c0 -> Octanoyl_CoA_c0 -> c0 Main, Helper
34 rxn20012_c0 Decanoate_c0 -> Decanoyl_CoA_c0 -> c0 Main, Helper
35 rxn20013_c0 Pentanol_c0 -> Pentanoyl_CoA_c0 -> c0 Main, Helper
36 rxn20014_c0 Heptanol_c0 -> Heptanoyl_CoA_c0 -> c0 Main, Helper
37 rxn20015_c0 Nonanol_c0 -> Nonanoyl_CoA_c0 -> c0 Main, Helper
38 rxn20016_c0 Undecanol_c0 -> Undecanoyl_CoA_c0 -> c0 Main, Helper
39 rxn20017_c0 Undecanol_c0 -> Undecanoyl_CoA_c0 -> c0 Main, Helper
40 rxn20018_c0 Dodecanoyl_CoA_c0 -> Decanoyl_CoA_c0 -> c0 Main, Helper
41 rxn20019_c0 Decanoyl_CoA_c0 -> Octanoyl_CoA_c0 -> c0 Main, Helper
42 rxn20020_c0 Octanoyl_CoA_c0 -> Hexanoyl_CoA_c0 -> c0 Main, Helper
43 rxn20021_c0 Hexanoyl_CoA_c0 -> Butanoyl_CoA_c0 -> c0 Main, Helper
44 rxn20022_c0 Butanoyl_CoA_c0 -> Acetyl_CoA_c0 -> c0 Main, Helper
45 rxn20023_c0 Undecanoyl_CoA_c0 -> Nonanoyl_CoA_c0 -> c0 Main, Helper
46 rxn20024_c0 Nonanoyl_CoA_c0 -> Heptanoyl_CoA_c0 -> c0 Main, Helper
47 rxn20025_c0 Heptanoyl_CoA_c0 -> Pentanoyl_CoA_c0 -> c0 Main, Helper
48 rxn20026_c0 Pentanoyl_CoA_c0 -> Propionyl_CoA_c0 -> c0 Main, Helper
49 rxn20028_c0 Ribulose -> Glyceraldehyd_3_phosphat -> c0 Main, Helper
50 rxn20029_c0 PQQH2 -> PQQ -> c0 Main, Helper
51 rxn20031_c0 Maltose -> 2Glucose -> c0 Main, Helper
52 rxn20032_c0 kynurenine -> anthranilat -> c0 Main, Helper
53 rxn20033_c0 Anthralinat _> Catechol -> c0 Main, Helper
54 EX_Glucose_e0 -> EX_Glucose_e0 -> e0 Main, Helper
55 rxn20039_c0 Butanoyl_CoA -> Succinat -> c0 Main, Helper
56 rxn20045_c0 Poly_ethylene -> Butanoyl_CoA -> c0 Main, Helper

By adding the terephthalic metabolism, as introduced in the wetlab, to the preexisting model and linking it to the wild-type beta-ketoadipate pathway, a new model was created simulating our main strain. Following this pathway, terephthalic acid is introduced to the TCA cycle. The helper strain was established using the same base model and the PET and PE depolymerization pathways were added, as well as the AlkB mediated oxidation into n-alkanols. The introduced alkanols are naturally linked to the fatty acid metabolism of the helper strain.

The completed helper and main strain models were optimized for maximum biomass flux and cleaned of infeasible fluxes by implementing loopless FBA. The non-growth associated ATP-maintenance reaction was set to a minimal reaction rate of 46.3 mmol/gDW/h to maintain the physiologically necessary ATP usage to keep the cell alive (Chen et al., 2011). It was assumed that the flux observed in E. coli during aerobic growth is the same for P. fluorescens due to lacking data.

Model Performance Evaluation

The comparison of measured and predicted growth rates on minimal media with one carbon was chosen as a reliable and easy method to assess models performance. The carbon sources tested are metabolites from the citric acid cycle and central carbon metabolism.

Measuring growth data

The maximum growth rate of P. fluorescens was determined in 3 biological replicates as a baseline. 18 different carbon sources with a concentration of 1% (w/v) were tested on M9 minimal media to cover a broad spectrum of pathways and to eliminate potential diauxic growth effects.

Once the pre-cultures in LB medium reached an OD600 of 0.1, the cells were washed with M9 minimal salts without carbon sources removing excess LB medium. The washed preculture was inoculated 1:10 in 200µl and plated on 96-well plates with M9 minimal medium. OD600 was measured every 10 min for 24 h using a TECAN Spark spectrometer and the maximum growth rate was determined using the Growthcurver package in R 4.1.3 (Sprouffske & Wagner, 2016).

In silico growth predictions

The in silico growth rates were produced by setting all fluxes of exchange reactions with metabolites not contained in the medium to zero and optimized for biomass. A direct comparison of the predicted growth rates rpredicted and the measured rates rmeasured was not possible, as the two values are proportional but not identical. By implementing a proportionality constant termed scaling factor k to account for this difference, the two datasets can be compared:

The equation was rearranged and k was adjusted so the rescaled mean sum of squares converged to zero. This step was repeated twice with the upper and lower values of the 1 sigma interval of the measured growth rates to estimate the 1σ interval of the scaling factor.

The carbon sources were also chosen to avoid the overfitting of particular degradation pathways.

Established FBA workflow

All following simulations were obtained by following this workflow: First, all exchange fluxes of components not present in the simulated medium were set to zero. Flux constraining and biomass optimization was performed. Using the scaling factor k, the whole model was scaled down to yield physiologically correct results. The scaled predicted fluxes were used for yield and degradation prediction.

Digital Twin Creation

cFBA implementation (co-culture simulation)

The co-culture was simulated using the two metabolic models for the helper and the main strain. A shared medium was created including the limited amounts of nutrients to be observed. By subtracting the consumed nutrient values of both strains from the initial medium values, biomass growth was linked to depletion of nutrients. The growth of the biomass is terminated when all energy sources in the medium are consumed. Predictions with cFBA followed the same workflow described for FBA, however, the objective function was altered to optimize for the sum of the biomass fluxes of helper and main strain.

dFBA implementation (time resolved simulation)

For dFBA simulations, an initial medium composition of 200 mmol carbon source (20 mmol for polymer degradation) was set. Other necessary non-carbon metabolites (e.g oxygen, carbondioxide, water, etc.) were set to 10 mol to ensure that the growth simulation is not limited by them. The fluxes are calculated and multiplied by a small, constant time interval to give the approximate change in media composition over the time frame. A time step of 0.1 h was found to give good resolution with moderate calculation time. Subtracting this change from the initial medium composition gives a new initial medium composition. The process is repeated until the substrate concentration in the medium reaches negative values. At each timestep, a vector of media concentration is stored for plotting and analysis.

Digital twin implementation - coculture dynamic flux-balance analysis(cdFBA)

dFBA of the coculture model was performed similarly to the dFBA of a singular model. However, the optimization step at each time point included two optimizations on the same medium composition. One optimization calculates the biomass and the substrate effects for the main strain, the other for the helper strain. The media compositions were updated accordingly and stored for later analysis.

Discovery of potential growth control genes for the helper strain

To identify the metabolic rations with the highest impact on the helper strains growth, its biomass flux was simulated on PET/PE minimal medium and in an unlimited medium containing all possible extracellular metabolites. Then each reaction was knocked-out individually and the effect on the PET/PE and full medium growth recorded. Knockouts that resulted in less than 1% of initial growth rate on both media were deemed essential and present good candidates for controlling the helper strain’s growth.

Simulating industrial scale-up

To assess the future potential applications of the ReMixHD system, we developed a function to compute up-scaling processes. This method is capable of simulating a chemostat (also known as a continuous bioreactor), which maintains constant media composition and cell density, allowing for consistent cell growth at maximum growth rate.

Our up-scaling model utilizes the bioreactor's volume, elapsed time, and substrate properties to predict the degradation of various substrates, the production of biomass, emitted CO2, and polyhydroxyalkanoates (PHA). This extrapolation function is based on the following assumptions: A cell of P. fluorescens has a comparable dry weight to E. coli, estimated at 3 x 10-13 grams (Neidhardt et al., 1996). In the bioreactor, a dominant cell density of OD = 1.2 is observed. Assuming that 0.1 OD units correspond to 7 x 105 cells/ml, this results in a cell count of 8.4 x 109 cells/l.

Next, the total weight of dry biomass in the bioreactor can be calculated. From there, the determination of the biomass produced during a degradation process depends on the growth rate and duration. The generated PHA mass is derived from the potential storage capacity of PHA in P. fluorescens, which is reported to be 15% of the cell dry weight according to (Vladu et al., 2019).

Subsequently, the quantity of the degraded substrate was calculated. This was achieved by multiplying the substrate's uptake rate with CDWbioreactor, t, and the molecular weight of our substrate.

To compute the emitted CO2, the model uses the molecular weight of CO2 as a multiplier for its flux. This value is multiplied by CDWbioreactor and t.

This model serves as an extrapolation for industrial-scale predictions of degradation rates, product yields, and CO2 emissions. It is extensible to include the status of co-cultures, such as the proportion of helper and main strain.

Resources

Chen, X., Alonso, A. P., Allen, D. K., Reed, J. L., & Shachar-Hill, Y. (2011). Synergy between 13C-metabolic flux analysis and flux balance analysis for understanding metabolic adaption to anaerobiosis in E. coli. Metabolic Engineering, 13(1), 38–48. https://doi.org/10.1016/j.ymben.2010.11.004

Ebrahim, A., Lerman, J. A., Palsson, B. O., & Hyduke, D. R. (2013). COBRApy: COnstraints-Based Reconstruction and Analysis for Python. BMC systems biology, 7, 74. https://doi.org/10.1186/1752-0509-7-74

Huang, X. and Lin, Y.-H. (2020). Reconstruction and analysis of a three-compartment genome-scale metabolic model for Pseudomonas fluorescens. Biotechnology and Applied Biochemistry, 67, 133-139. https://doi.org/10.1002/bab.1852

Sprouffske, K., Wagner, A. (2016). Growthcurver: an R package for obtaining interpretable metrics from microbial growth curves. BMC Bioinformatics 17, 172. https://doi.org/10.1186/s12859-016-1016-7

Neidhardt et al. (1996). Escherichia coli and Salmonella typhimurium, 2nd edition..

Vladu, M.-G., Petrescu, M., Stefaniu, A., Savoiu, G., Spiridon, M., Eremia, M., & Irina, L. (2019). Studies on polyhydroxyalkanoates biosynthesis by some Pseudomonas spp. strains. Romanian Biotechnological Letters, 24, 388-394. https://doi.org/10.25083/rbl/24.3/388.394