From industrial fermentor to CFD-guided downscaling
what have we learned?
Haringa, Cees; Mudde, Robert F.; Noorman, Henk J.
DOI
10.1016/j.bej.2018.09.001
Publication date
2018
Document Version
Final published version
Published in
Biochemical Engineering Journal
Citation (APA)
Haringa, C., Mudde, R. F., & Noorman, H. J. (2018). From industrial fermentor to CFD-guided downscaling:
what have we learned? Biochemical Engineering Journal, 140, 57-71.
https://doi.org/10.1016/j.bej.2018.09.001
Important note
To cite this publication, please use the final published version (if applicable).
Please check the document version above.
Copyright
Other than for strictly personal use, it is not permitted to download, forward or distribute the text or part of it, without the consent of the author(s) and/or copyright holder(s), unless the work is under an open content license such as Creative Commons. Takedown policy
Please contact us and provide details if you believe this document breaches copyrights. We will remove access to the work immediately and investigate your claim.
This work is downloaded from Delft University of Technology.
‘You share, we take care!’ – Taverne project
https://www.openaccess.nl/en/you-share-we-take-care
Otherwise as indicated in the copyright section: the publisher
is the copyright holder of this work and the author uses the
Dutch legislation to make this work public.
Contents lists available atScienceDirect
Biochemical Engineering Journal
journal homepage:www.elsevier.com/locate/bej
Regular article
From industrial fermentor to CFD-guided downscaling: what have we
learned?
Cees Haringa
a,b,⁎, Robert F. Mudde
b, Henk J. Noorman
a,caDSM Biotechnology Center, Alexander Fleminglaan 1, 2613 AX Delft, Netherlands
bTransport Phenomena, Chemical Engineering Department, Delft University of Technology, Lorentzweg 1, 2628 CJ Delft, Netherlands cBioprocess Engineering, Biotechnology Department, Delft University of Technology, Lorentzweg 1, 2628 CJ Delft, Netherlands
H I G H L I G H T S
•
A comprehensive overview of the Euler–Lagrange bioreactor simulation approach.•
Application of Euler–Lagrange CFD to three different case studies.•
Different strategies to design scale-down experiments using CFD data are discussed.•
Approach selection chart based on hydrodynamic characteristics of modeled reactor.•
The potential of combining Euler–Lagrange CFD with microfluidics is discussed. A R T I C L E I N F O Keywords: CFD Euler–Lagrange Fermentation Downscaling Metabolic modeling A B S T R A C TEuler–Lagrange computational fluid dynamics simulations offer great potential for the integration of transport dynamics and metabolic dynamics in fermentation systems. Since the seminal work of Lapin et al.[1,2], progress has been made, mainly in the analysis of CFD data and translation to laboratory setup designs. Different large-scale processes require different analysis methods; in this paper we discuss which analysis methods are best suited for given reactor types, by reviewing prior simulation cases as well as introducing new test cases. Fur-thermore, we address challenges in the translation from Euler–Lagrange simulations to laboratory scale systems, and propose methods to work around these shortcomings. Based on the current state of the art, we propose guidelines for the selection of data analysis methods, and we discuss the design of rational scale-down simu-lators. We conclude with a brief discussion regarding the requirements and possibilities of next-generation scale-down simulators, such as microfluidic single-cell analysis, and possible ways to approximate cellular lifelines from invasive intra-cellular measurements.
Many things are known about scale-up. No longer are Rushton impellers the answer. No longer is our concern only in maintaining the same kLa.
Environmental stress due to poor mixing and“hidden” auxotrophy are two factors not fully addressed nor appreciated on scale-up. As a con-sequence, scale-up is still an art not a science.
Arthur Humphrey: Shake Flask to Fermentor: What Have We Learned? (1998)[3]
1. Introduction
In the past 20 years,“from shake flask to fermentor” (scale-up) has been shifting to scale-down: from industrial-scale to lab-scale. The broth in industrial reactors may be heterogeneous in substrate con-centration, dissolved oxygen (DO), shear rate, etcetera. Humphrey, from a scale-up philosophy, hence posed:“we need to learn how to achieve better mixing in large-scale fermentors” [3]. With larger re-actors this is increasingly unachievable; heterogeneity manifests due to the competition between mixing and reaction (expressed in the
https://doi.org/10.1016/j.bej.2018.09.001
Received 30 April 2018; Received in revised form 27 August 2018; Accepted 3 September 2018
Abbreviations: CFD, computational fluid dynamics; CRD, computational reaction dynamics; DO, dissolved oxygen; DRW, dynamic random walk; EL, Euler–Lagrange; FACS, fluorescence activated cell sorting; MAMS, micro array mass spectroscopy; PBM, population balance model; (A)PFR, (axially dispersed) plug flow reactor; RTD, residence time distribution; SD, scale-down; STR, stirred tank reactor; UDM, user defined memory; UDF, user defined function
⁎Corresponding author at: DSM Biotechnology Center, Alexander Fleminglaan 1, 2613 AX Delft, Netherlands.
E-mail address:cees.haringa@dsm.com(C. Haringa).
Available online 08 September 2018
1369-703X/ © 2018 Elsevier B.V. All rights reserved.
Damköhler number Da =τcirc/τrxn), with mixing inherently
scale-de-pendent whereas reactions are microbe-descale-de-pendent. The scale-down philosophy instead asks how broth heterogeneity will affect process performance, taking the effect of heterogeneity into account early in process development. This approach utilizes “scale-down simulators” (SD-simulators), lab-scale reactors that deliberately produce a hetero-geneous environment[4,5]. While developed in the eighties[6–9], SD-simulators gained recent popularity [4] as omics advances enable deeper quantification of the impact of extra-cellular variations on micro-organisms.
A major challenge in downscaling is to devise a lab-scale environ-ment that reflects industrial broth heterogeneity, with only limited knowledge of said heterogeneity. Detailed local substrate/DO mea-surements are generally lacking for existing processes, and even if flow-following probes enable collection of such data in the future, this is of little help during novel process design. Consequently, SD-simulators are often calibrated without industrial reference [10–12], or based on correlations[6,13]. Such experiments provide valuable insight in the response of micro-organisms to extra-cellular variations, but do not necessarily represent large-scale process performance. Computational Fluid Dynamics (CFD) allows to simulate the hydrodynamics of large-scale reactors, and may provide detailed in-silico insight in the microbial environment within fermentors[14,5]. The scale of industrial fermen-tors prohibits full resolution of turbulence and gas–liquid hydro-dynamics, hence approximative models are required[15]. Still, keeping these limitations in mind, CFD offers useful approximation of the large-scale environment, in more detail than offered by correlations or in-dustrial-scale measurements.
CFD is, amongst others, applied to study mixing [16,17], mass transfer[18,19]and substrate gradients[14]. While such studies can supply input for downscaling, the aim of SD-simulators is to capture temporal environmental fluctuations as observed by micro-organisms. Euler–Lagrange CFD simulations, in which a large number of virtual particles (parcels) is tracked, are especially suitable: they allow to track extra-cellular environment of each parcel in time (called “lifelines”)
[1,2]. In our view, these lifelines form a preferential basis for SD-design
as they provide detailed dynamic information that is unavailable from purely Eulerian simulations. In a Sino-Dutch collaboration, we explored the use of Euler–Lagrange CFD for SD-design[20]; here we provide an overview of the lessons learned during this project. These lessons in-clude Euler–Lagrange reaction coupling, lifeline analysis methods and their relation to reactor configuration, and scale-down design as well limitations therein. To conclude, we outline some considerations for future scale-down simulators, and present a conceptual design of a micro-fluidic scale-down simulator that may overcome limitations ob-served for bench scale systems.
2. Bioreactor hydrodynamics
First, the Eulerian model– resolving broth (and air) flow – must be developed. The following aspects are to be considered:
•
Macromixing (impeller & gas/liquid interaction)•
Mesomixing (feed mixing)•
Micromixing (e.g.film & intra-pellet mass transfer)•
Aeration (phase interaction, bubble size, etc.)•
Solid-liquid interaction (turbulent motion, settling, etc.)•
Gas–liquid mass transfer (kla models)•
Heat transfer•
Rheology (& turbulence/aeration interaction)•
Liquid-microbe transfer/reactions (uptake/excretion)•
TurbulenceThe Euler–Lagrange approach is compatible regardless of the choices made regarding the above. Inclusion of all aspects leads to demanding models, while data for modeling and validation may be scarce. Often, one or more aspects are negligible or absent; a timescale/ force balance analysis then allows to discard irrelevant processes. Macro-mixing[14], gas–liquid interaction[21,22]and microbial reac-tions are often considered dominant[14]. Meso-mixing can be relevant near the feed inlet point, where concentrated glucose blobs may be segregated from the broth (hence not observed by biomass)[14,23,24]. Nomenclature
A area (general), m2
Cs substrate concentration (mole based), mol/kg
Cext generic extracellular concentration, mol/L
Cint generic intracellular concentration, mol/gx
cs substrate concentration (mass based), mg/L
Cx Biomass concentration, g/kg
CD drag coefficient, –
Dr dilution rate, h−1
m diffusion coefficient, m2/s
t Turbulent diffusion coefficient, m2/s
F feed rate (general), kg/s
Fs Substrate feed rate (general), mol/s
kLa overall mass transfer coeff., h−1
Ks affinity constant for substrate, mol/kg
Ns agitation rate, s−1
Nc total number grid cells,–
Np total number particles,–
Np,c Number particles in cell c,–
qp specific production rate of product, molp/Cmolx/h
qs specific uptake rate of substrate, mols/Cmolx/s
qs,max max. specific uptake rate of substrate, mols/Cmolx/s
qref reference qs/qs,max,–
Rp reaction rate of s, parcel-based, m
Ss source term of s, mol/kg/s
Sij average rate-of-strain (component), s−1
t time (general), s
Ul liquid velocity, m/s
Ug superficial gas velocity, m/s
V tank volume, m3
Vc cell volume, m3
Vp parcel associated volume (VT/NP), m3
VT total volume (general), m3
α gas fraction,–
αg geometry constant for no. parcels,–
γ surface tension, N/m
γ shear rate, s−1
Δt timestep size, s
ϵ turbulent energy dissipation rate, m3/s2 θ95 dimensionless mixing time,–
λ inter-phase mass imbalance (Euler–Lag.), %
μ specific growth rate, h−1
μl effective Viscosity, dynamic, Pa s
μl,eff viscosity, dynamic, Pa s
μt turbulent viscosity, dynamic, Pa s
ρ density, kg/m3
τ95 mixing time, s
τarc arc-time, s
τcirc circulation timescale, s
τrxn uptake timescale of substrate, s
ϕi circulationflowrate, kg/s
χ heterogeneity in domain,–
Based on the analysis of Linkes et al., eddy micro-mixing andfilm dif-fusion are estimated to be rarely limiting[23], although this may differ
in highly viscous processes. These aspects may be included via subgrid concentration distribution and film diffusion models, respectively. A concern is that potential mass transfer limitations at the cell surface are rarely checked in experiments, which may compromise the determi-nation of kinetic parameters[25,24]. Even if rarely an issue in practice, it is wise to adopt practices from catalysis engineering to rigorously exclude- or account for them.
Non-Newtonian rheology may be highly relevant [26,27], but is often omitted due to modeling difficulties. In highly turbulent flows, assuming an overall effective Newtonian viscosity μl,effmay be allowed,
but in transitional flows effective viscosities may differ locally. Moi-lanen et al. found large stagnant caverns in simulations with a Her-schel–Bulkley fluid[27], yielding simulated mixing times in excess of 10 s, compared to ca. 3 s experimentally. While cavern formation may appear in yield-stressfluids, the stagnation observed in these simula-tions does not appear to reflect reality. They based the local μl,effon the
meanflow-shear rate γ =|∂Ui/∂xj|; it can be argued that the shear rate
based on energy dissipation,γ= ϵ /ρ μl,eff, may be more representative. In exploratory simulations we found reasonable single-phase results with this approach, but multiphase simulations diverged (Appendix A). For single cells, the Stokes number St < < 1 meaning they can be considered massless flow-followers. For cell agglomerates, however, settling effects may become relevant. Turbulent motion of parcels has to be considered. The default way to describe this is via the discrete random walk model (DRW), which may suffer from artificial accumu-lation of parcels in low-turbulence regions, e.g. near walls [2,28]. Continuous random walk (CRW) models[28]or a Fokker–Planck type
treatment[2], may yield superior results. Practically, we observed little accumulation in stirred tanks due to the low surface/bulk gridcell ratio
[29]. Possible intra-pellet diffusion limitations can be evaluated via reaction-diffusion models, although determining the intra-particle dif-fusion coefficient may be challenging. The computational burden makes large eddy simulations unsuitable for (routine) application in bioprocess simulation; cruder Reynolds-Averaged Navier Stokes (RANS) simulations have to be relied on. In reactors lacking significant flow transients (including stirred tanks with the MRF-impeller model), a steady-state solution may be a reasonable approximation, reducing computation time[30]. In the end, the choice of models is dictated by the relevant physics, the required level of detail, and the availability of verification data.
3. Reaction coupling 3.1. Coupling methods
Reactions can be coupled to CFD simulations via regular Eulerian models, population balance models (PBMs), and Euler–Lagrange (EL) models (which are our primary interest). The current focus is on me-tabolic models [1,2,31,32], and metabolic models including key en-zymes [33]. Similar strategies can be applied to include a broader spectrum of cellular responses (transcription and translation dynamics
[34–36], protein formation[37], cell cycles[38], and more), provided the resolved timescales are sufficient to capture such dynamics.
Regular and PBM coupling. Regular reaction models compute rates based on the local extra-cellular conditions. This assumes instantaneous equilibrium between intra- and extracellular conditions, which is valid for uptake processes lacking metabolic feedback (for example lacking intra-cellular glucose inhibition) and with constant enzyme capacity, which is justified when considering “process snapshots” of a few mixing times. Directly linking long-term intra-cellular (enzymatic) adaptation to local extra-cellular conditions may lead to gross errors; for example coupling a penicillin production model (enzyme adaptation time τad≈ 20 h) to a CFD simulation led to an improbable yield loss of 85% [39]. PBMs account for local non-equilibrium and parcel history by
incorporating a transportable distribution of the intra-cellular compo-sition, typically with growth rateμ as a representative parameter. Im-proved acetate overflow predictions for E. coli were found with the approach [31], and and simulations of the full process duration are feasible if (pseudo-) steady hydrodynamics can be assumed [40,41]. Capturing the distribution by traditional class-based methods is tract-able for narrow composition ranges (bubble size distributions for ex-ample[21]), but becomes cumbersome for broad distributions. Moment methods offer an alternative. Initially, a comparison between 100 and 400 classes and several moment methods found mostly gains in memory use[42], more recently also time-gains were observed[43]. Most stu-dies assume one-dimensional heterogeneity, captured byμ; this may be insufficient in some cases. Consider a model with enzyme X (adaptation time τX≈ 1 h ≫ τcirc) subject to metabolic control by M
(τM≈ 30 s ≈ τcirc). This requires accounting for heterogeneity at two
timescales; X may be distributed in the population, but the long turn-over time makes the distribution spatially homogeneous. The distribu-tion in M will vary locally, making X and M uncorrelated. Yn class-interactions may occur for n pools with Y classes, posing a potential computational challenge. Models with two compositional variables have been reported [44]. How the PBM framework performs with higher-dimensional metabolic models is a development to be followed with interest.
Lagrangian coupling. In EL methods Npindividual parcels are tracked.
EL simulations come in two forms. EL simulations with Eulerian reac-tion coupling solely employ the parcels to probe extra-cellular condi-tions over time, registering what organisms “see” [32,39,22,38]. However, this approach allows for only limited interaction between the extra- and intra-cellular domain[39,33](see Section 3.2). Reactions can also be coupled to the Lagrangian phase. In this case, the cellular composition is tracked by assigning a composition vector Xn with n
pools to each parcel. Each pool is quantified by a single value, making pool interactions straightforward. A metabolic model, taking the extra-cellular concentrationfield(s) of parcel p as a boundary condition, de-scribes pool and environmental interactions. In contrast to PBM, there are no pseudo-steady states for the Lagrangian phase, requiring a transient solution with shortΔt to capture parcel motion and pool dy-namics[29]. Typically O(> 105) parcels are needed for a smooth
bio-mass distribution[29], making full fermentation simulations infeasible with full Lagrangian coupling. These requirements make EL most sui-table to simulate multi-dimensional heterogeneity on minute-hour timescales; the unique parcel perspective provides “lifelines” of the observed environment in time. We do emphasize these models re-present an average response to extra-cellular variations. At the single-cell level, intrinsic heterogeneity[45]andfluctuations in transcription levels[46]will induce variations in this response, that are not captured by typical metabolic models.
Computational advances are required to enable full fermentation and (faster than) real-time simulations. With innovations in highly parallel and GPU-accelerated CFD giving rise to applications like ANSYS Discovery Live, the prospects are positive, underlining our preference for the EL method. In the next section, practical coupling implementations are discussed.
3.2. Lagrangian coupling: practical aspects
Lagrangian reaction coupling comes with practical challenges, for example regarding Np. Whether or not full Lagrangian coupling is
re-quired depends on how the inter-phase (uptake) reactions are con-trolled by the cellular composition, where we distinguish three situa-tions:
•
C1: no rapid or slow intracellular control.•
C2: only slow intracellular control.The cases are graphically summarized inFig. 1. In C1, uptake re-actions are completely uncoupled from the Lagrangian phase, meaning straightforward Eulerian reaction coupling can be employed and par-cels are only used to register extra-cellularfluctuations[32,39,22,38]. The number of parcels Npneeds to be sufficient to converge fluctuations
statistics (typically 100–104, depending on resolvedflow-time), but is
otherwise unrestricted. The assumption of no intra-cellular control (e.g. simple Monod kinetics with constant parameters) is valid if there are no rapid turnover intra-cellular pools (τpool<τmix) controlling the rate
(or: instantaneous equilibrium can be assumed), while changes in slow pools (τpool≫ τmix) are insignificant (growth and enzyme dynamics can
be neglected in the chosen timeframe). Intra-cellular responses can be calculated in the CFD software or during post-processing, as there is no feedback to the Eulerian phase.
In C2, growth and/or enzyme level changes have to be accounted for, and the population may exhibit significant heterogeneity in enzyme levels. However, sinceτpool≫ τmix, the population distribution will be
the same everywhere in the domain. This allows a hybrid
Euler–Lagrange model: uptake reactions are coupled to the Eulerian phase, but reaction parameters like qs,maxare determined based on on
population average enzyme levels[33]in the Lagrangian phase. In this case Npmust suffice to capture enzyme pool distributions. In C3, local
uptake/excretion rates are affected by intra-cellular heterogeneity in rapid turnover pools, which requires full Lagrangian coupling to cap-ture[1,2]. In this case, the distribution of parcels determines the dis-tribution of biomass, and Np must be sufficient to produce a (near)
homogeneous biomass distribution, in order to ensure realistic con-centrationfields. This puts a criterion on Np[29]that is treated in the
next section.
3.2.1. C3: interphase coupling
In Lagrangian coupling, inter-phase species transfer is key; uptake by parcels must equal field depletion. The pseudocode inFig. 2 de-scribes an (explicit) implementation of parcel-bound kinetics in ANSYS FLUENT. The workflow in FLUENT means DPM scalars, and hence source terms Ss, are updated after the field iterations: what parcels
consume in step t, is taken from the extra-cellular domain in step t+1. Aside from challenges with explicit integration, an issue is that the mass balance is violated if uptake exceeds availability (Ss(c, t)Δt >
Vc(c)Cs(c))[29]; source term‘clipping’ will occur. As parcels do not
“see” each other, this may also occur when multiple parcels reside in the same gridcell, competing for the same substrate. Re-calculating the exchange terms in each iteration accounts for the presence of multiple parcels (uptake now fully takes place in timestep t), but increases computational expense, and it was not found to yield improved accu-racy[29]. In fact computing qs once per timestep performed slightly
better, possibly due to exploiting sub-flowtime time-stepping for parcel updating. In open-source or dedicated codes, an implementation that combines sub-flowtime stepping with implicit integration may be more
Fig. 1. Coupling approaches used in this work. Left: C1: instantaneous intra-cellular response or constant intra-cellular environment, leading to a homogeneous population. Middle: C2: slow intra-cellular variations, leading to a heterogeneous population, but without spatial variation. Right: C3: rapid intra-cellular variations, leading to a population distribution that is also spatially heterogeneous. The shown distributions are hypothetical distributions on the single-cell level. In practice, the resolved number of parcels in current Euler–Lagrange simulations does not suffice to locally reconstruct such instantaneous distributions of cellular composition.
f
f f
Fig. 2. Pseudocode for an (explicit) implementation of parcel-bound kinetics in ANSYS FLUENT. Updating the reaction-rates and intracellular pools is done using a DPM scalar UDF. Uptake rates are transferred to used defined memory (UDM) and included in Eulerian computations in the next timestep via source term UDFs. An execute at end UDF is used for data storage, and to reset the UDM.
straightforward to implement than in commercial codes.
3.2.2. C3: simulation criteria
To resolveflow dynamics for many parcels requires a proper bal-ance between computational burden and accuracy; Npand timestepΔt
have a large influence. It is clear from the pseudocode that if Np< < Nc, Ss= 0 in many gridcells, i.e. no biomass will be present,
resulting in unphysical qsand Csgradients. The volume of gridcells in a
locally refined grid may range over orders of magnitude, meaning some cells contain many parcels, while others are mostly empty even if Np≫ Nc; ensuring homogeneity by increasing NPis unrealistic. An
al-ternative, distributing qsin a parcels‘effective volume’ Vp= Vdomain/Np
requires cumbersome listing of surrounding gridcells. However, if tur-bulent distribution of substrate is sufficiently fast, artificial gradients will be minor even with inhomogeneous parcel distributions. By linking the magnitude of unphysical gradients to the ratio of timescales of re-action and turbulent dispersion, and includingfluctuations in parcels-per-cell Np,c, a criterion for Np can be developed. Substituting the
timescales with the underlying hydrodynamic and kinetic parameters this gives (for Monod and 1st order kinetics)[29]:
∑
⎜ ⎜ ⎟⎟ = ⎡ ⎣ ⎢ ⎢ ⎛ ⎝ ⎜ ⎛ ⎝ − ⎛ ⎝ ⎞ ⎠ ⎞ ⎠ ⎞ ⎠ ⎟ ⎤ ⎦ ⎥ ⎥ N α q C K χ V V V V 1 p g s x s T N c t c c T 2 ,max 7/6 , 1/2 2 c (1) Withχ an accuracy parameter and α a geometry-dependent con-stant. For typical conditions, this gives Np= 104–106ifχ = 0.05. Withthis, running full fermentation simulations (O(> 10 h) on desktop computers is infeasible, but several mixing times (O(min–h)) to gain insight in metabolicfluctuations are manageable. A timestep size of min (1/(10Ns), τrxn/10) is advised for stirred tanks, with Ns the stirring
speed in s−1andτrxnthe timescale of the fastest exchange reaction[29].
Thefirst of these criteria considers flow trajectory calculation and may differ per reactor type. The pioneering studies by Lapin et al. did not explicitly note the effect of parcel distribution on uptake; with 105
parcels in 443gridcells[1,2], their work was likely in accordance with these criteria. All criteria are based on 1st order accurate explicit in-tegration for the exchange terms; the presence of stiff intracellular re-actions may change the criteria and favor the use of higher order ex-plicit or imex-plicit methods, at least for the intra-cellular equations.
3.2.3. Model structure
Dynamic kinetic models should be designed with CFD in mind, meaning they should be stable towards turbulentfluctuations in extra-cellular concentrations Cs. Using a model for penicillin production[47],
two potential issues were identified: (1) noisy rates when extra-cellular control was assumed over intra-cellular rates and (2) non-zero rates at zero pool size[33]. Concerning thefirst, storage of glycolitic inter-mediates to polymers was non-linearly controlled using extra-cellular substrate as a signal molecule, strongly affecting rates in the range of turbulent fluctuations. Rates hence oscillated rapidly, inducing in-sufficient release of storage carbohydrates during starvation, which reduced ATP availability. Such issues (and numerical issues noted above) are more likely present in systems with low Ks, as small Cs
variations may lead to strong qsvariations. Low XATPquenches
depo-lymerization, while the ATP-consuming storage process continued (a non-zero storage rate at low XATPwas no realistic scenario in model
development), eventually driving XATPnegative.
Numerical instability was excluded as a cause by varying integrators (implicit, explicit, RK4), timestep sizes (Δt = 0.0005–0.5 s), and en-forcing noise in MATLAB. The issues related to the non-linear differ-ential-algebraic formulation of the model, not to the solution algorithm. The problem was bypassed by correlating XATPwith the glycolitic pool
XGLY, but storage rate oscillations remained[33]. More generally, direct
extra-cellular rate control over intra-cellular rates should be avoided (instead, use intra-cellular substrate, ATP or redox factors as a signal), or a signaling timescale should otherwise be introduced, for example using an “activated enzyme pool”, to buffer rapid extra-cellular dy-namics. Furthermore, ensure all rates quench when their reactants de-plete, even if depletion is unlikely. Most dynamic metabolic models described in literature are based on ideal-mixed systems with smooth dynamics, not anticipating CFD application. For example, Lei et al.[48]
eliminated all intra-cellular pools for S. cerevisae, as intra-cellular changes could be considered instantaneous compared to slow (fed-) batch dynamics. Imposing rapid noise on this model led to similar issues as described above. We did not conduct rigorous stability analysis, which may be a future avenue to formalize model requirements.
3.2.4. Compartment-based Lagrangian simulations
The high Nprequired to simulate the microbial phase induces a large
computational burden. With C3 coupling and unsteady hydrodynamics,
Fig. 3. Lagrangian analysis methods used in this work. Left: Regime analysis, based on regions in qs-space with a consistent metabolic response. Most suitable for
compartmentalizedflows (often found in stirred tanks), and for application with multi-compartment downscaling. Middle: Arc-analysis, based on the duration between subsequent crossings of a reference value in qref. Most suitable for non-compartmentalizedflows (smooth extra-cellular variations), with a wide residence
time distribution (some stirred tanks, bubble columns) (Right: Fourier analysis, seeking the dominantfluctuation frequencies. Most suitable for flows with narrow residence time distributions, such as airlift-loop or pipe reactors. Here, methods were illustrated with qsvariations, but they are broadly applicable.
simulating several mixing times can take weeks on a desktop computer; in contrast, a full fermentation period is feasible in 1–2 weeks assuming steady state flow and C1–C2 reaction coupling. Using the
hydro-dynamics of compartment models using CFD, and modeling
(Lagrangian) reactions in these compartment models, can strongly re-duce computation time by lowering both Npand grid size Nc, without
severely compromising accuracy [49,50]; parcel tracking in such a model can be done, for example, using a statistical (Markov-chain) rather than a trajectory-resolved approach[51–54].
4. Lifeline analysis methods
Lifelines represent cellular data as a function of time, and can be subjected to signal/statistical analysis to report extra-cellular fluctua-tions experienced by organisms. One of the main goals of lifeline ana-lysis is to provide guidance in the design of scale-down experiments. The methods described below aim to capture statistics describing the duration and magnitude of extra-cellularfluctuations, which can then be used to set design parameters for scale-down simulators (as de-scribed in Section5). To our knowledge, such a direct CFD-based scale-down design has not yet been used, at least in public literature. There are, however, examples where scale-down design was aided by CFD
results, for example in setting the order of magnitude in fluctuation frequency, shear-rate orfluctuation amplitude[55–58].
Three methods were addressed: Fourier analysis, Arc-analysis and Regime-analysis. A graphical overview of methods is provided inFig. 3. The preferred method depends on:
•
Flow behavior (plug-like or well circulating).•
Compartment formation.•
Operating conditions.•
Application (downscaling method).Low-passfiltering of lifelines to remove rapid turbulent variations is advised, as these are unlikely to have metabolic influence and may affect the statistics of large scale fluctuations. Smoothing using the timescale of turbulentfluctuations as a window is the most straight-forward approach[39].
Fourier analysis. The Fourier method is fast and straightforward, giving most insight with periodicfluctuations, e.g. narrow circulation time distributions. For broad distributions, data interpretation is chal-lenging[39]. First, the mean is subtracted from the lifeline (as a re-ference value) which is then multiplied with a window function and fed to an FFT algorithm. The real parts of the transforms are subsequently
Fig. 4. Analysis selection diagram. Top: some considerations for SD-selection. Gradient-induced population heterogeneity means different cells have significantly different lifelines (In single-vessel scale-down simulators, some difference between individual cell experiences may exist by micro-mixing, but all cells experience the same explicitly induced variations). Operational challenges are linked to pumping issues for shear-sensitive and complex rheology systems. The selection chart gives the preferred analysis method; other analysis methods may still be applicable.
summed to a composite spectrum to reveal dominant frequencies. Amplitudes are retrieved by visual inspection of the lifelines; extracting them from spectra is more challenging.
Arc-analysis. Arc-analysis records fluctuations compared to a re-ference qref(e.g. the mean) in the time-domain[39,22,33]. The interval
between consecutive qrefcrossings gives thefluctuation duration, akin
to a circulation time measurement in reactant rather than physical space. The magnitude follows from the extreme lifeline value between crossings; correlating magnitude and duration provides insight in ty-pical circulation trajectories. The three points (start, extreme and end) form an“arc” describing the fluctuation. The time at which the max-imum is recorded can be used to study iffluctuation trajectories are symmetric. To remove negligible oscillations, it is advised to only register events that exceeding qrefby a certain threshold. A drawback
arc-analysis is that secondary motions are not captured; the method is most applicable for relatively simple circulation behavior, not for strongly compartmentalized reactors.
Regime analysis. In regime analysis, the domain is divided into zones based on a characteristic quantity such as a certain metabolic response (e.g. Cs> Cs,overflow), or a relatively homogeneous concentration (a
compartment around an impeller) [39,22]. The fluctuation duration follows from the duration of individual regime visits. Again, it's advised to set a crossing threshold (a“fuzzy boundary”[29]) to remove minor fluctuations. The amplitude is replaced by a per-regime mean value; if oscillations around this mean are relevant, arc-analysis can be con-ducted within a regime. Regime analysis works well with more complex systems, like compartment-forming multi-impeller reactors.
The role of downscaling. Fourier and arc-analysis are best suited for single-vessel, fluctuating input designs. Both provide the duration (distribution) of observedfluctuations, and the amplitude compared to a reference value. Regime analysis directly links to multi-vessel scale-down, each vessel representing a regime, operating at the regime mean. This means that the number of regimes may be practically constrained
[22]. For a combination of feedfluctuations and multiple vessels the combination of regime and arc-analysis seems most applicable. The preferred downscaling method follows from lifeline structure (smooth versus discrete jumps) and practical considerations. In a single vessel SD-simulator, all organisms see the same; the effect of rare, extreme fluctuations cannot be captured. A multi-vessel simulator accounts for the duration distribution by design, but it does not account for the
Fig. 5. Analysis of Case I. A: Example lifeline, the red and blue dashed line providing regime boundaries. B: Residence time distributions (RTDs) per regime (solid lines),fitted ideal stirred tanks RTDs in a 3-vessel SD-simulator (dashed). The volumetric regime distribution is reported inFig. 3A; the inset shows conceptual lifeline discretization. C: 3-vessel SD-simulator based on the regime division, RTDs, and regime-mean operating conditions. D: Example of a SD-simulator lifeline. Image reproduced from[22]with permission.
amplitude distribution; this is better captured influctuating input sys-tems. Single-vessels are also more straightforward to operate and ana-lyze. Sweere et al.[7–9]and Wang et al.[55]are among the few to compare two downscaling approaches. Both observed significantly different responses at equal mean fluctuation time, including hints at different substrate transporter expression, showing that an increased understanding of the metabolic response to extra-cellular variations is still required. InFig. 4we present a selection chart for the preferred analysis method depending on the problem statement and downscaling approach, built on our experience to date.
5. Case studies
Three different case-studies were explored in this project, with different lifeline structures and analysis methods.
5.1. Case I: stirred yeast fermentation
A top-fed stirred, aerated yeast fermentation in a 22 m3 research fermentor was simulated[14,59–61,56]. Details regarding computation and validation are found in[22]. The analysis considers a short time-frame (ca. 1800 s). Monod kinetics for glucose consumption were as-sumed, using C1-coupling; dissolved oxygen and ethanol were not si-mulated. The observed concentration gradient is in decent agreement with experimental data[14]. Eight circulation loops form around the 4 impellers, yielding strongly compartmentalizedflow; regime-analysis is the preferred analysis method. The reactor circulation timeτcirc≈ 40 s
is close to the characteristic reaction time Ks/Cxqs,max= 38 s which
suggests a global gradient. The per-loop circulation time isτloop≈ 40/
8≈5 s, making the individual compartments well mixed.
7 qsranges can be discerned[22]. While a 7-regime analysis can be
conducted, the qsjump between compartments is small, and the average
residence time impractically short. A 3-regime division is therefore selected, in line with state of the art downscaling[62]. The boundary between regime 1 and 2 is set to qs/qs,max= 0.05 to capture the large,
nearly homogeneous bottom region. The regime 2–3 boundary is set to
qs/qs,max= 0.2, the expected transition to Crabtree ethanol
fermenta-tion[63]. The regime layout is provided in Fig. 3A, the lifelines in
Fig. 5A. From the acquired regime residence time distributions (RTDs) (Fig. 5B), it is possible to design a multi-vessel SD-simulator using the 5 degrees of freedom posed by Noorman[64]:
•
Number of regimes→ number of reactors.•
Per-regime volume fraction→ Reactor volume ratio.•
RTD→ Reactor circulation pattern.•
Mean RTD, regime volume→ Exchange flowrates.•
Mean regime conditions→ Feed/aeration rates.Despite the complex flow patterns, each regime RTD shows the exponential decay of a recirculating flow (Fig. 5B); 3 well-mixed compartments are hence selected. The exchangeflowrates are found by fitting the regime RTDs with perfect stirred-vessel RTDs (Ignoring the short term peaks, and under the restriction that∑(ϕcirc) = 0). Knowing
theflowrates, MATLAB optimization is used to specify feed and drain rates, leading to the design inFig. 5C, with an example lifeline depicted in Fig. 5D. In order to set the operating qs, the SD-system needs to
operate at or above Cx= 15 g/kg; the underlying reason is further
discussed in[39,22].
5.1.1. Case study II: airlift yeast fermentation
The forcedflow structure in an airlift loop reactor leads to a more narrow RTD[66,67], which is revealed in the lifelines. Lacking an in-dustrial case, we simulated a hypothetical fermentation based on the geometry and conditions of Simcik et al.[65]. A bubble size of 5 mm was assumed. The realizable k− ϵ turbulence model, Universal drag model and Simonin et al. turbulent dispersion model were used[68]. 2nd upwind discretization was used except for volume fraction, for which 1st order upwind was used to ensure stability. Validation results are inTable 1. We employed the yeast kinetics used in case I with a biomass concentration of 30 gdw/kg. The liquid phase properties
equaled water, for air ρ = 1.2 g/L, μl= 18.6 × 10−6Pa s were set.
Profiles of steady-state holdup and liquid velocity are shown inFig. 6A, B.
The feed Fs= 0.00194 mol/s, included as a source term at the
bottom of the draft tube, at the high-qsspot inFig. 6C. The feed was set
such that under ideal mixing conditions Ks= Cs. The non-ideal
simu-lated conditions result in a gradient ranging from 0.45 < qs/
qs,max< 0.8. 7500 particles were tracked for 200 s.Fig. 7shows a
ty-pical track, showing periodicity, in contrast to those from stirred ves-sels. Still there is duration variability, from two sources: 1) dispersion in the riser/down-comer and 2) recirculation in the headspace. Two ‘types’ of peaks are distinguished in the lifelines; peaks with qs/
Table 1
Validation of airlift loop simulations, compared with experimental data of Simcik et al.[65].
Parameter Simulation Experimental ul(m/s) (draft tube) 0.7 0.75 ul(m/s) (downcomer) 0.18 0.18 α (draft tube) 0.06 0.07
α (downcomer) 0 0.02
Fig. 6. Eulerian profiles for the Airlift loop reactor. A: Gas holdup profile. B: Liquid velocity profile. Vectors show direction, colors show magnitude. C: Uptake rate profile. Substrate is fed right above the sparger, leading to a high-uptake hotspot at this location.
qs,max> 0.6 from feed point passages, peaks with 0.52 < qs/
qs,max< 0.6 from down-comer passes and headspace recirculation; a
high peak followed by a low peak is a direct full circulation, a high followed by multiple lows signifies repeated headspace recirculation. The riser peaks are uniform, reflecting the uniform velocity profile; the recirculation/downcomer peaks show more variability. The periodic nature makes the lifelines more suitable for Fourier analysis than prior examples. The combined spectrum of 7500 lifelines (Fig. 7) shows two distributed peaks centered around 0.05 Hz and 0.1–0.12 Hz, the first reflecting the interval between feed-passings, the second reflecting the secondary circulations.
Three straightforward scale-down strategies can be devised: (1) a variable feed reactor, total cycle time of 20 s (0.05 Hz), imposing qs/
qs,max≈ 0.65 at the cycle start, optionally using a secondary injection at
t = 10 s to reach qs/qs,max≈ 0.54. (2) A STR-PFR combination with the
STR representing the headspace/downcomer and the PFR as the riser, containing the feed point (Fig. 8A). (3) A 2-PFR combination, one near-perfect PFR for the riser and a more dispersed PFR for the headspace/ downcomer (Fig. 8B). Putting the “riser” above the “downcomer” compartment has the benefit of allowing selective aeration in the “riser” section, if combined glucose and oxygen gradients are desired. We present conceptual designs for systems (2) and (3), byfitting the qs/
qs,maxdistribution between the CFD results and idealized systems using
the MATLAB function fmincon, assuming the system in chemostat op-eration with a dilution rate Dr= 0.05 h−1and Cx= 30 g/kg. Thefitted
parameter values are reported in Fig. 8A, B, for concept 2 and 3,
Fig. 7. A: Single lifeline in the airlift-loop reactor. B: Composite Fourier spectrum of 7500 tracks, peaking at 0.05 and 0.11 Hz, associated with reactor and headspace circulations respectively.
Fig. 8. Conceptual design of multi-vessel SD-simulators for the airlift case. A, B: Sketches of the PFT-CSTR and PFR-PFR con-figuration, respectively, with design parameters as determined byfitting the qs/qs,maxdistributions. Note the desired axial
dis-persion depends on the reactor cross-sectional area A C: Frequency spectra for the lifelines in both designs, compared with the CFD result. D: qs/qs,maxdistribution for both designs,
compared with the CFD result. E: Examples of lifelines for the scale down simulator designs proposed in A, B, generated using idealized CSTR/axially-dispersed PFR models in MATLAB. These lifelines are designed to replicate the airlift-loop lifelines re-ported inFig. 7A.
respectively.Fig. 8C–E shows the Fourier spectrum, qs/qs,max
distribu-tion and example lifelines, respectively. Note the Fourier spectrum was not fitted – by requiring a similar qs/qs,max distribution a decent
fre-quency spectrum match emerged automatically. Due to rapid dilution near the feed point in CFD simulations, the peak in qs/qs,maxcannot be
reproduced explicitly.
5.1.2. Case study III: stirred penicillin fermentation
A broad study with both C1- and C2-coupling was conducted in a 54 m3 stirred fermentor (described in [39,33]), using the metabolic model by Tang et al.[47]with adaptations to ensure stability[33]. The effect of viscosity was omitted due to divergence. The mixing time over-estimated experimental data by about 25% with aeration included, excluding aeration gave a similar under-estimation. Withτrxn= 0.32 s
and τcirc≈ 26 s, a strong gradient appeared, which manifests itself
completely in the top impeller loop, with the vessel bottom operating under starvation conditions. For end-of-fermentation snapshot simula-tions (Cx= 55 g/L, C1-coupling), a penicillin yield Ysploss up to 40% is
found, depending on the circulation time. A linear relation between Ysp
and Damköler numberτcirc/τrxnis observed in the range, allowing for
interpolation. The single-loop gradient gives smoothfluctuations with a wide duration distribution (Fig. 9A), suitable for arc-analysis (Fig. 9B,C).
Setting qs/qs,max= 0.05 as qref, only arcs above qref have to be
analyzed in magnitude, showing a clear relation between duration and magnitude (Fig. 9D,E). Below qref assuming qs→ 0 and registering
duration suffices. From this data, an SD-simulator with variable inter-vals can be designed. Using inverse transform sampling, a string of random intervals can be generated that abides the duration pdfs (Fig. 9D), alternating between positive and negative (Fig. 9F,G). For positive arcs, the amplitude is set from the duration-magnitude relation, computing the feed from the mass balance. In-silico assessment of the SD-system assuming ideal mixing shows good agreement in yield loss compared to large scale CFD.
In-silico assessment allows to study simplifications; faithfully re-producing the arc structure (Fig. 9F), with gradual feeding followed by consumption, requires the industrial Cx= 55 g/kg. This may be
chal-lenging to operate withfilamentous fungi at SD scales. It is found that administering the pulse instantaneously and reducing the consumption rate by halving Cx(Fig. 9G) does not affect the predicted cellular
re-sponse, while theoretically reducing the viscosity by a factor 5.6[26]. Furthermore, withfinite mixing times of 2.2 and 13.2 s the predicted metabolic response still equals ideal mixing, implying the SD-protocol is realizable in practice. If lower Cxand/or regularly spaced pulse
inter-vals are required, faithful downscaling based on extra-cellular varia-tions is not possible for the simulated condivaria-tions. An alternative is to conduct“intra-cellular downscaling”: finding a set of inputs that re-plicates the metabolic model output (qp,μ and/or intra-cellular pool
Fig. 9. A: Lifeline for case III. B, C: Arc analysis approach with qref= 0.05qs,max. D: Arc time distributions based on all lifelines. E: Relation between arc-duration and
magnitude. The solid line gives the mean magnitude as function of duration. F: SD-lifeline generated from the arc-time distributions and magnitude. Feeding is gradual, with Cx= 55 g/kg. G: Alternative feed protocol with instantaneous pulses allows to set Cx= 27.5 g/kg. A comparison with the metabolic model of Tang
sizes) within set constraints (fixed feed period, pre-set Cx). This can be
done by feeding intra-cellular target values and constraints to a genetic algorithm for optimization[69], for brevity a further discussion is in Appendix B.
Besides downscaling, CFD-CRD (Computational Reaction Dynamics) for design optimization and full fed-batch simulations are explored. Moving the feed from the top to the impeller discharge stream reduces the projected yield loss from 31% to 9%. To conclude, fed-batch si-mulations with C2-coupling were conducted; the predicted yield-drop agrees with experimental observations. The model predicts significant heterogeneity in the number of glucose transporters depending on the cellular history. This observation has not yet been experimentally as-sessed, but it shows the capacity of simulations to generate hypothesis for experimental follow-up. All-together, this case highlights the spec-trum of applications of coupled CFD-CRD[33].
5.2. Other Euler–Lagrange studies
The major aspect not addressed in the Sino-Dutch collaboration is C3-coupling with metabolic models; guidelines for C3-coupling were explored with a Monod-model, while the P. chrysogenum model lacked rapid feedback and did not warrant the computationally expensive approach. The seminal work of Lapin et al. remains the sole example of C3-coupling, exemplified by replicating experimentally observed NADH oscillations in a yeast culture [1] and gradients in the phosphoe-nolpyruvate to pyruvate ratio in E. coli[2]. C1-tracking has been used more frequently; McClure et al. quantified the duration of substrate fluctuations in a bubble column[32]. Kuschel et al. simulated an E. coli cultivation and used regime analysis for transitions between replication regimes; they linked the transition frequency to experimentally tran-scriptome changes[38]. Liu et al.[70]used Euler–Lagrange simulations
to study shear exposure of plant cells, linking the death rate to peak shear exposures in a lifeline extension of the EDCF function[71]. De-lafosse et al. used a CFD-based compartment model, to compare tra-jectories with experimentally recorded tracks[54].
6. Towards novel scale down simulators
Current generation downscaling approaches have inherent limita-tions: all cells observe the same in single-vessel systems, and both in single and multi-compartment systems, a limited range of frequencies/ amplitudes is available, with no straightforward way to reduce glucose concentration other than by consumption alone [39,22](for oxygen, there is the option of stripping with an inert gas). While current scale-down simulators operating at industrial Cxmay capture dynamics at the
average level, more rapid dynamics at the individual level are un-attainable. Formulating the desires for a new generation of SD-simu-lators:
•
Decouple rate-of-change from consumption.•
Full range of amplitudes can be imposed.•
Full range of durations can be imposed.•
Arbitrary duration (circulation) patterns.•
Different cells undergo different experiences.It would additionally be desirable if individual or small sub-popu-lations of cells could be analyzed, to gain insight in individual experi-ences rather than population averages. It may be challenging to satisfy all the above, but there are promising developments in microfluidic development[69]
6.1. Towards microfluidic downscaling
Microfluidic devices allow to conduct small-population or even single-cell cultivation. They have been applied to study growth under controlled conditions[72–76], and to study enzyme dynamics during slow extra-cellular variations [77,78], but not to impose rapid fluc-tuations (to our knowledge). Convectiveflow can renew volumes in less than a second in typical micro-cultivations, and with a diffusion coef-ficient = ×6 10− m /s
m 10 2 (glucose in water), diffusive transport over
a typical 50μm lengthscale takes around a second. Theoretically, all but the most rapid turbulent variations can be imposed. By mixing an oxygen saturated and oxygen-free stream, DOfluctuations may also be imposed. Laminarflow means there will be axial dispersion, which has to be quantified.
A 2D CFD simulation of a feed pulse in a microfluidic system was done as a proof-of-principle, the geometry being a simplification of a cell-trap design[79], shown inFig. 10A. A 40μm diameter was set for the feed- and mainflow channels. Both inlets had a liquid inlet velocity of 5–20 × 10−4m/s, making the main channel velocity 1–4 mm/s. The
path length was approx. 3 mm/s. First, theflow was solved in steady state with purefluid B. At t = 0, inlet 1 is switched to pure A, with
properties equal to B and a molecular diffusion coefficient
= ×6 10− m /s
m 10 2 .Fig. 10A shows a snapshot of the transition from B to A + B. After 3 s, inlet 1 switched back to B. The concentration of A at the cell cluster was measured. While the imposed step is indeed dis-persed somewhat, a change from 5 to 95% saturation takes just 0.46–1.15 s (Fig. 11), for the highest and lowest simulated velocity. Still, if required, backmixing may be further reduced by using droplet-based microfluidics to enhance plug flow (Fig. 10B), albeit at the drawback that cells grow in compartments separate from the mainflow, requiring diffusion-controlled substrate transfer [80]. In both cases, with a transfer function describing the device dynamics, it can be checked if convolution of imposed dynamics by device dynamics affects the metabolic response. Deconvolution of the device dynamics may be possible if lifeline dynamics are slower than device dynamics. The above conceptual analysis concludes that inducing rapidfluctuations using microreactors is feasible. The challenge likely lies in quantifying the dynamic metabolic response.
6.1.1. Single-cell analysis: towards experimental lifelines?
Several microfluidic studies have focused on visual measurements,
Fig. 10. A: 2-D CFD simulation of a simple single-phase micro-reactor (channel diameter 40μm) designed to impose rapid extra-cellular variations on micro-organisms, physically trapped in a cell chamber. Contour plot shows a switch from purefluid B to fluid A fed via feed 1, while feed 2 constantly feedsfluid B. B: concept layout of a 2-phase microreactor de-signed to reduce backmixing in the substrate flow by introducing Taylor-flow.
such as cell division frequency orfluorescent labeling[81,73,76]. These techniques could be employed to study the response of (repetitive) extra-cellular perturbations on the single cell level. Still, many dimen-sions of the cellular response cannot (currently) be measured on-line, while knowledge about the dynamics of these responses is relevant for modeling purposes, where the overall cellular performance (e.g. growth
rate, enzyme expression, production rates) is computed based on the dynamics of a number of intra-cellular pools (metabolites, redox fac-tors, enzymes, etcetera).
Cell-sorting techniques, such as Fluorescence Activated Cell Sorting
[82,83](FACS), offer the option to segregate cells in a (micro)popula-tion based. For example, this can be used to segregate a cell-sample in
Fig. 11. Response at the cell-trap (black) to a pulse of A in-serted via one inlet. The pulse profile (red) has been visualized as the response for a perfect plugflow with the same super-ficial flow-rate. A: inlet velocity of 2 mm/s, B: inlet velocity of 0.5 mm/s per channel.
Fig. 12. Conceptual outline of activity segregated lifeline ana-lysis, using three activity classes (called X, Y, Z), defined pre-perturbation. A: Microcultivation containing cells with a fast-re-spondingfluorescent label. Samples are taken at regular intervals; an extra-cellular perturbation is introduced between samples−1 and 0. B: The samples are sorted and segregated based on their activity class, e.g. using FACS. Due to the fast-responding label, thefluorescence-activity of each class may change in time; this is accounted for by on-linefluorescence monitoring of cells that were grouped in the same activity class. C: After sorting, the es-tablished sub-populations are analyzed with off-line techniques, for example to measure the abundance of key metabolites (or metabolite classes). The graph illustrates how the dynamic abundance of a key metabolite can be monitored, segregated on fluorescent activity. If the correlation between label activity and composition is strong, this allows to approximate quantitative experimental lifelines. A prerequisite is that samples can be sorted and quenched sufficiently fast to allow metabolite ana-lysis.
Fig. 13. Multi-scale bioreactor evaluation and design approach. Micro-fluidics allow single cell or small population analysis. Combined with omics-information from the other lab-scales, this provides a comprehensive modeling basis. The microtiter-scale and/or massive parallel microfluidics provide the environment for rapid parallel strain evaluation and selection, which may make the shake-flask obsolete. The regular bench-scale remains valuable for verification and population behavior. Large-scale (CFD) data feeds environmental information, and is the target for reactor optimization. For well-defined processes, the intermediate pilot scale may in time be omitted.
several classes of cells after a perturbation, based on theirfluorescent activity. A range of off-line techniques can then be used for an in-depth analysis of the sub-population characteristics (sub-population omics
[84,83]). This allows, for example, to correlate the availability of cer-tain key-metabolites after a perturbation with the strength of the on-line (e.g. fluorescent) response. Such data, providing links between metabolite levels and fluorescent response, is highly useful for mod-eling purposes. In principle, it is possible quantify certain compound classes or specific metabolites at the single-cell level [85], although their size limits such possibilities for microbes. In cases where the metabolite abundance is sufficient, the variance within a given sub-population could be assessed to comment on the intrinsic heterogeneity in such pools, as well as the strength of the correlation to the on-line measured response.
The above methodology can be applied in“classical” bench-scale cultivations, as the activity at time of sorting can be correlated with subsequent measurement data. However, such cultivations do not allow to follow the response of individual cells on-line, meaning that we can only link metabolite levels between timesteps if thefluorescent activity of sub-populations is static: for example, how do metabolite levels after a substrate pulse develop in relation to the abundance of a certain (labeled) enzyme, the level of which does not change in the studied timeframe. If thefluorescent activity itself changes in response to the pulse, for example using fast-responding bio-sensors to monitor the availability of a single key metabolite or signal molecule[86,45,87], it is not possible to reconstruct dynamics: we don’t know how the activity of a cell at timestep t relates to that at step t− 1. In micro-fluidic cul-tivations establishing such correlations is possible, as we can track in-dividual cells in time, thereby allowing to construct dynamic“activity classes”. Again, we can take (off-line, destructive) samples of sorted sub-populations and correlate the levels of measured metabolites with thefluorescent response, but now we can also use the on-line single-cell responses to monitor how sub-populations develop in time, in terms of theirfluorescent activity and (by correlation) other relevant metabolic pools. If population heterogeneity is already visible prior to a pertur-bation, it is possible to monitor how different sub-populations respond to a perturbation (how strong is the response, and is it consistent be-tween cells in a class?), in relation to their pre-perturbation composi-tion. Graphically, this approach is outlined inFig. 12. Of course, we are still typically dealing with sub-population averages, and sub-popula-tions will still be intrinsically heterogeneous, hence individual lifelines of cells in the sub-population may still differ. Still, being able to measure and dynamically model extrinsic heterogeneity between sub-popula-tions[45]would be a large step forward from the current situation, where models are based on whole-population averages. Significant experimental developments are required to realize experimental lifeline acquisition as described above, but the insights might well be worth pushing for, and we’ve only just begun to explore the opportunities. 7. Future perspective
Euler–Lagrange CFD simulations, especially with C3 coupling, re-quire days to weeks of computation time to resolve several mixing times offlow time. This limits its applications; routine usage in process design and optimization is not yet feasible. Computational advances such as highly parallelized simulations and GPU-assisted CFD are promising: real-time Euler–Lagrange simulations are already feasible if the flow itself can be considered steady-state, and ANSYS Discovery live already delivers real-time capabilities for transient single-phaseflows. If this is extended to multi-phase and reactive flows, application for routine design and optimization are truly within reach. Even if crude meshes are applied, the ability to assess the qualitative impact of design changes and subsequently case the most promising configurations in detail is of great interest. Next-generation microfluidic SD-simulators may then impose lifelines from such simulations directly on organisms, while advances in micro-population omics allow for ever-increasing
insight in the cellular response. Our current understanding of extra-cellularfluctuation effects is insufficient to go directly from lab to full scale, but with increasing understanding the need for piloting may re-duce. Automated microwell devices like the Biolector[88]are replacing shake-flasks, allowing for parallelized strain testing, and in due time milli-wells may be replaced by micro-droplet cultivations, converging the capacity of microfluidic micro-population and/or even single-cell analysis with highly parallel strain testing and controlled population heterogeneity studies. The combined information from massive parallel micro-cultivations and bench-scale cultivations refines metabolic un-derstanding and associated CFD-CRD models; these simulations in turn provide environmental input for cultivations at all scales, illustrated in
Fig. 13. Real-time simulations also have implications for process con-trol. Imagine a“digital twin” of an industrial reactor, predicting the effect of process disturbances on the metabolic composition; with vi-sualization techniques like virtual reality, operators may even“look” into running processes on-site. While this may be some longer term speculation, the potential is certainly there. Advances in simulation, miniaturization and visualization will impact the field of bioprocess engineering in the years to come, and the biggest advances are likely still ahead.
Acknowledgements
We gratefully acknowledge the input from our colleagues at the DSM biotechnology center, Delft university of Technology, and colla-borators at East China University of Technology. We are especially thankful to Wenjun Tang, Guan Wang, Wouter van Winden, Walter van Gulik, Amit Deshmukh, Ju Chu, Jianye Xia, Matthias Reuss and Sef Heijnen for their frequent discussions and contributions to the colla-borative project that was the origin of this work. This work has been conducted within a multi-party research project, among DSM-Sinochem Pharmaceuticals, TU Delft, East China University of Science and Technology and Shanghai Guojia Ltd., funded by NWO and MoST (2013DFG32630).
Appendix A. Supplementary data
Supplementary data associated with this article can be found, in the online version, athttps://doi.org/10.1016/j.bej.2018.09.001.
References
[1] A. Lapin, D. Müller, M. Reuss, Dynamic behavior of microbial populations in stirred bioreactors simulated with Euler–Lagrange methods: traveling along the lifelines of single cells, Ind. Eng. Chem. Res. 43 (2004) 4647–4656.
[2] A. Lapin, J. Schmid, M. Reuss, Modeling the dynamics of E. coli populations in the three-dimensional turbulentfield of a stirred-tank bioreactor—a structured-segre-gated approach, Chem. Eng. Sci. 61 (2006) 4783–4797.
[3] A. Humphrey, Shakeflask to fermentor: what have we learned? Biotechnol. Progr. 14 (1998) 3–7.
[4] P. Neubauer, S. Junne, Scale-down simulators for metabolic analysis of large-scale bioprocesses, Curr. Opin. Biotechnol. 21 (2010) 114–121.
[5] G. Wang, W. Tang, J. Xia, J. Chu, H.J. Noorman, W.M. van Gulik, Integration of microbial kinetics andfluid dynamics toward model-driven scale-up of industrial bioprocesses, Eng. Life Sci. 15 (2015) 20–29.
[6] N.M. Oosterhuis, N.W. Kossen, Dissolved oxygen concentration profiles in a pro-duction-scale bioreactor, Biotechnol. Bioeng. 26 (1984) 546–550.
[7] A.P.J. Sweere, Y.A. Matla, J. Zandvliet, K.C.A.M. Luyben, N.W.F. Kossen, Experimental simulation of glucosefluctuations, Appl. Microbiol. Biotechnol. 28 (1988) 109–115.
[8] A.P.J. Sweere, J.R. Mesters, L. Janse, K.C.A.M. Luyben, N.W.F. Kossen, Experimental simulation of oxygen profiles and their influence on baker's yeast production: I. One-fermentor system, Biotechnol. Bioeng. 31 (1988) 567–578. [9] A.P.J. Sweere, L. Janse, K.C.A.M. Luyben, N.W.F. Kossen, Experimental simulation
of oxygen profiles and their influence on baker's yeast production: II. Two-fer-mentor system, Biotechnol. Bioeng. 31 (1988) 579–586.
[10] A. Lemoine, M.H. Limberg, S. Kästner, M. Oldiges, P. Neubauer, S. Junne, Performance loss of Corynebacterium glutamicum cultivations under scale-down conditions using complex media, Eng. Life Sci. (2016).
[11] A.-L. Heins, R. Lencastre Fernandes, K.V. Gernaey, A.E. Lantz, Experimental and in silico investigation of population heterogeneity in continuous Sachharomyces
