Tryptophan Production Maximization in a Fed-Batch Bioreactor with Modified E. coli Cells, by Optimizing Its Operating Policy Based on an Extended Structured Cell Kinetic Model

Hybrid kinetic models, linking structured cell metabolic processes to the dynamics of macroscopic variables of the bioreactor, are more and more used in engineering evaluations to derive more precise predictions of the process dynamics under variable operating conditions. Depending on the cell model complexity, such a math tool can be used to evaluate the metabolic fluxes in relation to the bioreactor operating conditions, thus suggesting ways to genetically modify the microorganism for certain purposes. Even if development of such an extended dynamic model requires more experimental and computational efforts, its use is advantageous. The approached probative example refers to a model simulating the dynamics of nanoscale variables from several pathways of the central carbon metabolism (CCM) of Escherichia coli cells, linked to the macroscopic state variables of a fed-batch bioreactor (FBR) used for the tryptophan (TRP) production. The used E. coli strain was modified to replace the PTS system for glucose (GLC) uptake with a more efficient one. The study presents multiple elements of novelty: (i) the experimentally validated modular model itself, and (ii) its efficiency in computationally deriving an optimal operation policy of the FBR.


Introduction
Over the last few decades, there has been a continuous trend to develop more and more effective bioreactors [1,2] "to industrialize important biosyntheses for producing fine chemicals used in the food, pharmaceutical, or detergent industry, by using freesuspended or immobilized cell cultures (or enzymes) in suitable bioreactors (or enzymatic reactors)", as reviewed by Maria [3]. The batch (BR), semi-batch (fed-batch, FBR), a serial sequence of BRs [4], and the continuously operated fixed-bed or three-phase fluidized-bed bioreactors (with immobilized biocatalyst) are successfully used to conduct biosyntheses aimed at replacing complex chemical and energetically intensive processes, as well as those generating toxic wastes [5,6].
Bioreactors with microbial/animal cell cultures have been developed in simple or complex constructive/operating alternatives as underlined in reviews by Maria [3] or by [9,10]. In spite of their larger volumes, the continuously mixing aerated tank reactors (CSTR), operated in BR or FBR modes, are preferred for bioprocesses requiring a high oxygen transfer and rigorous temperature/pH control. For these reasons, an effective FBR was used in the approached case study of TRP production.  [51] (i.e., the wild cell, including the PTS-system). Fluxes characterizing the membrane transport [Metabolite(e) ↔ Metabolite(c)] and the exchange with environment are omitted from the plot (see [38] for details and explanations regarding the numbered reactions). Notations: [e] = environment; [c] = cytosol. Adapted from Maria et al. [38] with the 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. [38]. The left rectangle indicates the chemical node inducing glycolytic oscillations [33,36]. Notations [+] and [−] denote the feedback positive and negative regulatory loops, respectively. GLC = glucose". See the abbreviation list for species names; V1-V6 = lumped reaction rates discussed in the Section 3.1.3.
At this point, it is worth underlining that the cell metabolism is highly sophisticated, involving 10 3-4 components, 10 3-4 transcription factors (TF-s), activators, inhibitors, and at least one order of magnitude higher number of (bio)chemical reactions, all ensuring a fast adaptation of the cell to the changing environment through complex genetic regulatory circuits (GRC-s) [50]. The cell is highly responsive to the environmental stimuli and highly evolvable by self-changing its genome/proteome and metabolism, mediating the stoichiometry and the reaction rates (fluxes) of the enzymatic reactions to get an optimized and balanced growth by using minimum resources (nutrients/substrates).  [51] (i.e., the wild cell, including the PTS-system). Fluxes characterizing the membrane transport [Metabolite(e) ↔ Metabolite(c)] and the exchange with environment are omitted from the plot (see [38] for details and explanations regarding the numbered reactions). Notations: [e] = environment; [c] = cytosol. Adapted from Maria et al. [38] with the 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. [38]. The left rectangle indicates the chemical node inducing glycolytic oscillations [33,36]. Notations [+] and [−] denote the feedback positive and negative regulatory loops, respectively. GLC = "glucose". See the abbreviation list for species names; V1-V6 = lumped reaction rates discussed in the Section 3.1.3. Development of extended CCM dynamic models on a deterministic basis to adequately simulate in detail the cell metabolism self-regulation, cell growth, and replication for such an astronomical cell metabolism complexity is practically impossible due to the lack of structured and comprehensive information, as well as computational limitations. Reviews of some trials were presented by Styczynski and Stephanopoulos [67] and by Maria [31,32,50].
In spite of such tremendous modeling difficulties, the development of structured reduced deterministic (rather than stochastic) models [31] able to adequately reproduce the dynamics of some CCM complex metabolic syntheses [48,67,68], as well as the dynamics of the genetic regulatory systems [50] tightly controlling the metabolic processes, has reported significant progress over the last few decades [69,70]. Even if they are rather based on sparse information from various sources, unconventional statistical identification, and lumping algorithms [31,41,45,50], such structured reduced deterministic kinetic models have been proven to be extremely useful for in silico analysis and characterization of the CCM, as well as for the design of a novel GRC-s conferring new properties/functions to the mutant cells [31,50,71]. This paper is aimed at proving the feasibility and advantage of using this novel concept to couple an extended cell structured deterministic kinetic model with a bioreactor macroscopic dynamic model. The resulting hybrid dynamic model was successfully used for engineering evaluations. The applied example involves the optimization of the FBR used for TRP synthesis.
"L-tryptophan is a high-value aromatic amino acid with important applications in food and pharma industry. TRP is an aromatic nonpolar α-amino acid essential in humans, which is used in the cell biosynthesis of proteins, being also a precursor to the neurotransmitter serotonin, the melatonin hormone, and vitamin PP" [72].
This paper uses a hybrid dynamic model built up by Maria [35] by linking a CCMbased structured kinetic model with an FBR simple dynamic model. The resulting hybrid FBR model was used to computationally determine the optimal (time stepwise) feeding policy of the FBR used by Chen et al. [73] to study TRP synthesis using a modified E. coli T5 strain culture. The thus obtained optimal operating policy of the FBR has proven to be very effective, by ensuring maximization of TRP production involving a few key control variables (i.e., the feed flow-rate and the feeding GLC concentration), and it reported better performance compared to the non-optimally operated FBR of Maria et al. [34,35] or of Chen [74].
The structured modular kinetic model of Maria [35,48] used in this numerical analysis includes modules characterizing the dynamics of the concerned cell pathways involved in TRP synthesis, i.e., glycolysis, ATP recovery system, TRP operon expression, and biomass growth. This bioprocess model was experimentally identified and checked over extensive experiments conducted by several authors, i.e., [33,48,52,68] for glycolysis, and by Chen et al. [73,74] and Maria [35] for TRP synthesis. Experimental data of Chen [74] for TRP synthesis were also used to compare the derived predictions of the hybrid model.
The present study presents multiple elements of novelty: (i) although production of TRP by engineered E. coli has been extensively studied, "the need of multiple precursors for its synthesis and the complex regulations of the biosynthetic pathways make the achievement of a high product yield still very challenging" [35]. This engineering problem was solved here by using a model-based (in silico) approach, completed with a biological improvement of the used E. coli cell culture; (ii) the derived optimal operating policy of the FBR is given in time intervals (the so-called "time-arcs") of equal length, with a reduced number, to be easily implemented. The control variables present optimal but constant levels over each time-arc (different between time-arcs) during the FBR operation; (iii) the used biomass culture refers to a modified E. coli T5 strain. The characteristics of this strain were reflected in the rate constants estimated by Maria [35]. This T5 strain was produced by Chen [75] and Chen et al. [73] to increase TRP production in their bench-scale FBR. They performed genetic modifications of the TRP-producing "wild" strain S028. Basically, "they removed the PTS import-system of GLC, by replacing it with a more effective one based on the galactose permease/glucokinase (GalP/Glk) uptake system, by modulating the gene expression of GalP/Glk. The resulting T5 strain showed an increase of the specific TRP production rate in a nonoptimal FBR by 52.93% (25.3 mg/gDW biomass /h) compared to the initial strain" [73] and by ca. 70% if the used FBR was optimally operated (this paper); (iv) the results reveal the close link between the cell key metabolites and the FBR operating conditions; (v) the used hybrid bilevel kinetic model is complex enough to adequately represent the dynamics of the FBR state variables (i.e., the biomass growth, the GLC depletion, and the excreted TRP and PYR in the bulk phase), as well as the dynamics of the cell key species involved in the concerned reaction pathway modules, i.e., (a) glycolysis, (b) ATP recovery system, and (c) TRP operon expression. Comparison between the reduced schemes for GLC import systems into the cell linked to the TRP synthesis. Adapted from [73,74] (see the acknowledgement). (A) The wild E. coli model of Chassagnole et al. [52] and Maria [48] uses the phosphoenolpyruvate/sugar phosphotransferase (PTS) system for the GLC uptake. (B) The modified E. coli T5 strain of Chen et al. [73] and Chen [74], studied in this paper, uses the more efficient GLC uptake system based on galactose permease/glucokinase (GalP/Glk). The numbers on arrows indicated the relative metabolic fluxes at QSS predicted by Chen [74,75]. The same authors predicted a maximum theoretical yield of 0.23 g Trp/g glucose for the wild E. coli strain and of 0.45 g Trp/g glucose for the modified T5 strain.

Experimental Bioreactor and the Recorded Kinetic Data
To estimate the rate constants of the hybrid structured kinetic model for the studied TRP synthesis using the modified E. coli T5 strain (Figure 3), Maria [35] used the experimental kinetic data of Chen [74] obtained in a lab-scale three-phase FBR operated under the so-called "nominal" (nonoptimal) conditions displayed in (Table 1). The completely automated FBR of 1.5 L capacity includes a large number of facilities described in detail by Chen [74]. The nominal operation of this bioreactor by Chen [74] requires addition of a controlled constant feed flow rate of substrate solution (GLC) of a constant concentration, together with nutrients, antibodies, etc. in recommended amounts (Table 1) along the entire batch. A reduced FBR scheme can be found in the upper left corner of Figure 4.
To obtain the necessary kinetic data, samples were taken from the FBR bulk during the batch (63 h), with a certain frequency (2 to 5 h), thus determining the concentration dynamics of the key species of interest, i.e., X (biomass), GLC, TRP, and PYR. These recorded data are presented in Table 1 (see also the blue points in Figures 4-8). The reader Figure 2. Comparison between the reduced schemes for GLC import systems into the cell linked to the TRP synthesis. Adapted from [73,74] (see the acknowledgement). (A) The wild E. coli model of Chassagnole et al. [52] and Maria [48] uses the phosphoenolpyruvate/sugar phosphotransferase (PTS) system for the GLC uptake. (B) The modified E. coli T5 strain of Chen et al. [73] and Chen [74], studied in this paper, uses the more efficient GLC uptake system based on galactose permease/glucokinase (GalP/Glk). The numbers on arrows indicated the relative metabolic fluxes at QSS predicted by Chen [74,75]. The same authors predicted a maximum theoretical yield of 0.23 g Trp/g glucose for the wild E. coli strain and of 0.45 g Trp/g glucose for the modified T5 strain.

Experimental Bioreactor and the Recorded Kinetic Data
To estimate the rate constants of the hybrid structured kinetic model for the studied TRP synthesis using the modified E. coli T5 strain (Figure 3), Maria [35] used the experimental kinetic data of Chen [74] obtained in a lab-scale three-phase FBR operated under the so-called "nominal" (nonoptimal) conditions displayed in (Table 1). The completely automated FBR of 1.5 L capacity includes a large number of facilities described in detail by Chen [74]. The nominal operation of this bioreactor by Chen [74] requires addition of a controlled constant feed flow rate of substrate solution (GLC) of a constant concentration, together with nutrients, antibodies, etc. in recommended amounts (Table 1) along the entire batch. A reduced FBR scheme can be found in the upper left corner of Figure 4.
To obtain the necessary kinetic data, samples were taken from the FBR bulk during the batch (63 h), with a certain frequency (2 to 5 h), thus determining the concentration dynamics of the key species of interest, i.e., X (biomass), GLC, TRP, and PYR. These recorded data are presented in Table 1 (see also the blue points in Figures 4-8). The reader interested in the analytical methods used to obtain the experimental data, as well as in the details related to the bench-scale bioreactor operation and to the data acquisition system, is referred to the PhD thesis of Chen [74] (see also the Acknowledgement). the details related to the bench-scale bioreactor operation and to the data acquisition system, is referred to the PhD thesis of Chen [74] (see also the Acknowledgement). Simplified structured reaction pathway in E. coli for glycolysis (after Maria [48]), and for the TRP synthesis (the gray area) (after Maria et al. [35,37,48]). This reaction pathway was used by Maria et al. [34,37] to derive a TRP synthesis kinetic model. Connection of the TRP synthesis to glycolysis is realized through the PEP node [33,37]. The modular model structure also includes the synthesis of adenosine cometabolites ATP, ADP, and AMP, as part of the ATP recovery system (the pink rectangle in the figure). Notations: GLC(ex)= glucose in the cell environment. Species abbreviations are given in the abbreviations list. Species in parenthesis are not explicitly included in the glycolysis model. Italic letters denote the enzymes. Squares include notations of enzymatic reactions V1-V6 included in the glycolysis model (Table 2; Table 3). Adapted from [48] with the courtesy of CABEQ Jl., and completed according to the Maria [35] kinetic model. Simplified structured reaction pathway in E. coli for glycolysis (after Maria [48]), and for the TRP synthesis (the gray area) (after Maria et al. [35,37,48]). This reaction pathway was used by Maria et al. [34,37] to derive a TRP synthesis kinetic model. Connection of the TRP synthesis to glycolysis is realized through the PEP node [33,37]. The modular model structure also includes the synthesis of adenosine cometabolites ATP, ADP, and AMP, as part of the ATP recovery system (the pink rectangle in the figure). Notations: GLC(ex)= glucose in the cell environment. Species abbreviations are given in the abbreviations list. Species in parenthesis are not explicitly included in the glycolysis model. Italic letters denote the enzymes. Squares include notations of enzymatic reactions V1-V6 included in the glycolysis model (Tables 2 and 3). Adapted from [48] with the courtesy of CABEQ Jl., and completed according to the Maria [35] kinetic model.           Table 1, with a constant feed flow rate, and with a constant GLC concentration in the feed. Both cases use the same modified E. coli T5 strain. (Bottom curves). Model-based simulated trajectories (-) of glucose (GLC) in the bioreactor bulk for the FBR operated in two alternatives: (i) (2, black) optimal operation derived in this paper (variable fed [GLC] and variable feed flow rate), and (ii) (1, blue) experimental data (•, blue) of Chen [71] derived under nominal, nonoptimal operation of Table 1, with a constant fed [GLC] and a constant feed flow rate.  [74] to collect the kinetic data of the TRP synthesis using a suspended culture of genetically modified E. coli cells (T5 strain). More experimental details were given by Chen [74].

Species Mass Balance Auxiliary Relationships and Estimated Rate Constants
glc,j = control variables to be optimized; j = 1, . . . , N div (equal time-arcs) c ext glc (t = 0) is given in (Table 1) for the nominal FBR of Chen [74] For the optimal FBR with adopted N div = 5, the feeding policy is (Footnote a): [33,86,87] (ii) c adp results from solving the thermodynamic equilibrium relationship c atp c amp = Kc 2 adp , i.e., The initial values of cell species concentrations are given in Table 1 (see also footnote (b)) (v) The lump c tca of Figure 3 includes species belonging to the TCA cycle; there are no measurements on this lump, so it was excluded from data fitting (vi) The adopted value for y trp by Maria [35] is y trp = r syn,trp /r syn,pep = 1/43.63 (at QSS) [88]; y trp was re-estimated from experimental data by Maria [35], resulting in y trp = 0.467 (vii) See Table 3 for the V 1 − V 6 flux expressions Species inside the cell Table 1; j = 1, . . . , N div (equal time-arcs) (viii) For the adopted N div = 5, the feeding policy is (see footnote (a)) (Table 1) (ix) The biomass growth inhibition corresponds to a modified Contois model [85] The estimated rate constants by Maria [35] are Table 3. Reaction rate expressions V1-V6 of the hybrid model of Table 2, describing the dynamics of the cellular glycolytic species according to the kinetic model of Maria [35,48] and of Chassagnole et al. [52]. In the present study, this glycolysis kinetic model was modified by replacing the PTS system (V1 flux) for the GLC uptake with those of the mutant T5 E. Coli strain tested in this paper. The model rate constants were estimated by Maria [46] to fit the experimental data of Chen [71] (presented in Table 1 and Figures 4-8). Species abbreviations are given in the abbreviation list.

Rate Expressions Estimated Rate Constants (Units in mM, min)
GLC import system glc + pep → f6p + pyr pyr + atp → pep + adp + h glc + atp → f6p + adp +h Modification for the T5 strain (i) Termonia and Ross [86,87] indicated experimental evidence of a very fast reversible reaction catalyzed by AKase, with the equilibrium being quickly reached (ii) The k6 constant takes values according to the microorganism phenotype (related to the gene encoding the enzyme ATPase that catalyzes this reaction) (iii) c amp + c adp + c atp = c amdtp = constant [86,87] (iv) c adp results from solving the following thermodynamic equilibrium relationship:    [74] to collect the kinetic data of the TRP synthesis using a suspended culture of genetically modified E. coli cells (T5 strain). More experimental details were given by Chen [74].

The FBR Initial Conditions Parameter
Nominal (Initial) Value Obs.

The Structured Hybrid Kinetic Model of Maria
Being a metabolite of high practical importance, intense efforts have been invested to decipher the synthesis regulation mechanism of TRP in various microorganisms, for deriving an adequate dynamic model of its QSS or oscillatory synthesis to be used for engineering purposes. Some results include the deterministic kinetic models of [37,81], while other studies [63] rather focused on determining correlations among flux distribution, flux control, and the optimized enzyme amount distribution, by employing a reduced kinetic model, not able to simulate most CCM reaction pathways.
As the TRP synthesis regulation is a very complex process, a significant number of simplified kinetic models with lumped terms (species and/or reactions) have been proposed in the literature (see the review of Maria et al. [34,37]). Kinetic modeling of this complex process is even more difficult because, as proven by Xiu et al. [82,83], Chen et al. [84], and Maria [33,36,37], under certain FBR operating conditions, TRP synthesis can become an oscillatory process. Oscillations in the TRP synthesis are produced due to the concomitant activation and high-order repression of the TRP operon expression, together with a nonlinear demand for end product, making its expression cyclic. The cell growth and dilution rates (related to the cell cycle and the liquid residence time in a (semi-)continuous bioreactor) strongly influence the TRP system stability, as proven in silico by Maria [33] and Maria et al. [34].
The adopted hybrid kinetic model is that of Maria [35] built up using the kinetic data of Chen [74] collected in an FBR operated under the nominal (nonoptimal) conditions of ( Table 1), using the T5 strain of E. coli. This complex structured kinetic model (presented in Tables 1, 2 and 4) is a deterministic one. The CCM-based model core is the glycolysis dynamic model of Maria [48] validated using literature data.  [81] were modified by Maria et al. [37] to better fit the experimental data, as follows: (i) PEP (from glycolysis) is the substrate of TRP synthesis and the node coupling this synthesis with the glycolysis [a] module; (ii) a novel model for the TRP synthesis inhibition was proposed and identified from experiments. The model rate constants were estimated by Maria [35] to fit the experimental data of Chen [74] (Figures 4-8) collected in the FBR using the modified E. coli T5 strain, under "nominal" operating conditions (Table 1). Species notations (TRP, OR, OT, mRNA, and E) are given in the abbreviation list. QSS = quasi-steady state.

Rate Expression
Kinetic Model Parameters (Units in mM, µM, min) (a) The adopted modification for the TRP synthesis inhibition replaces the C3 variable of the Bhartiya et al. [81] model (not displayed here, see [35]) with a modified Contois model, including a power-law inhibition with TRP growth at the denominator. (b) The nitrogen source in the TRP synthesis is considered in excess and included in the model constants. (c) To be connected to the glycolysis kinetic model, the PEP species dynamics, generated by the glycolysis model, was explicitly included in the TRP synthesis rate as a substrate [35]. (d) The initial concentrations of the TRP operon species (OR, mRNA, E, and TRP) are given in Table 1.
To keep the bilevel hybrid model of Maria [33,35] adapted here with a reasonable extension, as well as to facilitate estimation of its rate constants, this dynamic model accounts for only the key species included in four linked cell reaction modules [a-c, X] responsible for the TRP synthesis. Three modules concern the following cell processes [35]: Module [a]-glycolysis with a modified GLC uptake system (due to the used modified E. coli T5 strain); The fourth kinetic module concerns the biomass [X] growth dynamics in the FBR bulk. This last module is connected to the cell processes, by influencing the GLC dynamics in the bulk phase through the X growth rate (Table 2), which, in turn, influences the GLC import flux V1 into the cell ( Table 3). The dynamic model is hybrid (bilevel) because it connects the macro state variables of the FBR (biomass X, GLC, TRP) with the cell nano-level key variables (GLC, F6P, FDP, PEP, PYR, and ATP; Tables 2 and 3) of the glycolysis and those (TRP, OR, OT, and mRNA) of the TRP operon expression (Table 4). All four kinetic modules are linked to the macroscopic FBR dynamic model through the formulated mass balances (Tables 1, 2 and 4).
The thus obtained hybrid extended kinetic model includes a large number (49) of rate constants. To facilitate the estimation rule and to avoid suboptimal estimates (i.e., rate constants, locally valid in the operating parametric space), only (27) independent rate constants were accounted for in the estimation step [35]. Moreover, a step-by-step estimation methodology was applied by Maria [35], by decomposing the estimation problem in successive subproblems of smaller dimensions, by taking advantage of the modular construction of the cell FBR hybrid model. The estimated rate constants by Maria [35] were validated by fitting the experimental key species kinetic curves of Chen [74] recorded in the FBR of ( Table 1) over a long batch (63 h), using the novel E. coli T5 strain of Chen et al. [73], as well as by comparison with the literature data reviewed by Maria [34,36,37,48,49]. In short, the methodology used by Maria [35] to estimate the adopted bilevel modular dynamic model consisted of a sequence of trial-and-error steps, by adjusting the literature information (reaction rate expressions and constants characterizing the dynamics of cell metabolic species of interest) to fit the available experimental kinetic data recorded from the above-described FBR. The sequence of computational steps is summarized below.

The Biomass [X] Growth
The cell culture in the bioreactor is considered to be homogeneous and introduced as a lump "X" in the FBR model (Table 2). A modified Contois model, modified by considering a power-law inhibition with the first-order growing biomass at the denominator [85], was proven to be the most adequate vs. the experimental data (Figure 8c). To overcome the absence, at this modeling stage, of the predicted values of [X] and [GLC] (coming from the FBR coupled with the glycolysis dynamic models), simulations of the biomass dynamics over the batch were performed using the experimentally recorded [X] and [GLC] species trajectories of Chen [74], interpolated with the cubic splines functions (INTERP1 facility of Matlab TM package). The estimated kinetic model of the biomass is given in Table 2.

The FBR Dynamic Model
The FBR ideal model of Maria [35] was adopted to describe the key species dynamics during the batch at a macroscopic level (in the liquid bulk phase). The bioreactor initial conditions and the time stepwise dynamics of the control variables (added GLC substrate solution concentration, and the feed flow rate F L ) were further explored to obtain the desired optimum operation of the studied FBR. The bioreactor ideal model main assumptions were as follows [2]: (i) isothermal, iso-pH, iso-DO operation; (ii) it is self-understood that nutrients, additives, antibiotics, and pH-control compounds are added initially and during FBR operation to ensure the optimal grow of the biomass, as indicated by Chen [74]; (iii) oxygenation with pure oxygen in excess over the batch to ensure an optimal biomass maintenance, and to contribute to the liquid homogeneity; (iv) perfectly mixed liquid phase (with no concentration gradients, see Table 1), of a volume increasing according to the liquid feed flow rate time-varying policy; (v) the limits of the liquid feed flow rate (F L,j in Table 2) are adjusted to not to exceed the bioreactor capacity (Max(V L ) in Table 1); (vi) negligible mass resistance to the transport of oxygen and compounds into the liquid and biomass flocks (if any); (vii) GLC substrate is initially added in the bioreactor and during the batch according to an optimal feeding policy to be determined; (viii) the feed flow rate during the batch F L,j is varied according to an optimal feeding policy to be determined for every "time-arc" index J in Equation (5).
From a mathematical point of view, in a general form, the FBR dynamic hybrid model (Tables 1, 2 and 4) translates to a set of 12 differential mass balances (ODE set) written as below for the key species of the FBR.: Species in the bulk-phase: where index "i" denotes species present in the FBR bulk (GLC, TRP, PYR, X); index "j" denotes the FBR feeding time-arcs; j = 1, . . . , N div . Key species inside cells: where index "i" denotes species inside cells, i.e., GLC, F6P, FDP, PEP, PYR, and ATP for glycolysis and OR, mRNA, E, and TRP for TRP operon expression. Biomass in the bulk phase: Liquid volume dynamics: In Equation (1), c inlet,i,j refers to the concentration of the species index "i" in the feeding solution, constant over the time interval index "j" (j = 1, . . . ,N div ). In the present case, only GLC is fed into the FBR during the batch. The reaction rate r i expressions together with the associated rate constants and other details are given in Tables 1 and 4. In Equations (1)-(3), c is the vector of species concentrations, c o is the initial value of c (at time t = 0) given in (Table 1), and is the vector of the model rate constants. The reactor content dilution (determined by the increasing V L in Equation (4) is due to the continuously added F L term.
Feed flow-rate policy: Similarly, for the adopted N div = 5 equal time-arcs, the feeding policy for the GLC solution concentration is To not complicate the engineering calculus, the main assumption in Equations (5) and (6) is the following: on each time stepwise "arc", index j = 1, . . . , N div , the control variables  (Table 1). Moreover, they are also the same, i.e., F L,0 = F L,1 = F L,2 = F L,3 = F L, 4 , and FBR optimal operating conditions. By contrast, under the optimal conditions studied in this paper, the suitable time stepwise values F L,0 − F L, 4 and those of c f eed glc,0 − c f eed glc,4 are to be determined together (simultaneously) to reach the optimum of an objective function (maximum of TRP production here). Multi-objective FBR optimization is also possible (see [89,90]) but is beyond the scope of this paper.

Module [a] Glycolysis and Module [b] ATP Recovery System
Once a dynamic model able to predict the biomass dynamics in the FBR was adequately fitted (Section 3.1.1), two other kinetic modules were considered for rate constant estimation, i.e., glycolysis *module [a]) and ATP recovery system (module [b]). Their reduced reaction pathways are given in Figure 3, while the reaction rate expressions are given in Tables 2 and 3. The model of Maria [35] was adapted on the basis of the literature information [33,34,37,48,52], by performing only a few modifications in the reaction rate expressions to reflect the modified GLC import system of the modified E. coli T5 strain schematically represented in Figure 2B , and [X]) kinetic models to fit the experimental kinetic data (i.e., the species dynamic trajectories recorded by Chen [74] in the FBR of Table 3) under the "nominal" operating conditions. The results are presented in Tables 1 and 2. For supplementary details, the reader is referred to Maria [35].
In short, glycolysis module [a] is a determined sequence of 10 enzyme-catalyzed reactions (see the reduced pathways of Figures 2 and 3 with only six lumped reactions) that convert 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 several enzymatic syntheses in the cell [91]. Adequate modeling of the glycolysis dynamics is important because the glycolytic intermediates provide entry/exit points to/from glycolysis. Thus, most of the monosaccharides, such as fructose or galactose, can be converted to one of these intermediates, further used in subsequent pathways. For example, PEP is the starting point for the synthesis of essential amino acids (AAs) such as tryptophan, cysteine, arginine, and serine [37,52,88,92].
Due to the tremendous importance of glycolysis in simulating the cell CCM, intense efforts have been made both in its experimental study and in modeling the dynamics of this process specifically in bacteria (short reviews [33,48,93]). The large number of glycolysis reduced or extended kinetic models proposed in the literature (review [48]) present a complexity ranging from 18-30 species, including 48-52 reactions, with a total of 24-300 or more rate constants. Most of these models are, however, too complex to be easily identified from (often) few available kinetic data and too complex to be further used for engineering calculations. Moreover, with a few exceptions, most of them cannot satisfactorily reproduce the occurrence of glycolytic oscillations on a mechanistic basis [33,36].
The adopted glycolysis kinetic model of Maria [33,48] even if of a reduced form, by accounting only for nine key species in lumped reactions including 17 easily identifiable rate constants belonging to V1-V6 metabolic fluxes ( Figure 3, and Tables 2 and 3) has been proven to adequately reproduce the cell glycolysis under steady-state, oscillatory, or transient conditions according to (i) the defined glucose concentration level/dynamics in the bioreactor bulk (liquid) phase, (ii) the total A(MDT)P cell energy resources, and (iii) the cell phenotype characteristics related to the activity of enzymes involved in the ATP utilization and recovery system (here denoted as module [b]). Detailed discussions about the operating conditions leading to glycolytic oscillations were extensively presented by Maria [33,36,37]. For this reason, the FBR and the glycolysis dynamic models have to be considered together (Tables 2 and 3) when simulating the dynamics of the [GLC] in the FBR bulk phase and of the metabolites of interest (F6P, FDP, PEP, PYR, and ATP) into the cell. The adopted rate expressions for the glycolysis main fluxes V1-V6 presented in Tables 2 and 3 are those of the basic model, except those of the GLC import system (V1), modified to match the T5 E. coli strain kinetic data [35]. It is worth mentioning that, even if not the case here, under certain conditions (i.e., external/environmental and internal/genomic factors), glycolysis and TRP synthesis can become oscillatory processes [34,36,37,48,49]. According to the experimental data, the produced TRP (module [c]) is excreted (Figure 3) through a process described by Chen [74]. The PYR key metabolite concentration in the cell is regulated through a complex mechanism [94,95], with the excess being excreted, as experimentally proven by Chen [74].
As revealed by the reactions in the pink square of Figure 3, the efficiency and the dynamics of the ATP recovery system are essential for the reaction rates of the whole CCM, as long as ATP plays a catalytic-chemical energy provider role. As underlined by Maria et al. [33,36,37], among the involved parameters, an essential factor is the k6 reaction rate (determined by the ATPase characteristics in Figure 3), included in the glycolysis model of Tables 2 and 3. The involved enzymes characteristics are directly related to the cell phenotype (i.e., cell genomic) controlling the [AMDTP] total energy resources level. To not complicate the simulations, the [AMDTP] level was kept unchanged in the present analysis at an average value given in Table 1, as suggested by Chassagnole et al. [52]. The adopted kinetic model for the glycolysis (i.e., the V1-V6 reaction rates of Figure 3 and Tables 2 and 3) and the equilibrium relationships for the ATP-ADP-AMP system given in Tables 2 and 3 were imported from the literature [33,35,48]. This kinetic model was proven by Maria [35], according to experimental checks to fairly represent the dynamics and the thermodynamics of the internal modules [a,b] in the modified E. coli T5 strain.

Module [c] TRP Synthesis
The adopted in silico evaluation of the TRP synthesis of Maria [35] is based on a simplified pathway displayed in Figure 3, derived from various studies reviewed by Maria et al. [37]. Modeling the TRP synthesis using a deterministic (mechanism-based) approach is difficult because this cellular process is known as being, under certain conditions, a QSS or an oscillatory one [33,81,83]. However, to avoid extended models, difficult to be estimated and used, most of the reduced dynamic models from the literature do not distinguish the process components from the regulatory components, and lumped reactions/species are considered instead, with the regulatory performance being included via adjustable model parameters and terms. Kinetic models trying to reproduce the TRP operon expression self-regulation [82,83] are too extended for our engineering evaluation purposes. Due to the process complexity, some modeling approaches [63] instead focused on determining correlations among flux distribution, flux control, and the optimized enzyme activity distribution, by employing a reduced kinetic model, not able to simulate most CCM key modules.
For such reasons, in the present analysis, simulations of the TRP synthesis were performed using the reduced CCM-based kinetic model of Maria et al. [33,37]. The adopted dynamic model of Maria [35] for the TRP synthesis (TRP operon expression) is given in Table 4. This kinetic model is a modification of those proposed by Bhartiya et al. [81]. The operon expression regulation terms (C1,C2) were kept unchanged. Only the TRP mass balance was changed according to the below reasons. The rate constants of the considered OR, mRNA, TRP, and E key species mass balances were re-estimated using the experimental data of Chen [74] given in Figures 4-8. The TRP mass balance of the Bhartiya et al. [81] model was modified and re-estimated step by step as follows: i.
An explicit connection of the TRP module to the glycolysis module [a] pathway was introduced through the PEP precursor sharing node (in Figure 3). Consequently, PEP is included as a substrate in the TRP mass balance (dc TRP /dt in Table 4), while the PEP consumption term is also considered in the PEP balance of the glycolysis model according to the recommended fluxes ratios of Stephanopoulos and Simpson [88], as a first guess ( Table 2). Analysis of this model suggests that intensifying TRP synthesis clearly depends on the glycolysis intensity (average levels of glycolytic species) and its dynamics (QSS or oscillatory). In fact, as remarked by Li et al. [78] and by Chen and Zeng [76], the PEP precursor is the limiting factor for TRP synthesis. This is why intense efforts have been made to increase its production through glycolysis intensification [33,34]. This can be realized by optimizing the FBR operating policy (as in this paper) and/or by using (also in this paper) the modified E. coli T5 strain culture of Chen et al. [73] and Chen [74]. ii. The TRP synthesis model of Bhartiya et al. [81] (Table 4) includes two terms for the product inhibition, i.e., the C3 term (of allosteric-type) plus a Michaelis-Menten term. Our tests proved that these terms do not adequately fit the TRP experimental kinetic data of Figure 4. Accordingly, the product inhibition term in the TRP balance of  Table 4, an activation inhibition term was considered by bringing together the substrate (PEP) and the first key enzyme (anthranilate synthase, E) that trigger TRP synthesis [35]. Such an approach was proven to better fit the experimental data of Figure 4, i.e., c trp (t u ), u = 1, . . . , n (where n = 17 denotes the number of experimental points) and to confer more flexibility to the derived model. The estimated negative g constant, of a small negative value, reflects the slight inhibition of TRP synthesis with the substrate PEP, as also suggested in the literature [35].

Rate Constant Estimation by Maria (2021)
In total, the developed hybrid structured kinetic model includes 49 rate constants to be estimated from the experimental kinetic curves of four observed species (GLC, TRP, PYR, and X), with each species time trajectory including 17 uniformly distributed recorded points (Figures 4-8). This estimation problem is equivalent to a nonlinear programming one (NLP) of high difficulty [41] due to its dimension, the high nonlinearity of the model, and its associated constraints.
To avoid unfeasible local estimates of the NLP problem, Maria [ [X]. To reduce the problem size, only 27 independent model rate constants were accounted during estimation (from the total of 49 rate constants). A number of rate constants were adopted from the literature [34,37]. Eventually, all rate constants were refined by Maria [35], as presented in Tables 1, 2 and 4. The thus identified FBR hybrid structured dynamic model fit the experimental data very well, as indicated by Figures 4-8.
The multimodal NLP estimation problem solved by Maria [35] is a difficult one, being highly nonlinear, including nonlinear constraints defining a nonconvex domain. For such large nonconvex estimation problems, the usual optimization routines usually encounter difficulties in reaching the feasible global solution with an acceptable reliability. This is why a very effective NLP solver was used instead, i.e., the adaptive random search MMA of Maria [96] implemented on the Matlab TM numerical calculus platform. The NLP solution was checked using several (randomly generated) initial guesses for the rate constants. A stiff integrator (ODE15S routine of Matlab TM package) was used to solve the ODE dynamic model with a high accuracy.
A comparison of the model-estimated rate constants for the modified T5 E. coli strain using the FBR experimental data of Chen [74] with those of the same model but estimated from experiments using the wild E. coli strain was presented by Maria [35]. As expected, most of the estimated rate constants presented similar values for some reaction steps. However, due to the mentioned modifications of the used E. coli T5 strain in the present kinetic model, important differences were reported for (i) the rate expression and parameters of the GLC import system (V1 in Tables 2 and 3, and Section 3.1.3), (ii) the biomass growing dynamics (Table 2), and (iii) the TRP synthesis module [c], in terms of both parameters and rate expressions ( Table 3). As another observation, for the nominal (nonoptimal) FBR experimental conditions of Table 1 used by Chen [74], the species dynamics belonging to inside the cell and to the external liquid phase tend to reach a quasi-steady state (QSS) that corresponds to a balanced cell growth (homeostasis) in the bioreactor [35].

Ways to Intensify the TRP Production in the FBR
As revealed by the concerned literature [34][35][36][37]74], intensifying the TRP synthesis strongly depends on (a) the glycolysis intensity (GLC uptake flux, and average levels of glycolytic species of module [a]), transmitted to the TRP synthesis module [c] via the shared PEP intermediate, and (b) on the glycolysis dynamics (QSS, or oscillatory behavior) [33,34,37]. More specifically, as pointed out in the literature by Maria et al. [36] and Maria [33], the glycolysis intensity is controlled by both cell internal and external factors, as follows: i. The GLC import system efficiency (V1 in Figure 3) is regulated and triggered by the external concentration of glucose and by the subsequent PEP and PYR synthesis (Tables 2 and 3). The regular GLC uptake system, i.e., the PTS translocation system in the wild strain (of a complex reaction rate expression discussed by [35,48,52]) was replaced in the present studied E. coli T5 strain, as mentioned in Section 2, with a more efficient one ( Figure 2B) able to accelerate the GLC uptake flux into the cell at least twofold [74]. Such a modified GLC import was modeled by simple Michaelis-Menten kinetics in the model of Table 3 by accounting for the well-known GLC substrate inhibition. ii. The quick import of GLC and its conversion to the precursor PEP require important amounts of regenerable ATP and a rapid enough ATP-to-ADP conversion rate, as well as its quick regeneration. The re-estimated rate constants of the kinetic module [b] (pink rectangle in Figure 3, and Section 3.1.3), concomitantly with those of the kinetic module [a] from the experimental data coming from the FBR operated with modified E. coli cells implicitly ensure the requirement that the A(MDT)P energy system is able to support the cell glycolysis (see V2, V4, and V6 expressions in Table 3 and the ATP mass balance in Table 2). On the other hand, limited A(MDT)P energy resources which exist in the cell slow down the GLC import if the ATP use/regeneration is not working fast enough [97]. Such an A(MDT)P resource is linked to the microorganism phenotype. Here, the total A(MDT)P was adopted (Tables 1 and 3) at the average level recommended by Chassagnole et al. [52]. iii. Additionally, due to the enzyme ATPase and AKase characteristics related to the bacteria genome and cell phenotype (Figure 3), a limited ATP conversion rate can sustain the glycolytic reactions, while the ATP recovery rate is limited by the enzymes partici-pating in the A(MDT)P interconversion reactions (i.e., the K and k6 rate constants in the kinetic model of Table 3). This is why the k6 rate constant was re-estimated here to fit the experimental data, as suggested by Maria et al. [36,49]. iv. At the same time, as glycolysis is a systemic process, with a complex regulatory structure, its dynamics (oscillatory, transient, or QSS) is also related to the rate constants of all involved reactions. Consequently, all these rate constants were considered in the final estimation step of the whole FBR hybrid kinetic model. Similarly, Silva and Yunes [98] found that glycolysis (QSS or oscillatory) is only possible if the external concentration of GLC and the maximum reaction rates controlled by the enzymes PFKase and GKase (which control the V1 and V2 reactions of Figure 3) are within specific intervals. Due to the same reason, the rate constants related to the GLC uptake system in the modified E. coli cell (modified V1 flux in Table 3) were re-estimated to match the experimental kinetic data. v.
As a corollary of the issue (iv), Maria [33,34,36,37] determined the operating conditions leading to glycolytic oscillations or QSS by varying the external factor [GLC]ext and some internal factors such as the total [AMDTP] level and the k6 rate constant of Table 3. Such an investigation was not necessary here, because no oscillatory process was identified in the present operating case. vi. Simulations by Maria [33,35] revealed that the TRP synthesis efficiency is also strongly influenced by external factors, related to the FBR operating regime, namely, (a) the cell dilution (taken into account as "µ" in the approached hybrid kinetic model of Table 2, (b) the GLC concentration in the external (bulk) phase (c ext glc in Table 2), and (c) the optimal operating policy for the control variables. In this paper, such an operating policy will correspond to the time stepwise variation of the feed flow-rate (F L,j in Equation (5)) and of the GLC feeding concentration (c f eed glc,j in Equation (6)).

Preliminary Considerations
To support further engineering calculations, a reasonable extended hybrid modular approach was adapted from literature [35], by expressing the macroscopic main state variable species dynamics (i.e., biomass X, GLC, and TRP) governing the FBR performance, as a function of intracellular species dynamics related to the cell CCM metabolic fluxes responsible for the TRP synthesis. This inner cell environment link is realized by means of model key species (GLC, X, PEP, ATP) (Section 3.1). The main modification in this paper of this adopted hybrid dynamic model refers to the introduction of a variable FBR feeding both in the feed flow rate F L (Equation (5) and Table 2) and in the GLC feeding solution concentration (Equation (6) and Table 2).
The reasonable compromise between the hybrid model details (number of intracellular species and reaction pathways accounted for) and its predictive value was realized by using only the cell key-modules  1.2 and 3.1.3). The fair adequacy of the resulted dynamic model (Tables 1, 2 and 4) vs. the experimental data was proven by Maria [35]. Consequently, this hybrid model becomes suitable for further engineering evaluations of the reactor and process efficiency, as is the case here.
The optimal FBR operation derived in this paper is more complex than the simple nonoptimal ("nominal") operation of Chen [74] (Table 1). Mainly, the feed flow rate and GLC concentration in the feeding solution are no longer kept constant. In contrast, (i) the batch time is divided in N div (equal "time-arcs") of equal lengths, and (ii) the control variables are kept constant only over every "time-arc" at optimal values for each timearc determined from solving an optimization problem (i.e., maximization of the TRP production in this case). The time intervals of equal lengths ∆t = t f /N div are obtained by dividing the batch time t f into N div parts t j−1 ≤ t ≤ t j , where t j = j∆t are switching points (where the reactor input is continuous and differentiable). Time intervals for the present case study with an adopted N div = 5 are shown in the "liquid volume dynamics" row of Table 2 and its footnote (a).

Selection of the FBR Control Variables
By analyzing the FBR hybrid model of Table 2, completed by Table 3; Table 4, the natural option is to choose as control variables those with a high influence on the biological process, which are easily to handle. In the present case, according to the discussion of Section 3.3, two control variables were chosen related to the bioreactor feeding: (a) The substrate c f eed glc,j (j = 1, . . . , N div ) whose concentration plays a major role in the cell glycolysis and TRP production; (b) The liquid feed flow rate F L,j (j = 1, . . . , N div ), with a GLC solution directly linked to the GLC feeding, responsible for the reactor content dilution.
In the present optimization strategy, each control variable is kept constant over each time-arc (index "j"). Of course, they are not necessarily equal between different time-arcs. For N div = 5, in total there are 5 × 2 = 10 unknowns in Equation (7) to be determined by optimization, under certain constraints ( Table 2): The FBR initial state is given in Table 1 for both inside cell and bulk-phase species. Those of the control and bulk phase variables, i.e., the initial liquid flow rate and the substrate initial concentration (as shown in Table 2, and Equations (5) and (6)) are included as unknown variables in the FBR optimization, i.e., in Equation (5), and in Equation (6).

Objective Function (Ω) Choice
By considering the mentioned control variables Equation (7), the FBR optimization consists of determining its optimal initial load simultaneously with its feeding policy for every time interval during the batch, eventually leading to maximization of the [TRP] production during the batch.
The control variables values of Equations (7)-(9) to reach Max Ω were identified, where The [TRP](t) dynamics in Equation (10) was evaluated in silico by solving the ODE dynamic model of the FBR (Equations (1)-(6)) over the whole batch time (t) ∈ [0, t f ].

Optimization Problem Constraints
The optimization problem in Equation (10) was subjected to the following multiple constraints: (a) The FBR model in Equations (1)- (6) including the bioprocess kinetic model (Tables 1, 2 and 4); (b) The FBR initial condition from Table 1, except F L,0 and c f eed glc,0 which were determined from solving the optimization problem (the initial guess was taken from Table 1); (c) To limit the excessive consumption of substrate and to prevent the hydrodynamic stress due to the limited reactor volume, feasible searching ranges were imposed on the control/decision variables, i.e., As an observation, the imposed ranges for the control variables were related to not only the implementation facilities, but also economic reasons, achieving minimum substrate consumption, reduced dilution of the reactor content, and an effective bioreactor control.

N div and Operating Alternatives Choice
The adopted FBR operating policy alternative of Section 4.2.1 is one of the simplest variable operating modes. It implies a time stepwise variable feeding of the bioreactor, over an adopted (N div = 5 here) equal time-arc that covers the whole batch time. Each time-arc "J" (j = 1, . . . , N div ) is characterized by optimal levels of the feed flow rate F L,j and of the GLC concentration c f eed glc,j (see Equations (7)- (9)). This type of FBR operation, despite its simplicity and easy implementation, still includes enough degrees of freedom to offer a wide range of operating facilities that, in principle, might be investigated, for instance (see also the discussion of Maria [3]), (a) by choosing unequal time-arcs, of lengths to be determined by the optimization rule, (b) by considering the whole batch time as an optimization variable, (c), by increasing the number of equal time-arcs (N div ) to obtain a more refined and versatile FBR operating policy, but keeping the same nonuniform feeding policy (of the two control variables here), (d) by considering the search min/max limits of the control variables as unknown (to be determined), or (e) by feeding the bioreactor with a variable feed flow rate, but with a GLC solution of an uniform concentration over a small/large number (N div ) of time-arcs. All alternatives (a-e) are not approached here for the reasons discussed below.
Alternatives (a-c) are not good options, because, as N div increases, the necessary computational effort grows significantly (due to a considerable increase in the number of searching variables), thus hindering the quick (real-time) implementation of the derived FBR operating policy. Additionally, multiple optimal operating policies can exist for the resulting overparameterized constrained optimization problem of a high nonlinearity, thus increasing the difficulty in quickly locating a feasible globally optimal solution of the FBR optimization problem.
Additionally, as the N div increases, the operating policy is more difficult to implement, since the optimal feeding policy requires a larger number of stocks with feeding substrate solutions of different concentrations, separately prepared to be fed for every time-arc of the FBR operation (an overly expensive alternative). Moreover, the NLP optimization problem is more difficult to solve because of the multimodal objective function, leading to multiple solutions difficult to discriminate and evaluate. This is the case, for instance, of an obtained infeasible optimal policy requiring a very high [X], difficult to be ensured due to limitations in keeping the necessary levels of the related running parameters of the bioreactor (i.e., dissolved oxygen, nutrients, pH-control substances, antibodies, etc.). Furthermore, FBR operation using a larger number of small time-arcs N div can raise special operating problems when including PAT (process analytical technology) tools [99].
A brief survey of the FBR optimization literature [100][101][102] reveals that a relatively small number N div < 10 is commonly used for such an FBR due to the abovementioned reasons. In fact, the present numerical analysis does not intend to exhaust all the possibilities of the approached FBR optimization. Thus, an extended analysis of the operating alternatives (a-d) of the FBR operation or the influence of the parametric uncertainty deserves a separate investigation, beyond the scope of this paper. To not complicate the computational analysis, only N div = 5 equal time-arcs are tested here, with equal time-arc lengths of t f /N div = 63/5 h.
The alternative (d) is unlikely because it might indicate unrealistic results, as explained in point (c) of Section 4.2.3. In our numerical analysis, carefully documented upper bounds of control variables were tested to ensure the practical implementation of the optimal operating policy.
Alternative (e) is also not feasible, even if a larger N div is used. That is because it is well known that the variability of the FBR feeding over the batch time-arcs is the main degree of freedom used to obtain FBR optimal operating policies of superior quality [3,6,89,101]. By neglecting the variable feed flow rate and substrate concentration, suboptimal FBR operating policies will be obtained of low performance.

The Used Numerical Solvers
The prediction of the species concentration time evolution inside the cell and in the bulk phase was obtained by solving the FBR dynamic model in Equations (1)-(6) with the initial condition of C j,0 = C j (t = 0) of Table 1 for the inside cell species, except the bulk [GLC] 0 to be determined from the FBR optimization, as indicated by Equations (7) and (9). The imposed batch time t f and the optimal medium conditions are those of Table 1. The dynamic model solution was obtained with a high precision, using the high-order stiff integrator ("ode15s") of the MATLAB™ numerical calculus platform, with suitable quadrature parameters to keep the integration error very low.
Because the FBR hybrid model structure in Equations (1)-(6), its reaction rate terms (Tables 1, 2 and 4), and the problem constraints from Equations (11)-(14) (Section 4.2.3) are all highly nonlinear, the formulated problem in Equations (7)-(10) translates into a nonlinear optimization problem (NLP) with a multimodal objective function and a nonconvex searching domain. To obtain the global feasible solution with enough precision, the multimodal optimization solver MMA of Maria [41,96,103] was used, proven in previous studies to be more effective compared to the common (commercial) algorithms. The computational time was reasonably short (minutes) using a common PC, thus offering a quick implementation of the obtained FBR optimal operating policy [96,103].

The Problem Solution Particularities
The obtained optimal operating policy of the FBR, for the optimization problem formulated in the Section 4.2.2, with the control variables of Section 4.2.1, the constraints of Section 4.2.3, and adopted N div in Section 4.2.4, is given in Figure 7 for the feeding policy of the GLC concentration c f eed glc,j (j = 1, . . . , 5) and in Figure 8a for the feed flow rate F L,j (j = 1, . . . , 5). It is to be observed that, due to the above formulated engineering problem, the FBR optimal operating policy is given for every time interval (of equal length) uniformly distributed throughout the batch time.
Such an optimal time stepwise variable feeding of the bioreactor presents advantages and inherent disadvantages. The advantages are related to the higher flexibility of the FBR operation, leading to a higher productivity in TRP, as proven in Section 5. Furthermore, the imposed limits of the control variables prevent excessive substrate consumption or an excessive reactor content dilution.
As a disadvantage, FBRs with such time-variable control are more difficult to operate than simple BRs, as long as the time stepwise optimal feeding policy requires different stocks of feeding substrate solutions of different concentrations to be used over the batch. This is the price paid for achieving the best performance of an FBR. This need to previously prepare different substrate stocks to be fed for every "time-arc" (i.e., a batch-time division in which the feeding is constant) is offset by the net higher productivity of FBR compared to that of BR as discussed below and pointed out in the literature [6,89,90,101,104]. In fact, the best operating alternative (FBR vs. BR) is related to many other economic factors (operating policy implementation costs, product cost compared to production costs, product price fluctuation, etc.), not discussed here.

Optimization Results and Discussion
The obtained optimization problem solution (of the type discussed in Section 4.2.6) is given in Figure 7 (top, curve 2) for the GLC feeding concentrations and in Figure 8a (curve 2) for the feed flow rate. The optimally operated FBR displays the bulk [TRP] dynamics of Figure 4 (curve 2). The corresponding dynamics of cell glycolytic species during the batch is presented in Figure 5, while that belonging to the TRP operon expression is presented in Figure 6. The dynamics of species present in the reactor liquid phase are presented in Figure 7 for GLC and in Figure 8c for the biomass (X). In these figures, the species dynamics plotted for the optimal FBR operation (black curve 2, i.e., the model predictions) are compared to those corresponding to the nominal, nonoptimal FBR operation (blue curve 1 of Maria [46]) and with the experimental results (blue points) of Chen [71]. Both operating policies (optimal 1 and nonoptimal 2) are obtained using the same modified E. coli T5 strain of Chen [70,71].
By analyzing the resulting FBR optimal operating policy (plot no. 2 in Figures 4-8) compared to the suboptimal (nominal) operation of Chen [71] (plot no. 1 in Figures 4-8), several observations can be derived, as follows: By using the same FBR operated under nominal (nonoptimal) conditions of Table 1, the modified E. coli T5 strain reported a higher GLC uptake rate and a much higher TRP production compared to the "wild" strain, as revealed by the analysis given in Table 5. Table 5. Efficiency of the modified E. coli T5 strain for GLC uptake and for the TRP production in the tested FBR of ( (*) By following the same optimal feeding policy, a higher productivity can be obtained for larger batch times (not presented here).
The efficiency of the optimally operated FBR (this paper) in TRP production is significantly higher (ca. 20%) compared to the same FBR but suboptimally (nominal) operated (Table 5), even if the same modified E. coli T5 strain is employed in both cases. The same conclusion also results by comparing the TRP final concentrations in the FBR bulk given in Figure 4 for the two operating policies.
The optimal FBR operation reported a similar dilution of the reactor content, as revealed by (Figure 8b).
The substrate (GLC) consumption in (Table 5) was computed using the following relationship: As expected, a higher TRP productivity requires a higher GLC consumption, as is the case when using a modified E. coli T5 strain instead of the "wild" type. As revealed by (Table 5), the GLC consumption is influenced by the FBR operating mode, even if the same cell strain is used. As indicated by our present analysis given in Table 5, the GLC overall consumption for the optimal (variable feeding) FBR operation is roughly similar to that of a nonoptimal (uniform feeding) FBR operation. Not surprisingly, the optimal operating mode requires a slightly lower GLC consumption (ca. 6%) because of its better use during the batch.
The comparative analysis of the glycolytic species dynamics in Figure 5 reveals close trajectories (even quasi-identical for F6P, FDP species), without any accumulation tendency, for both nominal (nonoptimal, curve 1) and optimal (curve 2) FBR operation. By contrast, the intermediate PEP intermediate species is formed in high amounts but then quickly consumed in the subsequent TRP synthesis, thus tending to reach a QSS. Such a more intensive GLC import for the optimal FBR operation (curve 2) and its successive transformation over the glycolysis pathway and TRP-operon expression are reflected by a higher ATP consumption compared to the nonoptimal FBR operation. The PYR metabolite is consumed in the TCA cycle and excreted in the bulk phase (fairly predicted by our kinetic model matching the experimental data).
The comparative analysis of the TRP operon expression species dynamics in Figure 6 reveals very close trajectories, except for the excreted TRP, for both nominal (nonoptimal, curve 1) and optimal (curve 2) FBR operation. Such a result can be explained by the operon expression mechanism, involving a tight control via its inhibition terms presented in Table 4.
The comparative plots of the GLC concentration dynamics in the FBR bulk phase are presented in Figure 7. They indicate similar decreasing trajectories for both investigated FBR operating alternatives, i.e., (i) nominal (nonoptimal, curve 1) and optimal (curve 2). Such a result can be explained by the same GLC uptake mechanism of the modified E. coli T5 strain. In the optimal case (curve 2), the GLC consumption is higher, due to a higher TRP productivity. The curve 2 unevenness is linked to the variable feeding with GLC of the optimally operated FBR (see the feeding plots in the top part of (Figure 7)).
The comparative plots of the biomass dynamics in the FBR bulk phase are presented in Figure 8c. They reveal similar increasing trajectories for both investigated FBR operating alternatives, i.e., (i) nominal (nonoptimal, curve 1) and (ii) optimal (curve 2). In the optimal operation case, the biomass growth is more intense, due to a significantly higher GLC uptake and a better GLC use during the batch, thus offering more favorable conditions for the biomass growth.
The TRP concentration dynamics in the bulk phase is plotted in Figure 4 for both investigated FBR operating alternatives, i.e., (i) nominal nonoptimal operation of Table 1 (curve 1) and the experimental data (•, blue) of Chen [71], and (ii) optimal operation (curve 2). The TRP higher final concentration leads to a higher productivity for the optimally operated FBR (see observation 2 above). Such a result proves that the optimal time stepwise FBR feeding (i.e., the GLC feeding curve 2 in Figure 7) and the feed flow rate policy of Figure 8a) is superior to the nonoptimal uniform feeding, leading to a better GLC use, even if the overall GLC consumption (see observation 4 above) is similar for both nominal and optimal FBR operation. The better GLC use for the optimal FBR operation is also proven by the less produced secondary metabolite PYR in Figure 5 (curve 2) and by the smaller QSS concentration of the PEP intermediate ( Figure 5, curve 2), quickly transformed into the final product TRP.

Conclusions
The extended bilevel (hybrid) kinetic model adopted in this paper was proven by Maria [46] to adequately represent the dynamics of an experimentally studied FBR under a nominal (uniform feeding) operating policy, for both macroscopic state variables and for the cell key species of the CCM reaction modules related to the TRP production in the FBR, i.e., [a] glycolysis, [b] ATP recovery system, [c] TRP operon expression, and biomass [X] growth. The hybrid structured model, linking the macro state variables to the nano cell-scale variables, was validated using the recorded data from the lab-scale FBR over a long batch time (63 h).
By adopting this adequate kinetic model, the paper exemplifies how the use of reduced CCM-based hybrid kinetic models, of modular construction, including the inter-connected complex metabolic pathways of interest, is a continuously challenging subject when developing structured cell simulators for various engineering applications, such as (a) metabolic flux analysis under variable operating conditions, (b) target metabolite synthesis optimization by optimizing the bioreactor operation, and/or by modifying the cell strain, (c) in silico reprogramming of the cell metabolism to design GMOs (not approached here), (d) a quick analysis of the cell metabolism, leading to an evaluation of substrate utilization, oscillation occurrence, and reactor QSS conditions or structured interpretations of the metabolic changes in modified cells or in direct connection to the bioreactor operation mode, (e) bioreactor/bioprocess optimization (the present study), (f) to derive simple lumped models, locally valid (in the operating parameters domain), and (g) to allow more robust extrapolations of the bioprocess behavior (not tested here).
The engineering evaluations developed in this paper can be further extended, for instance, by deriving a multi-objective optimization of the FBR operating policy, by accounting for not only maximization of the product (TRP), but also minimization of the substrate consumption and of the batch time. The engineering evaluations using such extended bilevel hybrid models present a higher prediction accuracy compared to the simple global (empirical) dynamic models.
Simulations and the experimental checks proved the advantage of using the modified E. coli T5 strain culture to improve the TRP production. The obtained results also proved that, in addition to the cell phenotype characteristics (linked to GLC uptake and glycolysis), the FBR operation mode is the major factor determining the TRP synthesis efficiency. Funding: This work did not receive any kind of funding under any circumstances from any source.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
Restrictions apply to the availability of the experimental data. These data were obtained from a third party (see the acknowledgement) and are available from the authors of the works [70][71][72][73][74] with their permission.

Acknowledgments:
The authors are deeply grateful to An-Ping Zeng and to their research group (Minliang Chen, Lin Chen, and Chengwei Ma) from the Hamburg University of Technology (Institute of Bioprocess and Biosystems Engineering, Germany) for providing on a free basis and without any claim of any kind (including no claim of authorship) experimental data from their published studies on TRP synthesis in a fed-batch bioreactor using a genetically modified E. coli (T5 strain) cell culture [70][71][72][73][74]. In particular, the authors want to thank An-Ping Zeng for the constructive comments and discussions on the manuscript.