The holistic approach to convert biology from an empirical into an axiomatic science is called systems biology. Advancements in this field have enabled researchers to create genome-scale metabolic models (GEMs) which gave invaluable insights into both the structure of metabolic networks as well as the effects of gene perturbations after metabolic engineering (Feist & Palsson 2008), (Lewis et al. 2012), (Zhang & Hua 2016). The approach of constrained-based reconstruction and analysis (COBRA) makes it possible to examine metabolic models based on pure assumptions and in silico simulations (Lewis et al. 2012), (Kabimoldayev et al. 2018).

Flux Balance Analysis

The most abundant method to achieve this is Flux Balance Analysis (FBA). Based on the assumption that eons of evolution guided organisms towards optimality in their network and growth, FBA assumes that an optimal steady state in metabolite concentration is given at any time in a cells growth (Raman & Chandra 2009). While the concentrations themselves cannot be deduced, it is however possible to get a flux distribution v after solving the optimization problem S × v = 0 (Goryanin 18.10.2020), (Orth et al. 2010), (Rajvanshi & Kareenhalli 2013). S is the stoichiometric matrix consisting of the stoichiometric coefficients of the reactions – therefore, only knowledge of the genomics and biochemistry of the target organism is needed for simulation (Reed 2017). However more constraints need to be imposed, including the reaction directionality based on thermodynamics, maximal transport capacities, but most importantly an objective function (Kauffman et al. 2003), (Orth et al. 2010), (Raman et al. 2009), (Rajvanshi & Kareenhalli 2013). This objective function is the mathematical expression of the main assumption of FBA and makes it possible to obtain a solution, which is a physiologically senseful flux distribution (Kauffman et al. 2003), (Lewis et al. 2012). While basic FBA does not take into account flux dynamics, metabolite concentrations or regulatory mechanisms, its purposes are somewhat limited, especially quantitatively (Zhang & Hua 2016). Still, its simplicity makes it a powerful tool for both foundational science and synthetic biology alike and modified versions like Flux Variability Analysis (FVA) (Segre et al. 2002), parsimonious enzyme usage FBA (pFBA) (Lewis et al. 2010) and geometric FBA (gFBA) (Smallbone & Simeonidis 2009) enhance the base approach while maintaining the simple nature of not requiring experimental, kinetic data.

Our approach

To our knowledge, no metabolic model of Physcomitrella patens exists as of yet that is not automatically generated. Those models are too large to manually curate them, especially on a limited home computer, and therefore we decided to approximate the metabolism of such a plant with the majority of the reactions from a metabolic minimal model of Chlamydomonas rheinhardtii (Kliphuis et al. 2012) and the biomass objective function (BOF) from an Arabidopsis thaliana model (Yuan et al. 2016).


We simplified the model of C. rheinhardtii in some ways like foregoing the photorespiration and the lipid metabolism, while including biosynthesis reactions for some cofactors. The model contains four compartments: extracellular, cellular, stroma (chloroplast) and mitochondrion (matrix).

The troubleshooting of the model did take a considerable amount of time so that we do not have gene annotations for any of the reactions as of now. The main reaction relevant for our approach is a reaction for synthesizing polyphosphate, the sum formula of which we approximated with PO3, and a corresponding sink reaction. Moreover, sink reactions for stomatal protons and mitochondrial oxygen were needed to make the model viable – we do not know why those metabolites over accumulate but we believe that those minor metabolites are not that relevant especially for the phosphate flux.

A Visualization of the model can be seen in Figure 1.

Fig. 1: Visualization of the main parts of thje metabolic network of our Physcomitrella patens model. The major subsystems are shown in their respective colors in the legend. The flux distribution is shown as well and was obtained by applying FBA on the model with the biomass reaction as the objective function.

Analysis and Results

Default FBA was applied on the model.

First, the model was optimized towards polyphosphate production to get the maximal flux through the polyphosphate synthesis reaction. The resulting flux was ~ 28.5. The model was then optimized towards biomass production and the flux through the polyphosphate synthesis set to increase in 100 intervals from a lower bound of 0 to the maximally possible. Both the biomass and phosphate import flux were plotted against the polyphosphate synthesis flux and the results can be seen in Figure 2.

Fig. 2: Depicted is the flux of the biomass reaction (in blue) and the phosphate import reaction (in red) on the y-axis against the polyphosphate reaction flux on the x-axis. The flux unit is mmol ✕ g-1 Dry Weight ✕ h-1.

What can first be seen is that logically the growth decreases the more the system focuses on polyphosphate synthesis.

Moreover, although the system does not incorporate any regulatory functions, the one effect of the phosphate starvation response shows itself in the increased phosphate import flux: that being the overexpression of phosphate importers (Plaxton & Tran 2011), (Wang et al. 2008). This makes sense because the synthesis of polyphosphate drains phosphate from the rest of the metabolism, and thus, causes internal phosphate deprivation (Wei 2020).

Furthermore, we also tried to consider the possibility of magnesium starvation due to complexation with cytosolic polyphosphate (Rudat et al. 2018). For this we added another sink reaction for polyphosphate, where two polyphosphate “molecules” vanish together with one magnesium atom (as the main polyphosphate unit has a charge of -1 and magnesium is a divalent cation). Then we repeated the steps from before, but this time we first optimized for this new sink reaction und used its maximal flux of ~ 14.25 as the upper bound for the optimization for growth. Also, the flux of the magnesium import reaction was plotted. The results can be seen in Figure 3 - as a comparison, the results of the exact same procedure as for Fig. 2 with the added magnesium flux are given in Figure 4.

Fig. 3: Depicted is the flux of the biomass reaction (in blue), the phosphate import reaction (in red) and the magnesium import reaction (in cyan) on the y-axis against the polyphosphate reaction flux on the x-axis. The system was optimized towards biomass production with a forced flux trough the polyphosphate sink reaction with magnesium starvation. The flux unit is mmol ✕ g-1 Dry Weight ✕ h-1.
Fig. 4: Depicted is the flux of the biomass reaction (in blue), the phosphate import reaction (in red) and the magnesium import reaction (in cyan) on the y-axis against the polyphosphate reaction flux on the x-axis. The system was optimized towards biomass production with a forced flux trough the polyphosphate synthesis reaction. The flux unit is mmol ✕ g-1 Dry Weight ✕ h-1.

The trend is the same with the one observable in Fig. 2 but with the difference, that the overall flux through the polyphosphate synthesis reaction is lower. This means that magnesium starvation can be coped with the same way as phosphate starvation, i. e. an overexpression of magnesium transporters, but only to a lower degree. As FBA does not incorporate internal metabolite concentration, a maximum magnesium capacity of polyphosphate cannot be simulated (which should be observed in vivo), as the results are just relative, not absolute. It would therefore be interesting to see whether this can be implemented and what the exact implications of this would be.


In the beginning, we had two assumptions about our model regarding the project: 1) the growth will be greatly inhibited due to the phosphate flux redirection towards polyphosphate and 2) the increasing phosphate import increases the allowable polyphosphate flux. While 2) has been confirmed because of the reason stated above, 1) has been disproven. It seems that the enhanced phosphate import balances the withdrawal of phosphate from the main metabolism. What can be deduced from the results as well, is that phosphate can not be significantly taken from the core metabolism. This is likely due to the optimality assumption, as the phosphate demand in the metabolism may be already at its minimum. This is confirmed when applying pFBA to the model. pFBA finds the lowest overall flux through the system for the same flux value of the objective (Lewis et al. 2010). pFBA delivers the same results as FBA.

What can be done as a next step is implementing the lipid remodeling as a response to phosphate starvation (Abida et al. 2015), (Wang et al. 2008). This can be done by extending the lipid metabolism and varying the composition of the LIPIDS metabolite (representing the lipidome in our model). Perhaps even simple regulatory mechanisms can be included, like a more intricate and direct phosphate starvation response.

Our model is a foundation for future analyses like these and others in the organism P. patens.

Abida, H., Dolch, L. J., Meï, C., Villanova, V., Conte, M., Block, M. A., ... & Rébeillé, F. (2015). Membrane glycerolipid remodeling triggered by nitrogen and phosphorus starvation in Phaeodactylum tricornutum. Plant physiology, 167(1), 118-136.

Feist, A. M., & Palsson, B. Ø. (2008). The growing scope of applications of genome-scale metabolic reconstructions using Escherichia coli. Nature biotechnology, 26(6), 659-667.

Goryanin, I. Flux Balance Analysis. Retrieved from Last date of access: 18.10.2020.

Kabimoldayev, I., Nguyen, A. D., Yang, L., Park, S., Lee, E. Y., & Kim, D. (2018). Basics of genome-scale metabolic modeling and applications on C1-utilization. FEMS microbiology letters, 365(20), fny241.

Kauffman, K. J., Prakash, P., & Edwards, J. S. (2003). Advances in flux balance analysis. Current opinion in biotechnology, 14(5), 491-496.

Kliphuis, A. M., Klok, A. J., Martens, D. E., Lamers, P. P., Janssen, M., & Wijffels, R. H. (2012). Metabolic modeling of Chlamydomonas reinhardtii: energy requirements for photoautotrophic growth and maintenance. Journal of applied phycology, 24(2), 253-266.

Lewis, N. E., Hixson, K. K., Conrad, T. M., Lerman, J. A., Charusanti, P., Polpitiya, A. D., ... & Weitz, K. K. (2010). Omic data from evolved E. coli are consistent with computed optimal growth from genome‐scale models. Molecular systems biology, 6(1), 390.

Lewis, N. E., Nagarajan, H., & Palsson, B. O. (2012). Constraining the metabolic genotype–phenotype relationship using a phylogeny of in silico methods. Nature Reviews Microbiology, 10(4), 291-305.

Orth, J. D., Thiele, I., & Palsson, B. Ø. (2010). What is flux balance analysis?. Nature biotechnology, 28(3), 245-248

Plaxton, W. C., & Tran, H. T. (2011). Metabolic adaptations of phosphate-starved plants. Plant physiology, 156(3), 1006-1015.

Rajvanshi, M, Kareenhalli, V. V. (2013). Flux Balance Analysis. In Dubitzky, W., Wolkenhauer, O., Yokota, H., & Cho, K. H. (Eds.) Encyclopedia of systems biology (pp. 749 – 752). Springer Publishing Company, Incorporated.

Raman, K., & Chandra, N. (2009). Flux balance analysis of biological systems: applications and challenges. Briefings in bioinformatics, 10(4), 435-449.

Reed, J. L. (2017, July). Genome-scale metabolic modeling and its application to microbial communities. In The Chemistry of Microbiomes: Proceedings of a Seminar Series. National Academies Press.

Rudat, A. K., Pokhrel, A., Green, T. J., & Gray, M. J. (2018). Mutations in Escherichia coli polyphosphate kinase that lead to dramatically increased in vivo polyphosphate levels. Journal of bacteriology, 200(6).

Segre, D., Vitkup, D., & Church, G. M. (2002). Analysis of optimality in natural and perturbed metabolic networks. Proceedings of the National Academy of Sciences, 99(23), 15112-15117.

Smallbone, K., & Simeonidis, E. (2009). Flux balance analysis: a geometric perspective. Journal of theoretical biology, 258(2), 311-315.

Wang, Y., Secco, D., & Poirier, Y. (2008). Characterization of the PHO1 gene family and the responses to phosphate deficiency of Physcomitrella patens. Plant Physiology, 146(2), 646-656.

Wei, R., Wang, X., Zhang, W., Shen, J., Zhang, H., Gao, Y., & Yang, L. (2020). The improved phosphorus utilization and reduced phosphorus consumption of ppk-expressing transgenic rice. Field Crops Research, 248, 107715.

Yuan, H., Cheung, C. Y., Hilbers, P. A., & van Riel, N. A. (2016). Flux balance analysis of plant metabolism: the effect of biomass composition and model structure on model predictions. Frontiers in plant science, 7, 537.

Zhang, C., & Hua, Q. (2016). Applications of genome-scale metabolic models in biotechnology and systems medicine. Frontiers in physiology, 6, 413.