*Article* **Predicting By-Product Gradients of Baker's Yeast Production at Industrial Scale: A Practical Simulation Approach**

#### **Christopher Sarkizi Shams Hajian 1, Cees Haringa 2, Henk Noorman 2,3 and Ralf Takors 1,\***


Received: 21 October 2020; Accepted: 25 November 2020; Published: 27 November 2020 -

**Abstract:** Scaling up bioprocesses is one of the most crucial steps in the commercialization of bioproducts. While it is known that concentration and shear rate gradients occur at larger scales, it is often too risky, if feasible at all, to conduct validation experiments at such scales. Using computational fluid dynamics equipped with mechanistic biochemical engineering knowledge of the process, it is possible to simulate such gradients. In this work, concentration profiles for the by-products of baker's yeast production are investigated. By applying a mechanistic black-box model, concentration heterogeneities for oxygen, glucose, ethanol, and carbon dioxide are evaluated. The results suggest that, although at low concentrations, ethanol is consumed in more than 90% of the tank volume, which prevents cell starvation, even when glucose is virtually depleted. Moreover, long exposure to high dissolved carbon dioxide levels is predicted. Two biomass concentrations, i.e., 10 and 25 g/L, are considered where, in the former, ethanol production is solely because of overflow metabolism while, in the latter, 10% of the ethanol formation is due to dissolved oxygen limitation. This method facilitates the prediction of the living conditions of the microorganism and its utilization to address the limitations via change of strain or bioreactor design or operation conditions. The outcome can also be of value to design a representative scale-down reactor to facilitate strain studies.

**Keywords:** scale-up; scale-down; computational fluid dynamics; *Saccharomyces cerevisiae*; mechanistic kinetic model; bioreactor; concentration gradients; digital twin; bioprocess engineering

#### **1. Introduction**

Bioprocesses are applied for the production of a vast spectrum of commodities, from food and pharmaceuticals to bioplastic and biofuel. Although different from their chemical counterparts, transferring promising lab approaches to industrial applications is a major challenge too. The problem lies within the different scales of lab, pilot, and industrial bioreactors. Whereas, ideally, mixed homogenous conditions are easily realized at lab scale, economic and physical constraints prevent the establishment of such ideal conditions in industrial tanks. As a result, gradients of limiting substrate concentrations, by-products, pH, temperature, and shear rates are formed inevitably. Circulating microorganisms in stirred and gassed large-scale tanks respond to the permanently changing microenvironmental conditions, finally causing uncertainty of process performance, possibly deteriorating key TRY criteria (titer, rate, yield) [1–6].

Different approaches to addressing such issues have been studied by bioreactor design experts [7–12]. The investigation of cellular interaction with substrate gradients has been the basis of numerous studies. Often, deteriorating TRY values have been reported [13,14] but there are also occasional observations of improved performance values [15]. In particular, *Saccharomyces cerevisiae* was found to respond to fluctuating conditions with improved viability [16], adapted cAMP-mediated metabolism [17–19], and global regulation [20,21]. However, most of these investigations were performed at laboratory scale, mimicking industrial-scale conditions. Real industrial-scale data, important for validation, are rare. Accordingly, researchers have been employing computational fluid dynamics (CFD) combined with metabolic models with different resolutions to shed light on gradients in the bioreactor that take place at the interface of various physical and biological phenomena [6,22–25]. It is worth mentioning that shear gradients may have a significant effect on shear sensitive hosts. Simulating shear fields and calculating the frequencies of cellular exposure may be a highly valuable tool for future applications to investigate this particular large-scale impact [26].

*Saccharomyces cerevisiae* strains are applied for a wide range of processes, from the food industry [27–31] and bulk chemical production [32–37] to the pharmaceutical industry [15]. These products are manufactured in large bioreactor scales, where gradients cannot be avoided. As one of the workhorses of industrial biotechnology, this yeast is well investigated [38,39], and key traits of the "Crabtree positive" strains [40–45] are thoroughly studied. One feature is the overflow production of ethanol under aerobic conditions, mirroring how *S. cerevisiae* consumes more glucose than it can metabolize. [46–48]. Under anaerobic conditions *S. cerevisiae* is known to consume ethanol [49]. Apparently, such traits may gain importance under varying microenvironmental conditions affecting TRY values and effecting population differentiation [50]. Interesting short- and long-term Crabtree effects have been observed [44] that create population responses at different timescales. Furthermore, stress exposure may cause growth phenotypes different from the well-reported Monod-type kinetics using substrate supply as the growth limiting impact. To this end, obtaining information regarding the surroundings of the cell will help to identify the triggers that could set off the stress responses [51,52], finally yielding further improved prediction of large-scale performance of the yeast.

The mathematical framework requires the joint use of CFD with microbial kinetics. Whereas the first is applied to predict hydrodynamics and mass transfer in large tanks, the latter describes the cellular phenotype which is dominated by metabolic models [53–57]. The effects interact and are both needed to predict gradients in large-scale reactors. However, recently, efforts have been made to integrate multiple levels of cellular regulation linking metabolism with enzyme activities [9] and gene expression [6]. The integration of these hierarchical control levels expands the timescales of cellular response severely [6,21,58–61], which causes extra computational burden, challenging conventional simulation capacities. With this in mind, such biokinetic models should be linked to CFD that describe the targeted phenotype with the least computational effort.

Nowadays, CFD is one of the must-have tools for process development and troubleshooting [5,62–65]. It provides the possibility to incorporate the main aspects of the process and yield further insights into the conditions inside the bioreactor. Because of the dynamic environment faced by the cells in large-scale bioreactors, different responses of cells will give rise to a heterogeneous population, which will result in heterogeneities in substrate/by-product gradients and/or in the cell population [46,48,66,67]. Setting up a simulation by itself needs to be purpose-driven for engineering applications. The computational resources and strategies should be allocated in a way that contributes the most to the goal of the project at hand.

Currently, Reynolds-averaged Navier–Stokes (RANS) methods are the most common way to model hydrodynamics including turbulence. However, more delicate methods are surfacing [68–70], with drastic changes in software and hardware requirements. Most of the bioproducts are produced in aerated fermenters. Moreover, multiphase problems are addressed with numerous methods. On the one hand, Euler–Euler approaches are used for considering the impact of the ever-changing environment on cells [71]. On the other hand, one may implement the Euler–Lagrange approach [72], for instance, if individual microorganisms are incorporated, typically as massless particles.

In this work, CFD is coupled to a mechanistic biokinetic model to predict concentration gradients in a large-scale bioreactor. The process of interest is the production of baker's yeast, *Saccharomyces cerevisiae*. A detailed assessment in terms of concentration gradients is undertaken. Aspects of biokinetics, mass transfer, and thermodynamics of the bioreactor are well characterized, and fundamentals are well established. Together with experimental validation, they prove to be an effective tool to shed light on the gradients within a bioreactor at industrial scale. By efficiently combining CFD with simple yet practical models, the assumption of non-limiting dissolved *O*<sup>2</sup> (*dO*2) is evaluated. In addition to this, dissolved *CO*<sup>2</sup> (*dCO*2) inhibition is known to occur at large scales [4,73] and is addressed within this work. Large-scale processes are suspected to lack real starvation in the case of overflow metabolism products that are formed near the feeding point and which are then re-consumed at the other end of the tank [4].

#### **2. Materials and Methods**

All pre-processing, solution, and post-processing were carried out with Ansys® Academic Research Workbench (2019 R1, Ansys, Canonsburg, PA, USA).

#### *2.1. Geometry and Mesh*

The reactor dimensions were based on the schematics available in the literature [25,74], as illustrated in Figure 1. For CFD simulation, a mesh with ~2.5 million hexahedron cells was generated and evaluated for turbulent kinetic energy (k) and dissipation of turbulent kinetic energy (ε) and velocity profile at the impeller heights, as discussed in the work of Haringa et al. [25].

**Figure 1.** Stavanger bioreactor configuration and dimensions.

#### *2.2. CFD Setup*

A Eulerian model with two phases was used. Turbulence was modeled with realizable k-ε using mixture formulation for calculating turbulence of the gas phase [75]. To investigate the concentration of chemical species, namely glucose, *O*2, *CO*2, ethanol, every phase was defined as a mixture. The physical properties of the liquid phase were approximated as water with volume weighted mixing law for density calculation. The gas phase was assumed as ideal gas considering an average bubble diameter 0.009 m as calculated by Haringa et al. [25]. With this assumption of single bubble size, the swarm coefficient needed to be set to −1.2 in the grace drag model [76] to reproduce the gas flow regime at the bottom impeller. Surface tension was set to that of the air–water system and was 0.072 N/m. Other interphase forces were assumed not to impose significant impacts according to Scargiali et al. [77].

Boundary conditions for walls were set with no-slip conditions for liquid phase and free slip for gas phase other than the impellers, where the no-slip conditions also apply for the gas phase to improve the reproduction of the vortices behind the impellers, as suggested by Haringa et al. [25]. Gas enters the bioreactor via a sparger through a mass flow inlet with 0.231 kg/s and leaves the system on top via a degassing boundary condition. Operating pressure was set close to the boundary to 130,710 Pa [4,25,74] and operating density was set to 0, as suggested by the Fluent manual [78]. The rotation of the impeller was modeled using multiple reference frames (MRF) at 2.22 1/s, as mentioned in previous works on this bioreactor [25]. The operating conditions are summarized in Table 1.


**Table 1.** Operating conditions for the fermentation process.

After setting up the phenomena describing the system, the solution methods and strategies were included. Phase coupled SIMPLE were chosen for pressure-velocity coupling and temporal (transient formulation) and spatial discretization scheme were set to first-order upwind for the first few hundred iterations to achieve solution stability and then were set to QUICK for velocity, turbulence, and volume fraction, and temporal discretization was set to bounded second-order implicit. The residuals were set to 10−<sup>6</sup> and a time step of 0.001 s with 50 iterations was chosen.

Simulations were qualified as "accomplished" once flow velocities and turbulent kinetic energies converged to pseudo-steady states with ±5% variation. Further validation was achieved by comparing mixing time (τ95) and integral mass transfer coefficient (*kla*) with published values. Flow fields served as basis for implementing mass transfer and biokinetics in subsequent steps.

#### *2.3. Biokinetics*

Sonnleitner and Käppeli [39] introduced a black-box model to describe the substrate uptake, growth, and by-product formation of *S.cerevisiae*. Glucose is considered as substrate whereas ethanol

may serve as substrate or product depending on metabolic and environmental conditions. The model assumes the respiratory capacity of the yeast as key metabolic bottleneck. If glucose uptake exceeds respiratory limits, remaining electrons are channeled via reductive pathways, leading to the secretion of ethanol. Notably, the model also allows ethanol uptake under aerobic conditions. Under anaerobic conditions, ethanol is considered as dominating product. Details are as follows:


Uptake rates are assumed to follow Monod kinetics (1)–(3).

$$q\_{\rm s} = q\_{\rm s,max} \frac{\mathbb{C}\_{\rm s}}{K\_{\rm s} + \mathbb{C}\_{\rm s}} \tag{1}$$

$$q\_o = q\_{o\,\mu\text{max}} \frac{\mathbb{C}\_o}{\mathbb{K}\_o + \mathbb{C}\_o} \tag{2}$$

$$q\_{\varepsilon} = q\_{c, \max} \frac{\mathbb{C}\_{\varepsilon}}{K\_{\varepsilon} + \mathbb{C}\_{\varepsilon}} \tag{3}$$

To distinguish between these cases, the catabolic capacity to metabolize glucose aerobically serves as the threshold. In essence, if biomass specific glucose uptake *qs* exceeds the related oxygen demand for oxidation, *Ys o* ·*qo* meaning (4):

*qs* > *Ys o* ·*qo* (4)

Aerobic ethanol formation starts. Acetaldehyde, upstream of ethanol in the fermentation pathway, serves as electron acceptor under such conditions. Accordingly, "anaerobic" carbon flux occurs and equals the remainder of the total substrate uptake (5) and (6), which will be metabolized to ethanol.

$$q\_{\text{safe}} = \min \left( \mathbf{Y}\_{\frac{\theta}{o}} \times q\_{\theta'}, q\_{\mathbf{s}} \right) \tag{5}$$

$$
\eta\_{\text{sav}} = \eta\_s - \eta\_{\text{sac}} \tag{6}
$$

To shed light on ethanol dynamics, its consumption under aerobic conditions is also considered, prioritizing glucose [79]. Given that (7) holds true, ethanol uptake rates *qeae* are calculated as shown in (8). In essence, the min modulator compares whether oxygen demands for ethanol oxidation after glucose consumption (*qo*−*Yo s* ·*qsae*) *Yo e* exceed the Monod-type ethanol uptake kinetics *qe*.

$$
\eta\_s \prec \mathcal{Y}\_{\frac{4}{9}} q\_o \tag{7}
$$

$$q\_{\rm enc} = \min\left(\frac{(q\_o - \mathcal{Y}\_{\frac{\rm g}{s}} q\_{\rm scc})}{\mathcal{Y}\_{\frac{\rm g}{s}}}, \ q\_c\right) \tag{8}$$

A graphical illustration of respiratory bottleneck is shown in Figure 2. The concentrations are calculated and, using Equations (1)–(3), the uptake rates are calculated, upon which the rest of the model is based. Yields for by-products are stoichiometrically approximated for every reaction.

To reveal the stoichiometry between substrates and products, elemental balances are applied to the process reaction (9), which then results in three different scenarios (10)–(12).

Process reaction

$$\text{sCH}\_2\text{O} + \text{oO}\_2 + n\text{NH}\_3 \rightarrow \text{CH}\_{1.79}\text{O} \\ \text{o.56} \\ \text{N}\_{0.15} + \text{cCH}\_3\text{O} \\ \text{o.5} + \text{cCO}\_2 + \text{wH}\_2\text{O} \tag{9}$$

Process reaction for aerobic growth on glucose (*rsae*)

$$\rm{sCH}\_2O + \rm{oO}\_2 + nNH\_3 \rightarrow CH\_{1.79}O\_{0.56}N\_{0.15} + \rm{cCO}\_2 + wH\_2O \tag{10}$$

Process reaction for anaerobic growth on glucose (*rsan*)

$$\text{sCH}\_2\text{O} + n\text{NH}\_3 \rightarrow \text{CH}\_{1.79}\text{O}\_{0.56}\text{N}\_{0.15} + \text{cCH}\_3\text{O}\_{0.5} + \text{cCO}\_2 + w\text{H}\_2\text{O} \tag{11}$$

Process reaction for aerobic growth on ethanol (*reae*)

$$\text{cCH}\_3\text{O}\_{0.5} + \text{oO}\_2 + \text{nNH}\_3 \rightarrow \text{CH}\_{1.79}\text{O}\_{0.56}\text{N}\_{0.15} + \text{cCO}\_2 + \text{wH}\_2\text{O} \tag{12}$$

**Figure 2.** Graphical representation of the bottleneck concept and the interplay of substrate availability in the kinetic model for ethanol consumption under aerobic conditions and overflow metabolism based on [39].

The conservation of mass results in (13). Since all underlying phenomena do not disturb mass or elemental conservation, it is safe to conclude that there is no net conversion of elements. For this purpose, the matrices for elemental composition "*E*" (14), stoichiometry coefficients "*S.C.*" (15), and reaction rates "*r*" (16) are set up.

$$\frac{d(E \cdot C \cdot V)}{dt} = E \cdot r \cdot V \text{ (Volume of the element)} + E \cdot MTR \text{ (Mass Transfer Rate)} \tag{13}$$

$$\begin{array}{c|cccccc} & - & - & - & - & - & - & - \\ & - & 1 & 0 & 0 & 1 & - & 1 & - \\ & & - & 2 & 0 & 3 & 1.79 & 3 & 0 & 2 \\ & & & 0 & 1 & 2 & 0 & 0.56 & 0.5 & 2 & 1 \\ & & & & 0 & 0 & 1 & 0.15 & 0 & 0 & 0 \\ \end{array}$$

$$r = \begin{bmatrix} r\_{S\_i} \\ r\_{O\_i} \\ r\_{R\_i} \\ r\_{T\_i} \\ r\_{C\_i} \\ r\_{C\_i} \\ r\_{T\_i} \end{bmatrix} \tag{15}$$

$$S.C.C\_{ox} \quad S.C.C\_{ox} \quad S.C\_{ox}$$

$$\begin{array}{ccccc} \cdot & S.C\_{ox} & S.C\_{ox} & S.C\_{ox} \\ & -s\_{ox} & -s\_{ox} & 0 \\ & & \\ \alpha & -o\_{ox} & 0 & -o\_{ox} \\ & n & -n\_{ox} & -n\_{ox} & -n\_{ox} \\ & \text{x} & \text{x}\_{ox} & \text{x}\_{ox} & \text{x}\_{ox} \\ & \text{c} & 0 & \text{c}\_{ox} & -\text{c}\_{ox} \\ & \text{c} & \text{c}\_{ox} & \text{c}\_{ox} & \text{c}\_{ox} \\ & \text{w} & \text{w}\_{ox} & \text{c}\_{ox} & \text{c}\_{ox} \\ & \text{w} & \text{w}\_{ox} & W\_{\alpha\alpha} & W\_{\alpha\alpha} \end{array} \tag{16}$$

Elemental conservation results in the following system of Equation (17):

$$E \cdot r = 0\tag{17}$$

By solving this system based on growth rate and carbon source uptake rate, other rates are calculated for each reaction in the models (18)–(21).

$$\text{Aeroicic growth on glucose}: \left\{ \begin{array}{l} r\_{\text{\ $}\_{\text{sur}}} = \theta r\_{\text{\$ }\_{\text{sur}}} - r\_{\text{\ $}\_{\text{sur}}}\\ r\_{\text{\$ }\_{\text{sur}}} = \theta r\_{\text{\ $}\_{\text{sur}}} - 1.05 r\_{\text{\$ }\_{\text{sur}}} \end{array} \tag{18}$$

$$\text{Anaerodic growth on glucose}: \left\{ \begin{array}{l} r\_{\varepsilon\_{\text{sun}}} = 2r\_{\varepsilon\_{\text{sun}}} - 0.3r\_{\chi\_{\text{sun}}}\\ r\_{\varepsilon\_{\text{sun}}} = r\_{\varepsilon\_{\text{sun}}} - 0.7r\_{\chi\_{\text{sun}}} \end{array} \right. \tag{19}$$

$$\text{Aeroicic growth on ethanol}: \begin{cases} \quad r\_{\varepsilon\_{\text{arc}}} = 2r\_{\varepsilon\_{\text{arc}}} - r\_{\chi\_{\text{arc}}}\\\ r\_{\varepsilon\_{\text{box}}} = 3r\_{\varepsilon\_{\text{arc}}} - 1.05r\_{\chi\_{\text{box}}} \end{cases} \tag{20}$$

The carbon source consumption rate is calculated using Equations (5)–(8) to solve the remaining equation for a known biomass concentration. The growth rate is calculated using the yield coefficient taken from the literature [80].

#### *2.4. Mass Transfer*

The mass transfer coefficient between the two phases was considered for *O*2, *CO*2, and ethanol (21). A common assumption is to estimate the mass transfer close to the equilibrium state. The gas phase is considered well mixed (zero resistance); hence, using the film theory, one can assume the mass transfer resistance to be on the liquid side of the interface. Moreover, for dilute gases in liquid, Henry's law (22) is implemented to calculate the equilibrium concentration "*C*\* " on the interface [81,82].

$$MRR = \; k\mu \; \Delta \mathcal{C} \tag{21}$$

$$C\_i^\* = \frac{P\_i}{H} \tag{22}$$

*H* is Henry's constant for the gas component at fermentation temperature (30 ◦C). *kl* is modeled by the surface renewal approach [83], *a* is the bubble surface, and both are known to be dependent on flow characteristics [84]. Mass transfer driving force (Δ*C*) differs for *O*<sup>2</sup> and *CO*<sup>2</sup> since the direction of the transport is different, which leads to (23) and (24).

$$
\Delta \mathbf{C}\_{\mathbf{c}} = \left( \mathbf{C}\_{\mathbf{c}\_{liq}} - \mathbf{C}\_{\mathbf{c}}^{\*} \right) \tag{23}
$$

$$
\Delta \mathbb{C}\_o = \left( \mathbb{C}\_o^\* - \mathbb{C}\_{o|\_{\mathbb{H}}} \right) \tag{24}
$$

For ethanol stripping, an approach by Löser et al. [85] is implemented in which ethanol transfer to gas phase is investigated. In this way, a partition coefficient (25) linking ethanol concentration in both phases to each other is available, which allows calculation of the mass transfer driving force for ethanol.

$$K\_{\dot{\mathfrak{k}}} = \frac{\mathbb{C}\_{\mathfrak{c}\_{\text{liq}}}}{\mathbb{C}\_{\mathfrak{c}\_{\text{gas}}}} \tag{25}$$

#### **3. Results**

#### *3.1. Flow Field Validation*

For the validation of the represented flow field, several criteria were evaluated. τ<sup>95</sup> was reproduced with a virtual pulse of glucose above the top impeller and by reading its concentration at the probe location, as disclosed in [74], at 0.97 m distance from the bottom. The estimated value in the current work is 186 s, which is in agreement with the results of previous investigations [4,25,74].

The simulated gas holdup of 19% slightly overpredicts experimental measurements (17.1%) [86] and previous numerical studies (17.6%) [25] but still falls within an acceptable range. In addition to the average holdup, gas distribution of the gas phase plays a crucial role in *kla* calculations (Figure 3).

**Figure 3.** Gas holdup distribution can reproduce the loading regime at the bottom impeller.

Under said conditions, *kla* of 190 h−<sup>1</sup> is estimated, which agrees well with the experimental measurements [86].

#### *3.2. Scenario I: Experimental Fermentation*

For this scenario, conditions are considered as explained by previous investigations [25,86]. Accordingly, model suitability could be checked. *dO*<sup>2</sup> concentrations are expected to show high values at the bottom of the tank for mainly three reasons: first, higher hydrostatic pressure increases the solubility of *dO*2. Second, lower metabolic activity of cells should occur due to substrate scarcity, and third, higher fraction of *O*<sup>2</sup> in gas phase close to the bottom should be observed. In contrast, opposite trends are found close to the feeding point, where the lowest *dO*<sup>2</sup> of ~3.7 <sup>×</sup> 10−<sup>5</sup> M (~1.2 ppm) is estimated (Figure 4a). This is an order of magnitude larger than the critical value of ~4.6 <sup>×</sup> <sup>10</sup>−<sup>6</sup> M [87]. Furthermore, the related volume is only a negligible fraction of the entire stirred tank reactor, which gives rise to the fair assumption of "no oxygen limitation" in the bioreactor for given conditions [4,25,86].

**Figure 4.** Concentration profiles for scenario I with 10 g/L biomass for (**a**) *dO*2, (**b**) *dCO*2, (**c**) ethanol in liquid phase, (**d**) stripped ethanol mole fraction in the gas phase.

*dCO*<sup>2</sup> concentrations are represented in Figure 4b. Interestingly, only minor *dCO*<sup>2</sup> gradients occur, varying no more than ±5% from average. Notably, CO2/HCO3 − creates a buffering system that consists of around 99% *CO*<sup>2</sup> at the operational pH 5 [73,88]. Accordingly, the simplifying assumption was made that inorganic carbon only encompasses *dCO*2. The current study snapshots a pseudo-steady state of late phase yeast fermentation [82,89]. Consequently, the reasonable postulation was made that the liquid phase is saturated with respect to *dCO*2. In the gas phase, *CO*<sup>2</sup> is estimated to reach mole fractions between 2.1 and 2.9%.

Ethanol gradients are more pronounced than those of *dCO*<sup>2</sup> but less so than *dO*2. The highest values are observed proximate to the feeding point, reflecting the highest cellular product formation and reduced stripping. The lowest titer is found at the bottom of the tank (3.25 <sup>×</sup> <sup>10</sup>−<sup>5</sup> M, Figure 4c,d).

#### *3.3. Scenario II: Protein Production*

To place the model into a more industrially relevant context, biomass was increased to 25 g/L to imitate protein production [15]. For simplification, putative impacts of increased biomass concentration on the viscosity were neglected [90]. The scenario shows that a significant volume is exposed to oxygen limitation (approximately 0.37 m3) above the top impeller, considering that 10% of saturation *dO*2. *dO*<sup>2</sup> levels below 4.6 <sup>×</sup> <sup>10</sup>−<sup>6</sup> M were observed in a volume of 0.04 m3, which is below the *dO*2,*crit* according to the available literature [87] (Figure 5a). Notably, elevated viscosity values would have even deteriorated the oxygen supply. Increasing biomass concentrations also increased microbial substrate consumption and product formation rates. As aeration and the energy input of the bioreactor remained equal, gradients for substrates and by-products became more pronounced.

**Figure 5.** Concentration profiles with 25 g/L biomass concentration for (**a**) *dO*<sup>2</sup> and (**b**) ethanol.

In turn, this affected the ethanol gradient twofold. First, the drop in *dO*<sup>2</sup> resulted in higher production of ethanol around the feeding point. Second, higher biomass concentration in the tank increased the volumetric ethanol consumption, causing lower ethanol concentrations at the bottom of the tank, as observed in Figure 5b.

#### **4. Discussion**

#### *4.1. Glucose Gradient*

Figure A1 (Appendix A) shows the anticipated heterogeneous glucose distribution, disclosing a hotspot of high glucose concentrations close to the inlet and low values at the bottom of the bioreactor. As expected, the resulting gradients are more pronounced for the "25 g/L biomass" case than for 10 g/L. Interestingly, studies [4,73] provided experimental values sampled from the top, middle, and bottom regions of the bioreactor (Table A1). Notably, sampling was performed at the wall of the bioreactors, only giving a very restricted resolution of local conditions. On the contrary, simulated values of scenario I and II indicate average concentrations of total reactor slices calculated at the same height. Consequently, the comparison of simulated predictions with experimental values is intrinsically biased. Nevertheless, the comparison shows that scenario II comes closest to the measurements. At the top, high glucose levels were equally predicted by simulation and measurements. Notably, each value indicates saturated glucose uptake. The strongest deviations are found for the bottom region, where simulations overestimate the glucose consumption. Consequently, model refinements should be considered in the next generation of metabolic models by implementing the co-consumption of intracellular buffers (such as trehalose) as an additional, not yet considered, carbon source in nutrient-limiting regions.

#### *4.2. Ethanol Gradient*

To the best of our knowledge, this study represents the first example of considering ethanol formation and re-consumption in a CFD-linked large-scale bioreactor simulation. Figure 4c indicates the well-distributed presence of ethanol in the entire reactor, giving rise to the assumption that ethanol-based growth should be possible in large parts of the bioreactor. Based on the results, growth on ethanol is expected to take place in 97% bioreactor. The conclusion is in accordance with the work of Noorman [4], who anticipated that no real "starvation" zone might exist because of the occurrence of ethanol. The finding does have implications for the design of proper scale-down approaches [63,65,91] as suitable settings ideally should consider the co-substrate ethanol too. For scenario II, the average ethanol concentration was approximately 26% lower (3.17 <sup>×</sup> <sup>10</sup>−<sup>5</sup> M) compared to scenario I (5 <sup>×</sup> <sup>10</sup>−<sup>5</sup> M). Nevertheless, more than 90% of the tank may offer sufficient ethanol uptake within seconds according to a radiocarbon study [92]. Such levels might be enough to prevent an actual starvation scenario.

#### *4.3. Oxygen Gradient*

A conventional approach for estimating the occurrence of gradients is the comparison of critical timescales τ for substrate supply τ*ssupply* versus substrate consumption τ*scons*. Whereas the first may be approximated by the mixing time τ*mix* or circulation time τ*circulation*, the latter resembles the quotient of average substrate concentration divided by volumetric substrate consumption rates (26). Regarding *dO*2, scenarios 1 and 2 anticipate the occurrence of gradients because τ*ocons*,1 and τ*ocons*,2 showing ~28 s and ~30 s are smaller than τ*circulation* ≈ 47 s (τ*mix* = 186 s). Indeed, the expectation is met by the CFD simulations.

$$\tau\_{o\_{\rm{cws}s}} = \frac{\mathbb{C}\_{\mathcal{o}}}{\left(\mathbb{Y}\_{\frac{\mathfrak{e}}{s}} \cdot q\_{\rm{sac}} + \mathbb{Y}\_{\frac{\mathfrak{e}}{s}} \cdot q\_{\rm{cac}}\right) \times \mathbb{C}\_{\mathfrak{x}}} \tag{26}$$

Assuming average values, this approach theoretically indicates that assuming the non-limiting role of *dO*<sup>2</sup> for scenario I is a reasonable approximation for the whole tank, other than a small region around the feeding point, where the *dO*<sup>2</sup> concentration is slightly above 3.67 <sup>×</sup> <sup>10</sup>−<sup>5</sup> M. The threshold for aerobic growth *dO*2,*crit* for *S. cerevisiae* is given as 4.6 <sup>×</sup> <sup>10</sup>−<sup>6</sup> M [87]. This allows us to make a distinction between the ethanol production caused by overflow metabolism or *dO*<sup>2</sup> limitation. Accordingly, no *dO*<sup>2</sup> limitation is observed in scenario I as mentioned; hence, all ethanol production in this case is attributed to overflow metabolism, which occurs in 1.63 m3 of the fermenter. However, this is not the

case for scenario II, where 10% of the ethanol production takes place in regions with *dO*<sup>2</sup> below the critical value. The total volume associated with ethanol production in scenario II is 0.3 m3.

#### *4.4. Carbon Dioxide Gradient*

Although the *dCO*<sup>2</sup> gradient is practically absent when compared to those for *dO*2, ethanol, and glucose (Appendix A), the key observation is the generally high level inside the bioreactor due to overpressure applied in the headspace plus the hydrostatic pressure from the liquid column. This should be accounted for in experimental scale-down. It should be noticed that by using a black-box model, some inherent flaws of such models affect the results. While such models could prove to be insightful for a specific case, the assumptions upon which the model is founded limit the generalization. In this case, the process reaction (6) considers *dCO*<sup>2</sup> only as a product of a single reaction. However, multiple decarboxylating reactions exist in the cellular metabolism, showing variable activity [93,94]. This intrinsic feature needs to be included if one wishes to reproduce the respiratory quotient (RQ). Despite this, the average ethanol consumption rate is an order of magnitude smaller than the average glucose consumption rate and, as a result, an RQ value of around 1.1 is achieved. Adding another layer of detail to the biokinetic model by including lumped reactions and metabolite pools [95] might be an interesting step forward. This might be possible by combining multi-reaction models like the one used in this work with lumped metabolic models [95,96]. From another perspective, *dCO*<sup>2</sup> creates a carbonate system in the fluid and within the cell which alters the cytosolic pH and hence induces stress and increases the cellular maintenance [64,73,97] or alters the metabolism based on gas composition [98]. Based on the actual process, one can decide to include some or all of the equilibrium reaction, but this does not fall within the scope of this work since, at pH 5, more than 99% is in the form of *dCO*2. Nevertheless, using a comparatively simple approach, the results indicate that the *dCO*<sup>2</sup> gradient is rather weak compared to other species. At such levels, *dCO*<sup>2</sup> inhibition inevitably takes place at industrial scale and impacts the transcription according to recent findings [37]. Our results indicate that while fluctuations in other concentrations might be experienced by cells on short timescales, the same does not hold true for *dCO*2, where cells are exposed to high *dCO*<sup>2</sup> for long timescales. The latter requires different experimental set-ups for scale-down tests.

#### **5. Conclusions**

This work suggests that in the case of baker's yeast production, ethanol production is inevitable around the feeding point—in this case, positioned at the top of the vessel. This causes lower growth rates above the top impeller and hence hinders the overall growth rate over the bioreactor volume and is not desirable when the final product is the biomass itself. It is possible to distinguish the ethanol production due to overflow metabolism (Crabtree effect) from *dO*<sup>2</sup> limitation (Pasteur effect). Such information can prove crucial for process optimization. *dCO*<sup>2</sup> gradients might not be as pronounced as the other species, but the fact that, in both scenarios, it reaches saturation levels hints at *CO*<sup>2</sup> stripping under real industrial conditions [4]. This points out the fact that, unlike glucose, ethanol, and *dO*2, where fluctuations might trigger a stress response, *dCO*<sup>2</sup> stress is different in nature and should be evaluated by long-term scale-down experiments. The results further suggest that a real starvation region in the lower parts of the tank might not exist because of the presence of ethanol compensating for glucose shortage. Accordingly, scale-down experiments should consider this impact, even investigating the putative benefits for long-term protein formation [15].

**Author Contributions:** Funding acquisition, R.T.; Investigation, C.S.S.H.; Methodology, C.S.S.H.; Resources, R.T.; Software, R.T.; Supervision, R.T.; Writing—original draft, C.S.S.H.; Writing—review & editing, C.S.S.H., C.H., H.N. and R.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** Authors are thankful to Höchstleistungsrechenzentrum Stuttgart (HLRS) for large-scale computations and to the University of Stuttgart for providing the required funds. This work was supported by the German Federal Ministry of Education and Research (BMBF), grant number: FKZ 031B0629. C.S.S.H. is supported by ERA CoBioTech/EU H2020 project (grant 722361) "ComRaDes", a public–private partnership between the University of

Stuttgart, TU Delft, University of Liege, DSM, Centrient Pharmaceuticals and Syngulon. C.H. and H.N. are paid employees of the DSM Biotechnology Center.

**Acknowledgments:** Authors acknowledge the insightful comments of Wouter van Winden in the process of manuscript preparation.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Nomenclature**


#### **Abbreviations**


#### **Appendix A. Glucose Gradients**

**Figure A1.** Glucose gradients with (**a**) 10 g/L and (**b**) 25 g/L biomass (logarithmic colormap).

As expected, a strong gradient for glucose as the main substrate exists in both cases, as shown in Figure A1, and also aligned with previous efforts [4,25,74]. Using a similar approach to Section 4.2. for glucose (A1) gives a timescale for substrate consumption of 7 and 17 s for 25 /L and 10 g/L biomass concentration, respectively. This still falls short of the circulation time τ*circulation* of approximately 47 s, meaning that supply is slower than the demand.

$$\pi\_{s\_{\text{cons}}} = \frac{\mathbb{C}\_{\text{s}}}{(q\_{\text{sinc}} + q\_{\text{sinc}}) \times \mathbb{C}\_{\text{x}}} \tag{A1}$$

It is worth noticing the glucose concentration predicted in this work is in the same order of magnitude (Table A1), but it drops to concentrations that are below measured quantities. This could be an interesting topic for further investigations.

**Table A1.** Comparison of glucose concentration at probe location (Top: 6.35 m, Mid.: 3.9 m, Bot.: 0.97 m from the bottom) from experimental (light blue 10 g/L biomass concentration) and simulation (light orange 10 g/L biomass concentration, orange 25 g/L biomass concentration) results from the literature and this work.


#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
