A Modular Mathematical Model of Fermentation in an Industrial-Scale Bioreactor

1. Introduction

The growing challenges of climate change demand effective and sustainable solutions to meet the needs of a rising global population while minimizing environmental impact [1]. Biomanufacturing is the biotechnological production of valuable chemicals [2] which has emerged as a key approach to address these challenges. By relying on renewable alternatives, biomanufacturing helps reduce dependence on fossil-based resources [2,3,4].

A continuous culturing is widely used in biomanufacturing, as it enables the controlled growth of microorganisms that produce target products. Advances in biological systems engineering have enabled laboratory-scale production of a wide range of compounds, including pharmaceuticals and precursors for materials and fuels, by utilizing microbial chassis at least in small amounts as a proof-of-concept [4,5,6,7,8,9]. There are also examples of biomanufacturing at an industrial scale, including the production of both high-value products (e.g., therapeutics, enzymes, antibiotics) and low-volume commodities (e.g., bioethanol, lactic acid, amino acids) [10]. However, the successful transition from a laboratory scale (with cultivation volumes ranging from milliliters to tens of liters) [11,12] to an industrial scale (with cultivation volumes reaching hundreds of cubic meters) [13,14] remains challenging and is often associated with significant physical, chemical, and biological limitations [2,12,15]. These constraints often prevent promising developments from being successfully translated into large-scale production [16,17].

One of the key limiting factors during scale-up is the emergence of spatial heterogeneity within large bioreactors. In laboratory-scale bioreactors, rapid mixing (on the order of seconds) results in a well-mixed, homogeneous environment. In contrast, industrial-scale bioreactors may require several minutes to reach ideal mixing conditions that lead to the formation of zones with concentration gradients of substrate, oxygen, and metabolic products [18,19]. For example, cells may shift toward anaerobic metabolism in zones with limited oxygen availability, leading to the accumulation of acetate – a byproduct that can inhibit microbial growth [20,21]. On the other hand, in substrate-depleted zones, cultures experience nutrient starvation, which reduces overall bioprocess efficiency [21,22].

Mathematical modeling plays a key role in analysis and optimization of industrial bioprocesses, as well as in mitigation of risks during scale-up from laboratory to industrial scale. Four main types of bioreactor models are commonly distinguished: kinetic, compartment, computational fluid dynamics (CFD), and hybrid models. Kinetic models provide a detailed description of metabolic pathways based on chemical kinetics, while generally neglecting hydrodynamic properties of the reactor [23]. In contrast, compartment and CFD models simplify the metabolic representation and instead focus on predicting the formation of mixing gradients. Compartment models consider the reactor as multiple ideally mixed zones (compartments), simulate liquid exchange flows between adjacent zones, and capture concentration gradients of substrate, oxygen, and other key compounds [19,24,25,26]. CFD models, on the other hand, calculate fluid dynamics by solving the Navier–Stokes equations, treating the liquid phase as a continuous medium [27,28]. While CFD models provide high spatial resolution, they require significant computational resources, making compartment models a computationally inexpensive alternative when moderate detail is sufficient. Finally, hybrid models combine kinetic models with either compartment or CFD models, in order to predict the spatial distribution of compound concentrations under realistic hydrodynamic conditions [29].

The objective of this study was to develop a hybrid mathematical model of an industrial bioreactor that simulates the growth of E. coli while accounting for key processes in central carbon metabolism and the spatial environment heterogeneity induced by mixing. To this end, we reproduced results from two previously published models on the BioUML platform: a simplified four-compartment model of a large stirred bioreactor [30] and a detailed kinetic model for the central carbon metabolism of E. coli [31]. These models were then combined into a unified modular model.

The structure of the paper is as follows. Section 2 provides a brief overview of the BioUML platform used for the model development and numerical simulations. In Section 3, the modular structure and main equations of the model are detailed. Section 4 presents the simulation results and their interpretation. Section 5 discusses the biological relevance of the results and outlines potential directions for further model refinement. Supplementary Materials (Tables S1–S3) provide the complete description of model equations, parameter values, and initial conditions.

2. Methods

2.1. Modeling Software

Model construction and numerical simulations were carried out using BioUML (Biological Universal Modeling Language, https://sirius-web.org/bioumlweb/), an open-source platform under active development [32,33]. BioUML is a Java-based software environment designed for modeling complex biological systems and analyzing biomedical data. It provides tools for visual representation of modeled biological systems and supports a wide range of systems biology standards, including Systems Biology Markup Language (SBML) [34] and Systems Biology Graphical Notation (SBGN) [35].

Model development in BioUML begins with representation of the modeled system as a compartmentalised graph. Based on this graph, the platform automatically generates Java code that defines the corresponding system of differential equations. Then this code is compiled and executed using internal numerical solvers to simulate the system’s dynamics.

To simplify the construction of models for complex systems, BioUML supports a modular approach, which enhances visual clarity and reduces the risk of technical errors. This approach, previously used to build various mathematical models [36,37], allows decomposition of any complex system into elementary components called modules (or submodels/blocks). Each module is defined independently and can follow its own formalism and level of abstraction. Interactions between modules are defined through variable connections, which, in the context of the bioreactor model, represent processes such as fluid exchange between adjacent compartments.

2.2. Model Parameters

The modular bioreactor model integrates two previously published models, using the same parameter values as reported in the original studies where they were originally described and validated [30,31].

2.3. Numerical Simulations

The unit of model time corresponds to one hour of real time. Numerical simulations were performed using the CVODE solver [38], ported to Java and integrated into the BioUML platform. All simulations presented in this paper were conducted for a duration of up to 1 hour.

3. Modular Model

3.1. Bioreactor Compartmentalization

The bioreactor compartment model utilizes the geometry of a 90-cubic-meter stirred bioreactor with four Rushton turbines, as described in [14,39]. The total reactor volume is divided into four vertical compartments using three horizontal planes, positioned equidistantly between the turbines. Three of the lower compartments have a volume of 20 m³ each, while the top compartment has a volume of 30 m³. A schematic representation of the bioreactor compartmentalization is shown in Figure 1.

3.2. Modular Model Structure

The mathematical model of a large stirred bioreactor, which describes the growth of E. coli and accounts for key processes in its central carbon metabolism as well as concentration gradients induced by mixing, was developed based on the principles of modular modeling for complex biological systems using the BioUML platform (see Section 2.1).

The BioUML implementation of the model consists of 9 modules: 4 metabolic modules (Metabolic_i) that calculate the dynamics of E. coli metabolism in the corresponding reactor zones; 4 compartment modules (Compartment_i) that calculate the liquid exchange flows between neighboring zones; and an auxiliary module (Feed_rate) that calculates the glucose feed rate to the top zone. The modular structure of the model in BioUML is shown in Figure 2.

The model includes 166 variables, 57 differential equations, 377 parameters, and 96 assignment operations. Hereafter, it is referred to as the modular model of the bioreactor.

3.3. Main Equations of the Modular Model

This section provides a brief overview of the main equations underlying the modular model of the bioreactor (Equations (1) – (23)). A complete description of all equations, parameter values, and initial conditions is given in Supplementary Materials (Tables S1 – S3).

The metabolic modules are based on a previously published kinetic model of the central carbon metabolism of E. coli [31]. This model includes the glycolytic pathway, tricarboxylic acid cycle, pentose phosphate pathway, Entner-Doudoroff pathway, anaplerotic pathway, glyoxylate shunt, oxidative phosphorylation, phosphotransferase and non- phosphotransferase systems, and alternative metabolic pathways. Regulation of metabolic activity is also incorporated through four key transcription factors: the cAMP receptor protein, the catabolite repressor/activator, the pyruvate dehydrogenase complex repressor, and the isocitrate lyase regulator.

The dynamics of each metabolic module are governed by the following set of equations (Equations (1) – (4)):

$${\displaystyle \frac{dX}{dt}}=\mu \left(x_1,x_2,y,p\right)X-DX,$$

(1)

$${\displaystyle \frac{dy}{dt}}=D\left(y_{feed}-y\right)-g\left(x_1,x_2,y,p\right)X,$$

(2)

$${\displaystyle \frac{d{x}_1}{dt}}=f\left(x_1,x_2,y,p\right),$$

(3)

$$x_2=h\left(x_1,y,p\right),$$

(4)

where $X$ denotes the biomass concentration in the compartment; ${\mathbf{x}}_1$ is the vector of time-dependent intracellular variables (e.g., concentrations of metabolites, proteins and enzymes); ${\mathbf{x}}_2$ is the vector of auxiliary variables, including intracellular regulatory proteins (such as transcription factors and their complexes with metabolites); $\mathbf{y}$ is the vector of extracellular concentrations of glucose and acetate; the vector ${\mathbf{y}}_{feed}$ contains the concentrations of glucose and acetate in the feeding medium; $\mathbf{p}$ is the vector of the kinetic parameters of the model; and $D$ is the dilution rate. The function $\mu \left({\mathbf{x}}_1,{\mathbf{x}}_2,\mathbf{y},\mathbf{p}\right)$ defines the specific growth rate; $\mathbf{g}\left({\mathbf{x}}_1,{\mathbf{x}}_2,\mathbf{y},\mathbf{p}\right)$ corresponds to the specific uptake rates of glucose and acetate; and $\mathbf{f}\left({\mathbf{x}}_1,{\mathbf{x}}_2,\mathbf{y},\mathbf{p}\right)$ represents the system of intracellular mass balance equations.

To incorporate mixing dynamics, the system was extended with compartment modules based on the simplified compartment model of a bioreactor [31]. The total reactor volume was divided into four compartments (see Section 3.1). To simulate the mixing process, additional exchange terms were introduced into the right-hand side of the differential equations describing the extracellular concentrations of glucose, acetate, and biomass in the metabolic modules (Equations (5) – (7)):

$$S_{i,flow}={\displaystyle \frac{1}{V_{L,i}}}\left({\displaystyle \sum F_{ex,j\leftrightarrow i}\cdot S_j}-{\displaystyle \sum F_{ex,i\leftrightarrow j}\cdot S_i}\right),$$

(5)

$$A_{i,flow}={\displaystyle \frac{1}{V_{L,i}}}\left({\displaystyle \sum F_{ex,j\leftrightarrow i}\cdot A_j}-{\displaystyle \sum F_{ex,i\leftrightarrow j}\cdot A_i}\right),$$

(6)

$$X_{i,flow}={\displaystyle \frac{1}{V_{L,i}}}\left({\displaystyle \sum F_{ex,j\leftrightarrow i}\cdot X_j}-{\displaystyle \sum F_{ex,i\leftrightarrow j}\cdot X_i}\right),$$

(7)

where $S$, $A$ and $X$ correspond to the extracellular concentrations of glucose, acetate, and biomass, respectively, and $F_{ex,j\leftrightarrow j}$ represents the liquid exchange flow between neighboring compartments i and j.

The calculation of liquid exchange between compartments was based on the volumetric pumping capacity of Rushton turbines under aerated conditions. It was assumed that half of the liquid flow generated by a turbine is directed downwards and the other half upwards [40,41]. Thus, the liquid exchange flow between compartments i and j ($F_{ex,j\leftrightarrow j}$) was defined as the average of the two volumetric flows meeting at the interface between these compartments (Equation (8)):

$$F_{ex,i\leftrightarrow j}={\displaystyle \frac{{\theta }_{P,g,i}+{\theta }_{P,g,j}}{4}},$$

(8)

with the following expressions (Equations (9) – (15)) used to compute intermediate parameters:

$${\theta }_{P,g,i}=0.2\cdot N_P\cdot {\left({\displaystyle \frac{P_g}{P_O}}\right)}_i\cdot N\cdot D^3,$$

(9)

$${\left({\displaystyle \frac{P_g}{P_O}}\right)}_i=0.0312\cdot F{r}^{-0.16}\cdot R{e}^{0.064}\cdot {\left({\displaystyle \frac{Q_{g,i}}{N\cdot D^3}}\right)}^{-0.38}\cdot {\left({\displaystyle \frac{T}{D}}\right)}^{0.8},$$

(10)

$$Re={\displaystyle \frac{{\rho }_L\cdot N\cdot D^2}{\eta }},$$

(11)

$$Fr={\displaystyle \frac{N^2\cdot D}{g}},$$

(12)

$$Q_{g,i}={\displaystyle \frac{F_g}{{\rho }_{g,i}}},$$

(13)

$${\rho }_{g,i}={10}^5\cdot {\displaystyle \frac{p_{abs,i}\cdot M{W}_G}{R\cdot Temp}},$$

(14)

$$p_{abs,i}={\displaystyle \frac{{\rho }_L\cdot g\cdot h_i}{{10}^5}}+p_{head}+p_{atm},$$

(15)

where $Q_{P,g,i}$ is the volumetric pumping capacity of a single Rushton turbine (with a power number $N_P$) in compartment i under aeration; $P_g/P_O$ is the ratio of aerated to unaerated power input; $Re$ and $Fr$ are the Reynolds and Froude numbers, respectively; $D$ and $N$ are the impeller diameter and rotation speed; $T$ is the reactor diameter; ${\rho }_L$ is the liquid density; $\eta$ is the dynamic viscosity; $g$ is the gravitational acceleration; $Q_{g,i}$ and $F_g$ denote the volumetric and mass flow rates of the gas, respectively; ${\rho }_{g,i}$ is the gas density; ${MW}_G$ is the molecular weight of air; $p_{abs,i}$ is the absolute pressure; $R$ is the universal gas constant; $Temp$ is the absolute temperature; $h_i$ is the height of the liquid above the impeller; $p_{head}$ is the excess pressure in the headspace of the reactor; and $p_{atm}$ is the atmospheric pressure.

The original kinetic model of E. coli central metabolism [31] did not account for changes in oxygen concentration. Therefore, in the modular model of the bioreactor developed in this study, oxygen dynamics were incorporated by adapting the gas–liquid mass transfer equations from the compartment model of the bioreactor [30].

The gas hold-up (i.e., the volume fraction of the gas relative to the total volume) in each compartment was calculated using the following expressions (Equations (16) – (18)):

$${\left(1-\epsilon \right)}_i=0.13\cdot {\left({\displaystyle \frac{P_{s,i}}{V_i}}\right)}^{{\displaystyle \frac{1}{3}}}\cdot {\left(v_{g,i}\right)}^{{\displaystyle \frac{2}{3}}},$$

(16)

$$P_{s,i}=N_P\cdot {\left({\displaystyle \frac{P_g}{P_O}}\right)}_i\cdot {\rho }_L\cdot N^3\cdot D^5,$$

(17)

$$v_{g,i}={\displaystyle \frac{Q_{g,i}\cdot 4}{\pi \cdot T^2}},$$

(18)

where ${\epsilon }_i$ represents the volume fraction of the bulk liquid; $V_i$ is the compartment volume; $P_{g,i}$ is the power input by a single impeller under aeration; and $v_{g,i}$ is the superficial gas velocity. Due to the gas hold-up, the total liquid volume across all compartments is slightly reduced compared to the full bioreactor volume of 90 m³.

Oxygen transfer from the gas phase to the bulk liquid in each compartment was modeled under the assumption that air bubbles do not coalesce. The following expressions (Equations (19) – (21)) were used:

$$OT{R}_i={\left(k_{L}A\right)}_i\cdot \left(O_i^{*}-O_i\right),$$

(19)

$${\left(k_{L}A\right)}_i=0.002\cdot {\left({\displaystyle \frac{P_{g,i}}{V_{L,i}}}\right)}^{0.7}\cdot {v_{g,i}}^{0.2},$$

(20)

$$O_i^{*}=y_{O,i}\cdot p_{abs,i},$$

(21)

where ${OTR}_i$ is the oxygen transfer rate; $k_{La}$ is the oxygen mass transfer coefficient; $O_i$ is the extracellular oxygen concentration; $O_i^{*}$ is the oxygen saturation concentration; $V_{L,i}$ is the bulk liquid volume; and $y_{O,i}$ is the mole fraction of oxygen in the gas phase.

The following expressions (Equations (22),(23)) were incorporated into the metabolic modules to compute the extracellular oxygen concentration:

$${\displaystyle \frac{d{O}_i}{dt}}=O_{i,flow}+OT{R}_i,$$

(22)

$${\displaystyle \frac{d{y}_{O,i}}{dt}}={\displaystyle \frac{\left(ga{s}_{flow,j\to i}\cdot \left(y_{O,j}-y_{O,i}\right)-OT{R}_i\cdot V_{L,i}\right)}{n_{g,i}}},$$

(23)

where ${gas}_{flow,j\leftrightarrow i}$ is the gas flow rate and $n_g$ is the number of moles of gas in the compartment.

The resulting extracellular oxygen concentration is consumed by the cells and subsequently utilized in oxidative phosphorylation, with the corresponding reaction rate contributing to the calculation of the specific growth rate. It is worth to note that the original metabolic model considers adaptation to oxygen availability according to an assumption of the oxygen excess in the medium, which allows for the instantaneous recovery of ATP, while the rate of ATP consumption is used in the rate of biomass production (Equations (24) – (28)):

$${\displaystyle \frac{dX}{dt}}=\mu X,$$

(24)

$$\mu =k_{ATP}{v}_{ATP},$$

(25)

$$\begin{array}{c}{v}_{ATP}=O{P}_{NADH}+O{P}_{FADH2}-v_{Glk}-v_{Pfk}+v_{Gapdh}+v_{Pyk}\\ -v_{Pps}+v_{Ack}-v_{Acs}+v_{\alpha kgdh}-v_{Pck}-v_{AceK}-v_{Cya},\end{array}$$

(26)

$$O{P}_{NADH}=\left(v_{Gapdh}+v_{Pdh}+v_{\alpha kgdh}+v_{Mdh}\right)\cdot \left(P/O\right),$$

(27)

$$O{P}_{FADH2}=v_{Sdh}{\left(P/O\right)}^{*},$$

(28)

where $X$ is the biomass concentration; $\mu$ is the specific growth rate; $v_{ATP}$ is the total ATP production flux; $k_{ATP}$ is the adjustable constant; ${OP}_{NADH}$ and ${OP}_{FADH2}$ are the specific ATP production fluxes through oxidative phosphorylation; $\left(P/O\right)$ and ${\left(P/O\right)}^{*}$ indicate the phosphate/oxygen ratios for NADH and FADH2, respectively under aerobic conditions; and $v_{Glk}$, $v_{Pfk}$, $v_{Gapdh}$, $v_{Pyk}$, $v_{Pps}$, $v_{Ack}$, $v_{Acs}$, $v_{\alpha kgdh}$, $v_{Pck}$, $v_{AceK}$, $v_{Cya}$, $v_{Pdh}$, $v_{Mdh}$, $v_{Sdh}$ are reaction rates for glucokinase, phosphofructokinase, glyceraldehyde 3-phosphate dehydrogenase, pyruvate kinase I, phosphoenolpyruvate synthase, acetate kinase, acetyl coenzyme A synthetase, α-keto-D-gluconate dehydrogenase, phosphoenolpyruvate carboxykinase, isocitrate dehydrogenase phosphatase/kinase, adenylate cyclase, pyruvate dehydrogenase, malate dehydrogenase, succinate dehydrogenase, respectively.

While in the proposed model, ATP is restored depending on the accumulated NADH and FADH2 and available oxygen (Equations (29) – (30)):

$$O{P}_{NADH}={\displaystyle \frac{v_{OP,max}\cdot NADH\cdot O}{\left(K_{OP,NADH}+NADH\right)\left(K_{OP,O}+O\right)}}\cdot \left(P/O\right),$$

(29)

$$O{P}_{FADH2}={\displaystyle \frac{v_{OP,max}\cdot FADH2\cdot O}{\left(K_{OP,FADH2}+FADH2\right)\left(K_{OP,O}+O\right)}}\cdot {\left(P/O\right)}^{*},$$

(30)

where $NADH$, $FADH2$ and $O$ are concentrations of NADH, FADH2 and intracellular oxygen, respectively; $v_{OP,max}$ is the maximum oxidative phosphorylation rate; $K_{OP,NADH}$ and $K_{OP,FADH2}$ are Michaelis–Menten constants.

Nevertheless, the process of oxidative phosphorylation is presented in a simplified form, taking into account the production of NADH, FADH2 and ATP. The absence of a complete oxidative phosphorylation process can become critical, especially in conditions of changing oxygen concentrations in the compartments of the bioreactor. This constraint is a jumping-off point of further development of the metabolic part of the modular model.

4. Results

4.1. Preliminary Information

The results of the numerical simulations for the modular model of the bioreactor are presented for a total simulation time of 1 hour, which was sufficient for the system to reach a steady state. It should be noted that reaching the stationary growth phase in bacterial cultures typically requires a longer timescale (5–10 hours) [43–47]. This highlights the need for further investigation of the factors and conditions influencing the long-term bioprocess dynamics, which will be addressed in future developments of the presented modular model.

Two different bioreactor mixing regimes were used: one corresponded to low mixing intensity (also referred to as “low settings”) with parameter values $N=80$ rpm, $F_g=0.4$ kg/s, $p_{head}=0$ bar; and another corresponded to high mixing intensity (“high settings”) with $N=140$ rpm, $F_g=0.79$ kg/s, $p_{head}=1$ bar. Additional simulations were performed by varying the mixing parameters between these two boundaries in order to identify an optimal mixing regime.

In all simulations, the specific glucose feed rate $q_{S_{feed}}$ was set to 0.08 grams of substrate per gram of biomass per hour (gS/gX/h).

4.2. Simulation Results for the Boundary Mixing Regimes

Figure 3 shows the dynamics of extracellular metabolite concentrations and biomass across different compartments of the bioreactor under low and high mixing intensities.

At low mixing intensity (Figure 3, left column), oxygen depletion occurs in the upper compartments of the reactor because the upward transport of oxygen from the bottom zones is insufficient to sustain active metabolism. In contrast, under high mixing intensity (Figure 3, right column), the oxygen concentration remains sufficiently high throughout all compartments, thereby preventing acetate accumulation. Since glucose feeding occurs only in the top compartment, substrate depletion is observed in the lower zones. However, under high mixing intensity, glucose distribution becomes slightly more uniform across compartments, and, combined with high oxygen availability, this leads to an accelerated biomass growth rate.

The obtained results qualitatively agree with experimental observations reported in the published data regarding the influence of hydrodynamic conditions on substrate distribution, oxygen availability, and by-product formation in industrial-scale bioreactors [30].

E. coli can simultaneously assimilate substrates through both aerobic and anaerobic pathways. Under oxygen-limited conditions, the proportion of anaerobic metabolism varies, leading to the development of distinct metabolic regimes depending on the ratio between anaerobic and total substrate uptake rates.

Table 1 and Table 2 summarize the simulated substrate uptake rates and the corresponding metabolic regimes under different mixing conditions. The classification criteria for metabolic regimes were adopted from the original compartment model [30]: substrate starvation was defined by a total substrate uptake rate below 0.04 gS/gX/h; oxygen limitation was defined when the anaerobic substrate uptake rate exceeded 0.00001 gS/gX/h; in all other cases, substrate limitation was assumed.

Three metabolic regimes were identified from the simulations: oxygen limitation, substrate limitation, and substrate starvation. Under oxygen limitation, the oxygen transfer rate into the liquid phase is insufficient to fully oxidize the available substrate, forcing cells to partially shift to anaerobic metabolism. Substrate starvation occurs when the substrate uptake rate is too low to sustain cell growth, while metabolism remains fully aerobic. In the substrate limitation regime, growth is constrained by the availability of substrate, and all consumed substrate is completely oxidized without switching to anaerobic pathways.

It can be seen that adjusting the mixing parameters primarily influences whether glucose is consumed aerobically or anaerobically in the upper part of the reactor, while the lower compartments still experience significant substrate depletion. Interestingly, increasing the mixing intensity did not reduce the overall extent of substrate starvation across the reactor. This observation is consistent with the well-known fact that Rushton turbines promote efficient radial mixing but do not contribute to better axial (vertical) mixing [19,42]. As a result, glucose remains unevenly distributed along the height of the reactor, and enhanced mixing primarily improves oxygen transfer rather than homogenizing substrate concentrations.

4.3. Optimal Mixing Regime Identification

In addition to the two boundary regimes (with low and high mixing intensity), a series of simulations was conducted by varying three mixing-related parameters within the range defined by their respective boundary values using linear interpolation:

$$par=pa{r}_{ls}+\alpha \cdot \left(pa{r}_{hs}-pa{r}_{ls}\right),$$

(31)

where $par\in \left\{N,F_g,p_{head}\right\}$ is the new value of the interpolated parameter; ${par}_{ls}$ and ${par}_{hs}$ are its values under low and high mixing settings, respectively; and $\alpha$ is the interpolation coefficient. A value of $\alpha =0$ corresponds to the low settings, while $\alpha =1$ corresponds to the high settings.

The aim of this analysis was to determine the minimum value of $\alpha$ at which the steady-state acetate concentration in all compartments stabilizes at zero. This state was interpreted as an optimal mixing regime, in which anaerobic metabolism is fully suppressed at a lower mixing intensity than in the high settings scenario. This regime avoids excessive energy consumption stirring while maintaining optimal conditions for biomass growth.

For each value of $\alpha \in \left[\mathrm{0,1}\right]$ with a step size of 0.01, a numerical simulation was performed, and the maximum acetate concentration across all compartments was recorded at the final simulation time point. Additionally, it was verified that the system had reached steady state. The resulting dependency is shown in Figure 4. The dynamics of extracellular metabolite and biomass concentrations for the identified optimal mixing regime are shown in Figure 5.

The acetate concentration was found to reach zero for the first time at $\alpha \approx 0.58$, which corresponds to a mixing regime with the following parameter values: $N=114.8$ rpm, $F_g=0.6262$ kg/s, $p_{head}=0.58$ bar.

Thus, the simulation results identified a minimum combination of mixing-related parameters at which acetate accumulation was suppressed in all compartments. This transition to fully aerobic substrate oxidation represents an optimal operating regime that sustains biomass growth rates comparable to those under high mixing conditions, while requiring lower energy demands for mixing.

5. Discussion

Fermentation processes in industrial-scale bioreactors are the focus of extensive experimental and computational research [18,48,49,50]. In this study, we developed a modular mathematical model of an industrial-scale stirred bioreactor capable of simulating the dynamics of metabolic processes in E. coli within different bioreactor zones, while explicitly accounting for concentration gradients induced by mixing. The main advantage of the model lies in its integration of a detailed kinetic description of the central carbon metabolism with a simplified hydrodynamic representation of the cultivation environment. This hybrid modeling approach provides an improved capability to capture how mixing conditions and spatial concentration gradients affect the metabolic responses of cell cultures.

The simulation results demonstrated that mixing intensity has a substantial impact on the overall dynamics within the reactor. Under low mixing intensity, oxygen limitation was observed in the upper reactor compartments, forcing cells to partially switch to anaerobic metabolism and leading to acetate accumulation. In contrast, under high mixing intensity, both substrate and oxygen are more evenly distributed, preventing acetate formation and enhancing biomass growth. Additionally, computational analysis allowed identification of an optimal mixing regime, characterized by zero steady-state acetate concentration at more than 40% lower mixing intensity compared to the high-mixing scenario. This result highlights the model’s potential in identifying energy-efficient bioreactor operating conditions, which simultaneously improve oxygen availability for aerobic metabolism and reduce risks of mechanical stress on cells associated with excessively high mixing rates [51,52].

However, the current implementation of the modular model relies on a simplified four-compartment representation of hydrodynamics, which does not account for potential intra-compartmental gradients of substrate and oxygen concentrations. This simplification leads to artificially steep transitions between metabolic regimes. In realistic systems, it is reasonable to expect an intermediate region between zones of oxygen limitation and substrate starvation – where substrate is still available but no longer sufficient to meet the full metabolic demands of the cells. The existence of such transitional regions has been computationally confirmed using CFD-based models [30].

Therefore, future improvements of the presented model should focus on incorporating more refined spatial concentration distributions, either by increasing the complexity of compartmentalization or by integrating CFD approaches.

6. Conclusions

In this study, we developed a modular mathematical model of an industrial-scale stirred bioreactor, integrating key aspects of E. coli metabolism with hydrodynamic properties of the cultivation environment. The model combines a simplified four-compartment model of a large stirred bioreactor with a detailed kinetic model for the central carbon metabolism of E. coli.

The compartment model divides the reactor volume into distinct zones and calculates local concentrations of glucose, oxygen, acetate, and biomass, as well as liquid exchange flows between adjacent compartments. The kinetic model incorporates detailed metabolic pathways, including glycolysis, the tricarboxylic acid cycle, oxidative phosphorylation, and accounts for shifts between aerobic and anaerobic metabolism based on oxygen availability. Compared to existing hybrid models that combine metabolic kinetics with hydrodynamic effects, the presented model provides the most detailed kinetic component, enabling more accurate simulation of metabolic pathways.

Numerical simulation results demonstrated the substantial impact of mixing parameters on spatial distributions of compounds and corresponding metabolic regimes within the bioreactor. Additionally, computational analysis identified an optimal mixing regime that eliminates anaerobic metabolism while reducing energy consumption for stirring. The developed model provides a computational tool for simulating cultivation processes in industrial bioreactors and analyzing metabolic dynamics under various hydrodynamic conditions.