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

: 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 ﬂuid dynamics equipped with mechanistic biochemical engineering knowledge of the process, it is possible to simulate such gradients. In this work, concentration proﬁles 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 overﬂow 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.


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][2][3][4][5][6].
Different approaches to addressing such issues have been studied by bioreactor design experts [7][8][9][10][11][12]. The investigation of cellular interaction with substrate gradients has been the basis 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 2 (dO 2 ) is evaluated. In addition to this, dissolved CO 2 (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].

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

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].
Processes 2020, 8, x FOR PEER REVIEW 3 of 20 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 ( ) is evaluated. In addition to this, dissolved ( ) 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].

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

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].

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. 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 −6 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 (k l a) with published values. Flow fields served as basis for implementing mass transfer and biokinetics in subsequent steps.

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: 1.
To distinguish between these cases, the catabolic capacity to metabolize glucose aerobically serves as the threshold. In essence, if biomass specific glucose uptake q s exceeds the related oxygen demand for oxidation, Y s o ·q o meaning (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.
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 q eae are calculated as shown in (8). In essence, the min modulator compares whether oxygen demands for ethanol oxidation after glucose consumption Y o e exceed the Monod-type ethanol uptake kinetics q e .
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 Process reaction for aerobic growth on glucose (r sae ) Process reaction for anaerobic growth on glucose (r san ) Process reaction for aerobic growth on ethanol (r eae ) 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 for aerobic growth on glucose ( ) Process reaction for anaerobic growth on glucose ( ) Process reaction for aerobic growth on ethanol ( ) 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. 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.   . (16) Elemental conservation results in the following system of Equation (17): (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) 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].

Mass Transfer
The mass transfer coefficient between the two phases was considered for , , and ethanol (21). A common assumption is to estimate the mass transfer close to the equilibrium state. The gas (16) Elemental conservation results in the following system of Equation (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).
Aerobic growth on glucose : Anaerobic growth on glucose : r c san = 2r s san − 0.3r x san r e san = r s san − 0.7r x san (19) Aerobic growth on ethanol : 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].

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].
H is Henry's constant for the gas component at fermentation temperature (30 • C). k l 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 2 and CO 2 since the direction of the transport is different, which leads to (23) and (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.

Flow Field Validation
For the validation of the represented flow field, several criteria were evaluated. τ 95 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 k l a calculations (Figure 3). Under said conditions, of 190 h −1 is estimated, which agrees well with the experimental measurements [86].

Scenario I: Experimental Fermentation
For this scenario, conditions are considered as explained by previous investigations [25,86]. Accordingly, model suitability could be checked.
concentrations are expected to show high values at the bottom of the tank for mainly three reasons: first, higher hydrostatic pressure increases Under said conditions, k l a of 190 h −1 is estimated, which agrees well with the experimental measurements [86].

Scenario I: Experimental Fermentation
For this scenario, conditions are considered as explained by previous investigations [25,86]. Accordingly, model suitability could be checked. dO 2 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 2 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 2 of~3.7 × 10 −5 M (~1.2 ppm) is estimated (Figure 4a). This is an order of magnitude larger than the critical value of~4.6 × 10 −6 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].  Ethanol gradients are more pronounced than those of but less so than . 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 × 10 −5 M, Figure 4c,d).

Scenario II: Protein Production
To place the model into a more industrially relevant context, biomass was increased to 25 g/L to   Figure 4b. Interestingly, only minor dCO 2 gradients occur, varying no more than ±5% from average. Notably, CO 2 /HCO 3 − creates a buffering system that consists of around 99% CO 2 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 2 is estimated to reach mole fractions between 2.1 and 2.9%. Ethanol gradients are more pronounced than those of dCO 2 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 × 10 −5 M, Figure 4c,d).

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 m 3 ) above the top impeller, considering that 10% of saturation dO 2 . dO 2 levels below 4.6 × 10 −6 M were observed in a volume of 0.04 m 3 , 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 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. In turn, this affected the ethanol gradient twofold. First, the drop in 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. Figure A1 (Appendix A) shows the anticipated heterogeneous glucose distribution, disclosing a   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.

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 × 10 −5 M) compared to scenario I (5 × 10 −5 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.

Oxygen Gradient
A conventional approach for estimating the occurrence of gradients is the comparison of critical timescales τ for substrate supply τ s supply versus substrate consumption τ s cons . 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 τ o cons,1 and τ o cons,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.
Assuming average values, this approach theoretically indicates that assuming the non-limiting role of dO 2 for scenario I is a reasonable approximation for the whole tank, other than a small region around the feeding point, where the dO 2 concentration is slightly above 3.67 × 10 −5 M. The threshold for aerobic growth dO 2,crit for S. cerevisiae is given as 4.6 × 10 −6 M [87]. This allows us to make a distinction between the ethanol production caused by overflow metabolism or dO 2 limitation. Accordingly, no dO 2 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 m 3 of the fermenter. However, this is not the case for scenario II, where 10% of the ethanol production takes place in regions with dO 2 below the critical value. The total volume associated with ethanol production in scenario II is 0.3 m 3 .

Carbon Dioxide Gradient
Although the dCO 2 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 2 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 2 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 2 gradient is rather weak compared to other species. At such levels, dCO 2 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 2 for long timescales. The latter requires different experimental set-ups for scale-down tests.

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 2 limitation (Pasteur effect). Such information can prove crucial for process optimization. dCO 2 gradients might not be as pronounced as the other species, but the fact that, in both scenarios, it reaches saturation levels hints at CO 2 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 2 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]. 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.  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 of approximately 47 s, meaning that supply is slower than the demand. 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.

Nomenclature
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.