Supporting Information A: Modified Kinetic Model by Douma Et Al

Douma et al. [1] based several kinetic parameters for the model equations (section 2.1) on work by van Gulik et al. [2]. More recently, De Jonge et al. re-evaluated the kinetic parameters, with improved methods for determining the residual í µí° ¶ í µí± [3], yielding í µí± í µí± ,í µí±í µí±í µí±¥ and í µí°¾ í µí± as reported in table S.1. For consistency, this means í µí°¾ í µí± in eq. 2 must be refitted, yielding í µí°¾ í µí± = 8.36 ⋅ 10 −6 mols/kg The updated model parameters are shown in table S.1. In eq. 2, í µí± í µí± is proportional to the availability of limiting enzyme; the first term represents enzyme synthesis, the second represents decay and dilution. In the Herbert-Pirt equation (eq. 3) í µí¼ is negative if í µí± í µí± < í µí± í µí± í µí± í µí± í µí± ⁄ + í µí± í µí± which can be interpreted as consumption of stored carbon, or cell death due to starvation; fundamentally, a negative í µí¼ is not problematic. However, í µí¼ < 0 is a problem in eq. 2. Synthesis becomes negative for í µí¼ < 0, which is not deemed realistic and can lead to a negative í µí± í µí± giving numerical instabilities. In the second term of eq. 2, í µí¼ < 0 gives rise to an unrealistic enzyme increase. We propose two simple modifications to prevent this behaviour. The first term is set to 0 for í µí¼ < 0 to prevent negative synthesis. In the second term, í µí¼ is replaced by |í µí¼| such that í µí¼ < 0 leads to enzyme destruction rather than concentration. We stress there is no experimental backing for these choices, but we deem them defendable for the goal of our investigation. The modified version of eq. 2 is eq. S.1: í µí±í µí± í µí± í µí±í µí±¡ = í µí»½⋅max(0,í µí¼) 1 + (í µí° ¶ í µí± í µí°¾ í µí±) í µí± − (í µí°¾ í µí±í µí°¸+ |í µí¼|)í µí± í µí± (S.1) Table S.1: Parameters of the black box kinetic model applied in this study [1-3].


Supporting information A: Modified kinetic model by Douma et al.
Douma et al. [1] based several kinetic parameters for the model equations (section 2.1) on work by van Gulik et al. [2]. More recently, De Jonge et al. re-evaluated the kinetic parameters, with improved methods for determining the residual [3], yielding , and as reported in to the availability of limiting enzyme; the first term represents enzyme synthesis, the second represents decay and dilution. In the Herbert-Pirt equation (eq. 3) is negative if < ⁄ + which can be interpreted as consumption of stored carbon, or cell death due to starvation; fundamentally, a negative is not problematic. However, < 0 is a problem in eq.
2. Synthesis becomes negative for < 0, which is not deemed realistic and can lead to a negative giving numerical instabilities. In the second term of eq. 2, < 0 gives rise to an unrealistic enzyme increase. We propose two simple modifications to prevent this behaviour.
The first term is set to 0 for < 0 to prevent negative synthesis. In the second term, is replaced by | | such that < 0 leads to enzyme destruction rather than concentration. We stress there is no experimental backing for these choices, but we deem them defendable for the goal of our investigation. The modified version of eq. 2 is eq. S.1:  , the angle with respect to the positive direction. This is a simplification as the real feed tube is too small to resolve. Preliminary work showed that the size of the feed region, provided it is much smaller than the tank volume, has little effect on the shape and size of the metabolic regime distribution due to rapid initial dilution. The baffles and impellers were modelled as 0-thickness walls [4]. All walls have a no-slip boundary condition, the top was a no-shear surface. The shaft was a moving (no slip) wall with an angular velocity equal to the impeller rotation speed.

Fluid dynamics model
Here, is the particle relaxation time, the eddy length-scale and a random number between 0 and 1. The constant is linked to , and is set to 0.45 rather than the default 0.15 to reflect the change in . is a normally distributed random number. Turbulence is assumed to be isotropic. A known issue with the DRW is that near walls → 0, resulting in particles artificially 'sticking' to walls. This can be solved by using an alternative calculation for [10] or using an improved random walk model [11]. In our simulations, sticking was of minor concern as is reflected in the regime distributions (table 2) and improvements over the DRW were deemed unnecessary. 2 order upwind discretization was used for all equations [4]. The SIMPLE algorithm was used for pressure-velocity coupling. For transient simulations, second order implicit discretization in time was applied. For particle tracking, a trapezoidal tracking scheme with an automatically adapting particle timestep < Δ was applied. Subroutines to store the observed substrate concentration and track the development of and per particle were added via user defined functions (UDFs).

Mass transfer limitations:
The results of Linkes et al. [12] indicate that sub-grid concentration variations due to turbulence lead to only small fluctuations in , even for cells in the vicinity of . As such, basing the reaction rate on , ̅̅̅̅̅ leads to no significant error with respect to turbulent concentration variations either (definitely not compared to the inherent approximations and limitations of turbulence modeling). For simplicity, we currently omitted mass transfer limitations from the bulk to the organism, which may in practice lead to lower local concentrations than the local bulk concentration. If desired, such effects can be straightforwardly included in the reaction term via common closure relations.

Solution strategy and convergence criteria
First, the hydrodynamics were solved; convergence was declared when | | ̅̅̅̅ was stable within 0.01% and the residuals were below 10 −5 . Next the flowfield was frozen and the mixing time dissipation near the impeller tip [15], but decent agreement in the bulk. It is established that the − model cannot properly capture the flow close to the impeller [16]. As the impeller region is relatively small, this will have limited effect on the overall mixing. The − model is strongly preferred over more elaborate turbulence models, to keep the calculation time tractable.

Regime map sensitivity
In our work we simulated the fluid as single phase water, rather than the true aerated similarly shaped for longer times, hence only the -type is shown for brevity (at short times, there were some changes due changes around the top impeller). The distribution becomes slightly wider for case N1.17, but then remains relatively unchanged for N0.7. To explain this, we refer back to section 3.7 of the paper, where we argued that from the microbial point of view, the rate of change in followed from consumption plus local mixing. The higher , the lower the contribution of mixing, and the closer to ideal plug flow the system operates.
Consequently, the mean residence time for approaches the time required to move from / , = 0.95 to / , = 0.05 in the limit → ∞, which is 8.13 s in this case.
Indeed, the mean -residence time was found to be , 1.17 ≈ 8.3 s and , 0.7 ≈ 8.0 s which are in agreement with this reasoning. In the starvation zone, not only broadening of the distribution occurs, but the changes close to the top impeller lead to a slightly different shape altogether. For N1.17 the peak at 1 second is reduced, and slightly longer trajectories are more common because the impeller disc and regime boundary no longer overlap. The lower broth velocity changes the slope of the long-time distribution, although the limited number of counts lead to significant noise in this region. The weaker slope for the long-term distribution is even more pronounced for N0.7, where the shortterm peak is nearly completely removed due to the top impeller being fully in starvation. Also, reducing the agitation rate reduces power input which influences particle movement. As such, the comparison shown here is to illustrate possible changes only. The dominant conclusions, however, remain unaltered: concentration variations occur in the order of magnitude of the circulation time, and to properly assess the rate-of-change in conditions, ideal simulators should operate at industrial biomass concentrations.

Filtering of turbulent fluctuations
The turbulent movement of particles leads to rapid variations in / , , with rates of change in an order of magnitude faster than can be replicated in an ideal. These rapid / , variations will introduce features in the regime residence time distributions and limitation regime dynamics, that cannot be addressed feasibly with conventional scale down simulators.
Hence, these turbulent variations must be filtered to attain tractable distributions. We filtering on the overall distribution is included in table 2; with both filters applied the deviation from raw data is less than 5% for all regimes, which we deem acceptable.

Statistics of turbulent fluctuations
In the main text, we worked under the assumption that rapid fluctuations in caused by turbulence would not significantly affect the metabolism; in any case, they'll be difficult to reproduce in scale down simulators. Hence, we did not address an analysis of turbulent movements in the main text. For the interested reader, we here present briefly a methodology to assess the statistics of turbulent movements. The magnitude of turbulent movements can be determined by subtracting the smoothed signal (called , ℎ ) from the raw signal. The data displayed in figure S.4 is based on the analysis of 25000 particles. We determine the probability distribution of the turbulent movements, , around the mean that follows from the smoothed data (figure S.4 A) which is narrow near the boundaries of , ℎ and has a maximum width at , ℎ = 0.5 , . This behavior is reflected in the distribution standard deviation and is consistent with the saturation of kinetics at both extremes in . The skewness is close to 0 at low , ℎ , and increases with due to the shift from 1 to 0 ℎ order kinetics. The increase in asymmetry with increasing , ℎ is also visible when fitting the turbulent distribution with a symmetric trial function; the symmetric Gaussian-and Laplace-distribution are used. Over the entire range, the sharper peak of the Laplace distribution provides a better approximation of the observed data (figure S.4 E) which is in accordance with motion on an exponentially distributed random time (which applies for the DRW) following this distribution. To account for the effect of the shifting kinetics, an asymmetric Laplace distribution may give a better fit.