Application of (bio) chemical engineering concepts and tools to model genetic regulatory circuits, and some essential central carbon metabolism pathways in living cells. Part 4. Applications in the design of some Genetically Modified Micro-Organisms (GMOs

Figure 1: Time-dependent mercury transfer between the E. coli cell and the three-phase fluidized bed bioreactor (TPFB), according to the proposed hybrid dynamic model HSMDM of Maria [59-62]. Enzyme concentrations in E. coli determine the apparent reaction rates in the bioreactor model (especially PT lumped permease, and PA lumped reductase), while the reactor state-variables (e.g. [Hg²⁺]L, nutrients) determine the cell metabolism and mer-operon expression adaptation.
E. coli cell model notations: The simplified GRC pathway of the mer-operon expression for the mercury ion uptake in E. coli cells includes 7 gene expression lumped regulatory modules (GERMs), by which 4 of type [G(PP)1] (see the Part-1 of this work) that is: 2 modules for mediated transport of [Hg²⁺]env into the cytosol (catalysed by PT) and its reduction (catalysed by PA); 5 regulatory modules of mer operon expression including successive synthesis of the corresponding proteins, that is: PR [ the transcriptional activator of other protein syntheses; its synthesis being triggered by the import of the mercury ions immediately linked as Hg(SR)2 into the cell]; the lumped PT permease; the lumped PA reductase, and of the control protein PD. One additional GERM regulatory module of the simplest type [G(P)1] (see Part 1 of this work) deals with the lumped proteome P and genome G replication into the cell (by thus mimicking the cell ”ballast”). The mer-operon GRC is placed in a growing cell, by mimicking the homeostasis, and cell response to stationary and dynamic perturbations in the environmental [Hg²⁺]env. The reductant NADPH and RSH are considered in excess of the cell. Notations: P = lumped proteome; G = lumped genome; NutG, NutP = lumped nutrients used for gene and protein synthesis, respectively; P• = proteins; G• = genes; RSH = low molecular mass cytosolic thiol redox buffer (such as glutathione). Perpendicular arrows on the reaction path indicate the catalytic activation, repressing, or inhibition actions. The absence of a substrate or product indicates an assumed concentration invariance of these species; ⊕ / Θ positive or negative feedback regulatory loops.

Figure 2: Simplified scheme of the three-phase fluidized bed reactor (TPFB). Notations: 1 - mercury ion fed solution including nutrients (C/saccharides, N/ureea, P/salts, mineral sources, pH-buffer additives, anti-bodies, etc.); 2 - liquid outlet; 3 - Hg vapour in air flow; 4 - sterile air input; 5 - immobilized bacteria on porous support; 6 - air bubbles.
The bioreactor performance depends on the following running parameters: suspended biomass concentration and efficiency; porous support size, feed flow rate, inlet [Hg²⁺], nutrients, additives, aeration rate, pH, etc. In turn, the biomass efficiency (mercury overall reduction rate) depends on the mer-operon enzymes, especially PT, and PA concentration (dependent on the cell resources/phenotype, and environmental [Hg²⁺]).

Figure 3: The present case study - In-silico design of a cloned E. coli with a maximized capacity of mercury uptake from wastewaters [9,59-62,88]. [Down-right] Prof. G. Maria and late Prof. W. Deckwer at the 2nd Croatian-German Conference on Enzyme Reaction Engineering, 21-24 Sept. 2005, Dubrovnik (Croatia), sharing opinions about the bioprocess of mercury removal from wastewaters by using cultures of cloned E. coli cells.

Figure 4: The reduced cell model, accounting for the mer-operon expression (5 GERMs), coupled with 2 kinetic modules referring to the main enzymatic reactions. See also the (Figure 1) caption. Adapted from [5,6,59-62].

Figure 5: Some applications of GRC models [1,2,5,6,92].

Figure 6: The macroscopic model of the three-phase fluidized-bed bioreactor (TPFB) with suspended immobilized E. coli on pumice beads. The Michaelis-Menten mercury reduction model rate constants depend on the mer-plasmid level [59-62].

Figure 7: Three-phase fluidized bed (TPFB) reactor parametric sensitivity related to the variations in the inlet [ $H{g}_{{}^L}^{2+}\underset{x\to \infty }{\lim }$ ]in, leading to variations of the following operating variables: outlet [ $H{g}_{{}^L}^{2+}$ ](top-left); outlet [ $H{g}_{{}^L}^0$ ](down-left); conversion (top-right); and outlet [ $H{g}_{{}^L}^{2+}$ ](down-right). Nominal conditions: [Gmer] = 3 nM; $F_L$ = 0.02 L/min; biomass $c_X$ = 1 g/L; particle size $d_p$ = 1 mm. The notation “in” refers to the bioreactor inlet. Simulations are made using the HSMDM structured kinetic model of the TPFB bioreactor.

Figure 8: Three-phase fluidized bed (TPFB) reactor parametric sensitivity related to inlet liquid flow rate F_L variations, leading to variations of: outlet [ $H{g}_{{}^L}^{2+}$ ](top-left); outlet [ $H{g}_{{}^L}^0$ ](down-left); $H{g}_{{}^L}^{2+}$ conversion (top-right); outlet [ $H{g}_{{}^G}^0$ ](down-right). Nominal conditions: [Gmer] = 3 nM; [ $H{g}_{{}^L}^{2+}$ ]in = 20 mg/L; biomass C_X = 1 g/L ; particle size d_p = 1 mm. Simulations are made using the HSMDM structured kinetic model of the TPFB bioreactor.

Figure 9: Three-phase fluidized bed (TPFB) reactor parametric sensitivity related to biomass concentration C_X variations, leading to variations of: outlet [ $H{g}_{{}^L}^{2+}$ ](top-left); outlet [ $H{g}_{{}^L}^0$ ](down-left); $H{g}_{{}^L}^{2+}$ conversion (top-right); outlet [ $H{g}_{{}^G}^0$ ](down-right). Nominal conditions: [Gmer] = 3 nM; F_L = 0.02 L/min; [ $H{g}_{{}^L}^{2+}$ ]in = 20 mg/L; particle size d_p = 1 mm. Simulations are made using the HSMDM structured kinetic model of the TPFB bioreactor.

Figure 10: Three-phase fluidized bed (TPFB) reactor sensitivity related to particle size d_p variations, leading to variations of: outlet [ $H{g}_{{}^L}^{2+}$ ](top-left); outlet [ $H{g}_{{}^L}^0$ ](down-left); $H{g}_{{}^L}^{2+}$ conversion (top-right); outlet mercury in the outlet-gas [ $H{g}_{{}^G}^0$ ](down-right). Nominal conditions: [Gmer] = 3 nM; F_L = 0.02 L/min; [ $H{g}_{{}^L}^{2+}\frac{1}{2}$ ]in = 20 mg/L; biomass C_X= 1 g/L. Simulations are made using the HSMDM structured kinetic model of the TPFB bioreactor.

Figure 11: The nano-scale cellular enzymatic process in the immobilized E. coli bacteria: mer-operon (4 gene lumps) expression (Figure 4), and the main enzymatic reactions (mercury permeation, and its reduction). Adapted from [5,59-62].

Figure 12: The main characteristics of the E. coli (K-12 strain) used in the present case study. Data from EcoCyc [103].

Figure 13: Typical evolution of relevant species concentrations predicted by the E. coli cell model, after a “step” perturbation in the TPFB bioreactor inlet from []s = 0.1 to 10 µM (ca. 2 mg/L), for the case of cell cloned with [Gmer]= 3 nM (___), or with [Gmer]= 140 nM (• • • • • ). The arrow indicates the quick and vigorous response of the two key mer-enzymes: PT = the mercury ions permease, and PA = the mercury ions reductase.

Table 1: Unstructured (apparent, reduced, global) kinetic model of Michaelis-Menten type for mercury ions reduction by E. coli (after Philippidis, et al. [88-90]). Notations: substrate S = $H{g}_{{}^L}^{2+}$ = $H{g}_{env}^{2+}$ ; PT = lumped permease for membranar transport of $H{g}_{{}^{env}}^{2+}$ ; PA = reductase of cytosolic mercury ions; RSH = low molecular-mass thiol redox cytosolic buffers; NADPH = nicotinamide adenine dinucleotide phosphate; ‘env’ = environmental; ‘cyt’ = cytosol. Adapted after [59-62].

Table 2: Nominal operating conditions and characteristics of the semi-continuous (SCR), three-phase fluidized bed (TPFB) reactor used for mercuric ions reduction (uptake) by the immobilized E. coli cells on pumice porous support (Footnote a). Notations are as provided in the Abbreviations and Notations chapter.
Footnotes:
(a) fluid physical properties correspond to those of water and air at the operating conditions.
(b) see the calculation rule of Table 3.
(c) particle tortuosity was evaluated by using the Shen and Chen [111] formula .

To simulate the (semi-)continuous TPFB bioreactor performance, a simple unstructured dynamic ideal model was considered (Table 7), as a first approach [98,99], by assuming homogeneous (perfectly mixed) liquid and gas phases, and a uniform distribution of the solid particles (of uniform characteristics) in the G-S-L fluidized bed. This reduced TPFB bioreactor dynamic model accounts for only the apparent (Michaelis-Menten) mercury uptake kinetics (denoted by MM1, MM2, PHM) of Deckwer et al. [93] (Table 8). For such global kinetics, the apparent mercury reduction rate (r_app) necessary in the bioreactor model (Table 7) is evaluated following the (Table 8) rule. The apparent rate $H{g}_{{}^{cyt}}^0$ was evaluated by solving on every integration time-increment the quasi-steady-state equality of mass fluxes at the solid-liquid (S-L) external interface (on the liquid side), by also including the external diffusion coefficient KSaS [100]. This rule also includes the diffusional resistance of the substrate (mercury ions) / product (dissolved metallic mercury) transport through the support pores.

If the apparent (Michaelis-Menten) mercury uptake kinetics of Philippidis et al. [88-90] (Table 1) are introduced in the TPFB bioreactor model, then a better unstructured global model is obtained (Table 3). By also including the [Gmer]-dependent rate constants, it results in the apparent mercury reduction rate (r_app) by the immobilized bacteria of concentration c_X in (Table 1) in the spherical solid carrier [of fraction εs in Table 2 (pumice)] for different levels of [Gmer] plasmids. The mercury mass balance in the liquid and gas phases is presented in (Table 7). This reduced (unstructured) TPFB reactor dynamic model (Table 7, and Table 8) includes terms referring to the mass balance in the bulk (liquid, L) phase, the gas phase (G), the interphase L-G transport (rtrans) of the volatile mercury, and the overall/apparent bioprocess of mercury reduction inside the solid particles (r_app). In such a way, the apparent mercury reduction rate also includes the diffusional resistance of the substrate (mercury ions) / product (dissolved metallic mercury) transport through the pumice support pores.

A hybrid reduced/unstructured dynamic model of the TPFB can by constructed by linking the macro0scale state-variables of the TPFB (Table 7), with the Michaelis-Menten kinetics of Philippidis et al. [88-90] (Table 1). This apparent model (Figure 6) can only give a rough idea of the bioreactor dynamics, but it is unable to describe the biomass adaptation to environmental changes, that is variations in both inlet feed flow-rate and inlet mercury load in the influent. To offer a prediction to such engineering requirements, a structured kinetic model of the mercury uptake in the E. coli bacteria at a cellular level is necessary. The next chapter describes the complex HSMDM structured model proposed by Maria [59-61], by considering the macro-level TPFB model linked to the cell biological process included by using the WCVV modelling approach (Part 2 of this work) to simulate the dynamics of the cell mer-operon expression self-regulation in the wild, or modified E. coli under various environmental (bioreactor operating) conditions.

The resulting unstructured TPFB global dynamic model of (Table 7), while keeping the apparent M-M kinetics of (Table 8) only for the mercury cytosolic reduction step, allows a rough simulation of the transient operating conditions of the TPFB bioreactor. This apparent model of (Figure 6, plus Table 8), even if suitable for quick/rough engineering evaluations, is unable to describe the biomass adaptation to environmental changes, that is variations in the both inlet feed flow rate and inlet mercury load in the influent, or the effect of cloning E. coli cells on the bioreactor efficiency.

Extended HSMDM model of the TPFB bioreactor

To avoid the above-mentioned limitations of the unstructured / reduced dynamic model of the TPFB bioreactor, an extended and complex HSMDM structured model was elaborated by also including the effect of cloning E. coli with various concentrations of mer-plasmids [Gmer], by also including the Michaelis-Menten [Gmer]-dependent kinetics of Philippidis et al. [88-90] in the cytosolic mercury reduction step (Table 4).

To offer a prediction to such engineering requirements, an extended HSMDM structured kinetic model of the mercury uptake in the E. coli bacteria at a cellular level is necessary (of WCVV type, see Part-2 of this work). And, of course, this cell-state-variables (nano-scale) part of the HSMDM should be linked to the TPFB bioreactor state-variables (macro-scale), because they are interrelated.

This section describes the WCVV-type structured cell model proposed by Maria [59-61] to simulate the dynamics of the mer-operon expression induction by the presence of environmental mercury, and the expression self-regulation in the “wild”, or in-silico design GMO E. coli under various environmental (bioreactor operating) conditions.

This HSMDM structured and modular cell kinetic model has been developed by Maria [59-61]. Roughly, the part of the extended HSMDM structured model referring to the mer-operon expression WCVV model includes a GRC of 7 GERMs (Figure 1) linked by following the rules described in this section, and by accounting for the few experimental information of Philippidis et al. [88-90], and of Barkay et al. [91]. The derived GRC model can simulate the dynamics of the mer-operon expression, and the process self-control at a molecular level under isothermal and isotonic conditions.

The extended HSMDM structured bioreactor bi-level model includes two inter-connected parts: (a) The macro-level state-variables of the TPFB bioreactor (Table 3), that is the mass balance of $H{g}_L^{2+}$ , $H{g}_L^0$ in the liquid-phase, and of $H{g}_G^0$ in the gas-phase; (b) The cell nano-level state-variables, that is the mass balance of the cell key-species included in the WCVV-GRC model of (Table 4), completed with the [Gmer]-dependent apparent Michaelis-Menten kinetics of (Table 1) for the mercury reduction reaction in the cell cytosol. All the cell parameters included in the HSMDM model correspond to the E. coli characteristics (Table 5), and (Figure 12).

The macro-scale state variables

More specifically, the mercury differential mass balance, in the TPFB model of (Table 3) includes the following terms:

1. The apparent mercury reduction rate, evaluated at the solid interface;
2. The substrate (that is S =  $H{g}_{{}^L}^{2+}$ = $H{g}_{env}^{2+}$ ) diffusional transport in the particle, expressed by the effectiveness factor (h) evaluated using the Thiele modulus for a Michaelis-Menten type reaction [101], the effective diffusivity (DS,ef) accounting for the molecular diffusion (DS,L), the particle porosity (ep) and the tortuosity (t) (other resistances being neglected [102]);
3. The apparent rate rj, app was evaluated by solving on every integration time-increment the quasi-steady-state equality of mass fluxes at the solid-liquid (S-L) external interface (on the liquid side), by also including the external diffusion coefficient kSaS [100]. This rule also includes the diffusional resistance of the substrate (mercury ions) / product (dissolved metallic mercury) transport through the support pores.
4. The liquid-to-gas transport rate rtrans of the metallic mercury ( $H{g}_{{}^L}^0\text{to }H{g}_{{}^G}^0$ ), evaluated from the quasi-steady-state equality of mass fluxes at bubbles G-L interface, by accounting for the mass transfer coefficients kLaL (on the liquid film side; experimental value adopted), and kGaG (on the gas film side; evaluated from the Sharma’s relationship given by Trambouze et al. [99] ).

The cell-scale state variables

This paragraph describes the WCVV structured cell model proposed by Maria [59-61] to simulate the dynamics of the mer-operon-induced expression, and its self-regulation under various environmental conditions. This model was constructed and validated by using the experimental data of Philippidis et al. [88-90], and the experimental information of Barkay et al. [90] on the mer-operon characteristics. The WCVV cell model was built up by accounting for the GERMs library, their P.I.-s, and the linking rules presented in Part 2 of this work.

The proposed E. coli cell model by Maria [59-61] includes the GRC responsible for the control of the mer-operon expression and the whole process of mercury ions removal. The proposed GRC includes 4 lumped genes (denoted by GR, GT, GA, GD in (Figure 1, Figure 2, Figure 3, and Figure 4) of individual expression levels induced and adjusted according to the level of mercury and other metabolites into the cytosol. As it follows from these figures, the mer-operon expression process is induced by the presence of a small concentration of mercury ions in the environment, leading to the appearance in the cytosol as Hg(SR)₂ compounds. The whole process is tightly cross- and self-regulated to hinder the import of large amounts of mercury into the cell, which eventually might lead to the blockage of cell resources (RSH, NADPH, other metabolites, and proteins), thus compromising the whole cell metabolism. The GRC model includes four GERMs of simple but effective [G(PP)1] type (see the discussion of Part-2 of this work) as follows (see Figure 11, and Table 4):

1. A GERM to regulate the Hg²⁺ transport across the cellular membrane, mediated by three proteins (PmerP, PmerT, and PmerC) from the periplasmatic space, considered as a lumped permease PT in the model. Phillippidis et al. [88-90] found this transport step to be energy-dependent and the rate-determining step for the whole mercury uptake process. Once the mercuric ion complex arrives in the cytosol, thiol redox buffers (such as glutathione of millimolar concentrations) form a dithiol derivative Hg(SR)2. Instantly, the GT lumped gene expression is induced by the regulatory protein PR, and easily ‘smoothed’ by the large ‘ballast’ effect of the proteome lump P. (see rule 6 of chapter. 2.3.6 of Part-2 of this work).
2. A GERM to control the expression of the PR protein that induces and controls the whole mer-operon expression in the presence of cytosolic  $H{g}_{{}^{cyt}}^{2+}$ (even if they are present in traces, that is low nM concentrations, (Table 5) and Figure 12). This GERM acts as an amplifier of the mer-expression leading to a quick (over ca. 30 s) cell response by starting the mer-enzymes production. First, the expression of the GR gene leads to obtaining the encoded PR protein. The process magnitude is also controlled by the protein PD, present in small amounts in the cell.
3. A GERM to control the expression of PA enzyme responsible for the Hg(SR)2 (that is $H{g}_{cyt}^{2+}$ ) reduction to metallic mercury ( $H{g}_{{}^{cyt}}^0$ ) into the cytosol. The metallic mercury is relatively non-toxic for the cell, being easily removable through membranar diffusion into the bioreactor bulk liquid phase ( $H{g}_{{}^L}^o$ ). From here, metallic mercury ( $H{g}_{{}^L}^o$ ) passes into the gas phase (air bubbles) $H{g}_{{}^G}^o$ , being entrained and continuously removed by the sparged air as mercury vapours ( $H{g}_{{}^G}^o$ ), to be later recovered. The GA gene expression is induced and controlled by the PT protein, whose expression is controlled by the PR level, which in turn is controlled by the cytosolic mercury and PD levels.
4. A GERM controlling the protein PD synthesis. This protein has a complex role, by maintaining a certain level of GR expression even when the mercury is absent in the cytosol [91]. It also hinders the over-expression of GR and GT when too large concentrations of mercury  $H{g}_{{}^{cyt}}^{2+}$ are present in the environment.
5. A GERM controlling the replication of the lumped cell proteome (P) and genome (G) (of concentrations 107 nM, and 4500 nM, respectively, to mimic the effect of “cell ballast”, see Part-2 of this work) in the immobilized E. coli cells. These data are based on the Ecocyc [103] databank (Table 5, and (Figure 12), thus mimicking the cell ’ballast’ effect on the cell genes expression, and on all considered reactions. The need to include the cell content lump (the so-called ’cell ballast’) in the WCVV model is legitimate by the possibility offered by such a structured cell model to reproduce the smoothing effect of perturbations leading to more realistic transient times (compared to a cell with a ‘sparing’ content), the synchronized response to certain inducers, and the ‘secondary perturbation’ effect transmitted via the cell volume to which all cell components contribute (see the discussion in the Part-2 of this work).
6. Thus, the whole mer-operon expression consists of a controlled expression in cascade (Figure 11) of all 4 mer-genes {GR, GT, GA, GD in (Figure 1, Figure 2)}.

In total, the cell GRC dynamic model describing the mer-operon expression includes only 26 individual or lumped cellular species involved in 33 reactions (Figure 1, Figure 11, and Table 4). The structured cell WCVV model is presented in the (Table 4). The cell model is coupled with the SCR –TPFB dynamic model (Table 3) through the rj,app link to the individual reactions $r_{j,(Cj,S)}$ occurring inside the cell. The WCVV cell model includes not only the reactions and the dynamics of the mer-GRC, but also the enzymatic reactions directly responsible for the environmental mercury $H{g}_{env}^{2+}$ (also denoted by $H{g}_L^{2+}$ ) import into the cell as $H{g}_{cyt}^{2+}$ , and for its reduction to cytosolic metallic mercury $H{g}_{cyt}^0$ in the E. coli cells cloned with a defined [Gmer] concentration of mer-plasmids.

All reactions in the cell model of (Table 4) are considered elementary, except some of them for which extended experimental information exists, that is [59-61]:

1. A Michaelis-Menten rate expression for the mercuric ion permeation through the membrane into the cell;
2. A Michaelis-Menten rate expression for the mercuric ion reduction in the cytosol;
3. A Hill-type quick induction of the GR expression that can rapidly initiate the production of permease PT (through the control protein PR) when mercuric ions are present in large amounts.
4. Dimmerization reactions of TF-s are considered to be much more rapid than the enzyme synthesis, while equal concentrations of active gene (G) and inactive (GPP) forms of the generic gene G are considered at homeostasis to maximize the GERM efficiency (see Part-2 of this work) (Figure 1, and Figure 11).
5. The lumped proteome P, present in a large amount, is included in all gene expression rates, thus leading to a more realistic evaluation of the GERM regulatory efficiency indices P.I.-s (see Part-2 of this work) [1,2,3,4,57,58].
6. The WCVV model equations for the mercury uptake in E. coli, together with the general hypotheses of the WCVV model are presented in (Table 4). The cell is considered an open system, of uniform content and negligible inner gradients. To not increase the number of parameters, the structured model includes GERM of the simplest form [G(PP)1] (see Part-2 of this work) for all mer-genes, by using dimeric TF-s (of Protein: Protein type) to increase the GERM regulatory efficiency, as experimentally proved by [1,2,3,38,55,58]. The resulting HSMDM model includes 26 individual or lumped cellular species and 33 reactions. All reactions are considered elementary, except some of them for which extended experimental information exists, that is the Michaelis-Menten rate expression for mercury permeation into the cell, and its reduction in cytosol. A Hill-type induction of the GR expression is adopted to rapidly amplify the mer-operon expression when mercuric ions are present in significant amounts inside the cell. Dimerization reactions of TF-s are considered to be much more rapid than the enzyme synthesis, while equal concentrations of active (G(i)) and inactive (G(i): TFTF) forms of the generic gene G(i) are considered at homeostasis to maximize the GERM efficiency [1,2,3,4,55]. The homeostatic characteristics of E. coli cells (belonging to a uniform culture) from the reactor, and the adopted species concentrations are presented in (Table 5).
7. The model rate constants are estimated from solving the cell stationary mass balances for nominal concentrations of observable species (see Part-2 of this work), but also from optimizing the GERM regulatory indices (see Part-2 of this work), for instance: (1) by adjusting the optimum TF level of gene expression to obtain the minimum recovering times after a 10% dynamic perturbation in one of the key species, and (2) realize the smallest sensitivities of the homeostatic species concentrations vs. external perturbations (see Part-2 of this work for an extended discussion). The E. coli cell characteristics and the considered cell key-species homeostatic concentrations used in the rate constants estimation are given in (Table 5). Exceptions are the Michaelis-Menten (M-M) rate constants for the mercury transport and its reduction in cytosol adopted from the Phillipidis et al. [88-90] kinetic data. By adopting this M-M reduced model, the rate constants of these two metabolic steps depend on the amount of [Gmer] plasmids in the cloned E. coli cells (Figure 11, and Figure 6). Simple correlations are used to include this essential aspect in the model.

HSMDM structured model advantages

As extensively discussed in Part 2 of this work, and proved in this section, such a cell WCVV structured model presents a large number of advantages, being able to:

1. Simulate the cell metabolism adaptation when the environmental mercury level changes. Such a reconfiguration of the levels of mer-genes and mer-proteins is presented in (Figure 13) as a step response after a ‚step’-like perturbation in the mercury level from []s = 0.1 µM to 10 µM (ca. 2 mg/L), in a cloned E. coli cell with mer-plasmids of [Gmer]= 140nM, compared to a cell cloned with only [Gmer]= 3 nM. The transient state toward the cell’s new homeostasis of the “adapted” mer-gene/protein levels stretches over 15-20 cell cycles (of ca. 0.5 h each) as long as the environment stationary step perturbation is maintained. An E. coli cell with a higher content of mer-plasmids reacts much strongly to the environmental perturbations, by quickly starting to produce the enzymes responsible for the mercury removal.
2. Because the Hg²⁺ reduction rate constants are dependent on the mer-plasmid level in the cell, the WCVV GRC model can predict the maximum level of mer-plasmids that can be added to the cell genome for improving its mercury uptaking capacity without exhausting the internal cell resources, thus putting in danger the cell survival. Consequently, it follows that this cell MSDKM model allows the in-silico design of modified E. coli cloned with a suitable amount of mer-plasmids to improve its efficiency in cleaning wastewater by improving the mercury uptaking capacity. As an example, (Figure 14 presented the cell key-species stationary levels, and the mercury concentration in the TPFB bioreactor bulk-phase for two GMOs: An E. coli cloned with [Gmer] = 67 nM, and another culture of E. coli cells cloned with [Gmer] = 140 nM. Simulations of Maria [59-62] revealed that as the mer-plasmids level increases, the mercury uptake capacity also increases. However, an upper limit exists (around 140 nM) over which the cell resources will be exhausted, putting its metabolism in danger.
3. The simulated state-variables plotted in (Figure 14), for using the unstructured global model (Table 7 and Table 8) comparatively with using the structured extended HSMDM model, reveal a higher quality (vs.- experimental data), and more detailed/refined predictions for a larger number of state-variables.
4. By coupling the structured MSDKM cell model of (Table 4) with the three-phase (TPFB) continuous bioreactor model of Table 3 (with immobilized E. coli cells on pumice beads; see also the bioreactor model in Figure 6), Maria et al. [59-62] have been able to determine the optimal operating policies of the bioreactor in relationship to the culture of cloned E. coli cells. Similar studies are reported by [7,59-62,88].

Quick comparison between the structured HSMDM and the unstructured/global dynamic model of the TPFB reactor:

The superiority of the structured extended HSMDM model (Table 3 and Table 4) against the unstructured reduced model (Table 7 and Table 8) was used to simulate the TPFB dynamics over a wide range of time.

Thus, when simulating the TPFB reactor performance (of nominal conditions given in Table 2, for the biomass immobilized on pumice support) by employing a classical unstructured reactor model, the biomass adaptation to variable input loads is not accounted for, that is the maximum rates v_m,t, v_m,P, and the Michaelis inhibition constants K_mt, K_mP and K_iP of (Table 1) are kept constant (eventually depending only on the Gmer plasmid level). In fact, v_m,t and v_m,P rate constants depend on the cell enzyme levels of PA (reductase) and PT (permease) (Figure 1), which vary during the bacteria adaptation in the bioreactor under stationary or transient conditions.

To solve this problem by also offering a more detailed and robust prediction on the system dynamics/cell metabolism evolution, more detailed (structured) dynamic models are necessary for the TPFB bioreactor, linked to the structured kinetic model of mer-GRC in E. coli cell. Thus, it results in the here-described hybrid structured HSMDM model. The two linked differential models are solved together simultaneously, by applying a mutual exchange of input/output parameters {[ $H{g}_{env}^{2+}$ )], [ $H{g}_{env}^0$ ], PT, PA} on every small time increment throughout the solution (integration) of the HSMDM model, as graphically represented in (Figure 1).

As described in the below paragraph dedicated to solving the HSMDM model of (Table 3 and Table 4)., due to its structure, only the HSMDM model can simulate the E. coli cell response to dynamic or stationary perturbations from the environment (i.e. the TPFB bulk-phase). And also the influence of [Gmer] plasmids on the bioreactor performances. For instance, (Figure 13) displays the E. coli cell adaptation after a ‘step’-like perturbation in the environmental [ $H{g}_{env}^{2+}$ ]s (that is in the bioreactor bulk-phase) from the background level of 0.1 µM to 10 µM (ca. 2 mg/L), for the case of cell cloned with [Gmer]= 3 nM (full line), or with [Gmer]= 140 nM (dash line) mer-plasmids. The transient state toward the cell’s new homeostasis usually lies over 15-20 cell cycles as long as the environmental stationary perturbation is maintained. The simulated species dynamic trajectories plotted in (Figure 13) reveal a vigorous response of the cell mer-species (especially of PT, and PA) to the perturbation in the environmental [ $H{g}_{env}^{2+}$ ]s (i.e. liquid bulk-phase). Such cell responses to environmental perturbations are impossible to reproduce by the simple unstructured reduced model (Table 7 and Table 8). A more detailed comparison between the structured HSMDM model, and the unstructured/global one, for this case study is made by Maria [3,4].

The link between the macro-, and cell-scale state variables when solving the HSMDM model

The link between the two parts of the macro-sate-variables and the cell-state-variables of the HSMDM is made by the [ $H{g}_{L,env}^{2+}$ ] controlling the [ $H{g}_{L,cyt}^{2+}$ ], which, in turn, will control the whole mer-operon dynamics and the mercury reduction process through the cascade expression of the key-proteins {PR, PT, PA, PD} (see Table 4). In turn, the cell bioprocess induced and intensified by the [ $H{g}_{L,env}^{2+}$ ] level, controls the TPFB bioreactor output [ $H{g}_G^0$ ] through [ $H{g}_L^0$ ], via [ $H{g}_{cyt}^0$ ].

The link between the bioreactor macro-state variables and the cell nano-scale-state variables can be better pointed out when describing the HSMDM DAE model solution, that is its integration over a large time domain.

Due to its high complexity, the solving rule of the HSMDM model involves successive integration steps with an adopted time-interval equal to the cell cycle (ca. 30 min), enough to obtain the steady-state of the TPFB reactor on every time-interval, as follows: (i) The rule starts with solving the extended hybrid HSMDM E. coli cell model by using the known initial condition (the cell variables’ initial state, or those from the end of the previous integration cycle), and by considering the current concentration of [ $H{g}_{env}^{2+}$ ] in the bioreactor liquid phase (the nutrients are considered in excess and of constant levels). Thus, the cell species dynamics over one cell cycle are obtained from solving the WCVV cellular model. (ii) Then, the TPFB reactor model is solved over the current time-interval by using the known initial conditions (i.e. the reactor state variables from the end of the previous time-interval), and by considering the enzymes PT and PA concentrations resulting from the E. coli model solution. The [PT] and [PA] are necessary for setting the maximum reaction rates $v_{m,t}=k_7{c}_{PT}\text{ and }{v}_{m,P}=k_8{c}_{PA}$ in the reactor model. The biomass level on the support is taken constant in the simulated case study, but an additional mass balance can be easily added to the reactor model if a kinetic model about biomass (X) growth is available.

The procedure is repeated a large number of times, over hundreds of cell cycles. For instance, (Figure 13) displays the E. coli cell adaptation after a ‘step’-like perturbation in the environmental [ $H{g}_{env}^{2+}$ ]s (that is in the bioreactor bulk-phase) from the background level of 0.1 µM to 10 µM (ca. 2 mg/L), for the case of cell cloned with [Gmer]= 3 nM (full line), or with [Gmer]= 140 nM (dash line) mer-plasmids. The transient state toward the cell’s new homeostasis usually lies over 15-20 cell cycles as long as the environmental stationary perturbation is maintained. The simulated species dynamic trajectories plotted in (Figure 13) reveal a vigorous response of the cell mer-species (especially of PT, and PA) to the perturbation in the environmental [ $H{g}_{env}^{2+}$ ]s (i.e. liquid bulk-phase). Such cell responses to the environmental (dynamic or stationary) perturbations are impossible to reproduce by the simple unstructured reduced model (Table 7 and Table 8).

Some simulations with the HSMDM extended dynamic model

The resulting hybrid dynamic HSMDM model of (Table 3), including the WCVV kinetic model of the GRC responsible for the mer-operon expression of (Table 4) allows simulating with more accuracy the dynamics of the TPFB bioreactor at both macro- and (call) nano-scale level. By also including in (Table 4) the Michaelis-Menten kinetics of Philippidis et al. [88-90], the above-described extended HSMDM model can also be used to study the parametric sensitivity of the studied TPFB. Thus, the dynamics of the TPFB (with the characteristics of (Table 2(pumice)) have been simulated under the nominal operating conditions of Table 2, but every time-varying one operating parameter, that is:

1. The inlet [ $H{g}_{{}^L}^{2+}$ ]in concentration from the background pollution level (0.02 mg L-1) to successively 1, 5, or 10 mg L-1 levels (Figure 7);
2. The inlet liquid flow rate FL of 0.01, 0.02, and 0.04 L/min (Figure 8);
3. The biomass load cX on the solid support taken as 0.1, 0.25, 1.0 gX/L (referred to as the liquid volume) for a constant fraction of the solid in the reactor (Figure 9);
4. The use of particle average diameter dp of 1 mm, or 4 mm respectively (Figure 10).

Simulations of the TPFB with this extended HSMDM model over ca. 50 min running time reveal that the liquid residence time in the reactor (related to the inlet liquid flow rate, F_L) and the biomass content (c_X) are the most influential operating parameters, being directly responsible for the realised mercury uptake conversion. For instance, by doubling the feed flow rate (from 0.02 to 0.04 L min-1), the uptake conversion can be reduced by ca. 5%. The particle size is also important, with an increase in the average d_p leading to a higher resistance to the diffusional transport in pores and to a particle effectiveness diminishment (from 0.51 for dp = 1 mm to 0.15 for dp = 4 mm under nominal conditions). On the other hand, the biomass average load in the reactor can be adjusted by employing a continuous purge-renewal system of the solid particles. The mercury content in the input flow, [ $H{g}_{{}^L}^{2+}$ ]in in the range of 1-40 mg/L, has a significant influence on the mercury content of the output gas (Figure 7). As further proved, this last parameter is also of tremendous importance for adaptation of bacteria metabolism when using “wild” E. coli cells, or those cloned with an increased copy-numbers of mer-plasmids (denoted by [Gmer]).”

In-silico design of a cloned (GMO) E.coli.

Simulations with the extended HSMDM model of the TPFB dynamics using several E.coli cultures, that is E.coli cells cloned with various levels of mercury-plasmids, reveal interesting conclusions. As plotted in Figure 15, the cellular uptake process (i.e. the stationary reduction rate, r_Hg ) may be improved by increasing the mer-permease (PT) content into the cell (GT/PT dynamics in Figure 14-right), up to a limit of ca. [Gmer] = 80-140 nM mer-plasmids, higher than the permeation rate remains unchanged, thus preventing exhaustion of the cell metabolic resources. Such a conclusion was experimentally confirmed by Philippidis et al. [88-90]. As revealed by the present simulations, for the same reason, the cell regulatory system maintains an upper limit for the membranar transport rate (rt) of environmental $H{g}_{env}^{2+}$ into the cell, irrespectively to the $H{g}_{env}^{2+}$ concentration in the environment, the cytosolic $H{g}_{cyt}^{2+}$ concentration never exceeding a 3500-5000 nM level (Figure 13).

Fitting the HSMDM model parameters from experimental data

To overcome the high computational effort necessary to fit the large number of rate constants of this extended HSMDM model, and to increase their confidence and physical meaning, a three-step procedure was employed based on the available experimental data and additional information from literature, as followings.

1. Mass transport parameters of the TPFB reactor, that is interfacial partial transfer coefficients (k_Sa_S, k_La_L, k_Gp_G), effective diffusivity (D_ef) and particle effectiveness factor (h) have been estimated by using the experimental data of Deckwer et al. [93], and based on common correlations from the chemical engineering literature (Table 3) evaluated for the specified reactor operating conditions of (Table 2- pumice support).
2. Cell model parameters are estimated by using the Philippidis et al. [88-90] kinetic data obtained from separate batch experiments with “wild” (not-cloned) E. coli cells. By using the defined cell nominal characteristics of (Table 5) (some key notations are: tC = cell cycle time; Vcyt,O = born cell volume; C_X = biomass concentration in the bioreactor; P = cell lumped proteome; G = cell lumped genome; NutG, NutP = lumped nutrients for G and P synthesis, respectively; [Gmer] = mer-plasmid levels into the cloned cell. The stationary levels of the essential cell mer-proteins (PA, PR, PD) are taken from the literature data. The involved TFs (PRPR, PTPT, PAPA, PDPD) intermediate concentrations resulted by maximizing the GERMs regulatory efficiency P.I.-s (see Part-2 of this work) [1,2,3,4,21,22]. The resulting rate constants of (Table 4) have been estimated for the most severe experimental conditions of [ $H{g}_{env}^{2+}$ ]s = 120 µM, and [Gmer] = 140 nM. The used first guess of Hill-induction rate constants (nPD = 1, nPR = -0.5, nH = 2, a = 3) have been adopted at values recommended in the literature, by similarity with the Hill-induction of gene expression in genetic switches (see the remarks included in Table 4), while the Michaelis-Menten rate constants of mercury membranar transport (r_max,t,K_mt) and its reduction rates’ constants (r_max,P, K_mP, K_iP) have been kept at the fitted values of Philippidis et al. [88-90] (Table 1). The reference concentration CHg2cy,ref was adopted at the average cytosolic level of mercuric ions detected by Philippidis et al. [88-90]. When applying the model estimation rule, the GERM regulatory indices have been kept at their optimized levels, which corresponds to: (i) Equal concentrations of catalytically active/inactive forms [G(j)]s = [G(j)TFn]s adopted at steady-state to ensure GERM maximum regulatory efficiency vs. perturbations (i.e. smallest sensitivities of the homeostatic levels vs. external perturbations [1-4], and (Part-2 of this work); (ii) adjustable optimum [TF]s level (ca. 4 nM here; Table 5) involved in the gene expressions to get the minimum recovering times after a dynamic perturbation in the key species [1,2,21,22,51,58], (see Part-2 of this work). Other adjustable parameters, such as the cell concentration in biomass (C_cell), are tuned to fit the experimental cellular mercury reduction rate of Philippidis et al. [88-90].
3. The model validity is extended over a wider experimental range, by covering the input [ $H{g}_L^{2+}$ ] = 0-100 mg/L, and plasmid levels of [Gmer] = 3-140 nM. All these have been realised by adjusting some key parameters of the extended HSMDM model, and by applying two concomitant fitting criteria, that is: (1) the cytosolic mercury reduction rate should fit the experimental values of Philippidis et al. [88-90], obtained for ‚permeabilized’ cells (i.e. Min $\Vert r_{Hg,\exp }-r_{8,\mathrm{mod}el}\Vert$ , see the r8 expression of Table 4); (2) the TPFB reactor predicted output CHg,Out by the unstructured model (Table 1, Table 7, and Table 8) with using the Philippidis et al. [88-90] Michaelis-Menten parameters and the mass transfer terms, to fit with those of the structured (cell + reactor) extended HSMDM model of (Table 3, and Table 4). The Hill parameter b=2a⁴ was adjusted by using the following approximate linear dependence on the inlet mercury load:

$a=a_{\min }+(a_{\max }-a_{\min })(c_{Hg,in}-c_{Hg,in}^{\min })\text{ / }(c_{Hg,in}^{\max }-c_{Hg,in}^{\min })$

(see the fitted values a_min and a_max in Table 4 – footnote d). The above-simplified correlation tries to account for the influence of the environmental mercury on the induction characteristics of the mer-operon expression, that is a slow induction in the presence of low levels of Hg²⁺ inducer and a sharp (sigmoidal) mer-operon expression response to high levels of [ $H{g}_{cyt}^{2+}$ ] inducer. When a certain saturation level is reached, the limited cell resources will impose ’flattening’ the metabolic mer-rates and the mer-protein levels, irrespectively to the increased amount of mercury in the environment [88-90].

The extended HSMDM model predictions are in satisfactory agreement with the experimental data. Thus, the model predictions for r_Hg in (Figure 15, A-C) (curves ’2’) roughly fall within the confidence band [‚Eup’, ‚Elow’] of the experimental curves (‚E’) of Philippidis et al. [88-90] for three different cloned cell cases (confidence curves being plotted by taking constant the reported maximum relative error of ca. 19%). The unstructured model (chap.2.3, Table 1, Table 7), and Table 8) predictions (that is curves ’1’) reported apparent r_Hg of low adequacy, that is too low values. Such a result for the unstructured reduced model of the TPFB bioreactor is due to the rough bioprocess representation by the M-M model, and due to the inclusion of the mass transport resistance between the three phases in contact. The extended HSMDM model adequacy in terms of standard error ’Std’ / average observed value ’Obs’ ratio is evaluated for every cloned cell culture, leading to Std/Obs of 22.3%, 16.7%, 20.4%, 19.9%, 18.8% for [Gmer] levels of 3, 67, 78, 124, and 140 nM respectively, that is acceptable values when compared to the maximum experimental error (ca. 25-33% in Figure 15). The predicted structured vs. unstructured model outputs, in terms of outlet concentration of mercury ( $c_{Hg,out}$ ) are also in fair agreement (curves not displayed here), that is Std/Obs values of 3.0%, 4.7%, 3.7%, and 12% for [Gmer] levels of 67, 78, 124, and 140 nM respectively. The model’s poor adequacy for the [Gmer] = 3 nM data set and low [ $H{g}_L^{2+}$ ] in the environment (input) might be explained by the use of a less adapted E. coli cell than probably those studied by Phillipidis et al. [88-90], reflected by a smaller [PA] initial level of 600 nM (see the PA-curve in Figure 15-D). Once a higher level of Gmer-plasmids is introduced into the cell, then, when higher $H{g}_L^{2+}$ stimulus levels are present in the environment, the cytosolic PA level is tripling.

The extended HSMDM structured model cellular rate constants (see the above point (b), and Table 4) are in good agreement with the reported data from the literature, as remarked in the same (Table 4). The fitted rate constants multiplied by the reactant lead to reaction rates of the same order of magnitude as those reported in the literature for similar genetic processes, such as the TF (repressor monomer) dimerization, the TF binding to gene operator, or the mRNA (genes) synthesis reactions (Table 4) [59-61]. This observation sustains the physical meaning of model parameters, thus increasing the HSMDM model robustness.

As a final observation, by extending the detailing degree of the bioreactor dynamic model at a cellular level, the resulting structured HSMDM model not only preserves but also extends the adequacy of the unstructured model, adding the possibility to predict the cell species/fluxes dynamics over dozens of cell cycles. By offering details on the cell metabolism adaptation, the ’intrinsic’ reduction rate, and the possibility to in-silico predict the modified (cloned) cell response to various stimuli, such an HSMDM model presented in this paper is superior compared to the unstructured (apparent) bioreactor dynamic models. Eventually, the extended HSMDM models are worth the supplementary experimental and computational effort to derive/identify them.

Conclusions

Meritorious structured deterministic CCM kinetic models have been reviewed by Maria [2]. Deterministic kinetic models using continuous variables have been developed by Maria [3] for the glycolysis, and by [77,104-107] for the CCM in bacteria of industrial interest. Such models can adequately reproduce the cell response to continuous perturbations, the cell model structure and size being adapted based on the available –omics information. Even if such extended structured models are currently used only for research purposes, being difficult to identify, it is a question of time until they will be adapted for industrial / engineering purposes in the form of reduced HSMDM models.

In other words, this work presents a holistic “closed loop” approach that facilitates the control of the in vitro through the in silico development of dynamic models for living cell (biological) systems [108], by deriving deterministic modular, structured cell kinetic models (MSDKM), with continuous variables, and based on cellular metabolic reaction mechanisms.

The ever-increasing availability of experimental (qualitative and quantitative) information, at the cell metabolism level, but also about the bioreactors’ operation necessitates the advancement of a systematic methodology to organize and utilize these data. The resulting HSMDM was proved to successfully solve difficult bioengineering problems. In such HSMDM models, the cell-scale model part (including nano-level state variables) is linked to the biological reactor macro-scale state variables for improving the both model prediction quality and its validity range. The case study presented and discussed here proves this engineering aspect. An alternative compromise is to use hybrid models that combine unstructured with structured process characteristics to generate more precise predictions (see the review of Maria [74]). Hybrid models use a two-level hierarchy: the bioreactor macroscopic state variables linked with the nano-scale variables describing the cell key metabolic processes, and those of practical interest.

As proved by the case study presented in this paper, and by some additional ones mentioned in the Introduction section of this work, the use of MSDKM and HSMDM models (developed under the novel WCVV modelling framework) to simulate the dynamics of the bioreactor and, implicitly, the dynamics of the key cellular metabolic processes (CCM-based, related GRC-s, and target metabolites syntheses) occurring in the bioreactor biomass, presents multiple advantages, such as: (i) a higher degree of accuracy and the prediction detailing for the bioreactor dynamic parameters (at both macro- and nano-scale level); (ii) the prediction of the biomass metabolism adaptation over tens of cell cycles to the variation of the operating conditions in the bioreactor; (iii) prediction of the CCM key-species dynamics, by also including the metabolites of interest for the industrial biosynthesis; (iv) prediction of the CCM stationary reaction rates (i.e. metabolic fluxes) allow to in-silico design GMO of desired characteristics.

The superiority of structured HSMDM models is proved by several case studies approached by the author and briefly mentioned in the Introduction section of this paper. For details, the reader is asked to consult the above-indicated references.

However, to avoid the intensive experimental and computational efforts necessary to develop an extended structured HSMDM model, unstructured dynamic models of bioprocesses continue to be used for various engineering purposes. Even if of low precision, the unstructured/global models for the bioreactors and the conducted bioprocesses continue to be largely used in the engineering practice, for a quick design, optimization, or control of the industrial bioprocesses. For instance, for the case study mentioned in this paper, a method was presented to optimize the operating policy of a SCR bioreactor by using the reduced unstructured hybrid model. Maria [4,62,92].