Lagrangian Trajectories to Predict the Formation of Population Heterogeneity in Large-Scale Bioreactors

Successful scale-up of bioprocesses requires that laboratory-scale performance is equally achieved during large-scale production to meet economic constraints. In industry, heuristic approaches are often applied, making use of physical scale-up criteria that do not consider cellular needs or properties. As a consequence, large-scale productivities, conversion yields, or product purities are often deteriorated, which may prevent economic success. The occurrence of population heterogeneity in large-scale production may be the reason for underperformance. In this study, an in silico method to predict the formation of population heterogeneity by combining computational fluid dynamics (CFD) with a cell cycle model of Pseudomonas putida KT2440 was developed. The glucose gradient and flow field of a 54,000 L stirred tank reactor were generated with the Euler approach, and bacterial movement was simulated as Lagrange particles. The latter were statistically evaluated using a cell cycle model. Accordingly, 72% of all cells were found to switch between standard and multifork replication, and 10% were likely to undergo massive, transcriptional adaptations to respond to extracellular starving conditions. At the same time, 56% of all cells replicated very fast, with µ ≥ 0.3 h−1 performing multifork replication. The population showed very strong heterogeneity, as indicated by the observation that 52.9% showed higher than average adenosine triphosphate (ATP) maintenance demands (12.2%, up to 1.5 fold). These results underline the potential of CFD linked to structured cell cycle models for predicting large-scale heterogeneity in silico and ab initio.


Introduction
The physiological state of bacterial cells is strongly dependent on the surrounding conditions. As outlined in Müller et al. [1], external stress is a key factor inducing the formation of population heterogeneity, which differs according to growth phenotypes and cell cycle patterns. Moreover, concentration fluctuations occurring under large-scale mixing conditions have a measurable influence on growth and production yield [2][3][4]. Accordingly, homogeneity of the bacterial population may be affected, yielding subpopulations that co-exist next to each other [1]. Makinoshima et al. [5] observed five and ten cell populations of Escherichia coli during exponential growth and the subsequent stationary phase, respectively. For Pseudomonas putida, steady-state chemostat cultivation revealed that industry-like stress conditions induced changes in the cell cycle process. Under stress, deoxyribonucleic acid (DNA) replication was accelerated in a dose-dependent manner, yielding subpopulations with different DNA contents [6].

Cell Cycle Model
Flow cytometry data ranging from µ = 0.1 h −1 to 0.6 h −1 for P. putida KT2440 were obtained by Lieder et al. [6] and processed as shown in Figure 1. The data were channeled and displayed as the frequency distribution of DNA content. The durations of cell cycle phases C (DNA replication) and D (period between replication and completed cell division) were determined iteratively by minimizing the deviation between experimental and theoretical DNA histograms. The theoretical DNA content of an asynchronous, ideal culture in which all cells have equal growth parameters was derived from the age distribution according to Cooper and Helmstetter [7]. Using this probability density function for cells of a specific cell age, Cooper and Helmstetter further calculated the theoretical chromosome content per cell at a specific cell age. This model was extended by Skarstad et al. [12] to calculate the frequency of the occurrence of a specific DNA content in an interval of ongoing DNA synthesis. The durations of phases C and D are decisive for the distribution of DNA content. Different values for C were obtained to fit the experimental histograms for various growth rates. Based on the work of Lieder [22], a function for C-phase duration, dependent on the growth rate of P. putida KT2440, was derived. A correlation for C proposed by Keasling et al. [23] was used.
where C is the length of the C period, C min is the minimal length of the C period, µ represents the growth rate and a and b are parameters that fit the experimental data. Based on the experimental data of Lieder et al. [6], the parameter estimation resulted in C min = 0.77 h, a = 1.83, and b = 4.88.

Cell Cycle Model
Flow cytometry data ranging from µ = 0.1 h −1 to 0.6 h −1 for P. putida KT2440 were obtained by Lieder et al. [6] and processed as shown in Figure 1. The data were channeled and displayed as the frequency distribution of DNA content. The durations of cell cycle phases C (DNA replication) and D (period between replication and completed cell division) were determined iteratively by minimizing the deviation between experimental and theoretical DNA histograms. The theoretical DNA content of an asynchronous, ideal culture in which all cells have equal growth parameters was derived from the age distribution according to Cooper and Helmstetter [7]. Using this probability density function for cells of a specific cell age, Cooper and Helmstetter further calculated the theoretical chromosome content per cell at a specific cell age. This model was extended by Skarstad et al. [12] to calculate the frequency of the occurrence of a specific DNA content in an interval of ongoing DNA synthesis. The durations of phases C and D are decisive for the distribution of DNA content. Different values for C were obtained to fit the experimental histograms for various growth rates. Based on the work of Lieder [22], a function for C-phase duration, dependent on the growth rate of P. putida KT2440, was derived. A correlation for C proposed by Keasling et al. [23] was used.
where C is the length of the C period, Cmin is the minimal length of the C period, µ represents the growth rate and a and b are parameters that fit the experimental data. Based on the experimental data of Lieder et al. [6], the parameter estimation resulted in Cmin = 0.77 h, a = 1.83, and b = 4.88.

Geometry and Reactor Setup
In order to generate a pseudostationary glucose gradient of an industrial fed batch fermentation, a large-scale stirred tank bioreactor was chosen. Precise dimensions and information about the inner geometry can be found in Appendix A ( Figure A1 and Table A1). The main geometry was derived from Haringa et al. [21] and only slightly modified for the purpose of this study. With an H/D ratio of 2.57, the total volume was about 54,000 L. The reactor setup included four baffles and a stirrer with two Rushton agitators. The lower stirring unit was equipped with eight blades, and the middle unit with six blades. With a stirring rate of 100 rpm, a tip speed of 5-8 m s −1 was reached. The impeller Reynolds number was 1.8 × 10 6 , the power number 13.15, and the needed power was 226 kW.
The feeding rate was set as half of the maximum uptake rate q s,max of P. putida with 0.738 kg glc ·kg −1 CDW ·h −1 . Aeration, gas transfer, and oxygen uptake were neglected in this study. Therefore, no gassing system was installed. A cell concentration of 10 kg CDW ·m −3 was assumed, and a simple Monod-like kinetic was used to simulate the substrate uptake q s : where q s,max is the maximum uptake rate, c s is the glucose concentration, and the approximated substrate specific uptake constant K s with 10 mg·L −1 . The maximum uptake rate was calculated with the biomass substrate yield Y XS = 0.40 g s ·g −1 CDW and the maximum growth rate µ = 0.59 h −1 [22,24].

Simulation Setup
For the numerical simulation, the commercial calculation tool ANSYS Fluent version 17.0 was used. Using this finite volume-based fluid dynamic analysis program, the virtual geometry was built, and spatial discretization was performed. A total of 445,000 numerical cells yielded the same circulation time as achieved by Haringa et al. [21]. The flow field was approximated by solving the Reynolds-averaged Navier-Stokes (RANS) equations in combination with the standard k-ε model for turbulence. All surfaces were set as slip boundaries, except for the frictionless top area, which implied the reactor filling height. Both impeller units were set to sliding mesh motion to generate a more realistic flow field.
For glucose feed, a separate volume at the top of the reactor was defined, and a constant mass flow was set. The feed was inserted as mass percentage, with constant pressure and volume. The hydrodynamic and kinetic was calculated every 10 ms until the overall glucose concentration was constant and a pseudostationary gradient was reached. Finally, an average flow field and glucose gradient were obtained over 150 s. In further simulations, the hydrodynamic and glucose gradient were set as frozen.
Bacteria lifelines were simulated as massless Lagrangian particles with a discrete random walk (DRW) model passing through the flow field. Every 30 ms, the position and glucose concentration for each bacterium were recorded. In total, 120,000 bacterial cells were tracked over 260 s. According to the ergodic theorem, the same average values are obtained by tracking 1,560,000 bacteria for 20 s (the approximate circulation time). The simulation would yield even more precise statistical evaluations by increasing the number of lifelines.

Statistical Evaluation
All bacterial lifelines were evaluated statistically and grouped according to the regime borders. The growth rate was calculated for each bacterial cell and each time interval. The regimes were classified as follows: standard forked replication S for µ ≤ 0. The second capital letter always indicates the area in which the retention time τ was measured. Before the bacterial lifelines were grouped in regimes, a moving-average filter was applied to filter unrealistic, turbulent fluctuations caused by the standard DRW model (see Appendix B). A second one-dimensional (1D) filter was conducted to erase rapid sequential regime transitions smaller than 0.09 s. Both filtering steps caused deviations from the raw data of less than 5%.
The distribution of the growth rates was derived by calculating the mean growth rate for the whole reactor and the mean growth rate for 20 s for each bacterium. This distribution combined with the cell cycle approach resulted in a distribution of different C-phase durations using Equation (1). Additionally, the energy level distribution was obtained based on Pirt's law [25]: with the Pseudomonas putida properties of nongrowth-associated maintenance m ATP = 3.96 mmol ATP ·g −1 CDW ·h −1 and the growth-associated maintenance Y XATP = 1 85 g CDW ·mmol −1 ATP [24].

Results and Discussion
In order to investigate heterogeneity in large-scale bioreactors, a pseudostationary glucose gradient occurring during fed batch fermentation of P. putida was simulated. Therefore, a biomass of 10 kg·m −3 was assumed, which remained constant within the time observed. For higher biomass concentrations, stronger gradients can be expected.

Gradient and Flow Field
In a 54,000 L stirred tank reactor, a pseudostationary glucose gradient was obtained with CFD simulations. The average glucose concentration was monitored until no further changes could be observed. The residual steady state glucose concentration was 20.7 mg·L −1 . The theoretical growth rate for every numerical cell was computed (Eulerian approach), resulting in an average growth rate of µ = 0.294 h −1 . Ideal mixing was assured by comparing the average growth rate in the reactor (Eulerian approach) and the expected growth rate for the set feed rate µ = 0.295 h −1 . In the fed batch fermentation, the feeding rate amounted to half the maximum uptake rate of P. putida. The objective of the simulation was to generate a realistic glucose gradient with concentrations for which theoretical growth rates ranging from 0.0 h −1 to 0.59 h −1 could be approximated. Moreover, the distribution of bacteria that were introduced from different vertical positions in the reactor at the start of the simulation is displayed.
In Figure 2, three reactor cross sections are depicted to describe (A) the growth rate regimes (see also Section 2.3), (B) the flow field, and (C) the bacterial distribution. Due to asymmetric reactor geometry (see Section 2.2.1), the mean flow field and mean glucose gradient showed periodic changes. Accordingly, the averages of the flow field and gradients over 150 s were computed to track the bacteria ( Figure 2C) as lifelines. Bacteria moved faster when approaching the stirrer. This clearly indicated zones with different residence times. However, tracking the bacterial paths showed that they evenly crossed every part in the reactor.
The underlying gradient was not expected to perfectly reflect the real experiment. Several assumptions had to be made. For simplicity, bubbling flow and oxygen transfer were neglected. The kinetic reaction of substrate consumption following a Monod-like kinetic was assumed to take place in every numerical cell. This implied that the bacterial cells were distributed homogeneously at each time step, which is only a simplified scenario ( Figure 2C). However, to examine the effects of cell history or lag phases of the bacteria on the gradient itself, an existing gradient had to be installed with the stated simplifications. In the following sections, a detailed statistical analysis is provided to study the influence of the gradient on the bacteria and reverse in a realistic manner.
Bioengineering 2017, 4, 27 6 of 12 kinetic reaction of substrate consumption following a Monod-like kinetic was assumed to take place in every numerical cell. This implied that the bacterial cells were distributed homogeneously at each time step, which is only a simplified scenario ( Figure 2C). However, to examine the effects of cell history or lag phases of the bacteria on the gradient itself, an existing gradient had to be installed with the stated simplifications. In the following sections, a detailed statistical analysis is provided to study the influence of the gradient on the bacteria and reverse in a realistic manner.

Lagrangian Trajectory
For 260 s, 120,000 bacteria were tracked on their paths crossing different substrate concentrations. Figure 3 depicts growth rate profiles of two organisms for 20 s, referred to as lifelines L1 and L2. Figure 3C shows the related paths.
According to the regime thresholds (see Section 2.3 and Figure 3A, dashed lines), the growth rate trajectories could be transferred to replication modus curves, as described in Figure 3B). The lifeline L1 revealed high variations in glucose concentrations that were likely to induce strong metabolic changes. In contrast, environmental shifts along L2 were moderate, and there were no effects on metabolism or the cell cycle. The first lifeline L1 gave information regarding five regime transition strategies (STS, TST, STM, TMT, and MTS) and the individual residence times. Lifelines L1 and L2 started from different positions in the reactor and were unequal in length because they moved according to the predominant velocity field. Within 20 s, L2 did not approach the feed zone, remaining in an area of reduced substrate concentration and increased shear stress, owing to the higher velocity of L2.
As shown in Figure 3B,C, within a defined timescale, bacteria completely sensed different environmental conditions. Whereas L2 seemed to remain in the same environment, L1 passed different glucose concentrations and performed several replication strategies. Each metabolic

Lagrangian Trajectory
For 260 s, 120,000 bacteria were tracked on their paths crossing different substrate concentrations. Figure 3 depicts growth rate profiles of two organisms for 20 s, referred to as lifelines L1 and L2. Figure 3C shows the related paths.
According to the regime thresholds (see Section 2.3 and Figure 3A, dashed lines), the growth rate trajectories could be transferred to replication modus curves, as described in Figure 3B). The lifeline L1 revealed high variations in glucose concentrations that were likely to induce strong metabolic changes. In contrast, environmental shifts along L2 were moderate, and there were no effects on metabolism or the cell cycle. The first lifeline L1 gave information regarding five regime transition strategies (STS, TST, STM, TMT, and MTS) and the individual residence times. Lifelines L1 and L2 started from different positions in the reactor and were unequal in length because they moved according to the predominant velocity field. Within 20 s, L2 did not approach the feed zone, remaining in an area of reduced substrate concentration and increased shear stress, owing to the higher velocity of L2.
As shown in Figure 3B,C, within a defined timescale, bacteria completely sensed different environmental conditions. Whereas L2 seemed to remain in the same environment, L1 passed different glucose concentrations and performed several replication strategies. Each metabolic adjustment will cost energy and could have an impact on the production yield.

Regime Transition Frequency
All bacterial lifelines were scanned for regime transitions and retention times in order to obtain the frequency distributions as a function of . Thus, six transition strategies were evaluated in a statistical manner to gain insights into cell histories and possible cell behaviors (see also Section 2.3). Figure 4 shows the counts for each regime transition at a certain retention time. All regime transition statistics, except the TST transition, exhibited a decay after at least 10 s. Bacteria starting from the transition regime T could remain in an area of low concentration for up to 73.5 s (data not shown), where they could grow regularly (standard forked S), before changing back to the T regime. This could be explained by the flow field and gradient pictured in Figure 2A,B. The critical concentrations representing possible growth rates for the regime transition (μ ≥ 0.3 h −1 and μ > 0.4 h -1 ) were located in the upper half of the reactor. Rushton turbines usually cause flow patterns moving away from the blades to the wall, where they circulate up or down, thereby forming large eddies for each stirrer set ( Figure 2B). Consequently, cells will often circulate in this segment and do not pass other areas of the reactor. The lower part of the reactor, which does not provoke a regime transition and, therefore, badly supplies the organisms with substrate, consisted of three segments. As a result, the average retention time in the TST transition was the longest ( ̅ TST = 8.54 s). All other average and maximum retention times are listed in Table 1. The shapes of the distributions follow a Poisson distribution. The maximal retention time was defined as the limit, within which 99% of the values were located.

Regime Transition Frequency
All bacterial lifelines were scanned for regime transitions and retention times in order to obtain the frequency distributions as a function of τ. Thus, six transition strategies were evaluated in a statistical manner to gain insights into cell histories and possible cell behaviors (see also Section 2.3). Figure 4 shows the counts for each regime transition at a certain retention time. All regime transition statistics, except the TST transition, exhibited a decay after at least 10 s. Bacteria starting from the transition regime T could remain in an area of low concentration for up to 73.5 s (data not shown), where they could grow regularly (standard forked S), before changing back to the T regime. This could be explained by the flow field and gradient pictured in Figure 2A,B. The critical concentrations representing possible growth rates for the regime transition (µ ≥ 0.3 h −1 and µ > 0.4 h −1 ) were located in the upper half of the reactor. Rushton turbines usually cause flow patterns moving away from the blades to the wall, where they circulate up or down, thereby forming large eddies for each stirrer set ( Figure 2B). Consequently, cells will often circulate in this segment and do not pass other areas of the reactor. The lower part of the reactor, which does not provoke a regime transition and, therefore, badly supplies the organisms with substrate, consisted of three segments. As a result, the average retention time in the TST transition was the longest (τ TST = 8.54 s). All other average and maximum retention times are listed in Table 1. The shapes of the distributions follow a Poisson distribution. The maximal retention time was defined as the limit, within which 99% of the values were located.  Lifeline statistics provide insights into the frequency of regime transitions and residence times. Depending on the cell history, i.e., the concentrations of bacteria encountered before the bacteria passed the actual concentration, the cells will adapt accordingly. Although metabolic adaptation is known to be very rapid, the initiation of regulatory programs involving transcriptional changes is slower. Investigating the impact of large-scale conditions for E. coli, Löffler et al. [26] showed that fundamental transcriptional programs were initiated after 70 s of glucose shortage. After 30 s, metabolic consequences were measured, and the first transcriptional changes were detected. In total, about 600 genes were found to be up-or downregulated repeatedly, indicating a strong adaption.
Considering this finding during the regime analysis, it is assumed that all cells travelling from high (M) to low (S) substrate availability should be influenced. Being prepared for multifork replication in M, the cells must adapt to standard replication (S). By analogy, this also includes travelers from T to S. Such cells can have a growth rate of about 0.4 h −1 before they adapt to growth rates of less than 0.3 h −1 . During the observation window of 260 s, 72.6% of all cells were expected to carry out this move at least once and to linger more than 30 s in regime S. About 14.7% of all cells were expected to stay more than 70 s in regime S after experiencing higher glucose concentrations in regime T. Furthermore, if a regime transition from maximal to moderate growth conditions (MTS) with the retention time in regime T and S is assumed, 55.5% of all cells performed this move for more than 30 s. A retention time of 70 s was calculated for 10.4% of all cells. The time scales of 30 s and 70 s were shown to significantly influence the transcriptional response of E. coli [26], leading to the  Table 1. Average and maximal retention time in a specific regime. For the six regimes (STS, TST, TMT, MTM, STM, and MTS), the average (τ) and maximal retention times (τ max ) are displayed in seconds. The maximum τ was defined as the limit, within which 99% of the values were located. Lifeline statistics provide insights into the frequency of regime transitions and residence times. Depending on the cell history, i.e., the concentrations of bacteria encountered before the bacteria passed the actual concentration, the cells will adapt accordingly. Although metabolic adaptation is known to be very rapid, the initiation of regulatory programs involving transcriptional changes is slower. Investigating the impact of large-scale conditions for E. coli, Löffler et al. [26] showed that fundamental transcriptional programs were initiated after 70 s of glucose shortage. After 30 s, metabolic consequences were measured, and the first transcriptional changes were detected. In total, about 600 genes were found to be up-or downregulated repeatedly, indicating a strong adaption.

Regime Transition
Considering this finding during the regime analysis, it is assumed that all cells travelling from high (M) to low (S) substrate availability should be influenced. Being prepared for multifork replication in M, the cells must adapt to standard replication (S). By analogy, this also includes travelers from T to S. Such cells can have a growth rate of about 0.4 h −1 before they adapt to growth rates of less than 0.3 h −1 . During the observation window of 260 s, 72.6% of all cells were expected to carry out this move at least once and to linger more than 30 s in regime S. About 14.7% of all cells were expected to stay more than 70 s in regime S after experiencing higher glucose concentrations in regime T. Furthermore, if a regime transition from maximal to moderate growth conditions (MTS) with the retention time in regime T and S is assumed, 55.5% of all cells performed this move for more than 30 s. A retention time of 70 s was calculated for 10.4% of all cells. The time scales of 30 s and 70 s were shown to significantly influence the transcriptional response of E. coli [26], leading to the assumption that changes in adenosine triphosphate (ATP) and guanosine triphosphate (GTP) levels of P. putida KT2440 could also be expected.

Energy and C-Phase Duration Distribution
For the observation window of 260 s, the growth rate profiles of 120,000 bacteria were calculated. Given the set feed rate, the average µ of 0.295 h −1 was expected. Using the Lagrangian approach, an average growth rate of µ = 0.269 h −1 was computed, indicating an adequate deviation of 8.5% compared to the Eulerian approach with µ = 0.294 h −1 (see Section 3.1).
The distribution of the ATP consumption rate q ATP is presented in Figure 5A. The growth rate µ and q ATP were not evenly distributed compared to the mean value, but exhibited individual distributions according to the gradient. The ATP consumption rate was calculated applying Pirt's law (see Equation (3)). While only 6.3% of all cells had a mean ATP consumption rate of q ATP,mean = 29.31 ± 2 mmol ATP ·g −1 CDW ·h −1 , 40.8% showed a reduced consumption rate of less than 27.31 mmol ATP ·g −1 CDW ·h −1 , and 52.9% showed an increased energy demand of 31.31 mmol ATP ·g −1 CDW ·h −1 in comparison to the average consumption rate. Moreover, 12.2% show an energy demand that was more than 1.5 times that of the mean value in the reactor. For the observation window of 260 s, the growth rate profiles of 120,000 bacteria were calculated. Given the set feed rate, the average µ of 0.295 h −1 was expected. Using the Lagrangian approach, an average growth rate of µ = 0.269 h −1 was computed, indicating an adequate deviation of 8.5% compared to the Eulerian approach with µ = 0.294 h −1 (see Section 3.1).
The distribution of the ATP consumption rate qATP is presented in Figure 5A. The growth rate µ and qATP were not evenly distributed compared to the mean value, but exhibited individual distributions according to the gradient. The ATP consumption rate was calculated applying Pirt's law (see Equation (3)). While only 6.3% of all cells had a mean ATP consumption rate of qATP,mean = 29.31 ± 2 mmolATP· g CDW −1 ·h −1 , 40.8% showed a reduced consumption rate of less than 27.31 mmolATP·g CDW −1 ·h −1 , and 52.9% showed an increased energy demand of 31.31 mmolATP·g CDW −1 ·h −1 in comparison to the average consumption rate. Moreover, 12.2% show an energy demand that was more than 1.5 times that of the mean value in the reactor. The distribution will differ if increased nongrowth-associated maintenance m ATP is considered. As outlined by Löffler et al. [26], m ATP increases by 40-50% when cells are exposed to large-scale substrate gradients.
The individual growth profiles of the cells are the basis for deducing cell cycle patterns using the cell cycle model (see Section 2.3). Distributions of the C-length (encoding DNA replication) could be derived for the population of 120,000 bacteria. Figure 5B shows the average duration of replication of 1.21 h and the frequency of cells with a C-phase duration ranging from Cmin = 0.86 h to Cmax = 2.05 h. Clearly, the bacteria were not evenly distributed according to the mean value, and there was a large heterogeneity in the reactor. Although only 22.3% of all cells had a replication phase of 1.21 ± 0.2 h, about 30% possessed a C-period of more than 1.41 h. In contrast, 47.7% displayed a shorter replication phase than the average time for replication (less than 1.01 h). Moreover, approximately 56.1% of the cells were rapidly replicating cells with a growth rate higher than µ = 0.3 h −1 . For these cells, it can be The distribution will differ if increased nongrowth-associated maintenance m ATP is considered. As outlined by Löffler et al. [26], m ATP increases by 40-50% when cells are exposed to large-scale substrate gradients.
The individual growth profiles of the cells are the basis for deducing cell cycle patterns using the cell cycle model (see Section 2.3). Distributions of the C-length (encoding DNA replication) could be derived for the population of 120,000 bacteria. Figure 5B shows the average duration of replication of 1.21 h and the frequency of cells with a C-phase duration ranging from C min = 0.86 h to C max = 2.05 h. Clearly, the bacteria were not evenly distributed according to the mean value, and there was a large heterogeneity in the reactor. Although only 22.3% of all cells had a replication phase of 1.21 ± 0.2 h, about 30% possessed a C-period of more than 1.41 h. In contrast, 47.7% displayed a shorter replication phase than the average time for replication (less than 1.01 h). Moreover, approximately 56.1% of the cells were rapidly replicating cells with a growth rate higher than µ = 0.3 h −1 . For these cells, it can be assumed that they already started to completely adjust their metabolism to achieve multifork replication. As shown in Figure 5B, the bioreactor population was strongly heterogeneous, characterized by a nonequal distribution of bacteria in different cell cycle states. Three different growth phenotypes are shown: C-phase durations of (i) 0.94 ± 0.08 h, (ii) 1.68 ± 0.1 h, and (iii) a transition state of C-phases ranging from 1.1 to 1.5 h. Previously, subpopulations resulting from chemostat experiments have been categorized in populations containing one, two, or more than two chromosomes [27]. With this simulation setup, a model-based superposition of subpopulations containing different growth rates to mimic the scenario in a (fed)batch fermentation was shown. For the underlying gradient, new categories of subpopulations according to the C-phase durations mentioned above can be formulated.

Conclusions
The existence of population heterogeneity in industrial fermenters has been demonstrated, but it still not completely understood. Improvements in fermenter operation, reactor design, and strain engineering can be achieved as more information of cell behaviors during large-scale production becomes available. In this study, the formation of heterogeneity by combining CFD with a cell cycle model of P. putida was investigated. With this method, heterogeneity can be interpreted from the bacterial point of view, particularly with respect to the growth phase durations and energy demands of the cell.
Average and maximum residence times for each transition strategy have been approximated and can be linked to scale-down experiments using STR-PFR setups. Moreover, distributions of growth rates, ATP consumptions, and C-phase durations could be generated. Such findings provide important insights into the intracellular mechanisms that determine growth phenotypes. These mechanisms may become a crucial part of strain and process engineering to predict ab initio and in silico whether and how large-scale performance will meet expectations. Realistic large-scale cultivation can be simulated by investigating the "subpopulations" individually. Specifically, it may be possible to elucidate whether the total drop in production performance during large-scale production is caused by all cells or by individual "subpopulations" that underperform.
To further investigate such problems, heterogeneity studies need to be coupled with single-cell product kinetics. Moreover, research will need to focus on the quantitative measurement of the impact of stress intensity on the m ATP level. This will enable prediction of the total energy demand for a given setup. Table A1. Dimensions of the reactor setup pictured in Figure A1.  Figure A1. Schematic diagram of reactor geometry derived from Haringa et al. [21]. The stirred tank reactor contains four baffles and two Rushton turbines with eight blades (bottom) and six blades (top). Dimensions indicated by capital letters are explained in Table A1.

Appendix B
The standard moving average filter of MATLAB is a linear filter (low pass filter), which removes high frequency components such as fluctuations caused by the DRW model. It is formulated as: m(t) = ∑ y t+j q < t < N − q q j=−q (A1) with: where N is the total number of measured time points and ̅ the filter timescale.  Figure A1. Schematic diagram of reactor geometry derived from Haringa et al. [21]. The stirred tank reactor contains four baffles and two Rushton turbines with eight blades (bottom) and six blades (top). Dimensions indicated by capital letters are explained in Table A1.

Appendix B
The standard moving average filter of MATLAB is a linear filter (low pass filter), which removes high frequency components such as fluctuations caused by the DRW model. It is formulated as: where N is the total number of measured time points and τ the filter timescale.