A kinetic metabolic study of lipid production in Chlorella protothecoides under heterotrophic condition

Abstract

Background

Microalgae have been proposed as potential platform to produce lipid-derived products, such as biofuels. Knowledge on the intracellular carbon flow distribution may identify key metabolic processes during lipid synthesis thus refining culture/genetic strategies to maximize cell lipid productivity. A kinetic metabolic model simulating cell metabolic behavior and lipid production was first applied in the microalgae platform Chlorella protothecoides under heterotrophic condition. It combines both physiology and flux information in a kinetic approach. Cell nutrition, growth, lipid production and almost 30 metabolic intermediates covering central carbon metabolism were included and simulated.

Results

Model simulations were shown to adequately agree with experimental data, which is suggesting that the proposed model copes with Chlorella protothecoides cells’ biology. The dynamic metabolic flux analysis using the model showed a reversible starch flux from accumulation to decomposing when glucose reached depletion, while net lipid flux shows a quasi-constant rate. The sensitive flux parameters on starch and lipid metabolism suggested that starch synthesis is the major competing pathway that affects lipid accumulation in C. protothecoides. Flux analysis also demonstrated that high lipid yield under heterotrophic condition is accompanied with high lipid flux and low TCA activity. Meanwhile, the dynamic flux distribution also suggests a relatively constant ratio of glucose distributed to biomass, lipid, starch, nucleotides as well as pentose phosphate pathway.

Conclusion

The model described not only experimental data, but also unraveled intracellular carbon flow distribution and identify key metabolic processes during lipid synthesis. Most of the metabolic kinetics also showed statistical significance for metabolic mechanism. Therefore, this study unravels the mechanisms of the glucose impact on the dynamic carbon flux distribution, thus improving our understanding of the links between carbon fluxes and lipid metabolism in C. protothecoides.

Background

Microalgae have been proposed as potential platform to produce lipid-derived products, such as biofuels. Previous studies have shown that most of microalgae cells accumulate lipids when the cell division is blocked or inhibited while carbon can still continue to be fixed (such as nitrogen shortage stress conditions), which results in reduced biomass growth and in tum constraints the total lipid yield [1,2,3,4]. Thus, there is a contradiction between lipid content and algae growth, which limits high lipid yield. Knowledge on the intracellular carbon flow distribution may identify key metabolic processes during lipid synthesis thus refining culture/genetic strategies to maximize cell lipid productivity [5,6,7]. The flux balance analysis (FBA) approach, which is based on pseudo steady-state approximation, has been applied to uncover the black-box of intracellular metabolism under steady state [8, 9]. It has been used in microalgae biosystems such as Arghrospia platensis, Synechocystis sp. PCC 6803 [10], Chlamydomonas reinhardtii [11, 12], Chlorella protothecoides [13] and Chlorella sp. [14]. For instance, in Chlorella sp., a shift in intracellular flux distribution was predicted during transition from nutrient sufficient phase to nutrient starvation phase of growth [14]. Another appealing modeling approach, which allows simulating a culture’s dynamics, is based on a kinetic transient-type approach [15, 16]. Such kinetic metabolic models describe cell dynamic behavior by inducing enzyme kinetics. An underestimated potential output of these kinetic metabolic models relies in their capacity to perform dynamic metabolic flux analysis from which key metabolic processes can be examined while assessing in silico hypothesis of genetic engineering and/or culture conditions management strategies [17]. However, to the best of our knowledge, kinetic metabolic model application in microalgae is relatively new. There are very few studies reporting dynamic metabolomics data of microalgae and even less on the development of mathematical models to describe cell metabolic dynamics [18].

In the present work, a kinetic metabolic model describing Chlorella protothecoides cellular metabolism was developed to describe heterotrophic culture mode. Unlike most dynamic models coping with algae physiology or steady-state metabolic level models (FBA), the current model combines both physiology and flux information in a kinetic approach. Cell growth, lipid production and almost 30 metabolic intermediates covering glycolysis, pentose phosphate pathway and TCA cycle and energetic metabolism were included and simulated. Multiple Michaelis–Menten equations are used to introduce metabolic kinetics of reaction rates and the Monod equation is used to describe the cell specific growth state, mass balances for each intermediate were considered in a dynamic profile. It can thus be used as an in silico platform for characterizing the cell lines as well as to search for ‘‘optimal’’ culture strategy through identify key metabolic processes during lipid synthesis.

Materials and methods

Algae stain and culture conditions

Details about algae species and culture conditions can be found in a previous work [19]. Briefly, Chlorella protothecoides (the Culture Collection of Alga at the University of Texas) culture in the dark was carried out in 2.8-L glass flasks with 10 g L⁻¹ glucose as the carbon source and the modified basal medium (MBM), thus imposing a strict heterotrophic metabolism. Glucose concentration in the medium was analyzed by a biochemistry analyzer (YSI Life Science, 2700 select, Ohio, USA). Intracellular metabolites extraction was performed as described in previous work [19] and their quantification was carried out by UPLC/MS/MS system (1290 model, Agilent Technologies, Santa Clara, CA, USA), starch analysis was performed using a starch assay kit (Sigma-Aldrich, St. Louis, MO, USA). Total lipid quantification was done according to Drochiou’s method and described in previous work [19].

Model development

Model structure

A kinetic metabolic model was developed to describe the central carbon metabolism of a microalgae platform, including glycolysis, TCA (tricarboxylic acid) cycle, pentose phosphate pathway, total lipid synthesis, starch synthesis, amino acids metabolism, energy metabolism and biomass synthesis. The metabolic network (Fig. 1) was first built according to databases such as KEGG, MetaCyc, DiatomCyc, BioCyc as well as from literature [18, 20, 21]. In this work, Chlorella protothecoides cells were considered as a unique compartment with no specific intracellular compartments such as mitochondria, chloroplast, vacuoles, vesicles and nucleus. Energy metabolism was considered as a global reaction where de novo synthesis and substrate level phosphorylation were combined in a unique pathway. Reversible reactions involving storage carbon such as starch and lipid catabolism were described. The stoichiometry of the biochemical reactions of the network is based on the flux balance analysis on Chlorella protothecoides [18]. A full list of the model reactions and reactions stoichiometry is listed in Table 1.

Fig. 1

The model metabolic network for heterotrophic Chlorella protothecoides

Table 1 Reactions of a metabolic network

The Michaelis–Menten kinetic equation is used to describe each flux rate (Table 2). The cells specific growth rate (equation no. 30 in Table 2), accounting for biomass synthesis from precursors of the major cell constituents such as RX (R5P and X5P), G6P, PYR, AcCOA and total lipid. RX is normally used to synthesize nucleotides, DNA and RNA; G6P leads to organic phosphates providing energy for maintenance and metabolism; PYR is feeding amino acids metabolism which leads to protein formation; AcCOA is the precursor of fatty acids; while the lipid pool is the main contributor to cell mass accumulation in the algae platform.

Table 2 Kinetic equations of the metabolites fluxes in the model

For each metabolite in the model, a mass balance equation (Eq. 1) included the sum of all the input and output fluxes minus cell dilution effect from the cell division phenomenon. The specific consumption of precursors to cell mass synthesis was also considered in the precursors’ mass balances.

$$\frac{{\left. {d[M_{i} } \right]}}{{{\text{d}}t}} = \left( {\mathop \sum \limits_{m = 1}^{a} V_{input_i} - \mathop \sum \limits_{n = 1}^{b} V_{output_j} - \left[ {M_{i} } \right] \times V_{growth} - c_{i} \times V_{growth} } \right)$$

(1)

where M_i is the concentration of each metabolite at time t, a is the input flux number, b is the output flux number, V_input and V_output are the flux rates at each metabolite node. $c_{i}$ is the stoichiometric coefficient for biomass precursors, V_growth is the specific growth rate. So $\left[ {M_{i} } \right] \times V_{growth}$ is the cell dilution term and $c_{i} \times V_{growth}$ is the growth contribution term for biomass precursors (G6P, RX, PYR, Lipid, AcCOA).

All the ordinary differential equations of metabolites were provided as in Additional file 1: Table S1. The ordinary equation for biomass was described by $\frac{dX}{dt} = V_{growth} *X$, listed in Additional file 1: Table S1 (Eq. 22). Where $V_{growth}$ is the specific growth rate, X is the biomass concentration. The cell specific growth rate $V_{growth}$ was described by Eq. 30 in Table 2. A program generating automatically the Matlab code and differential equations for a specific metabolic network of reactions was developed and used (S. Peres and M. Jolicoeur, unpublished). The ordinary differential equations system were performed using Matlab (the MathWorks Inc., Natick, MA, USA) with the “ode23” solver to get the model simulations result.

Model parameters estimation

The model has 77 parameters, which include 34 maximum flux rates, 38 enzyme half-saturated constants, and 5 growth coefficients for the 5 growth precursors contributing to biomass synthesis. Initial metabolite concentrations (i.e. at t = 0) were taken from experimental data, which include 24 intracellular metabolites distributed in 8 pathways covering 30 metabolic reactions [19] or from literature (Table 3).

Table 3 State variables description and initial conditions

Initial kinetic parameter values (V_max, K_m) for each flux and enzymes in the model were taken within ranges found in the enzyme database BRENDA (http://www.brenda-enzymes.org), and the respective units (mmol/L) were converted to comply with the model (mmol/gDW) by dividing 10 gDW/L biomass obtained in our culture. First estimates of maximal flux rates (V_max) have been calculated from experimental data [19], or from BRENDA. Model parameter values were determined following the method proposed in Rizzi et al. [22]. Briefly, the time course of each metabolite with experimental concentration data were defined as fixed mathematical functions, enabling the procedure for parameter values optimization to focus first on non-measured metabolites. In the present case of a high number of parameters, this approach allows accelerating parameter values identification. An objective function (Eq. 2), defined as the weighted sum of squared residues between experimental data and simulated values for each metabolites m at time k, where the weight is the experimental data for each state variable, was used to quantify simulation error.

$${ \text{min} }\left( {\mathop \sum \limits_{i = 1}^{N} \mathop \sum \limits_{t = 1}^{T} \left( {\frac{{M_{i,t}^{exp} - M_{i,t}^{sim} }}{{M_{i,t}^{exp} }}} \right)^{2} } \right)$$

(2)

Based on this objective function, a sensitivity analysis of model parameters was performed to identify the sensitive ones in order to avoid over-parameterization, by then keeping constant non-sensitive parameters. Sensitivity analysis was performed by changing each parameter from − 70 to + 150% one at a time while holding others constant. From the Matlab optimization toolbox, the “linsqurfit” sub-routine was used to identify model parameter values. This process of parameter calibration was continued until minimizing the objective function, i.e. the simulated results closely following experimental data. Final parameter values of the model are shown in Table 4. Confidence intervals of estimated parameters were evaluated using the Matlab sub-routine “nlparci.m” (Table 4). It is clear there is no unique solution for parameter values in such an underdetermined system.

Table 4 Parameter values and 95% confidence intervals of the highly sensitive parameters

Results and discussion

Model simulates algae cell behavior under heterotrophic condition

The final calibrated model adequately simulates the experimental data (Fig. 2). Cell growth, as well as extracellular metabolites such as glucose and glycine are closely simulated. More importantly, total lipids and starch as the main products were also simulated adequately. The model simulation also followed closely the dynamics of intracellular metabolites, which distributed in glycolysis, PPP pathway, TCA cycle as well as energy metabolism. These results thus confirm the model structure as well as its calibrated kinetic parameters to simulate algae cells metabolism and products accumulation dynamics. Indeed, in this work, both the experimental data and model simulations show glycolysis and PPP pathways being more affected by glucose supply while TCA metabolism, which is fed by both carbon and nitrogen metabolisms, seems more robust to perturbations such as extracellular glucose depletion. The intracellular nutrition storage pool in the form of TCA cycle metabolites seemed to maintain biomass in the later growth phase. From both simulation and experimental data, algae biomass still accumulates while glycine and other amino acids pool (AA) reached values under the detection limits. However, this phenomenon was only observed for nitrogen sources since cell biomass growth stopped simultaneously to glucose depletion. This intracellular nutrients management phenomenon has also been modeled and proved in phytoplankton and plant cells [23]. Where the model premises were based on observations that cell growth continued after the exhaustion of external nitrogen pool, being then supported by the consumption of intracellular nitrogen pools such as chlorophyll molecules.

Fig. 2

Simulation result versus experimental data for Chlorella protothecoides heterotrophic behavior (experimental data: open squares, simulated data: solid line. Experimental data were taken from a previous work [19], where the error bars represent the standard error of triplicates data.)

Dynamic metabolic flux analysis reveals lipid and cellular metabolic behavior in Chlorella protothecoides

Considering all the above, it is thus clear that the model structure allows simulating heterotrophic Chlorella protothecoides cell behavior. The model was then taken as an in silico tool and perform a dynamic metabolic flux analysis estimating flux distribution. A dynamic metabolic flux analysis was performed from model simulation (Fig. 3). Looking at glucose flux (V_HK), glycine flux (V_GHMT) as well as cell specific growth rate (V_growth) (Fig. 3a), it is clear that cell growth proceeds simultaneously to carbon source uptake, but not proportionally to nitrogenous source uptake. Interestingly and as previously discussed for glycine concentration, glycine flux (V_GHMT) ceased more than 24 h prior to growth cessation. Fluxes of PPP pathway and starch synthesis (Fig. 3c, d) originate from G6P are partially affected in some extent by glucose flux (Fig. 3a). For instance, model simulation V_PGM flux showed being reversible from accumulation to decomposing at around day 3, where glucose reached depletion. This suggests that starch, which is an intracellular carbon storage pool, rapidly responds to a low carbon source level threshold, contributing providing continuous carbon flow feeding cell metabolism and maintenance. However, as an alternative carbon storage pool, net lipid flux shows a quasi-constant rate, composed of a synthesis flux (V_FASN) that was slightly affected at glucose depletion and two catabolic fluxes (V_Lipase and V_GPAT) which stayed quite constant (11.87–12.04 mmol gDW⁻¹ day⁻¹ and 0.05 mmol gDW⁻¹ day⁻¹ respectively) (Fig. 3e). Interestingly, TCA cycle fluxes (V_ISOD, V_SDH) (Fig. 3f) exhibited a minimum value at glucose depletion, for increasing thereafter. As previously mentioned, the TCA cycle is closely related to carbon and nitrogen metabolism, it seems after carbon source depletion, some carbon and nitrogen dependent compounds (such as pigment) stopped synthesis, which squeezed the intracellular nitrogen source flux goes to TCA cycle. As in heterotrophic culture, C. protothecoides represents yellowish because of carotene content is higher than chlorophyll [24]. However, after glucose depletion, we found the color of culture turns from yellow to green, and the carotene gets to degrade. Moreover, CS flux dynamics closely follows the lipid synthesis flux although it’s quite low compared with lipid synthesis. As CS is competing the same substrate from FASN, the flux of these two enzymes are quite dependent on the concentration of AcCOA, which is in agreement from model prediction.

Fig. 3

Flux rates of all the reactions in the model system. a Nutrition fluxes and growth rate; b glycolysis fluxes; c PPP pathway fluxes; d starch synthesis fluxes; e lipid synthesis fluxes; f TCA cycle fluxes. All flux units were in mmol gDW⁻¹ day⁻¹

A closer view of flux rates were estimated at 48 h before glucose depletion in the exponential growth phase. For comparison purposes, all the flux values were normalized to an uptake flux of 100 mmol g⁻¹DW h⁻¹ glucose (Fig. 4). Flux results agree with that reported in the previous reports [13, 18]. For example, in [18], who performed a flux balance analysis at steady state for Chlorella sp. under heterotrophic condition, with a GPI flux of 66.28 mmol g⁻¹ DW h⁻¹ (leading to glycolysis), G6PDH of 12.15 (leading to PPP pathway) and PGM of 13.83 (leading to starch), compared to 49.8 mmol g⁻¹ DW h⁻¹, 32.04 and 17.3 respectively in our work. The total flux to G6P obtained from our model is of 92.26 mmol g⁻¹ DW h⁻¹ compared that of 99.22 in literature. The net flow From F6P to GD (PFK minus FBPase) was of 73.47 mmol g⁻¹DW h⁻¹ compared to 70.35 mmol g⁻¹ DW h ⁻¹ (from F6P to GAP), and the flux from GD to PEP was of 150.51 mmol g⁻¹ DW h⁻¹ in our model versus 148.25 (from GAP to PEP) in literature. The fluxes of nucleotides synthesis (from RX to ADP) was of 1.36 mmol g⁻¹DW h⁻¹ compared to 0.67 mmol g⁻¹DW h⁻¹ (from PRPP to DNA and RNA). Biomass synthesis rate was of 9.19 compared to 7.36 in [18].

Fig. 4

Flux distribution under heterotrophic cultivation at exponential phase (48 h) and lipid peak time point (72 h). All the flux values were normalized to an uptake flux of 100 mmol g⁻¹DW h⁻¹ glucose. Red arrows represents the flux direction for the four reversible reactions

Furthermore, downstream fluxes to AcCOA, the sum of the downstream lipid and TCA cycle flux was of 74.01 mmol g⁻¹DW h⁻¹ compared to 86.14 g⁻¹DW h⁻¹ in literature. Although a similar total flux around the TCA cycle was obtained, with 73.93 mmol g⁻¹DW h⁻¹ at lipid branch and 0.08 at TCA branch, different results were reported in [18] with 81.21 mmol g⁻¹DW h⁻¹ at TCA cycle and 4.82 mmol g⁻¹DW h⁻¹ at lipid branch. This discrepancy may rely on a high lipid level (13.13% DW) in our cell culture compared to that in [18] (1% DW). Differences in culture conditions may be involved as well. Meanwhile, this difference of high lipid synthesis flux and low TCA cycle fluxes were also found in Wu’s work, where a ¹³C metabolic flux analysis was accompanied with flux balance analysis in Chlorella protothecoides [13]. Thus, from both the dynamic and steady state flux analysis, high lipid content in Chlorella is mainly due to low TCA split-flow. In our result, lipid was fully accumulated at 72 h, so the metabolic fluxes at 72 h was also extracted from the dynamic flux profile. Along with the glucose deleption, most of the metabolic fluxes decreased in different extent, some even decreased up to 78.87% (EMP pathway fluxes). However, the re-arrangement of flux distribution from starch and energy catabolism to lipid synthesis was obvious. At 48 h, starch and ATP were accumulating (V_PGM is 13.83 mmol g⁻¹DW h⁻¹) along with glucose assimilation. However, after glucose deleption at 72 h, starch was catablising (V_PGM is − 23.25 mmol g⁻¹DW h⁻¹) as a storage pool and ATP was also consuming (V _PPRiBP is − 0.01 mmol g⁻¹DW h⁻¹) to maintain other metabolism. Lipid synthesis flux (V_FASN), decreased a little bit from 73.93 to 63.89 mmol g⁻¹DW h⁻¹ (decreased 13.51%), which was less impacted by glucose deleption compared with other fluxes on EMP pathway. This may due to a constant feeding flux converting from starch and ATP catabolism. Therefore, the energy stored in starch and ATP seams to be converted to a more stable storage pool as lipid after glucose deleption.

We also looked at the major carbon distribution before glucose depletion. Model simulations show that the glucose uptake rate (V_HK) and the glycolytic fluxes went down to a very low level after day 2.6, we have thus analyzed their related flux ratios only before glucose depletion (< 2.6 days). First, we evaluated that 6% of the glucose flux contributes to biomass synthesis and growth (V_growth-to-V_PK ratio) (Fig. 5a), a value comparable to the literature with 3.9% [25]. Within the same range, 8% (8.323–8.325%) of the glucose flux feed lipid synthesis (V_FASN–V_Lipase to V_PDH ratio) (Fig. 5b). However, as a main product contributing to biomass, the lipid catabolism-to-biomass synthesis and growth ratio (V_FASN –V_Lipase –V_GPAT to V_growth) shows two successive constant values at around 60% increasing at 80% at mid-exponential growth phase (1.5 d) (Fig. 5e). Model simulations also suggest that around 1% of the glucose flux goes to starch synthesis (V_PGM to V_HK) (Fig. 5c), and that 15% to 7% of the glucose flux feed nucleotides synthesis (V_PPRiBP to V_G6PDH) (Fig. 5d). Concerning the PPP pathway activity, around 12% of the glucose uptake flux flow into the pentose phosphate pathway (V_G6PDH to V_HK) (Fig. 5f). Therefore, the dynamic metabolic flux analysis give our lights of the carbon distribution in Chlorella protothecoides.

Fig. 5

Metabolic flux ratios between different pathways. a Glucose contribution to biomass, b to lipids, c to starch and d to nucleotides; e lipid contribution to biomass; f PPP pathway activity

Parameter sensitivity showed biological significance on metabolic kinetics

A sensitivity analysis on model parameters showed flux maximum rate constants (V_max,i) to be more sensitive than affinity constants (K_m,i). For the final calibrated model 21 parameters, 15 maximum flux rates and 6 enzyme affinity constant (Fig. 6), out of 77 revealed greater sensitivity, defined as affecting the objective function of more than 10% when applying a − 70% to + 150% parameter value change around its optimized value.

Fig. 6

Sensitivity analysis on model parameters. Vertical axis value represents percentage change in the objective function for parameter change from − 70 to + 150% around the optimized value. Parameters not shown have percentage changes less than 10%

The most sensitive parameters are V_max,HK and V_max,GHMT, which are both at the entrance of the major carbon and nitrogen sources; V_max,GPI, V_max,PGM and V_max,PDH, which refer to fluxes at the intersection of glycolysis, starch and lipid metabolisms are also highly sensitive, while PPP pathway (V_max,TK) and TCA cycle parameters, showed a low sensitivity level. Since the major intersections from glycolysis including starch synthesis, PPP pathway which finally leads to nucleic acids, and TCA cycle which is related to protein synthesis, and the fatty acids synthesis which leads to lipids. The sensitive flux parameters on starch and lipid metabolism suggested that starch synthesis is the major competing pathway that affect lipid accumulation in C. protothecoides. It has also been reported in Chlamydomonas reinhardtii, when starch biosynthesis is blocked (sta6 mutant), the lipid content could be greatly boost, some can reach up to 30-fold [26]. Except for starch flux sensitivity, V_max,FASN, V_max,Lipase and V_max,GPAT are also sensitive. As V_max,FASN and V_max,Lipase are related to lipid synthesis and degradation, they are responsible for the balance of cellular lipid pool. Meanwhile, V_max,GPAT is in charge of providing glycerone-phosphate as the neutral lipid skeleton. The sensitivity of these fluxes gave us light on the genetic strategy for lipid yield promotion. Interestingly, there are two highly sensitive affinity constants (km_{,growth_lipid} and km_{,Lipase_lipid}), referring to the importance of lipid for cell biomass growth. Algae cell is a great platform accumulating lipids, some algae species could accumulate lipids up to 70% of their biomass. In C. protothecoides, the lipid content could reach 36% under heterotrophic condition, although in our work lipid content only reached 13% DCW, the sensitive of lipid affinity constant to growth suggest a huge potential to optimize the enzyme activity. Some reactions or pathways (i.e. their kinetic parameters) such as the maximum specific growth rate, PPP pathway and TCA cycle, showed a low sensitivity level, which suggest these are robust pathways. Final parameter values and the 95% confidence intervals for the sensitive parameters are shown in Table 4. They are all within ranges found in the BRENDA databank.

Conclusion

A model simulating Chlorella protothecoides cell metabolic behavior under heterotrophic condition and describing metabolic network flux kinetics and energetic states has been developed and calibrated. Simulation results show adequate fit with experimental data. Flux analysis is also in high agreement with literature data. A sustained high lipid synthesis metabolic activity was further confirmed from model simulations with higher lipid flux and lower TCA activity. The model was also used to analyze the dynamic distribution of the carbon source to the main carbon pathways, such as PPP pathway, starch synthesis, lipid synthesis and nucleotides synthesis. In addition, as the model included a high number of parameters, it described not only experimental data, but also most of the metabolic kinetics that showed statistical significance. It can thus be used as an in silico platform for characterizing the cell lines as well as to search for ‘‘optimal’’ culture strategy either by management through rational adjustment of the main nutrient concentrations that affect glucose and/or glycine concentration with time or by genetic manipulation of certain predicted critical enzymes. However, much work remains to be done: it would be of interest to add more metabolic reactions from extracellular multiple nutrients, like ions; and separating the lipids pools to more interest classes; get larger data sets including both extra- and intracellular experimental data, to test and validate the platform as a predictive tool.

Abbreviations

Abbreviations of metabolites names and parameter symbols were detailed in Tables 1, 2, 3 and 4.