Open Access

Simulating cyanobacterial phenotypes by integrating flux balance analysis, kinetics, and a light distribution function

  • Lian He1,
  • Stephen G. Wu1,
  • Ni Wan2,
  • Adrienne C. Reding3 and
  • Yinjie J. Tang1Email author
Microbial Cell Factories201514:206

Received: 7 October 2015

Accepted: 12 December 2015

Published: 24 December 2015



Genome-scale models (GSMs) are widely used to predict cyanobacterial phenotypes in photobioreactors (PBRs). However, stoichiometric GSMs mainly focus on fluxome that result in maximal yields. Cyanobacterial metabolism is controlled by both intracellular enzymes and photobioreactor conditions. To connect both intracellular and extracellular information and achieve a better understanding of PBRs productivities, this study integrates a genome-scale metabolic model of Synechocystis 6803 with growth kinetics, cell movements, and a light distribution function. The hybrid platform not only maps flux dynamics in cells of sub-populations but also predicts overall production titer and rate in PBRs.


Analysis of the integrated GSM demonstrates several results. First, cyanobacteria are capable of reaching high biomass concentration (>20 g/L in 21 days) in PBRs without light and CO2 mass transfer limitations. Second, fluxome in a single cyanobacterium may show stochastic changes due to random cell movements in PBRs. Third, insufficient light due to cell self-shading can activate the oxidative pentose phosphate pathway in subpopulation cells. Fourth, the model indicates that the removal of glycogen synthesis pathway may not improve cyanobacterial bio-production in large-size PBRs, because glycogen can support cell growth in the dark zones. Based on experimental data, the integrated GSM estimates that Synechocystis 6803 in shake flask conditions has a photosynthesis efficiency of ~2.7 %.


The multiple-scale integrated GSM, which examines both intracellular and extracellular domains, can be used to predict production yield/rate/titer in large-size PBRs. More importantly, genetic engineering strategies predicted by a traditional GSM may work well only in optimal growth conditions. In contrast, the integrated GSM may reveal mutant physiologies in diverse bioreactor conditions, leading to the design of robust strains with high chances of success in industrial settings.


Glycogen Multiple-scale modeling Photobioreactors Photosynthesis efficiency Self-shading Synechocystis 6803


In photobioreactors (PBRs), light penetration depth at high cell density can be as short as a few centimeters [1]. Thus, during large-size PBR cultivation, cyanobacteria move continuously between the “light zone” (where light is sufficient) and the “dark zone” (where light is substantially shaded). As a consequence, cyanobacterial metabolism in PBRs is spatially and temporally dependent: cells have autotrophic growth in the light zone, and they perform heterotrophic growth in the dark zone by consuming energy-storage compounds. Moreover, PBR performances are also affected by the efficiency of CO2 gas–liquid transfer. To enhance mass transfer, people often use CO2-enriched air in combination with high intensity mixing. Many models have been developed to understand how cyanobacterial physiological dynamics are influenced by the light intensity, CO2 supply, temperature, and geometry of PBRs [26]. Those kinetic and reactor studies are useful in optimizing PBR design and operations. However, bioprocess modeling is unable to provide an understanding of intracellular enzyme functions and metabolic fluxes in cyanobacteria. To improve engineered microalgae strains’ metabolisms in large-size PBRs, it is necessary to link process models to metabolic models.

On the other hand, metabolic flux analyses (MFA) can quantify in vivo enzyme reaction rates, and thus allow us to investigate the flux phenotypes resulting from complicated gene-protein-metabolite regulations. 13C-MFA measures carbon fluxes through the central metabolism via 13C labeling experiments. Alternatively, genome-scale flux balance analysis (FBA) can generate a holistic intracellular flux distribution map [7] owing to its extended coverage of genomic information [8]. Computational platforms, such as COBRA [9] and OptForce [10], can predict genetic targets and guide rational designs of engineered strains. FBA can also be integrated with constraint-based elementary flux mode analysis to identify optimal pathways for bio-productions [11]. However, an inherent limitation of traditional GSM is that it predicts only flux distributions that result in maximal yields in an optimal culture condition. They cannot forecast mutant strains’ production titers and rates in dynamic and heterogeneous bioreactors.

In this study, the major goal is to demonstrate multiple-scale modeling approaches by linking cell metabolisms to PBR environmental fluctuations. Specifically, the modeling efforts focus on Synechocystis 6803, a most widely used cyanobacterial biorefinery. Appealing traits of this species include amenability to genetic modifications, well-studied genomics, and native genes for biosynthesis of alkanes/alkenes and hydrogen [1214]. To predict cyanobacterial growth and metabolic flux phenotypes in PBR settings, we integrated a genome-scale cyanobacteria model, iJN678 [15], with growth kinetics, cell movements based on reported PBR hydrodynamics, and a heterogeneous light distribution (Fig. 1). The model assumption is that heterogeneous PBR conditions affect cyanobacteria, leading to heterogeneous cell metabolisms in different sub-populations. Such an approach can provide biological information ranging from the intracellular domain to the PBR domain, and fill the gaps between systems biology and the PBR process. The multiple-scale modeling is useful for estimating mutant strains’ potentials to achieve the production metrics required for commercialization.
Fig. 1

Algorithm for simulating cyanobacterial growth and intracellular flux distribution by integrating the flux analysis model, kinetics, and a light distribution function


Simulation of cyanobacterial optimal growth in a cylindrical PBR

The integrated GSM was first applied to predict cyanobacterial growth in a cylindrical PBR, which was assumed to have a radius of 60 mm and a constant surface light intensity of 50 µE/m2/s. Although the maximal photosynthetic efficiency in photosynthetic species can reach 4.6–6 % [16], not all incident radiation in PBRs can be efficiently used by cyanobacteria, thereby resulting in a lower conversion efficiency [17]. Hence, we chose a photosynthesis efficiency of 1.5 %, which was within a reasonable range of actual photosynthesis efficiencies of microalgae [18]. Based on a previous study, the mass transfer rate of CO2 was assumed to be 10 h−1 [19]. Under such a condition, cyanobacterial biomass concentration could increase from 0.1 to 5 g/L in 3 weeks, provided that other mineral nutrients are supplied continuously (Fig. 2a and Additional file 2: Fig. S1). The modelling results also showed continuous decreases in the growth rate (Fig. 2b) and intracellular fluxes in the central metabolism (Fig. 2c–e), which was caused by a continuous decrease in local light intensity over time (Fig. 2g). As the ‘dark zone’ expanded, some cyanobacteria switched from autotrophic growth to heterotrophic growth in the late growth phase, and eventually became resting cells (Fig. 2f). The expanding ‘dark zone’ also led to a gradual reduction in glycogen content per gram of biomass, which was the same when all the cells were located in the light zone (Fig. 2b). This prediction agrees with two previous studies [20, 21].
Fig. 2

Simulations of dynamics of cyanobacterial performance in a 60 mm-radius cylindrical PBR under 50 µE/m2/s surface light intensity. a biomass concentration and dissolved CO2/HCO3 concentration; b glycogen content in biomass and average specific growth rate; c flux through the PEP carboxylase reaction; d flux through the RuBP carboxylase reaction; e fluxes through the TCA cycle (negative flux means the flux direction is reversed); f ratios of autotrophic, heterotrophic, and resting sub-populations; and g light distribution in the bioreactor as a function of time and relative distance r/R, where r is the local position and R is the radius of the bioreactor. The parameters used for simulation are given in Table 1. Abbreviations of metabolites: 3PG 3-phosphoglycerate, AceCoA acetyl-CoA, CIT citrate, MAL malate, OAA oxaloacetate, PEP phosphoenolpyruvate, and RuBP ribulose-1,5-bisphosphate

Next, we tested the sensitivity of biomass production to the mass transfer rate, PBR surface light intensity, and PBR diameter (Fig. 3) . With a light intensity of 100 µE/m2/s and a moderate mass transfer rate of 15 h−1, small PBRs (30 mm radius) could produce 20 g/L of biomass in 21 days. Although such productivity has been experimentally observed in small PBRs [22], it can be hardly achieved in large-size PBRs. As shown by the model, the biomass productivity is highly sensitive to the surface-to-volume ratios of the PBRs, and increasing the PBR diameter reduces biomass productivity dramatically. Hence, to improve biomass production in PBRs, one needs to reduce the surface-to-volume ratios, increase the culture mixing and air flow rate [23], and maintain a sufficient surface light intensity.
Fig. 3

Test of biomass growth performance sensitivity to the mass transfer rate (K La in h−1), surface light intensity (I 0 in µE/m2/s), and the bioreactor geometry (R in mm). Each three-dimensional figure shows biomass growth as a function of time and surface light intensity

Finally, the simulations also demonstrate that, due to random cell movements in PBRs, single cell fluxome may show stochastic changes (Additional file 2: Fig. S2 b, c). Additionally, we tested the sensitivity of biomass growth to circulation time. The model indicates that perturbing the circulation speeds of cell subpopulations did not affect total biomass production by PBRs, provided that the CO2 mass transfer and surface light were constant (Additional file 2: Fig. S3).

Simulation of cyanobacterial oxidative pentose phosphate pathway in a cylindrical PBR

In cyanobacteria, the oxidative pentose phosphate (OPP) pathway and the Calvin cycle operates in opposite directions: The former generates CO2 and NADPH, while the latter consumes CO2 and NADPH. Figure 2d shows that the Calvin cycle had a strong flux in the early growth phase, while the OPP pathway remained silent under light-sufficient conditions (Fig. 4a). In the late growth stage, active fluxes through the OPP pathway appeared (Fig. 4a) due to the self-shading effect. The activity of the OPP pathway increased concurrently with the glycogen consumption rate in darkness (Fig. 4a and Additional file 2: Fig. S4). Thus, an active OPP flux in photoautotrophic cultures is the metabolic response to light deficiency in PBRs. Recent 13C-flux measurements also showed positive OPP fluxes in Synechocystis 6803 PBR cultures [24, 25]. To further confirm our model predictions, we examined the labelling patterns of histidine by growing Synechocystis 6803 with NaH13CO3 and [1-13C] glucose. When glucose was metabolized via the OPP pathway, non-labeled ribose-5-phosphate was generated from [1-13C] glucose [26], which is a precursor to histidine. Therefore, an active OPP pathway was expected to reduce the 13C-enrichment of proteinogenic histidine. Figure 3b shows that the 12C-concentration of histidine was high under low light conditions, supporting the model prediction that light deficiency leads to an active OPP pathway for C6 sugar utilizations.
Fig. 4

The oxidative pentose phosphate (OPP) pathway in cyanobacteria. a Comparison of simulated fluxes through the OPP pathway between wild-type (black bars) and glycogen-rich (white bars) cyanobacteria strains at different biomass concentrations. ‘Asterisk’ means the flux is zero. Compared to the wild type strain, the model assumes that the glycogen-rich strain accumulates five times more glycogen during autotrophic growth and consumes glycogen five times faster during heterotrophic growth. b Relative abundance of histidine labeling profile under different Synechocystis growth phases. Black bars (low light/early growth phase): light intensity of ~50 µE/m2/s; white bars (low light/late growth phase): light intensity of ~50 µE/m2/s; grey bars (high light/late growth phase): light intensity of ~100 µE/m2/s

Investigation of cyanobacterial photosynthesis efficiency in shake flasks

Next, we used the integrated GSM to determine the photosynthesis efficiency of Synechocystis 6803 by minimizing the sum of squared errors between experimental and simulated averaged specific grow rates. We simplified the geometry of the shake flasks into a two-dimensional rectangle (Additional file 2: Fig. S5), and made the local light intensity dependent on the vertical distance from a cell to the light source. The CO2 mass transfer rates in shake flasks were calculated based on Eq. (9). As a consequence, a photosynthesis efficiency of 2.7 % (Additional file 2: Fig. S6) resulted from the best fit of specific growth rates under shake flask cultures (Diamond and circle markers in Fig. 5). Furthermore, this photosynthesis efficiency was used to simulate the growth of Synechococcus elongatus UTEX 2973 (a fast-growing cyanobacterium species) in a column PBR (with 3 % CO2 and under 500 µE/m2/s light intensity) [27]. The model predicted slightly lower specific growth rate than the experimental value (Square marker in Additional file 2: Fig. S5). This difference is possibly due to an increased photosynthesis activity under high CO2 concentrations [28].
Fig. 5

Comparison of experimental and simulation results of cyanobacterial growth rates in shake flasks. Diamond markers represent cultures with different volumes under the same light intensity of ~50 μE/m2/s. The volumes are 50, 100, and 150 mL, respectively, corresponding to diamond markers from top to bottom. Circle markers represent 15 mL cultures growing under different light conditions. The light intensities are ~15, ~25 and ~35 μE/m2/s, respectively, corresponding to circle markers from top to bottom. The specific growth rates were calculated based on OD730 values in the early exponential growth phase. The square marker represents a reported cyanobacterial growth rate in a mini-PBR [27]. The parameters used for simulating cyanobacterial growth are given in Additional file 2: Table S1

Model-based investigation of lactate production by engineered cyanobacteria in PBRs

We further applied the integrated GSM to predict the growth and volumetric d-lactate productivity of engineered cyanobacterial strains, in which a mutated glycerol dehydrogenase was overexpressed for producing optically pure d-lactate [29]. The MOMA algorithm was applied to simulate the metabolism in engineered strains (See Methods). Growth-associated lactate production was assumed (i.e., lactate production was proportional to biomass synthesis). First, we tested the relationship between lactate efflux (v lac ) and specific growth rate (µ) using only the FBA model. Figure 6a shows v lac and µ as functions of the ratio v lac /µ, which denotes the amount of lactate produced per gram of biomass (or mmol lactate/g biomass). Within a wide v lac /µ range, from 0.01 to 100 mmol lactate/g biomass, µ decreased with increasing v lac /µ, but v lac showed a parabolic tendency, peaking at 0.3 mmol/g/h (Fig. 6a). Next, we used the integrated GSM to simulate the cyanobacterial growth and D-lactate production in PBRs at different v lac /µ ratios (Fig. 6b, c). As a result, increasing the v lac /µ ratio led to lower biomass production, which, however, did not necessarily improve the overall d-lactate production. For example, when v lac /µ was increased from 1 to 10 mmol/g, lactate production in PBRs remained the same, but biomass production was significantly diminished. Enhancing the lactate pathway (i.e., increase v lac /µ ratio) can improve lactate production, but excessive overexpression of this pathway may sacrifice biomass growth and impair overall lactate productivity. To resolve this problem, it is desirable to induce the lactate synthesis pathway at late biomass growth phase.
Fig. 6

Simulation results of d-lactate producing cyanobacteria performance. a FBA simulations of d-lactate flux and growth rate as functions of v lac/µ (mmol lactate/g biomass). White markers wild-type strain, black markers glycogen-knockout strain. b, c Simulation of biomass growth (b) and d-lactate production (c) of wild-type cyanobacteria at different v lac/µ ratios in a cylindrical PBR. d, e Simulation of biomass growth (d) and d-lactate production (e) of glycogen-knockout cyanobacteria at different v lac/µ ratios in a cylindrical PBR

Deleting carbon storage in cyanobacteria is one strategy to redirect carbon flux to product synthesis [30]. However, the integrated GSM shows that such a strategy may not offer significant improvements in final lactate productivity in PBRs (Fig. 6a, d, e). This finding is consistent with two recent reports: (1) Glycogen knockout did not enhance lactate productivity under nutrient-sufficient growth conditions [30]. (2) Removal of glycogen in an isobutanol-producing cyanobacterium yielded no benefit in production titer and rate [31]. Possibly, glycogen serves as the carbon and energy reservoir to store the energy and carbon excess flow in the light zone, and this carbon and energy reservoir can maintain redox homeostasis under stressed growth conditions or in darkness [32]. Therefore, deleting glycogen or other carbon storage may impair cyanobacterial survival as well as its resistance to environmental stresses and contaminations.


In this study, a genome-scale FBA model was integrated with information on kinetics, light distribution, and cell movement. Using the integrated GSM, one can simultaneously learn both intracellular information (e.g., flux distributions as functions of time) and extracellular information (e.g., growth curve and nutrient changes in the medium) simultaneously.

In the extracellular domain, the integrated GSM can describe changes in nutrient concentrations, biomass accumulation, and local light intensities. As demonstrated by Fig. 3, cyanobacteria are intrinsically capable of reaching high biomass concentration in PBRs, however, their performance is usually limited by low light availability and low mass transfer rate. To reduce cell self-shading, high surface-to-volume ratio PBRs equipped with thin panel or hollow fibers have been employed [23]. To improve the mass transfer rate, enhancing aeration rates has proved efficient [19]. Better mixing conditions not only lead to better gas transfer rates, but also help maintain more homogeneous conditions for both cells and nutrients.

In the intracellular domain, we observed continuously changing fluxes in the cyanobacterial metabolic network, which were mostly affected by reduced energy and carbon inputs. One interesting finding is the OPP pathway activation as cell cultures get denser. It is a traditional point of view that the Calvin cycle and the OPP pathway are separate systems because the two pathways are reciprocally regulated [33]. In this study, the integrated GSM describes that these two pathways could be employed by two different subpopulations in PBRs at the same time. This simulation explains that a measureable flux through the OPP pathway could be observed in both wild type and engineered cyanobacterial strains via 13C-based flux analysis [24, 25]. In addition, the elementary modes analysis shows that the Calvin cycle and the OPP pathway may function in complementary ways in photoautotrophs, since an active OPP pathway ensures a maintainable flux to triose phosphate synthesis from carbohydrate degradation in low light or darkness [33]. Furthermore, we noticed that an active OPP pathway was always present in the d-lactate producing strain, and it became stronger with increased d-lactate production (Additional file 2: Fig. S7). This indicates that the enhanced OPP activity benefits bio-production by providing more reducing power. In summary, the plasticity of the OPP pathway endows cells with high vitality and energy flexibility [34, 35].

The traditional FBA model usually describes the optimal growth condition, and thus it may not be suitable for suboptimal and heterogeneous cultivation conditions. Our model, on the other hand, is integrated with growth kinetics and a heterogeneous light distribution in PBRs. Hence, the model can not only predict the production yield, titer and rate, but also offer insights into how cells adjust their internal metabolisms to survive under different growth conditions and genetic manipulations. Moreover, the integrated GSM may give more accurate predictions of mutant physiology than GSM alone in bioreactor conditions. For example, the integrated GSM correctly indicates that glycogen knockout may not be an effective strategy to improve PBR lactate production. Lastly, the integrated GSM can reveal real-time variations/dynamics in metabolisms of different subpopulation cells, and thus improve understandings of cellular responses to large-size PBRs.

Nevertheless, our model still has limitations. For example, previous studies have shown that glycogen synthesis could be connected with unknown regulations affecting cyanobacterial viability under stress conditions [32, 36]. However, the integrated GSM may not give the same prediction without further constraints from knowledge of genetic regulations. Additionally, it has been demonstrated that cyanobacteria have circadian behaviors (i.e., their metabolism exhibits day and night rhythms) [37], while our model does not include this property. Finally, some inhibition factors may also influence cyanobacterial growth (e.g., effects of crowding), which are not included in the model. In the future, this model platform should be further improved via additional multi-scale modeling approaches.


This study demonstrates a genome-scale FBA model integrated with kinetics, cell movements, and a light distribution function. With constraints obtained from bioprocess variables, the integrated GSM can not only simulate the dynamic metabolisms in sub-population cells but also predict PBR overall productivity under light and CO2 conditions. The integration of GSMs with PBR modeling can facilitate the development of new cyanobacterial strains for industrial settings.


Cell cultivation

Synechocystis PCC 6803 was cultivated in a modified BG-11 medium [26] at 30 °C and 180 rpm. We first tested the cyanobacterial growth in different culture volumes. In brief, 50, 100, and 150 mL of cell suspensions were cultivated in 250 mL shake flasks under continuous illumination of ~50 µE/m2/s. We also tested the cyanobacterial growth under different light conditions. Specifically, 15 mL of cultures were grown in 150 mL shake flasks under different light intensities (from ~15 to ~35 µE/m2/s). OD730 was used to measure biomass density, and the relationship between the biomass concentration and OD730 was 0.45 × OD730 = Biomass (g/L) [26]. We made duplicate cultures of each condition (n = 2).

13C-Labelling experiment

13C-labeling experiments were performed to determine histidine labeling as evidence of OPP pathway activity under different light conditions. We grew photomixotrophic cultures in BG-11 medium supplied with 2.5 g/L [1-13C] glucose and 4 g/L NaH13CO3 (tracers were purchased from Sigma-Aldrich, Saint Louis, USA). The TBDMS (N-tert-butyldimethylsilyl-N-methyltrifluoroacetamide) method [26] was used to analyze the labeling patterns of proteinogenic histidine. In brief, cells were harvested by centrifugation, and cell pellets were hydrolyzed in 6 mol/L HCl solution at 100 °C for 24 h. The amino acid solution was air-dried and then derivatized by TBDMS (Sigma-Aldrich, USA) at 70 °C for 1 h. A gas chromatograph (GC) (Hewlett-Packard model 7890A; Agilent Technologies, CA) equipped with a DB5-MS column (J&W Scientific, Folsom, CA) and a mass spectrometer (model 5975C; Agilent Technologies, CA) were used for analyzing amino acid labeling profiles. The GC–MS fragment [M-57]+ contains the complete amino acid backbone, and MS data M + 0, M + 1, and M + 2 represent isotopomers with zero, one, and two 13C atoms, respectively.

Flux balance analysis model

The FBA model was modified from the cyanobacterial model iJN878 [15], which has 843 reactions, including photosynthesis and the central carbon metabolism. A complete list of reactions is provided in Additional file 1. The iJN878 model contains a recently discovered γ-aminobutyrate shunt [38] which converts 2-oxoglutarate to succinate in Synechocystis 6803. In our model, two new reactions were added, namely ‘glycogen storage → glycogen [c]’ and ‘d-lactate [c] → d-lactate [external]’, which were respectively used to simulate glycogen storage/consumption and d-lactate production by an engineered cyanobacterial strain [29]. The mathematical description of our FBA model is as follows:
$$\left[ \begin{array}{l} {\text{maximize }}\mu \hfill \\ {\text{subject to }}S \cdot v = 0 \hfill \\ \quad \quad \quad \quad lb \le v \le ub \hfill \\ \quad \quad \quad \quad v_{{CO_{2} }} \le f_{1} \left( {K_{{La}} ,[{\text{CO}}_{2} ],Km} \right) \hfill \\ \quad \quad \quad \quad v_{{photon}} \le f_{2} \left( {l,X,v_{{photon,0}} } \right) \hfill \\ \end{array} \right],$$
where µ represents the specific growth rate, S is the stoichiometric matrix, v represents a vector of flux distribution, and lb and ub represent vectors of the lower and upper boundaries, respectively. Further, f 1 is a function of the mass transfer coefficient K La , dissolved CO2 concentration [CO2], and half-saturation constant for dissolved CO2, Km; f 2 is a function of the cell’s local position l, biomass concentration X, and photon influx on the PBR surface v photon,0 . The linear optimization problem was solved by the MATLAB® (2012b) built-in function ‘linprog’ using the ‘simplex’ algorithm. To estimate the flux distribution in engineered cyanobacterial strains, we used the MOMA (minimization of metabolic adjustment) algorithm [39], which was solved by the MATLAB built-in function ‘quadprog’ using the ‘interior-point-convex’ algorithm:
$$\left[ \begin{array}{l} \text{minimize} \, \frac{1}{2}v^{T} Hv - f_{opt}^{T} v \hfill \\ {\text{subject to }}S \cdot v = 0 \hfill \\ \quad \quad \quad \quad lb \le v \le ub \hfill \\ \quad \quad \quad \quad v_{{CO_{2} }} \le f_{1} \left( {K_{La} , [ {\text{CO}}_{2} ],Km} \right) \hfill \\ \quad \quad \quad \quad v_{photon} \le f_{2} \left( {l,X,I_{0} } \right) \hfill \\ \end{array} \right],$$
where H is a unit matrix, and f opt is the optimal flux distribution of wild-type cyanobacteria. The remaining notations have the same meanings as above.

Moreover, we considered three growth states for cyanobacteria in the FBA model: (1) the ‘light condition’: an autotrophic sub-population in the light zone, (2) the ‘dark condition’: a heterotrophic sub-population in the dark zone, where the photon influx is below 0.4 mmol/g/h (under which the cyanobacterial growth rate is lower than the heterotrophic growth rate in darkness) and glycogen is consumed at a rate of 0.01 mmol/g/h [40] to maintain minimal growth, and (3) the ‘resting condition’: a glycogen-depleted sub-population with no active fluxes in the dark zone. To improve the calculation efficiency, we built a database containing all the flux distributions in response to different photon influxes (Additional file 1). By having such a database, we could directly use pre-calculated fluxome from the database according to culture conditions in PBRs. Thereby, we did not need to redo flux calculations at each time interval during new simulations.

Simulation of cyanobacterial growth via integrating FBA, kinetics, and cell movements

Figure 1 shows our modeling algorithm. To simulate biomass growth as a function of time, we divided the entire time period into finite intervals of 0.002 h (Additional file 2: Fig. S1 shows that further decreasing the interval period did not change the simulation results). In each time interval, a simplified sinusoid equation [41] was used to estimate the cell location in a well-mixed PBR:
$$l = \frac{r}{2} - \frac{r}{2}\cos (\frac{2\pi }{{f_{r} }}t),$$
where l is the shortest distance between the PBR surface and the cell local position, in mm; r is the radius or thickness of the PBR, in mm; f r represents the cyanobacteria circulation frequency, in h; and t is time, in h. Because cell circulation frequencies in PBRs vary from cell to cell, stochastic effects are induced on a single cell’s metabolism. In fact, the random movements of cells in PBRs have been measured and simulated, and, in the present study, are described by a probability distribution function [42]. In our model platform, we distinguished cell populations with different circulation times (Additional file 2: Fig. S2-3). Thus, the whole culture was considered to be comprised of twelve sub-populations instead of a plethora of cyanobacterial cells. Based on cell locations and the biomass concentrations, we calculated the local light intensity [43]:
$$\frac{{v_{photon} }}{{v_{photon,0} }} = \left( {\frac{1}{{(0.0216 \cdot l + 1)^{1.54} (0.130 \cdot X \cdot + 1)^{1.18} }}} \right),$$
where v photon and v photon,0 represent the local and surface photon influxes, respectively, in mmol/g/h; and X is biomass concentration, in g/L.
The CO2 uptake flux was described by a Michaelis–Menten equation:
$$v_{{{\text{CO}}_{2} }} = v_{{{\text{CO}}_{2} \,,\hbox{max} }} \frac{{[{\text{CO}}_{2} ]}}{{K_{m} + [{\text{CO}}_{2} ]}},$$
where v CO2,max is the maximum uptake rate of dissolved CO2/HCO3 , K m is the half-saturation constant, and [CO2] represents dissolved CO2 concentration. This study assumed that pH was constant at 8.0, and that the dissolved CO2 and cell culture were homogeneous in PBRs. Because we assumed that cell metabolism was pseudo-steady in each interval [44], the FBA model could use linear optimization to profile the intracellular fluxes constrained by light and carbon input fluxes (Eqs. 4 and 5). The FBA model then predicted the growth rates, glycogen synthesis rates, and CO2 uptake rates of cell populations with different circulation times in PBRs. Those values were averaged based on the probability distribution function (Additional file 2: Fig. S2a):
$$\begin{aligned} \mu _{{app}} & = \sum\limits_{{i = 1}}^{n} {P_{i} \mu _{i} } \\ v_{{{\text{CO}}_{2} ,app}} & = \sum\limits_{{i = 1}}^{n} {P_{i} v_{{{\text{CO}}_{2} ,i}} } \\ v_{{glycogen,app}} & = \sum\limits_{{i = 1}}^{n} {P_{i} v_{{glycogen,i}} } , \\ \end{aligned}$$
where P i is the fraction of ith cell population (Additional file 2: Fig. S2a), µ app is the apparent specific growth rate in PBRs, v CO2,app is the apparent CO2 uptake rate, and v glycogen,app is the overall glycogen production rate.
For the kinetic model, we used ordinary differential equations (ODEs) to describe changes in biomass production, glycogen accumulation, dissolved CO2/HCO3 concentrations, and so forth.
$$\begin{aligned} \frac{dX}{dt} &= \mu_{app} \cdot X - K_{d} \cdot X \hfill \\ \frac{d[glycogen]}{dt} &= v_{glycogen,app} \cdot X - \beta \cdot K_{d} \cdot X \hfill \\ \frac{{d[{\text{CO}}_{2} ]}}{dt} &= K_{La} \cdot ([{\text{CO}}_{2} ]^{*} - [{\text{CO}}_{2} ]) + v_{{{\text{CO}}_{2} ,app}} \cdot X, \hfill \\ \end{aligned}$$
The ODEs were resolved in their numerical discrete form (Euler-like integration scheme):
$$\begin{aligned} X_{i + 1} = X_{i} + \mu_{app} \cdot X_{i} \cdot \Delta t - K_{d} \cdot X_{i} \cdot \Delta t \hfill \\ [glycogen]_{i + 1} = [glycogen]_{i} + v_{glycogen,app} \cdot X_{i} \cdot \Delta t - \beta \cdot K_{d} \cdot X_{i} \cdot \Delta t \hfill \\ [{\text{CO}}_{2,i + 1} ] = [{\text{CO}}_{2,i} ] + K_{La} \cdot ([{\text{CO}}_{2} ]^{*} - [{\text{CO}}_{2,i} ]) \cdot \Delta t + v_{{{\text{CO}}_{2} ,app}} \cdot X_{i} \cdot \Delta t \hfill \\ \end{aligned}$$
where i and i + 1 represent the current and next intervals, respectively; Δt is the time interval (0.002 h); K d is the death rate, in h−1; β is the glycogen composition in the biomass, in mmol/g; [glycogen] is the overall glycogen concentration in the PBR, in mmol/L; K La is the mass transfer rate of CO2, in h−1; [CO2] represents the dissolved CO2 and HCO3 concentrations, in mmol/L; [CO2]* is the combined concentrations of dissolved CO2 and HCO3 in equilibrium with atmospheric CO2 (0.039 %, v/v), in mmol/L; and µ app , v glycogen,app , and v CO2,app are fluxes determined previously. The updated values of the biomass concentration, dissolved CO2 concentration, etc., were then used to constrain the FBA model in the next interval. The kinetic parameters are given in Table 1. For shake flask conditions, K La was determined by the following equation: [45]
Table 1

List of parameters used to simulate the growth and metabolic fluxes of cyanobacteria growing in a cylindrical PBR



Value (range)



K La

Mass transfer rate of CO2

10a (3–15)



K d

Death rate





Radius of PBR



Similar to the reactor used in Reference [42]

I 0

Surface light intensity



Equivalent to ~50 µE/m2/sb


Half-saturation constant of CO2 uptake rate





Medium pH



BG-11 medium


Initial biomass concentration



Equivalent to an OD730 of ~0.2


Initial concentration of dissolved CO2 and HCO3 (in equilibrium with air)




aThe value is used in Figs. 2, 3 and 5; ba photosynthesis efficiency of 1.5 % is assumed in Figs. 2, 3 and 5; and ccalculation is based on experimental conditions (Additional file 2)

$${\it{\text{K}}}_{\text{La}} { \,= 0} . 0 3 2 \times {\text{N}}\left( {\frac{\text{V}}{\text{L}}} \right)^{ 0. 8 4 5} ,$$
where N is the rotation speed, in rpm; V is the shake flask volume, in mL; and L is the culture volume, in mL. MATLAB code of the integrated GSM is provided in Additional file 3 and 4.


Author’s contributions

YJT and LH initialized this project. LH and SW built the integrated GSM. LH, NW, and AR performed the experiments. All authors read and approved the final manuscript.


We would like to thank Dr. Le You for valuable discussions for this project. We are also grateful for the funding supports from a US DOE grant (DESC0012722) and an NSF Grant (CBET 1438125).

Competing interests

The authors declare that they have no competing interests.

Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated.

Authors’ Affiliations

Department of Energy, Environmental and Chemical Engineering, Washington University
Department of Mechanical Engineering and Materials Science, Washington University
Department of Biochemistry and Molecular Biology, College of Wooster


  1. Janssen M, Tramper J, Mur LR, Wijffels RH. Enclosed outdoor photobioreactors: light regime, photosynthetic efficiency, scale-up, and future prospects. Biotechnol Bioeng. 2003;81:193–210.View ArticleGoogle Scholar
  2. Wu X, Merchuk JC. A model integrating fluid dynamics in photosynthesis and photoinhibition processes. Chem Eng Sci. 2001;56:3527–38.View ArticleGoogle Scholar
  3. Béchet Q, Shilton A, Guieysse B. Modeling the effects of light and temperature on algae growth: state of the art and critical assessment for productivity prediction during outdoor cultivation. Biotechnol Adv. 2013;31:1648–63.View ArticleGoogle Scholar
  4. Cornet JF, Dussap CG. A simple and reliable formula for assessment of maximum volumetric productivities in photobioreactors. Biotechnol Prog. 2009;25:424–35.View ArticleGoogle Scholar
  5. Takache H, Christophe G, Cornet JF, Pruvost J. Experimental and theoretical assessment of maximum productivities for the microalgae Chlamydomonas reinhardtii in two different geometries of photobioreactors. Biotechnol Prog. 2010;26:431–40.Google Scholar
  6. Cornet J, Dussap C, Dubertret G. A structured model for simulation of cultures of the cyanobacterium Spirulina platensis in photobioreactors: I. Coupling between light transfer and growth kinetics. Biotechnol Bioeng. 1992;40:817–25.View ArticleGoogle Scholar
  7. Orth JD, Thiele I, Palsson BO. What is flux balance analysis? Nat Biotechnol. 2010;28:245–8.View ArticleGoogle Scholar
  8. Lerman JA, Chang RL, Hyduke DR, Palsson BØ. Genome-scale models of metabolism and gene expression extend and refine growth phenotype prediction. Mol Syst Biol. 2013;9:693–705.Google Scholar
  9. Becker SA, Feist AM, Mo ML, Hannum G, Palsson BØ, Herrgard MJ. Quantitative prediction of cellular metabolism with constraint-based models: the COBRA Toolbox. Nat Protoc. 2007;2:727–38.View ArticleGoogle Scholar
  10. Ranganathan S, Suthers PF, Maranas CD. OptForce: an optimization procedure for identifying all genetic manipulations leading to targeted overproductions. PLoS Comput Biol. 2010;6:e1000744.View ArticleGoogle Scholar
  11. Teusink B, Wiersma A, Jacobs L, Notebaart RA, Smid EJ. Understanding the adaptive growth strategy of Lactobacillus plantarum by in silico optimisation. PLoS Comput Biol 2009; 5:e1000410.View ArticleGoogle Scholar
  12. Schirmer A, Rude MA, Li X, Popova E, Del Cardayre SB. Microbial biosynthesis of alkanes. Science. 2010;329:559–62.View ArticleGoogle Scholar
  13. Bandyopadhyay A, Stöckel J, Min H, Sherman LA, Pakrasi HB. High rates of photobiological H2 production by a cyanobacterium under aerobic conditions. Nat Commun. 2010;1:139.View ArticleGoogle Scholar
  14. Berla BM, Saha R, Immethun CM, Maranas CD, Moon TS, Pakrasi HB. Synthetic biology of cyanobacteria: unique challenges and opportunities. Front Microbiol. 2013;4:246.View ArticleGoogle Scholar
  15. Nogales J, Gudmundsson S, Knight EM, Palsson BO, Thiele I. Detailing the optimality of photosynthesis in cyanobacteria through systems biology analysis. Proc Natl Acad Sci USA. 2012;109:2678–83.View ArticleGoogle Scholar
  16. Blankenship RE, Tiede DM, Barber J, Brudvig GW, Fleming G, Ghirardi M, Gunner M, Junge W, Kramer DM, Melis A. Comparing photosynthetic and photovoltaic efficiencies and recognizing the potential for improvement. Science. 2011;332:805–9.View ArticleGoogle Scholar
  17. Zhu X-G, Long SP, Ort DR. What is the maximum efficiency with which photosynthesis can convert solar energy into biomass? Curr Opin Biotechnol. 2008;19:153–9.View ArticleGoogle Scholar
  18. Melis A. Solar energy conversion efficiencies in photosynthesis: minimizing the chlorophyll antennae to maximize efficiency. Plant Sci (Amsterdam, Neth). 2009;177:272–80.Google Scholar
  19. Kazim SA. Experimental and empirical correlations for the determination of the overall volumetric mass transfer coefficients of carbon dioxide in stirred tank bioreactors. The University of Western Ontario. 2012. p. 48.Google Scholar
  20. Aikawa S, Izumi Y, Matsuda F, Hasunuma T, Chang J-S, Kondo A. Synergistic enhancement of glycogen production in Arthrospira platensis by optimization of light intensity and nitrate supply. Bioresour Technol. 2012;108:211–5.View ArticleGoogle Scholar
  21. Hasunuma T, Kikuyama F, Matsuda M, Aikawa S, Izumi Y, Kondo A. Dynamic metabolic profiling of cyanobacterial glycogen biosynthesis under conditions of nitrate depletion. J Exp Bot 2013;64(10):2943–54.View ArticleGoogle Scholar
  22. Pirt SJ, Lee YK, Walach MR, Pirt MW, Balyuzi HH, Bazin MJ. A tubular bioreactor for photosynthetic production of biomass from carbon dioxide: design and performance. J Chem Technol Biotechnol. 1983;33:35–58.View ArticleGoogle Scholar
  23. Lee C-G, Palsson BØ. High-density algal photobioreactors using light-emitting diodes. Biotechnol Bioeng. 1994;44:1161–7.View ArticleGoogle Scholar
  24. Young JD, Shastri AA, Stephanopoulos G, Morgan JA. Mapping photoautotrophic metabolism with isotopically nonstationary 13C flux analysis. Metab Eng. 2011;13:656–65.View ArticleGoogle Scholar
  25. Xiong W, Morgan JA, Ungerer J, Wang B, Maness P-C, Yu J. The plasticity of cyanobacterial metabolism supports direct CO2 conversion to ethylene. Nat Plants 2015; 1.Google Scholar
  26. You L, Berla B, He L, Pakrasi HB, Tang YJ. 13C-MFA delineates the photomixotrophic metabolism of Synechocystis sp. PCC 6803 under light- and carbon-sufficient conditions. Biotechnol J. 2014;9:684–92.View ArticleGoogle Scholar
  27. Yu J, Liberton M, Cliften PF, Head RD, Jacobs JM, Smith RD, Koppenaal DW, Brand JJ, Pakrasi HB. Synechococcus elongatus UTEX 2973, a fast growing cyanobacterial chassis for biosynthesis using light and CO2. Sci Rep 2015; 5.Google Scholar
  28. Levitan O, Rosenberg G, Setlik I, Setlikova E, Grigel J, Klepetar J, Prasil O, Berman-Frank I. Elevated CO2 enhances nitrogen fixation and growth in the marine cyanobacterium Trichodesmium. Glob Chang Biol. 2007;13:531–8.View ArticleGoogle Scholar
  29. Varman AM, Yu Y, You L, Tang YJ. Photoautotrophic production of d-lactic acid in an engineered cyanobacterium. Microb Cell Fact. 2013;12:117.View ArticleGoogle Scholar
  30. van der Woude AD, Angermayr SA, Puthan Veetil V, Osnato A, Hellingwerf KJ. Carbon sink removal: Increased photosynthetic production of lactic acid by Synechocystis sp. PCC6803 in a glycogen storage mutant. J Biotechnol. 2014;184:100–2.View ArticleGoogle Scholar
  31. Li X, Shen C, Liao J. Isobutanol production as an alternative metabolic sink to rescue the growth deficiency of the glycogen mutant of Synechococcus elongatus PCC 7942. Photosynth Res. 2014;120:301–10.View ArticleGoogle Scholar
  32. Gründel M, Scheunemann R, Lockau W, Zilliges Y. Impaired glycogen synthesis causes metabolic overflow reactions and affects stress responses in the cyanobacterium Synechocystis sp. PCC 6803. Microbiology. 2012;158:3032–43.View ArticleGoogle Scholar
  33. Poolman MG, Fell DA, Raines CA. Elementary modes analysis of photosynthate metabolism in the chloroplast stroma. Eur J Biochem. 2003;270:430–9.View ArticleGoogle Scholar
  34. Wu S, Huang A, Zhang B, Huan L, Zhao P, Lin A, Wang G. Enzyme activity highlights the importance of the oxidative pentose phosphate pathway in lipid accumulation and growth of Phaeodactylum tricornutum under CO2 concentration. Biotechnol Biofuels. 2015;8:78.View ArticleGoogle Scholar
  35. Wasylenko TM, Ahn WS, Stephanopoulos G. The oxidative pentose phosphate pathway is the primary source of NADPH for lipid overproduction from glucose in Yarrowia lipolytica. Metab Eng. 2015;30:27–39.View ArticleGoogle Scholar
  36. Xu Y, Tiago Guerra L, Li Z, Ludwig M, Charles Dismukes G, Bryant DA. Altered carbohydrate metabolism in glycogen synthase mutants of Synechococcus sp. strain PCC 7002: Cell factories for soluble sugars. Metab Eng. 2013;16:56–67.View ArticleGoogle Scholar
  37. Kondo T, Strayer CA, Kulkarni RD, Taylor W, Ishiura M, Golden SS, Johnson CH. Circadian rhythms in prokaryotes: luciferase as a reporter of circadian gene expression in cyanobacteria. Proc Natl Acad Sci USA. 1993;90:5672–6.View ArticleGoogle Scholar
  38. Xiong W, Brune D, Vermaas WF. The γ-aminobutyric acid shunt contributes to closing the tricarboxylic acid cycle in Synechocystis sp. PCC 6803. Mol Microbiol. 2014;93:786–96.View ArticleGoogle Scholar
  39. Segre D, Vitkup D, Church GM. Analysis of optimality in natural and perturbed metabolic networks. Proc Natl Acad Sci USA. 2002;99:15112–7.View ArticleGoogle Scholar
  40. Knoop H, Gründel M, Zilliges Y, Lehmann R, Hoffmann S, Lockau W, Steuer R. Flux balance analysis of cyanobacterial metabolism: the metabolic network of Synechocystis sp. PCC 6803. PLoS Comput Biol. 2013;9:e1003081.View ArticleGoogle Scholar
  41. Wu X, Merchuk JC. Simulation of algae growth in a bench-scale bubble column reactor. Biotechnol Bioeng. 2002;80:156–68.View ArticleGoogle Scholar
  42. Luo H-P, Al-Dahhan MH. Verification and validation of CFD simulations for local flow dynamics in a draft tube airlift bioreactor. Chem Eng Sci. 2011;66:907–23.View ArticleGoogle Scholar
  43. Katsuda T, Arimoto T, Igarashi K, Azuma M, Kato J, Takakuwa S, Ooshima H. Light intensity distribution in the externally illuminated cylindrical photo-bioreactor and its application to hydrogen production by Rhodobacter capsulatus. Biochem Eng J. 2000;5:157–64.View ArticleGoogle Scholar
  44. Feng X, Xu Y, Chen Y, Tang YJ. Integrating flux balance analysis into kinetic models to decipher the dynamic metabolism of Shewanella oneidensis MR-1. PLoS Comput Biol. 2012;8:e1002376.View ArticleGoogle Scholar
  45. Nikakhtari H, Hill GA. Modelling oxygen transfer and aerobic growth in shake flasks and well-mixed bioreactors. Can J Chem Eng. 2005;83:493–9.View ArticleGoogle Scholar
  46. Kim HW, Vannela R, Zhou C, Rittmann BE. Nutrient acquisition and limitation for the photoautotrophic growth of Synechocystis sp. PCC6803 as a renewable biomass source. Biotechnol Bioeng. 2011;108:277–85.View ArticleGoogle Scholar
  47. Badger MR, Andrews TJ. Photosynthesis and inorganic carbon usage by the marine cyanobacterium, Synechococcus sp. Plant Physiol. 1982;70:517–23.View ArticleGoogle Scholar


© He et al. 2015