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

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. Electronic supplementary material The online version of this article (10.1186/s12934-019-1163-4) contains supplementary material, which is available to authorized users.


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

Open Access
Microbial Cell Factories  [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.

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 −1 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 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, Meta-Cyc, 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.
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.
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. (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 [M i ] × V growth is the cell dilution term and c i × 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 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 Fig. 1 The model metabolic network for heterotrophic Chlorella protothecoides 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).
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. brend a-enzym es.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.
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 subroutine "nlparci.m" (Table 4). It is clear there is no unique solution for parameter values in such an underdetermined system.

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.

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 −1 day −1 and 0.05 mmol gDW −1 day −1 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.
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 −1 DW h −1 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 −1 DW h −1 (leading to glycolysis), G6PDH of 12.15 (leading to PPP pathway) and PGM of 13.83 (leading to starch), compared to 49.   [18]. Furthermore, downstream fluxes to AcCOA, the sum of the downstream lipid and TCA cycle flux was of 74.01 mmol g −1 DW h −1 compared to 86.14 g −1 DW h −1 in literature. Although a similar total flux around the TCA cycle was obtained, with 73.93 mmol g −1 DW h −1 at lipid branch and 0.08 at TCA branch, different results were reported in [18] with 81.21 mmol g −1 DW h −1 at TCA cycle and 4.82 mmol g −1 DW h −1 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

Table 4 Parameter values and 95% confidence intervals of the highly sensitive parameters
The flux-rates' units are in mmol gDW −1 day −1 , except for the maximum specific growth rate (V max, growth ) which is in day −1 . Enzymes affinity constants' units are in mmol gDW −1 , except for HK and GHMT, which are in mmol L −1 flux and low TCA cycle fluxes were also found in Wu's work, where a 13 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 −1 DW h −1 ) along with glucose assimilation. However, after glucose deleption at 72 h, starch was catablising (V PGM is − 23.25 mmol g −1 DW h −1 ) as a storage pool and ATP was also consuming (V PPRiBP is − 0.01 mmol g −1 DW h −1 ) to maintain other metabolism. Lipid synthesis flux (V FASN ), decreased a little bit from 73.93 to 63.89 mmol g −1 DW h −1 (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.

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.
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.