Development of dynamic compartment models for industrial aerobic fed-batch fermentation processes

Inhomogeneities in key cultivation variables (e.g., substrate and oxygen concentrations) have been shown to affect key process metrics in large-scale bioreactors. Being able to understand these gradients is hence of key interest from both an industrial and academic perspective. One of the main shortcomings of current modelling approaches is that volume change is not considered. Volume increase is a key feature of fed-batch fermentation processes. Existing models are restricted to simulating snapshots (hundreds of seconds) of industrial processes, which can last several weeks. This study presents a novel methodology that overcomes this limitation by constructing dynamic compartment models for the simulation of fed-batch fermentation processes. This strategy is applied to an industrial aerobic fed-batch fermentation process (40–90 m3) with Saccharomyces cerevisiae. First, it has been validated numerically that the compartmentalization strategy used captures the mixing performance and fluid dynamics. This was done by comparing the mixing times and the local concentration profiles of snapshot fermentation process simulations calculated with both CFD and compartment models. Subsequently, simulations of the entire process have been performed using the dynamic compartment model with kinetics. The simulation allows the spatio-temporal characterization of all process variables (e.g., glucose and DO concentrations), as well as the quantification of the metabolic regimes that the cells experience over time. This strategy enables the rapid characterization and assessment of the impact of gradients on process performance in industrial (aerobic) fed-batch fermentation processes and can be readily generalized to any type of bioreactor and


Introduction
Industrial biotechnology is increasingly being used for the production of a wide range of compounds, including pharmaceuticals, enzymes, food products and commodity chemicals [1][2][3]. Understanding the local environment of large-scale fermentation processes is key to the success of any biotechnological process, as the conditions experienced by the microorganisms directly affect the yield and productivity of the process, as well as the product quality [4][5][6][7]. Many industrially relevant processes are operated as aerobic fed-batch processes where the final volume is often 2-3 times the initial volume, meaning the conditions inside the reactor can undergo significant change throughout the course of the batch.
Large-scale bioreactors for aerobic fed-batch fermentation processes have volumes of the order 20 -1000 m 3 [8] and at this scale gradients of substrate and dissolved oxygen have been shown to occur [6,9,10]. These gradients can lead to organisms being exposed to fluctuating environmental conditions, which can impact their physiology and hence process performance [5,6,[11][12][13][14]. Therefore, substantial work aiming to understand gradients in large-scale bioreactors has been performed [12,[15][16][17]. Large-scale experimental work typically involves the installation of multiple probes which can be used for sampling and measuring the concentration of key metabolites. Such an approach has the advantage of providing direct measurements in realistic operational settings [12,14,18]. However, the spatial resolution of such measurements is often relatively low as for practical reasons only a small number of probes are used, typically located close to the wall of the reactor. Additionally, there is a cost involved in the risk of process disturbances related to malfunctioning of the sensor in feedback loops, and of contamination. For these reasons, relatively few data from industrialscale reactors are available in the open literature [6,9,12,14,19]. An alternative approach is the utilization of free-floating sensor particles which can be introduced to the reactor [20]. These particles record the local values for measurements of interest (e.g., dissolved oxygen concentration and pH) as they move throughout the vessel. The significant technical issues involved in the deployment of such particles (inclusion of sensors with low response time, ease of recovery, ability to be sterilized, sealing and robustness) means they are currently under development and are yet to be fully commercialized.
Computational methods offer an alternative approach to quantifying the conditions experienced by cells in industrial bioreactors. Here, a model of the hydrodynamics is combined with a kinetic model to predict the extent of gradients. Computational Fluid Dynamics (CFD) modelling has been used to model the hydrodynamics inside a range of bioreactor designs [16,[21][22][23][24][25]. CFD modelling provides a high degree of spatial and temporal resolution in the results, and it also allows for the examination of operational conditions far from the norm. It also provides the advantage of being able to visualize the processes occurring inside industrial-scale equipment in a way that is extremely challenging to replicate experimentally. The major limitation of CFD modelling approaches is the high computational demand involved, which is caused by both the need to model the system as a transient with short time-steps and the large mesh sizes needed to correctly model the complex internal geometry found in some bioreactor designs (i.e., presence of coils and baffles). In practice, it often takes several weeks to simulate several hundreds of seconds of operation, meaning that CFD modelling can only be used to model snapshots of a fed-batch process. Without substantial reductions in model run time, simulating the entire duration of a fedbatch fermentation process is currently not feasible.
Compartment models can be used as an alternative to CFD based models to simulate the hydrodynamics. These models consist of several ideally mixed volumes (of the order 2-70) [26,27] with the connections and flows between them being set to mimic the flow behavior of the system of interest. One of the main challenges in the construction of compartment models is the definition of the compartment volumes, flows and connections between them [28]. Many approaches are available for compartment model design [29], with hybrid CFD-compartment models being one of the most widely-used. In hybrid CFD-compartment model development, the fluid flow is first spatially divided in compartments. Subsequently, non-homogeneous flow parameters (e.g. velocity field) are used for the calculation of the exchange flows between compartments [30,31]. Due to the much smaller number of volumes and a simpler set of equations involved, run times when combined with microbial kinetic models are typically several seconds [32], which are at least 10 5 times faster than an equivalent CFD model in combination with microbial kinetics. Compartment models have been combined with kinetic models and used to simulate fixed-volume large-scale yeast (Saccharomyces cerevisiae) [33] and Escherichia coli [26] fermentation processes.
Current models (both CFD and compartment models) have the limitation that they are only applicable to fixed volume processes. This means that they cannot be used to model industrial fed-batch processes in their entirety, where feed is added. Furthermore, volume reductions, e.g. related to broth withdrawal, fill-and-draw-operations and evaporation, cannot be taken into account either. Although volume addition can be implemented in CFD models, the simulated time needs to be of the magnitude of hours to see a meaningful change in the total volume at realistic feed rates. This practice is not feasible with the current CFD models owing to the large computational cost. Existing compartment models for industrial bioreactors do not allow volume change. Addressing this limitation would allow industrial fermentation processes to be simulated in their entirety. The relatively short run time of the compartment model means it could also be used to evaluate different design and operating conditions. Hence, in this study, a numerically validated dynamic compartment model strategy which accounts for both fluid dynamics and volume change is constructed and used to describe gradient development over the course of an entire large-scale aerobic fed-batch fermentation process.

Computational fluid dynamics (CFD) modelling set-up
Computational fluid dynamics (CFD) was used to solve the hydrodynamics using the Euler-Euler two fluid approach and provided the data to be used for the compartmentalization procedure. Interphase momentum transfer was considered to arise from drag and turbulent dispersion. The Favre averaged drag model developed by Burns et al. [34] was used to model turbulent dispersion. Drag was modelled using the Ishii-Zuber correlation [35]. The standard k-ε model was used to model liquid-phase turbulence; the source terms developed by Yao and Morel [36] were included to account for bubble-induced turbulence. The dispersed phase zero equation option was used to model turbulence in the gas phase.
Here, a stirred-tank fermenter equipped with four Rushton turbines was simulated, at three different volumes (40, 60 and 90 m 3 ) which correspond to the start, middle and end of an industrial aerobic fedbatch fermentation process. A schematic of the reactor and the computational mesh used is given in Fig. 1. For simulation purposes, the mesh was split into two domain types: the tank domain (which also contained the baffles, the cooling coils, the shaft and the sparger), and the rotating domain containing the impellers and a small section of the shaft (to maintain continuity). The meshes used contained 3.41 × 10 6 , 5.19 × 10 6 and 7.10 × 10 6 elements for the 40, 60 and 90 m 3 cases, respectively.
The bubble behaviour in the model is determined via an Eulerian approach, so unlike in the Lagrangian approach which is good for low volume fractions of bubbles, explicit bubble trajectories are not followed. Instead, the bubbles are described by average volume fraction and velocity fields throughout the domain. The bubble motion depends on the forces due to drag, buoyance and turbulent dispersion. Bubbles were modelled as having a single size (d b ), which was varied to account for changes in pressure (Eq. (1)): where d b,out is the bubble diameter at the outlet, P abs is the absolute pressure and P ref is the reference pressure (1 atm). A value of 5 mm was used for d b,out based on experimental measurements [37] of the bubble size in fermentation medium. Changes in bubble size due to break-up and coalescence were not modelled. While it is possible to model such effects [38], implementation of such models leads to substantial increases in computational demand while not necessarily providing improved predictions, with a detailed discussion of such issues available elsewhere [39,40]. The liquid density used was 1050 kg m − 3 [41] , while the viscosity used was 0.692 (viscosity of water at 37 • C), 2.08 and 4.14 mPa s for the 40, 60 and 90 m 3 cases, respectively, the value being increased to mimic the presence of higher cell densities at later stages of the fermentation process [42]. Physical properties of air at 25 ℃ were used for the gas phase, except for the density which was calculated using the ideal gas law, this being done to account for the change in volume due to differences in hydrostatic pressure throughout the reactor. Air was introduced at the top face of the sparger at a fixed mass flow rate (0.40 kg s − 1 ) for all simulations using an inlet boundary condition. The top of the reactor was modelled using the degassing boundary condition, while all other surfaces were modelled as walls (having free-slip for the gas-phase and no-slip for the liquid). The oxygen transfer rate (OTR) was calculated according to: where k L is the liquid film mass transfer coefficient, a is the interfacial area per unit volume, O * is the saturation oxygen concentration and O is the dissolved oxygen concentration. Here k L was calculated as a function of the local energy dissipation rate using the correlation of Lamont and Scott [43], where a value of 2.42 × 10 -9 m 2 s − 1 was used for the diffusivity of oxygen in water [44]. The interfacial area (a) was calculated based on the local values of the bubble size and gas volume fraction (α), assuming the bubbles were spherical: Henry's law was used to calculate the value of O * : where y is the mole fraction of oxygen in the gas and H O is the Henry's law constant (790.6 atm L mol − 1 ) [44]. In this work, the gas phase oxygen mass balance was not solved, meaning that the value of y has been fixed (0.2095). This assumption was made in order to simultaneously run multiple scenarios (top and bottom feeding). Consequently, this led to an over-estimate of the OTR. This has been estimated to be up to 15%, assuming that the mole fraction of oxygen at the outlet is 0.15 and using the logarithmic mean concentration difference.
Here, the average volumetric power input (P i /V L ) was fixed at 0.9 kW m − 3 for all three volumes simulated, and this corresponded to rotational speeds of 60, 70 and 80 rpm for the 40, 60 and 90 m 3 cases, respectively. The average power input per volume was calculated as shown in Eq. (5): where T Imp corresponds to the torque of each impeller, N is the agitation speed, V T is the total volume and α is the gas volume fraction.
Its value was time-averaged for 90 s, during the transient averaging of the hydrodynamic variables. If case studies with time-varying volumetric power inputs were simulated, the flows between compartments should be adjusted as a function of the power input per volume. The transient rotor-stator model was used to model the interface between the two domains. A time step size of 10 ms was used, to correctly capture the behaviour of the system.
Mixing times were determined by introducing scalars at two locations corresponding to the top and bottom feed points ( Fig. 6A-C). Three tracers were introduced at each feed point for a period of 1 s at a flow rate of 312 kg s − 1 . A value of 8.5 × 10 -10 m 2 s − 1 was used for the kinematic diffusivity of the tracer (corresponding to the value for glucose in water [44]). The tracer concentration was monitored at two points (as shown in Fig. 6A-C) and the mixing time was defined as the time for the tracer concentration to reach ± 5% of the final, equilibrium value. Mixing times were also calculated with the compartment models at 40, 60 and 90 m 3 . The same flow rate and criteria as with CFD modelling were used. The tracer addition and monitoring compartments are shown in Fig. 6A-C.
The CFD simulations used in this work were set-up and solved with Ansys CFX 19.2. The high-resolution advection scheme was used, with a first order scheme being used for the turbulence equations, the secondorder backward Euler scheme was utilized for the transient terms. Here, we have used the coupled solver with the volume fraction equations included. Generally, 5-10 coefficient loops were needed to reduce the RMS residuals below the target value (1 × 10 -4 ). Simulations were performed for 30 s. After this time, transient averaging was started, and the simulations were run for an additional 90 s. These transient average results were used in the compartmentalization process. After a total of 120 s, microbial kinetics and the tracers were implemented, and the simulations were run for an additional 30 s. Subsequently, transient averages of the glucose, dissolved oxygen and ethanol concentrations were taken for 60 s. When the mixing time was longer than 90 s, the CFD models with tracers were run until the target homogeneity level (±5% of equilibrium value) was reached. To provide an idea of the elapsed real time needed, for the smallest mesh (40 m 3 ), an average run time of 1 day was needed to simulate 3 s when kinetics were implemented using 24 CPU cores (Intel Xeon Gold 6126, 12 core, 2.60 GHz).
When plotting the CFD results with contour plots, the data for the impeller domain have been excluded because the transient averages were taken in a rotational frame of reference. Then, when visualizing the results in a stationary frame of reference, a sharp interface between the tank and impeller domains was observed. Thus, the impeller data were excluded to avoid the misrepresentation arising from this modelling artifact.

Kinetic modelling
A kinetic model for Saccharomyces cerevisiae with model parameters from the literature was implemented in the CFD and compartment models. The model parameters are summarized in Table 1.
The model structure follows the approach of the Sonnleitner and Käppeli model for S. cerevisiae [46]. The model assumes the occurrence of five metabolic regimes (glucose starvation, oxidation, overflow, oxygen limitation and glucose starvation and oxygen limitation) depending on the values of the specific uptake rates of glucose (q G ) and oxygen (q O ). The distribution of regimes with the specific rate thresholds is depicted in Fig. 2. The thresholds between glucose or oxygen sufficient and limited regimes are the maintenance levels for glucose (m G ) and oxygen (m O ), respectively. The boundary between oxidation and overflow regimes is the critical specific rate for glucose uptake (q G,crit ), which corresponds to the value of q G when cells are grown at the critical specific growth rate (μ crit ). The switches between metabolic regimes are assumed to be instantaneous and have been modelled using a hyperbolic tangent approach.
The specific growth rate on glucose (μ G ) and ethanol (μ BP ) follow Monod kinetics as shown in Eqs. (6) and (7), respectively: The expression for the specific rates for glucose and dissolved oxygen uptake, and for ethanol formation and assimilation are summarized in Table 2.

Implementation of reaction kinetics to CFD modelling
The kinetic model for S. cerevisiae was implemented in the CFD models. Concentrations of ethanol, dissolved oxygen and glucose were modelled as scalars, with source and sink terms corresponding to transfer/production and consumption being included. A linearization coefficient (C Si ) was specified for sink terms to improve convergence, whose mathematical expression is shown in Eq. (21): where S i corresponds to the source term of a variable i. Biomass growth was neglected due to the relatively short simulation time scale (90 s). Glucose was used as feed, and it was introduced at two locations as different case studies (i.e. top and bottom feed) ( Fig. 6A-C). All concentrations were calculated per liquid volume. The diffusivity values for glucose, oxygen and ethanol used were 8.5 × 10 -10 , 2.42 × 10 -9 and 1.49 × 10 -9 m 2 s − 1 , respectively [44]. The biomass concentrations used were 20, 58 and 79.4 kg m − 3 for the 40, 60 and 90 m 3 simulations, respectively, and the glucose feed rate was 226, 278 and 408 kg h − 1 for the 40, 60 and 90 m 3 simulations, respectively. These were calculated Table 1 Values of the model parameters of S. cerevisiae with their references. The parameters marked with an asterisk (*) have been calculated based on the references.  ), three regimes can arise depending on the value of the specific glucose uptake rate (q G ). If glucose is insufficient, i. e. when the specific glucose uptake rate is lower than the maintenance coefficient for glucose (q G < m G ), the cells will undergo glucose starvation.

Table 2
Equations of specific rates for growth (μ), glucose uptake (q G ), ethanol formation and re-assimilation (q BP ) and oxygen uptake (q O ) of the kinetic model used.

Regime Conditions
Oxygen limitation Glucose starvation and oxygen limitation with the dynamic simulations assuming ideal mixing (section 2.5.3).

Dynamic compartment model set-up
The approach used for the set-up and solution of the dynamic compartment models is outlined schematically in Fig. 3.
In the first step of the process, the CFD results at different volumes (40, 60 and 90 m 3 ) are compartmentalized based on the axial and radial velocities using a previously developed [28] automated methodology. The methodology uses the data of half a 2-D plane to define the compartment volumes, flows and connections of the compartment model. The 2-D plane selected in this study corresponds to that shown in Fig. 1D. By using this methodology, it is assumed that a 3-D geometry can be simplified to a 2-D plane. This is reasonable for symmetric systems around the central vertical axis, as is the case for those simulated in this study. Furthermore, the velocity data used correspond to transient average results, which are also symmetric around the central vertical axis.
The second step is the development of a library of compartment maps for those volumes not directly simulated by CFD modelling (intermediate volumes). It was assumed that the flow behaviour would not change significantly for volume increases below 10%. Thus, each intermediate volume was approximately 10% larger than the previous one, meaning the fed-batch fermentation process was divided into eight phases. For each intermediate volume, the CFD information used to develop the compartment map corresponded to that of the next higher volume from which CFD data were available. For example, if a compartment map for a volume of 53 m 3 is being developed, the CFD simulation results that are used for the compartmentalization are those of the 60 m 3 simulation, while if a compartment map for 66 m 3 is being built, the CFD simulation results correspond to those of the 90 m 3 simulation. For the development of compartment maps of different volumes other than those of the CFD simulations, the CFD extracted data corresponding to a larger volume than that to be re-compartmentalized were excluded. The target parameter used to select the excluded data was the height of the two-phase mixture. With the utilization of this methodology, flow particularities between different intermediate volumes (i.e. within each fed-batch phase) may be simplified. For example, the compartment model may not be able to capture the flow behaviour at the surface when a new impeller is partially flooded. These simplifications may lead to the occurrence of flow artefacts in some compartments. These potential issues have not been investigated in this publication, as the main flow behaviour was kept across the collection of compartment models, and no inconsistent hydrodynamic or kinetic results were observed. They may need to be considered if significant mismatches with validation data are found, e.g. by increasing the number of intermediate volumes simulated with CFD modelling.
Subsequently, to simulate the increase of mass in the compartments, the volume of the top compartments is increased. The kinetic model state variables (e.g. glucose mass) at the top compartments are assumed to get diluted according to the volume increase. The compartmentalization strategy used [28] initially divides the geometry of a twodimensional surface into a grid, where the surface used is shown schematically in Fig. 3, Step 3. The number of elements in the vertical direction is variable and depends on the height of the two-phase mixture (ranging between 32 and 61 elements), while the horizontal direction is fixed (15 elements) as the radius of the vessel is constant. The extension of the top compartments is performed by assuming that they increase in the same compartment proportion as that of the top fluid surface (Fig. 3, Step 3).
The final step is the calculation of the overall mass transfer coefficient (k L a) and the dissolved oxygen concentration at saturation (O * ) in each compartment over each fed-batch phase. To this end, an empirical approach which calculates the k L a and O * values in each compartment as a function of the total volume has been used. The resulting k L a and O * values in each compartment can only be used for the specific case study for which the CFD simulations are performed. Thus, if the equipment is changed (e.g. addition of an impeller, utilization of different reactor configurations), the CFD simulations need to be re-done. Then, new compartment maps need to be constructed, and different k L a and O * values in each compartment need to be calculated. First, the volumeweighted average values of k L a and O * of the CFD results at all different volumes are calculated. The resulting values are an average k L a of 300 ± 23 h − 1 taking the three volumes into account and a linear correlation with the total volume for O * values (slope = 4 × 10 -5 g kg − 1 m − 3 , intercept = 0.0082 g kg − 1 ). Subsequently, for all volumes of a fedbatch phase, the fraction (β) of each variable in a certain volume (V T ) with respect to its value for the initial volume (V T0 ) is calculated following the approach from Eqs. (22) and (23): Then, the volume-weighted average k L a and O * values in each compartment at the beginning of each fed-batch phase from the CFD plane data (k L a Comp (V T0 ) and O * Comp (V T0 )) are multiplied by the fractions calculated in Eqs. (22) and (23) to give the k L a and O * values in each compartment over the fed-batch phase (k L a Comp and O * Comp ). Moreover, a correction is performed to ensure that the volume-weighted average of the k L a and O * values of all compartments correspond to those calculated considering the entire CFD volume rather than only the data from the plane. The mathematical expression of the correction performed is shown in Eqs. (24) and (25):

Simulation of aerobic fed-batch fermentation processes 2.5.1. Dynamic compartment models
The dynamic compartment models with the kinetic model of S. cerevisiae with different top and bottom feeding positions ( Fig. 6A-C) were implemented in MATLAB® 2019a as a set of ordinary differential equations (ODEs), where the solver ode15s was used for solving the mass balances of biomass (M X,Comp ), glucose (M G,Comp ), ethanol (M BP,Comp ) and dissolved oxygen (M O,Comp ) in each compartment. The absolute and relative tolerance values were set to 1 × 10 -5 and 1 × 10 -6 , respectively. The mathematical expressions describing the mass changes of these state variables are summarized in Eqs. (26)- (30). For the description of the glucose mass in each compartment, Eq. (27) was used for the feeding compartment, while in the other compartments it was described with Eq. (28). and q O,Comp are the specific rates for glucose and oxygen consumption, respectively, q BP,Comp is the specific rate for ethanol uptake and production, M L,Comp is the liquid mass in a compartment, and k L a Comp and O * Comp correspond to the overall mass transfer coefficient and the oxygen concentration at saturation in a compartment, respectively. F Feed corresponds to the feed rate calculated assuming ideal mixing (section 2.5.3) and G Feed is the glucose concentration in the feeding solution, which has a value of 305 g kg − 1 .
To account for the change in total liquid mass (M L ), an additional ODE was implemented as shown in Eq. (31): Hence, the rate of change in the liquid mass depends on the feed rate, the evaporation rate (F Evap ), and the summation of the flow of oxygen transferred from the gas to the liquid phase (F O2 ) and of carbon dioxide removed (F CO2 ) in all compartments.
The evaporation rate was calculated by estimating the mass of water at the inlet and at the outlet of the tank, following the approach from Mears et al. [41]. The mathematical description of F Evap is shown in Eq. (32).
where Q G corresponds to the volumetric air flow rate (NL h − 1 at 0 • C and 1 bar), p * out and p * in correspond to the saturated vapour pressure (in Pa) at the outlet and at the inlet, respectively, M w,H2O is the molecular weight of water, R is the ideal gas constant and T out and T in are the temperature at the outlet (37 • C) and at the inlet (40 • C), respectively. The values of p * in and p * out are calculated using the approach from Bolton [48] (Eqs. (33) and (34)): where RH corresponds to the relative humidity (79%) [49]. The flow of oxygen transferred from the gas to the liquid phase in each compartment (F O2 ) was calculated from the oxygen transfer rate, following the approach from Eq. (35).
The carbon dioxide removed in each compartment (F CO2 ) was estimated from the oxygen uptake rate. It was assumed that for each mole of oxygen consumed, one mole of carbon dioxide was produced. This assumption is valid under glucose starvation and oxidation regimes (Fig. 2). Nevertheless, it is applied for all cases because these regimes are the most experienced by the cells based on the simulation results (Figs. 11 and 12). Hence, it is not expected that this calculation will have a significant impact in the simulation output when the other metabolic regimes are expressed. The mathematical expression is shown in Eq. (36): where M w,CO2 and M w,O2 correspond to the molecular weight of carbon dioxide and of oxygen, respectively.

Fixed-volume compartment models
The kinetic model for S. cerevisiae was implemented to the fixedvolume CFD-based compartment models of 40, 60 and 90 m 3 developed in the first step of the dynamic compartment model development strategy (Fig. 3) for their comparison with the CFD results with kinetics. Both top and bottom feeding position case studies were simulated ( Fig. 6A-C). In that case, Eqs. (27)- (30) were implemented in MATLAB® 2019a to describe the change of glucose, dissolved oxygen and ethanol masses in each compartment. The ODE describing the change of biomass in each compartment was not included to ensure results were consistent with the CFD simulation set-up. The biomass concentration and glucose feed rates simulated are equal to those used with CFD modelling.

Simulations assuming ideal mixing
Dynamic simulation of aerobic fed-batch fermentation processes assuming ideal mixing were used for calculating the glucose feed rate and the biomass concentration values that would be used to perform CFD, fixed-volume and dynamic compartment model simulations with kinetics. The methodology followed is summarized in Fig. 4.
To this end, the ODEs describing the changes in biomass, glucose, dissolved oxygen, ethanol and liquid broth mass were implemented to MATLAB® 2019a as in Eqs. (26), (27), (29)-(31) assuming the presence of only one compartment. The values of F Evap , F O2 and F CO2 were also calculated as in Eqs. (32)- (36). The feed rate (F Feed ) was calculated in two different parts. First, it was set as an exponential feed (F Exp , Eq. (37)) to keep the specific growth rate at a constant level (μ set ). The μ set value was chosen to be 0.15 h − 1 to prevent the onset of overflow metabolism, which is triggered at a specific growth rate value of 0.25 h − 1 (μ crit ) [45]. When the DO level reached 20% of its saturation value, the control strategy changed. A proportional-integral (PI) controller that manipulated the feed rate (F PI , Eq. (38)) to keep the DO concentration at 20% (O set ) was then activated.
In Eq. (38), F Bias corresponds to the feed rate value at the end of the exponential phase, K c is the gain of the PI controller and I is the integral component. When the PI controller is activated, an additional differential equation for the integral component I was implemented as shown in Eq. (39): where τ corresponds to the time constant of the PI controller. The values for K c and τ correspond to 10 6 kg 2 g -1 h − 1 and 0.1 h, respectively.
The values of k L a and O * used were calculated from the CFD results of hydrodynamics after 90 s of time-averaging by calculating their volumeweighted average values at the three different volumes simulated. As previously underlined, a constant value of 300 h − 1 was used for the k L a and a total volume-dependent correlation (slope = 4 × 10 -5 g kg − 1 m − 3 , intercept = 0.0082 g kg − 1 ) for the O * . The initial conditions had the same values as those used in the dynamic compartment model simulations.

Fluid dynamics and mixing performance
First, the fluid dynamics and mixing performance of the collection of compartment maps developed at other volumes than those simulated in CFD modelling were examined to learn how the physical part of the process would influence the simulation of aerobic fed-batch fermentation processes with dynamic compartment models. The collection of compartment maps generated for the case simulated here is shown in Fig. 5.
The total number of compartments in each compartment map ranged between 14 and 27, linearly increasing with the total volume. This observation is in line with the impeller configuration, which leads to a series of radial re-circulation loops adjacent to the impeller [50]. With increasing volume and number of impellers, the number of recirculation loops and compartments consequently increases. It is also observed that the distribution and size of compartments does not change significantly between all compartment models at different volumes, including those developed from different CFD results (at 40, 60 and 90 m 3 ). From a practical perspective, the distribution and size of compartments at the same coordinates are kept across volumes and, with volume addition, new compartments are added at the top. This observation strengthens the strategy used to develop compartment models of Fig. 4. Workflow of the methodology used to calculate the feed rate and the biomass concentration for the CFD and compartment model simulations combined with kinetic models. other volumes than those simulated in CFD modelling as no significant changes are observed in compartment distribution and size between the different volumes.
For the dynamic compartment model to offer meaningful predictions of the process performance it must first give reasonable predictions of the mixing and oxygen mass transfer. This is because the microbial kinetics depend on both physical processes and if their underlying description is incorrect, it is not likely the model will offer meaningful predictions. Ideally, comparison would be made with data from an industrial-scale process. However, no relevant measurements at such a scale are available in the open literature. Hence, comparison has been made with the CFD models. These have already shown good agreement  the Table). The top feeding position is shown with a red sphere and a blue arrow for the CFD and CMs, respectively. The bottom feeding position is shown with a green sphere and a grey arrow for the CFD and CMs, respectively. The monitoring points are shown with yellow spheres or circles. The mixing times estimated for the top (TMP) and bottom (BMP) monitoring points are summarized in subplots D and E, respectively. The error bars show the standard deviation of the tracer triplicates in CFD modelling. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) with mixing time predictions of pilot-and production-scale bioreactors [16,18,51]. Such a comparison between CFD and compartment model data is also interesting as it allows for an increased understanding of the trade-off between predictive ability and computational demand involved in the two modelling approaches.
Oxygen transfer is not discussed because the CFD data were already utilized as an input for the calculations of k L a and O * in each compartment. Thus, the following discussion focuses on comparing the mixing times calculated using CFD and compartment models, whose results are shown in Fig. 6. The results agree with a correlation from the literature (36 -181 s, Eq. (40)) [52] describing the mixing time in stirred tanks with multiple RTD impellers: (40) where N is the agitation speed, Po is the power number (corresponding to 5.5 for RTD impellers [53]), D is the impeller diameter and H is the height of the two-phase mixture. It was found that there was good agreement between the predictions of the CFD and CFD-based compartment models for volumes of 40 and 60 m 3 , while for the 90 m 3 volume the compartment models underpredicted the mixing time by a factor of approximately two. Differences in the mixing time as a function of the feeding location were only observed for the 40 m 3 case, and the location of the monitoring point did not appear to have large impact on the calculated mixing time. It is encouraging that the compartment models were able to predict mixing times of the same order of magnitude, given the considerable amount of simplification involved in the development of such models. Whether the difference in mixing time values between CFD and compartment models is sufficiently high to create a miss-match when predicting the magnitude and occurrence of gradients needs to be assessed by extending both models with microbial kinetics and comparing their results. This is done in the following section.

Comparison between fixed-volume CFD and compartment model simulations
The next step to evaluate the predictive ability of compartment models resembling the fluid dynamics was their extension with microbial kinetics to model the local concentrations of glucose, dissolved oxygen and ethanol. These predictions were subsequently compared with results from CFD modelling with kinetics at different points of the fermentation process. The predicted local concentrations with both modelling approaches are shown in Fig. 7 for the largest volume (90 m 3 ). The same assessment was done for the other volumes (40 and 60 m 3 ) and yielded similar results in terms of agreement between CFD and compartment model predictions.
As expected, the highest local glucose concentration was found in the feeding locations ( Fig. 7A and B). It was observed that a relatively large portion of the reactor volume had relatively low (<0.01 kg m − 3 ) local glucose concentrations for both top and bottom feeding locations. These results suggest that glucose starvation may be the most significant metabolic regime occurring in the reactor. Plots of the dissolved oxygen concentration are given in Fig. 7C and D. It was found that the lowest dissolved oxygen concentrations were also in the feeding points, this being a result of the kinetic approach used where the oxygen uptake rate is proportional to the glucose uptake rate (Ref to Eq. (19)). Finally, the local ethanol concentration values were evaluated. An ethanol concentration gradient with highest local values (up to 0.16 kg m − 3 ) in the feeding location was observed. The values of the local ethanol concentration were also correlated with those of the local glucose concentration (Ref to Eqs. (17) and (18)).
The level and location of the local concentrations of glucose, dissolved oxygen and ethanol suggest a good agreement of the metabolic activities that the cells would undergo between CFD and CM modelling strategies. The low glucose concentration levels indicate that glucose starvation is the most significant metabolic regime occurring in the reactor, at levels of 40 and 60% v/v and 68 and 84% v/v for CM and CFD simulations with bottom and top feeding positions, respectively. Then, oxidation and oxygen limitation have the higher volume fraction levels, with higher oxidation levels with CFD modelling in comparison with compartment modelling. With both modelling approaches, low levels (<5% v/v) of both overflow and combined glucose starvation and oxygen limitation are found. Thus, compartment models overestimate the degree of metabolic heterogeneity in comparison with CFD models. Despite the differences, both approaches lead to the prediction of the same metabolic regimes with strong similarities in the metabolic regime ratios calculated.
Encouragingly, the predicted profiles of the local glucose, dissolved oxygen and ethanol concentrations and the metabolic regime expression levels calculated using both CFD and CM approaches were in good agreement. Examination of Fig. 7 also illustrates the tradeoffs involved with both modelling approaches. The CFD predictions have a high degree of spatial resolution, however as previously noted these results have a substantial computational demand. Contrastingly, the compartment models have lower spatial resolution, but have negligible computational demand in comparison with CFD models. The compartment models offer a reduction of up to six orders of magnitude of the run time needed to simulate a snapshot of the process in comparison with CFD modelling. As previously noted, the compartment model predictions are in good agreement with those from the CFD model, suggesting that the key characteristics of the hydrodynamics are being captured by the compartmentalization approach employed in this work. Hence, the dynamic compartment model was used to quantify the performance of two full industrial aerobic fed-batch fermentation processes.

Dynamic compartment model simulations
The overall predicted biomass, glucose, ethanol and dissolved oxygen concentrations for a case assuming ideal mixing (IM), as well as the dynamic compartment model with top (TF) and bottom (BF) feed locations are shown in Fig. 8. More detailed contour plots showing the concentration of glucose, ethanol and dissolved oxygen in each compartment are presented in Fig. 9 for top feeding and in Fig. 10 for bottom feeding. A plot showing the metabolic regime experienced in each compartment is given in Fig. 11 , while the volume fractions of each regime over the course of the fermentation processes simulated here are shown in Fig. 12.
As expected the gradients in glucose and dissolved oxygen concentrations tended to increase with time (see Figs. 9 and 10), this being due to both the higher biomass concentration (leading to higher uptake rates) and the increased liquid volume (leading to increased transport times, see Fig. 6). These gradients in turn led to increased heterogeneity in the metabolic regimes experienced by the cells (Figs. 11 and 12).
As shown in Figs. 11 and 12, at the start of the process the entire volume experienced oxidation and as the process progressed, the volume experiencing glucose starvation, as well as the volume experiencing oxygen limitation, increased. The final oxidation volume fractions were 26 and 13% for the simulations with top and bottom feeding positions, respectively. The occurrence of starvation was also reported in a CFD study of an industrial (22 m 3 ) S. cerevisiae aerobic fed-batch fermentation process [16]. The occurrence of glucose starvation is in line with the results shown in Figs. 9 and 10, where low glucose and dissolved oxygen concentrations were found in significant portions of the reactor volume, particularly as the fermentation process progressed and the liquid volume increased. As shown in Figs. 8-10C, there was a high local ethanol concentration (>0.1 kg m − 3 ). The model did not predict the occurrence of overflow metabolism, the production of ethanol being primarily due to oxygen limitation. These results are consistent with those previously discussed where the dissolved oxygen concentration was found to be The predicted biomass concentration for the compartment model was lower than for the ideally mixed case, with the bottom feed having a slightly lower biomass concentration than the top (Fig. 8A). Results from Fig. 8 can be used to calculate relevant process metrics. For example, the yield of the fermentation process can be determined. In this instance it was calculated that the overall biomass yield on substrate was 0.33 and 0.32 g biomass per gram of glucose for the top and bottom feed positions respectively. A yield of 0.39 g biomass per gram of glucose was calculated for an ideally mixed case using the conditions applied to the compartment model case. The calculated yield for the dynamic compartment model case was 82-85% of the ideally mixed case, corresponding to a yield reduction of 15-18%. George et al. [5] compared both laboratory and production scale baker's yeast fermentations, finding that the biomass yield on substrate was reduced by approximately 7% at the production scale. These results are in reasonable agreement with the predictions of this work, particularly given the different reactor configurations used.
The results shown in Figs. 9-12 clearly illustrate the advantage of the dynamic compartment model approach as it allows for the gradients in the reactor to be visualized throughout the entire duration of the fermentation process. Existing approaches are only capable of providing snapshots of the process, whereas the current approach can model the fed-batch process in its entirety. Another advantage of the dynamic compartment model approach is that it has a relatively short run time (4-7 min for the cases studied with 1 CPU core), and hence can be rapidly used to examine different reactor designs (e.g. varying the substrate feed location) or operating conditions (e.g. control systems). While this work has examined the production of baker's yeast the approach outlined can also be readily generalized to other microorganisms. , ethanol (C) and dissolved oxygen (D) concentration evolution over the dynamic simulation of the aerobic fed-batch fermentation process with Saccharomyces cerevisiae ranging between 40 and 90 m 3 using dynamic compartment modelling with top (TF, red) and bottom (BF, blue) feeding, and assuming ideal mixing (IM, black). The concentration values were calculated by adding the masses of each variable in all compartments and dividing them by the total liquid broth volume. The dissolved oxygen concentration was normalized with respect to its volume-weighted average saturation concentration calculated from the CFD results. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Conclusions
The aim of this work was to develop a dynamic compartment modelling approach which could be used to simulate fed-batch fermentation processes in their entirety taking both the fluid dynamics and volume change into account. Such an approach is necessary as many industrially relevant fermentation processes undergo significant volume change throughout the process and existing modelling approaches are not capable of accounting for this. In this work a dynamic compartment modelling approach was developed and the predictions from the resulting constructed models were compared with results from CFD simulations. Both models used kinetics for the growth of S. cerevisiae, a widely used industrial organism.
It was found that the compartmentalization approach used gave predicted mixing times (16-119 s) in good agreement with those predicted by the CFD models (27-198 s) for the range of operating conditions examined. Similarly, predictions of the glucose, dissolved oxygen and ethanol concentrations were in good agreement for both modelling approaches. Results from the CFD model had a much higher spatial resolution, however the tradeoff involved was substantially longer run times. For example, for the least computationally-demanding case study (40 m 3 ), 90 s of CFD with kinetics simulation time corresponded to approximately 30 days of real-world time with 24 CPU cores, while solving the compartment model took 2 s with 1 CPU core. The compartment model was then used to simulate a fed-batch process from beginning to end, it was found that the biomass yield decreased by 15-18% compared with the ideally mixed case. Using the dynamic compartment model, it was possible to visualize the metabolic regimes experienced by the microorganisms throughout the duration of the fedbatch process, something which is impossible to achieve using alternative approaches. While the approach here has been applied to S. cerevisiae and a stirred tank with Rushton turbine disk impellers it can be readily generalized to other microorganisms and reactor configurations.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.