A General Process-Based Model for Describing the Metabolic Shift in Microbial Cell Cultures

The metabolic shift between respiration and fermentation at high glucose concentration is a widespread phenomenon in microbial world, and it is relevant for the biotechnological exploitation of microbial cell factories, affecting the achievement of high-cell-densities in bioreactors. Starting from a model already developed for the yeast Saccharomyces cerevisiae, based on the System Dynamics approach, a general process-based model for two prokaryotic species of biotechnological interest, such as Escherichia coli and Bacillus subtilis, is proposed. The model is based on the main assumption that glycolytic intermediates act as central catabolic hub regulating the shift between respiratory and fermentative pathways. Furthermore, the description of a mixed fermentation with secondary by-products, characteristic of bacterial metabolism, is explicitly considered. The model also represents the inhibitory effect on growth and metabolism of self-produced toxic compounds relevant in assessing the late phases of high-cell density culture. Model simulations reproduced data from experiments reported in the literature with different strains of non-recombinant and recombinant E. coli and B. subtilis cultured in both batch and fed-batch reactors. The proposed model, based on simple biological assumptions, is able to describe the main dynamics of two microbial species of relevant biotechnological interest. It demonstrates that a reductionist System Dynamics approach to formulate simplified macro-kinetic models can provide a robust representation of cell growth and accumulation in the medium of fermentation by-products.


INTRODUCTION
Glucose is the main carbon and energy source for microbial metabolism. Glucose uptake supplies the glycolytic process producing different intermediates, with pyruvate representing a central catabolic hub, followed by the respiratory or fermentative pathway, depending on oxygen availability.
Respiration is able to maximize ATP production and consequently biomass yield. However, despite the fully aerobic conditions, in several microbial species when glucose concentration is high, the respiratory metabolism is replaced by a fermentative one, which produces partially oxidized products (Molenaar et al., 2009;Goel et al., 2012).
Such metabolic shift between two different ATP producing metabolisms, respiration and fermentation, is a widespread phenomenon in the biological world (Molenaar et al., 2009;Goel et al., 2012). In yeast it is known as Crabtree effect (De Deken, 1966) recognizing its similarity with the respiration/fermentation shift occurring in mammalian cells where it is commonly reported as Warburg effect (Warburg, 1956), and considered a hallmark of cancer (Hanahan and Weinberg, 2011). Moreover, some yeasts are recognized to be Crabtree-positive such as S. cerevisiae, others are Crabtreenegative (De Deken, 1966), the difference mainly relies on the extent of the glycolytic flux which, in turn, depends on the glucose uptake rate (Huberts et al., 2012). Recently, it has been shown how the overexpression of a single transcription factor (the ortholog of S. cerevisiae GAL4) in Komagataella phaffii results in a switch of the Crabtree phenotype from negative to positive with an increase in specific glucose uptake (Ata et al., 2018).
The fitness advantage associated to the metabolic shift and, more in general its significance, has been largely debated (Pfeiffer and Morley, 2014;Liberti and Locasale, 2016). Recently, a review and clarification of the process dynamics beyond this phenomenon has been proposed (de Alteriis et al., 2018).
The metabolic shift has been attributed to an "overflow metabolism, " caused by the saturation of the limited respiratory capacity of the cell, leading to an overflow reaction at pyruvate level, as first shown for the yeast Saccharomyces cerevisiae (Sonnleitner and Käppeli, 1986). It is now established that a complex interplay of molecular mechanisms is also responsible for the phenomenon, with the ascertained role of regulatory systems referred to as either carbon catabolite or glucose FIGURE 1 | Model diagram of microbial growth. Simplified cell metabolism with explicit representation of the major metabolic pathways. (1) Glucose uptake; (2) respiration; (3a) fermentation; (3b) acetate production by mixed fermentation; (4) acetate respiration; (5) secretion of inhibitory compounds; (6) inhibitory effects; (7) reserves accumulation; (8) cell death.
As known, the prokaryotic cell factories Escherichia coli, and to a lesser extent Bacillus subtilis, together with the eukaryotic

Feeding
All symbols are described in Tables 2-4.

Description Formula
Initial feed rate

Glycolysis products concentration
Medium volume in the reactor Other parameters are described in Tables 3, 4.
Frontiers in Microbiology | www.frontiersin.org unicellular fungus Saccharomyces cerevisiae are the prevalent microbial platforms for biotechnological applications (Öztürk et al., 2016;Sanchez-Garcia et al., 2016). The phenomenon of the metabolic shift with the consequent production of fermentative products has been widely described for these species, representing one of the problems which may limit the achievement of high cell densities and productivities (Riesenberg et al., 1991;Sandén et al., 2003), for both non-recombinant and recombinant microbial strains (Lee, 1996;Riesenberg and Guthke, 1999;Porro et al., 2005;Shiloach and Fass, 2005;Öztürk et al., 2016). In unrestricted growth conditions, the fermentative microbial metabolism generally leads to a main end-product and other by-products: E. coli and B. subtilis predominantly form acetate, but also lactate and propionate respectively, while in the case of S. cerevisiae a production of ethanol and, to lesser extent, acetate is observed.
In the case of E. coli, aerobic acetate production is very detrimental for growth and productivity, and for this reason it has been largely investigated (Wolfe, 2005;De Mey et al., 2007;Bernal et al., 2016). Several strategies have been proposed to avoid acetate production, from technological approaches (fedbatch cultures with controlled glucose supply, removal of acetate from culture medium, use of alternative carbon sources such as glycerol or mannose) to genetic approaches aimed to obtain strains with low propensity to acetate formation (Sandén et al., 2003;Shiloach and Fass, 2005;Eiteman and Altman, 2006).
The outstanding importance of E. coli in biotechnological processes supported the development of several mathematical models aimed to describe strain performance in different cultural conditions, optimizing their cell/product density and avoiding acetate overproduction. First attempts to use mechanistic models to simulate the kinetics of E. coli population growth followed different approaches, from simplified representations of batch cultures (Corman et al., 1986), based on classic biomass-resource model (Monod, 1949) to more detailed models of the main metabolic fluxes by an optimization method (Majewski and Domach, 1990;Ko et al., 1993Ko et al., , 1994. More recently, a processbased kinetic model first developed by Xu et al. (1999) was further improved to study the growth of E. coli W3110 strain in batch and fed-batch cultures, explicitly including the inhibitory effect of acetate accumulation on glucose and oxygen consumption (Lin et al., 2001;Neubauer et al., 2003).
The recent increasing studies on E. coli metabolism improved the understanding of the acetate production on one hand and, on the other hand, the co-assimilation of both acetate and glucose in sugar-limited conditions (Wolfe, 2005;Lara et al., 2008;Peebo et al., 2015;Basan et al., 2015;Bernal et al., 2016). Such new findings were integrated in new macro-kinetic models    (Anane et al., 2017;Retamal et al., 2018). In particular, Anane et al. (2017) developed a mechanistic model based on previous works (Xu et al., 1999;Lin et al., 2001;Neubauer et al., 2003) with two major improvements: (i) a mathematical formulation deriving a set of tractable and continuously differentiable equations leading to better computational performance and allowing the use of gradient-based optimization methods and (ii) the inclusion of a continuous process of production and re-assimilation of intracellular acetate even under non-overflow conditions, as recently highlighted in proteomic and systems biology studies (Valgepea et al., 2010;Basan et al., 2015;Peebo et al., 2015). The model proposed by Retamal et al. (2018), based on the overflow metabolism assumption (Sonnleitner and Käppeli, 1986), assumed that the critical glucose uptake rate responsible for the activation of the metabolic overflow is not constant, but decreases with increasing acetate concentrations. Following the consideration that the metabolic shift in S. cerevisiae is controlled by both limited respiratory capacity (Sonnleitner and Käppeli, 1986) and repression of respiration (Westergaard et al., 2007), our group developed a novel macro-kinetic model based on the System Dynamics approach (Forrester, 1961), capable to reproduce the growth of the budding yeast in both batch and high-cell-density cultures (Mazzoleni et al., 2015). The main assumption of this model was that the glycolytic intermediates represent the central metabolic hub regulating the shift between respiratory and fermentative pathways.
In this work, considering the similarity of the metabolic shift between S. cerevisiae and prokaryotic cells as determined by the level of the glycolytic intermediates, we extend the model by Mazzoleni et al. (2015) to simulate the growth of different strains of E. coli and B. subtilis cultured in batch and fed-batch bioreactors under aerobic conditions, with glucose as carbon and energy source.

MATERIALS AND METHODS
The model developed to simulate the growth behavior of a generic microbial cell cultured in a bioreactor is presented. Figure 1 shows a schematic diagram of the implemented processes, the limited number of which was achieved by a top-down approach, i.e., selecting the essential elements sufficient to reach a robust representation of the system behavior.
The resulting model is composed of a set of 7 ordinary differential equations representing glucose in the growth medium (G), glycolysis intermediates from glucose-6-phosphate to pyruvate (P), acetate produced by fermentation (A), cellular components produced by either fermentation or respiration (C M ), reserve compounds (R), growth-associated inhibitory byproducts (I), and dead cells (D).
Glucose (G) is provided according to the feeding conditions of the bioreactor described for the simulated experiment. Glucose is assimilated by microbial cells, and then converted into the different intermediate products of glycolysis, from glucose-6-phosphate to pyruvate (P). These are used for the construction of new cellular material (C M ), either through respiration or fermentation. In the case of fermentation, acetate (A) is the main end-product, which can also be used as carbon source for the respiratory pathway when glucose is limiting. The essential assumption of the model is the key role of the glycolytic products (P) in the regulation of the metabolic shift between respiration/fermentation and, in general, cell metabolism. Therefore, high levels of P are assumed to be responsible for (i) the activation of aerobic fermentation due to overflow metabolism, (ii) the repression of respiration ("glucose effect"), (iii) the accumulation of reserve materials (R), and (iv) the induction of mortality with accumulation of dead cells (D) (de Alteriis et al., 2018).
Moreover, considering that acetate is not the only endproduct of bacterial fermentations even in aerobic conditions (Park et al., 1992;Kim et al., 2015), in the presented model we also assumed that the allocation toward secondary byproducts proceeds in parallel with acetate. In particular, their production is assumed to increase with the acetate outflow from the cells.
The model also considers growth-associated inhibitory byproducts (I), as already described in Mazzoleni et al. (2015), whose production is related to anabolic pathways, hence it is expressed as a proportion of the respiration and fermentation fluxes. Both the inhibitors and acetate are assumed to separately exert a negative feedback on cell growth in a concentrationdependent way.
The model is formulated with the following mass-balance equations: The equations of the model are described in detail in Tables 1, 2, while fixed and calibrated parameters are described in Tables 3, 4. The mathematical equations were integrated using MATLAB R2018b (the MathWorks) with a variable order solver (ode15s).
The model calibration was performed by minimizing the sum of the squared errors (SSE) where n 1 , n 2 , n 3 are the number of samples per observed outputs, C Mi , G i , A i , are the values of the ith measured outputs and C * Mi , G * i , A * i , are the values of the ith outputs predicted by the model. The minimization was performed by using the fminsearch MATLAB routine which implements a Nelder-Mead simplex algorithm (Lagarias et al., 1998).
Furthermore, a sensitivity analysis was implemented to analyze the model behavior under parameters perturbations. Using a local sensitivity analysis (Morris, 1991;Norton, 2015), the following normalized sensitivity index was calculated by changing each parameter by ± 5% one-at-a time while keeping the rest constant: where, SSE i, is the Standardized elementary effect of the parameter p i with (± 5%) perturbation on the model outputs; X j (P) represents the simulation values of the state variables Microbial mass, Glucose and Acetate without any parameter perturbation; max(X j (p))−min(X j (p)) is the standardization factor referring to the values of the baseline simulation; k is the number of parameters.

RESULTS
Model simulations were compared to experiments of growth in bioreactor of two strains of E. coli among the most used in biotechnological applications, namely W3110 (Anane et al., 2017) and TG1 (Riesenberg et al., 1991;Korz et al., 1995). For E. coli TG1, also experiments of two recombinant strains were considered (Hellmuth et al., 1994;Rinas and Hoffmann, 2004). In the case of B. subtilis, model simulation was compared to the experiment by Huang and co-workers (Huang et al., 2004).
In the experiments selected for model simulations, microbial growth was carried out in bioreactors initially operating in batch mode and then fed with a glucose-based inlet stream (fed-batch). The possible metabolic shift occurred according to the value of the specific growth rate of the population, determined by the specific feeding rate (SFR) applied to the bioreactor during the fed-batch phase (Enfors and Häggström, 1998). Figure 2 shows the model simulation reproducing a twophase (exponential and constant feeding regimes) fed-batch culture of E. coli W3110, performed by Anane et al. (2017). The initial batch culture, characterized by a maximum specific growth rate value (µ MAX ) of 0.31 h −1 , presents the typical exponential growth behavior of an E. coli population growing on glucose and displaying a fermentative metabolism with acetate production. When glucose is completely depleted in the medium, fermentative metabolism is replaced by a respiratory one, with a short period of acetate consumption by respiration observed during the batch phase. Then, an exponential increasing feeding regime is activated (after 13 h from the beginning of the experiment) at a Specific Feeding Rate (SFR) value of 0.22 h −1 , equal to the population specific growth rate, µ. This value is high enough (75% µ MAX ) to switch back the population to a fermentative metabolism with accumulation of acetate up to 0.3 g l −1 . Three hours later (16 h from the beginning of the experiment), a constant feeding is applied to the bioreactor beginning at a SFR value of 0.11h −1 , so that the bacterial cells are allowed to display a respiratory metabolism with no accumulation of acetate. During this phase, however, three glucose pulses are performed inducing temporary acetate production (Figure 2, lower panel).   (Riesenberg et al., 1991;Korz et al., 1995). The first simulated experiment (Figure 3, left column) is characterized by an initial batch culture lasting 12 h when glucose in the medium is completely depleted, followed by a two-phase fed-batch, carried out at an exponential increasing feeding regime (up to 30 h from the beginning of the experiment) corresponding to a SFR value of 0.11 h −1 and then a constant feeding regime (from 30 h onward), starting at 0.11 h −1 . In the first phase of exponential feeding, the simulated microbial population follows the observed growth corresponding to the feeding regime, while in the second phase (constant feeding) the model properly describes the declining growth rate due to selfproduced inhibitory compounds. In turn, the reduced growth rate compared to the glucose feeding induces the metabolic switch reactivating acetate production. In the second simulated experiment (Figure 3, right column) the initial batch phase is followed by a single-phase fed-batch, carried out at a SFR value of 0.17 h −1 for the first 3 h and then reduced to 0.14 h −1 . The simulated dynamics of both experiments are very similar, showing an exponential growth of the cell population in the first fed-batch phase at the imposed µ value. In these conditions, no acetate is produced, since the set-point µ values (0.11 and 0.17 h −1 respectively) are below the threshold value for acetate production reported for the strain (Korz et al., 1995).
In brief, in both experiments of Figure 3, a reduction in the growth rate was observed near the end of the run. Such growth reduction appeared when microbial mass achieves a value around 100 g l −1 , and as explained above is modeled as ascribed to the production of growth-linked inhibitory compounds, different from acetate, consistent with previous findings in yeast (Mazzoleni et al., 2015). Figure 4 shows the simulations of two fed-batch cultures of recombinant strains of E. coli TG1 (Hellmuth et al., 1994;Rinas and Hoffmann, 2004) carried out at SFR of 0.13 and 0.12 h −1 , respectively. The behavior of the recombinant strains results very similar to the non-recombinant ones of Figure 3, with accumulation of acetate, as indicative of a fermentative metabolism, observed only in the initial batch phase.
The last simulation presented in Figure 5, shows the growth dynamics of B. subtilis cultured in fed-batch carried out at an exponential feeding regime of SFR = 0.12 h −1 . Also in this case, there is an adequate fit between the model and the experimental data. The culture achieves a very low density (about 17 g l −1 ) at the end of the simulated experiment, and it is characterized by a continuous, although very low, production of acetate. Figure 6 provides a summary of the model simulation performance, showing very good agreement between measured and simulated values of microbial mass for all the selected strains. The results of the sensitivity analysis are presented in   The generally low response of the model outcomes to the variation of the parameters shows that the model formulation is robust. In fact, a ± 5% change in each parameter induces significant variations only in a few cases. In particular, two parameters related to glucose uptake (v G and η G ) affected all strains (Figure 7), reflecting the relevance of this process in the model formulation. Differently, it is interesting to notice that each strain showed a specific sensitivity to different parameters. In particular, E. coli W3110 appeared to be sensitive to the fermentation process (v F and η FA ); E. coli TG1 showed higher responsiveness to the secretion rate of inhibitory compounds (ρ) and the respiration parameters (v R and η RP ); B. subtilis showed to specifically respond to the metabolic shift between respiration and fermentation (b 1 ) and the efficiency of biomass production by fermentation (η FP ).

DISCUSSION
The presented model is capable to reproduce the dynamic behavior of several Escherichia coli strains, as well as of Bacillus subtilis growing both in batch and fed-batch cultures on glucose as carbon as energy source. The highly significant agreement between experimental data and simulations obtained for different microbial species and strains demonstrates how the processbased System Dynamics approach, already followed in the case of the yeast Saccharomyces cerevisiae (Mazzoleni et al., 2015), can be successful to develop a general model of microbial growth in bioreactors, despite the extremely simplified representation of the main physiological functions limited to very few, but fundamental metabolic processes.
Indeed, as for yeasts, also in bacterial cells the dynamic levels of pyruvate, as end metabolic hub of the glycolytic process, play a central role in the control of metabolism. At high glucose concentrations, and consequently at high concentration of glycolytic intermediates, the differential rates of reactions along the fermentative pathway progressively trigger the activation of the overflow metabolism and the consequent repression of respiration (Mazzoleni et al., 2015;de Alteriis et al., 2018).
This essential assumption characterizing the presented model proved indispensable to predict the occurrence of the metabolic shift between respiration and fermentation in two species of prevalent biotechnological interest such as E. coli and B. subtilis, as shown by the very good agreement between simulations and experimental data (Figure 6). Therefore, acetate is produced in both batch and unrestricted fed-batch cultures, such as the first fed-batch phase of the experiment presented in Figure 2, when the imposed µ value is higher than the critical one for the examined strain. On the contrary, when glucose supply to the bioreactor is controlled, such as in the other simulated experiments of Figures 3, 4, an oxidative metabolism is ensured. Concerning the latter experimental setups of the E. coli TG1 strains, the supplied culture medium was the same between recombinant and non-recombinant strains (i.e., mineral medium supplemented with trace elements containing glucose as carbon and energy source). The obtainment of the desired products (beta-galactosidase and human growth factor) by the recombinant strains was achieved by shifting the temperature to 42 • C or by addition of 0.5 mM IPTG after 22 h fed-batch phase. In both cases, induction of the products did not affect the dynamics of growth, in terms of glucose consumption and acetate formation, as shown by the simulations.
Noteworthy, in the case of B. subtilis (Figure 5), even though the glucose supply was fairly low (SFR = 0.12 h −1 ), a limited production of acetate was observed during all the culture run, showing a less clear-cut metabolic shift for this microbial species. This is also reflected by the high sensitivity of this strain to the parameter related to the metabolic switch (b 1 ). A further assumption of the proposed model is related to glucose transporters which have been considered as constant within each microbial strain. Clearly, this is a strong simplification since it is known that transporters can be modulated according to glucose availability in the media and future specific studies could address this point more in depth.
Moreover, in the previous S. cerevisiae model presented by Mazzoleni et al. (2015), secondary products of fermentation were not considered, being ethanol largely predominant when budding yeast is in conditions of overflow metabolism. Differently, in the case of bacteria, acetate production is followed by a significant production of other partially oxidized products even in aerobic conditions. This is explicitly represented in the model by the description of a mixed fermentation which is assumed to increase with the acetate outflow from the cell. Due to lack of experimental data on such secondary by-products it was impossible to explicitly represent their accumulation in our simulations even though their production is accounted for. The model also considers that the produced acetate can be reassimilated through the respiratory metabolism at the same time as glucose if both substrates are available in the growth medium, but for simplicity, secondary pathways for acetate consumption were not considered (Anane et al., 2017).
The mechanistic models developed by Pham et al. (1998) to reproduce the growth of S. cerevisiae in aerobic fed-batch cultures and later extended to describe the fermentation dynamics of E. coli by Xu et al. (1999); Lin et al. (2001), and Neubauer et al. (2003) are probably the most studied microbial macrokinetic models. Very recently, Anane et al. (2017), refined the previous model formulation to improve the mathematical analyzability and include the latest knowledge on the acetate metabolism of E. coli. Our proposed model has a similar level of simplification compared to the abovementioned works, although it was designed to be more general in order to reproduce different microbial species and strains rather than being applied to a single one. Moreover, our model is the only one considering the inhibitory effect on growth and metabolism of self-produced toxic compounds different from fermentation products which is relevant to reproduce the late phases of high cell-density cultures when this phenomenon becomes significant in limiting the growth rate. The phenomenon of self-toxicity regulating cell proliferation is evident only in prolonged fed-batch cultures, and it was clearly pointed out in the case of both wild-type and auxotrophic yeast strains cultured in fed-batch reactors (Mazzoleni et al., 2015). However, in this paper the cell densities achieved in E. coli experiments are sufficiently far from the theoretical value of maximum cell density for bacterial cells (200 g d.w. l −1 ) (Lee, 1996). For this reason, the growth decline due to self-toxicity was observed only in the simulated experiment presented in Figure 3, reproducing the experiments of Riesenberg et al. (1991), where the cell density achieved at the end of the run was higher than 100 g d.w. l −1 . In this case, the final growth decline was clearly visible (Figure 3, left panels) and the model was perfectly capable to reproduce such behavior.
Following a completely different approach, metabolic flux analysis models specifically focused on the central carbon metabolism of E. coli, show a detailed description of the metabolic pathways (e.g., Chassagnole et al., 2002;Lemuth et al., 2008). Recently, Millard et al. (2017) developed a detailed kinetic model linking the internal metabolism to the environment and cell proliferation through the description of the dynamics of 62 metabolites, and 68 reactions divided into 3 compartments (environment, periplasm and cytoplasm). This model has been validated using 226 experiments from different sources, allowing the authors to conclude that the self-regulating capabilities of the E. coli central metabolism are far more important than expected, also undermining the relevance of gene regulation to explain these dynamics. In typical metabolic flux analyses of Systems Biology, all the measured processes are considered and eventually reduced by selection techniques based on their relevance (bottom-up approach). On this point, Chassagnole et al. (2002) declare: "Because the many biochemical details of the metabolic networks appear overwhelming at first sight, there is a demand for decreasing the enormous complexity of the problem." As an example, Erdrich et al. (2015) used a combination of "pruning" and "compression" procedures, dramatically reducing the number of reactions from 2384 to 88.
Instead, our modeling procedure directly aims at the identification of a minimal number of processes sufficient to simulate the emergent properties of a complex system (topdown approach) and it is based on logical reasoning on existing knowledge of the system. Then, in this work only those variables relevant for the growth of microbial populations on glucose and their metabolic shift were mathematically described, whereas the many secondary pathways, also directly or indirectly affecting the selected variables, were not taken into account. In summary, our approach is based on the principle of parsimony and has the advantage of keeping the mathematical formulation simple and robust with a reduced number of parameters.

CONCLUSION
In conclusion, the results demonstrate how a reductionist System Dynamics approach can be used to formulate simplified macrokinetic models, still capable to accurately capture the dynamics of biomass growth, glucose consumption and accumulation of fermentation by-products in the medium. The conceptual base of the model, similar to that already proposed for S. cerevisiae (Mazzoleni et al., 2015), suggests a unifying theoretical view for all microbial species, with a key role of the metabolic shift phenomenon despite differences in specific molecular mechanisms.
Moreover, the robustness of the results supports a future potential application of the model as a tool for optimization and control of microbial fermentation processes of the main species of biotechnological importance.