Model-based Identification of Some Conditions Leading to Glycolytic Oscillations in E . coli Cells

Autonomous oscillations of glycolytic intermediate concentrations reflect the dynamics of the control and regulation of this major catabolic pathway, and this phenomenon has been reported in a broad range of bacteria. Understanding glycolytic oscillations might therefore prove crucial for the general understanding of the regulation of cell metabolism with immediate practical applications, allowing in silico design of modified cells with desirable ‘motifs’ of practical applications in the biosynthesis industry, environmental engineering, and medicine. By using a kinetic model from literature, this paper is aiming at in silico (model-based) identification of some conditions leading to the occurrence of stable glycolytic oscillations in the E. coli cells.


Introduction
Autonomous oscillations of the glycolytic intermediates' concentrations reflect the dynamics of the control and regulation of this major catabolic pathway, and the phenomenon has been reported in a broad range of cell types. 1 Understanding glycolytic oscillations might therefore prove crucial for our general understanding of the regulation of metabolism and the interplay among different parts of metabolism, as illustrated, for instance, by the hypothesis that glycolytic oscillations play a role in complex pulsatile insulin secretion, 2 or in the amino-acid synthesis. 3In this context, the key question concerns the mechanism(s) of the oscillations but, "despite much work over the last 40 years, it remains unsettled". 1 Glycolysis is an essential part of the cell metabolism.5][6] The CCM model is the 'core' part of any systematic and structured model-based analysis of the cell metabolism with immediate practical applications, such as target metabolite synthesis optimization, insilico re-programming of the cell metabolism to design new microorganisms for industrial bioprocess optimization, etc. [4][5][6][7][8][9][10][11][12] However, to cope with the astronomic complexity of cellular processes, of low observability, involving O(10 3 -10 4 ) number state variables (species conc.),O (10 3 ) gene expression transcription factors TF, and O(10 4 -10 5 ) reactions, versatile models of 'building-blocks' like modular constructions, including individual and lumped species and reactions have been developed over decades. 5,7[10] Consequently, understanding and simulation of the cell characteristics and environmental conditions leading to an oscillating glycolysis is an old subject, but still of high interest. 4To simulate the glycolysis in bacteria, a large number of glycolysis models, of a reduced or extended form, have been proposed over decades (review of Maria 4 ).Recently, Maria 4 proposed a reduced glycolysis model (denoted by mTRM in the E. coli cells, including only 9 species, 7 lumped reactions, and 17 estimable parameters.This model was identified using experimental dynamic data from literature, 6,18 and has been proved to adequately reproduce the cell glycolysis under steady state, oscillatory, or transient conditions according to the defined glucose input flux, environmental conditions, the total A(MDT)P cell energy resources, and cell phenotype characteristics (related to the activity of enzymes involved in the ATP utilization and recovery system).This paper is aiming at using the mTRM model of Maria 4 to simulate some conditions leading to glycolytic oscillations in the E. coli cells.Because this model is one of good adequacy, relevant results are expected.

Kinetic model of glycolysis in the E. coli prokaryotic bacteria
Glycolysis (from an older term with the meaning of glucose degradation) is the metabolic pathway that converts glucose (C 6 H 12 O 6 ) into pyruvate (CH 3 COCOO − + H + ).The free energy released by the subsequent tricarboxylic acid cycle (TCA) originating from pyruvate is used to form the high-energy molecules ATP (adenosine triphosphate), and NADH (reduced nicotinamide adenine dinucleotide) that support the glycolysis and numerous enzymatic syntheses into the cell. 13,14lycolysis is a determined sequence of ten enzyme-catalyzed reactions (Fig. 1).The intermediates provide entry points to glycolysis.For example, most monosaccharides, such as fructose or galactose, can be converted to one of these intermediates.The intermediates may also be directly useful.For example, the intermediate dihydroxyacetone (DHAP, an intermediate in the reaction of f6p conversion to g3p in Fig. 1) is a source of the glycerol that combines with fatty acids to form fat. Also, NADPH is formed by the pentose-phosphate pathway (PPP), which converts glucose into ribose, which is used in the synthesis of nucleotides and nucleic acids.Pep is also the starting point for the synthesis of essential aminoacids such as tryptophan, cysteine, arginine, serine, etc. 15,16 Due to the tremendous importance of this metabolic process in simulating the cell CCM, intense efforts have been made both in the experimental study, and in modeling the glycolysis dynamics in Escherichia coli, [17][18][19][20] and other cell types.
On the other hand, to model in detail the dynamics and regulation of glycolysis, which is one of the most complex cellular processes, is a very difficult task, if not impossible.Consequently, a large number of glycolysis reduced or extended kinetic models have been proposed over decades (review of Maria 4 ), of a complexity ranging from 18-30 species, included in 48-52 reactions, with a total of 24-150 parameters.Most of these models are however too complex to be easily used and identified.Besides, with a few exceptions, most of the models cannot satisfactorily reproduce the glycolytic oscillations on a mechanistic basis.
However, starting from an extended reaction pathway and model, and by applying chemical engineering lumping techniques, 12,21,22 Maria 4 pro-posed a valuable reduced dynamic model of glycolysis (denoted by mTRM), accounting for 9 species, 7 lumped reactions, and including 17 easily identifiable parameters.The rate constants of this model have been identified using the kinetic experimental data of Chassagnole et al., 23 and Visser et al. 6 The mTRM model is presented in Table 3.The model has been proved to adequately reproduce the cell glycolysis under oscillatory, or transient conditions according to the defined glucose concentration in the bioreactor, the total A(MDT)P cell energy resources, and the cell phenotype characteristics (concerning the ATPase enzyme activity, this essential enzyme being involved in the ATP utilization and recovery system).This is why the bioreactor (of model presented in Table 2) and the glycolysis dynamic models have to be considered together (Tables 2-3) when simulating the cell CCM.

How glycolytic oscillations occur
Oscillations in chemical systems represent periodic state variable (i.e.species concentrations) transitions in time.
According to Franck, 24 spontaneous occurrence of self-sustained oscillations in chemical systems is due to the coupled actions of at least two simultaneous processes.Oscillations sourced in a so-called "oscillation node" (that is, a chemical species, or a reaction), on which concomitant rapid positive (perturbing) and slow negative (recovering) regulatory loops act.Because the coupling action between the simultaneous processes is mutual, the total coupling effect actually forms closed feedback loops for each kinetic variable involved.There exists a well-established set of essential thermodynamic and kinetic prerequisites for the occurrence of spontaneous oscillations.In short, according to Franck, 24 these are the following: 1) Sustained oscillations can only occur in thermodynamically open systems far from equilibrium; 2) Oscillatory systems always consist of more than one degree of kinetic freedom, i.e. the description of their temporal behaviour requires a corresponding set of simultaneous differential equations; 3) Oscillations occur as a result of simultaneous feedback effects; 4) There exist extremely nonlinear relationships between the involved driving forces and driving fluxes or reactions; 5) Oscillatory systems always contain unstable states; 6) Oscillations are the result of mutual kinetic coupling between processes being otherwise independent from each other; 7) Once an oscillation occurs, it propagates in the whole reaction pathway; 8) Feedback occurs when a process acts kinetically upon itself; it therefore consists basically of a closed chain of action which causes the well-known effects of self-enhancement in the case of "positive feedback" (denoted by ), and self-inhibition in the case of "negative feedback" (denoted by ) respectively, in a non-systemic, or systemic feedback (i.e.rate constants depend on their own reaction products or reactants); 9) In general, there are four possible modes of coupling control loops leading to oscillatory processes, the positive ( ) or the negative ( ) feedbacks simultaneously acting on the synthesis and consumption of the oscillating "node species", that is: i) positive, and negative feedbacks; ii) positive feedback, and positive feedforward; iii) negative feedback and negative feedforward; iv) negative feedforward and positive feedforward.Such an oscillation "engine" of type (iv) is displayed in Fig. 2 for the glycolysis case; 10) Chemical oscillations exhibit positive and negative feedback simultaneously; according to the "principle of antagonistic feedback of chemical oscillators".Oscillations are understood as a consequence of an antagonistic interaction of a relatively fast acting positive feedback of labilizing tendency and a slower acting negative feedback of a recovering tendency for stabilizing the system; 11) Oscillations occurrence and characteristics depend not only upon the presence of both kinds of feedbacks but also upon the correct ratios of the time parameters (rate-constants) of the feedback loops involved.
In the glycolysis system case, extensive experiments have revealed that self-sustained oscillations are reported in a broad range of cell types 1 (see for instance in Fig. 3 the plotted experimental glycolytic oscillations measured by Schaefer et al. 25 in E. coli, and plotted by Chiarugi et al. 26 ).As revealed by Termonia and Ross, 27,28 Vermeer, 29 and Chiarugi et al. 26 glycolytic oscillations occurrence is due to the antagonistic action of two processes on regulating the V 2 reaction rate that converts F6P into FDP (see reaction scheme in Fig. 2).Valuable contributions to explain and model, on an experimental basis, the glycolytic oscillations "engine" and their self-control, have been done on E. coli, 27,28,[30][31][32][33] or on yeast. 34lycolytic oscillation occurrence and characteristics (period) are influenced by both external (environmental) and internal (genomic) factors, that is: 35,36 I) From one side it is the glucose (Glc) import driving force through the phosphotransferase (PTS)-system (Fig. 1) regulated and triggered by the external concentration of glucose c ext GLC = [Glc]ext F i g . 2 -Chemical node inducing glycolytic oscillations (after 27,28 ).Notations , and denote the feedback positive or negative regulatory loops, respectively.Glc = glucose; F6P= fructose-6-phosphate; FDP = fructose-1,6-biphosphate; V 1 -V 3 = reaction rates in Fig. 1.
and by the PEP and PYR levels (see the V 1 rate expression in Table 3).II) However, the Glc import and conversion to PYR requires important amounts of regenerable ATP, and a sufficiently rapid ATP to ADP conversion rate, as well as its quick regeneration (Fig. 1 and rate expressions in Table 3).III) On the other hand, the limited A(MDT)P cell energy resources exist in the cell, slowing-down the Glc import if the ATP use/regeneration is not working "fast enough". 36V) Additionally, due to the enzymes ATP-ase and AKase characteristics related to the bacteria genome and cell phenotype, a limited ATP conversion rate can sustain the glycolytic reactions, while the ATP recovery rate is limited by the enzymes participating in the A(MDT)P inter-conversion reactions (K and k 6 , constants in Table 3).
V) At the same time, glycolysis being a systemic process with a complex regulatory structure, oscillations are also related to the rate constants of all the involved reactions.Similarly, Silva and Yunes 37 found that oscillations are only possible if the Glc concentration, and the maximum reaction rates controlled by the PFKase and GKase enzymes (involved in the V 1 reaction, i.e., the PTS import system) are within specific intervals.VI) Consequently, among the glycolytic oscillation factors, it is natural to approach in this study in the first place the influence of the external factor [Glc]ext, and of some internal factors, such as the k 6 rate constant (belonging to the A(MDT)P inter-conversion system) on the glycolytic oscillations.The total [AMDTP] level was kept constant in this study to highlight the influence of the other mentioned factors.
Roughly, the same conclusions have been underlined by Silva and Yunes 37 (even if in the S. cerevisiae case): "It appears that glycolytic oscillations are focused on the maintenance of energy levels in the cell (negative regulation of PFKase by ATP) and thus the ability to limit the conversion into energy in situations where it is not needed.Therefore, it would be more advantageous to store it or deviate the flux towards other cell cycle activities such as cell division."Consequently, mutant cells with modified enzymes activity (especially PFKase, PKase, ATPase, AKase, GKase) will lead to a noticeable modification in the cell metabolism.
Here it is important to mention the works 30,38 dealing with explaining specific regulation of the glucose influx by PTS system in E. coli glycolytic, modelled in detail by Chassagnole et al., 23 and Visser et al., 6 and also the biochemical interactions among the PTS system, and the components of the ATP regeneration pathway (e.g.PFKase, PKase), used to evaluate the dynamic behavior of the glycolytic pathway of E. coli under steady-state conditions.1][32][33] The advantage of the mTRM model of Maria 4 is its capacity to reproduce glycolytic oscillations using a reduced kinetic model but still preserving the essential glycolytic and environmental parameters with a major influence on the process (see the next chapter on oscillation conditions).

Simulation of some oscillation conditions
By adopting the glycolysis kinetic model of Maria 4 , one can determine, by repeated simulations, the cell external and internal conditions leading to glycolytic oscillation occurrence.In simulations, one considers the E. coli cell growing conditions of the semi-continuous bioreactor of Chassagnole et al. 23 given in Table 1 (using sparging air in excess, and necessary nutrients for a cell culture equilibrated growth).The main mass balance equations of the bioreactor and glycolysis dynamic model are presented in Table 2. To obtain the model solution with enough precision, a low-order stiff integrator ("ODE23S" routine) of the Matlab™ package was used.
Simulations were made for the cell culture conditions given in Table 1, and for cells with [AM-DTP]total = 5.82 mM. 4,23Following the discussion in the previous chapter on oscillations occurence, the influence of two main factors is studied here, that is: F i g .3 -Experimentally measured glycolytic oscillations of fructose-6-phosphate (F6P) and fructose-1,6-bisphosphate (FDP)."In spite of white gaussian noise, the two plots have a clear constant period and amplitude, showing a stationary oscillatory pattern" in E. coli (adapted after Chiarugi et al. 26 similar results have been obtained by Schaefer et al. 25 ).Time axis in minutes.Concentrations are in mM.

i) [Glc]ext (related to the bioreactor operating conditions);
ii) k 6 reaction rate (determined by the ATP-ase characteristics, related to the cell phenotype); iii) all other reaction rate constants, and [AM-DTP] level were kept unchanged during simulations of values given in Tables 1, 3. Culture dilution rate (F L /V L ) 1.667•10 -3 min -1 (adjusted to be identical to D)  27,28 iii) c ADP results from solving the thermodynamic equilibrium relationship

-Glycolytic stationary oscillation domains (thick lines) in E. coli in the plane [Glc]ext (at t = 0), and k 6 for the bioreactor operating conditions in Table 1 [AMDTP] = 5.82 mM
Ta b l e 3 -Glycolysis kinetic model mTRM of Maria 4 and its parameters (the units are in mM, min)  ( / )   be quasi-constant in the present case study), are determined by the ATP to ADP conversion rate, and ATP regeneration rate (reflected here by k 6 , and K constants of Table 3).II) Oscillations occur for low [Glc]ext but with a slow Glc import, due to relatively low k 6 constant values (i.e., a cell with a slow ATP conversion to ADP and ATP recovery).III) By contrast, high levels of [Glc]ext, triggering high rate import, produce glycolytic oscillations for larger values of k 6 , due to the limited ATP recovery rate (k 6 being also related to the K constant governing the AMDTP pathway).Eventually, for too small, or too large k 6 values, the glycolysis reaches its steady-state.
IV) The glycolytic oscillation domains in Fig. 4, plotted in terms of k 6 and [Glc]ext, are very nar-row, revealing their high sensitivity with respect to the inducing factors, and their poor stability (as proved by the limit cycles plotted in Figs.5-6).As expected, such a result indicates that oscillations stability is also dependent on the microorganism characteristics.For instance, by contrast, the glycolytic oscillations in yeast have been proved 39,40 to be very robust even in the presence of environmental noise, "oscillations being a side-effect of the tradeoffs between robustness and regulatory efficiency of the feedback control of the autocatalytic reaction network".
V) As may be observed in Table 4, the increasing values of k 6 have, as a consequence, a slight decrease in the oscillation period until oscillation disappearance.This effect is better illustrated in Fig. 6 obtained for [Glc]ext = 1 mM (at t = 0), and k 6 = 16 min -1 , with an oscillation period of ca.0.5 min.These oscillations are amortized and, eventually disappear (plots for larger operating times are not presented here) due to the sharp decline of [Glc]ext from the initial 1 mM level (Fig. 6 down-left), due to its consumption by the cells, and its washout.VI) For comparison, the simulation result plotted in Fig. 5, which was obtained for [Glc]ex = 0.0557 mM; k 6 = 12 min -1 , presents a smaller oscillation period of ca.0.33 min, and a higher stability due to a smaller Glc environmental level and a higher ATP use/recovering rate.Table 4 also reveals that oscillation period takes values in the range of 0.4-0.9min, being smaller as k 6 is bigger, and [Glc]ext is smaller.For comparison, experimentally determined glycolytic oscillations present periods of ca.0.2 min, 1 or 2-100 s, 35 15 s, 37 or 1-20 min, 36 up to 3 h, 41 or 0.2 min to hours. 30ven if the checked time-interval, oscillation period, and bioreactor conditions of our simulations given in Fig. 5 are different of those experimentally checked in Fig. 3 (not mentioned by the authors), the theoretical curve shapes for FDP and F6P are similar to the experimental ones, even if the timescale of the abscissa is very different (400 min in Fig. 3, compared to 10 min in Fig. 5).The slight increase in the amplitude of oscillations of FDP and F6P are similar to the simulation results of Selkov, 33 Bier et al., 36 Elias, 42 de la Fuente, 43   in the two operating cases.As expected, when the initial concentration of Glc is small (0.0557 mM), the high Glc level in the feed (200 mM) ensures a relatively quick fulfilment of the reactor steadystate conditions (of 0.188 mM in Fig. 5 bottom-left).By contrast, when the initial concentration of Glc is higher (1 mM), the reactor transition toward the steady-state (of 0.28 mM in Fig. 6 bottom-left) is slower.
As for all in silico analyses, the results are strongly dependent on the used model quality.
To summarize, simulations varying the considered search variables in Table 4, and Fig. 4, clearly reveal the overwhelming importance of the environmental level of Glc, and of cell phenotype (that is, cell genomic and phenotype factors controlling the [AMDTP] total energy resources level, the glycolytic synthesis and regulation reactions, and especially the activity of enzymes involved in the A(MDT)P inter-conversion system).At the same time, glycolysis being a systemic process, with a complex regulatory structure, oscillations are also related to the rate constants of all the involved reactions and their appropriate ratios.
For the defined input data and model hypotheses, one can conclude that the derived results are conclusive enough.That is mainly because the used glycolysis mTRM model presents very good adequacy vs. independent experimental data used for its identification by Maria. 4 Consequently, one may conclude that the in-silico analysis results present satisfactory confidence and relevance.Of course, other variables, not accounted for in the model (cell characteristics reflected in the model constants) can influence the location of the problem solution.Subsequent experimental checks can validate the estimated glycolytic stationary oscillation domains in E. coli and, eventually, in the case of inconsistencies, they will lead to the model updating/completions for correcting its adequacy in order to perform future in silico analyses.

Conclusions
The use of reduced kinetic models describing the dynamics of complex metabolic pathways is a continuous challenging subject when developing structured cell simulators for various applications (flux analysis, target metabolite synthesis optimization, in silico re-programming of the cell metabolism for microorganism design purposes, bioreactor and bioprocess optimization).As exemplified by the E. coli glycolysis case study, the reduced mTRM model, of simple and easily adaptable structure to various cell cultures, can be used in quick analyses of the cell metabolism, such as substrate utilization, oscillation occurrence, or structured interpretation of metabolic changes in modified cells.
Reduced structured glycolysis models of satisfactory adequacy for the key-species are preferred to other semi-empirical or very extended models, being easily included in modular cell simulation platforms used for solving various problems, such as: analysis of cell adaptation to certain environmental conditions; simulation of genetic regulatory circuits controlling the synthesis of some target metabolites; simulation of metabolic flux distribution and its dynamics under transient regimes; in silico reprogramming of some metabolic pathways to design new microorganisms. 5,7erivation of reduced kinetic structures to describe some parts of the cell core metabolism is worth the associated identification effort, due to the considerable reduction in the model parameterization (e.g., 17 parameters in the mTRM model vs. 127 in the Chassagnole et al. 23 model), while preserving a fair adequacy over a wide experimental domain.Besides, when cell characteristics change significantly, the reduced model can be upgraded quickly using the available experimental informa-Ta b l e 4 -Some cell external and internal conditions leading to glycolytic oscillations occurrence in the E. coli cells.Total [AMDTP] = 5.82 mM [nominal conditions of Chassagnole et al. 23 and Maria 4 presented in Ta b l e 1 ].NO = evolution to quasi-steady-state with no oscillations.
[Glc]ext (at t = 0), (mM) k 6 (min -1 ) Approx.oscillation period (min.)tion.Thus, the cell metabolic process complexity appears to be described by a succession of "locally" reduced models "enfolded" on the real process long time-interval dynamics. 4,11eing quite versatile, the reduced mTRM model includes enough information to reproduce not only the cell energy potential through the total A(MDT)P level, but also the role played by ATP/ ADP ratio as a glycolysis driving force.The model can also reproduce the oscillatory behaviour depending on the external/internal conditions, as well as situations when homeostatic conditions are not fulfilled.
The glycolysis core model can be easily extended by including any complex synthesis and regulatory pathway deriving from the main carbon uptake stream (e.g.CCM, succinate, SUCC, amino-acids production; Maria 5,16 ), without necessarily complicating the 'core' model with too many species and parameters of less importance for the target metabolite production.
Simulations were conducted in an exhaustive way, by covering the ranges of the initial [Glc]ext = [0.01-1]mM (at t = 0), and k 6 = [0.1-20]min -1 .The results are presented in Table 4, and plotted in Fig. 4.These simulation results lead to several model-based conclusions: I) Oscillations are basically determined by the external level of [Glc] (triggering the glucose import into the cell) but also, for a certain [AMDTP], total energy resources level in the cell (assumed to Ta b l e 1 -Operating conditions of the Chassagnole et al. 23 semi-continuous bioreactor with suspended E. coli cell culture used to simulate the glycolytic oscillation occurrence Parameter Value Biomass concentration (C x ) 8.7 gDW L -1 culture volume Cell content dilution rate (D) 1.667•10 -3 min -1

e 2 -
Bioreactor and glycolysis mass balance eqns.for the kinetic model of Maria4 species initial concentrations are those measured by Chassagnole et al.,23 that is (in mM): ( 0) ext GLC c t = = tried reference value of 0.0557 mM, or 1 mM, 6 formation from PYR has been neglected from this model; v) biomass concentration (C x ) is assumed to be quasi-constant.
Termonia and Ross27,28 indicated experimental evidence of a very fast reversible reaction catalysed by AKase, the equilibrium being quickly reached.Obs.: Other values of k 6 are also possible (to be investigated), according to the microorganism phenotype (characteristics of the gene encoding the enzyme ATPase that catalyse this reaction).
etc. VII) Also noteworthy are the dynamics of the [Glc]ext in the bioreactor bulk phase (a semi-continuous bioreactor with suspended E. coli cell culture)

A
b b r e v i a t i o n s a n d n o t a t i o n s