Computational fluid dynamics simulation of an industrial P. chrysogenum fermentation with a coupled 9-pool metabolic model: Towards rational scale-down and design optimization

a r t i c l e i n f o

Computational fluid dynamics simulation of an industrial P. chrysogenum fermentation with a coupled 9-pool metabolic model: Towards rational scale-down and design optimization

Introduction
Due to the presence of gradients in substrate concentration (Enfors et al., 2001) and Kossen, 1984) and other process variables in industrial bioreactors, organisms are subject to temporal variations in their environment. Such variations impose stresses on these organisms (Lara et al., 2006;Neubauer and Junne, 2010;Wang et al., 2014), which may in turn affect the process yield (de Jonge et al., 2011). There are cases where extra-cellular variations appear to be advantageous (Enfors et al., 2001), but typically the impact is negative as the process is driven away from the conditions set for yield optimization (de Jonge et al., 2011;Wang et al., 2015). Being related to mixing behavior, these gradients may occur in any reactor type, and are expected to amplify upon scale-up, which may hence come with a yield loss that should be considered when judging scale-up economics. Furthermore, knowledge on the impact of bioreactor heterogeneity can be used to guide design changes to the reactor and, with genetic engineering, the micro-organism itself.
Previously, we used Euler-Lagrange computational fluid dynamics (CFD) to study the environmental fluctuations experienced by micro-organisms (called lifelines) (Haringa et al., 2016), and showed how fluctuation statistics can be acquired from such simulations to guide scale-down (SD) simulator design (Haringa et al., 2017a). These works focused on simulation and fluctuation quantification using the substrate uptake (q s ) lifeline, and did not quantitatively consider the metabolic response. When a dynamic metabolic model is available for the studied organism, coupled metabolic-hydrodynamic simulations can be used to evaluate the expected metabolic impact (Lapin et al., 2004(Lapin et al., , 2006. Combined with experiments in representative scale-down simulators, such a coupled hydrodynamic-metabolic approach can be used for: (1) scale-down verification: does a scale-down simulator result in the same metabolic response as observed in the large-scale CFD simulation? and (2) design optimization: what is the expected impact of reactor design changes or metabolic modifications based on numerical assessment? The most promising changes can then be experimentally tested in representative scale-down simulations, offering a powerful approach to rational bioreactor design and scale-up (Wang et al., 2014(Wang et al., , 2015. We numerically study five topics, outlined in Fig. 1, highlighting the different aspects of the CFD-based scale-down workflow. A penicillin production process is used as a case-study. Part I considers the coupled hydrodynamic-metabolic simulation of a 54 m 3 industrial P. chrysogenum fermentation (Haringa et al., 2016), focusing on mixing dynamics and neglecting slow processes such as biomass growth. We study the impact of mixing on metabolic variations using a 9-pool metabolic model (Tang et al., 2017). Part II focuses on the design of a representative lab-scale SD-simulator for the 54 m 3 reactor. In part III, we perform numerical verification of the proposed SD-simulator performance, first assuming ideal mixing, and second by a CFD simulation of a 3 L reactor with dynamic feed. In part IV, we discuss process optimization and propose a simple reactor alteration to improve the penicillin yield. To conclude, in part V we simulate 60 h of a fed-batch fermentation for comparison with industrial data. With this, we explore various aspects of the use of coupled hydrodynamic-metabolic modeling for process evaluation and optimization.

Methodology
All CFD simulations were conducted in ANSYS FLUENT 15:7, MATLAB 8:6:0 was used for post-processing and ideal mixing simulations.

Metabolic model
The 9-pool metabolic model for P. chrysogenum developed by Tang et al. (2017) contains 5 intra-cellular metabolite and 4 enzymatic pools, and couples to the extra-cellular substrate concentration C s and phenylacetic acid (PAA) concentration C PAA . The metabolite pools are: Glycolytic intermediates (X gly ), Amino acids (X AA ), Storage polymers (X sto ), ATP (X ATP ) and intra-cellular PAA (X PAA ), all reported in lmol=g dw with g dw being the dry biomass weight. Three dimensionless enzyme pools influence metabolic rates: X E;11 (the substrate uptake capacity), X E;32 (PAA export capacity) and X E;4 (storage capacity). The 4 th enzyme pool controls the biomass specific penicillin production rate q p (Douma et al., 2010) and is reported in mol p =Cmol x =h (van Gulik et al., 2000;Douma et al., 2010;de Jonge et al., 2011). For brevity, the mathematical model formulation is provided in Supplementary material A, together with additional information regarding FLUENT coupling. The effect of oxygen limitations has not been studied sufficiently to be included currently model (Haringa et al., 2016;Tang et al., 2017). Hence we currently assume sufficient oxygen supply in all cases; oxygen limitations will be considered in future extensions. For a model overview, we refer to Tang et al. (2017).
Model simplifications. Tang et al. developed and validated their model against a range of experimental data (Tang et al., 2017;van Gulik et al., 2000;de Jonge et al., 2011) including 360 s feast-famine cycles (de Jonge et al., 2011). These results provide confidence that the model is able to capture the impact of circulation-timescale substrate variations. However, instabilities in X ATP were encountered in our CFD simulations, which resulted from the sensitivity of the storage pool fluxes to turbulence- induced C s fluctuations on the sub-second timescale, which were not accounted for in model development (for details see supplementary material A). A structural solution of this issue requires deeper analysis of the signaling mechanism behind storage dynamics. As we currently lack the information to develop such improvements, we instead opted for a patch solution by assuming the ATP pool is in quasi-steady state, meaning the fluxes in-and out of the ATP pool balance, giving dX ATP =dt % 0 (Nikerel et al., 2012). This converts the dynamic ATP-balance in an algebraic expression: For the current non-linear kinetics, Eq. (1) was evaluated for 100; 000 randomly generated sets of intra-cellular pools. Subsequent correlation showed X ATP can be modeled as gly Þ. The model response was deemed satisfactory under all tested conditions. Further details on the approach and verification against experimental data are reported in supplementary material A).

54 m 3 reactor setup
We use the 54 m 3 reactor simulation (Haringa et al., 2016) with simplified single-phase hydrodynamics as the industrial base-case. We furthermore simulate the same case including aeration, with a superficial gas velocity of U g ¼ 0:05 m=s, measured under STP conditions. The headspace pressure in the reactor was 1:85 bar, which gives an air density of 2:4 kg=m 3 based on the log-mean pressure; the gas flowrate in the vessel was adjusted accordingly. The total domain height H t ¼ 11 m to account for broth expansion upon gassing, the gas-filled headspace is removed during parcel tracking (Haringa et al., 2017a). A discrete population balance (8 bins, 0:5-12:7 mm) with the kernels of Luo and Svendsen (1996) was employed to capture the bubble size distribution. A Sauter mean diameter d b ¼ 7 mm was observed; we lack experimental data to verify this, unfortunately. Furthermore we used the standard k À model (dispersed turbulent formulation), multiple-reference frame impeller modeling, and the universal drag model for interphase momentum exchange. Other inter-phase forces were neglected (Khopkar et al., 2003;Gunyol et al., 2009;Haringa et al., 2017a). Simulations using Casson rheology (Roels et al., 1974) diverged in volume fraction a. For simplicity, we hence set the broth rheology equal to water (Newtonian, l l ¼ 0:001). We realize this is a strong deviation from reality; we defend this assumption by observing that the measured air-broth circulation time lies in between the circulation times in pure water and airwater (Table 1), and capturing the range suffices for the current purpose. The air-water surface tension r ¼ 0:072 N=m, the turbulent Schmidt number was set to Sc t ¼ 0:2. Both single-phase and aerated simulations were conducted in a mesh with 180 periodicity. 235 Á 10 3 hexahedral grid cells were used for single-phase cases, 923 Á 10 3 hexahedral cells for aerated cases.
The gas flow number Fl ¼ Q g =ND 3 ¼ 0:1 implies the fermentor operates at the boundary of the 3-3 cavity regime and recirculation regime, where the mixing time s 95 is equal to or above that of single phase-flow, respectively (van der Lans and van't Riet, 2011). Available industrial data on the circulation time (Haringa et al., 2016) (s circ % s 95 =4, Noorman, 2011) suggests the latter; the circulation time s circ is compared to simulation results in Table 1. The single-phase and two-phase simulation under-and over-estimate s circ for aerated broth with 30%, respectively. Note the experimental value is based on a single measurement and hence comes with a significant margin of error; furthermore, transient effects may lead to a natural variability in recorded mixing times (McClure et al., 2015), introducing additional uncertainty. With the present industrial data, it is unfortunately not possible to quantify this uncertainty. We regard the single-phase and aerated simulation as a lower and upper bound mixing time scenario, with the true mixing behavior in the range. This level of accuracy suffices for our current demonstration purposes, but we stress the need for further investigation into modeling true aerated, non-Newtonian fermentation broths, and associated with that, a wider range of large-scale validation data (gas hold-up, local mixing curves and preferably local DO/substrate concentrations). Table 1 shows the gas-holdup is over-estimated compared to both air-water and air-broth experiments. This is likely an effect of model approximations, such as omitting inter-phase forces (except drag) and the empirical nature of inter-phase models/pop. balance kernels. For broth, the simplified rheology and effect of surfactants and anti-foam on the broth-water surface tension r play an additional role. Currently we are not directly interested in gasholdup, but in case oxygen dynamics are included, this aspect requires further study.

3-l laboratory reactor setup
A round-bottom vessel with a working volume of 3 L (Tang et al., 2017) is simulated for scale-down verification (452 Á 10 3 hexahedral grid cells. Geometric parameters are reported in supplementary information B. The gas flowrate applied in prior scaledown experiments is 2 L=min (0:66VVM) Tang et al., 2017)), giving Fl ¼ 0:009 with an agitation rate of N ¼ 10 s À1 (600 RPM) (Tang et al., 2017). This value is outside of the range probed in mixing experiments (van der Lans and van't Riet, 2011), but implies s 95 is similar to or slightly higher than for single-phase conditions. For simplicity, we hence ignore the effect of gas flow and model single-phase water. All walls were no-slip while the top surface had a no-shear free surface condition. Computational mixing simulations at 600 RPM yield a dimensionless mixing time Ns 95 ¼ 22, in excellent agreement with experiments (supplementary information B); the dimensionless circulation time s circ % s 95 =4 (Noorman, 2011). At high C x , the high effective liquid viscosity l l may practically lead to transitional flow, possibly increasing h 95 significantly. Previous non-Newtonian simulations of aerated lab-scale reactors did not produce realistic mixing results due to stagnant zones (Moilanen et al., 2007), and preliminary work using a low À Re n/a = not applicable; n/m = not measured.
k À model with l l ¼ 0:15 Pa s led to parcel tracking issues, with parcels sticking in the impingement point of the impeller discharge stream. We hence opted to decrease the agitation rate N to 1:67 s À1 to assess the effect of mixing time on the performance of a labscale scale-down simulator, and again assume a Newtonian fluid with l l ¼ 0:001. This approach suffices for our current interest in the qualitative effect of a significant change in s 95 ; we do stress that for predictive quantitative modeling a more realistic rheology model is required. For 600 and 100 RPM, s circ ¼ 0:55 and 3:3 s, respectively. Experimental evaluation of mixing behavior in real fermentation broths is required to comment on whether this range of s circ represents lab-scale practice.

Metabolic model coupling
The 9-pool metabolic model (Tang et al., 2017) is coupled to the Lagrangian (parcel) phase to study the response of microorganisms to environmental variations (Lapin et al., 2004(Lapin et al., , 2006Haringa et al., 2017b). In the 9-pool model the glucose uptake rate q s is subject to transporter control, where the availability of transporter (X E;11 ) is controlled by growth rate l ðh À1 Þ. This means that strictly speaking 2-way coupling is required to resolve the substrate environment, which requires simulating long timespans (OðhÞ) due to the long transporter adaptation time, and is therefore computationally expensive (see supplementary material C).
The long adaptation time allows for the assumption that the average transport capacity X E;11 is homogeneous in the fermentor.
As X E;11 ¼ f ðlÞ, its value can be estimated based on growth rate under ideally mixing conditions, l id . For the applied model, the average growth rate under dynamic conditions l was typically close to l id , and the estimated X E;11 was similarly close. A-priori estimation of X E;11 allows to use 1-way coupling, as was done in earlier work (Haringa et al., 2016(Haringa et al., , 2017aMcClure et al., 2016), which means the number of tracked parcels N p does not influence the substrate gradient and can be freely chosen. This simplification does not hold when intra-cellular dynamics affect q s at short timescales (% s circ ) (Lapin et al., 2004(Lapin et al., , 2006Haringa et al., 2017b), or when X E;11 under dynamic conditions differs strongly from the ideal-mixing assessment.
The above 1-way coupled approach was used to study mixingtimescale dynamics, assuming constant C x ; X E;11 , feed rate F and liquid-filled height H. This practically represents a chemostat cultivation, where the dilution rate D r is equal to the mean growth rate l. Parcel tracking for both 1-and 2-way coupling is conducted in FLUENT, but segregating the extra-cellular and intra-cellular reactions allows 1-way coupling to be executed after rather than during the FLUENT simulation, using MATLAB to perform the metabolic computations. Additional information regarding the practical implementation of the metabolic computations is provided in supp. mat. A. The (statistical) steady state allows to simulate Oð10Þ mixing times to acquire fluctuation statistics; lifelines of 80 h are subsequently generated to study the adaptation of q p to mixing-time dynamics (with constant C x ; X E;11 by construction) by joining together individual lifelines, exploiting the statistically-steady extra-cellular nature.
For the fed-batch simulation we use 2-way coupling to include temporal changes in C x and X E;11 , meaning the metabolic computations are conducted in FLUENT as part of the simulation. The long variation time of both parameters allows the assumption that C x and X E;11 are spatially homogeneous (X E;11 may be heterogeneous within the population, but to a same degree at every spatial location). This means that each timestep C x and X E;11 can be calculated as the parcel population ensemble average, and the local uptake rate can be computed from the Eulerian framework as Eq. (2): This simplified 2-way coupling requires the parcel number to be sufficient to capture overall heterogeneity, for which N p ¼ Oð10 3 Þ typically suffices (Haringa et al., 2016;McClure et al., 2016); full 2-way coupling would require N p ¼ Oð10 5 Þ À Oð10 6 Þ (Haringa et al., 2017b). 1-and 2-way coupling require similar computation time per hour flow-time, but 2-way coupling does require the full fermentation time to be simulated to account for changes in C x and X E;11 . A comparison of the assumptions between 1 and 2-way coupling is given in Table 2.

Overview of cases
We provide an overview of all simulations (Table 3), both conducted with CFD (FLUENT) and with the ideal or instantaneous mixing assumption (MATLAB), including the made assumptions and sections where these simulations are conducted. There is some variability in the applied timestep size Dt in FLUENT; in all cases it was ensured the particle trajectories were completed within the default accuracy settings. In all cases, glucose concentration C s was variable, and the PAA concentration was fixed at C PAA ¼ 3 mmol=kg.
As noted in Section 2.2.3, the uptake capacity q s;max ¼ k 11 Á X E;11 in the 9-pool model depends on the growth rate l. In the chemostat simulations, we aimed at l % 0:03 h À1 to maximize q p ; at this value of l, the 9-pool model predicts q s;max % 1:13 mmol=g dw =h under well mixed conditions, which is markedly lower than the q s;max ¼ 1:6 mmol=g dw =h reported by de Jonge et al. (2011), mea- MU-A were conducted with q s;max ¼ k 11 Á X E;11 ¼ 1:13 mmol=g dw =h and K s ¼ 9:8 lmol=kg.
The scale-down analysis and associated lab-scale CFD simulations (part II and III) were conducted before the 9-pool model was available, which meant we had to rely on the kinetic parameters of de Jonge et al. (2011), as in our previous work where we solely considered glucose uptake (Haringa et al., 2016). For consistency, we hence report a set of CFD simulations (TU-B, TG-B, MU-B) which use the 9-pool model, but with the uptake kinetics as published by De Jonge et al., K s ¼ 7:8 lmol=kg and q s;max ¼ 1:6 mmol=g dw =h. We note that the fluctuations in q s and the intra-cellular pools are too strong in these cases. The purpose of these simulations is to show that the intra-cellular response predicted between the industrial and lab-scale simulations matches; not to predict the metabolic response in the absolute sense.
Part I: Model response study. Part I focuses on TU-A (1-phase hydrodynamics, top feed) and TG-A (2-phase hydrodyn., top feed), to study the metabolic response to extra-cellular variations in an industrial-scale reactor with a statistically steady extra-cellular environment. As in our earlier work, a late fermentation stage Table 2 Comparison of the assumptions between 1-way and 2-way coupling method used in this work. 1-way coupling is here used for chemostat cultivation, and 2-way coupling for fed-batch cultivation. was modeled, with C x ¼ 55 g=kg and substrate feed rate F ¼ 1:23g=m 3 s (Haringa et al., 2016). The 1-way coupling approach means X E;11 remains unchanged in time. All other pools were variable, and initialized based on ideal mixing results. For consistency with part II, III, TU-B and TG-B are also reported here. The results are compared with a CFD simulation coupled with the dynamic gene regulation model of Douma et al. (1-phase, top feed, case TU-1P), and ideal-mixing simulations with both the dynamic gene regulation (ID À 1P) and 9-pool (ID-9P) model.
Part II: Scale-down design. In part II we show how a representative single-vessel SD-simulator with dynamic feed can be designed from the lifelines gathered in part I, using CFD-case TU À 1B as a basis. Two designs are proposed, with biomass concentrations C x ¼ 55 g=kg and C x ¼ 27:5 g=kg, respectively. As noted above, the uptake kinetics of de Jonge et al. (2011) were used. As in previous work (Haringa et al., 2017a), the default SD protocol is based on matching q s -lifelines between the scales.
Part III: Scale-down verification. First the performance of the scale-down protocols from part II is assessed assuming ideal mixing (cases ID-SD-27 and ID-SD-55). Next, CFD simulations of the 3 L lab scale reactor were conducted with the C x ¼ 27:5g dw =kg scale-down protocol, to study the effect of non-ideal mixing on SD performance. Instantaneous feed pulse injection was assumed in a small volume near the top surface. The hydrodynamics were frozen, but the substrate field was updated every timestep. The feed pulse scheme was supplied to FLUENT via a user defined function coupled to a lookup table. The fast mixing required time resolutions of Dt ¼ 0:002 s for N ¼ 600 RPM (case CFD-SD-600), Dt ¼ 0:01 s for N ¼ 100 RPM (case CFD-SD-100); this limited the resolved flow-time to 650 s, in which 42 feed pulses were applied. This number is too small for a proper replication of the industrialscale fluctuation statistics; therefore, scale-down performance was judged by comparing the model performance with the idealmixing response for the same 42 pulses.
Part IV: Design optimization. Industrial-scale CFD simulations were conducted with the substrate feed directly in the top impeller discharge stream (1-phase hydrodynamics), referred to as MU-A and MU-B.
Part V: Full-scale fed-batch Verification. We simulated a 60 h timespan of a fed-batch fermentation (top feed, 1-phase hydrodyn.) which was conducted in the current 54 m 3 geometry, named TU À FB, to verify model performance with industrial data which was kindly provided by the DSM biotechnology center. The simulation was started at t ¼ 10 h after the batch start C x ¼ 14 g=L. All model parameters are initialized based on the ideally-mixed 9pool model outcome for the given starting conditions. In the industrial fermentation the total broth mass increased from 36 to 46 Á 10 3 kg over the simulated timespan. However, explicitly modeling the volume change is computationally costly. As an approximation, kept the volume constant at 54 m 3 , with the hydrodynamics of MU À 1; as both impellers are submerged at all times, the change in s circ over the course of the fermentation is assumed to be minor. To compensate for the higher volume, the provided feed profile (reported in Fig. 6) was adjusted to ensure an equal feed in g=kg=s between the simulation and industrial fermentation at all times. Experimental data for q p and l were used to evaluate model performance for TU À FB, as well as an ideal-mixed simulation with the model of Douma et al., case ID À 1P À FB.

Results and discussion
3.1. Part I: Model response study

CFD simulations
We study the long-term adaptation of P. chrysogenum exposed to a strong substrate gradient. The most notable difference between TU-A/B and TG-A/B is the higher s circ for the latter, as discussed in Section 2.2.1, yielding q s fluctuations of longer duration. As q s is locally saturated in all cases, the fluctuation amplitude hardly differs. Examples of single lifelines for TU-A and TG-A are shown in Fig. 2, top panel. Fig. 2 shows the pool dynamics over an 80 h period for TU-A, TG-A, TU-B, TG-B. All cases show qualitatively similar behavior, but the higher X E;11 for TU/TG-B has a clear negative impact on q p . This illustrates the error introduced by taking kinetic parameters directly from literature, without accounting for the adaptation of q s;max to l.
Practically, q p is controlled by X gly : high X gly inhibits synthesis of penicillin producing enzyme, but it increases growth rate l which enhances enzyme synthesis. The first effect scales with X 6 gly (Tang et al., 2017), meaning that high values of X gly are highly repressive, but below-average values of X gly are hardly influential. This explains the large difference in q p between the cases, even though all cases have a nearly equal average X gly . The cases with the highest X gly buildup show the biggest q p loss. For aerated cases,the higher s circ translates to prolonged exposures to excess conditions, resulting in strong X gly accumulation. Similarly, the higher trans- Table 3 Overview of all the simulations, both CFD and ideal/instantaneous-mixing based (IDM), conducted in this work. All cases were conducted as chemostats, except for TU=ID À FB, which are a fed-batch simulations. Naming convention: T = top feed. M = mid feed (impeller discharge stream). U = ungassed. G = gassed. FB = fed-batch (2-way coupled). ID = instantaneously mixed. SD = scale-down. 9 À P indicates the 9-pool model of Tang et al. (2017) is used for metabolic coupling, 1 À P indicates the Dynamic Gene Regulation model of Douma et al. is used (Douma et al., 2010). A and B indicate which kinetic parameter values are used. SD À 100 and SD À 600 indicate agitation rates of 100 and 600 RPM, respectively. IDM 1 À P I n/a n/a n/a 1:6 7 :8 55 n/a st. st. ID-9P IDM 9 À P I n/a n/a n/a 1:6 7 :8 55 n/a st. st. ID-SD-27 IDM 9 À P II/III n/a n/a n/a 1:6 7 :8 27 n/a 0:03 ID-SD-55 IDM 9 À P II/III n/a n/a n/a 1:6 7 :8 55 n/a 0:03 CFD-SD-100 CFD 9 À P III no 1-way top 1:6 7 :8 27 5000 0:01 CFD-SD-600 CFD 9 À P III no 1-way top 1:6 7 :8 27 5000 0:002 ID-1P-FB IDM 9 À P V n/a n/a n/a var. 9:8 var. port capacity for TU/TG-B causes increased glycolytic accumulation.
The effect of both kinetics and s circ is summarized in the Damköhler number Da ¼ s circ =s rxn , where we take s rxn ¼ K s =ðq s;max C x Þ, the limit for C s ! 0. This definition for s rxn does not require specifying a value of C s , which makes it straightforward to evaluate for both experimental and CFD cases. Including the impeller-fed cases MU-A/B (part IV), a linear trend between the penicillin yield Y sp (Table 4) and Da is observed: Y sp ¼ 0:3417 À 0:0015Da (R 2 ¼ 0:97), graphically shown in supp. material D.
Within the range of fluctuations, the effect of X gly on l, while non-linear in nature (Tang et al., 2017), can be reasonably linearized. Hence, the effect of high and low X gly values on l nearly averages out: lðX gly Þ % lðX gly Þ. Only the most extreme case (TG-B) deviates from this; the very lengthy exposures to starvation conditions leads to a lower l. The data clearly shows that the duration of exposures to excess-and starvation conditions strongly impacts the metabolic response. Since these time periods are highly distributed, there is considerable heterogeneity in X gly at any given location. This feature is clearly visible in supplementary videos (available online), and is inherently not captured by black-box models that instantaneous adaptation of the intra-cellular to the extra-cellular domain.  Table 4 Comparing yields and productivity between experimental data of van Gulik et al. (2000), the black box (BB) model of Douma et al. (2010) and the 9-pool model (Tang et al., 2017) with ideal mixing assumption and the 9-pool (9P) of Tang  Video 1 X E;32 and X E;4 are hardly affected in case TU/TG-B where X E;11 was preconditioned for l ¼ 0:03 h À1 , whereas the higher uptake for TU/ TG-B causes some changes in these pools. X AA is hardly affected in all cases. The value of these pools is homogeneous within the population (supp. material E).

Experimental data and yields
The CFD results are compared with experimental chemostat data (van Gulik et al., 2000) and ideal-mixing simulations using both the model of Douma et al. (2010) (ID-1P) and the model of Tang et al. (2017) (ID-9P) in Table 4. Both models are known to under-predict q p around l ¼ 0:03 h À1 compared to steady-state experiments. As shown in Table 4, the CFD simulations show a yield loss between 18% (TU-A) and 46% (TG-B) compared to the 9-pool model with ideal mixing. The real circulation time for the 54 m 3 reactor lies in between the extremes simulated here; based on the Da-correlation a yield loss of 22% is expected for s circ ¼ 25:7 s, using X E;11 value for l ¼ 0:03 h À1 .
For demonstration, we have also coupled the model of Douma directly to FLUENT (TU-1P), which yields an extreme 85% decrease in Y sp and strong increase in Y sx (discussed in detail in Haringa et al. (2016)). These results are deemed unrealistic; the model of Douma was not designed to cope with rapid substrate concentration fluctuations, and the results show that applying the model in a situation where such fluctuations are present leads to extreme results. Returning now to the 9-pool model of Tang et al.; Although the chemostat assumption used here introduces some simplifications, we are confident the overall trends hold, making the outlined method suitable for a quick assessment of the impact of design changes on the fermentation process. The most promising cases can subsequently be studied in more detail with 2-way coupling and experimental scale-down assessment.

Part II: Scale-down design
A scale-down design analysis is conducted for TU-B. Feed protocols for a single-vessel, fluctuating feed scale-down simulation with variable pulse duration were designed based on the arcanalysis methodology proposed in Haringa et al. (2016). In contrast to earlier work, we did not divide the lifelines in regimes first; the arc analysis method was directly applied to the full (smoothed) lifelines, using a reference value q ref ¼ 0:05q s;max . This means the lifelines are divided in feast-arcs (q s =q s;max > 0:05) and faminearcs (q s =q s;max < 0:05). The arc duration s arc is registered as the time between two consecutive crossings of q ref , as graphically indicated in Fig. 3A, B. The distribution in s arc is reported in Fig. 3C. For famine arcs, we can assume negligible magnitude: q s % 0, regardless of duration. For feast arcs, the maximum q s =q s;max for each arctrajectory, called X s;max , is recorded. This gives a correlation between magnitude X s;max and duration s arc (Fig. 3D). The rationale behind q ref ¼ 0:05 follows from the results; the famine arcs show a complex distribution in s arc , but with negligible amplitude. For the feast arcs, the s arc distribution is comparatively simple, and a clear correlation between X s;max and s arc exists. Together, these statistics quantify q s -lifeline fluctuations and form a basis for representative scale-down simulation.
Representative profiles of alternating feast-famine arcs are generated from the s arc -distributions by inverse transform sampling; for each feast event, the maximum q s is retrieved from the mean s arc -X s;max correlation, Determining the feed rate F is straightforward from the mass balance, assuming an instantaneously mixed lab-scale reactor. During famine intervals, F ¼ 0 and q s % 0 by construction. The most truthful approach is to feed gradually over a period of 0:5s arc , such that the arc-shape is symmetric (Fig. 3E); this requires the lab-scale to operate at the industrial biomass concentration C x ¼ 55 g dw =kg (case ID-SD-55). Applying instantaneous feed pulse administration (Fig. 3F) relaxes this to C x ¼ 27:5 g dw =kg (ID-SD-27); the rate-of-change in q s is reduced as q s decreases over the entire period s arc . The lower C x leads to a reduced effective viscosity (Roels et al., 1974) which may facilitate practical operation. However, it must be ensured the change in rate-of-change does not result in a different metabolic response. Further decreasing C x inherently compromises either the fluctuation duration or magnitude, and thereby the representation of q s -lifelines.

Part III: Scale-down verification
In this section we assess the scale-down protocols of part II, first assuming instantaneous mixing and second using lab-scale CFD simulations. Note that instantaneous/ideal mixing in this context means the feed is immediately spatially distributed; due to the pulsed feed nature, there are temporal variations in q s .

Instantaneous mixing
Both for ID-SD-55 and ID-SD-27, 5 statistically representative lifelines were generated and analyzed. Table 5 lists the metabolic response in q p and l compared to TU-B. Additionally, we conduct a regime analysis (using the definitions of Haringa et al. (2016)) on the generated lifeline to determine the exposure to excess (E), limitation (L) and starvation (S) conditions. Case ID-SD-55 slightly over-estimates exposure to excess conditions. This results in a higher l, mildly lower q p and minor offsets in the intra-cellular pool sizes (reported in supp. material E), but overall we conclude that both cases excellently represent the large-scale simulation. The good performance of ID-SD-27 follows from the notion that the total uptake within a pulse of length t; R t 0 q s dt, is equal between the two pulse administration methods, and the turnover time of X gly is sufficiently slow to yield similar responses in X gly (see supp. material E). If the turnover time of X gly was well below s arc , the metabolic response is expected to differ between the cases, and lowering C x might not be allowed. We hence regard the possible reduction in C x as a case-depended effect, and it should be evaluated as such. Furthermore, operating at industrial C x whenever possible may avoid unforeseen responses, not captured by the metabolic model. In case no predictions regarding the metabolic response are available, a scale-down simulator should in any case aim to produce the best possible replication of the extra-cellular environment (q s -lifelines), and no compromises in C x should be made.

CFD verification
In many cases, lab-scale fermentors can be assumed as ideally mixed. However, the combination of a very short s rxn (due to a low K s ) and mixing issues due to rheological issues (Moilanen et al., 2007) could lead to spatial heterogeneity in lab-scale fermentors. To assess whether this impacts scale-down performance, CFD simulations of a SD-simulator were conducted with C x ¼ 27:5 g=kg to probe the possible impact of non-instantaneous mixing. Spatial heterogeneity relates to the Damköhler number Da ¼ s circ =s rxn .
Here, s circ % 0:5-3:3 s (for 600 and 100 RPM, respectively). As the C s field is now dynamic, we employ a more general definition of the reaction time, s rxn ¼ C s =R s ¼ ðC s þ K s Þ=ðq s;max C x Þ with C s the volume average substrate concentration. Right after pulse administration, C s ) K s and Da ( 1: this implies the pulse will be mixed before C s and thereby s rxn drop significantly, leading to a homogeneous broth and equal experiences by all micro-organisms in the domain. This is reflected in the model response for both case CFD-SD-600 and CFD-SD-600. The q s lifelines in Fig. 4B (600 RPM) and C (100 RPM) show evidence of spatial heterogeneity directly following pulse administration, which for case CFD-SD-600 rapidly wears off, meaning the lifeline under the instantaneous mixing assump-tion is retrieved (Fig. 4A). The heterogeneous period lasts longer for CFD-SD-100, but eventually the population synchronizes, and the metabolic response is hardly affected (Fig. 4D). To comment on the role of non-ideal mixing in (aerated) SD-simulators with a high liquid viscosity l l , experimental measurements are required, Table 5 Comparison of the instantaneous mixing cases ID-SD-27 and ID-SD-55 with CFD simulation TU-B. Inst. ¼ instantaneous feed, Grad. ¼ gradual feed. Dq p is reported with respect to the ideal mixing benchmark ID-9P. The last three columns report the exposure (in %) to E (excess), L (limitation) and S (starvation) conditions, based on the definitions of Haringa et al. (2016).

Case
Cx ½g=kg but the results for CFD-SD-100 imply very poor mixing is required to yield significant heterogeneity in the population, and to yield a different metabolic response compared to the pulse-profile under the assumption of instantaneous mixing. This stems positive for practical application of fluctuating-feed SD-simulators.

Part IV: Design optimization
Part I revealed that reducing the frequency of q s variations reduces the amplitude of X gly fluctuations, which reduces inhibition of q p . Cronin et al. reduced s 95 by a factor 2-2:5 by placing the feed point just below the top impeller (Cronin et al., 1994;Vrábel et al., 1999;van der Lans and van't Riet, 2011). We find s 95 ¼ 23 s (1-phase hydrodyn., MU-A/B) when the feed is placed in the top-impeller discharge stream, a 2:7-fold reduction in s 95 compared to the top feed. This exceeds expectations and may be excessively low for a true penicillin fermentation when rheology and aeration are accounted for, but we accept this result for the sake of demonstration. The pool response for simulations MU-A and MU-B is reported in Fig. 5.
Compared to TU-A, the q s -lifelines for MU-A show a lower fluctuation amplitude, and strong reduction in fluctuation duration (Fig. 5, top). This translates to much milder X gly variations that directly relate to a higher q p for MU-A/B cases (Table 6). Again, X gly and hence l remains virtually equal between the cases. The q p loss is reduced to 8:6% (with respect to ID-9P), where the topfeed case with equal s circ , TU-A, showed a yield loss of 17%. The reduced exposure to starvation conditions furthermore is observed Table 6 Comparing yields and productivity between experimental data of van Gulik et al. (2000), the black box (BB) model of Douma et al. (2010) and the 9-pool model (Tang et al., 2017) with ideal mixing assumption and the 9-pool (9P) of Tang et al. (2017).

Model
Case  to yield a higher X sto for MU À 1 cases. An alternative process improvement may be to modify increase K s by modifying the glucose transporter, thereby reducing sensitivity to C s fluctuations.
Within the current metabolic model, this also requires altering the sensitivity of the storage/release process to C s ; we did not further pursue this option within the scope of this work.

Part V: Industrial-scale fed-batch simulation
The long-term metabolic response in an industrial fed-batch reactor is simulated; the dynamic feed profile that was supplied to the simulation is reported in Fig. 6A. C x and l are well captured (( Fig. 6B and C, resp.), although an ideal-mixed simulation with model of Douma et al. (2010) (ID-FB) better captures the final 20 h. The 9-pool CFD simulation, however, performs superior in predicting the gradual reduction in q p (Fig. 6D). The initial offset results from the lower peak q p prediction by the 9-pool model The trends in intra-cellular pools (Fig. 6E) reveal major temporal changes in the pool averages (solid lines), as well as the emergence of significant heterogeneity within the population; the dashed lines in Fig. 6E represent the pool size standard deviation over 2500 tracks. The decreasing trend in all enzyme pools is a consequence of the reduction in l to 0:01 h; the drop in X E32 reduces the PAA export capacity, giving rise to a strong PAA build-up. Similarly, a buildup in X sto is observed. As before, the AA pool is least sensitive, although it undergoes some changes in later stages. The strong rise in population heterogeneity roughly coincides with the switch to a constant feed rate F % 1:6 kg=m 3 =h.
For brevity, figures further detailing the onset of and degree of population heterogeneity are reported in supp. material F. The high degree of heterogeneity in the enzyme pools may be surprising at first glance; their adaptation timescale strongly exceeds s circ , and all parcels are expected to observe highly similar C s fluctuations during the cultivation. The link between l and X E;11 plays a key role; a parcel residing in a famine zone (l % 0) for a prolonged time undergoes a reduction in X E;11 . This reduces subsequent substrate uptake q s with respect to the population average, further decreasing l and hence X E;11 , thereby amplifying the original disturbance.
A deeper analysis in supp. material F shows that parcels with a below-average X E;11 early on end in the bottom of the X E;11 distribution. A prolonged exposure to excess conditions could reverse the disturbance, but the results show that this practically rarely occurs. The further the deviation from the population average, the more unlikely recovery becomes. We do recognize that the observation that starving organisms lower their uptake capacity appears counter-intuitive; we stress this is a model prediction, that should be verified experimentally, thereby showing how coupled simulations may generate new targets for experimental investigation. The variation in all other intra-cellular pools eventually stems from the variations in X E;11 ; in the chemostat simulations of part I, where X E;11 was necessarily fixed, no population heterogeneity was observed (supp. mat. E).
The parcels with high X E;11 are the fastest growers; some acquire double the population average biomass over the cultivation time, whereas for the poorest growers l % 0 in the late process stage. As a low l has a negative effect on q p , the fastest growers are also among the best penicillin producers, whereas the poor growers mostly accumulate storage material (supp. mat. F). Whether or not the predicted degree of heterogeneity is realistic requires an experimental scale-down study where population heterogeneity is probed on the single-cell level (Zenobi, 2013;Delvigne and Goffin, 2014). The simulations predict notable heterogeneity enzyme levels, which may provide suitable targets for fluorescent marking for experimental quantification. Besides bench-scale scale-down, the use of microfluidic tools (Grünberger et al., 2014;Dusny and Schmid, 2015) with highly controllable substrate feed rates may be a promising route towards studying the effects of substrate variations on enzyme expression and population heterogeneity.

Concluding remarks
We reported on the use of coupled hydrodynamic-metabolic simulations to assess large-scale fermentation processes in five parts: (I) industrial-scale metabolic response analysis, (II) scaledown design, (III) scale-down verification, (IV) design optimization and (V) industrial-scale fed-batch analysis. Combined, these steps provide a methodology for the analysis, scale-down and optimization of large-scale fermentation processes. Combining the 9-pool metabolic model for P. chrysogenum of Tang et al. (2017) with CFD simulations of a 54 m 3 fermentor (Haringa et al., 2016) (part I), we report a predicted penicillin yield loss of 18-45%, which correlated linearly with the Damköhler number, assuming a chemostat cultivation with 1-way metabolic coupling to simplify the simulation. The yield loss resulted from level of glycolytic intermediates, relating the circulation time and substrate uptake capacity of the organism. These observations provide targets for reactorand metabolic optimization.
The arc analysis methodology of (Haringa et al., 2016) was used to design a representative single-vessel, dynamic-feed scale-down simulator, based on the q s -lifelines (part II). Numerical evaluation (part III) showed that the proposed scale-down design is predicted to accurately reflect the metabolic response recorded in the industrial reactor. Capturing the rate-of-change experienced by microorganisms on the industrial scale requires operating the lab-scale at the industrial biomass concentration C x ¼ 55 g=kg. The 9-pool model response shows that it is possible to compromise the rate of change without changing the metabolic response to some degree, allowing for a factor 2 reduction in C x . This would facilitate operation, but may induce metabolic responses for which the model does not account. Hence, we do emphasize that operating the scale-down simulator at industrial C x is preferred, especially when no metabolic response prediction is available, to ensure the best possible replication of q s -lifelines. CFD simulations of the proposed scale-down simulator with pulsed feeding showed that noninstantaneous mixing at the lab scale (assessed for circulation times of 0:55 and 3:3 s) did not compromise the metabolic response, which gives confidence in the practical application of the proposed simulator. This operational window may depend on the organism and geometry, and should be evaluated per-case.
Changing the substrate feed location in the industrial-scale fermentor to improve substrate distribution reduced the yield loss from 18:4% to 8:6% (part IV). This showcases the prospects for in silico design optimization. To conclude, we present a 60 h fedbatch study (part V) with 2-way metabolic coupling, showing good agreement in l and q p compared to industrial data, while significant intra-cellular heterogeneity was observed due to the interplay between l and the glucose transport capacity X E;11 . The results illustrate the importance of simulating fed-batch dynamics including 2-way coupling to capture population heterogeneity. We do stress this does not imply that the 1-way coupled approach is futile; it is preferred for a rapid assessment of the metabolic response to design changes. We do, however, advise that the most promising chemostat cases are subsequently simulated with 2way coupling (and/or experimentally assessed) to verify their performance when population heterogeneity is included.
Altogether, we outlined the different roles of coupled hydrodynamic-metabolic modeling in the assessment and improvement of large-scale fermentor designs. In future work, the proposed scale-down simulators are to be tested to verify model predictions; the predicted yield loss and population heterogeneity provide clear targets for assessment and model verification. There is room for improvement in both the CFD models and dynamic metabolic models, which would greatly benefit from a broader availability of industrial-scale data for verification. Such improvements act towards increasing the accuracy and reliability of the here-shown coupled CFD approach, but will not influence the methodology in itself. We believe the here-presented methodology, combined with practical scale-down simulation, opens up a new approach towards rational fermentor design and scale-up, accounting for the effect of large-scale reactor heterogeneity.