Dynamic Metabolic Flux Analysis Demonstrated on Cultures Where the Limiting Substrate Is Changed from Carbon to Nitrogen and Vice Versa

The main requirement for metabolic flux analysis (MFA) is that the cells are in a pseudo-steady state, that there is no accumulation or depletion of intracellular metabolites. In the past, the applications of MFA were limited to the analysis of continuous cultures. This contribution introduces the concept of dynamic MFA and extends MFA so that it is applicable to transient cultures. Time series of concentration measurements are transformed into flux values. This transformation involves differentiation, which typically increases the noisiness of the data. Therefore, a noise-reducing step is needed. In this work, polynomial smoothing was used. As a test case, dynamic MFA is applied on Escherichia coli cultivations shifting from carbon limitation to nitrogen limitation and vice versa. After switching the limiting substrate from N to C, a lag phase was observed accompanied with an increase in maintenance energy requirement. This lag phase did not occur in the C- to N-limitation case.


Introduction
Material balances of reactors describe how certain compounds go in the reactor and are transformed into other compounds. Metabolic models consider the cell as the reactor and, additionally to the exchange rates, also include information on the different steps (reactions) on how the input compounds are transformed to the output compounds. This should allow to gain insight in the internal cellular fluxes. Typically the number of equations (material balances of the different compounds) obtained is less than the number of unknowns (exchange rates and reaction rates) one wants to calculate.
Techniques for solving metabolic models can be divided into two main groups: data-driven techniques and knowledge/assumption-driven techniques.
The prime example of a data-driven technique is metabolic flux analysis (MFA) where intracellular fluxes are calculated based on measured exchange rates. The mathematical technique used was initially developed for black box modeling [1,2], but can easily be applied to stoichiometric models [3,4]. Exchange rates are not always sufficient to solve a stoichiometric model. In such cases, the option exists to measure intracellular fluxes via 13 C enrichment analysis [5].
The counterpart of MFA is flux balance analysis (FBA) which uses linear optimisation to solve the stoichiometric model [6] and does not directly use measurements (although measurements can be included [7]). Whereas FBA typically yields one optimal solution, the calculation of elementary modes or extreme pathways gives a full view of the metabolic capabilities of an organism [8,9]. An overview of the different extensions applied to MFA and FBA is given in [10].
Both MFA and FBA have the fundamental assumption that the cell is in a pseudo-steady state, thus that there is no accumulation or depletion of intracellular metabolites. This implies that metabolic models can only be used for cells in pseudo-steady state. Models that fully describe the cellular dynamics exist, but are complex and cumbersome to implement and/or very limited in scope [11,12].
However, some examples of FBA applied in nonstationary cases can be found. Varma and Palsson [13] used flux balance analysis to describe fed batches by iteratively solving the model for maximal biomass production. At each time point, the available substrate is calculated from the results of the FBA in the previous step. This way the time profiles of cell density, glucose, and by-products could be quantitatively predicted. This concept of dynamic FBA (DFBA) is further developed by Mahadevan et al. [14]. They formalise the methodology of Varma and Palsson [13] and name it static optimisation-based DFBA. They also introduce a new methodology, called dynamic optimisation-based DFBA, in which the optimisation is done over the entire time period of interest to obtain time profiles of fluxes. The methodology of dynamic optimisation-based DFBA is combined with the concept of minimisation of metabolic adjustment (MOMA) to model the myocardial energy metabolism [15]. Minimisation of metabolic adjustment (MOMA) is an alternative goal function to better predict the behaviour of mutated strains [16]. Using this technique, it was shown that the cellular metabolism does not always remain optimal during transient perturbations [16]. Lee et al. [17] introduced integrated dynamics FBA (idFBA) that dynamically simulates cellular phenotypes arising from integrated networks. These networks combine metabolic models with genetic, regulatory, and intracellular signalling information.
In this contribution, not FBA but MFA is extended so that it can be used for modelling cultures during transient conditions. This is reasonable as intracellular pseudo-steady state can also be assumed under certain dynamic conditions, because of the relatively small time constants of cellular processes, for example, mass action and metabolic adaption to novel conditions, in comparison with processes affecting the observed environmental conditions. Indeed, perturbation experiments performed in E. coli [11,[18][19][20] showed that intracellular metabolite pools reach a new pseudo-steady state after 20 seconds. Hence, an intracellular pseudo-steady state can be assumed during the transient.
Pseudo-steady state conditions are typically encountered in chemostat operated cultures [21,22]. Chemostats are mostly run under the concept that only one substrate is limiting. However, conditions do exist where multiple substrates can be limiting at the same time. Practical applications of such dual limitations have been described by Zinn et al. [23]. Egli [24] has described concentration zones in which both glucose and nitrogen are limiting. In this contribution, however, the limiting compound in the feed medium of steady state cultures is switched from glucose to ammonia and vice versa. No abrupt changes are thus sensed by the cells (the reactor broth is only gradually replaced), an ideal test case for applying dynamic FBA. The exchange fluxes of metabolites are determined based on their concentration in the reactor broth and are subsequently used to solve an overdetermined metabolic model, resulting in the determination of the intracellular flux values during the transient period. As exchange fluxes are based on the derivative of the concentration data, the latter ones are first smoothed to avoid amplification of the noise in the exchange fluxes.  [25].

Cultivations.
Culture conditions and inoculation procedure are described in [25].
The experiment in which the limiting substrate was changed from glucose to ammonia was conducted at a dilution rate of 0.155 h −1 . The experiment in which the limitation was changed from nitrogen to carbon was conducted at a dilution rate of 0.142 h −1 .
For each experiment, steady state samples were taken: optical density (OD), cell dry weight (CDW), and HPLC analyses. These samples were taken after waiting at least five residence times after the batch phase. To ensure that the perturbations caused by sampling did not influence steady state, switching the limiting nutrient was only done after another five residence times. The transition between the two limitations was frequently sampled for OD and metabolite data (HPLC). Sampling of the second steady state was performed again after five residence times. Furthermore, during the culture, dissolved oxygen, pH, stirrer rate, temperature, airflow and oxygen and carbon dioxide content of the offgas were continuously logged. Two balances, one under the influent and one under the effluent barrel, were coupled to the system to allow precise measurement of the dilution rate.

Measurements.
The sampling and measurement methodologies used are described in [25].

Data Analysis.
Due to the amount of data that had to be combined for each experiment, a custom-made program was written in Python using the SciPy scientific library [26,27] to process the data and apply the algorithm as explained in Section 3 (polynomial fitting, extracting the derivatives, and calculating the fluxes).
Each point in the time series of the transient experiments is only measured once. This is deemed sufficient, as nonsystematical variations are smoothed out when approximating the data with polynomials (as described in Section 3.2). In this polynomial approximation, error propagation calculations are not performed. For the MFA calculations, however, a covariance matrix is needed. This matrix is used to rescale the measurements so that fluxes that are known to be inaccurate do not weight as much during the data reconciliation as other, more trustworthy fluxes [1]. Therefore, a "mean" covariance matrix was calculated from data presented in [25] and used for the MFA part of this analysis. The experiments described therein are similar to those performed for this publication (also using the same equipment) and it is assumed that the utilised covariance matrix is a good approximation of the real one. The standard deviation flags in the figures are based on this covariance matrix.
The MFA modelling was performed as described in [4], although error analysis on the reconciliated fluxes was not done. The model used is based on the one in [4] and is the same metabolic model as described in [25]. It contains 136 reactions (Table 1) and 150 metabolites (Table 2), of which 12 are exchangeable. There are 142 independent equations and 136 + 12 = 148 unknowns. Thus at least 6 measurements have to be performed in order to fully solve the model. Ten exchange metabolites were measured and used for solving the stoichiometric models (CO 2 , O 2 , NH 3 , PiOH, acetate, lactate, pyruvate, succinate, glucose, and biomass) giving 4 redundant measurements that can be balanced. The methods for balancing measurements were initially developed for black box models [1,2,28], but can easily be applied on stoichiometric models [3,4].
Although it is known that the biomass composition is different under C-and N-limitation [29,30], the same model was used in both cases thus implying the same biomass composition for cells grown under N-and C-limitation. Limited testing showed that the influence of a different protein content of the biomass on the flux fluctuations was minimal. This approach of Provost and Bastin [31] can be formalised and extended. Instead of limiting oneself to cases where a straight line is sufficient to adequately capture the dynamics of the measured metabolites during the transients, one can take the derivative in each point, without approximating the whole interval with a straight line. This derivative is, by definition, a flux.

Dynamic Metabolic Flux Analysis
The differential equation that governs the change of a component in the reactor broth can be written as where C is the concentration of the component, C in the concentration of the component in the influent, r p the production rate, and D the dilution rate. From (1) the net reaction rate r p can be calculated using the following approximation for the derivative in a point n of the time series: Eventually, more data around the working point might be used to perform the linear regression. However, because of the noisiness of the data, this approach does not give satisfying results. Therefore, approximation techniques have to be used. In this work, a polynomial approximation was chosen.

Polynomial
Fitting. The polynomial fitting algorithms found in SciPy [26] were used to approximate the data. The whole time series were covered by different polynomials. To this end, a moving window approach was used ( Figure 1). The points inside the outer window, W 1 , are used for fitting the polynomial. The derivatives of the polynomial in the points of interest within in the inner window W 2 are used for dynamic MFA. Making the inner window smaller than the outer window ensures a smooth transition between the different polynomials. However, if the inner window is too small compared to the outer window, the noise on the data is not filtered sufficiently, yielding nonrealistic fluctuations in the derivatives. Some of the available data series were impossible to smooth with polynomials. They behaved like a typical logistic curve. Hence, the logistic curve was added to the set of functions that could be used to approximate the data: The parameters W 1 and W 2 and the maximal degree of the polynomials (or logistic curve) were manually selected by visually inspection of the generated curves. The actual degree of the polynomials could vary for each polynomial fitted (thus the curve for a time series of a certain metabolite consists of polynomials of different degrees chained together): for each W 1 window the polynomial giving the least amount of outliers was selected. Data points were identified as outliers when the difference between the measured value and  Journal of Biomedicine and Biotechnology    Time series Figure 1: A moving window is run through a time series of data. W 1 is the polynomial fitting window, W 2 is the interpolation window.  the estimated value was larger than 5% of the range of the response variable. Table 3 summarized the fitting parameters used for each metabolite. The polynomial or logistic functions yield a smoothed time profile of the concentrations and enable to calculate the derivative of the concentrations to time. However, when calculating this derivative directly by symbolically deriving the function and filling in the time for the point of interest, one sometimes obtains unrealistically high values (D 1 in Figure 2). Therefore, the derivative was approximated by taking the slope of the point of interest and the previous point on the polynomial (D 2 in Figure 2).
The derivative can subsequently be used in (1) to calculate the cellular production (or consumption) fluxes: Once r p has been calculated, classical MFA techniques (as described in [4]) can be applied at each time point. The results of those models, evaluated at the different time points, can then be combined to obtain a time profile of internal fluxes.

Growth Rate Calculation.
During the transient period, the assumption that the growth rate is equal to the dilution rate does not hold. However, (4) can explicitly be rewritten for the following biomass: from which the growth rate can be calculated as

Polynomial Fitting.
The polynomials fitted through the data (Figures 3 and 4) are used to obtain the concentration, C, and the derivative of the concentration to the time, dC/dt, in each point. From those values, the conversion rate can be calculated (4). Biomass concentrations were calculated from OD values, after generating a calibration curve between cell dry weight (CDW) and OD values. This calibration curve was generated at the end of each experiment. As the curves for both experiments were very close to each other, no distinct calibration was used for the C-limited and N-limited conditions. The steady state concentration of NH 3 during Nlimitation in the C-to N-limitation experiment is around 0.2 g/L (data are not shown on Figure 3, as this value falls out the time frame shown), while it is 0.1 g/L in the N to C-limitation case. These values are high compared to what is found in the literature for the residual nitrogen concentration in similar strains [29] and the nitrogen measurement kit used is probably not suitable for accurately measuring very low NH 3 concentrations.

Growth Rate and Biomass Production.
At time zero, the influent medium was simply switched. No cells were extracted and resuspended in the other medium, as done in [32]. Thus at the switch of the limiting compound, there is an abundance of both, as the medium limited in carbon is saturated with nitrogen while the medium limited in nitrogen is saturated with glucose. In the C to N-limitation experiment, it seems as if N-limitation is reached not much later than when carbon abundance starts (Figure 3). In the Journal of Biomedicine and Biotechnology N to C-limitation case, nitrogen abundance is reached before all the glucose is depleted: after 5 hours, both glucose and nitrogen are not limiting. This can be seen in Figure 4, in the glucose profile. The glucose concentration becomes only zero after 12 hours while the nitrogen concentration already reaches steady state after 5 hours. At the medium switch, as there is excess nitrogen, the growth rate could increase, but apparently E. coli needs some time to adapt to the new conditions. The biomass flux even decreases as nitrogen is added and becomes almost zero (top of Figure 5). Only seven hours after the medium switch, the growth rate increases sharply to attain almost 0.6 h −1 (right part of Figure 5). This increased growth rate coincides with the uptake of all glucose left in the broth (Figure 4). Figure 5 depicts the different growth rates during the transient conditions occurring when switching the limiting substrate. The same trends as in the biomass flux (top of Figure 5) can be observed. For the C-to N-limitation case, a sharp increase in growth rate is observed in the beginning, while for the N-to C-limitation the growth rate decreases to almost zero before increasing to 0.6 h −1 , not far from the maximal growth rate of 0.7 h −1 . In the carbon to nitrogenlimitation case, the maximal possible growth rate is never achieved and after a sharp increase in the beginning to 0.35 h −1 the growth rate decreases.
The sharp increase in growth rate of the cells at the beginning of the medium switch in the C-to N-limitation experiment (left part of Figure 5) is not accompanied by an increased flux through the Krebs cycle ( Figure 8, the first hour after the steady state point, labeled with a triangle). An increase in biomass flux does not imply that intracellular fluxes relative to the biomass flux should increase as well. On the contrary, in the beginning the cells are more energyefficient than during the carbon-limited steady state (left subfigure of Figure 7 shows that the ATP hydrolysis flux relative to the biomass flux is lower in the beginning implying that the maintenance energy requirement is lower). However, very quickly the flux through the glycolysis increases.
During the switch of limiting substrate from carbon to nitrogen, less oxygen is consumed and less CO 2 is produced than during steady state ( Figure 6). In the beginning of the medium switch, a sharp decline in O 2 consumption and CO 2 production is observed. This correlates with the increased growth rate and biomass production rate ( Figure 5). For the experiment where nitrogen limitation is switched to carbon limitation, two parts can be observed: in the first period (until five hours), the O 2 consumption rate and CO 2 production rate increase. Then they suddenly decrease (around 7 hours), to steadily increase again until they reach the steady state level. In both experiments, the steady state value of oxygen consumption and carbon dioxide production is around 16 g/mol Biomass/h. This is in accordance with values obtained with a similar strain [4].
The ATP hydrolysis reaction (Figure 7) lumps all the maintenance requirements and futile cycles (as including all the futile cycles in the metabolic model would yield parallel pathways and the system of equations would not be solvable without measuring intracellular fluxes). It can be seen that during the lag phase of the experiment switching nitrogen to carbon-limitation (right part of Figure 7), the ATP maintenance requirement is around 0.6 mol ATP/mol Biom/h, which is high compared to the values found in a similar E. coli [4]. It seems even more unusually high because there is minimal growth during this period. This high ATP requirement in the beginning of the lag phase correlates with a high uptake of glucose (PTS reaction in Figure 9) and an increased flux in the TCA cycle ( Figure 9). The somewhat higher steady state ATP flux at the end of the C-to Nlimitation (left part of Figure 7) may be due to multiple physiological steady states [33].

Central Carbon Metabolism.
During the transient in the C-to N-limitation experiment, the fluxes through the glycolysis increase while those through the Krebs cycle decrease. Also an increase in acetate excretion is observed (Figure 8). The opposite is observed in the N-to Climitation experiment: fluxes through the glycolysis and acetate production decrease while the fluxes through the Krebs cycle increase (Figure 9). The decrease is not linear: at 15 hours the fluxes through the glycolysis reach a minimum to increase again until the steady state values are obtained after the subsequent 5 hours. In the Krebs cycle, the fluxes gradually increase the first five hours (coinciding with the lag phase observed in the biomass production, Figure 5) after which a sharp decline is observed, immediately followed by a gradual increase, peaking at around 15 hours, and then decrease to steady state level. The succinate flux sharply decreases in the beginning of the switch to gradually increase during the lag phase (5 hours) at which point it reaches a value that does not change anymore. The PEPCB flux, which is responsible for replenishing to the Krebs cycle, follows (or predates) this succinate excretion rate. The excretion of pyruvate decreases during the lag phase and then sharply increases again to reach a lower steady state than under N limitation. The pentose phosphate pathway is entirely driven by the demands in biomass components, and as each flux is expressed relatively to the biomass, to remove the influence of biomass on the flux values, the fluxes of the pentose phosphate pathway do not change, which results in seemingly huge error bars.

Fluxes.
When going from an nitrogen-limited medium to a carbon-limited one, a lag phase was observed ( Figure 5). A hypothesis could be that when more nitrogen becomes available, the cell starts to break down the accumulated glycogen which is redirected to energy metabolism for protein production. As the biomass composition is assumed to be constant in the model, this increased respiration is measured in the off-gas (O 2 and CO 2 fluxes, Figure 6) and propagated in the model to the ATPHY flux ( Figure 7). However, the question remains why the cell would degrade glycogen while there is still plenty of carbon available in the environment. Furthermore, the carbon balance data (not shown) do not reflect an increase in carbon output compared to the carbon input. Another explanation is that ammonia abundance for cells previously grown under ammonia limiting conditions causes a futile cycle of ammonia import and export. Such cycles have been described previously in micro-organisms, although not under the experimental conditions described in this work [34]. Normally bacterial cells are highly resistant to ammonium, and the negative effects of high NH + 4 concentrations are due to an enhanced osmolarity or increased ionic strength of the medium and not caused by ammonium itself [35]. But when changing the medium from nitrogenlimitation to carbon-limitation, the environment of the bacteria switches from low NH + 4 concentrations to high ones, and it could be that E. coli needs a certain adaptation period. The scavenging active during N-limitation could still be active in the beginning of the transition, resulting in a high influx of ammonia. To counter this, active efflux is needed, draining ATP and explaining the high energetic cost in the beginning of the addition of ammonium (right part of Figure 7).
Such a mechanism has been observed in plants, where NH + 4 toxicity is due to the inability of root cells to limit the influx of ammonium. The high cytosolic NH + 4 concentration activates the high-capacity and energy demanding ammonium efflux systems. This ammonium efflux can constitute as much as 80% of primary influx, resulting in a futile cycle of nitrogen across the plasma membrane of root cells. This futile cycle imposes a high energetic cost on the cell that is independent of N metabolism and is accompanied by a decline in growth. Plants that are resistant to high ammonium concentrations (e.g., Oriza sativa) limit the influx of ammonium by lowering the polarisation of the cellular membrane. This lowering of membrane polarisation is not observed in Hordeum vulgare, sensitive to ammonium [36].   The lack of growth resulting from adding ammonium to cells that are adapted to ammonium limitation has previously only been observed in cells lacking glnE (giving the gene product ATase, responsible for activating/inactivating glutamine synthetase) while the wild-type E. coli did not show an impaired growth when shifted from ammonium limitation to excess [37,38]. Furthermore, the lack of growth in glnE mutants could be alleviated in cells constitutively expressing only low levels of glutamine synthetase [37], suggesting that the toxicity of ammonium in this case is not due to NH + 4 itself, but maybe to accumulation of glutamine and/or glutamate. Normally, when shifting from ammonium limitation to carbon-limitation, the amount of active glutamine synthetase decreases instantly [39]. Possibly due to the lack of ATase, the enzymatic glutamine synthetase regulation is lost, and only genetic regulation is present. If this would be the case for the strain used in this study, it could explain the long lag phase after the medium switch. The cells survive until enough glutamine synthetase is degraded before they are able to grow again.
Thus it seems as E. coli MG1655, the strain used in this study, has some anomalies in its nitrogen metabolism. It is known that multiple strains are denominated under this name [40]. The strain used in this work originated from the Coli Genetic Stock Center (CGSC) (personal communications from the Netherlands Culture Collection of Bacteria). The NCCB acquired this strain from CGSC prior to 1991. At that time, the strain was sent out as being f nr − . Later it was found that the strain was a mixture of f nr + and f nr − strains (personal communication from Mary Berlyn). For the experiments in this work, it turned out that the f nr − strain was used. More research is needed to find out why this strain does not grow well when NH + 4 is added to a nitrogenlimiting environment.
The lag phase found in the N-limitation to C-limitation experiment stands in sharp contrast with the experiment in which the carbon-limiting medium is replaced with the nitrogen-limiting one and thus carbon is added in the beginning of the switch while nitrogen is not yet limiting. The growth rate increases instantly (top of Figure 5). Death and Ferenci [41] demonstrated that during carbon-limitation, sugar regulons are upregulated by internally synthesised sugars. Lequeux et al. [4] observed that under glucose limitation ptsG is upregulated. The results of the experiments described here show that this up-regulation permits the cells to almost instantly increase their growth rate when carbon is added to a carbon-limited culture. The cells grow as fast as possible, until all nitrogen is consumed and thus the glucose concentration in the reactor broth increases around the same time as the nitrogen decreases ( Figure 3). Hence the peak of succinate (and the smaller peak of pyruvate) occurring around 5 hours after the medium switch ( Figure 8) could be due to the depletion of ammonium ( Figure 3). The high flux through the glycolysis is not stopped quickly enough and the carbon is excreted as pyruvate or diverted to PEPCB to be excreted as succinate.

Dynamic MFA.
The extension of metabolic flux analysis to allow to predict also intracellular fluxes for cells growing under nonsteady state conditions was successfully applied to transient cultures. As raw measurement data are too noisy, smoothing had to be performed prior to the calculation of the exchange rates, needed for solving the metabolic model. An improvement for this smoothing process would be to incorporate more than two points in the calculation of the derivative, to avoid steep derivatives (and thus unrealistic fluxes) caused by the quirks of the higher-order polynomials.
A possible improvement to the model is to allow fluctuations of the biomass composition. However, this should be accompanied with the appropriate measurements such as protein, RNA, DNA, lipid... content of the biomass. The biomass reaction can then be removed from the model and instead the different biomass compounds can be considered as excreted by the cell.

Conclusion
This work introduced the concept of dynamic metabolic flux analysis. We demonstrated how to transform extracellular measurement data from dynamic experiments to flux values. It is not always clear in transient experiments whether the decrease in concentration of the extracellular metabolites is due to the cells stopping production and the remaining product diluting out or whether the cells are actively taking up the metabolites. Transforming the extracellular time series of concentrations to flux values can give the answer to that question. Furthermore, those flux values can help to get insight into the intracellular reactions.
The transformation of time series of concentration measurements to flux values is based on differentiation (in the mathematical sense: finding the derivative) of those time series. Differentiation typically amplifies the noise on the data; therefore a noise-reducing step is needed prior to the differentiation. In this work, this was done by polynomial filtering. Selecting the right parameters for this filtering was done manually and produced satisfying results for the case presented. However, a statistically more sound approach may be desirable.
We subsequently applied dynamic metabolic flux analysis on E. coli cultures in which the limiting compound in the medium was changed from nitrogen to glucose or from glucose to nitrogen.
A lag phase occurred when cells adapted to a nitrogenlimiting environment were supplied with nitrogen. During this lag phase of several hours, the ATP demand was high. No clear physiological reason could be given. The only similar  case described in literature for E. coli is a glnE knock out, in which case the toxicity was not due to ammonium itself, but to the accumulation of glutamine/glutamate because the down-regulation of glutamine synthetase is not operative (glutamine synthetase is down-regulated by the gene product of glnE). A possible explanation for the long lag phase observed in the N-to-C-limitation experiment could be that sufficient glutamine synthetase has to be degraded so that the accumulation of glutamine/glutamate does not impair growth anymore.
No lag phase occurred when supplying cells adapted to carbon limitation, with carbon. Instead the growth rate increased immediately. However, the highest growth rates were attained in the N to C-limitation experiment and not in the C-to N-limitation one.