In silico Determination of Some Conditions Leading to Glycolytic Oscillations and Their Interference With Some Other Processes in E. coli Cells

Introduction

“Autonomous oscillations of species levels in the glycolysis express the self-control of this essential cellular pathway belonging to the CCM, and this phenomenon takes place in a large number of bacteria.

The study of glycolytic oscillations might, therefore, prove crucial for the general understanding of the cell metabolism regulation and the connections among different parts of metabolism. The key question in this context is the mechanism of the oscillations but, despite much work over the last 40 years, it remains unsettled (Wierschem and Bertram, 2004; Madsen et al., 2005).

A model able to simulate the dynamics of the cell CCM must include linked modules relating to (i) the glycolysis (Figures 1, 2); (ii) the phosphotransferase (PTS) system for GLC import into the cell (Figure 1); (iii) the pentose-phosphate pathway (PPP) to generate NADPH and pentoses (5-carbonsugars), as well as ribose 5-phosphate (R5P, a precursor for the synthesis of nucleotides); (iv) the tricarboxylic acid cycle (TCA); (v) the ATP recovery system, and several other pathways” (Palsson, 2005; Kadir et al., 2010; KEGG, 2011; Maria, 2018).

FIGURE 1

Figure 1. Simplified representation of the CCM pathway in E. coli of Edwards and Palsson (2000). Fluxes characterizing the membranar transport [Metabolite(e) ↔ Metabolite(c)] and the exchange with environment have been omitted from the plot (see Maria et al., 2011 for details, and for explanations regarding the numbered reactions). [e], environment; [c], cytosol. Adapted from Maria et al. (2011, 2018c) courtesy of CABEQ Jl. The considered 72 metabolites, the stoichiometry of the 95 numbered reactions, and the net fluxes for specified conditions are given by Maria et al. (2011). The pink rectangle indicates the chemical node inducing glycolytic oscillations (after Termonia and Ross, 1981a,b; see also Maria et al., 2018c). Notations + and − denote the feedback positive or negative regulatory loops, respectively. Glc, glucose; F6P, fructose-6-phosphate; FDP, fructose-1,6-biphosphate; V1–V3, reaction rates indicated in Figure 2.

FIGURE 2

Figure 2. Simplified representation of the structured reaction pathway of glycolysis (Maria, 2014), and of the TRP synthesis (in the yellow area) (Maria et al., 2018a) in E. coli used by Maria et al. (2018a,b,c) to derive the process kinetic model and the operating conditions of the fed-batch bioreactor that maximize the TRP synthesis. Connection of the TRP synthesis to glycolysis is realized through the PEP node (Maria et al., 2018c; Mihalachi and Maria, 2019). The modular model structure also includes the adenosine co-metabolites ATP, ADP, and AMP synthesis as part of the ATP recovery system (the pink rectangle). Notations: Species in parenthesis are not explicitly included in the glycolysis model. Italic letters denote the enzymes. Squares include notations of enzymatic reactions V1–V5 included in the glycolysis model. Adapted from Maria (2014; 2018b) courtesy of CABEQ Jl. GLC(ex), glucose in the cell environment; G6P, glucose-6-phosphate; F6P, fructose-6-phosphate; HK-ASE, hexokinase; PFK-ASE, phosphofructokinase; ATP-ASE, ATP monophosphatase; ADP, adenosine diphosphate; ATP, adenosine triphosphate; AMP, adenosine monophosphate; AK-ASE, adenylate kinase; Pi, phosphoric acid; FDP, fructose-1,6-biphosphate; G3P,GAP, glyceraldehyde-3-phosphate; 13DPG=PGP, 1,3-diphosphoglycerate; 3PG, 3-phosphoglycerate; 2PG, 2-phosphoglycerate; PEP, phosphoenolpyruvate; PYR, pyruvate; SUCC, succinate; NAD(P)H, nicotinamide adenine dinucleotide (phosphate) reduced; CIT, citrate; ACCOA, acetyl-coenzyme A; LAC, lactate; ETOH, ethanol; AC, acetate.

“Modeling bacteria CCM, or parts of CCM, is a subject of very high interest, because the CCM is the essential part of any systematic and structured (model-based) analysis of the cell metabolism with immediate applications, such as biosynthesis optimization, metabolic fluxes evaluation (Stephanopoulos et al., 1998), model-based design of GMO with target characteristics of various applications in the industry, medicine, etc. (Hatzimanikatis et al., 1996; Visser et al., 2004; Styczynski and Stephanopoulos, 2005; Hempel, 2006; Maria and Luta, 2013; Maria, 2018).” Modeling of the CCM kinetics in close connection to the bioreactor environment dynamics falls at the border of several research areas, that is: (bio)chemical physics, computational biochemistry, biochemical engineering, and physical (bio)chemistry.

“To cope with the very high complexity of cell metabolic processes, involving ca. 10⁴ species concentrations, 10³ gene expression transcription factors, and ca. 10⁵ enzymatic reactions, adaptable reduced dynamic models, of ‘building-blocks' like modular construction, have been developed over the last decades (reviews of Maria, 2017a,b), with including individual/lumped species and/or reactions.

Modeling the glycolysis dynamics is of particular interest, because most of its intermediates are starting points for the cell production of various metabolites of industrial/medical use (e.g., amino acids, succinate, citrate, etc.; Maria, 2018).

By using two adequate dynamic models validated in previous studies (Tables 2, 3), this paper exemplifies how the model-based analysis can be used (i) to predict some of the internal/environmental conditions inducing glycolytic oscillations in the Escherichia coli culture grown in a fed-batch bioreactor (FBR), and (ii) to simulate the influence of the glycolytic oscillations on the TRP oscillatory synthesis by means of the key-species PEP shared by the two oscillatory processes. Industrial applications are immediate seeking for the TRP production optimization. The glycolytic oscillations occurrence will be analyzed vs. external ([GLC]ex in the bulk-phase), and internal factors (that is the ATP recovery rate, dependent on the cell phenotype) (Mathews et al., 1999; Maria, 2014).” The main elements of novelty of this paper consists in the following: (i) Testing and highlighting the influence of the connection of the two metabolic pathways (glycolysis and TRP synthesis), of the external factors ([GLC]ex in the bulk phase, cell dilution rate (μ) adjusted to be equal to the bioreactor dilution rate D), and of the cell internal factors (activity of some enzymes in the ATP recovery system) on the two approached oscillatory processes in E. coli cells. Several operating policies have been checked in this respect (ii) to eventually suggest how to modulate these most influential factors ([GLC]ex, k6, D = μ, and others) before deriving a future optimal operating policy of the FBR (by obtaining, for instance, a timestep-wise GLC feeding policy, not approached here) leading to maximization of the TRP synthesis; (iii) accomplishment of the sensitivity analysis of the studied FBR, which is the preliminary step for any future TRP bioreactor production maximization.

The Tested FBR Bioreactor

“The in silico study of the glycolytic oscillations connected to the TRP synthesis dynamics is performed by considering a FBR with a suspended E. coli cell culture, operated with the initial/nominal conditions given in Table 1. It is worth mentioning that the bioreactor includes an excess of sparged air and necessary nutrients for a balanced growth of the cell culture. This FBR was used by Chassagnole et al. (2002) to develop experimental kinetic studies to validate their CCM model (denoted here by CHASSM). The same experiments have also been used by Maria (2014) to validate his reduced kinetic model of glycolysis (denoted here by MGM).”

TABLE 1

Table 1. The nominal operating conditions of Chassagnole et al. (2002) and Maria et al. (2018a) for the FBR with suspended E. coli cell culture used to simulate the glycolytic and TRP synthesis processes.

The adopted FBR model of Table 2 “is a classical one (Moser, 1988), developed with the following main hypotheses: (i) the operation is isothermal, iso-pH, and iso-DO (dissolved oxygen); (ii) nutrients are added to the FBR, in recommended quantities, together with an aeration in excess for ensuring an optimal biomass maintenance; (iii) a perfectly mixed liquid phase (with no concentration gradients). The mass balance equations account for the main species in the FBR bulk and of the cellular ones referring to the glycolysis and TRP synthesis dynamics. The resulting ODE model of the FBR (Table 2) also includes the process kinetics (species dynamics). To obtain the species time trajectories with a satisfactory accuracy, by using Table 1 (or other) initial conditions, a low-order stiff integrator (ODE23S) of the Matlab™ software was used” (Maria et al., 2018b,c).

TABLE 2

Table 2. Species mass balance in the FBR model describing the dynamics of the cellular glycolysis species according to the MGM kinetic model of Maria (2014).

Dynamic Models for the Oscillatory Glycolysis and for the Oscillatory TRP Synthesis in the E. coli Cells

Glycolysis Model

“Glycolysis is a sequence of enzymatic reactions (Figures 1, 2) that converts glucose (GLC) into pyruvate (PYR). The free energy released by the subsequent TCA originating from PYR is used to form the high-energy molecules ATP, and NADH that support the glycolysis and the other enzymatic reactions into the cell (Mathews et al., 1999). Consequently, an adequate modeling/simulation of the glycolysis kinetics is of high importance because its intermediates are entry/exit points to/from glycolysis. For instance, most of the monosaccharides, such as fructose or galactose, can be converted to one of these intermediates. In turn, glycolytic intermediates are directly used in subsequent metabolic pathways. For example, DHAP (an intermediate in the F6P conversion to G3P in Figure 1) is a source of the glycerol that combines with fatty acids to form fat. In addition, NADPH is also formed by the PPP, which converts GLC into R5P, which is used in the synthesis of nucleotides and nucleic acids. PEP is, as well, the starting point for the synthesis of essential amino acids (AA) such as TRP, cysteine, arginine, serine, etc. (Calhoun and Swartz, 2006; Maria et al., 2018a).”

“Due to the huge importance of the glycolysis in simulating the CCM dynamics, intense efforts have been invested both in the experimental study and in modeling of its dynamics in various bacteria (Reeves and Sols, 1973; Bennett et al., 2009; Flamholz et al., 2013; Alberton et al., 2015).”

“However, modeling in detail the glycolysis kinetics and its regulation is a difficult task due to its high complexity. Despite these difficulties, a large number of extended/lumped kinetic models have been proposed (some of them being mentioned in Outline 1) of a complexity varying in the range of 18–30 key species and 48–52 key reactions, with a total of 24–300 or more rate constants. Most of these models are however too complex to be easy to use and identified. Besides, their adequacy is not always satisfactory. Thus, with few exceptions, most of the mentioned models cannot satisfactorily simulate the glycolytic oscillations on a deterministic basis” (Maria, 2014; Mihalachi and Maria, 2019).

TABLE OUTLINE 1

TABLE Outline 1. Some dynamic models of glycolysis from the literature.

Starting from the reaction pathway of Figure 1, from the CHASSM and other kinetic models (Outline 1), and by applying advanced lumping algorithms belonging to the physical chemistry (based on species lumping rules with preserving the chemical reaction invariants) (see details of Maria, 2004, 2005), Maria (2014) obtained a reduced kinetic model (that is the MGM) of the glycolysis. The MGM (adopted here) “is accounting for 9 key species, 7 lumped reactions, and includes 17 rate constants (Table 2). Its parameters have been estimated by using the experimental kinetic data of Chassagnole et al. (2002) and Visser et al. (2004).”

The MGM model proved that “it can satisfactorily simulate the dynamics of the glycolytic species concentrations (steady-state QSS, oscillatory, or transient) according to various internal/external regimes, related to (i) the GLC concentration level/dynamics in the bioreactor, (ii) the cell total energy resources in A(MDT)P, and (iii) the cell phenotype responsible for activity of the enzymes involved in the ATP utilization/recovery system. The MGM has been inserted in the bioreactor model template (Table 2) when simulating the dynamics of the [GLC] in the liquid phase simultaneously with that of the cell metabolites. A direct connection between the macro-scale (bioreactor bulk-phase) and the nano-scale (cellular) process variables is thus realized.”

According to Franck (1980), “oscillations in chemical systems represent periodic transitions in time of species concentrations. Thus, spontaneous occurrence of self-sustained oscillations in chemical systems is due 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 (see the discussion of Maria et al., 2018c on the glycolytic oscillation occurrence). 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 kinetics prerequisites for the occurrence of spontaneous oscillations, as well as their consequences, extensively discussed by Franck (1980) and Maria et al. (2018c).”

“In the glycolysis case, oscillations is due to the antagonistic action of two processes on regulating the V2 reaction rate (i.e., the oscillation node; Termonia and Ross, 1981a,b, 1982; Maria et al., 2018c). The V2 reaction converts F6P in FDP (see the pink rectangle of Figure 1, including a lumped representation of the glycolytic oscillations node). Glycolytic oscillations properties (period, amplitude) are determined by both external and internal (phenotype) factors. According to Maria (2014) and Maria et al. (2018a,b,c), the glycolysis dynamics [quasi-steady state (QSS) or oscillatory] depends on several factors:”

i) The glucose level in the liquid-phase {denoted by [GLC]ex}, which varies according to the FBR operating conditions;

ii) “The efficiency and the dynamics of the whole ATP recovery system. Among the involved parameters, an essential one is the k6 rate constant (related to the ATP-ase characteristics in Figure 2). The involved enzyme characteristics are determined by the cell phenotype (genom) controlling the total energy resources. To not complicate our simulations, the [AMDTP] level was kept unchanged in the present analysis at the value given in Table 1” (Maria, 2004).

iii) “As an important remark, the glycolysis is a systemic process, with a complex regulatory structure. Consequently, oscillations are also related to the rate constants of all the involved reactions, and their appropriate ratios (depending on the enzymes' activity of each microorganism)” (Maria et al., 2018c).

TRP Synthesis Model

“TRP is an aromatic non-polar α-amino-acid essential in humans that is used in the cell biosynthesis of proteins, being also a precursor to the neuro-transmitter serotonin, of the melatonin hormone, and of vitamin PP. Therefore, maximizing its production via model-based analyses is of particular industrial interest” (Slominski et al., 2002).

“The TRP operon of E. coli is one of the most extensively studied molecular regulatory systems (Hernandez-Valdez et al., 2010). The TRP synthesis is known as being an oscillatory process. However, due to the process high complexity, only reduced dynamic models involving lumped reactions/species are used, the regulatory performance being included in adjustable model terms and rate constants. For this reason, in the present analysis, the in silico simulation of the TRP synthesis was performed by using the lumped kinetic model of Maria et al. (2018a,b).

This kinetic model is based on the simplified TRP synthesis scheme displayed in Figure 2, derived from various studies reviewed by Maria et al. (2018a). The adopted model for the TRP synthesis, presented in Table 3, is a modification of the Bhartiya et al. (2006) model in order to better fit the experimental kinetic curves of the key species {OR, mRNA, T, E}. Besides, the model was explicitly connected to the glycolysis (as displayed in Figure 2), by including in the TRP mass-balance [i.e., (dc_T/dt) term in Table 3] a term accounting for the PEP precursor, while the PEP consumption term is included in the PEP balance of the MGM model (Table 2). Other dynamic models for the TRP synthesis module are reviewed by Maria et al. (2018a,b).”

TABLE 3

Table 3. The TRP synthesis kinetic model of Maria et al. (2018a) modified to be coupled with the glycolysis model.

Results and Discussion

Glycolytic Oscillations

Repeated simulations of the bioreactor dynamics using the FBR/MGM kinetic models, with the initial conditions of Table 1 and the parametric ranges of [GLC]ex ∈ [0.01–1.5] mM (at t = 0), and k6 ∈ [10–⁵-20] min⁻¹ lead to the following results (Maria et al., 2018a,b,c):

i) “Several glycolytic stationary oscillations domains exist in the E. coli cells, as indicated by the thick lines of Figure 3F plotted in the {[GLC]ex vs. k6} coordinates.

ii) As displayed in Figure 3F, glycolytic stationary oscillations occur for a slow GLC import not only due to a low [GLC]ex level in the environment but also due to small k6 constant values (that correspond to a low recovery rate of the ATP). Conversely, higher concentrations of GLC in the bioreactor will trigger higher GLC import rates. In this case, glycolytic oscillations are also possible if the k6 constant reported large values (for a certain K constant controlling the AMDTP pathway/equilibrium given in the Table 2). However, the ATP recovery rate is limited by the AMDTP resources and by the interconversion balance of the AMDTP system (Figures 1, 2 and Table 2). As reported by Maria et al. (2018a,b,c), in the cells with too small, or too large k6 values, the glycolysis often reaches its (non-oscillatory) steady state.”

iii) The “glycolytic oscillation domains plotted in Figure 3F are very narrow. Such a result reflects their high sensitivity vs. lots of external and internal factors. Besides, oscillations present a poor stability vs. internal/external factors, as proved by the plotted limit cycles (omitted here; see Maria et al., 2018a,b,c). Experiments in the literature have found that this stability is dependent on the metabolism characteristics of every microorganism. For instance, by contrast, the glycolytic oscillations in yeast have been proved to be very robust even in the presence of environmental noise, with oscillations being a side effect of the trade-offs between robustness and regulatory efficiency of the feedback control of the autocatalytic reaction network (Chandra et al., 2011; Gehrmann et al., 2011).”

iv) The numerical analysis results also indicated that “larger values of k6 lead to a slight decrease in the oscillation period and, eventually, the oscillation disappearance. This is due to the quick consumption of GLC by the cells following a more rapid ATP recovery system (Maria et al., 2018c).”

v) FBR dynamic simulations “have identified glycolytic oscillations with a period of 0.4–1 min, depending on the k6 value and on the [Glc]ex level. For comparison, various experiments in the literature have reported periods in a large range, that is: 0.2 min (Madsen et al., 2005), 2–100 s (Westermark and Lansner, 2003), 15 s (Silva and Yunes, 2006), 1–20 min (Bier et al., 1996), up to 3 h (Rapp, 1979), or 0.2 min to h (Diaz Ricci, 2000).”

vi) The simulated glycolytic oscillations of Figure 3 (FDP and F6P species) are similar to the experimentally recorded dynamics by Schaefer et al. (1999) and Chiarugi et al. (2006) and also similar to the dynamic simulations of Sel'kov (1968), Bier et al. (1996), Elias (2010), and de la Fuente (2010). In fact, Figures 3A–E display an incipient phase of the oscillation occurrence, when the species oscillation amplitude grows. However, over a longer time domain (not shown here), the oscillations stabilize and become stationary.

vii) “The simulated [GLC]ex dynamics in the FBR proved that, for a relatively high [GLC] = 200 mM in the feed, and for all the abovementioned ranges of internal/external operating conditions, the bioreactor evolution is always toward a steady state (QSS), with a faster or slower rate depending on the initial [GLC] in the bioreactor, irrespectively to the cell metabolism (stationary/homeostatic, or unbalanced) (Maria et al., 2018a,b,c).

viii) The factors influencing the glycolysis dynamics mentioned at the end of section Glycolysis Model are confirmed to have a major influence on the glycolysis dynamics as proved by the present analysis.”

FIGURE 3

Figure 3. (A–E) Simulated glycolytic stationary oscillations of the main glycolytic metabolites (PYR, F6P, FDP, ATP, and PEP) in E. coli for the bioreactor nominal operating conditions of Table 1 ([AMDTP] = 5.82 mmol/L, D = 0.001667·1/min), [GLC]ex = 0.0557 mM (at t = 0), and k6 = 12 1/min. The simulated FBR running time is 10 min. (F) Glycolytic stationary oscillation domains (thick lines) in E. coli plotted in the plane [Glc]ext (at t = 0), and k6, for the bioreactor nominal operating conditions of Table 1. The red point corresponds to the cell species dynamics plotted in (A–E). [Glc]ex, glucose concentration in the cell environment (bulk phase). Panel (F) was adapted from Maria et al. (2018c) courtesy of CABEQ Jl.

TRP Synthesis Oscillations

“Under certain conditions, the TRP synthesis presents oscillations (Bhartiya et al., 2006). Being strongly connected with the glycolysis (via PEP species), it is important to also study the influence of the glycolytic oscillations on the TRP synthesis dynamics. Such an analysis turns out to be of high practical interest in order to adapt the bioreactor operation to maximize the TRP production and, eventually, for the in silico design GMO with such an objective (by the modification of the cell metabolic fluxes) (Maria et al., 2018a,b).

In particular, the glycolysis intermediate PEP is the starting point for the synthesis of essential amino acids, including TRP (Stephanopoulos and Simpson, 1997; Lodish et al., 2000; Calhoun and Swartz, 2006).”

Having PEP as one of the precursors, maximization of TRP synthesis clearly depends on the glycolysis intensity (average levels of glycolytic species) and dynamics (QSS or oscillatory). “On the other hand, as previously discussed, glycolysis is controlled by cell internal and external factors, which indirectly will also influence the TRP synthesis and oscillations, as follows:

(A) The GLC import rate (flux 50 in Figure 1) triggered by the external [GLC]ex and by the PEP and PYR levels into the cell;

(B) The limited ATP energy resources and a slow recovery system can slow down the GLC import, the glycolysis and, implicitly, the all-metabolic syntheses, including the TRP production;

(C) The bacteria genome (cell phenotype) plays an essential role, because it determines the characteristics of the ATP-ase and AK-ase enzymes responsible for the ATP-to-ADP conversion and for the ATP recovery rate during the glycolytic reactions (accounted by means of the K and k6 rate constants in the MGM kinetic model of Table 2). In fact, the A(MDT)P interconversion system affects most of the metabolic reactions.”

(D) “Being a systemic process, inherently the glycolytic oscillations are also related to the rate constants of all the glycolysis reactions. As an example, Silva and Yunes (2006) found that oscillations are only possible if the [GLC]ex and the maximum reaction rates controlled by the PFKase and GKase are within specific intervals. The GKase is one of the enzymes controlling the V1 lumped reaction related to the PTS import system (GLC to G6P and then to F6P in Figure 2 and Table 2). The PFKase controls the V2 reaction (of Figure 2) responsible for the FDP synthesis.”

(E) The results of section Glycolytic Oscillations demonstrated that both external (like [GLC]ex level) and internal/metabolic parameters (like K, k6 constants) governing the glycolytic oscillations will influence the TRP synthesis dynamics and level.

“As has been proven experimentally by Bhartiya et al. (2006), the TRP synthesis is an oscillatory process with a complex engine. Model-based numerical analyses of Maria et al. (2018a,b) highlighted some of the essential factors on which the TRP process dynamics depends.”

As mentioned by Santillan and Mackey (2001); Mackey et al. (2004), and Hernandez-Valdez et al. (2010), “oscillations in the TRP synthesis are produced due to the concomitant activation and high-order repression of the TRP-operon expression, together with a non-linear demand for the end product, making its expression cyclic. Maria et al. (2018a,b) pointed out through model-based simulations that the cell dilution rate (related to the cell cycle), adjusted to be consistent with the liquid residence time in the bioreactor, also exerts a strong influence on the TRP system dynamics.”

The numerical analysis of the present paper has been performed by using the FBR reactor model, which includes the coupled glycolysis MGM model (section Glycolysis Model; Table 2) and the TRP synthesis kinetic model (section TRP Synthesis Model; Table 3). “The sensitivity analysis of the TRP production was performed by considering some of the most influential parameters checked in the range of [GLC]ex ∈ [0.01–1.5] mM (at t = 0); k6 ∈ [10⁻⁵-20] min⁻¹; D ∈ [10⁻⁴-10⁻²] min⁻¹; and initial [GLC]ex ∈ [0.005–5] mM (at t = 0). The simulation results for only two relevant operating conditions have been plotted in Figures 4A,B. This analysis leads to several results and conclusions:

i) The [GLC]ex, the constant k6, and the bioreactor dilution D (considered equal to the cell dilution) exert the highest influence not only on the glycolysis dynamics but also on the TRP synthesis dynamics and production (due to its close link to the glycolysis via PEP). Thus, under the initial FBR conditions of Table 1, for a low FBR dilution rate (D), and for conditions leading to a QSS glycolysis, the TRP synthesis also displays a stationary evolution (Figure 4A). By contrast, at higher dilutions, and when glycolysis meets the conditions necessary for an oscillatory process (of Figure 3F), the TRP synthesis also presents an oscillatory dynamics (Figure 4B). Consequently, the bioreactor dilution presents a strong influence on the QSS or oscillatory regime of the linked glycolysis and TRP synthesis. The TRP production (see its definition in the caption from Figures 4A,B) is influenced accordingly. A value of k6 = 12 min⁻¹ was considered in all the tested cases here.

ii) While glycolysis exerts a strong influence on the TRP synthesis dynamics, as proved by Figures 4A,B, the reverse influence is minor, as proved by disconnected glycolysis simulations (not reproduced here).

iii) For the high feeding rates (Table 1), the initial [GLC] in the FBR bioreactor does not quantitatively influence the TRP bioreactor performances (Maria et al., 2018b).

vi) Simulations of the only TRP synthesis, disconnected from the glycolytic process, but with employing various [PEP] average levels (Mihalachi and Maria, 2019), indicated that PEP average level has a huge influence on the dynamics and concentrations of the TRP synthesis species. In turn, [PEP] is controlled by the glycolysis dynamics, which, in turn, is controlled by the abovementioned external and internal factors.

v) It clearly appears that, beside cell phenotype (governing the TRP operon expression), glycolysis is one of the major factors influencing the TRP production. Thus, by ranging the FBR operating parameters, the TRP production can be maximized (Maria et al., 2018b).

vi) In all the FBR operating cases checked by Maria et al. (2018b) with the initial conditions of Table 1, but in the range of D ∈ [10⁻⁴-0.01] min⁻¹, simulations demonstrated that [GLC]ex in the liquid bulk always evolves toward its steady state irrespectively of the stationary or oscillatory dynamics of the cell metabolic processes.”

FIGURE 4

Figure 4. (A) Simulated dynamics of the coupled glycolysis and TRP synthesis main species in E. coli cells under the bioreactor nominal conditions of Table 1, for a cell phenotype with k6 = 12 1/min. The following operating conditions have been checked:
D = 0.0001 (1/min) (cell/bioreactor dilutions); [GLC]ex = 1 mM (at t = 0), when both glycolysis [left] and TRP synthesis [right] display a stationary behavior (quasi-steady-state, QSS). The realized TRP production over the batch time is 0.39 μM/min.
TRP production (μM/min) = {Dilution rate (F_L/V_L)} × {max [TRP] (t)}
Species dynamics was generated by using the coupled bioreactor/glycolysis/TRP models of Tables 2, 3. (B) Simulated dynamics of the coupled glycolysis and TRP synthesis main species in E. coli cells under the bioreactor nominal conditions of Table 1, for a cell phenotype with k6 = 12 1/min. The following operating conditions have been checked:
D = 0.01 (1/min) (cell/bioreactor dilutions); [GLC]ex = 0.0557 mM (at t = 0), when both glycolysis [left] and TRP synthesis [right] display stationary oscillations. The realized TRP production over the batch time is 0.1085 μM/min.
TRP production (μM/min) = {Dilution rate (F_L/V_L)} × {max [TRP] (t)}
The species dynamics were generated by using the coupled bioreactor/glycolysis/TRP models of Tables 2, 3.

Engineering Implications—Sensitivity Analysis

From an engineering point of view, the study suggests how to modulate the most influential factors, that is: [GLC]ex {via [GLC](feed), and [GLC](initial)}, D = μ, k6, and others (via cell phenotype), before deriving an optimal FBR operating policy (for instance, a timestep-wise GLC feeding policy) leading to maximization of the TRP synthesis.

An exhaustive or an adaptive model-based search can identify the FBR optimal operating policy that corresponds to a maximum of TRP production (an analysis not developed in this work). In the present study, a brief sensitivity analysis based on the above model has been performed, thus preceding the FBR future optimization. Numerical simulations revealed several interesting conclusions:

a.- FBR simulations using various operating parameters proved that FBR efficiency (TRP production) is not influenced by [GLC](initial) < 100 mM in the bioreactor, once [GLC](feed) > 100 mM.

b.- In all tested cases covering the ranges [GLC](feed) of 100–200 mM, [GLC](initial) of 1–50 mM, and D = 0.00001–0.01 (1/min), the FBR rapidly evolves toward its steady state, corresponding to a usually small [GLC](stationary) (below 1 mM), which proves the FBR efficiency. An example of GLC dynamics in the bioreactor is displayed in Figure 5A.

c.- The TRP productivity increases with D, as plotted in Figure 5C.

d.- The TRP productivity also increases with the [GLC]ex, as plotted in Figure 5D, where the external (bulk) concentration is given by the steady-state level of [GLC](stationary). The combined dependency of TRP productivity on the reactor dilution and the [GLC](stationary) is given in Figure 5B, confirming conclusions (c) and (d) for a wide range of the reactor dilutions.

FIGURE 5

Figure 5. The dynamics and sensitivity analysis of the FBR (simulated results). (A) [GLC] dynamics in the FBR for [GLC](feed) = 200 mM; [GLC](initial) = 5 mM; D = 0.001 (1/min). The same behavior for D = 0.00001–0.001 (1/min). (B) TRP production (μM/min) function of dilution and the [GLC](stationary) at the steady state. Fixed parameters of [GLC](feed) = 200 mM; [GLC](initial) = 5 mM. (C) TRP production (μM/min) function of reactor dilution. Fixed parameters: of [GLC](feed) = 200 mM; [GLC](initial) = 5 mM. (D) TRP production (μM/min) function of [GLC](stationary) at the steady state, generated for dilutions in the range of 0.00001–0.01 (1/min). Fixed parameters of [GLC](feed) = 200 mM; [GLC](initial) = 5 mM.

As proved in this work, the TRP productivity is also strongly dependent on the oscillatory characteristics of the glycolysis, determined by the abovementioned operating parameters, and on the activity of enzymes involved in the ATP recovery system (that is, the rate constants k6, K, and [AMDTP] of Table 2). Thus, from a biological point of view, as mentioned by Silva and Yunes (2006), “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 toward other cell cycle activities such as cell division. Consequently, GMO with modified enzyme activity (especially those related to the ATP use/recovery system of Figure 2) will lead to noticeable modifications in the metabolic species dynamics and concentrations.”

Conclusions

The in silico analysis of this paper demonstrates in a meaningful and relevant way the importance of using a detailed enough and adequate structured dynamic model linking the metabolic cellular processes and the bioreactor state variables for engineering purposes. Such a modular hybrid model can link the macro-scale (bioreactor liquid phase) process variables to the nano-scale (cellular) ones.

The structured model not only can be used for future FBR optimization, “but can also be a valuable tool to evaluate the cellular metabolic fluxes (i.e., the homeostatic metabolic reaction rates, not evaluated here), thus opening the possibility to in silico re-design the cell metabolism to obtain GMOs with industrial or medical applications (Hatzimanikatis et al., 1996; Stephanopoulos et al., 1998; Visser et al., 2004; Styczynski and Stephanopoulos, 2005; Maria, 2017a,b, 2018).

Thus, the large experimental and computational effort to validate such structured cell models is eventually fully justified through the practical advantages offered by such an engineering analysis.

It is also to be emphasized that such a modular and structured approach of the dynamic cellular models offers the possibility to study the interference of the CCM sub-process (e.g., glycolysis and AA synthesis here), together with the influence of the external conditions. Such a modular simulation platform presents the advantage to be easy to extend by the inclusion of new CCM modules.”

Data Availability Statement

All datasets generated for this study are included in the article/supplementary material.

Author Contributions

The author confirms being the sole contributor of this work and has approved it for publication.

Conflict of Interest

The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

References

Alberton, K. P. F., Alberton, A. L., Di Maggio, J. A., Estrada, V. G., Díaz, M. S., and Secchi, A. R. (2015). Simultaneous parameters identifiability and estimation of an E. coli metabolic network model. BioMed. Res. Int. 2015:454765. doi: 10.1155/2015/454765

CrossRef Full Text | Google Scholar

Bennett, B. D., Kimball, E. H., Gao, M., Osterhout, R., Van Dien, S. J., and Rabinowitz, J. D. (2009). Absolute metabolite concentrations and implied enzyme active site occupancy in Escherichia coli. Nat. Chem. Biol. 5, 593–599. doi: 10.1038/nchembio.186

CrossRef Full Text | Google Scholar

Bhartiya, S., Chaudhary, N., Venkatesh, K. V., and Doyle, F. J. (2006). Multiple feedback loop design in the tryptophan regulatory network of E. coli suggests a paradigm for robust regulation of processes in series. J. R. Soc. Interface 3, 383–391. doi: 10.1098/rsif.2005.0103

PubMed Abstract | CrossRef Full Text | Google Scholar

Bier, M., Teusink, B., Kholodenko, B. N., and Westerhoff, H. V. (1996). Control analysis of glycolytic oscillations. Biophys. Chem. 62, 15–24. doi: 10.1016/S0301-4622(96)02195-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Buchholz, A., Hurlebaus, J., Wandrey, C., and Takors, R. (2002). Metabolomics: quantification of intracellular metabolite dynamics. Biomol. Eng. 19, 5–15. doi: 10.1016/S1389-0344(02)00003-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Calhoun, K. A., and Swartz, J. R. (2006). Total amino acid stabilization during cell-free protein synthesis reactions. J. Biotechnol. 123, 193–203. doi: 10.1016/j.jbiotec.2005.11.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Chandra, F. A., Buzi, G., and Doyle, J. C. (2011). Glycolytic oscillations and limits on robust efficiency. Science. 333, 187–192. doi: 10.1126/science.1200705

PubMed Abstract | CrossRef Full Text | Google Scholar

Chassagnole, C., Noisommit-Rizzi, N., Schmid, J. W., Mauch, K., and Reuss, M. (2002). Dynamic modeling of the central carbon metabolism of Escherichia coli. Biotechnol. Bioeng. 79, 53–73. doi: 10.1002/bit.10288

PubMed Abstract | CrossRef Full Text | Google Scholar

Chiarugi, D., Chinellato, M., Degano, P., Lo Brutto, G., and Marangoni, R. (2006). “Feedbacks and oscillations in the virtual cell VICE,' in Proceedings of International Conference on Computational Methods in Systems Biology (Trento: CMSB), 93–107. doi: 10.1007/11885191_7

CrossRef Full Text | Google Scholar

Costa, R. S., Machado, D., Rocha, I., and Ferreira, E. C. (2010). Hybrid dynamic modeling of Escherichia coli central metabolic network combining Michaelis–Menten and approximate kinetic equations. BioSystems 100, 150–157. doi: 10.1016/j.biosystems.2010.03.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Costa, R. S., Rocha, I., and Ferreira, E. C. (2008). Model reduction based on dynamic sensitivity analysis: A systems biology case of study (Ph. D. grant report). University of Minho, Braga, Portugal.

Google Scholar

de la Fuente, I. M. (2010). Quantitative analysis of cellular metabolic dissipative, self-organized structures. Int. J. Mol. Sci. 11, 3540–3599. doi: 10.3390/ijms11093540

PubMed Abstract | CrossRef Full Text | Google Scholar

Degenring, D., Froemel, C., Dikta, G., and Takors, R. J. (2004). Sensitivity analysis for the reduction of complex metabolism models. J. Proc. Control 14, 729–745. doi: 10.1016/j.jprocont.2003.12.008

CrossRef Full Text | Google Scholar

Diaz Ricci, J. C. (2000). ADP modulates the dynamic behavior of the glycolytic pathway of Escherichia coli. Biochem. Biophys. Res. Commun. 271, 244–249. doi: 10.1006/bbrc.2000.2603

PubMed Abstract | CrossRef Full Text | Google Scholar

Edwards, J. S., and Palsson, B. O. (2000). The Escherichia coli MG1655 in silico metabolic genotype: its definition, characteristics, and capabilities. Proc. Natl. Acad. Sci. U.S.A. 97, 5528–5533. doi: 10.1073/pnas.97.10.5528

PubMed Abstract | CrossRef Full Text | Google Scholar

Elias, A. C. (2010). Escherichia coli: Dynamic Analysis of the Glycolytic Pathway (Ph.D. thesis). Facultad de Bioquímica y Farmacia, Universidad Nacional de Tucumán, San Miguel de Tucumán, Argentina.

Flamholz, A., Noor, E., Bar-Even, A., Liebermeister, W., and Milo, R. (2013). Glycolytic strategy as a tradeoff between energy yield and protein cost. Proc. Natl. Acad. Sci. U.S.A. 110, 10039–10044. doi: 10.1073/pnas.1215283110

PubMed Abstract | CrossRef Full Text | Google Scholar

Franck, U. F. (1980). Feedback kinetics in physicochemical oscillators. Ber. Bunsenges. Phys. Chem. 84, 334–341. doi: 10.1002/bbpc.19800840407

CrossRef Full Text | Google Scholar

Gehrmann, E., Glässer, C., Jin, Y., Sendhoff, B., Drossel, B., and Hamacher, K. (2011). Robustness of glycolysis in yeast to internal and external noise. Phys. Rev. E 84:021913. doi: 10.1103/PhysRevE.84.021913

PubMed Abstract | CrossRef Full Text | Google Scholar

Hatzimanikatis, V., and Bailey, J. E. (1997). Studies on glycolysis- I. Multiple steady states in bacterial. Chem. Eng. Sci. 52, 2579–2588. doi: 10.1016/S0009-2509(97)00074-2

CrossRef Full Text

Hatzimanikatis, V., Floudas, C. A., and Bailey, J. E. (1996). Analysis and design of metabolic reaction networks via mixed-integer linear optimization. AIChE J. 42, 1277–1292. doi: 10.1002/aic.690420509

CrossRef Full Text | Google Scholar

Hempel, D. C. (2006). Development of biotechnological processes by integrating genetic and engineering methods. Eng. Life Sci. 6, 443–447. doi: 10.1002/elsc.200620149

CrossRef Full Text | Google Scholar

Hernandez-Valdez, A., Santillan, M., and Zeron, E. S. (2010). Cycling expression and cooperative operator interaction in the trp operon of Escherichia coli. J. Theor. Biol. 263, 340–352. doi: 10.1016/j.jtbi.2009.12.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Kadir, T. A., Mannan, A. A., Kierzek, A. M., McFadden, J., and Shimizu, K. (2010). Modeling and simulation of the main metabolism in Escherichia coli and its several single-gene knockout mutants with experimental verification. Microb. Cell Fact. 9:88. doi: 10.1186/1475-2859-9-88

PubMed Abstract | CrossRef Full Text | Google Scholar

KEGG (2011). Kyoto Encyclopedia of Genes and Genomes, Kanehisa Laboratories, Bioinformatics Center of Kyoto University. Available online at: http://www.genome.jp/kegg/pathway.html (accessed September 10, 2020).

Lodish, H., Berk, A., Zipursky, S. L., Matsudaira, P., Baltimore, D., and Darnell, J. E. (2000). Molecular Cell Biology. New York, NY: W. H. Freeman & Co. Available online at: http://2009.igem.org/wiki/images/a/a7/The_tryptophan_biosynthetic_pathway.pdf (accessed May 30, 2009).

Mackey, M. C., Santillán, M., and Yildirim, N. (2004). Modeling operon dynamics: thetryptophan and lactose operons as paradigms. C. R. Biologies. 327, 211–224. doi: 10.1016/j.crvi.2003.11.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Madsen, M. F., Danø, S., and Sørensen, P. G. (2005). On the mechanisms of glycolytic oscillations in yeast. FEBS J. 272, 2648–2660. doi: 10.1111/j.1742-4658.2005.04639.x

CrossRef Full Text | Google Scholar

Maria, G. (2004). A review of algorithms and trends in kinetic model identification for chemical and biochemical systems. Chem. Biochem. Eng. Q. 18, 195–222.

Google Scholar

Maria, G. (2005). Relations between apparent and intrinsic kinetics of programmable drug release in human plasma. Chem. Eng. Sci. 60, 1709–1723. doi: 10.1016/j.ces.2004.11.009

CrossRef Full Text | Google Scholar

Maria, G. (2014). In silico derivation of a reduced kinetic model for stationary or oscillating glycolysis in Escherichia coli bacterium. Chem. Biochem. Eng. Q. 28, 509–529. doi: 10.15255/CABEQ.2014.2002

CrossRef Full Text | Google Scholar

Maria, G. (2017a). Deterministic Modelling Approach of Metabolic Processes in Living Cells - A Still Powerful Tool for Representing the Metabolic Process Dynamics. Newbury Park, CA: Juniper Publishers. Available online at: https://juniperpublishers.com/ebook-info.php (accessed June 15, 2017).

Maria, G. (2017b). A Review of Some Novel Concepts Applied to Modular Modelling of Genetic Regulatory Circuits. Newbury Park, CA: Juniper Publishers. Available online at: https://juniperpublishers.com/ebook-info.php (accessed March 15, 2017).

Maria, G. (2018). In-silico Design of Genetic Modified Micro-Organisms (GMO) of Industrial Use, by Using Systems Biology and (Bio)Chemical Engineering Tools. Simi Valley, CA: Juniper Publishers. Available online at: https://juniperpublishers.com/ebook-info.php (accessed March 21, 2018).

Maria, G., Gijiu, C. L., Maria, C., and Tociu, C. (2018a). Interference of the oscillating glycolysis with the oscillating tryptophan synthesis in the E. coli cells. Comp. Chem. Eng. 108, 395–407. doi: 10.1016/j.compchemeng.2017.10.003

CrossRef Full Text | Google Scholar

Maria, G., and Luta, I. (2013). Structured cell simulator coupled with a fluidized bed bioreactor model to predict the adaptive mercury uptake by E. coli cells. Comp. Chem. Eng. 58, 98–115. doi: 10.1016/j.compchemeng.2013.06.004

CrossRef Full Text | Google Scholar

Maria, G., Mihalachi, M., and Gijiu, C. L. (2018b). In silico optimization of a bioreactor with an E. coli culture for tryptophan production by using a structured model coupling the oscillating glycolysis and tryptophan synthesis. Chem. Eng. Res. Design. 135, 207–221. doi: 10.1016/j.cherd.2018.05.011

CrossRef Full Text | Google Scholar

Maria, G., Mihalachi, M., and Gijiu, C. L. (2018c). Model-based identification of some conditions leading to glycolytic oscillations in E. coli cells. Chem. Biochem. Eng. Q. 32, 523–533. doi: 10.15255/CABEQ.2017.1241

CrossRef Full Text | Google Scholar

Maria, G., Xu, Z., and Sun, J. (2011). Multi-objectiveMINLP optimization used to identify theoretical gene knockout strategies for E. coli cell. Chem. Biochem. Eng. Q. 25, 403–424. .

Google Scholar

Mathews, C. K., van Holde, K. E., and Ahem, K. G. (1999). Biochemistry. New Jersey, Prentice Hall.

Google Scholar

Mihalachi, M., and Maria, G. (2019). Influence of pep glycolytic precursor on tryptophan synthesis dynamics in E. coli cells. UPB. Sci. Bull. Series B 81, 29–36.

Moser, A. (1988). Bioprocess Technology - Kinetics and Reactors. Berlin: Springer Verlag. doi: 10.1007/978-1-4613-8748-0

CrossRef Full Text | Google Scholar

Palsson, B. O. (2005). Systems Biology – Properties of Reconstructed Networks. Cambridge: Cambridge University Press.

Google Scholar

Peskov, K., Mogilevskaya, E., and Demin, O. (2012). Kinetic modelling of central carbon metabolism in Escherichia coli. FEBS J. 279, 3374–3385. doi: 10.1111/j.1742-4658.2012.08719.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Rapp, P. E. (1979). An atlas of cellular oscillators. J. Exp. Biol. 81, 281–306.

PubMed Abstract | Google Scholar

Reeves, R. E., and Sols, A. (1973). Regulation of Escherichia coli phosphofructokinase in situ. Biochem. Biophys. Res. Commun. 50, 459–466. doi: 10.1016/0006-291X(73)90862-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Santillan, M., and Mackey, M. C. (2001). Dynamic behavior in mathematical models of the tryptophan operon. Chaos 11, 261–268. doi: 10.1063/1.1336806

PubMed Abstract | CrossRef Full Text | Google Scholar

Schaefer, U., Boos, W., Takors, R., and Weuster-Botz, D. (1999). Automated sampling device for monitoring intracellular metabolite dynamics. Analy. Biochem. 270, 88–96. doi: 10.1006/abio.1999.4048

PubMed Abstract | CrossRef Full Text | Google Scholar

Sel'kov, E. E. (1968). Self-oscillations in glycolysis. 1. a simple kinetic model. Eur J. Biochem. 4, 79–86. doi: 10.1111/j.1432-1033.1968.tb00175.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Silva, A. S., and Yunes, J. A. (2006). Conservation of glycolytic oscillations in Saccharomyces cerevisiae. Genet. Mol. Res. 5, 525–535.

PubMed Abstract | Google Scholar

Slominski, A., Semak, I., Pisarchik, A., Sweatman, T., Szczesniewski, A., and Wortsman, J. (2002). Conversion of L-tryptophan to serotonin and melatonin in human melanoma cells. FEBS Lett. 511, 102–106. doi: 10.1016/S0014-5793(01)03319-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Stephanopoulos, G., and Simpson, T. W. (1997). Flux amplification in complex metabolic networks. Chem. Eng. Sci. 52, 2607–2627. doi: 10.1016/S0009-2509(97)00077-8

CrossRef Full Text | Google Scholar

Stephanopoulos, G. N., Aristidou, A. A., and Nielsen, J. (1998). Metabolic Engineering. Principles and Methodologies. San Diego, CA: Academic Press.

Google Scholar

Styczynski, M. P., and Stephanopoulos, G. (2005). Overview of computational methods for the inference of gene regulatory networks. Comput. Chem. Eng. 29, 519–534. doi: 10.1016/j.compchemeng.2004.08.029

CrossRef Full Text | Google Scholar

Termonia, Y., and Ross, J. (1981a). Oscillations and control features in glycolysis: numerical analysis of a comprehensive model. Proc. Natl. Acad. Sci. U.S.A. 78, 2952–2956. doi: 10.1073/pnas.78.5.2952

CrossRef Full Text | Google Scholar

Termonia, Y., and Ross, J. (1981b). Oscillations and control features in glycolysis: analysis of resonance effects. Proc. Natl. Acad. Sci. U.S.A. 78, 3563–3566. doi: 10.1073/pnas.78.6.3563

PubMed Abstract | CrossRef Full Text | Google Scholar

Termonia, Y., and Ross, J. (1982). Entrainment and resonance in glycolysis. Proc. Natl. Acad. Sci. U.S.A. 79, 2878–2881. doi: 10.1073/pnas.79.9.2878

PubMed Abstract | CrossRef Full Text | Google Scholar

Visser, D., Schmid, J. W., Mauch, K., Reuss, M., and Heijnen, J. J. (2004). Optimal re-design of primarymetabolism in Escherichia coli using linlog kinetics. Metab. Eng. 6, 378–390. doi: 10.1016/j.ymben.2004.07.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Westermark, P. O., and Lansner, A. (2003). A model of phosphofructokinase and glycolytic oscillations in the pancreatic b-cell. Biophys. J. 85, 126–139. doi: 10.1016/S0006-3495(03)74460-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Wierschem, K., and Bertram, R. (2004). Complex bursting in pancreatic islets: a potential glycolytic mechanism. J. Theor. Biol. 228, 513–521. doi: 10.1016/j.jtbi.2004.02.022

PubMed Abstract | CrossRef Full Text | Google Scholar

Notations

Abbreviations