Complex Nonlinear Behavior in Metabolic Processes: Global Bifurcation Analysis of Escherichia Coli Growth on Multiple Substrates

The nonlinear behavior of metabolic systems can arise from at least two different sources. One comes from the nonlinear kinetics of chemical reactions in metabolism and the other from nonlinearity associated with regulatory processes. Consequently, organisms at a constant growth rate (as experienced in a chemostat) could display multiple metabolic states or display complex oscillatory behavior both with potentially serious implications to process operation. This paper explores the nonlinear behavior of a metabolic model of Escherichia coli growth on mixed substrates with sufficient detail to include regulatory features through the cybernetic postulate that metabolic regulation is the consequence of a dynamic objective function ensuring the organism's survival. The chief source of nonlinearity arises from the optimal formulation with the metabolic state determined by a convex combination of reactions contributing to the objective function. The model for anaerobic growth of E. coli was previously examined for multiple steady states in a chemostat fed by a mixture of glucose and pyruvate substrates under very specific conditions and experimentally verified. In this article, we explore the foregoing model for nonlinear behavior over the full range of parameters, γ (the fractional concentration of glucose in the feed mixture) and D (the dilution rate). The observed multiplicity is in the cybernetic variables combining elementary modes. The results show steady-state multiplicity up to seven. No Hopf bifurcation was encountered, however. Bifurcation analysis of cybernetic models is complicated by the non-differentiability of the cybernetic variables for enzyme activities. A methodology is adopted here to overcome this problem, which is applicable to more complicated metabolic networks.


Introduction
Historically, microorganisms have been utilized for the production of valuable products in our daily life, e.g., bread, vinegar, wine and beer.With the advent of recombinant DNA technology several decades ago, it is common practice to make genetic modifications to microbes for the industrial production of food, energy, medicine and other valuable products.Towards ensuring the economic competitiveness of those commercial processes, maximizing productivity is one of the goals to achieve.
It is a challenge to manipulate cellular metabolism due to its complexity.Metabolic systems often exhibit intricate nonlinear behaviors, such as steady-state multiplicity and dynamic oscillations.It is necessary to understand what triggers this breadth of behavior and to predict when and under what conditions they would occur.Such a study is also practically important, as nonlinear behavior should be avoided if it prevents stable operations [1] or may be induced if it can lead to higher productivity [2].
A basic source of nonlinearity in a metabolic system is the intrinsic kinetics of biochemical reactions.More importantly, however, nonlinear metabolic behavior becomes much more complex and diverse due to regulation that dynamically drives individual reactions in response to environmental changes.Dramatic shift between multitudes of metabolic pathways often arises in a dynamic environment as a consequence of metabolic regulation.For the nonlinear analysis of metabolic systems, therefore, it is essential to employ metabolic models that are able to appropriately account for dynamic regulation.Various modeling ideas have been developed for the analysis of metabolic systems, including metabolic pathway analysis [3,4], constraint-based approaches [5,6], kinetic models [7] and the cybernetic approaches [8].In the discussion of the conceptual distinctions and commonalities among different modeling frameworks, Song et al. [9] highlighted the essential need for dynamic modeling frameworks in a wide range of applications, such as the study of complex nonlinear behavior of metabolic processes.Our preference for cybernetic models has been based on its comprehensive accounting for dynamic regulation, not present in other dynamic approaches.
A full kinetic description of metabolic regulation requires detailed knowledge of its molecular mechanism, which is incomplete in most cases.Alternatively, the cybernetic approach [8] provides a rational description of regulation based on optimal control theory.The cybernetic description of metabolic regulation is based on the assumption that a cell is frugal in using its resources and optimally allocates them among a subset of enzymes to achieve a certain metabolic objective (such as the carbon uptake rate or growth rate).The resulting selective activation of reactions is realized by the cybernetic control variables without introducing additional parameters.
Cybernetic models have been successfully used to perform bifurcation analysis of metabolic systems, such as Klebsiella oxytoca [10], hybridoma cells [11] and Saccharomyces cerevisiae [12].While these analyses were made using lumped reaction networks, it is possible to consider a detailed network structure using the hybrid cybernetic modeling (HCM) framework [13][14][15].Recently, Kim et al. [16] built an HCM of the anaerobic growth of Escherichia coli on glucose and pyruvate.Using this model, they predicted three and five steady states in a chemostat and experimentally validated them.These predictions were made by generating hysteresis curves using continuation methods [17] only at a selected set of parameter values, however.
In this article, Kim et al.'s HCM is subjected to more comprehensive nonlinear analysis for the following two purposes.First, we revisit this model to construct global bifurcation maps over a wide range of parameter space.This will lead to the complete identification of domains where the model exhibits qualitatively different behavior.Second, we develop a practical method that facilitates the nonlinear analysis of the cybernetic models containing the non-differentiable max function.While examples for the rigorous nonlinear analyses of cybernetic models are available in the literature [10,12], we test a more practical method that can readily be realized using an automated software package, such as MATCONT [18,19].
This paper is organized as follows.In the subsequent sections, we provide a summarized description of the Kim et al.'s HCM of E. coli and discuss an idea of introducing an approximate function as a replacement of the non-differentiable max function.Using this idea, we perform comprehensive nonlinear analysis of the model to construct global bifurcation diagrams in a two-parameter space of dilution rate and feed composition.The effect of the total sugar concentration in the feed on bifurcation behavior is also discussed.

The HCM Framework
Dynamic mass balances of extracellular metabolites in a chemostat can be represented as follows: where t is time, c is the biomass concentration, x and x IN are the vectors of n x concentrations of extracellular components in the reactor and feed, respectively, including substrates, products and biomass, r is the vector of n x fluxes, S x is the (n x × n r ) stoichiometric matrix and D is the dilution rate.Under the quasi steady-state approximation, the flux vector, r, can be represented as non-negative (or convex) combinations of basic pathways, termed elementary modes (EMs) [20], i.e., where Z is the (n r × n z ) matrix composed of EMs as its columns and r M is the vector of n z fluxes through EMs.EMs may be viewed as metabolic pathways composed of a minimal set of reactions that can operate alone in steady state.Nonnegative combinations of EMs can represent any feasible metabolic state (i.e., flux distribution) in a network.The cybernetic approach assumes a certain metabolic objective, such as the maximization of the carbon uptake rate (or growth rate) for which metabolic reactions are optimally regulated.The HCM framework views EMs as metabolic options to achieve such an objective and describe metabolic regulation in terms of their optimal combinations.Flux through the jth EM is modeled as regulated by the control of enzyme level and its activity, i.e., , , , , where v M,j is the cybernetic variable controlling enzyme activity, r kin M,j is the kinetic term, and e rel M,j is the relative enzyme level to its theoretical maximum, i.e., e M,j /e max M,j .Enzyme level e M,j is governed by the following dynamic equation, i.e.: where u M,j is the cybernetic variable regulating the induction of enzyme synthesis, r kin ME,j is the kinetic part of the inducible enzyme synthesis rate, β M,j is the degradation rate and µ is the specific growth rate.The four terms of the right-hand side denote constitutive and inducible rates of enzyme synthesis and the decrease of enzyme levels by degradation and dilution, respectively.The cybernetic control variables, u M,j and v M,j , are computed from the Matching and Proportional laws [21,22], respectively: , , ; max( ) where the return-on-investment, ρ j , denotes the carbon uptake flux through the jth EM, i.e., f C,j e rel M,j r kin M,j , and f C,j denotes the factor converting EM flux to the carbon uptake rate.Dynamic shifts among different pathways are realized by two controlling variables, u M,j and v M,j .

HCM for Anaerobic E. coli Growth
Kim et al. [16] used the HCM framework to model the anaerobic growth of E. coli GJT001 on glucose and pyruvate.The metabolic network contains 14 reactions (one reversible and 13 irreversible) and 18 metabolites (eight extracellular and 10 intracellular).Among 49 EMs obtained using METATOOL 5.1 [23], four key modes that can represent yield data of fermentation products [15] were extracted for modeling.Each mode is associated with the consumption of different substrates, i.e., EM1 and EM2 with respective consumption of glucose and pyruvate, while EM3 and EM4 are with simultaneous consumption of both sugars.Model equations and parameters are summarized in Table 1.For a full description of the model, refer to [24].

Methods
Bifurcation analysis of cybernetic models requires special treatment of the non-smooth max function contained in the v M -variables.Among many possibilities, we discuss two ideas of handling this issue, i.e., the combinatoric approach used by Namjoshi and Ramkrishna [10] and the smooth approximation to the max function.

Rigorous Combinatoric Analysis
Namjoshi and Ramkrishna [10] proposed a strategy to enumerate all combinatorial cases, in each of which the model equations are fully differentiable.This leads to four cases by setting one of the v M variables to be 1, while the others are less than or equal to 1 (Table 2).Table 1.Model equations and parameter values.EM, elementary mode.

Equations or parameter values
Parameters and stoichiometric coefficients

Case
In each case, we force v M,j to be 1 by replacing the denominator of v M,j , i.e., , with f C,j e rel M,j r kin M,j , leading to four independent sets of model equations.Figure 1 shows the resulting four hysteresis curves in the D − c space with a fixed value of γ (i.e., 0.2), obtained from the analysis of Cases I to IV, respectively.Segments highlighted in color represent feasible branches satisfying the constraint, i.e., v M,j = 1 (j = 1 − 4), i.e., green (b), cyan (c) and magenta (d), respectively.Note that no such colored branch is found in Figure 1a, indicating that there exists no feasible solution satisfying v M,1 = 1 along the whole profile.Finally, we put together individual pieces of feasible branches of each case to obtain the hysteresis curve over the whole range of D (Figure 2a).Throughout this article, we use colors to distinguish one branch from others characterized with different dominant (i.e., most activated) modes.That is, blue, green, cyan and magenta lines imply that their dominant modes are EM1, EM2, EM3 and EM4, respectively.The black line, on the other hand, indicates the trivial solution with nonzero biomass (i.e., wash-out as marked with W).In Figure 2a, other than typical limit points (also called folds, turning points or saddle nodes), there are two sharp corners (C) (solid red circles), as well.These non-smooth folds represent catch-up points at which the maximally activated mode is overtaken by another.That is, around catch-up points in Figure 2a, the most dominant mode is switched between EM2 and EM4 (left) and between EM3 and EM4 (right).This clearly manifests the pathway shift by regulation.
The shape of the hysteresis curve becomes somewhat different at a higher fractional concentration of glucose, γ = 0.4 (Figure 2b).The dominant mode at lower values of D is EM1 (instead of EM2), the mode taking up glucose only.Interestingly, one of the two catch-up points (open red circle) does not correspond to a limit point.Thus, we differentiate this simple transition (T) (open red circle), which does not form a sharp limit point, from non-smooth catch-up points (C).The existence of simple transition points has not been reported in earlier studies using lumped network-based cybernetic models.As simple transitions are not bifurcation points, we do not trace them.

Smooth Approximation to the Max Function
While the combinatoric approach described above allows for rigorous bifurcation analysis in theory, it is ineffective in cases where the number of EMs is large.Alternatively, we may mollify the pain of handling non-smooth functions by making smooth approximations.L p -norm is considered as an accurate approximation to the max function when p is sufficiently large.That is, we may approximate Figure 3 shows the reproduction of the hysteresis curve using the L p -norm approximation with different p-values.No appreciable errors are found when p ≥ 30, while some deviations are observed when p-values are lower than that.
Figure 4 provides an enlarged view of two red windows in Figure 3 around the catch-up points.Approximate models progressively approach the rigorous solution (obtained with the combinatoric method described above) as the value of p increases.While small deviations are unavoidable regardless of how large p is, these tiny errors of below 1 percent are acceptable.Stable and unstable branches are also successfully reproduced using this approximate function.In all calculations hereafter, therefore, we use the L p -norm approximation with a p value of 70.
The usefulness of the smooth approximation depends on THE cases in consideration [25].In a number of studies, introduction of smooth approximation facilitated the bifurcation analysis by providing the system with global differentiability.On the other hand, approximate functions may become stiffer to integrate or may generate more complex bifurcation diagrams than the original function.Thus, it would be critical to have a previous check if smooth approximation yields any unexpected difficulties or errors.If the approximate representation is acceptable as in our case, nonlinear analysis of piecewise-smooth functions is greatly facilitated by using an automated software, such as MATCONT, a standard continuation software package [18].

Integration of Two Methods
As a compromise, we may integrate combinatoric enumeration and smooth approximation.That is, we can sketch a bifurcation diagram conveniently using the smooth approximation and refine non-smooth folds using rigorous computations based on the combinatoric approach, because they are only the regions where errors may occur.Catch-up points are readily identified from the hysteresis curve using the approximate function.This combined approach is more accurate than the approximate function alone and more convenient than the full combinatoric enumeration.

Results and Discussion
Among three methods discussed in the previous section, we use the L p -norm approximation (Section 3.2) to explore the nonlinear behavior of the HCM model by Kim et al. presented in Table 1.The smooth approximation is conveniently implementable with no appreciable errors in our case.The main parameters subject to variation include dilution rate (D) and the fractional molar concentration of glucose in the feed (γ), i.e., , , , where F is the volume flow rate of the feed, V is the culture volume and x IN,G and x IN,P are concentrations of glucose and pyruvate in the feed, respectively.The total sugar concentration (x IN,total ) is the sum of x IN,G and x IN,P .

Hysteresis Behaviors and Bifurcation Diagram
Figure 5 shows the concentration profiles of all components (including glucose, pyruvate, biomass, formate, acetate and ethanol) at a specific condition (i.e., γ = 0.4 and x IN,total = 50).The implication of different colors and solid and dashed lines is the same as before.This parameter set yields up to five steady states in a range of D between 0.325 and 0.335.A catch-up point is observed between EM3 and EM4.
To get a global bifurcation diagram, we explore the whole parameter space spanned by D and γ.In the comprehensive search of all possible bifurcation points using MATCONT, we ended up with only two different kinds of bifurcations: limit and catch-up points.No Hopf bifurcation was detected.That is, the nonlinear behaviors we could identify are limited to steady-state multiplicity, and no existence of metabolic oscillation is found.
Figure 6 shows a global map of multiplicity when x IN,total is fixed to 50 mM.It shows two closed curves in black and gold (left) and four pairs of lines highlighted in the same color, respectively (right).The gold curve represents the neutral saddles, equilibrium points characterized by two real eigenvalues with the opposite sign.Neutral saddles are, however, not bifurcation points of interest and have nothing to do with steady-state multiplicity.Solid lines (other than neutral saddle lines) represent typical limit points, while thick dotted lines, catch-up points.Therefore, inside each envelop, there exist three multiple steady states (i.e., domains I, II, III, IV and V), at least.In the region where two envelops overlap (i.e., domains VI, VII and VIII), five steady states exist.In the remaining region, a unique solution exists.

Bifurcation Diagram at a Higher Sugar Concentration in the Feed
The effect of the total sugar concentration in the feed (x IN,total ) on nonlinear behavior of the E. coli model is examined.When lowering the total sugar concentration from 50 to 25 mM, no qualitative change is observed in bifurcation behavior.Increasing x IN,total to 100 mL, on the other hand, leads to an additional domain not observed previously.
Figure 8 shows a global bifurcation diagram in the D − γ space at x IN,total of 100 mL.The implication of lines and colors is the same as before.Unlike the previous case, this condition leads to multiplicity regimes with up to seven steady states.That is, seven steady states emerge in the domain (IX) where three different envelops (i.e., red, orange and purple ones) are overlapped.We have highlighted this domain in the figure .Figure 9 shows hysteresis curves of biomass (a) and pyruvate (b) concentrations when γ = 0.83.We can see seven steady states in a small range of D around 0.34.The figure also provides zoomed-in views of seven steady states existing between two vertical dashed lines.For instance, pyruvate can have seven different concentrations in this domain, i.e., if enumerated from the top, one on the black line, one on the cyan, two on the magenta and three on the blue.

Experimental Validation
Kim et al. [16] have provided an experimental verification of the stable steady states for the foregoing two sets of conditions, i.e., at γ = 0.2 and x IN,total = 50, yielding a total of three steady states and five steady states, γ = 0.4 and x IN,total = 25, with a total of five steady states.
Through the comprehensive bifurcation analysis in this work, we could identify a new domain with seven steady states.Experimental verification would require precise control of conditions and concentration measurements, however.

Conclusions
Since the cybernetic variables for enzyme activity control are max functions and, therefore, non-smooth, nonlinear analysis of cybernetic models has had to rely on a suitably convenient methodology to confront this issue.The L p -norm approximation of the max function tested in this work is a practically useful idea, as it is applicable to general cases considering a large number of metabolic pathway options (i.e., EMs).Replacement of the max function with the L p -norm representation allows for accurate computation of bifurcation points.While slight errors around non-smooth folds (or catch-up points) are unavoidable, they are negligibly small in our case.When these errors are appreciable in certain cases, however, we can redo rigorous computation only for the non-smooth folds based on the combinatoric idea of Namjoshi and Ramkrishna [10].Such a combination of these two methods guarantees rigorous results at a minimal level of inconvenience, thus serving as a promising strategy.
Using the approximate function, we could construct global bifurcation diagrams on the D − γ space to identify various multiplicity domains, including the one with seven steady states.Considering the narrowness of that domain, however, there are some issues to be resolved for experimental validation.Despite such a comprehensive analysis performed in this work, dynamic nonlinear behaviors, such as metabolic oscillations, were not detected.This could be a consequence of the condensed set of elementary modes in the model.
Clearly, more detailed models comprising more EMs could produce a considerably greater number of steady states, which may be difficult to observe experimentally without accurate analytical measurements and precise control of experimental conditions.It is our premise that this is an area for extensive future exploration by researchers concerned with modeling metabolism.

Figure 3 .
Figure 3. Reproduction of the hysteresis curve of Figure 2a using the L p -norm approximation with different p-values.

Figure 4 .
Figure 4. Magnified views of two red windows around the catch-up points in the lower-right panel of Figure 3: (a) left upper window, (b) right lower window.

Figure 6 .
Figure 6.A global bifurcation diagram in the D − γ space when x IN,total = 50 mM.

Figure 7 .
Figure 7. Hysteresis of biomass concentration profiles with different γ values.

Figure 8 .
Figure 8.A global bifurcation diagram in the D − γ space when x IN,total = 100 mM.
: factor converting the EM flux (i.e., growth rate) to the carbon uptake rate, C-mmol/gDW (DW = dry weight) k F : rate constant for formate decomposition k max j : maximal rate constant for the jth EM flux, 1/h K F : Michaelis constant for formate decomposition, mM K G,j , K P,j : Michaelis constants for the jth EM flux, mM r F : specific rate of formate decomposition into CO 2 and H 2 , mmol/(gDWh) r M,j , r kin M,j : regulated and unregulated fluxes through the jth EM, mmol/(gDW h) r kin ME,j : kinetic part of inducible enzyme synthesis rate, 1/h s A,j , s E,j , s F,j , s G,j , s P,j : stoichiometric coefficients, mmol/gDW t: time, h u M,j : cybernetic variable regulating the enzyme induction v M,j : cybernetic variable regulating the enzyme activityx A , x E , x F , x G ,x P : concentrations of acetate, ethanol, formate, glucose and pyruvate, mM x IN,G , x IN,P : feed concentration of glucose and pyruvate, mM Notations c: biomass concentration, g/L D: dilution rate, 1/h e M,j , e max M,j : level of enzyme that catalyzes the jth EM flux and its maximal level f C,j

Table 2 .
Four combinatorial cases for Kim et al.'s E. coli model.